text:
module Day22

using AdventOfCode2019

function day22(input::String = readInput(joinpath(@__DIR__, "..", "data", "day22.txt")))
    instructions = parse_instructions(input)
    p1 = part1(instructions)
    p2 = part2(instructions)
    return [p1, p2]
end

function part1(instructions)
    n = 10_006
    i = 2019
    for (inst, N) in instructions
        if inst == :stack
            i = n - i
        elseif inst == :cut
            i = mod(i - N, n + 1)
        elseif inst == :increment
            i = mod(i * N, n + 1)
        end
    end
    return i
end

function part2(instructions)
    i = BigInt(2020)
    n = BigInt(119_315_717_514_047) - 1
    niter = BigInt(101_741_582_076_661)

    # Observation: every single shuffle instruction from the input
    # can be represented as `g(i) = a*i + b (modulo n + 1)`, so
    # the result of applying all the shuffle instructions from the
    # input can also be represented as `f(i) = a*i + b (modulo n + 1)`.
    # We need to find the integers `a` and `b`. Set `x = i`, `y = f(x)`
    # and `z = f(y)`. Then `y = a*x + b` and `z = a*y + b`,
    # so `a = (x - y)^{-1} * (y - z) (modulo n + 1)` where
    # `(x - y)^{-1}` is the modular inverse and `b = y - a*x (modulo n + 1)`
    x = i
    y = apply_reverse_instructions(instructions, x, n + 1)
    z = apply_reverse_instructions(instructions, y, n + 1)
    a = mod((y - z) * invmod(x - y, n + 1), n + 1)
    b = mod(y - a * x, n + 1)

    # Now we need to apply `f(i) = a*i + b` `k = niter` times.
    # Observe that
    # f^k(i) = a^k * i + b * (a^{k-1} + a^{k-2} + ... + a + a^0) mod n + 1
    #        = a^k * i + b * (a^k - 1) / (a - 1) mod n + 1
    # This allows us to calculate `f^k(i)` directly without a loop.
    pm = powermod(a, niter, n + 1)
    c = mod(mod(pm - 1, n + 1) * invmod(a - 1, n + 1), n + 1)
    i = mod(pm * i + b * c, n + 1)

    return i
end

function apply_reverse_instructions(instructions, i, decksize)
    # Apply all the instructions from the input in reverse order
    # to an index `i`
    n = decksize - 1
    for (inst, N) in reverse(instructions)
        if inst == :stack
            i = n - i
        elseif inst == :cut
            i = mod(i + N, n + 1)
        elseif inst == :increment
            i = mod(invmod(N, n + 1) * i, n + 1)
        end
    end
    return i
end

function parse_instructions(input::String)
    instructions = Array{Tuple{Symbol,Union{Int,Nothing}},1}()
    newStackRegex = r"deal\sinto\snew\sstack"
    cutRegex = r"cut\s(-?[0-9]+)"
    incrementRegex = r"deal\swith\sincrement\s([0-9]+)"
    for line in split(strip(input), "\n")
        m = match(newStackRegex, line)
        if !isa(m, Nothing)
            push!(instructions, (:stack, nothing))
            continue
        end
        m = match(cutRegex, line)
        if !isa(m, Nothing)
            val = parse(Int, m.captures[1])
            push!(instructions, (:cut, val))
            continue
        end
        m = match(incrementRegex, line)
        if !isa(m, Nothing)
            val = parse(Int, m.captures[1])
            push!(instructions, (:increment, val))
            continue
        end
    end
    return instructions
end

end # module

ext: jl


text:
"""
   HeistEvent is an object that represents a gallery event along with a counter that
   advances when an event is seen.

   You can make a HeistEvent by passing in a `HeistFiles` object that represents the
   files to process.

   If you are looping over events, you should create a `HiestEvents` object (plural) instead
   of this one.
"""
mutable struct HeistEvent
   galleryEvent::Cxx.CppValue
   nSeen::Int

   HeistEvent(hf::HeistFiles) = new( icxx"gallery::Event ev($(hf.theVector)); return ev;", 1)
end

"""
   nextEvent advances to the next event in the file or to the next file to process.
"""
nextEvent(he::HeistEvent) = ( @cxx he.galleryEvent->next() ; he.nSeen += 1; he )

"""
   atEnd returns true if there are no more events to process
"""
atEnd(he::HeistEvent) = @cxx he.galleryEvent->atEnd()


text:
# include("../les/custom_avg.jl")
# include("../les/output_to_data.jl")
# include("../les/get_les_data.jl")

function closure_free_convection_tke_full_evolution(parameters, N, Δt, les;
                                 subsample = 1, grid = 1)

     # set parameters
     # parameters = TKEMassFlux.TKEParameters( Cᴰ = Cᴰ )
     # Build the model with a Backward Euler timestepper
     constants = Constants(Float64; α = les.α , β = les.β, ρ₀= les.ρ, cP=les.cᵖ, f=les.f⁰, g=les.g)
     model = TKEMassFlux.Model(grid = UniformGrid(N, les.L), stepper=:BackwardEuler, constants = constants, tke_equation = parameters)
     # Get grid if necessary
     if grid != 1
         zp = collect(model.grid.zc)
         @. grid  = zp
     end

    # define the closure
    function free_convection()
        # get average of initial condition of LES
        T⁰ = custom_avg(les.T⁰, N)
        # set equal to initial condition of parameterization
        model.solution.T[1:N] = copy(T⁰)
        # Set boundary conditions
        model.bcs.T.top = FluxBoundaryCondition(les.top_T)
        model.bcs.T.bottom = GradientBoundaryCondition(les.bottom_T)
        # set aside memory
        if subsample != 1
            time_index = subsample
        else
            time_index = 1:length(les.t)
        end
        Nt = length(les.t[time_index])
        𝒢 = zeros(N, Nt)

        # loop the model
        ti = collect(time_index)
        for i in 1:Nt
            t = les.t[ti[i]]
            run_until!(model, Δt, t)
            @. 𝒢[:,i] = model.solution.T[1:N]
        end
        return 𝒢
    end
    return free_convection
end

function closure_free_convection_tke(parameters, N, Δt, les;
                                 subsample = 1, grid = 1)

     # set parameters
     # parameters = TKEMassFlux.TKEParameters( Cᴰ = Cᴰ )
     # Build the model with a Backward Euler timestepper
     constants = Constants(Float64; α = les.α , β = les.β, ρ₀= les.ρ, cP=les.cᵖ, f=les.f⁰, g=les.g)
     model = TKEMassFlux.Model(grid = UniformGrid(N, les.L), stepper=:BackwardEuler, constants = constants, tke_equation = parameters)
     # Get grid if necessary
     if grid != 1
         zp = collect(model.grid.zc)
         @. grid  = zp
     end

     # Set boundary conditions
     model.bcs.T.top = FluxBoundaryCondition(les.top_T)
     model.bcs.T.bottom = GradientBoundaryCondition(les.bottom_T)
     # set aside memory
     if subsample != 1
         time_index = subsample
     else
         time_index = 1:length(les.t)
     end
     Nt = length(les.t[time_index])
     𝒢 = zeros(N, n_steps+1)

     ti = collect(time_index)

    # define the closure
    function evolve_forward(; T⁰=T⁰, n_steps = 1)
        # get average of initial condition of LES
        T⁰ = custom_avg(T⁰, N)
        # set equal to initial condition of parameterization
        model.solution.T[1:N] = copy(T⁰)

        # loop the model
        for i in 1:n_steps+1
            t = les.t[ti[i]]
            run_until!(model, Δt, t)
            @. 𝒢[:,i] = model.solution.T[1:N]
        end
        return 𝒢
    end
    return evolve_forward
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using DifferentiableStateSpaceModels, Symbolics, Test

# Setup for test
@testset "Basic construction of model" begin
    ∞ = Inf
    @variables α, β, ρ, δ, σ, Ω_1, Ω_2
    @variables t::Integer, k(..), z(..), c(..), q(..)

    x = [k, z]
    y = [c, q]
    p = [α, β, ρ, δ, σ]

    H = [1 / c(t) - (β / c(t + 1)) * (α * exp(z(t + 1)) * k(t + 1)^(α - 1) + (1 - δ)),
         c(t) + k(t + 1) - (1 - δ) * k(t) - q(t), q(t) - exp(z(t)) * k(t)^α,
         z(t + 1) - ρ * z(t)]

    steady_states = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
                     c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                            δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
                     q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1))]

    steady_states_iv = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
                        c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                               δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
                        q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1))]

    n_ϵ = 1
    n_z = 2
    n_x = length(x)
    n_y = length(y)
    n_p = length(p)
    Γ = reshape([σ], n_ϵ, n_ϵ)
    η = reshape([0; -1], n_x, n_ϵ) # η is n_x * n_ϵ matrix

    Q = zeros(n_z, n_x + n_y)
    Q[1, 1] = 1.0
    Q[2, 3] = 1.0

    Ω = [Ω_1, Ω_2]

    model_name = "rbc_temp"
    verbose = true
    save_ip = true
    save_oop = false
    max_order = 2
    skipzeros = false
    fillzeros = false

    module_cache_path = make_perturbation_model(H; model_name, t, y, x, p, steady_states,
                                                steady_states_iv, Γ, Ω, η, Q,
                                                overwrite_model_cache=true, verbose,
                                                max_order, save_ip, save_oop, skipzeros,
                                                fillzeros)
    make_perturbation_model(H; model_name, t, y, x, p, steady_states, steady_states_iv, Γ,
                            Ω, η, Q, overwrite_model_cache=true)

    make_perturbation_model(H; model_name, t, y, x, p, steady_states, steady_states_iv, Γ,
                            Ω, η, Q, overwrite_model_cache=false, verbose, max_order,
                            save_ip, save_oop, skipzeros, fillzeros)
    #module_cache_path = join_path(default_model_cache_location(), "rbc_temp.jl")

    # Load the module
    include(module_cache_path)

    # Test the construction
    m = PerturbationModel(Main.rbc_temp)
    @test m.n_y == n_y
    @test m.max_order == max_order
    @test m.mod.n_z == n_z
    # Note that this is inherently dynamic and cannot be inferred, so @inferred PerturbationModel(Main.rbc_observables) would fail

    c = SolverCache(m, Val(2), [:a, :b, :c]) # the exact symbol names won't matter for inference
    @inferred SolverCache(m, Val(2), [:a, :b, :c])
    @inferred SolverCache(m, Val(2), [:a])  # less differentiated parameters shouldn't matter
    @inferred SolverCache(m, Val(1), [:a, :b, :c])  # lower order shouldn't break inference either

    # Test convenience macro.  Use dict or named tuple of args, without the model_name
    args = (; t, y, x, p, steady_states, steady_states_iv, Γ, Ω, η, Q,
            overwrite_model_cache=true, verbose, max_order, save_ip, save_oop, skipzeros,
            fillzeros)

    # If the model_nam.jl doesn't exist it creates it, inclues it, and uses that to construct
    # If it exists but the modeule name isn't loaded, then it includes it and uses the constructor
    # If the module name is loaded, it just calls the constructor
    m = @make_and_include_perturbation_model("rbc_conv2", H, args)
    m = PerturbationModel(Main.rbc_conv2) # After creation, can use module constructor or simply call this function again
    m = @make_and_include_perturbation_model("rbc_conv2", H, args) # should not regenerate
    m = @make_and_include_perturbation_model("rbc_temp", H, args) # already loaded from previous include
end

@testset "Loading Examples" begin
    # Tests for loading examples
    H, mod_vals = DifferentiableStateSpaceModels.Examples.rbc()
    make_perturbation_model(H; model_name="rbc_test_2", overwrite_model_cache=true,
                            verbose=true, mod_vals...)

    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc)
    @test m.n_y == 2

    # Load second one
    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc)
    @test m.n_y == 2

    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc_observables)
    @test m.n_y == 2

    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc_solve_steady_state)
    @test m.n_y == 2

    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc_multiple_shocks)
    @test m.n_y == 2

    m = @include_example_module(DifferentiableStateSpaceModels.Examples.rbc_observables_separate_variance)
    @test m.n_y == 2
end


text:
"""
Compute p(x) where p is the unique polynomial of degree length(xi),
such that p(x[i]) = y[i]) for all i.
Call with
```julia
aitken_neville!(ti::Vector, xi::Matrix, t::Number, x::Vector)
```
where
* `ti`: interpolation nodes
* `xi`: interpolation values
* `t`:  evaluation point
* `x`:  evaluation value
"""
function aitken_neville!(ti::Vector{TT}, xi::Matrix, t::TT, x::Vector) where {TT}
    @assert length(ti) == size(xi,2)
    @assert length(x)  == size(xi,1)

    for j in eachindex(ti)
        for i in 1:(length(ti)-j)
            for k in axes(x,1)
                xi[k,i] = xi[k,i+1] + (xi[k,i] - xi[k,i+1]) * (ti[i+j] - t) / (ti[i+j] - ti[i])
            end
        end
    end
    for k in eachindex(x)
        x[k] = xi[k,1]
    end
end


text:
module RadiiPolynomial

    using Printf, Reexport
    @reexport using IntervalArithmetic, LinearAlgebra, SparseArrays

include("utilities.jl")

# Sequence spaces

include("sequence_spaces/vector_spaces.jl")
    export VectorSpace, ParameterSpace, SequenceSpace, BaseSpace,
        TensorSpace, ⊗, Taylor, Fourier, Chebyshev,
        CartesianSpace, CartesianPower, CartesianProduct, ×

    export order, frequency, space, spaces, dimension, dimensions, indices,
        nb_cartesian_product

#

include("sequence_spaces/sequence.jl")
    export Sequence
include("sequence_spaces/linear_operator.jl")
    export LinearOperator, domain, codomain
include("sequence_spaces/broadcast.jl")

    export coefficients, eachcol, eachrow, eachcomponent, component

#

include("sequence_spaces/norm.jl")
    export Weights, GeometricWeights, geometricweights, AlgebraicWeights,
        algebraicweights, rate, weight, Norm, ℓᵖNorm, Weightedℓ¹Norm,
        CartesianPowerNorm, CartesianProductNorm, 𝐻ˢNorm, norm, opnorm

#

include("sequence_spaces/linear_algebra.jl")

#

include("sequence_spaces/arithmetic/add_conv_image.jl")
include("sequence_spaces/arithmetic/sequence.jl")
include("sequence_spaces/arithmetic/linear_operator.jl")
include("sequence_spaces/arithmetic/action.jl")
    export image, +̄, -̄, add!, radd!, ladd!, sub!, rsub!, lsub!
include("sequence_spaces/arithmetic/convolution.jl")
    export banach_rounding_order, banach_rounding!, banach_rounding_mul,
        banach_rounding_pow, *̄, banach_rounding_mul_bar, ^̄, banach_rounding_pow_bar
include("sequence_spaces/arithmetic/fft.jl")
    export fft_size, fft, fft!, ifft!, rifft!

#

include("sequence_spaces/special_operators/projection.jl")
    export project, project!
include("sequence_spaces/special_operators/multiplication.jl")
    export Multiplication
include("sequence_spaces/special_operators/calculus.jl")
    export Derivative, differentiate, differentiate!, Integral, integrate, integrate!
include("sequence_spaces/special_operators/evaluation.jl")
    export Evaluation, evaluate, evaluate!
include("sequence_spaces/special_operators/scale.jl")
    export Scale, scale, scale!
include("sequence_spaces/special_operators/shift.jl")
    export Shift, shift, shift!

# Radii polynomial

include("rpa/interval_existence.jl")
    export interval_of_existence
include("rpa/newton.jl")
    export newton, newton!

end


text:
using EquationsOfStateOfSolids
using Test

@testset "EquationsOfStateOfSolids.jl" begin
    include("collections.jl")
    include("Fitting.jl")
end


text:
using PyPlot
using MPI
using Printf

"""Gradient function. 2nd-order centred."""
function ∂x!(du::AbstractVector,u::AbstractVector)
    m = length(du)
    @boundscheck m+1 == length(u) || throw(BoundsError())

    @inbounds for i ∈ 1:m
        du[i] = one_over_dx*(u[i+1] - u[i])
    end
end

"""2-point interpolation."""
function Ix!(du::AbstractVector,u::AbstractVector)
    m = length(du)
    @boundscheck m+1 == length(u) || throw(BoundsError())

    @inbounds for i ∈ 1:m
        du[i] = 0.5*(u[i] + u[i+1])
    end
end

# Forcing
F(x,t) = -F0*sin.(8π*x/Lx .+ 200*t).*sin.(π*x/Lx).^2/rho

"""Computes the right-hand side of
        ∂u/∂t = -u*∂u/∂x -g*∂η/∂x + ν∂²u/∂x² + Fx
        ∂η/∂t = -∂(uh)/∂x. """
function rhs!(du,dη,u,η,h,u²,h_u,u²_h,dudx,p,dpdx,U,dUdx,t)

    @views h .= η .+ H      # layer thickness
    @views u² .= u.^2

    Ix!(h_u,h)
    Ix!(u²_h,u²)
    ∂x!(dudx,u)

    # Bernoulli potential + stress "tensor"
    @views p .= -.5*u²_h .- g*η[2:end] .+ ν*dudx

    # momentum equation
    ∂x!(dpdx,p)
    @views du[2:end-1] .= dpdx .+ F(x_up,t)./h_u[2:end]

    # continuity equation
    @views U .= u[1:end-1].*h_u         # volume flux
    ∂x!(dUdx,U)
    @views dη[2:end-1] .= -dUdx
end

"""Size of the subdomain per process."""
function subdomain_size(N::Int,size::Int)
    if mod(N,size) == 0
        return N ÷ size
    else
        throw(error("Domain size must be divisible by number of processes."))
    end
end

"""Rank of the neighbouring processes with wrap-around behaviour (periodic boundary conditions)."""
function neighbours(size::Int,rank::Int)
    # left, right neighbour
    lnb,rnb = mod(rank-1,size),mod(rank+1,size)
end

"""Add halo to the left and right of each subdomain for ghost point communication."""
function add_halo(u::AbstractVector,η::AbstractVector,rmsg::AbstractVector,lmsg::AbstractVector)
    u = cat(0.,u,0.,dims=1)
    η = cat(0.,η,0.,dims=1)

    ghost_points!(u,η,rmsg,lmsg)
    return u,η
end

"""Propagate the boundary conditions (BC) via ghost point communication.
For multiple processes, BC of subdomains are communicated via MPI."""
function ghost_points!(u::AbstractVector,η::AbstractVector,rmsg::AbstractVector,lmsg::AbstractVector)
    if size == 1
        # sequential, periodic BC
        u[1] = u[end-1]
        u[end] = u[2]

        η[1] = η[end-1]
        η[end] = η[2]

    else
        # parallel with communication
        MPI.Isend([u[2],η[2]],lnb,rank+10,comm)
        MPI.Isend([u[end-1],η[end-1]],rnb,rank+20,comm)

        rreq = MPI.Irecv!(rmsg,rnb,rnb+10,comm)
        lreq = MPI.Irecv!(lmsg,lnb,lnb+20,comm)

        MPI.Waitall!([rreq,lreq])

        u[end],η[end] = rmsg
        u[1],η[1] = lmsg
    end
end

"""Preallocate intermediate variables used in the right-hand side computation."""
function preallocate(u::AbstractVector,η::AbstractVector)
    # u, η already have the halo
    u0,η0 = zero(u),zero(η)
    u1,η1 = zero(u),zero(η)
    du,dη = zero(u),zero(η)

    u²,h = zero(u),zero(η)
    p,dpdx = zeros(length(η)-1),zeros(length(η)-2)

    # each centred gradient/interpolation removes one grid point
    h_u = zeros(length(η)-1)
    u²_h = zeros(length(u)-1)
    dudx = zeros(length(u)-1)

    U = zeros(length(u)-1)
    dUdx = zeros(length(u)-2)

    return u0,η0,u1,η1,du,dη,h,u²,h_u,u²_h,dudx,p,dpdx,U,dUdx
end

"""4th order Runge-Kutta time integration."""
function time_integration(Nt,u,η)

    # preallocate for parallel: message from right, left
    rmsg = Array{Float64}(undef,2)
    lmsg = Array{Float64}(undef,2)

    # pre-allocate memory
    u,η = add_halo(u,η,rmsg,lmsg)
    u0,η0,u1,η1,du,dη,h,u²,h_u,u²_h,dudx,p,dpdx,U,dUdx = preallocate(u,η)

    t = 0.

    # for output
    u_out = zeros(Nt+1,length(u)-2)     # -2 to remove ghost points
    η_out = zeros(Nt+1,length(η)-2)

    # store the initial conditions
    u_out[1,:] = u[2:end-1]
    η_out[1,:] = η[2:end-1]

    for i = 1:Nt

        ghost_points!(u,η,rmsg,lmsg)
        u1 .= u
        η1 .= η

        for RKi = 1:4
            if RKi > 1
                ghost_points!(u1,η1,rmsg,lmsg)
            end

            rhs!(du,dη,u1,η1,h,u²,h_u,u²_h,dudx,p,dpdx,U,dUdx,t)

            if RKi < 4 # RHS update for the next RK-step
                u1 .= u .+ RKβ[RKi]*dt*du
                η1 .= η .+ RKβ[RKi]*dt*dη
            end

            # Summing all the RHS on the go
            u0 .+= RKα[RKi]*dt*du
            η0 .+= RKα[RKi]*dt*dη

        end

        u .= u0
        η .= η0
        t += dt

        # store for output
        u_out[i+1,:] .= u[2:end-1]
        η_out[i+1,:] .= η[2:end-1]

    end
    u_out,η_out
end

MPI.Init()
comm = MPI.COMM_WORLD

const rank = MPI.Comm_rank(comm)
const size = MPI.Comm_size(comm)

const g = 10.
const H = 10.
const N = 500
const Nt = 1000
const cfl = 0.9
const rho = 1.
const ν = 0.006
const F0 = 500.
const Lx = 1.
const dx = Lx/N
const one_over_dx = 1/dx
const cph = sqrt(g*H)
const dt = cfl * dx / cph
const RKα = [1/6.,1/3.,1/3.,1/6.]
const RKβ = [0.5,0.5,1.]

# grid
const x_h = dx/2:dx:Lx
const x_u = dx:dx:Lx        # no u-point at x=0 but at x=L (periodicity)

# subdomains for parallel
const Np = subdomain_size(N,size)
const x_up = x_u[(rank*Np+1):(rank+1)*Np]
const lnb,rnb = neighbours(size,rank)

# initial conditions
u_ini = fill(0.,Np)
η_ini = fill(0.,Np)

u,η = time_integration(Nt,u_ini,η_ini)

# receive solutions from all processors
ηall = MPI.Gather(η,0,comm)
if rank == 0
    ηall = reshape(ηall,Nt+1,N)
end

MPI.Barrier(comm)
MPI.Finalize()

## PLOTTING
ioff()
if rank == 0
    fig,ax = subplots(figsize=(6,3))

    l1, = ax[:plot](x_h,ηall[1,:],lw=2)
    ax[:set_xlim](0,1)
    ax[:set_ylim](-1.2,1.7)
    ax[:set_yticks]([-1,0,1])
    ax[:set_xticks]([0,1])
    ax[:set_xlabel]("x")
    ax[:set_ylabel]("η")
    tight_layout()

    ax[:plot]([0,.504],[1.3,1.3],"C1")
    ax[:plot]([0.996,1.0],[1.3,1.3],"C1")

    ax[:plot]([0.496,1],[1.33,1.33],"C2")
    ax[:plot]([0.0,0.004],[1.33,1.33],"C2")

    ax[:text](0.1,1.4,"Process 0")
    ax[:text](0.8,1.4,"Process 1")

    #global iframe = 0
    for it = 2:2:1000
        pause(0.00001)
        l1[:set_data](x_h,ηall[it,:])
        #savefig("/home/kloewer/julia/vid/waveframes/frame"*Printf.@sprintf("%04d",iframe)*".png",dpi=75)
        #iframe += 1
    end

    pause(2)
    close(fig)
end

max_stars_repo_name: milankl/misc

max_issues_repo_name: milankl/misc

max_forks_repo_name: milankl/misc


text:
function ellk_bulirsch(k2::Real)
ca = sqrt(eps(k2))
kc=sqrt(1.0-k2)
h=1.0
m=1.0
while abs(h-kc) > (ca*h)
  h  = m
  m += kc
  kc = sqrt(h*kc)
  m  = 0.5*m
#  println(h,' ',m,' ',kc,' ',h-kc)
end
h=m
m+=kc
return pi/m
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
#-----------------------------------------------------------------------
#   TEST NULL ENFORCER FUNCTIONS
#-----------------------------------------------------------------------

X   = rand(4)

v       = zeros(4)
v[2:4]  = rand(3)
vnorm  = ηdot(v,v)

nvf = nullenforcerf( v , X , η )
nvp = nullenforcerp( v , X , η )

@test nvf[1] > 0.
@test ηdot(nvf,nvf) ≈ 0 atol=1e-15

@test nvp[1] < 0.
@test ηdot(nvp,nvp) ≈ 0 atol=1e-15

#-----------------------------------------------------------------------
#   TEST HAMILTONIAN
#-----------------------------------------------------------------------

Z = [0.;0.;0.;0.;1.;0.;0.;0.]
@test HamGeo( Z , η ) ≈ -0.5

Z = [0.;0.;0.;0.;0.;1.;0.;0.]
@test HamGeo( Z , η ) ≈ 0.5

v3 = rand(3)
Z[6:8] = v3
@test HamGeo( Z , η ) ≈ dot(v3,v3)/2

#-----------------------------------------------------------------------
#   TEST SYMPLECTIC OPERATOR
#-----------------------------------------------------------------------

Z = ones(8)

@test Jsympl(zeros(8)) ≈ zeros(8)
@test Jsympl(π*Z)[1:4] ≈ π*ones(4)
@test Jsympl(π*Z)[5:8] ≈ -π*ones(4)

#-----------------------------------------------------------------------
#   TEST HAMILTON EQUATIONS
#-----------------------------------------------------------------------

@test ZdotGeo( π*Z , η )[1]     ≈ -π
@test ZdotGeo( π*Z , η )[2:4]   ≈ π*ones(3)
@test ZdotGeo( π*Z , η )[5:8]   ≈ zeros(4)

#-----------------------------------------------------------------------
#   TEST SOLVEZ FUNCTION
#-----------------------------------------------------------------------

Z0 = [0.;0.;0.;0.;1.;1.;0.;0.]      # Index raised in p's here
ZF = [1.;1.;0.;0.;-1.;1.;0.;0.]     # Index lowered in p's here
@test solveZ( Z0 , η , 1e-9 , 1e-9 , AutoVern7(Rodas5()) ) ≈ ZF

Z0 = [0.;0.;0.;0.;1.;0.;1.;0.]      # Index raised in p's here
ZF = [1.;0.;1.;0.;-1.;0.;1.;0.]     # Index lowered in p's here
@test solveZ( Z0 , η , 1e-9 , 1e-9 , AutoVern7(Rodas5()) ) ≈ ZF

Z0 = [0.;0.;0.;0.;1.;0.;0.;1.]      # Index raised in p's here
ZF = [1.;0.;0.;1.;-1.;0.;0.;1.]     # Index lowered in p's here
@test solveZ( Z0 , η , 1e-9 , 1e-9 , AutoVern7(Rodas5()) ) ≈ ZF


text:
module TestPropagator

using Test
using PIMD, LinearAlgebra

ω, Δt, m = 2.3, 0.1, 1.4
A = [0 inv(m); - m * ω^2 0]

@testset "Propagator" begin
    @test exp(A * Δt) ≈ [cos(ω * Δt) sin(ω * Δt) * inv(ω * m); -sin(ω * Δt) * (ω * m) cos(ω * Δt)]
    @test sqrtcayley(A * Δt) ≈ [1 Δt / m; -ω^2 * Δt * m 1] / √(1 + ω^2 * Δt^2)
    @test cayley(A * Δt) ≈ [1 - ω^2 * Δt^2 / 4 Δt / m; -ω^2 * Δt * m 1 - ω^2 * Δt^2 / 4] / (1 + ω^2 * Δt^2 / 4)
end

end


text:
# ## Surface
# ### Surface (1)

#md # ![FILE_NAME.png](images/FILE_NAME.png)

f(x,y) = x^2 + y^2
x = -10:10
y = x
surface(x, y, f)


text:
"""
    AMatrix{T} <: AbstractMatrix{T}

This type is for testing `@isinsrc` and `@insrc`. We implement
enough of the `AbstractMatrix` interface to test these macros.
"""
struct AMatrix{T} <: AbstractMatrix{T}
    n::Int
end

# ::AMatrix is square with sides of length `n`
Base.size(a::AMatrix) = (a.n, a.n)

# elements of ::AMatrix are all equal to 1
function Base.getindex(a::AMatrix{T}, i, j) where {T}
    @boundscheck checkbounds(a, i, j)
    return one(T)
end

# Convert to an allocated `Matrix`
Base.Matrix(a::AMatrix{T}) where {T} = ones(T, a.n, a.n)

# We implement an efficient method for computing the sum of the elements.
Base.sum(a::AMatrix{T}) where T = T(a.n^2)

# We implement an efficient method for computing the product of the elements.
# BUT WE FORGET TO PREPEND `Base.`. So this method doe not extend the function in
# `Base` and is not visible unless qualified with the package name.
prod(a::AMatrix{T}) where T = one(T)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type InitializeStage
	InitializeStage()=begin
		NComp=outers.NComp
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components"
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Tray Temperature at Column Top",
				:Default=>300
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Tray Temperature at Column Bottom",
				:Default=>300
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Component Molar Fraction at Top"
			)),(NComp)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Component Molar Fraction at Bottom"
			)),(NComp)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Column holdup",
				:Default=>0.1
			)),
			[:NComp,:TopStageTemperature,:BottomStageTemperature,:TopStageComposition,:BottomStageComposition,:ColumnHoldup,],
		)
	end
	NComp::DanaInteger
	TopStageTemperature::temperature
	BottomStageTemperature::temperature
	TopStageComposition::Array{fraction}
	BottomStageComposition::Array{fraction}
	ColumnHoldup::positive
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export InitializeStage
function atributes(in::InitializeStage,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Initial conditions for Packed column section."
	drive!(fields,_)
	return fields
end
InitializeStage(_::Dict{Symbol,Any})=begin
	newModel=InitializeStage()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# snoop.jl
#
# Author: Tom Price
# Date: December 2018
#
# snoop on julia compilation
# see https://github.com/timholy/SnoopCompile.jl
#
using SnoopCompile

# Log the compiles
# This only needs to be run once to generate "/tmp/qpcr_compiles.csv"
SnoopCompile.@snoop "/tmp/qpcr_compiles.csv" begin
	cd("/home/vagrant/chaipcr/bioinformatics/QpcrAnalysis")
	push!(LOAD_PATH,pwd())
	using QpcrAnalysis
	include("../test/test_functions.jl")
	results = test_dispatch()
end

# Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts
data = SnoopCompile.read("/tmp/qpcr_compiles.csv")
pc = SnoopCompile.parcel(reverse!(data[2]))

# remove test functions not used in production mode
filter!(
    x -> (match(r"QpcrAnalysis\.verify_re",x) == nothing),
    pc[:QpcrAnalysis])
filter!(
    x -> (match(r"QpcrAnalysis\.\w+_test",x) == nothing),
    pc[:QpcrAnalysis
filter!(
    x -> (match(r"FactCheck",x) == nothing),
    pc[:Base])
filter!(
    x -> (match(r"FactCheck",x) == nothing),
    pc[:Core])
delete!(pc,:FactCheck)
delete!(pc,:Main)

# write output files to be included in QpcrAnalysis.jl
SnoopCompile.write("/tmp/precompile", pc)

filter(x -> (match(r"FactCheck",x) != nothing),p)

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#=##############################################################################
# DESCRIPTION
    Utilities

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

# Color pallete for plots
clrs = [
            # [17, 20, 76]/255,
            [58, 150, 144]/255,
            [250, 188, 96]/255,
            [225, 98, 98]/255,
            [194, 232, 206]/255,
            [148, 170, 42]/255,
            [42, 148, 170]/255,
            [170, 42, 84]/255,
            [113, 145, 146]/255,
            [252, 127, 178]/255
        ]

"""
    `plot_maneuver(maneuver::KinematicManeuver; ti::Real=0, tf::Real=1,
vis_nsteps=300, figname="maneuver", tstages=[])`

Plots the kinematics and controls of a `KinematicManeuver`.
"""
function plot_maneuver(maneuver::KinematicManeuver;
                        ti::Real=0, tf::Real=1, vis_nsteps=300,
                        figname="maneuver", tstages=[], size_factor=1.0)

    Vvhcl = maneuver.Vvehicle
    avhcl = maneuver.anglevehicle
    ts = range(ti, tf, length=vis_nsteps)

    # -------------------- Vehicle velocity history ----------------------------
    figure(figname*"-kinematics", figsize=[7*3, 5*1]*size_factor)
    suptitle("VEHICLE KINEMATICS")

    subplot(131)
    title("Velocity")
    Vs = Vvhcl.(ts)
    Vmax = max([maximum([V[i] for V in Vs]) for i in 1:3]...)
    Vmin = min([minimum([V[i] for V in Vs]) for i in 1:3]...)

    plot(ts, [V[1] for V in Vs], "-", label=L"V_x", alpha=0.8, color=clrs[1])
    plot(ts, [V[2] for V in Vs], "-", label=L"V_y", alpha=0.8, color=clrs[2])
    plot(ts, [V[3] for V in Vs], "-", label=L"V_z", alpha=0.8, color=clrs[3])

    for tstage in tstages
        plot(tstage*ones(2), [Vmin, Vmax], ":k", alpha=0.5)
    end

    xlabel("Non-dimensional time")
    ylabel("Non-dimensional velocity")
    legend(loc="best", frameon=false)

    # -------------------- Mission profile -------------------------------------
    Xinit = zeros(3)        # Initial position
    dt = ts[2]-ts[1]        # Time step

    Xs = [Xinit]
    for i in 1:(size(ts,1)-1)
        push!(Xs, Xs[end] + dt*Vs[i])
    end

    subplot(132)
    title("Position")
    Xmax = max([maximum([X[i] for X in Xs]) for i in 1:3]...)
    Xmin = min([minimum([X[i] for X in Xs]) for i in 1:3]...)

    plot(ts, [X[1] for X in Xs], "-g", label=L"x", alpha=0.8, color=clrs[1])
    plot(ts, [X[2] for X in Xs], "-b", label=L"y", alpha=0.8, color=clrs[2])
    plot(ts, [X[3] for X in Xs], "-r", label=L"z", alpha=0.8, color=clrs[3])

    for tstage in tstages
        plot(tstage*ones(2), [Xmin, Xmax], ":k", alpha=0.5)
    end

    xlabel("Non-dimensional time")
    ylabel("Non-dimensional position")
    legend(loc="best", frameon=false)

    # -------------------- Vehicle angle history -------------------------------
    subplot(133)
    title("Angles")
    as = avhcl.(ts)
    amax = max([maximum([a[i] for a in as]) for i in 1:3]...)
    amin = min([minimum([a[i] for a in as]) for i in 1:3]...)

    plot(ts, [a[1] for a in as], "-g", label=L"\theta_x", alpha=0.8, color=clrs[1])
    plot(ts, [a[2] for a in as], "-b", label=L"\theta_y", alpha=0.8, color=clrs[2])
    plot(ts, [a[3] for a in as], "-r", label=L"\theta_z", alpha=0.8, color=clrs[3])

    for tstage in tstages
        plot(tstage*ones(2), [amin, amax], ":k", alpha=0.5)
    end

    xlabel("Non-dimensional time")
    ylabel(L"Angle ($^\circ$)")
    legend(loc="best", frameon=false)

    # -------------------- Tilting system history ------------------------------
    angle_syss = [a.(ts) for a in maneuver.angle]    # Angles of every tilt sys
    RPM_syss = [rpm.(ts) for rpm in maneuver.RPM]    # RPM of every rotor system

    nplots = 3*(length(angle_syss)!=0) + 1*(length(RPM_syss)!=0)
    gdims = nplots==4 ? [2, 2] : nplots==3 ? [2, 2] :
            nplots==1 ? [1, 1] : nplots==0 ? [0, 0] : [2, 2]
    gnum = gdims[1]*100 + gdims[2]*10
    ploti = 1

    if length(angle_syss)!=0 || length(RPM_syss)!=0
        figure(figname*"-controls", figsize=[7, 5].*gdims * size_factor * (nplots==1 ? 3/6 : 2/3))
        suptitle("VEHICLE CONTROLS")
    end

    if length(angle_syss)!=0
        for i in 1:3
            subplot(gnum+ploti)
            ploti += 1

             # i-th angle of every tilting system
            a_syss = [[a[i] for a in angle_sys] for angle_sys in angle_syss]

            amax = max([maximum(a_sys) for a_sys in a_syss]...)
            amin = min([minimum(a_sys) for a_sys in a_syss]...)

            for (j, a_sys) in enumerate(a_syss)
                plot(ts, a_sys, "-", label="Tilt-sys #$j", alpha=0.8,
                                                color=clrs[(j-1)%length(clrs) + 1])
            end

            xlabel("Non-dimensional time")
            ylabel("Angle "*(i==1 ? L"\theta_x" : i==2 ? L"\theta_y" : L"\theta_z") * L" ($^\circ$)")
            legend(loc="best", frameon=false)
        end
    end

    # -------------------- Rotor systems history -------------------------------

    if length(RPM_syss)!=0
        subplot(gnum+ploti)
        ploti += 1

        RPMmax = max([maximum(RPM_sys) for RPM_sys in RPM_syss]...)
        RPMmin = min([minimum(RPM_sys) for RPM_sys in RPM_syss]...)

        for (j, RPM_sys) in enumerate(RPM_syss)
            plot(ts, RPM_sys, "-", label="Rotor-sys #$j", alpha=0.8,
                                            color=clrs[(j-1)%length(clrs) + 1])
        end

        xlabel("Non-dimensional time")
        ylabel("Non-dimensional RPM (RPM/RPMh)")
        legend(loc="best", frameon=false)
    end

end

"""
    `visualize_kinematics(sim::Simulation, nsteps::Int, save_path::String)`

Generate VTKs of the kinematics of this simulation.
"""
function visualize_kinematics(sim::Simulation{V, KinematicManeuver{N, M}, R},
                                nsteps::Int, save_path::String;
                                run_name="vis",
                                prompt=true,
                                verbose=true, v_lvl=0, verbose_nsteps=10,
                                paraview=true,
                                save_vtk_optsargs=[],
                              ) where {V<:AbstractVehicle, R<:Real, N, M}

    dt = sim.ttot/nsteps                # Time step

    # Create save path
    gt.create_path(save_path, prompt)
    strn = ""

    if verbose
        println("\t"^(v_lvl)*"*"^(73-7*v_lvl))
        println("\t"^(v_lvl)*"START $(joinpath(save_path,run_name))")
        println("\t"^(v_lvl)*"*"^(73-7*v_lvl))

        figure(run_name, figsize=[7*2, 5*1])
        subplot(131)
        xlabel("Simulation time")
        ylabel("Velocity")
        Vlbls = [L"V_x", L"V_y", L"V_z"]
        subplot(132)
        xlabel("Simulation time")
        ylabel(L"Angular velocity ($^\circ/t$)")
        Wlbls = [L"\Omega_x", L"\Omega_y", L"\Omega_z"]
        subplot(133)
        xlabel("Simulation time")
        ylabel(L"$O$ position")
        Olbls = [L"O_x", L"O_y", L"O_z"]
    end

    # Time stepping
    for i in 0:nsteps

        # if i!=0
            # Move tilting systems, and translate and rotate vehicle
            nextstep_kinematic(sim, dt)
            rotate_rotors(sim, dt)
        # end

        # Verbose
        if verbose && i%verbose_nsteps==0
            println("\t"^(v_lvl+1)*"Time step $i out of $nsteps")
        end
        if verbose
            for j in 1:3
                subplot(131)
                plot(sim.t, sim.vehicle.V[j], ".", label=Vlbls[j], alpha=0.8,
                                                                color=clrs[j])
                subplot(132)
                plot(sim.t, sim.vehicle.W[j], ".", label=Wlbls[j], alpha=0.8,
                                                                color=clrs[j])
                subplot(133)
                plot(sim.t, sim.vehicle.system.O[j], ".", label=Olbls[j], alpha=0.8,
                                                                color=clrs[j])
            end
            if i==0
                for j in 1:3
                    subplot(130+j)
                    legend(loc="best", frameon=false)
                end
            end
        end

        strn = save_vtk(sim, run_name; path=save_path, save_vtk_optsargs...)
    end

    # Tweak vtk string to be a time sequence
    if nsteps>1
        strn = replace(strn, ".$(nsteps)."=>"...")
    end

    # Call paraview
    if paraview
        if verbose; println("\t"^(v_lvl)*"Calling Paraview..."); end
        run(`paraview --data="$save_path/$strn"`)
    end

    return strn
end

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
indices = 0:49
rates = [2000, 4000, 6000, 8000, 10000, 14000, 18000, 22000, 26000, 30000, 34000, 38000, 42000, 46000, 50000]
max_iter = 2000

fw_arguments = (
    line_search=FrankWolfe.MonotonousNonConvexStepSize(),
    max_iteration=max_iter,
    print_iter=max_iter / 10,
    verbose=true,
    #lazy=true,
)

max_stars_repo_name: jmaces/fw-rde

max_issues_repo_name: ZIB-IOL/fw-rde

max_forks_repo_name: ZIB-IOL/fw-rde


text:
"""
Initialitation of model of 0-state synchronous (classic model) machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.BaseMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data
    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    #Assumption of Classical Machine: Xq = Xd_p
    Xd_p = PSY.get_Xd_p(machine)

    δ0 = angle(V + (R + Xd_p * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * (Vf0 - V_dq[2]) - R * V_dq[1])  #15.36
        i_q = (1.0 / (R^2 + Xd_p^2)) * (Xd_p * V_dq[1] + R * (Vf0 - V_dq[2])) #15.36
        Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q
        out[1] = τm - Pe #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
    end
    x0 = [δ0, τm0, 1.0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Not necessary to update Vf for AVR in Base Machine. Update eq_p:
        PSY.set_eq_p!(machine, sol_x0[3])
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
    end
end

"""
Initialitation of model of 2-state (One d- and One q-) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.OneDOneQMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data
    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)

    #States of OneDOneQMachine are [1] eq_p and [2] ed_p
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0, eq_p, ed_p
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        eq_p = x[4]
        ed_p = x[5]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d = (1.0 / (R^2 + Xd_p * Xq_p)) * (Xq_p * (eq_p - V_dq[2]) + R * (ed_p - V_dq[1]))  #15.32
        i_q =
            (1.0 / (R^2 + Xd_p * Xq_p)) * (-Xd_p * (ed_p - V_dq[1]) + R * (eq_p - V_dq[2]))  #15.32
        Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q
        out[1] = τm - Pe #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
        out[4] = -eq_p - (Xd - Xd_p) * i_d + Vf0 #∂(eq_p)/∂t
        out[5] = -ed_p + (Xq - Xq_p) * i_q #∂(ed_p)/∂t
    end
    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[2], V_dq0[1]]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update eq_p and ed_p for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.OneDOneQMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ed_p
    end
end

"""
Initialitation of model of 6-state (Marconato) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.MarconatoMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    T_AA = PSY.get_T_AA(machine)
    Td0_p = PSY.get_Td0_p(machine)
    γd = PSY.get_γd(machine)
    γq = PSY.get_γq(machine)

    #States of MarconatoMachine are [1] ψq, [2] ψd, [3] eq_p, [4] ed_p, [5] eq_pp and [6] ed_pp
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0, eq_p, ed_p
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        ψq = x[4]
        ψd = x[5]
        eq_p = x[6]
        ed_p = x[7]
        eq_pp = x[8]
        ed_pp = x[9]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d = (1.0 / Xd_pp) * (eq_pp - ψd)      #15.18
        i_q = (1.0 / Xq_pp) * (-ed_pp - ψq)     #15.18
        τ_e = ψd * i_q - ψq * i_d               #15.6
        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
        out[4] = R * i_q - ω0 * ψd + V_dq[2]                                    #15.9 ψq
        out[5] = R * i_d + ω0 * ψq + V_dq[1]                                    #15.9 ψd
        out[6] = -eq_p - (Xd - Xd_p - γd) * i_d + (1 - (T_AA / Td0_p)) * Vf0    #15.16 eq_p
        out[7] = -ed_p + (Xq - Xq_p - γq) * i_q                                 #15.16 ed_p
        out[8] = -eq_pp + eq_p - (Xd_p - Xd_pp + γd) * i_d + (T_AA / Td0_p) * Vf0       #15.16 eq_pp
        out[9] = -ed_pp + ed_p + (Xq_p - Xq_pp + γq) * i_q #15.16 ed_pp
    end

    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[1], V_dq0[2], V_dq0[2], V_dq0[1], V_dq0[2], V_dq0[1]]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.MarconatoMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #ψq
        machine_states[2] = sol_x0[5] #ψd
        machine_states[3] = sol_x0[6] #eq_p
        machine_states[4] = sol_x0[7] #ed_p
        machine_states[5] = sol_x0[8] #eq_pp
        machine_states[6] = sol_x0[9] #ed_pp
        #Update fluxes inner vars
        get_inner_vars(dynamic_device)[ψq_var] = sol_x0[4]
        get_inner_vars(dynamic_device)[ψd_var] = sol_x0[5]
    end
end

"""
Initialitation of model of 4-state (Simple Marconato) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.SimpleMarconatoMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)
    T_AA = PSY.get_T_AA(machine)
    Td0_p = PSY.get_Td0_p(machine)
    γd = PSY.get_γd(machine)
    γq = PSY.get_γq(machine)

    #States of SimpleMarconatoMachine are [1] eq_p, [2] ed_p, [3] eq_pp and [4] ed_pp
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0, eq_p, ed_p
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        eq_p = x[4]
        ed_p = x[5]
        eq_pp = x[6]
        ed_pp = x[7]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (Xq_pp * (eq_pp - V_dq[2]) + R * (ed_pp - V_dq[1]))      #15.25
        i_q =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (-Xd_pp * (ed_pp - V_dq[1]) + R * (eq_pp - V_dq[2]))      #15.25
        Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q
        out[1] = τm - Pe #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
        out[4] = -eq_p - (Xd - Xd_p - γd) * i_d + (1 - (T_AA / Td0_p)) * Vf0             #15.16 eq_p
        out[5] = -ed_p + (Xq - Xq_p - γq) * i_q                                         #15.16 ed_p
        out[6] = -eq_pp + eq_p - (Xd_p - Xd_pp + γd) * i_d + (T_AA / Td0_p) * Vf0        #15.16 eq_pp
        out[7] = -ed_pp + ed_p + (Xq_p - Xq_pp + γq) * i_q                              #15.16 ed_pp
    end
    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[2], V_dq0[1], V_dq0[2], V_dq0[1]]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update eq_p and ed_p for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.SimpleMarconatoMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ed_p
        machine_states[3] = sol_x0[6] #eq_pp
        machine_states[4] = sol_x0[7] #ed_pp
    end
end

"""
Initialitation of model of 6-state (Anderson-Fouad) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.AndersonFouadMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)

    #States of Anderson-Fouad are [1] ψq, [2] ψd, [3] eq_p, [4] ed_p, [5] eq_pp and [6] ed_pp
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0, eq_p, ed_p
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        ψq = x[4]
        ψd = x[5]
        eq_p = x[6]
        ed_p = x[7]
        eq_pp = x[8]
        ed_pp = x[9]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d = (1.0 / Xd_pp) * (eq_pp - ψd)      #15.18
        i_q = (1.0 / Xq_pp) * (-ed_pp - ψq)     #15.18
        τ_e = ψd * i_q - ψq * i_d               #15.6
        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
        out[4] = R * i_q - ω0 * ψd + V_dq[2]                                    #15.9 ψq
        out[5] = R * i_d + ω0 * ψq + V_dq[1]                                    #15.9 ψd
        out[6] = -eq_p + (Xd - Xd_p) * i_d + Vf0           #15.19 eq_p
        out[7] = -ed_p + (Xq - Xq_p) * i_q                 #15.19 ed_p
        out[8] = -eq_pp + eq_p - (Xd_p - Xd_pp) * i_d      #15.19 eq_pp
        out[9] = -ed_pp + ed_p + (Xq_p - Xq_pp) * i_q      #15.19 ed_pp
    end

    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[1], V_dq0[2], V_dq0[2], V_dq0[1], V_dq0[2], V_dq0[1]]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.AndersonFouadMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #ψq
        machine_states[2] = sol_x0[5] #ψd
        machine_states[3] = sol_x0[6] #eq_p
        machine_states[4] = sol_x0[7] #ed_p
        machine_states[5] = sol_x0[8] #eq_pp
        machine_states[6] = sol_x0[9] #ed_pp
        #Update fluxes inner vars
        get_inner_vars(dynamic_device)[ψq_var] = sol_x0[4]
        get_inner_vars(dynamic_device)[ψd_var] = sol_x0[5]
    end
end

"""
Initialitation of model of 4-state (Simple Anderson-Fouad) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.SimpleAFMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = PSY.get_Xq_pp(machine)

    #States of SimpleAndersonFouad are [1] eq_p, [2] ed_p, [3] eq_pp and [4] ed_pp
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))
    #To solve: δ, τm, Vf0, eq_p, ed_p
    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        eq_p = x[4]
        ed_p = x[5]
        eq_pp = x[6]
        ed_pp = x[7]

        V_dq = ri_dq(δ) * [V_R; V_I]
        i_d =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (Xq_pp * (eq_pp - V_dq[2]) + R * (ed_pp - V_dq[1]))      #15.25
        i_q =
            (1.0 / (R^2 + Xd_pp * Xq_pp)) *
            (-Xd_pp * (ed_pp - V_dq[1]) + R * (eq_pp - V_dq[2]))      #15.25
        Pe = (V_dq[1] + R * i_d) * i_d + (V_dq[2] + R * i_q) * i_q
        out[1] = τm - Pe #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_d + V_dq[2] * i_q) #Output Power
        out[3] = Q0 - (V_dq[2] * i_d - V_dq[1] * i_q) #Output Reactive Power
        out[4] = -eq_p + (Xd - Xd_p) * i_d + Vf0           #15.19 eq_p
        out[5] = -ed_p + (Xq - Xq_p) * i_q                 #15.19 ed_p
        out[6] = -eq_pp + eq_p - (Xd_p - Xd_pp) * i_d      #15.19 eq_pp
        out[7] = -ed_pp + ed_p + (Xq_p - Xq_pp) * i_q      #15.19 ed_pp
    end
    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[2], V_dq0[1], V_dq0[2], V_dq0[1]]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update eq_p and ed_p for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.SimpleAFMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ed_p
        machine_states[3] = sol_x0[6] #eq_pp
        machine_states[4] = sol_x0[7] #ed_pp
    end
end

function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{M, S, A, TG, P},
) where {
    M <: Union{PSY.RoundRotorQuadratic, PSY.RoundRotorExponential},
    S <: PSY.Shaft,
    A <: PSY.AVR,
    TG <: PSY.TurbineGov,
    P <: PSY.PSS,
}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_p = PSY.get_Tq0_p(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xq_p = PSY.get_Xq_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    Se = PSY.get_Se(machine)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)
    γ_q2 = PSY.get_γ_q2(machine)
    γ_qd = PSY.get_γ_qd(machine)

    ## Initialization ##
    ## Fluxes
    ψ_pp0 = V + (R + Xd_pp * 1im) * I
    ψ_pp0_ang = angle(ψ_pp0)
    ψ_pp0_abs = abs(ψ_pp0)
    Se0 = saturation_function(machine, ψ_pp0_abs)
    ## Angles
    _a = ψ_pp0_abs * (Se0 * γ_qd + 1.0)
    _b = (Xq_pp - Xq) * abs(I)
    θ_It = ψ_pp0_ang - angle(I)
    δ0 = ψ_pp0_ang + atan((_b * cos(θ_It)) / (_b * sin(θ_It) - _a))
    #Fluxes
    ψ_pp0_dq = ri_dq(δ0 + pi / 2) * [real(ψ_pp0); imag(ψ_pp0)]
    ψd_pp0 = ψ_pp0_dq[1]
    ψq_pp0 = -ψ_pp0_dq[2]
    ## Currents
    I_d0, I_q0 = ri_dq(δ0) * [real(I); imag(I)]
    ## Voltages and Fluxes
    V_d0, V_q0 = ri_dq(δ0) * [V_R; V_I]
    ψq_pp0 = V_d0 - R * I_d0 - Xq_pp * I_q0
    ψd_pp0 = V_q0 + I_d0 * Xd_pp + I_q0 * R
    # States
    ## External Variables
    τm0 = I_d0 * (V_d0 + I_d0 * R) + I_q0 * (V_q0 + I_q0 * R)
    Vf0 = I_d0 * (Xd - Xd_pp) + ψd_pp0 * (Se0 + 1.0)
    ψd0 = V_q0 + R * I_q0
    ψq0 = -V_d0 - R * I_d0
    ## States
    eq_p0 = I_d0 * (Xd_p - Xd) - Se0 * ψd_pp0 + Vf0
    ed_p0 = I_q0 * (Xq - Xq_p) - Se0 * γ_qd * ψq_pp0
    ψ_kd0 = -I_d0 * (Xd - Xl) - Se0 * ψd_pp0 + Vf0
    ψ_kq0 = I_q0 * (Xq - Xl) - Se0 * γ_qd * ψq_pp0
    #Auxiliar Variables
    Xad_Ifd0 =
        ed_p0 +
        (Xq - Xq_p) * (γ_q2 * ed_p0 - γ_q2 * ψ_kq0 - γ_q1 * I_q0) +
        Se0 * ψq_pp0 * γ_qd

    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf = x[3]
        eq_p = x[4]
        ed_p = x[5]
        ψ_kd = x[6]
        ψ_kq = x[7]
        Xad_Ifd_aux = x[8]

        V_dq = ri_dq(δ) * [V_R; V_I]
        ψq_pp = γ_q1 * ed_p + ψ_kq * (1 - γ_q1)
        ψd_pp = γ_d1 * eq_p + γ_d2 * (Xd_p - Xl) * ψ_kd
        ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)
        I_d =
            (1.0 / (R^2 + Xq_pp * Xd_pp)) *
            (-R * (V_dq[1] - ψq_pp) + Xq_pp * (-V_dq[2] + ψd_pp))
        I_q =
            (1.0 / (R^2 + Xq_pp * Xd_pp)) *
            (Xd_pp * (V_dq[1] - ψq_pp) + R * (-V_dq[2] + ψd_pp))
        Se = saturation_function(machine, ψ_pp)
        Xad_Ifd = eq_p + (Xd - Xd_p) * (γ_d1 * I_d - γ_d2 * ψ_kd + γ_d2 * eq_p) + Se * ψd_pp
        Xaq_I1q =
            ed_p +
            (Xq - Xq_p) * (γ_q2 * ed_p - γ_q2 * ψ_kq - γ_q1 * I_q) +
            Se * ψq_pp * γ_qd
        τ_e = I_d * (V_dq[1] + I_d * R) + I_q * (V_dq[2] + I_q * R)

        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_dq[1] * I_d + V_dq[2] * I_q) #Output Power
        out[3] = Q0 - (V_dq[2] * I_d - V_dq[1] * I_q) #Output Reactive Power
        out[4] = (1.0 / Td0_p) * (Vf - Xad_Ifd) #deq_p/dt
        out[5] = (1.0 / Tq0_p) * (-Xaq_I1q) #ded_p/dt
        out[6] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d) #dψ_kd/dt
        out[7] = (1.0 / Tq0_pp) * (-ψ_kq + ed_p + (Xq_p - Xl) * I_q) #deq_pp/dt
        out[8] = Xad_Ifd_aux - Xad_Ifd
    end
    x0 = [δ0, τm0, Vf0, eq_p0, ed_p0, ψ_kd0, ψ_kq0, Xad_Ifd0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = 1.0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in GENROU Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update Xad_Ifd for AVR in GENROU Machine
        get_inner_vars(dynamic_device)[Xad_Ifd_var] = sol_x0[8]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, typeof(machine))
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ed_p
        machine_states[3] = sol_x0[6] #ψ_kd
        machine_states[4] = sol_x0[7] #ψ_kq
    end
end

function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.SalientPoleQuadratic, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)

    ## Initialization ##
    E = V + (R + Xq * 1im) * I
    δ0 = angle(E)
    ## Currents
    I_d0, I_q0 = ri_dq(δ0) * [real(I); imag(I)]
    ## Voltages and fluxes
    V_d0, V_q0 = ri_dq(δ0) * [V_R; V_I]
    ψ_d0 = V_q0 + R * I_q0
    ψ_q0 = -V_d0 - R * I_d0
    ψd_pp0 = V_q0 + I_d0 * Xd_pp + I_q0 * R
    # States
    ψq_pp0 = -(V_d0 - I_q0 * Xq_pp - I_d0 * R)
    eq_p0 = (1 / (γ_d1 + γ_q1)) * (ψd_pp0 + I_d0 * (Xd_p - Xl) * γ_q1)
    ψ_kd0 = eq_p0 - I_d0 * (Xd_p - Xl)
    ## External Variables
    τm0 = ψ_d0 * I_q0 - ψ_q0 * I_d0
    Se0 = saturation_function(machine, eq_p0)
    Vf0 = eq_p0 + I_d0 * (Xd - Xd_p) + Se0 * eq_p0
    Xad_Ifd0 =
        eq_p0 +
        Se0 * eq_p0 +
        (Xd - Xd_p) * (I_d0 + γ_d2 * (eq_p0 - ψ_kd0 - (Xd_p - Xl) * I_d0))

    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf = x[3]
        eq_p = x[4]
        ψ_kd = x[5]
        ψq_pp = x[6]
        Xad_Ifd_aux = x[7]

        V_d, V_q = ri_dq(δ) * [V_R; V_I]
        #Additional Fluxes
        ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd

        #Currents
        I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d + ψq_pp) + Xd_pp * (ψd_pp - V_q))
        I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d + ψq_pp) + R * (ψd_pp - V_q))
        Se = saturation_function(machine, eq_p)
        Xad_Ifd =
            eq_p +
            Se * eq_p +
            (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
        τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)

        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_d * I_d + V_q * I_q) #Output Power
        out[3] = Q0 - (V_q * I_d - V_d * I_q) #Output Reactive Power
        out[4] = (1.0 / Td0_p) * (Vf - Xad_Ifd) #deq_p/dt
        out[5] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d) #dψ_kd/dt
        out[6] = (1.0 / Tq0_pp) * (-ψq_pp - (Xq - Xq_pp) * I_q) #dψq_pp/dt
        out[7] = Xad_Ifd_aux - Xad_Ifd
    end
    x0 = [δ0, τm0, Vf0, eq_p0, ψ_kd0, ψq_pp0, Xad_Ifd0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = 1.0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in GENSAL Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update Xad_Ifd for AVR in GENSAL Machine
        get_inner_vars(dynamic_device)[Xad_Ifd_var] = sol_x0[7]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, typeof(machine))
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ψ_kd
        machine_states[3] = sol_x0[6] #ψq_pp
    end
end

function initialize_mach_shaft!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.SalientPoleExponential, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}

    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Get parameters
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    Td0_p = PSY.get_Td0_p(machine)
    Td0_pp = PSY.get_Td0_pp(machine)
    Tq0_pp = PSY.get_Tq0_pp(machine)
    Xd = PSY.get_Xd(machine)
    Xq = PSY.get_Xq(machine)
    Xd_p = PSY.get_Xd_p(machine)
    Xd_pp = PSY.get_Xd_pp(machine)
    Xq_pp = Xd_pp
    Xl = PSY.get_Xl(machine)
    γ_d1 = PSY.get_γ_d1(machine)
    γ_q1 = PSY.get_γ_q1(machine)
    γ_d2 = PSY.get_γ_d2(machine)
    γ_qd = (Xq - Xl) / (Xd - Xl)

    ## Initialization ##
    ## Fluxes
    ψ_pp0 = V + (R + Xd_pp * 1im) * I
    ψ_pp0_ang = angle(ψ_pp0)
    ψ_pp0_abs = abs(ψ_pp0)
    Se0 = saturation_function(machine, ψ_pp0_abs)
    ## Angles
    _a = ψ_pp0_abs * (Se0 * γ_qd + 1.0)
    _b = (Xq_pp - Xq) * abs(I)
    θ_It = ψ_pp0_ang - angle(I)
    δ0 = ψ_pp0_ang + atan((_b * cos(θ_It)) / (_b * sin(θ_It) - _a))
    #Fluxes
    ψ_pp0_dq = ri_dq(δ0 + pi / 2) * [real(ψ_pp0); imag(ψ_pp0)]
    ψd_pp0 = ψ_pp0_dq[1]
    ψq_pp0 = -ψ_pp0_dq[2]
    ## Currents
    I_d0, I_q0 = ri_dq(δ0) * [real(I); imag(I)]
    ## Voltages and Fluxes
    V_d0, V_q0 = ri_dq(δ0) * [V_R; V_I]
    ψq_pp0 = V_d0 - R * I_d0 - Xq_pp * I_q0
    ψd_pp0 = V_q0 + I_d0 * Xd_pp + I_q0 * R
    ψ_d0 = V_q0 + R * I_q0
    ψ_q0 = -V_d0 - R * I_d0
    # States
    eq_p0 = (1 / (γ_d1 + γ_q1)) * (ψd_pp0 + I_d0 * (Xd_p - Xl) * γ_q1)
    ψ_kd0 = eq_p0 - I_d0 * (Xd_p - Xl)
    ## External Variables
    τm0 = ψ_d0 * I_q0 - ψ_q0 * I_d0
    Vf0 = eq_p0 + I_d0 * (Xd - Xd_p) + ψd_pp0 * Se0
    Xad_Ifd0 =
        eq_p0 +
        Se0 * ψd_pp0 +
        (Xd - Xd_p) * (I_d0 + γ_d2 * (eq_p0 - ψ_kd0 - (Xd_p - Xl) * I_d0))

    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf = x[3]
        eq_p = x[4]
        ψ_kd = x[5]
        ψq_pp = x[6]
        Xad_Ifd_aux = x[7]

        V_d, V_q = ri_dq(δ) * [V_R; V_I]
        #Additional Fluxes
        ψd_pp = γ_d1 * eq_p + γ_q1 * ψ_kd
        ψ_pp = sqrt(ψd_pp^2 + ψq_pp^2)

        I_d = (1.0 / (R^2 + Xd_pp^2)) * (-R * (V_d - ψq_pp) + Xq_pp * (-V_q + ψd_pp))
        I_q = (1.0 / (R^2 + Xd_pp^2)) * (Xd_pp * (V_d - ψq_pp) + R * (-V_q + ψd_pp))
        Se = saturation_function(machine, ψ_pp)
        Xad_Ifd =
            eq_p +
            Se * ψd_pp +
            (Xd - Xd_p) * (I_d + γ_d2 * (eq_p - ψ_kd - (Xd_p - Xl) * I_d))
        τ_e = I_d * (V_d + I_d * R) + I_q * (V_q + I_q * R)

        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_d * I_d + V_q * I_q) #Output Power
        out[3] = Q0 - (V_q * I_d - V_d * I_q) #Output Reactive Power
        out[4] = (1.0 / Td0_p) * (Vf - Xad_Ifd) #deq_p/dt
        out[5] = (1.0 / Td0_pp) * (-ψ_kd + eq_p - (Xd_p - Xl) * I_d) #dψ_kd/dt
        out[6] = (1.0 / Tq0_pp) * (-ψq_pp + (Xq - Xq_pp) * I_q - Se * γ_qd * ψq_pp) #dψq_pp/dt
        out[7] = Xad_Ifd_aux - Xad_Ifd
    end
    x0 = [δ0, τm0, Vf0, eq_p0, ψ_kd0, ψq_pp0, Xad_Ifd0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = 1.0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in GENSAL Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update Xad_Ifd for AVR in GENSAL Machine
        get_inner_vars(dynamic_device)[Xad_Ifd_var] = sol_x0[7]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, typeof(machine))
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #eq_p
        machine_states[2] = sol_x0[5] #ψ_kd
        machine_states[3] = sol_x0[6] #ψq_pp
    end
end

#=
"""
Initialitation of model of 5-state (Kundur) synchronous machine in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""
function initialize_mach_shaft!(
    device_states,
static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicGenerator{PSY.FullMachine, S, A, TG, P},
) where {S <: PSY.Shaft, A <: PSY.AVR, TG <: PSY.TurbineGov, P <: PSY.PSS}
    #PowerFlow Data

    P0 = PSY.get_active_power(static)
    Q0 = PSY.get_reactive_power(static)
    Vm = PSY.get_magnitude(PSY.get_bus(static))
    θ = PSY.get_angle(PSY.get_bus(static))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)

    #Machine Data
    machine = PSY.get_machine(dynamic_device)
    R = PSY.get_R(machine)
    R_f = PSY.get_R_f(machine)
    R_1d = PSY.get_R_1d(machine)
    R_1q = PSY.get_R_1q(machine)
    Xq = PSY.get_L_q(machine)'
    Xd = PSY.get_L_d(machine)
    Xad = PSY.get_L_ad(machine)
    inv_d_fluxlink = PSY.get_inv_d_fluxlink(machine)
    inv_q_fluxlink = PSY.get_inv_q_fluxlink(machine)

    #States of FullMachine are [1] ψd, [2] ψq, [3] ψf , [4] ψ1d, [5] ψ1q
    δ0 = angle(V + (R + Xq * 1im) * I)
    ω0 = 1.0
    τm0 = real(V * conj(I))

    function f!(out, x)
        δ = x[1]
        τm = x[2]
        Vf0 = x[3]
        ψd = x[4]
        ψq = x[5]
        ψf = x[6]
        ψ1d = x[7]
        ψ1q = x[8]

        V_dq = ri_dq(δ) * [V_R; V_I]
        #Obtain electric currents] using Flux Linkage equations
        i_df1d = inv_d_fluxlink * [ψd; ψf; ψ1d]  #11.18 in Machowski (3.127, 3.130 and 3.131 in Kundur)
        i_q1q = inv_q_fluxlink * [ψq; ψ1q]        #11.19 in Machowski (3.128 and 3.132 in Kundur)
        τ_e = ψd * i_q1q[1] - ψq * i_df1d[1]            #15.6 in Milano or 3.117 in Kundur
        out[1] = τm - τ_e #Mechanical Torque
        out[2] = P0 - (V_dq[1] * i_df1d[1] + V_dq[2] * i_q1q[1]) #Output Power
        out[3] = Q0 - (V_dq[2] * i_df1d[1] - V_dq[1] * i_q1q[1]) #Output Reactive Power
        out[4] = R * i_df1d[1] + ω0 * ψq + V_dq[1]                                    #15.9 ψd
        out[5] = R * i_q1q[1] - ω0 * ψd + V_dq[2]                                    #15.9 ψq
        out[6] = -R_f * i_df1d[2] + Vf0          #15.19 eq_p
        out[7] = -R_1d * i_df1d[3]                #15.19 ed_p
        out[8] = -R_1q * i_q1q[2]        #15.19 eq_pp
    end

    V_dq0 = ri_dq(δ0) * [V_R; V_I]
    x0 = [δ0, τm0, 1.0, V_dq0[2], -V_dq0[1], 1.0, 0.0, 0.0]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Synch. Machine failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        get_inner_vars(dynamic_device)[VR_gen_var] = V_R
        get_inner_vars(dynamic_device)[VI_gen_var] = V_I
        #Update δ and ω of Shaft. Works for every Shaft.
        shaft_ix = get_local_state_ix(dynamic_device, S)
        shaft_states = @view device_states[shaft_ix]
        shaft_states[1] = sol_x0[1] #δ
        shaft_states[2] = ω0 #ω
        #Update Mechanical and Electrical Torque on Generator
        get_inner_vars(dynamic_device)[τe_var] = sol_x0[2]
        get_inner_vars(dynamic_device)[τm_var] = sol_x0[2]
        #Update Vf for AVR in OneDOneQ Machine.
        get_inner_vars(dynamic_device)[Vf_var] = sol_x0[3]
        #Update states for Machine
        machine_ix = get_local_state_ix(dynamic_device, PSY.FullMachine)
        machine_states = @view device_states[machine_ix]
        machine_states[1] = sol_x0[4] #ψd
        machine_states[2] = sol_x0[5] #ψq
        machine_states[3] = sol_x0[6] #ψf
        machine_states[4] = sol_x0[7] #ψ1d
        machine_states[5] = sol_x0[8] #ψ1q
        #Update fluxes inner vars
        get_inner_vars(dynamic_device)[ψd_var] = sol_x0[4]
        get_inner_vars(dynamic_device)[ψq_var] = sol_x0[5]
    end
end
=#


text:
using DataAndPlotting
using Test
using HDF5
using Plots

filename = "array.h5"

arr = initialize()
rm(filename, force = true)
output(arr, filename)

arr2 = input(filename)

@test isequal(arr2, arr)

plotname = "wave.png"
arr = initialize()
makeplot(arr, plotname)

max_stars_repo_name: brodriguezmilla/DataAndPlotting

max_issues_repo_name: brodriguezmilla/DataAndPlotting

max_forks_repo_name: brodriguezmilla/DataAndPlotting


text:
mutable struct EqualityConstraint{T,N,Nc,Cs} <: AbstractConstraint{T,N}
    id::Int64
    name::String

    constraints::Cs
    pid::Union{Int64,Nothing}
    bodyids::SVector{Nc,Int64}

    s0::SVector{N,T}
    s1::SVector{N,T}

    function EqualityConstraint(data...; name::String="")
        jointdata = Tuple{Joint,Int64,Int64}[]
        for info in data
            if typeof(info[1]) <: Joint
                push!(jointdata, info)
            else
                for subinfo in info
                    push!(jointdata, subinfo)
                end
            end
        end

        T = getT(jointdata[1][1])# .T

        pid = jointdata[1][2]
        bodyids = Int64[]
        constraints = Joint{T}[]
        N = 0
        for set in jointdata
            push!(constraints, set[1])
            @assert set[2] == pid
            push!(bodyids, set[3])
            N += getNc(set[1])
        end
        constraints = Tuple(constraints)
        Nc = length(constraints)

        s0 = @SVector zeros(T, N)
        s1 = @SVector zeros(T, N)

        new{T,N,Nc,typeof(constraints)}(getGlobalID(), name, constraints, pid, bodyids, s0, s1)
    end
end

Base.length(::EqualityConstraint{T,N}) where {T,N} = N

function setForce!(mechanism, eqc::EqualityConstraint{T,N,Nc}, Fτ; K = mechanism.No) where {T,N,Nc}
    for i = 1:Nc
        setForce!(eqc.constraints[i], getbody(mechanism, eqc.pid), getbody(mechanism, eqc.bodyids[i]), Fτ[i], K)
    end
end

function setVelocity!(mechanism, eqc::EqualityConstraint{T,N,Nc}, vω) where {T,N,Nc}
    # TODO currently assumed constraints are in order which is the case unless very low level constraint setting and only joints
    n = Int64(Nc/2)
    body1 = getbody(mechanism, eqc.pid)
    for i = 1:n
        body2 = getbody(mechanism, eqc.bodyids[i])
        Δv = getVelocityDelta(eqc.constraints[i], body1, body2, vω[i])
        Δω = getVelocityDelta(eqc.constraints[i+1], body1, body2, vω[i+1])

        p1, p2 = eqc.constraints[i].vertices
        setVelocity!(mechanism, body1, body2; p1 = p1, p2 = p2, Δv = Δv, Δω = Δω)
    end
end

function setPosition!(mechanism, eqc::EqualityConstraint{T,N,Nc}, xθ) where {T,N,Nc}
    # TODO currently assumed constraints are in order which is the case unless very low level constraint setting and only joints
    n = Int64(Nc/2)
    body1 = getbody(mechanism, eqc.pid)
    for i = 1:n
        body2 = getbody(mechanism, eqc.bodyids[i])
        Δx = getPositionDelta(eqc.constraints[i], body1, body2, xθ[i])
        Δq = getPositionDelta(eqc.constraints[i+1], body1, body2, xθ[i+1])

        p1, p2 = eqc.constraints[i].vertices
        setPosition!(mechanism, body1, body2; p1 = p1, p2 = p2, Δx = Δx, Δq = Δq)
    end
end

@generated function minimalCoordinates(mechanism, eqc::EqualityConstraint{T,N,Nc}; K = mechanism.No) where {T,N,Nc}
    vec = [:(minimalCoordinates(eqc.constraints[$i], getbody(mechanism, eqc.pid), getbody(mechanism, eqc.bodyids[$i]), K)) for i = 1:Nc]
    :(svcat($(vec...)))
end

@generated function g(mechanism, eqc::EqualityConstraint{T,N,Nc}) where {T,N,Nc}
    vec = [:(g(eqc.constraints[$i], getbody(mechanism, eqc.pid), getbody(mechanism, eqc.bodyids[$i]), mechanism.Δt, mechanism.No)) for i = 1:Nc]
    :(svcat($(vec...)))
end

@inline function ∂g∂pos(mechanism, eqc::EqualityConstraint, id::Int64)
    id == eqc.pid ? ∂g∂posa(mechanism, eqc, id) : ∂g∂posb(mechanism, eqc, id)
end

@inline function ∂g∂vel(mechanism, eqc::EqualityConstraint, id::Int64)
    id == eqc.pid ? ∂g∂vela(mechanism, eqc, id) : ∂g∂velb(mechanism, eqc, id)
end

@generated function ∂g∂posa(mechanism, eqc::EqualityConstraint{T,N,Nc}, id::Int64) where {T,N,Nc}
    vec = [:(∂g∂posa(eqc.constraints[$i], getbody(mechanism, id), getbody(mechanism, eqc.bodyids[$i]), mechanism.No)) for i = 1:Nc]
    return :(vcat($(vec...)))
end

@generated function ∂g∂posb(mechanism, eqc::EqualityConstraint{T,N,Nc}, id::Int64) where {T,N,Nc}
    vec = [:(∂g∂posb(eqc.constraints[$i], getbody(mechanism, eqc.pid), getbody(mechanism, id), mechanism.No)) for i = 1:Nc]
    return :(vcat($(vec...)))
end

@generated function ∂g∂vela(mechanism, eqc::EqualityConstraint{T,N,Nc}, id::Int64) where {T,N,Nc}
    vec = [:(∂g∂vela(eqc.constraints[$i], getbody(mechanism, id), getbody(mechanism, eqc.bodyids[$i]), mechanism.Δt, mechanism.No)) for i = 1:Nc]
    return :(vcat($(vec...)))
end

@generated function ∂g∂velb(mechanism, eqc::EqualityConstraint{T,N,Nc}, id::Int64) where {T,N,Nc}
    vec = [:(∂g∂velb(eqc.constraints[$i], getbody(mechanism, eqc.pid), getbody(mechanism, id), mechanism.Δt, mechanism.No)) for i = 1:Nc]
    return :(vcat($(vec...)))
end


text:
using JSON, BinaryBuilder, Test

@testset "Meta JSON" begin
    mktempdir() do tmpdir
        meta_json_buff = IOBuffer()
        # Run autobuild() a few times to generate a moderately complex `meta.json`:
        build_output_meta = autobuild(
            tmpdir,
            "libfoo",
            v"1.0.0",
            [("https://julialang.org" => "123123"), "./bundled"],
            "exit 1",
            [Linux(:x86_64)],
            Product[ExecutableProduct("libfoo", :libfoo)],
            ["Zlib_jll"];
            meta_json_stream=meta_json_buff,
        )
        @test build_output_meta == Dict()

        build_output_meta = autobuild(
            tmpdir,
            "libfoo",
            v"1.0.0",
            [("https://github.com/JuliaLang/julia.git" => "5d4eaca0c9fa3d555c79dbacdccb9169fdf64b65")],
            "exit 0",
            [Linux(:x86_64), Windows(:x86_64)],
            Product[ExecutableProduct("julia", :julia)],
            [];
            meta_json_stream=meta_json_buff,
        )
        @test build_output_meta == Dict()

        # Now, deserialize the info:
        seek(meta_json_buff, 0)

        # Strip out ending newlines as that makes our while loop below sad
        meta_json_buff = IOBuffer(strip(String(take!(meta_json_buff))))
        objs = []
        while !eof(meta_json_buff)
            push!(objs, JSON.parse(meta_json_buff))
        end

        # Ensure that we get two JSON objects
        @test length(objs) == 2

        # Merge them, then test that the merged object contains everything we expect
        meta = BinaryBuilder.cleanup_merged_object!(BinaryBuilder.merge_json_objects(objs))

        @test all(haskey.(Ref(meta), ("name", "version", "script", "platforms", "products", "dependencies")))
        @test meta["name"] == "libfoo"
        @test meta["version"] == v"1.0.0"
        @test length(meta["sources"]) == 3
        @test all(in.(
              (
                  "https://julialang.org" => "123123",
                  "https://github.com/JuliaLang/julia.git" => "5d4eaca0c9fa3d555c79dbacdccb9169fdf64b65",
                  "./bundled",
              ), Ref(meta["sources"])))
        @test length(meta["products"]) == 2
        @test all(in.(("libfoo", "julia"), Ref(meta["products"])))
        @test length(meta["script"]) == 2
        @test all(in.(("exit 0", "exit 1"), Ref(meta["script"])))
    end
end


text:
"""
$(TYPEDEF)
"""
mutable struct SLAMWrapper
  fg::IncrementalInference.FactorGraph
  tree::NothingUnion{IncrementalInference.BayesTree}
  lndmidx::Int
  lastposesym::Symbol
  lastpose::SE3
  lbl2usrid::Dict{Symbol,Int}
  usrid2lbl::Dict{Int, Symbol}
  SLAMWrapper() = new()
  SLAMWrapper(a,b,c) = new(a,b,c,:x1, SE3(0), Dict{Symbol,Int}(), Dict{Int,Symbol}()) # TODO to be deprecated
end

function addOdoFG!(slaml::SLAMWrapper, odo::Pose3Pose3;
                  N::Int=100, ready::Int=1,
                  saveusrid::Int=-1)
  #
  vprev = getVert(slaml.fg, slaml.lastposesym)
  # vprev, X, nextn = getLastPose(fgl)
  npnum = parse(Int,string(slaml.lastposesym)[2:end]) + 1
  nextn = Symbol("x$(npnum)")
  vnext = addVariable!(slaml.fg, nextn, Pose2(labels=["POSE";]), N=N, ready=ready)
  # vnext = addVariable!(slaml.fg, nextn, getVal(vprev) ⊕ odo, N=N, ready=ready, labels=["POSE"])
  slaml.lastposesym = nextn
  fact = addFactor!(slaml.fg, [vprev;vnext], odo)

  if saveusrid > -1
    slaml.lbl2usrid[nextn] = saveusrid
    slaml.usrid2lbl[saveusrid] = nextn
  end
  return vnext, fact
end

function addposeFG!(slaml::SLAMWrapper,
      constrs::Vector{IncrementalInference.FunctorInferenceType};
      N::Int=100,
      ready::Int=1,
      saveusrid::Int=-1   )
  #
  vprev = getVert(slaml.fg, slaml.lastposesym)

  npnum = parse(Int,string(slaml.lastposesym)[2:end]) + 1
  nextn = Symbol("x$(npnum)")
  # preinit
  vnext = nothing
  if !haskey(slaml.fg.IDs, nextn)
    vnext = addVariable!(slaml.fg, nextn, Pose2, N=N, ready=ready)
    # vnext = addVariable!(slaml.fg, nextn, getVal(vprev), N=N, ready=ready, labels=["POSE"])
  else
    vnext = getVert(slaml.fg, nextn) #, api=localapi # as optimization, assuming we already have latest vnest in slaml.fg
  end
  slaml.lastposesym = nextn

  addsubtype(fgl::FactorGraph, vprev, vnext, cc::IncrementalInference.FunctorPairwise) = addFactor!(fgl, [vprev;vnext], cc)
  addsubtype(fgl::FactorGraph, vprev, vnext, cc::IncrementalInference.FunctorSingleton) = addFactor!(fgl, [vnext], cc)

  facts = Graphs.ExVertex[]
  PP = BallTreeDensity[]
  for cns in constrs
    fa = addsubtype(slaml.fg, vprev, vnext, cns)
    push!(facts, fa)
  end

  # set node val from new constraints as init
  val, = predictbelief(slaml.fg, vnext, facts, N=N)
  setVal!(vnext, val)
  IncrementalInference.dlapi.updatevertex!(slaml.fg, vnext)

  if saveusrid > -1
    slaml.lbl2usrid[nextn] = saveusrid
    slaml.usrid2lbl[saveusrid] = nextn
  end
  return vnext, facts
end

#


text:
include("test_A_mul_B.jl")
include("test_Ac_mul_B.jl")


text:
module DataType
using MLStyle
using MLStyle.Toolz: isCapitalized, ($), cons, nil
using MLStyle.MatchCore
using MLStyle.Qualification
using MLStyle.Infras
using MLStyle.Pervasives
using MLStyle.Render: render
using MLStyle.Record: def_record
using MLStyle.TypeVarExtraction

export @data
is_symbol_capitalized = isCapitalized ∘ string

"""
Codegen for `@data`, e.g.,

```
@data A begin              abstract type A end
    A1(Int, Int)     ≡     struct A1
                               _1 :: Int
end                            _2 :: Int
                           end

@data B{T}<:A begin        abstract type B{T} <: A end
    B1(T, Int)        ≡    struct B1{T}
end                            _1 :: T
                               _2 :: Int
                           end

@data C{T}<:B{T} begin          abstract type C{T} <: B{T} end
    C1{A, B} ::
    (A, B)=>C{Tuple{A, B}} ≡    struct C1{A, B} <: C{Tuple{A, B}}
end                                 _1 :: A
                                    _2 :: B
                                end

@data D{T}<:A begin        abstract type D{T} <: A end
    D1(a::T)          ≡    struct D1{T} <: D{T}
end                            a :: T
                           end
```

Additionally, `@data` macros setup up pattern matching for those datatypes.
```julia

@data S begin
    S1(Int)
end

@match S1(1) begin
    S1(x) => x == 1
end # => true
```
"""
macro data(typ, def_variants)
    data(typ, def_variants, :public, __module__) |> esc
end

macro data(qualifier, typ, def_variants)
    data(typ, def_variants, qualifier, __module__) |> esc
end

function get_tvars(t :: UnionAll)
    cons(t.var.name, get_tvars(t.body))
end

function get_tvars(t :: Base.DataType)
   nil()
end

function get_tvars(t :: Union)
   nil()
end

function data(typ, def_variants, qualifier_ast, mod)
    typename =
        @match typ begin
           :($typename{$(a...)})       => typename
           :($typename{$(a...)} <: $b) => typename
           :($typename <: $b)          => typename
           :($(typename :: Symbol))    => typename
        end

    mod.eval(:(abstract type $typ end))

    original_ty = getfield(mod, typename)
    tvars_of_abst = get_tvars(original_ty)

    for (ctor_name, pairs, each) in impl(original_ty, def_variants, mod)
        mod.eval(each)
        ctor = getfield(mod, ctor_name)
        qualifier = get_qualifier(qualifier_ast, mod)
        def_record(ctor, convert(Vector{Symbol}, map(first, pairs)), qualifier, mod)
    end
    nothing
end

function impl(t, variants :: Expr, mod :: Module)
    l :: LineNumberNode = LineNumberNode(1)
    abs_tvars = collect(get_tvars(t))
    defs = []
    abst() = isempty(abs_tvars) ? t : :($t{$(abs_tvars...)})
    VAR = mangle(mod)

    for each in variants.args
        @match each begin
            ::LineNumberNode => (l = each)
            :($case{$(generic_tvars...)} :: ($(params...), ) => $(ret_ty) where {$(implicit_tvars...)})                          ||
            :($case{$(generic_tvars...)} :: ($(params...), ) => $(ret_ty && Do(implicit_tvars=[])))                              ||
            :($case{$(generic_tvars...)} :: $(arg_ty && Do(params = [arg_ty])) => $ret_ty where {$(implicit_tvars...)})          ||
            :($case{$(generic_tvars...)} :: $(arg_ty && Do(params = [arg_ty])) => $(ret_ty && Do(implicit_tvars=[])))            ||
            :($(case && Do(generic_tvars = [])) :: ($(params...), ) => $(ret_ty) where {$(implicit_tvars...)})                   ||
            :($(case && Do(generic_tvars = [])) :: ($(params...), ) => $(ret_ty && Do(implicit_tvars=[])))                       ||
            :($(case && Do(generic_tvars = [])) :: $(arg_ty && Do(params = [arg_ty])) => $(ret_ty) where {$(implicit_tvars...)}) ||
            :($(case && Do(generic_tvars = [])) :: $(arg_ty && Do(params = [arg_ty])) => $(ret_ty && Do(implicit_tvars=[])))     ||

            :($case($((params && Do(generic_tvars=abs_tvars))...))) && Do(ret_ty = abst(), implicit_tvars=[]) => begin

              config = Dict{Symbol, Any}([(gtvar isa Expr ? gtvar.args[1] : gtvar) => Any for gtvar in implicit_tvars])

              pairs = map(enumerate(params)) do (i, each)
                 @match each begin
                    :($(a::Symbol && function (x) !is_symbol_capitalized(x) end)) => (a, Any, Any)
                    :($(a && function is_symbol_capitalized  end)) => (Symbol("_$i"), a, render(a, config))
                    :($field :: $ty)                => (field, ty, render(ty, config))
                 end
              end

              definition_body = [Expr(:(::), field, ty) for (field, ty, _) in pairs]
              constructor_args = [Expr(:(::), field, ty) for (field, _, ty) in pairs]
              arg_names = [field for (field, _, _) in pairs]
              getfields = [:($VAR.$field) for field in arg_names]
              definition_head = :($case{$(generic_tvars...), $(implicit_tvars...)})

              generic_tvars = collect(map(to_tvar, extract_tvars(generic_tvars)))
              implicit_tvars = collect(map(to_tvar, extract_tvars(implicit_tvars)))

              convert_fn = isempty(implicit_tvars) ? nothing : let (=>) = (a, b) -> convert(b, a)
                        out_tvars    = [generic_tvars; implicit_tvars]
                        inp_tvars    = [generic_tvars; [mangle(mod) for _ in implicit_tvars]]
                        fresh_tvars1 = []
                        fresh_tvars2 = []
                        n_generic_tvars = length(generic_tvars)
                        for i in 1 + n_generic_tvars : length(implicit_tvars) + n_generic_tvars
                            TAny = inp_tvars[i]
                            TCov = out_tvars[i]
                            push!(fresh_tvars2, :($TCov <: $TAny))
                            push!(fresh_tvars1, TAny)
                        end
                        arg2 = :($VAR :: $case{$(inp_tvars...)})
                        arg1_cov = :($Type{$case{$(out_tvars...)}})
                        arg1_abs = :($Type{$ret_ty})
                        casted = :($case{$(out_tvars...)}($(getfields...)))
                        quote
                            $Base.convert(::$arg1_cov, $arg2) where {$(generic_tvars...), $(fresh_tvars1...), $(fresh_tvars2...)} = $casted
                            $Base.convert(::$arg1_abs, $arg2) where {$(generic_tvars...), $(fresh_tvars1...), $(fresh_tvars2...)} = $casted
                        end
                    end

            def_cons =
                isempty(generic_tvars) && isempty(implicit_tvars) ?
                    !isempty(constructor_args) ?
                    quote
                        function $case(;$(constructor_args...))
                            $case($(arg_names...))
                        end
                    end                       :
                    nothing        :
                    let spec_tvars = [generic_tvars; [Any for _ in implicit_tvars]]
                        if isempty(generic_tvars)
                            quote
                                function $case($(constructor_args...), )
                                    $case{$(spec_tvars...)}($(arg_names...))
                                end
                            end
                        else
                            quote
                                function $case($(constructor_args...), ) where {$(generic_tvars...)}
                                    $case{$(spec_tvars...)}($(arg_names...))
                                end
                            end
                        end
                    end

              definition = @format [case, l, ret_ty, definition_head] quote
                 struct definition_head <: ret_ty
                     l
                     $(definition_body...)
                 end
                 $def_cons
                 $convert_fn
              end
              push!(defs, (case, pairs, definition))
            end
        end
    end
    defs
end

end # module