text:
module Interact

using Reexport

@reexport using InteractBase
import InteractBase: notifications
import Widgets: Widget, @layout, @nodeps
import Observables: @on, @map!, @map

@reexport using OrderedCollections
@reexport using Observables
@reexport using Knockout
@reexport using CSSUtil
@reexport using WebIO
@reexport using Widgets

struct Bulma<:InteractBase.WidgetTheme; end

const notebookdir = joinpath(@__DIR__, "..", "doc", "notebooks")

const bulma_css = joinpath(@__DIR__, "..", "assets", "bulma.min.css")
const bulma_confined_css = joinpath(@__DIR__, "..", "assets", "bulma_confined.min.css")

function InteractBase.libraries(::Bulma)
    bulmalib = InteractBase.isijulia() ? bulma_confined_css : bulma_css
    vcat(InteractBase.font_awesome, InteractBase.style_css, bulmalib)
end

function __init__()
    InteractBase.registertheme!(:bulma, Bulma())
    settheme!(Bulma())
    nothing
end

end


text:
#
# Subroutine readprm: Reads the parameter files to obtain force constants
#
# L. Martinez, Nov 27, 2013
# Institute of Chemistry - State University of Campinas - Brazil
#
# Tranlsated to Julia in April 15, 2019
#

# Methods to read only one of the components or to call without the flags

readprm!(parfiles, atoms :: Vector{Atom}) = readprm!(parfiles; atoms = atoms)
readprm!(parfiles, atoms :: Vector{Atom}, bonds :: Vector{Bond}) = readprm!(parfiles; atoms = atoms, bonds = bonds)
readprm!(parfiles, atoms :: Vector{Atom}, angles :: Vector{Angle}) = readprm!(parfiles; atoms = atoms, angles = angles)
readprm!(parfiles, atoms :: Vector{Atom}, dihedrals :: Vector{Dihedral}) = readprm!(parfiles; atoms = atoms, dihedrals = dihedrals)
readprm!(parfiles, atoms :: Vector{Atom}, impropers :: Vector{Improper}) = readprm!(parfiles; atoms = atoms, impropers = impropers)
readprm!(parfiles, atoms :: Vector{Atom},
                   bonds :: Vector{Bond},
                   angles :: Vector{Angle},
                   dihedrals :: Vector{Dihedral},
                   impropers :: Vector{Improper}) =
  readprm!(parfiles; atoms = atoms, bonds = bonds, angles = angles, dihedrals = dihedrals, impropers = impropers)

# General function

function readprm!( parfiles;
                   atoms = nothing,
                   bonds = nothing,
                   angles = nothing,
                   dihedrals = nothing,
                   impropers = nothing )

   if typeof(parfiles) == String
     parfiles = [ parfiles ]
   end

   dataname = [ "NONBONDED", "BONDS", "ANGLES", "DIHEDRALS", "IMPROPER", "CMAP", "HBOND" ]
   dataread = [    false,    false,    false,     false,       false,     false,  false  ]

   # Vectors for checking if parameters were found
   if atoms != nothing ; check_atoms = zeros(Bool,length(atoms)) ; end
   if bonds != nothing ; check_bonds = zeros(Bool, length(bonds)) ; end
   if angles != nothing ; check_angles = zeros(Bool, length(angles)) ; end
   if dihedrals != nothing ; check_dihedrals = zeros(Bool, length(dihedrals)) ; end

   for file in parfiles

     file = open(file,"r")
     lines = readlines(file)

     for line in lines

       # Remove comments
       comment = findfirst( isequal('!'), line )
       if comment != nothing
         line = strip(line[1:comment-1])
       else
         line = strip(line)
       end
       if length(line) == 0 ; continue ; end

       # Read data in line
       data = split(line)
       idata = findfirst( isequal(data[1]), dataname )
       if idata != nothing
         @. dataread = false
         dataread[idata] = true
         continue
       end

       #
       # Read atom nonbonded parameter data
       #
       if atoms != nothing && dataread[1]
         for iatom in 1:length(atoms)
           if data[1] == atoms[iatom].type
             eps = parse(Float64,data[3])
             rmin = parse(Float64,data[4])
             eps14, rmin14 = try
               arse(Float64,data[6]), parse(Float64,data[7])
             catch
                0., 0.
             end
             a = atoms[iatom]
             atoms[iatom] = Atom(a.index,a.residue,a.resid,a.name,a.resname,a.segname,a.type,a.element,a.chain,
                                 a.charge,a.mass,eps,rmin,eps14,rmin14,a.backbone,a.coor)
             check_atoms[iatom] = true
           end
         end
         continue
       end

       #
       # Read bond parameter data
       #
       if bonds != nothing && dataread[2]
         for ibond in 1:length(bonds)
           i = bonds[ibond].i
           j = bonds[ibond].j
           if ( ( data[1] == atoms[i].type && data[2] == atoms[j].type ) ||
                ( data[1] == atoms[j].type && data[2] == atoms[i].type ) )
             kb = parse(Float64,data[3])
             b0 = parse(Float64,data[4])
             b = bonds[ibond]
             bonds[ibond] = Bond(i,j,kb,b0)
             check_bonds[ibond] = true
           end
         end
         continue
       end

       #
       # Read angle parameter data
       #
       if angles != nothing && dataread[3]
         for iangle in 1:length(angles)
           i = angles[iangle].i
           j = angles[iangle].j
           k = angles[iangle].k
           if ( ( data[1] == atoms[i].type && data[2] == atoms[j].type && data[3] == atoms[k].type ) ||
                ( data[1] == atoms[k].type && data[2] == atoms[j].type && data[3] == atoms[i].type ) )
             ktheta = parse(Float64,data[4])
             theta0 = parse(Float64,data[5])
             kub, s0 = try
               parse(Float64,data[6]), parse(Float64,data[7])
             catch
               0., 0.
             end
             angles[iangle] = Angle(i,j,k,ktheta,theta0,kub,s0)
             check_angles[iangle] = true
           end
         end
         continue
       end

       #
       # Read dihedral parameter data
       #
       if dihedrals != nothing && dataread[4]
         for idihed in 1:length(dihedrals)
           i = dihedrals[idihed].i
           j = dihedrals[idihed].j
           k = dihedrals[idihed].k
           l = dihedrals[idihed].l
           if ( ( data[1] == atoms[i].type || data[1] == "X" ) &&
                ( data[2] == atoms[j].type ) &&
                ( data[3] == atoms[k].type ) &&
                ( data[4] == atoms[l].type || data[4] == "X" ) ) ||
              ( ( data[4] == atoms[i].type || data[4] == "X" ) &&
                ( data[3] == atoms[j].type ) &&
                ( data[2] == atoms[k].type ) &&
                ( data[1] == atoms[l].type || data[1] == "X" ) ) ||
              ( ( data[1] == atoms[i].type ) &&
                ( data[2] == atoms[j].type || data[2] == "X" ) &&
                ( data[3] == atoms[k].type || data[3] == "X" ) &&
                ( data[4] == atoms[l].type ) ) ||
              ( ( data[4] == atoms[i].type ) &&
                ( data[3] == atoms[j].type || data[3] == "X" ) &&
                ( data[2] == atoms[k].type || data[2] == "X" ) &&
                ( data[1] == atoms[l].type ) )

             mult = dihedrals[idihed].mult + 1
             kchi = Vector{Float64}(undef,mult)
             n = Vector{Int64}(undef,mult)
             delta = Vector{Float64}(undef,mult)

             if mult > 1
               kchi[1:dihedrals[idihed].mult] = dihedrals[idihed].kchi
               n[1:dihedrals[idihed].mult] = dihedrals[idihed].n
               delta[1:dihedrals[idihed].mult] = dihedrals[idihed].delta
             end

             kchi[mult] = parse(Float64,data[5])
             n[mult] = parse(Int64,data[6])
             delta[mult] = parse(Float64,data[7])

             dihedrals[idihed] = Dihedral(i,j,k,l,mult,kchi,n,delta)
             check_dihedrals[idihed] = true
           end
         end
         continue
       end

       #
       # Read improper dihedral parameter data
       #
       if impropers != nothing && dataread[5]
         for iimpr in 1:length(impropers)
           i = impropers[iimpr].i
           j = impropers[iimpr].j
           k = impropers[iimpr].k
           l = impropers[iimpr].l
           if ( ( data[1] == atoms[i].type || data[1] == "X" ) &&
                ( data[2] == atoms[j].type ) &&
                ( data[3] == atoms[k].type ) &&
                ( data[4] == atoms[l].type || data[4] == "X" ) ) ||
              ( ( data[4] == atoms[i].type || data[4] == "X" ) &&
                ( data[3] == atoms[j].type ) &&
                ( data[2] == atoms[k].type ) &&
                ( data[1] == atoms[l].type || data[1] == "X" ) ) ||
              ( ( data[1] == atoms[i].type ) &&
                ( data[2] == atoms[j].type || data[2] == "X" ) &&
                ( data[3] == atoms[k].type || data[3] == "X" ) &&
                ( data[4] == atoms[l].type ) ) ||
              ( ( data[4] == atoms[i].type ) &&
                ( data[3] == atoms[j].type || data[3] == "X" ) &&
                ( data[2] == atoms[k].type || data[2] == "X" ) &&
                ( data[1] == atoms[l].type ) )

             mult = impropers[iimpr].mult + 1
             kpsi = Vector{Float64}(undef,mult)
             psi0 = Vector{Int64}(undef,mult)

             if mult > 1
               kpsi[1:impropers[iimpr].mult] = impropers[iimpr].kpsi
               psi0[1:impropers[iimpr].mult] = impropers[iimpr].psi0
             end

             kpsi[mult] = parse(Float64,data[5])
             psi0[mult] = parse(Float64,data[7])

             impropers[iimpr] = Improper(i,j,k,l,mult,kpsi,psi0)
           end
         end
         continue
       end

     end
     close(file)
   end

   #
   # Report missing parameters
   #
   missing = false
   if atoms != nothing
     for iatom in 1:length(atoms)
       if ! check_atoms[iatom]
         println(" ERROR: Could not fing NONBONDED parameters for: $(atoms[iatom].type)")
         error = true
       end
     end
   end
   if bonds != nothing
     for ibond in 1:length(bonds)
       if ! check_bonds[ibond]
         i = bonds[ibond].i
         j = bonds[ibond].j
         println(" ERROR: Could not fing BOND parameters for: $(atoms[i].type) $(atoms[j].type)" )
         missing = true
       end
     end
   end
   if angles != nothing
     for iangle in 1:length(angles)
       if ! check_angles[iangle]
         i = angles[iangle].i
         j = angles[iangle].j
         k = angles[iangle].k
         println(" ERROR: Could not fing ANGLE parameters for: $(atoms[i].type) $(atoms[j].type) $(atoms[k].type)" )
         missing = true
       end
     end
   end
   if dihedrals != nothing
     for idihed in 1:length(dihedrals)
       if ! check_dihedrals[idihed]
         i = dihedrals[idihed].i
         j = dihedrals[idihed].j
         k = dihedrals[idihed].k
         l = dihedrals[idihed].l
         println(" ERROR: Could not fing DIHEDRAL parameters for: $(atoms[i].type) $(atoms[j].type) $(atoms[k].type) $(atoms[l].type)" )
         missing = true
       end
     end
   end
   if missing
     error("Missing parameters.")
   end

end

ext: jl


text:
include("Segment.jl")
include("PeriodicSegment.jl")
include("Ray.jl")
include("Circle.jl")
include("Line.jl")
include("Arc.jl")

include("UnionDomain.jl")
include("PiecewiseSegment.jl")
include("Curve.jl")

include("Point.jl")

const AffineDomain = Union{AbstractInterval,Segment,PeriodicSegment,Ray,Line}

# These are needed for spaces to auto-convert [a,b] to Interval
function convert(::Type{PeriodicDomain},d::ClosedInterval)
    a,b=d.left,d.right
    a,b = float(a),float(b)
    if isinf(norm(a)) && isinf(norm(b))
        PeriodicLine(d)
    elseif isinf(norm(a)) || isinf(norm(b))
        error("PeriodicRay not implemented")
    else
        PeriodicSegment(d)
    end
end

convert(::Type{Space},d::ClosedInterval) = Space(Domain(d))

#issubset between domains

issubset(a::PeriodicSegment, b::IntervalOrSegment) = Segment(endpoints(a)...)⊆b
issubset(a::IntervalOrSegment, b::PeriodicSegment) = PeriodicSegment(endpoints(a)...)⊆b
issubset(a::IntervalOrSegment{T}, b::PiecewiseSegment{T}) where {T<:Real} =
    a⊆Segment(first(b.points),last(b.points))
issubset(a::IntervalOrSegment, b::Line) = leftendpoint(a)∈b && rightendpoint(a)∈b
issubset(a::Ray{angle}, b::Line{angle}) where angle = leftendpoint(a) ∈ b
issubset(a::Ray{true}, b::Line{false}) = true
issubset(a::Ray{false}, b::Line{true}) = true

function intersect(a::Union{Interval,Segment,Ray},b::Line)
    @assert a ⊆ b
    a
end

function union(a::Union{Interval,Segment,Ray},b::Line)
    @assert a ⊆ b
    b
end

intersect(b::Line,a::Union{Interval,Segment,Ray}) = intersect(a,b)
union(b::Line,a::Union{Interval,Segment,Ray}) = union(a,b)

function setdiff(b::Line,a::Segment)
    @assert a ⊆ b
    if leftendpoint(a)>rightendpoint(a)
        b\reverseorientation(a)
    else
        Ray([leftendpoint(b),leftendpoint(a)]) ∪ Ray([rightendpoint(a),rightendpoint(b)])
    end
end

function setdiff(b::Segment,a::Point)
    if !(a ⊆ b)
        b
    elseif leftendpoint(b) == a.x  || rightendpoint(b) == a.x
        b
    else
        Segment(leftendpoint(b),a.x) ∪ Segment(a.x,rightendpoint(b))
    end
end

# sort

isless(d1::IntervalOrSegment{T1},d2::Ray{false,T2}) where {T1<:Real,T2<:Real} = d1 ≤ d2.center
isless(d2::Ray{true,T2},d1::IntervalOrSegment{T1}) where {T1<:Real,T2<:Real} = d2.center ≤ d1

## set minus
function Base.setdiff(d::AffineDomain,ptsin::UnionDomain{AS}) where {AS <: AbstractVector{P}} where {P <: Point}
    pts=Number.(elements(ptsin))
    isempty(pts) && return d
    tol=sqrt(eps(arclength(d)))
    da=leftendpoint(d)
    isapprox(da,pts[1];atol=tol) && popfirst!(pts)
    isempty(pts) && return d
    db=rightendpoint(d)
    isapprox(db,pts[end];atol=tol) && pop!(pts)

    sort!(pts)
    leftendpoint(d) > rightendpoint(d) && reverse!(pts)
    filter!(p->p ∈ d,pts)

    isempty(pts) && return d
    length(pts) == 1 && return d \ pts[1]

    ret = Array{Domain}(undef, length(pts)+1)
    ret[1] = Domain(leftendpoint(d) .. pts[1])
    for k = 2:length(pts)
        ret[k] = Domain(pts[k-1]..pts[k])
    end
    ret[end] = Domain(pts[end] .. rightendpoint(d))
    UnionDomain(ret)
end

function Base.setdiff(d::SegmentDomain,p::Point)
    x = Number(p)
    (x ∉ d || x ≈ leftendpoint(d) || x ≈ rightendpoint(d)) && return d
    DifferenceDomain(d,p)
end

# multivariate domainxs

include("multivariate.jl")
include("Disk.jl")


text:
# ============================================================
# Tracker Methods (see TrackerFlux.jl for details)
# ============================================================

track(m) = fmap(x -> x isa AbstractArray ? Tracker.param(x) : x, m)

untrack(m) = fmap(Tracker.data, m)

function Flux.Optimise.update!(opt, xs::Tracker.Params, gs)
    for x in xs
        Flux.Optimise.update!(opt, x, gs[x])
    end
end

function tracker_mode()
    @eval Flux.Optimise.update!(opt, x, x̄) =
        Tracker.update!(x, -Flux.Optimise.apply!(opt, Tracker.data(x), Tracker.data(x̄)))
    @eval Flux.gradient(f, args...) = Tracker.gradient(f, args...)
end

function zygote_mode()
    @eval Flux.Optimise.update!(opt, x, x̄) =
        Flux.Optimise.update!(x, -Flux.Optimise.apply!(opt, x, x̄))
    @eval Flux.gradient(f, args...) = Flux.Zygote.gradient(f, args...)
end


text:
struct CarRacingEnvParams{T}
    m::T
    Izz::T
    h_cm::T
    l_f::T
    l_r::T
    C_D0::T
    C_D1::T
    C_αf::T
    C_αr::T
    μ_f::T
    μ_r::T
    δ_max::T
    δ_dot_max::T
    Fx_max::T
    Fx_min::T
    λ_brake::T
    λ_drive::T
    β_limit::T
end

Base.show(io::IO, params::CarRacingEnvParams) = print(
    io,
    join(["$p=$(getfield(params, p))" for p in fieldnames(CarRacingEnvParams)], ","),
)

mutable struct CarRacingEnv{A,T,R<:AbstractRNG} <: AbstractEnv
    params::CarRacingEnvParams{T}
    action_space::A
    observation_space::Space{Vector{ClosedInterval{T}}}
    state::Vector{T}
    done::Bool
    t::Int
    dt::T
    δt::T
    track::Track
    rng::R
end

"""
    CarRacingEnv(;kwargs...)

# Keyword arguments
- `T = Float64`
- `m = 2000.0, 			                        # Mass (kg)`
- `I_zz = 3764.0, 		                    	# Moment of inertia in vertical direction (kg m^2)`
- `h_cm = 0.3,	 			                    # Height of CoM above ground (m)`
- `l_f = 1.53, 			                        # Distance from CoM to front axle (m) ["a" in paper]`
- `l_r = 1.23, 			                        # Distance from CoM to rear axle (m) ["b" in paper]`
- `C_D0 = 241, 		    	                    # Constant drag term (N)`
- `C_D1 = 25.1, 	                        	# Linear drag coefficient (N s/m)`
- `C_αf = 150000.0,     	                    # Front cornering stiffness (N/rad)`
- `C_αr = 280000.0,                         	# Rear cornering stiffness (N/rad)`
- `μ_f = 0.9, 			                        # Front tire friction`
- `μ_r = 0.9, 				                    # Rear tire friction`
- `δ_max = deg2rad(18),	                        # Steering angle limits (rad)`
- `δ_dot_max = deg2rad(90),                    	# Stearing angle rate limits (rad/s)`
- `Fxmax = 7200,                                # maximum force due to motor torque (N)`
- `Fxmin = 22500                                # maximum force due to breaking`
- `λ_brake = 0.6,                               # Fixed brake distribution`
- `λ_drive = 0.0,                               # Fixed drive distribution`
- `dt = 0.1`,                                   # Time step between actions
- `δt = 0.01,                                   # Time step used for integration`
- `β_limit = deg2rad(45)                        # Beta penalty limit`
- `track = string(@__DIR__, "/car_racing_tracks/curve.csv"),   # Track to load`
****** TO-DO Add lane_width
- `rng = Random.GLOBAL_RNG`
"""
function CarRacingEnv(;
    T = Float64,
    m = 2000.0,
    I_zz = 3764.0,
    h_cm = 0.3,
    l_f = 1.53,
    l_r = 1.23,
    C_D0 = 241,
    C_D1 = 25.1,
    C_αf = 150000.0,
    C_αr = 280000.0,
    μ_f = 0.9,
    μ_r = 0.9,
    δ_max = deg2rad(18),
    δ_dot_max = deg2rad(90),
    Fx_max = 7200.0,
    Fx_min = 22500.0,
    λ_brake = 0.6,
    λ_drive = 0.0,
    dt = 0.1,
    δt = 0.01,
    β_limit = deg2rad(45),
    track = string(@__DIR__, "/car_racing_tracks/curve.csv"),
    track_sample_factor = 20,
    rng = Random.GLOBAL_RNG,
)
    params = CarRacingEnvParams{T}(
        m,
        I_zz,
        h_cm,
        l_f,
        l_r,
        C_D0,
        C_D1,
        C_αf,
        C_αr,
        μ_f,
        μ_r,
        δ_max,
        δ_dot_max,
        Fx_max,
        Fx_min,
        λ_brake,
        λ_drive,
        β_limit,
    )
    return CarRacingEnv(params, T=T, dt=dt, δt=δt, track=track, track_sample_factor=track_sample_factor, rng=rng)
end

"""
    CarRacingEnv(params::CarRacingEnvParams;kwargs...)

# Keyword arguments
- `T = Float64`
- `dt = 0.1`,                                                   # Time step between actions
- `δt = 0.01,                                                   # Time step used for integration`
- `track = string(@__DIR__, "/car_racing_tracks/curve.csv"),   # Track to load`
- `rng = Random.GLOBAL_RNG`
"""
function CarRacingEnv(
    params::CarRacingEnvParams;
    T = Float64,
    dt = 0.1,
    δt = 0.01,
    track = string(@__DIR__, "/car_racing_tracks/curve.csv"),
    track_sample_factor = 10,
    rng = Random.GLOBAL_RNG,
)

    action_space = ClosedInterval{Vector{T}}(
        [-1.0, -1.0],
        [ 1.0,  1.0],
        )
    observation_space = Space([
        -Inf..Inf,                          # X position in XY plane (x = north, y = west)
        -Inf..Inf,                          # Y position in XY plane (x = north, y = west)
        -π..π,                              # yaw (rotation from x axis toward y axis [north to west])
        -Inf..Inf,                          # Longitudinal velocity
        -Inf..Inf,                          # Lateral velocity
        -Inf..Inf,                          # yaw rate
        -params.δ_max..params.δ_max,        # steering angle
        -1.0..1.0,                          # acceleration/brake amount [-1, 1]
        ])

    env = CarRacingEnv(
        params,
        action_space,
        observation_space,
        zeros(T,8),
        false,
        0,
        dt,
        δt,
        Track(track, sample_factor=track_sample_factor),
        rng,
        )

    reset!(env)
    env
end

CarRacingEnv{T}(; kwargs...) where {T} = CarRacingEnv(; T = T, kwargs...)

Random.seed!(env::CarRacingEnv, seed) = Random.seed!(env.rng, seed)

RLBase.action_space(env::CarRacingEnv) = env.action_space
RLBase.state_space(env::CarRacingEnv{T}) where {T} = env.observation_space
RLBase.is_terminated(env::CarRacingEnv) = env.done
RLBase.state(env::CarRacingEnv) = env.state

function within_track(env::CarRacingEnv)
    return within_track(env.track, env.state[1:2])
end
function calculate_β(env::CarRacingEnv)
    return atan(env.state[5], env.state[4])
end
function exceed_β(env::CarRacingEnv)
    if abs(calculate_β(env)) > env.params.β_limit
        return true
    end
    return false
end

"""
    reward(env::CarRacingEnv)

    - Reward for going fast, penalty for boundary violation,
        excess β and distance from center of track
    - Boundary violation    : -1000000
    - Excessive β           :  -5000
    - Distance from center  : -(x² + y²)^(1/2)
    - Speed                 : +2.0 * (Vx² + Vy²)^(1/2)
"""
function RLBase.reward(env::CarRacingEnv{T}) where {T}
    rew = 0.0
    within_tuple = within_track(env)
    if !within_tuple.within
        rew += -1000000.0
    end
    if exceed_β(env)
        rew += -5000.0
    end
    rew += -within_tuple.dist
    rew += 2.0*norm(env.state[4:5])
    return rew
end

function RLBase.reset!(env::CarRacingEnv{A,T}) where {A,T}
    ss_size = length(env.state)
    env.state = zeros(T, ss_size)
    env.state[3] = deg2rad(90)
    env.state[4] = 10.0
    env.t = 0
    env.done = false
    nothing
end

function RLBase.reset!(env::CarRacingEnv{A,T}, state::Vector{T}) where {A,T}
    env.state = state
    env.t = 0
    env.done = false
    nothing
end

"""
    a::Vector{Float64}
    size = (2,)
    a[1] = Turn angle [-max turn angle, max turn angle] (-1 right turn, +1 left turn)
    a[2] = Pedal amount (-1 = full brake, 1 = full throttle)
"""
function (env::CarRacingEnv{<:ClosedInterval})(a::Vector{Float64})
    a in env.action_space || error("Action is not in action space")
    _step!(env, a)
end

function (env::CarRacingEnv{<:ClosedInterval})(a::Vector{Int})
    env(Float64.(a))
end

function (env::CarRacingEnv{<:ClosedInterval})(a::Matrix{Float64})
    size(a)[2] == 1 || error("Only implented for one step")
    env(vec(a))
end

function calc_tire_fy(α, μ, C_α, fzt, fxt)
    fy_max = sqrt(max((μ * fzt)^2 - fxt^2, 1e-8))
    ta = tan(α)
    if abs(α) < atan(3 * fy_max / C_α)
        return -C_α*ta + (C_α^2/(3*fy_max))*abs(ta)*ta - ((C_α^3)/(27*fy_max^2))*ta^3
    else
        return -fy_max * sign(α)
    end
end

function calc_tire_fz(params::CarRacingEnvParams, fx, tire::Char)
    mass = params.m
    l_t = params.l_f
    h_cm = params.h_cm
    L = params.l_r + params.l_f
    if tire == 'f'
        l_t = params.l_r
        h_cm *= -1
    end
    return (mass * l_t * 9.81 + h_cm * fx) / L
end

"""
Planar single-track model with tire forces coming from a
single-friction-coefficient brush tire model.

M. Brown and J. C. Gerdes, Coordinating Tire Forces to Avoid Obstacles Using
Nonlinear Model Predictive Control, in IEEE Transactions on Intelligent
Vehicles, vol. 5, no. 1, pp. 21-31, March 2020, doi: 10.1109/TIV.2019.2955362.
"""
function _step!(env::CarRacingEnv, a::Vector{Float64})
    env.t += 1
    x = env.state[1]
    y = env.state[2]
    Ψ = env.state[3]
    Vx = env.state[4]
    Vy = env.state[5]
    Ψ_dot = env.state[6]
    δ = env.state[7]

    dt = env.dt
    δt = env.δt

    commanded_Δδ_rate = abs(a[1]*env.params.δ_max - δ) / dt # Commanded turn rate
    Δδ_rate = min(commanded_Δδ_rate, env.params.δ_dot_max) * sign(a[1]*env.params.δ_max - δ)
    pedal = a[2] # Assume instant torque for accel or brake

    integration_steps = round(Int, dt/δt)
    for _ in 1:integration_steps
        δ += Δδ_rate * δt

        # Slip angles for the tires
        α_f = atan((Vy + env.params.l_f*Ψ_dot) , Vx) - δ
        α_r = atan((Vy - env.params.l_r*Ψ_dot) , Vx)

        # Simple drag
        fx_aero = (env.params.C_D0 + env.params.C_D1*abs(Vx)) * sign(Vx)

        accel = env.params.Fx_max * max(pedal, 0.0)
        brake = env.params.Fx_min * min(pedal, 0.0) * sign(Vx) # Opposite direction of long velocity
        fx = accel + brake

        # Distribution of the forces on the front and rear tires
        fxf = (pedal <= 0 ? env.params.λ_brake : env.params.λ_drive) * fx
        fxr = (1 - (pedal <= 0 ? env.params.λ_brake : env.params.λ_drive)) * fx
        fzf = calc_tire_fz(env.params, fx, 'f')
        fzr = calc_tire_fz(env.params, fx, 'r')
        fyf = calc_tire_fy(α_f, env.params.μ_f, env.params.C_αf, fzf, fxf)
        fyr = calc_tire_fy(α_r, env.params.μ_r, env.params.C_αr, fzr, fxr)

        Ψ_ddot = (1/env.params.Izz) * (env.params.l_f * (fxf*sin(δ) + fyf*cos(δ))  - env.params.l_r*fyr)
        Vy_dot = (1/env.params.m) * (fyf*cos(δ) + fxf*sin(δ) + fyr) - Ψ_dot*Vx
        Vx_dot = (1/env.params.m) * (fxf*cos(δ) - fyf*sin(δ) + fxr - fx_aero) + Ψ_dot*Vy

        Ψ_dot += Ψ_ddot * δt                # Updated yaw rate
        Vx += Vx_dot * δt                   # Updated longitudinal velocity
        Vy += Vy_dot * δt                   # Updated lateral velocity
        Ψ += Ψ_dot * δt                     # Updated yaw (heading)
        Ψ = atan(sin(Ψ), cos(Ψ))
        x += (Vx*cos(Ψ) - Vy*sin(Ψ)) * δt   # Updated x position yaw heading is counterclockwise, x is north, y is west
        y += (Vx*sin(Ψ) + Vy*cos(Ψ)) * δt   # Updated y position yaw is counterclockwise, x is north, y is west
    end

    env.state[1] = x
    env.state[2] = y
    env.state[3] = Ψ
    env.state[4] = Vx
    env.state[5] = Vy
    env.state[6] = Ψ_dot
    env.state[7] = δ
    env.state[8] = pedal
    return env
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
@testset "Grid generation" begin

    @testset "Regular pattern" begin

        @test collect(MicroscopyTools.regular_pattern((4, 4), offset = -1, step = 2)) == [(-1, -1) (-1, 1); (1, -1) (1, 1)]

        @test collect(MicroscopyTools.regular_pattern((4, 4), step = 2)) == [(1, 1) (1, 3); (3, 1) (3, 3)]

        @test collect(MicroscopyTools.regular_pattern((3, 3))) == [(1, 1) (1, 2) (1, 3); (2, 1) (2, 2) (2, 3); (3, 1) (3, 2) (3, 3)]

    end

end

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
@testset "Narrow Picard group of non-maximal orders" begin
  D = [-97, -95, -94, -93, -91, -89, -87, -86, -85, -83, -82, -79, -78, -77, -74, -73, -71, -70, -69, -67, -66, -65, -62, -61, -59, -58, -57, -55, -53, -51, -47, -46, -43, -42, -41, -39, -38, -37, -35, -34, -33, -31, -30, -29, -26, -23, -22, -21, -19, -17, -15, -14, -13, -11, -10, -7, -6, -5, -3, -2, 1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, 17, 19, 21, 22, 23, 26, 29, 30, 31, 33, 34, 35, 37, 38, 39, 41, 42, 43, 46, 47, 51, 53, 55, 57, 58, 59, 61, 62, 65, 66, 67, 69, 70, 71, 73, 74, 77, 78, 79, 82, 83, 85, 86, 87, 89, 91, 93, 94, 95, 97]
  Qx, x = QQ["x"]
  for d in D
    OK = maximal_order(number_field(x^2 + d)[1])
    C = narrow_class_group(OK)[1]
    CC = Hecke.narrow_picard_group(OK)[1]
    @test is_isomorphic(C, CC)
  end

  K = number_field(x^4 + 2*x^3 - 35*x^2 - 36*x + 5, "a", cached = false)[1]
  E = equation_order(K)
  @test order(picard_group(E)[1]) == 1
  C, _dlog, _exp = Hecke.narrow_picard_group(E)
  for i in 1:10
    c = rand(C)
    @test _dlog(_exp(c)) == c
  end

  Qx, x = QQ["x"]
  K = number_field(x^2 - 4 * 9 * 13, "a", cached = false)[1]
  E = equation_order(K)
  C, _dlog, _exp = Hecke.narrow_picard_group(E)
  for i in 1:10
    c = rand(C)
    @test _dlog(_exp(c)) == c
  end

  f = x^2 - 136*x + 4590
  E = equation_order(f, cached = false)
  C, _dlog, _exp = Hecke.narrow_picard_group(E)
  @test length(unique(_exp.(C))) == 4
end


text:
# test assignment problem

using DocStringExtensions
using Distributions
using KernelDensityEstimate, KernelDensityEstimatePlotting
using Munkres
using TransformUtils

using Optim

# new cost
function S_min_kld(p::BallTreeDensity,q::BallTreeDensity)
  minimum(abs.([kld(p,q), kld(q,p)]))
end

objCost(cst::Array{Float64,2}, asgn::Vector{Int}) = sum(Float64[cst[i,j] for (i,j) in enumerate(asgn[asgn.>0])])

N=100

noiseOld = diagm([0.25;0.25].^2)
noiseNex = diagm([1.0;1.0].^2)

Tnoise = [randn(2); 0.1*randn()]
Tnoise = [0.5; 0.7; -0.1]

# old and new positions
old = [0.0 0 1; 1 1 1; 0 3 1; 0 4 1]'
nex = SE2(Tnoise)*old

# occlusion plus new
old2 = [0.0 0 1; 1 1 1; 2 3 1; 0 4 1]'
nex2 = SE2(Tnoise)*old2

OLD = [kde!(rand(MvNormal(old[1:2,1], noiseOld),N));
       kde!(rand(MvNormal(old[1:2,2], noiseOld),N));
       kde!(rand(MvNormal(old[1:2,3], noiseOld),N));
       kde!(rand(MvNormal(old[1:2,4], noiseOld),N)) ]

NEX = [kde!(rand(MvNormal(nex[1:2,1], noiseNex),N));
       kde!(rand(MvNormal(nex[1:2,2], noiseNex),N));
       kde!(rand(MvNormal(nex[1:2,3], noiseNex),N));
       kde!(rand(MvNormal(nex[1:2,4], noiseNex),N)) ]

NEX2 = [kde!(rand(MvNormal(nex2[1:2,1], noiseNex),N));
        kde!(rand(MvNormal(nex2[1:2,2], noiseNex),N));
        kde!(rand(MvNormal(nex2[1:2,3], noiseNex),N));
        kde!(rand(MvNormal(nex2[1:2,4], noiseNex),N)) ]

#
plotKDE(OLD, dims=[1;2], levels=2, c=["cyan"])
#
plotKDE(NEX2, dims=[1;2], levels=2, c=["red"])

"""
    $(SIGNATURES)

Return cost matrix with old elements as rows and new elements being assigned as columns.
"""
function assignmentCost(OLDl, NEXl, z::Vector{Float64})
  cost = zeros(length(OLDl),length(NEXl))
  xTo = inv(SE2(z))
  pts = ones(3,100)
  for i in 1:4, j in 1:4
    pts[1:2,:] = getPoints(NEX[j])
    nPts = xTo*pts
    nK = kde!(nPts[1:2,:])
    cost[i,j] = S_min_kld(OLDl[i], nK)
  end
  cost
end

"""
    $(SIGNATURES)

Determine base cost for new landmarks.
"""
function getNewElemCost(OLDl, NEXl;quant=0.75, xTo=zeros(3))
  c = assignmentCost(OLDl, NEXl, xTo)
  nelemcost = Float64[quantile(c[:,i],quant) for i in 1:length(NEXl)]'
  # nelemcost = Base.mean(c,1)
  repmat(nelemcost, length(NEXl),1) - 0.5*diagm(nelemcost')
end

function findAssignmentCost(OLDl, NEXl, z::Vector{Float64}; costNew::Union{Void, Array{Float64,2}}=nothing, quant=0.75)
  costNew = costNew == nothing ? getNewElemCost(OLDl, NEXl;quant=0.75, xTo=zeros(3)) : costNew
  cost = assignmentCost(OLDl, NEXl, z)
  cost = [cost; costNew]
  comp = munkres(cost)
  comp, objCost(cost, comp), cost
end

ep = 1e-5

costnew=getNewElemCost(OLD, NEX2;quant=0.5, xTo=zeros(3))
a1,c1,C1 = findAssignmentCost(OLD, NEX2, [0.0;0.0;-0.0], costNew=costnew, quant=0.01)
a2,c2,C2 = findAssignmentCost(OLD, NEX2, [0.0;0.0;-0.0+ep])

dfdx = (c2-c1)/ep

# base cost for assigning existing landmarks only
costnew=getNewElemCost(OLD, NEX2, quant=0.53, xTo=zeros(3))
# setup lambda for optmization
gg(x) = findAssignmentCost(OLD, NEX2, x, costNew=costnew)[2]
# find optimum solution
opm = optimize(gg, zeros(3))
# recover optimal assignment
@show opm.minimizer
oAsgn,oCost,oCostMat = findAssignmentCost(OLD, NEX2, opm.minimizer, costNew=costnew)

# from Munkres.jl

cost = rand(4,4)
best = [3;4;1;2]
for (i,j) in enumerate(best); cost[i,j]=0; end
comp = munkres(cost)
@assert best == comp

objCost(cost, comp)

# trying stuff

S_min_kld(OLD[1], NEX[4])

pts = randn(2,100); pts[2,1:50] += 4
testd = kde!(pts)
S_min_kld(NEX[4], testd)

#


text:
function dynmodel!(x::Vector, u::Vector, dt::Number, params::Vector)
    ## Parameters and constants
    # motor
    k1_left = exp(params[1]) # Motor constant: current -> torque
    k1_right = exp(params[2]) # Motor constant: current -> torque
    k2_left = exp(params[3]) # Motor constant: rot.speed -> induced voltage
    k2_right = exp(params[4]) # Motor constant: rot.speed -> induced voltage
    # segway
    Mb = 1.76 # kg # Masse des Körpers
    Mw = 0.147 # kg # Masse der Räder
    R = 0.07 # m #Radius der Räder
    cz = 0.077 # m #Schwerpunkt des Körpers Höhe über Radachse
    b = 0.1985/2 # m #Half the distance between the wheels (OO_w in paper)
    Ixx = (0.166^2+0.21^2)*Mb/12 + Mb*(cz-0.07)^2 # kg m^2 #Inertia of the body
    Iyy = (0.072^2+0.21^2)*Mb/12 + Mb*cz^2 # kg m^2 #Inertia of the body
    Izz = (0.072^2+0.166^2)*Mb/12 # kg m^2 #Inertia of the body
    Iwa = Mw*R^2/2 # kg m^2 #Inertia of the wheels about their axis
    Iwd = Mw*b^2 # kg m^2 #Inertia of the wheels about their axis
    seg = Segway_unmutable(Mb, Mw, R, cz, b, Ixx, Iyy, Izz, Iwa, Iwd)
    # friction
    fl_1 = exp(params[5]) # left wheel
    fl_2 = exp(params[6]) # left wheel
    fr_1 = exp(params[7]) # right wheel
    fr_2 = exp(params[8]) # right wheel

    ## Evaluate dynamics
    # get wheel rotation speeds
    v_left = x[6]-x[7]*b
    v_right = x[6]+x[7]*b
    omega_right = v_right/R
    omega_left = v_left/R
    # apply simplified motor model
    tau_right = k1_right*(u[2]-omega_right*k2_right)
    tau_left = k1_left*(u[1]-omega_left*k2_left)
    # apply simplified friction model
    tau_right = frictionmodel(tau_right, omega_right, fr_1, fr_2)
    tau_left = frictionmodel(tau_left, omega_left, fl_1, fl_2)
    # apply segway model
    tau = [tau_right, tau_left]
    segway_rk4!(x, tau, dt, seg)

    return x
end

function dynmodel(x::Vector, u::Vector, dt::Number, params::Vector)
    x2 = copy(x)
    dynmodel!(x2, u, dt, params)
    return x2
end

@inline function frictionmodel(tau, omega, param1, param2)
    # Viskoser Reibterm
    dtau_v = param1*omega

    # Haftreibung
    dtau_h = param2*sign(omega)

    return tau-dtau_v-dtau_h
end

"""
State vector is [x0, y0, phi, alpha, dalpha, v, dphi, ir, il]
"""
function segway_cost(xt::AbstractVector, u::AbstractVector)
    cost = (xt[4]^2 + 0.1*xt[5]^2
            + 10*xt[2]^2 + 10*xt[1]^2 # or squared here...
            + 0.1* xt[7]^2
            + 0.01*(u'u)) # punish for falling over

    reward = 0.0 #reward for moving along the x-axis

    if xt[4]^2 > 2 #check if lager than pi/2
        terminal = 1
        cost *= 1.2
    else
        terminal = 0
    end

    return cost-reward, terminal
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
@testset "show" begin
    x = CSTParser.parse("a + (b*c) - d")
    @test sprint(show, x) ===
    "  1:13  call\n  1:2    OP: -\n  3:12   call\n  3:4     OP: +\n  5:6     a\n  7:12    brackets\n  7:9      call\n  7:7       OP: *\n  8:8       b\n  9:9       c\n 13:13   d"
end


text:
module Morton

export cartesian2morton, cartesian3morton
export morton2cartesian, morton3cartesian
export tree2morton, tree3morton
export morton2tree, morton3tree
export tree2cartesian, tree3cartesian
export cartesian2tree, cartesian3tree

### https://fgiesen.wordpress.com/2009/12/13/decoding-morton-codes/

# "Insert" a 0 bit after each of the 16 low bits of x
function _Part1By1(x::Integer)
    x &= 0x0000ffff                   # x = ---- ---- ---- ---- fedc ba98 7654 3210
    x = (x | (x <<  8)) & 0x00ff00ff  # x = ---- ---- fedc ba98 ---- ---- 7654 3210
    x = (x | (x <<  4)) & 0x0f0f0f0f  # x = ---- fedc ---- ba98 ---- 7654 ---- 3210
    x = (x | (x <<  2)) & 0x33333333  # x = --fe --dc --ba --98 --76 --54 --32 --10
    x = (x | (x <<  1)) & 0x55555555  # x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
    x+1
end

# "Insert" two 0 bits after each of the 10 low bits of x
function _Part1By2(x::Integer)
    x &= 0x000003ff                   # x = ---- ---- ---- ---- ---- --98 7654 3210
    x = (x | (x << 16)) & 0xff0000ff  # x = ---- --98 ---- ---- ---- ---- 7654 3210
    x = (x | (x <<  8)) & 0x0300f00f  # x = ---- --98 ---- ---- 7654 ---- ---- 3210
    x = (x | (x <<  4)) & 0x030c30c3  # x = ---- --98 ---- 76-- --54 ---- 32-- --10
    x = (x | (x <<  2)) & 0x09249249  # x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
    x+1
end

"""
    cartesian2morton(c::Vector) -> m::Integer

Given a 2-D Cartesian coordinate, return the corresponding Morton number.

# Examples
```jldoctest
julia> cartesian2morton([5,2])
19
```
"""
function cartesian2morton(c::AbstractVector{T}) where T<:Integer
    (_Part1By1(c[2]-1) << 1) + _Part1By1(c[1]-1) - 2
end

"""
    cartesian3morton(c::AbstractVector) -> m::Integer

Given a 3-D Cartesian coordinate, return the corresponding Morton number.

# Examples
```jldoctest
julia> cartesian3morton([5,2,1])
67
```
"""
function cartesian3morton(c::AbstractVector{T}) where T<:Integer
    (_Part1By2(c[3]-1) << 2) + (_Part1By2(c[2]-1) << 1) + _Part1By2(c[1]-1) - 6
end

function _Compact1By1(x::Integer)
    x &= 0x55555555                   # x = -f-e -d-c -b-a -9-8 -7-6 -5-4 -3-2 -1-0
    x = (x | (x >>  1)) & 0x33333333  # x = --fe --dc --ba --98 --76 --54 --32 --10
    x = (x | (x >>  2)) & 0x0f0f0f0f  # x = ---- fedc ---- ba98 ---- 7654 ---- 3210
    x = (x | (x >>  4)) & 0x00ff00ff  # x = ---- ---- fedc ba98 ---- ---- 7654 3210
    x = (x | (x >>  8)) & 0x0000ffff  # x = ---- ---- ---- ---- fedc ba98 7654 3210
    x+1
end

function _Compact1By2(x::Integer)
    x &= 0x09249249                   # x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
    x = (x | (x >>  2)) & 0x030c30c3  # x = ---- --98 ---- 76-- --54 ---- 32-- --10
    x = (x | (x >>  4)) & 0x0300f00f  # x = ---- --98 ---- ---- 7654 ---- ---- 3210
    x = (x | (x >>  8)) & 0xff0000ff  # x = ---- --98 ---- ---- ---- ---- 7654 3210
    x = (x | (x >> 16)) & 0x000003ff  # x = ---- ---- ---- ---- ---- --98 7654 3210
    x+1
end

"""
    morton2cartesian(m::Integer) -> [x,y]

Given a Morton number, return the corresponding 2-D Cartesian coordinates.

# Examples
```jldoctest
julia> morton2cartesian(19)
2-element Array{Int64,1}:
 5
 2
```
"""
function morton2cartesian(m::Integer)
    m -= 1
    [_Compact1By1(m>>0), _Compact1By1(m>>1)]
end

"""
    morton3cartesian(m::Integer) -> [x,y,z]

Given a Morton number, return the corresponding 3-D Cartesian coordinates.

# Examples
```jldoctest
julia> morton3cartesian(67)
3-element Array{Int64,1}:
 5
 2
 1
```
"""
function morton3cartesian(m::Integer)
    m -= 1
    [_Compact1By2(m>>0), _Compact1By2(m>>1), _Compact1By2(m>>2)]
end

function _treeNmorton(t::AbstractVector{T}, ndim::Integer) where T<:Integer
    n=m=0
    it=length(t)
    ndim2=2^ndim
    while it>0
      m += (t[it]-1)*ndim2^n
      n += 1
      it-=1
    end
    m+1
end

"""
    tree2morton(t::AbstractVector) -> m::Integer

Given a quadtree coordinate, return the corresponding Morton number.

# Examples
```jldoctest
julia> tree2morton([2,1,3])
19
```
"""
tree2morton(t::AbstractVector{T}) where T<:Integer = _treeNmorton(t,2)

"""
    tree3morton(t::AbstractVector) -> m::Integer

Given a octree coordinate, return the corresponding Morton number.

# Examples
```jldoctest
julia> tree3morton([2,1,3])
67
```
"""
tree3morton(t::AbstractVector{T}) where T<:Integer = _treeNmorton(t,3)

function _mortonNtree(m::T, ndim::Integer) where T<:Integer
    t=T[]
    ndim2=2^ndim
    while true
        d,r = [divrem(m-1,ndim2)...]+[1,1]
        pushfirst!(t,r)
        d==1 && break
        m=d
    end
    t
end

"""
    morton2tree(m::Integer) -> t::AbstractVector

Given a Morton number, return the corresponding quadtree coordinate.

# Examples
```jldoctest
julia> morton2tree(19)
3-element Array{Any,1}:
 2
 1
 3
```
"""
morton2tree(m::Integer) = _mortonNtree(m,2)

"""
    morton3tree(m::Integer) -> t::AbstractVector

Given a Morton number, return the corresponding octree coordinate.

# Examples
```jldoctest
julia> morton3tree(67)
3-element Array{Any,1}:
 2
 1
 3
```
"""
morton3tree(m::Integer) = _mortonNtree(m,3)

function _treeNcartesian(t::AbstractVector{T}, ndim::Integer) where T<:Integer
    c = [((t[1]-1)>>b)&1+1 for b in 0:ndim-1]
    if length(t)>1
        cn = _treeNcartesian(t[2:end], ndim)
        m = 2^(length(t)-1)
        return [m*(c[i]-1)+cn[i] for i in 1:ndim]
    else
        return c
    end
end

"""
   tree2cartesian(t::AbstractVector) -> c::AbstractVector

Given quadtree coordinate, return the corresponding 2-D Cartesian coordinate.

# Examples
```jldoctest
julia> tree2cartesian([2,1,3])
2-element Array{Int64,1}:
 5
 2
```
"""
tree2cartesian(t::AbstractVector{T}) where T<:Integer = _treeNcartesian(t, 2)

"""
   tree3cartesian(t::AbstractVector) -> c::AbstractVector

Given octree coordinate, return the corresponding 3-D Cartesian coordinate.

# Examples
```jldoctest
julia> tree3cartesian([2,1,3])
3-element Array{Int64,1}:
 5
 2
 1
```
"""
tree3cartesian(t::AbstractVector{T}) where T<:Integer = _treeNcartesian(t, 3)

function _cartesianNtree(c::AbstractVector{T}, half, ndim::Integer) where T<:Integer
    t = 1
    for d=1:ndim
        t += 2^(d-1)*(c[d]>half)
    end
    if half>1
        return [t, _cartesianNtree(map(x->(x-1)%half+1,c), half>>1, ndim)...]
    else
        return [t]
    end
end

"""
   cartesian2tree(c::AbstractVector) -> t::AbstractVector

Given a 2-D Cartesian coordinate, return the corresponding quadtree coordinate.

# Examples
```jldoctest
julia> cartesian2tree([5,2])
3-element Array{Int64,1}:
 2
 1
 3
```
"""
cartesian2tree(c::AbstractVector{T}) where T<:Integer =
      _cartesianNtree(c, max(2,nextpow(2, widen(maximum(c))))>>1, 2)

"""
   cartesian3tree(c::AbstractVector) -> t::AbstractVector

Given a 3-D Cartesian coordinate, return the corresponding octree coordinate.

# Examples
```jldoctest
julia> cartesian3tree([5,2,1])
3-element Array{Int64,1}:
 2
 1
 3
```
"""
cartesian3tree(c::AbstractVector{T}) where T<:Integer =
      _cartesianNtree(c, max(2,nextpow(2, widen(maximum(c))))>>1, 3)

end # module


text:
"""
Abstract type used as base type for the type created by the `@resumable` macro.
"""
abstract type FiniteStateMachineIterator{R} end

"""
Mutable struct that contains a single `UInt8`.
"""
mutable struct BoxedUInt8
  n :: UInt8
end

"""
Implements the `iteratorsize` method of the *iterator* interface for a subtype of `FiniteStateMachineIterator`.
"""
Base.IteratorSize(::Type{T}) where T<:FiniteStateMachineIterator = Base.SizeUnknown()

"""
Implements the `eltype` method of the *iterator* interface for a subtype of `FiniteStateMachineIterator`.
"""
Base.eltype(::Type{T}) where T<:FiniteStateMachineIterator{R} where R = R

function Base.iterate(fsm_iter::T, state::UInt8=0x00) where T<:FiniteStateMachineIterator
  fsm_iter._state = state
  result = fsm_iter()
  fsm_iter._state == 0xff && return nothing
  result, fsm_iter._state
end


text:
module MonkeyLang

export start_repl, @monkey_eval_str, @monkey_vm_str

using Printf

const MONKEY_VERSION = v"0.2.1"
const MONKEY_AUTHOR = "Gabriel Wu"

# Tokens
include("token.jl")

# AST Nodes
include("ast.jl")

# Op Codes and Instructions
include("code.jl")

# Wrapped Objects
include("object.jl")

# Lexer
include("lexer.jl")

# Parser
include("parser.jl")

# Builtin Functions
include("builtins.jl")

# Evaluator
include("evaluator.jl")

# Auxiliary AST Modification Functions Used for Macros
include("ast_modify.jl")

# Symbol and Symbol Table
include("symbol_table.jl")

# Analyzer
include("analyzer.jl")

# Compiler
include("compiler.jl")

# VM
include("vm.jl")

# REPL
include("repl.jl")

# Transpilers
include("transpilers/transpilers.jl")

using .Transpilers.JuliaTranspiler: @monkey_julia_str
export @monkey_julia_str

# CLI
include("cli.jl")

end


text:
# using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters

# detect if using SLURM
if "SLURM_JOBID" in keys(ENV)
    SLURM_JOBID = ENV["SLURM_JOBID"]
    using ClusterManagers
else
    SLURM_JOBID = ""
end

using CountPlus, Distributed, JLD2
using Optim: minimizer, Options, BFGS, NelderMead
using ShaleDrillingLikelihood: value_function, EVobj, cost, serial_simloglik!
using SparseArrays: nonzeros
using Formatting: generate_formatter

num_nonzeros(x) = length(nonzeros(x))
intstr = generate_formatter("%'d")

# ------------------- number of simulations ----------------------

pargs = parse_commandline()
# using ArgParse: parse_args
# pargs = parse_args(["--noPar"], ShaleDrillingLikelihood.SDLParameters.arg_settings())
print_parsed_args(pargs)

DATAPATH = pargs["dataset"]

# estimation
M_cnstr = pargs["Mcnstr"]
M_full  = pargs["Mfull"]
DO_CNSTR = pargs["doCnstr"]
DO_FULL  = !pargs["noFull"]
MAXTIME_CNSTR = pargs["maxtimeCnstr"]
MAXTIME_FULL  = pargs["maxtimeFull"]
DO_PAR = !pargs["noPar"]

# parms
COMPUTE_INITIAL_VALUES = pargs["computeStarting"]
THETA0_FULL_OVERRIDE = pargs["theta"]

#payoffs
COST = pargs["cost"]
EXT = pargs["extension"]
REV = pargs["revenue"]
#
# THETA0_FULL_OVERRIDE=[-10.74244, -8.16438, -7.25647, -6.9641, -6.94902,     1.57742, -1.16717, -3.07728,    0.6, 0.355,   0.268, 0.261, 0.261, 0.148, 0.290, 0.261, 0.513, 0.638,     1.40348,             0.1, 0.6293, 1.1674, -1.6306, 0.1370,    4.2669, 0.8183, 1.3028, 1.3365, 1.0687,         -14.458, 0.09679, 0.275]
# REV = DrillingRevenue(Unconstrained(), TimeFE(2008, 2016), GathProcess() )
# COST = DrillingCost_TimeFE(2008, 2012)
REWARD = DrillReward(REV, COST, EXT)
ANTICIPATE = pargs["anticipateT1EV"]
DISCOUNT = pargs["discount"]

# approximation
PSI = PsiSpace(pargs["numPsi"])
NUM_P = pargs["numP"]
EXTEND_GRID = pargs["extendPriceGrid"]
MINP = pargs["minTransProb"]

# print some things
println("Model is\n\t$REWARD\n")
println("Firms anticipate T1EV shocks? $ANTICIPATE")
println("Psi-space is\n\t$PSI with length $(length(PSI))")

if cost(REWARD) isa DrillingCost_TimeFE
    HASRIGS = "no-rigs"
elseif cost(REWARD) isa DrillingCost_TimeFE_rigrate
    HASRIGS = "WITH-rigs"
elseif cost(REWARD) isa DrillingCost_TimeFE_costdiffs
    HASRIGS = "COSTDIFF"
else
    throw(error("don't have values for this cost fct"))
end

if "SLURM_JOBID" in keys(ENV)
    DATADIR = "/home/magerton/haynesville/intermediate_data"
else
    DATADIR = "E:/projects/haynesville/intermediate_data"
end

# --------------- create data ---------------

# load in data from disk
rdatapath = joinpath(DATADIR, DATAPATH)
println_time_flush("loading $rdatapath")
data_royalty    = DataRoyalty(REWARD, rdatapath)
data_produce    = DataProduce(REWARD, rdatapath)
data_drill_prim = DataDrillPrimitive(REWARD, rdatapath)
println_time_flush("Data imported")

thetarho0 = ThetaRho()
theta0_royalty = Theta(data_royalty; θρ = thetarho0)
theta0_produce = Theta(data_produce)
theta0_drill   = Theta(REWARD; θρ=thetarho0)

# transition matrices
zrng, ztrans = GridTransition(data_drill_prim, EXTEND_GRID, NUM_P; minp=MINP)
println_time_flush("Z range is\n\t$zrng with lengths $(length.(zrng))")
println_time_flush("Transition matrix drops probs < $MINP and has $(intstr(num_nonzeros(ztrans))) elements")

# create drilling dataset
wp = statespace(data_drill_prim)
ddm = DynamicDrillModel(REWARD, DISCOUNT, wp, zrng, ztrans, PSI, ANTICIPATE)
data_drill = DataDrill(data_drill_prim, ddm)

evdims = size(EVobj(value_function(ddm)))
println_time_flush("EV dimension is $evdims, implying $(intstr(prod(evdims))) states")

# full dataset
dataset_full = DataSetofSets(data_drill, data_royalty, data_produce, CoefLinks(data_drill))

# Solve unconstrained full model
theta0s = (theta0_drill, theta0_royalty, theta0_produce)
theta0_full = merge_thetas(theta0s, dataset_full)

leo = LocalEstObj(dataset_full, theta0_full)
reo = RemoteEstObj(leo, M_full)
ew = EstimationWrapper(leo, reo)

dograd = true

println_time_flush("Starting LL eval")
serial_simloglik!(ew, theta0_full, dograd; n=1)
println_time_flush("Done!")


text:
import Base: OneTo
import LatticeBoltzmann:
    deviatoric_tensor,
    velocity!,
    temperature,
    ConstantDensity,
    AnalyticalVelocityAndStress,
    AnalyticalEquilibrium,
    AnalyticalEquilibriumAndOffEquilibrium,
    IterativeInitializationMeiEtAl

struct InitializationTestProblem <: FluidFlowProblem
    u_max::Float64
    ν::Float64
    NX::Int64
    NY::Int64
    domain_size::Tuple{Float64, Float64}
end
InitializationTestProblem(N = 4, ν = 1.0 / 6.0) = InitializationTestProblem(
    1.0,
    ν,
    N,
    N,
    # (2.0 * pi, 2.0 * pi)
    (1.0, 1.0),
)

const U = 0.001

LatticeBoltzmann.density(
    q::Quadrature,
    problem::InitializationTestProblem,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real} = pressure(q, problem, x, y, timestep)
function LatticeBoltzmann.pressure(
    q::Quadrature,
    problem::InitializationTestProblem,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    return 1.0 -
           (1 / 4) * q.speed_of_sound_squared * problem.u_max^2 * U^2 * (cos(2x) + cos(2y))
end
LatticeBoltzmann.velocity(
    problem::InitializationTestProblem,
    x,
    y,
    timestep = 0.0,
) where {T <: Real} =
    problem.u_max * U * [
        cos(2pi * x) * sin(2pi * y)
        -sin(2pi * x) * cos(2pi * y)
    ]
function LatticeBoltzmann.velocity_gradient(
    problem::InitializationTestProblem,
    x::T,
    y::T,
    timestep::Real = 0.0,
) where {T <: Real}
    u_x = -sin(2pi * x) * sin(2pi * y)
    v_y = sin(2pi * x) * sin(2pi * y)
    u_y = cos(2pi * x) * cos(2pi * y)
    v_x = -cos(2pi * x) * cos(2pi * y)

    return problem.u_max * U * 2pi .* [u_x v_x; u_y v_y]
end

function LatticeBoltzmann.deviatoric_tensor(
    q::Quadrature,
    problem::InitializationTestProblem,
    x::T,
    y::T,
    time::Real = 0.0,
) where {T <: Real}
    a = LatticeBoltzmann.velocity_gradient(problem, x, y, time)
    ν = LatticeBoltzmann.viscosity(problem)

    σ = -ν * [
        2 * a[1, 1] a[1, 2] + a[2, 1]
        a[1, 2] + a[2, 1] 2 * a[2, 2]
    ]

    return σ * (problem.NX)
end

@testset "Initialization strategies with $q" for q in [D2Q9(), D2Q13(), D2Q17(), D2Q21()]
    ν = 1.0 / q.speed_of_sound_squared
    problem = InitializationTestProblem(1 * 8, ν)
    τ = 1.0
    range_x, range_y = range(problem)

    @testset "Velocity + constant density" begin
        strategy = ConstantDensity()
        # @show strategy
        f = initialize(strategy, q, problem)
        f_ = Array{eltype(q.weights)}(undef, size(f, 3))
        u = zeros(dimension(q))

        for x_idx in OneTo(problem.NX), y_idx in OneTo(problem.NY)
            x = range_x[x_idx]
            y = range_y[y_idx]
            f_ = f[x_idx, y_idx, :]

            ρ = density(q, f_)
            velocity!(q, f_, ρ, u)
            T = temperature(q, f_, ρ, u)
            p = pressure(q, f_, ρ, u)
            σ = deviatoric_tensor(q, τ, f_, ρ, u)
            expected_σ = deviatoric_tensor(q, problem, x, y)

            if (density(q, problem, x, y) ≉ 1.0)
                @test ρ ≉ density(q, problem, x, y)
                @test p ≉ pressure(q, problem, x, y)
            end

            @test u ≈ velocity(problem, x, y)
            @test σ ≉ deviatoric_tensor(q, problem, x, y)
        end
    end

    @testset "Velocity + pressure" begin
        strategy = AnalyticalEquilibrium()
        # @show strategy
        f = initialize(strategy, q, problem)
        f_ = Array{eltype(q.weights)}(undef, size(f, 3))
        u = zeros(dimension(q))

        for x_idx in OneTo(problem.NX), y_idx in OneTo(problem.NY)
            x = range_x[x_idx]
            y = range_y[y_idx]
            f_ = f[x_idx, y_idx, :]

            ρ = density(q, f_)
            velocity!(q, f_, ρ, u)
            T = temperature(q, f_, ρ, u)
            p = pressure(q, f_, ρ, u)
            σ = deviatoric_tensor(q, τ, f_, ρ, u)
            expected_σ = deviatoric_tensor(q, problem, x, y)

            @test ρ ≈ density(q, problem, x, y)
            @test p ≈ pressure(q, problem, x, y)

            @test u ≈ velocity(problem, x, y)
            @test σ ≉ deviatoric_tensor(q, problem, x, y)
        end
    end

    @testset "Velocity + stress initialization" begin
        strategy = AnalyticalVelocityAndStress()
        # @show strategy
        f = initialize(strategy, q, problem)

        f_ = Array{eltype(q.weights)}(undef, size(f, 3))
        u = zeros(dimension(q))

        for x_idx in OneTo(problem.NX), y_idx in OneTo(problem.NY)
            x = range_x[x_idx]
            y = range_y[y_idx]
            f_ = f[x_idx, y_idx, :]

            ρ = density(q, f_)
            velocity!(q, f_, ρ, u)
            T = temperature(q, f_, ρ, u)
            p = pressure(q, f_, ρ, u)
            σ = deviatoric_tensor(q, τ, f_, ρ, u)
            expected_σ = deviatoric_tensor(q, problem, x, y)

            # @warn deviatoric_tensor(q, problem, x, y) ./ σ

            if (density(q, problem, x, y) ≉ 1.0)
                @test ρ ≉ density(q, problem, x, y)
                @test p ≉ pressure(q, problem, x, y)
            end

            @test u ≈ velocity(problem, x, y)
            @test σ ≈ deviatoric_tensor(q, problem, x, y)
            # break;
        end
    end

    @testset "Velocity + pressure + stress" begin
        strategy = AnalyticalEquilibriumAndOffEquilibrium()
        # @show strategy
        f = initialize(strategy, q, problem)

        f_ = Array{eltype(q.weights)}(undef, size(f, 3))
        u = zeros(dimension(q))

        for x_idx in OneTo(problem.NX), y_idx in OneTo(problem.NY)
            x = range_x[x_idx]
            y = range_y[y_idx]
            f_ = f[x_idx, y_idx, :]

            ρ = density(q, f_)
            velocity!(q, f_, ρ, u)
            T = temperature(q, f_, ρ, u)
            p = pressure(q, f_, ρ, u)
            σ = deviatoric_tensor(q, τ, f_, ρ, u)
            expected_σ = deviatoric_tensor(q, problem, x, y)

            # @warn deviatoric_tensor(q, problem, x, y) ./ σ

            @test ρ ≈ density(q, problem, x, y)
            @test p ≈ pressure(q, problem, x, y)

            @test u ≈ velocity(problem, x, y)
            @test σ ≈ deviatoric_tensor(q, problem, x, y)

            # break;
        end
    end

    @testset "Mei et al" begin
        # problem = InitializationTestProblem(16, ν)
        strategy = IterativeInitializationMeiEtAl(τ, 1E-14)
        # @show strategy
        f = initialize(strategy, q, problem)

        f_ = Array{eltype(q.weights)}(undef, size(f, 3))
        u = zeros(dimension(q))

        for x_idx in OneTo(problem.NX), y_idx in OneTo(problem.NY)
            x = range_x[x_idx]
            y = range_y[y_idx]
            f_ = f[x_idx, y_idx, :]

            ρ = density(q, f_)
            velocity!(q, f_, ρ, u)
            T = temperature(q, f_, ρ, u)
            p = pressure(q, f_, ρ, u)
            σ = deviatoric_tensor(q, τ, f_, ρ, u)
            expected_σ = deviatoric_tensor(q, problem, x, y)

            @test_broken ρ ≈ density(q, problem, x, y)
            if p ≈ pressure(q, problem, x, y)
                @show (p, pressure(q, problem, x, y), q)
            else
                @test_broken p ≈ pressure(q, problem, x, y)
            end

            @test_broken u ≈ velocity(problem, x, y)
            @test_broken σ ≈ deviatoric_tensor(q, problem, x, y)
        end
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# What is included in the file `FMI2_md.jl` (model description)?
# - [Sec. 1a] the function `fmi2LoadModelDescription` to load/parse the FMI model description [exported]
# - [Sec. 1b] helper functions for the load/parse function [not exported]
# - [Sec. 2]  functions to get values from the model description in the format `fmi2Get[value](md::fmi2ModelDescription)` [exported]
# - [Sec. 3]  additional functions to get useful information from the model description in the format `fmi2Get[value](md::fmi2ModelDescription)` [exported]

using EzXML

using FMICore: fmi2ModelDescriptionModelExchange, fmi2ModelDescriptionCoSimulation, fmi2ModelDescriptionDefaultExperiment
using FMICore: fmi2ModelDescriptionReal, fmi2ModelDescriptionBoolean, fmi2ModelDescriptionInteger, fmi2ModelDescriptionString, fmi2ModelDescriptionEnumeration
using FMICore: fmi2ModelDescriptionModelStructure
using FMICore: fmi2DependencyKind
using FMICore: FMU2

######################################
# [Sec. 1a] fmi2LoadModelDescription #
######################################

"""
Extract the FMU variables and meta data from the ModelDescription

ToDo: New docstring format!
"""
function fmi2LoadModelDescription(pathToModellDescription::String)
    md = fmi2ModelDescription()

    md.stringValueReferences = Dict{String, fmi2ValueReference}()
    md.outputValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.inputValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.stateValueReferences = Array{fmi2ValueReference}(undef, 0)
    md.derivativeValueReferences = Array{fmi2ValueReference}(undef, 0)

    md.enumerations = []
    typedefinitions = nothing
    modelvariables = nothing
    modelstructure = nothing

    doc = readxml(pathToModellDescription)

    root = doc.root

    # mandatory
    md.fmiVersion = root["fmiVersion"]
    md.modelName = root["modelName"]
    md.guid = root["guid"]

    # optional
    md.generationTool           = parseNodeString(root, "generationTool"; onfail="[Unknown generation tool]")
    md.generationDateAndTime    = parseNodeString(root, "generationDateAndTime"; onfail="[Unknown generation date and time]")
    variableNamingConventionStr = parseNodeString(root, "variableNamingConvention"; onfail="flat")
    @assert (variableNamingConventionStr == "flat" || variableNamingConventionStr == "structured") ["fmi2ReadModelDescription(...): Unknown entry for `variableNamingConvention=$(variableNamingConventionStr)`."]
    md.variableNamingConvention = (variableNamingConventionStr == "flat" ? fmi2VariableNamingConventionFlat : fmi2VariableNamingConventionStructured)
    md.numberOfEventIndicators  = parseNodeInteger(root, "numberOfEventIndicators"; onfail=0)
    md.description              = parseNodeString(root, "description"; onfail="[Unknown Description]")

    # defaults
    md.modelExchange = nothing
    md.coSimulation = nothing
    md.defaultExperiment = nothing

    # additionals
    md.valueReferences = []
    md.valueReferenceIndicies = Dict{UInt, UInt}()

    for node in eachelement(root)

        if node.name == "CoSimulation" || node.name == "ModelExchange"
            if node.name == "CoSimulation"
                md.coSimulation = fmi2ModelDescriptionCoSimulation()
                md.coSimulation.modelIdentifier                        = node["modelIdentifier"]
                md.coSimulation.canHandleVariableCommunicationStepSize = parseNodeBoolean(node, "canHandleVariableCommunicationStepSize"   ; onfail=false)
                md.coSimulation.canInterpolateInputs                   = parseNodeBoolean(node, "canInterpolateInputs"                     ; onfail=false)
                md.coSimulation.maxOutputDerivativeOrder               = parseNodeInteger(node, "maxOutputDerivativeOrder"                 ; onfail=nothing)
                md.coSimulation.canGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUstate"                     ; onfail=false)
                md.coSimulation.canSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUstate"                     ; onfail=false)
                md.coSimulation.providesDirectionalDerivative          = parseNodeBoolean(node, "providesDirectionalDerivative"            ; onfail=false)
            end

            if node.name == "ModelExchange"
                md.modelExchange = fmi2ModelDescriptionModelExchange()
                md.modelExchange.modelIdentifier                        = node["modelIdentifier"]
                md.modelExchange.canGetAndSetFMUstate                   = parseNodeBoolean(node, "canGetAndSetFMUstate"                     ; onfail=false)
                md.modelExchange.canSerializeFMUstate                   = parseNodeBoolean(node, "canSerializeFMUstate"                     ; onfail=false)
                md.modelExchange.providesDirectionalDerivative          = parseNodeBoolean(node, "providesDirectionalDerivative"            ; onfail=false)
            end
        elseif node.name == "TypeDefinitions"
            md.enumerations = createEnum(node)

        elseif node.name == "ModelVariables"
            md.modelVariables = parseModelVariables(node, md)

        elseif node.name == "ModelStructure"
            md.modelStructure = fmi2ModelDescriptionModelStructure()

            for element in eachelement(node)
                if element.name == "Derivatives"
                    parseDerivatives(element, md)
                elseif element.name == "InitialUnknowns"
                    parseInitialUnknowns(element, md)
                elseif element.name == "Outputs"
                    parseOutputs(element, md)
                else
                    @warn "Unknown tag `$(element.name)` for node `ModelStructure`."
                end
            end

        elseif node.name == "DefaultExperiment"
            md.defaultExperiment = fmi2ModelDescriptionDefaultExperiment()
            md.defaultExperiment.startTime  = parseNodeReal(node, "startTime")
            md.defaultExperiment.stopTime   = parseNodeReal(node, "stopTime")
            md.defaultExperiment.tolerance  = parseNodeReal(node, "tolerance")
            md.defaultExperiment.stepSize   = parseNodeReal(node, "stepSize")
        end
    end

    # creating an index for value references (fast look-up for dependencies)
    for i in 1:length(md.valueReferences)
        md.valueReferenceIndicies[md.valueReferences[i]] = i
    end

    md
end

#######################################
# [Sec. 1b] helpers for load function #
#######################################

# Returns the indices of the state derivatives.
function getDerivativeIndices(node::EzXML.Node)
    indices = []
    for element in eachelement(node)
        if element.name == "Derivatives"
            for derivative in eachelement(element)
                ind = parse(Int, derivative["index"])
                der = nothing
                derKind = nothing

                if haskey(derivative, "dependencies")
                    der = split(derivative["dependencies"], " ")

                    if der[1] == ""
                        der = fmi2Integer[]
                    else
                        der = collect(parse(fmi2Integer, e) for e in der)
                    end
                end

                if haskey(derivative, "dependenciesKind")
                    derKind = split(derivative["dependenciesKind"], " ")
                end

                push!(indices, (ind, der, derKind))
            end
        end
    end
    sort!(indices, rev=true)
end

# Parses the model variables of the FMU model description.
function parseModelVariables(nodes::EzXML.Node, md::fmi2ModelDescription)
    numberOfVariables = 0
    for node in eachelement(nodes)
        numberOfVariables += 1
    end
    scalarVariables = Array{fmi2ScalarVariable}(undef, numberOfVariables)

    index = 1
    for node in eachelement(nodes)
        name = node["name"]
        valueReference = parse(fmi2ValueReference, node["valueReference"])

        scalarVariables[index] = fmi2ScalarVariable(name, valueReference)

        if !(valueReference in md.valueReferences)
            push!(md.valueReferences, valueReference)
        end

        if haskey(node, "description")
            scalarVariables[index].description = node["description"]
        end
        if haskey(node, "causality")
            scalarVariables[index].causality = fmi2StringToCausality(node["causality"])

            if scalarVariables[index].causality == fmi2CausalityOutput
                push!(md.outputValueReferences, valueReference)
            elseif scalarVariables[index].causality == fmi2CausalityInput
                push!(md.inputValueReferences, valueReference)
            end
        end
        if haskey(node, "variability")
            scalarVariables[index].variability = fmi2StringToVariability(node["variability"])
        end
        if haskey(node, "initial")
            scalarVariables[index].initial = fmi2StringToInitial(node["initial"])
        end

        # type node
        typenode = nothing
        typename = node.firstelement.name

        if typename == "Real"
            scalarVariables[index]._Real = fmi2ModelDescriptionReal()
            typenode = scalarVariables[index]._Real
            if haskey(node.firstelement, "quantity")
                typenode.quantity = node.firstelement["quantity"]
            end
            if haskey(node.firstelement, "unit")
                typenode.unit = node.firstelement["unit"]
            end
            if haskey(node.firstelement, "displayUnit")
                typenode.displayUnit = node.firstelement["displayUnit"]
            end
            if haskey(node.firstelement, "relativeQuantity")
                typenode.relativeQuantity = parseBoolean(node.firstelement["relativeQuantity"])
            end
            if haskey(node.firstelement, "min")
                typenode.min = parseReal(node.firstelement["min"])
            end
            if haskey(node.firstelement, "max")
                typenode.max = parseReal(node.firstelement["max"])
            end
            if haskey(node.firstelement, "nominal")
                typenode.nominal = parseReal(node.firstelement["nominal"])
            end
            if haskey(node.firstelement, "unbounded")
                typenode.unbounded = parseBoolean(node.firstelement["unbounded"])
            end
            if haskey(node.firstelement, "start")
                typenode.start = parseReal(node.firstelement["start"])
            end
            if haskey(node.firstelement, "derivative")
                typenode.derivative = parse(UInt, node.firstelement["derivative"])
            end
        elseif typename == "String"
            scalarVariables[index]._String = fmi2ModelDescriptionString()
            typenode = scalarVariables[index]._String
            if haskey(node.firstelement, "start")
                scalarVariables[index]._String.start = node.firstelement["start"]
            end
            # ToDo: remaining attributes
        elseif typename == "Boolean"
            scalarVariables[index]._Boolean = fmi2ModelDescriptionBoolean()
            typenode = scalarVariables[index]._Boolean
            if haskey(node.firstelement, "start")
                scalarVariables[index]._Boolean.start = parseFMI2Boolean(node.firstelement["start"])
            end
            # ToDo: remaining attributes
        elseif typename == "Integer"
            scalarVariables[index]._Integer = fmi2ModelDescriptionInteger()
            typenode = scalarVariables[index]._Integer
            if haskey(node.firstelement, "start")
                scalarVariables[index]._Integer.start = parseInteger(node.firstelement["start"])
            end
            # ToDo: remaining attributes
        elseif typename == "Enumeration"
            scalarVariables[index]._Enumeration = fmi2ModelDescriptionEnumeration()
            typenode = scalarVariables[index]._Enumeration
            # ToDo: Save start value
            # ToDo: remaining attributes
        else
            @warn "Unknown data type `$(typename)`."
        end

        # generic attributes
        if typenode != nothing
            if haskey(node.firstelement, "declaredType")
                typenode.declaredType = node.firstelement["declaredType"]
            end
        end

        md.stringValueReferences[name] = valueReference

        index += 1
    end

    scalarVariables
end

# Parses the `ModelStructure.Derivatives` of the FMU model description.
function parseUnknwon(node::EzXML.Node)
    if haskey(node, "index")
        varDep = fmi2VariableDependency(parseInteger(node["index"]))

        if haskey(node, "dependencies")
            dependencies = node["dependencies"]
            if length(dependencies) > 0
                dependenciesSplit = split(dependencies, " ")
                if length(dependenciesSplit) > 0
                    varDep.dependencies = collect(parse(UInt, e) for e in dependenciesSplit)
                end
            end
        end

        if haskey(node, "dependenciesKind")
            dependenciesKind = node["dependenciesKind"]
            if length(dependenciesKind) > 0
                dependenciesKindSplit = split(dependenciesKind, " ")
                if length(dependenciesKindSplit) > 0
                    varDep.dependenciesKind = collect(fmi2StringToDependencyKind(e) for e in dependenciesKindSplit)
                end
            end
        end

        if varDep.dependencies != nothing && varDep.dependenciesKind != nothing
            if length(varDep.dependencies) != length(varDep.dependenciesKind)
                @warn "Length of field dependencies ($(length(varDep.dependencies))) doesn't match length of dependenciesKind ($(length(varDep.dependenciesKind)))."
            end
        end

        return varDep
    else
        return nothing
    end
end

# ToDo: Comment
function parseDerivatives(nodes::EzXML.Node, md::fmi2ModelDescription)
    @assert (nodes.name == "Derivatives") "Wrong element name."
    md.modelStructure.derivatives = []
    for node in eachelement(nodes)
        if node.name == "Unknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                # find states and derivatives
                derSV = md.modelVariables[varDep.index]
                derVR = derSV.valueReference
                stateVR = md.modelVariables[derSV._Real.derivative].valueReference

                if stateVR ∉ md.stateValueReferences
                    push!(md.stateValueReferences, stateVR)
                end
                if derVR ∉ md.derivativeValueReferences
                    push!(md.derivativeValueReferences, derVR)
                end

                push!(md.modelStructure.derivatives, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        else
            @warn "Unknown entry in `ModelStructure.Derivatives` named `$(node.name)`."
        end
    end
end

# ToDo: Comment
function parseInitialUnknowns(nodes::EzXML.Node, md::fmi2ModelDescription)
    @assert (nodes.name == "InitialUnknowns") "Wrong element name."
    md.modelStructure.initialUnknowns = []
    for node in eachelement(nodes)
        if node.name == "Unknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                push!(md.modelStructure.initialUnknowns, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        else
            @warn "Unknown entry in `ModelStructure.InitialUnknowns` named `$(node.name)`."
        end
    end
end

# ToDo: Comment
function parseOutputs(nodes::EzXML.Node, md::fmi2ModelDescription)
    @assert (nodes.name == "Outputs") "Wrong element name."
    md.modelStructure.outputs = []
    for node in eachelement(nodes)
        if node.name == "Unknown"
            if haskey(node, "index")
                varDep = parseUnknwon(node)

                # find outputs
                outVR = md.modelVariables[varDep.index].valueReference

                if outVR ∉ md.outputValueReferences
                    push!(md.outputValueReferences, outVR)
                end

                push!(md.modelStructure.outputs, varDep)
            else
                @warn "Invalid entry for node `Unknown` in `ModelStructure`, missing entry `index`."
            end
        else
            @warn "Unknown entry in `ModelStructure.Outputs` named `$(node.name)`."
        end
    end
end

# Parses a Bool value represented by a string.
function parseBoolean(s::Union{String, SubString{String}}; onfail=nothing)
    if s == "true"
        return true
    elseif s == "false"
        return false
    else
        @assert onfail != nothing ["parseBoolean(...) unknown boolean value '$s'."]
        return onfail
    end
end

# parses node (interpreted as boolean)
function parseNodeBoolean(node, key; onfail=nothing)
    if haskey(node, key)
        return parseBoolean(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# Parses an Integer value represented by a string.
function parseInteger(s::Union{String, SubString{String}}; onfail=nothing)
    if onfail == nothing
        return parse(fmi2Integer, s)
    else
        try
            return parse(fmi2Integer, s)
        catch
            return onfail
        end
    end
end

# parses node (interpreted as integer)
function parseNodeInteger(node, key; onfail=nothing)
    if haskey(node, key)
        return parseInteger(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# Parses a real value represented by a string.
function parseReal(s::Union{String, SubString{String}}; onfail=nothing)
    if onfail == nothing
        return parse(fmi2Real, s)
    else
        try
            return parse(fmi2Real, s)
        catch
            return onfail
        end
    end
end

# parses node (interpreted as real)
function parseNodeReal(node, key; onfail=nothing)
    if haskey(node, key)
        return parseReal(node[key]; onfail=onfail)
    else
        return onfail
    end
end

# parses node (interpreted as string)
function parseNodeString(node, key; onfail=nothing)
    if haskey(node, key)
        return node[key]
    else
        return onfail
    end
end

# Parses a fmi2Boolean value represented by a string.
function parseFMI2Boolean(s::Union{String, SubString{String}})
    if parseBoolean(s)
        return fmi2True
    else
        return fmi2False
    end
end

# set the datatype and attributes of an model variable
function fmi2SetDatatypeVariables(node::EzXML.Node, md::fmi2ModelDescription, sv)
    type = fmi2DatatypeVariable()
    typenode = node.firstelement
    typename = typenode.name
    type.start = nothing
    type.min = nothing
    type.max = nothing
    type.quantity = nothing
    type.unit = nothing
    type.displayUnit = nothing
    type.relativeQuantity = nothing
    type.nominal = nothing
    type.unbounded = nothing
    type.derivative = nothing
    type.reinit = nothing
    type.datatype = nothing

    if typename == "Real"
        type.datatype = fmi2Real
        sv._Real = fmi2ModelDescriptionReal()
    elseif typename == "String"
        type.datatype = fmi2String
    elseif typename == "Boolean"
        type.datatype = fmi2Boolean
    elseif typename == "Integer"
        type.datatype = fmi2Integer
    elseif typename == "Enumeration"
        type.datatype = fmi2Enum
    else
        @warn "Unknown data type `$(type.datatype)`."
    end

    if haskey(typenode, "declaredType")
        type.declaredType = typenode["declaredType"]
    end

    if haskey(typenode, "start")
        if typename == "Real"
            type.start = parse(fmi2Real, typenode["start"])
            sv._Real.start = type.start
        elseif typename == "Integer"
            type.start = parse(fmi2Integer, typenode["start"])
        elseif typename == "Boolean"
            type.start = parseFMI2Boolean(typenode["start"])
        elseif typename == "Enumeration"
            for i in 1:length(md.enumerations)
                if type.declaredType == md.enumerations[i][1] # identify the enum by the name
                    type.start = md.enumerations[i][1 + parse(fmi2Integer, typenode["start"])] # find the enum value and set it
                end
            end
        elseif typename == "String"
            type.start = typenode["start"]
        else
            @warn "setDatatypeVariables(...) unimplemented start value type $typename"
            type.start = typenode["start"]
        end
    end

    if haskey(typenode, "min") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        if type.datatype == fmi2Real
            type.min = parse(fmi2Real, typenode["min"])
        else
            type.min = parse(fmi2Integer, typenode["min"])
        end
    end
    if haskey(typenode, "max") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        if type.datatype == fmi2Real
            type.max = parse(fmi2Real, typenode["max"])
        elseif type.datatype == fmi2Integer
            type.max = parse(fmi2Integer, typenode["max"])
        end
    end
    if haskey(typenode, "quantity") && (type.datatype == fmi2Real || type.datatype == fmi2Integer || type.datatype == fmi2Enum)
        type.quantity = typenode["quantity"]
    end
    if haskey(typenode, "unit") && type.datatype == fmi2Real
        type.unit = typenode["unit"]
    end
    if haskey(typenode, "displayUnit") && type.datatype == fmi2Real
        type.displayUnit = typenode["displayUnit"]
    end
    if haskey(typenode, "relativeQuantity") && type.datatype == fmi2Real
        type.relativeQuantity = convert(fmi2Boolean, parse(Bool, typenode["relativeQuantity"]))
    end
    if haskey(typenode, "nominal") && type.datatype == fmi2Real
        type.nominal = parse(fmi2Real, typenode["nominal"])
    end
    if haskey(typenode, "unbounded") && type.datatype == fmi2Real
        type.unbounded = parse(fmi2Boolean, typenode["unbounded"])
    end
    if haskey(typenode, "derivative") && type.datatype == fmi2Real
        type.derivative = parse(fmi2Integer, typenode["derivative"])
        if typename == "Real"
            sv._Real.derivative = type.derivative
        end
    end
    if haskey(typenode, "reinit") && type.datatype == fmi2Real
        type.reinit = parseFMI2Boolean(typenode["reinit"])
    end
    type
end

#=
Read all enumerations from the modeldescription and store them in a matrix. First entries are the enum names
-------------------------------------------
Example:
"enum1name" "value1"    "value2"
"enum2name" "value1"    "value2"
=#
function createEnum(node::EzXML.Node)
    enum = 1
    idx = 1
    enumerations = []
    for simpleType in eachelement(node)
        name = simpleType["name"]
        for type in eachelement(simpleType)
            if type.name == "Enumeration"
                enum = []
                push!(enum, name)
                for item in eachelement(type)
                    push!(enum, item["name"])
                end
                push!(enumerations, enum)
            end
        end
    end
    enumerations
end

################################
# [Sec. 2] get value functions #
################################

"""
Returns startTime from DefaultExperiment if defined else defaults to nothing.

    ToDo: update docstring format.
"""
function fmi2GetDefaultStartTime(md::fmi2ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.startTime
end

"""
Returns stopTime from DefaultExperiment if defined else defaults to nothing.

    ToDo: update docstring format.
"""
function fmi2GetDefaultStopTime(md::fmi2ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.stopTime
end

"""
Returns tolerance from DefaultExperiment if defined else defaults to nothing.

ToDo: update docstring format.
"""
function fmi2GetDefaultTolerance(md::fmi2ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.tolerance
end

"""
Returns stepSize from DefaultExperiment if defined else defaults to nothing.

ToDo: update docstring format.
"""
function fmi2GetDefaultStepSize(md::fmi2ModelDescription)
    if md.defaultExperiment == nothing
        return nothing
    end
    return md.defaultExperiment.stepSize
end

"""
Returns the tag 'modelName' from the model description.

ToDo: update docstring format.
"""
function fmi2GetModelName(md::fmi2ModelDescription)#, escape::Bool = true)
    md.modelName
end

"""
Returns the tag 'guid' from the model description.

ToDo: update docstring format.
"""
function fmi2GetGUID(md::fmi2ModelDescription)
    md.guid
end

"""
Returns the tag 'generationtool' from the model description.

ToDo: update docstring format.
"""
function fmi2GetGenerationTool(md::fmi2ModelDescription)
    md.generationTool
end

"""
Returns the tag 'generationdateandtime' from the model description.

ToDo: update docstring format.
"""
function fmi2GetGenerationDateAndTime(md::fmi2ModelDescription)
    md.generationDateAndTime
end

"""
Returns the tag 'varaiblenamingconvention' from the model description.

ToDo: update docstring format.
"""
function fmi2GetVariableNamingConvention(md::fmi2ModelDescription)
    md.variableNamingConvention
end

"""
Returns the tag 'numberOfEventIndicators' from the model description.

ToDo: update docstring format.
"""
function fmi2GetNumberOfEventIndicators(md::fmi2ModelDescription)
    md.numberOfEventIndicators
end

"""
Returns the number of states of the FMU.

ToDo: update docstring format.
"""
function fmi2GetNumberOfStates(md::fmi2ModelDescription)
    length(md.stateValueReferences)
end

"""
Returns true, if the FMU supports co simulation

ToDo: update docstring format.
"""
function fmi2IsCoSimulation(md::fmi2ModelDescription)
    return (md.coSimulation != nothing)
end

"""
Returns true, if the FMU supports model exchange

ToDo: update docstring format.
"""
function fmi2IsModelExchange(md::fmi2ModelDescription)
    return (md.modelExchange != nothing)
end

##################################
# [Sec. 3] information functions #
##################################

"""
Returns if the FMU model description contains `dependency` information.

ToDo: update docstring format.
"""
function fmi2DependenciesSupported(md::fmi2ModelDescription)
    if md.modelStructure === nothing
        return false
    end

    return true
end

"""
Returns if the FMU model description contains `dependency` information for `derivatives`.

ToDo: update docstring format.
"""
function fmi2DerivativeDependenciesSupported(md::fmi2ModelDescription)
    if !fmi2DependenciesSupported(md)
        return false
    end

    der = md.modelStructure.derivatives
    if der === nothing || length(der) <= 0
        return false
    end

    return true
end

"""
Returns the tag 'modelIdentifier' from CS or ME section.

ToDo: update docstring format.
"""
function fmi2GetModelIdentifier(md::fmi2ModelDescription; type=nothing)

    if type === nothing
        if fmi2IsCoSimulation(md)
            return md.coSimulation.modelIdentifier
        elseif fmi2IsModelExchange(md)
            return md.modelExchange.modelIdentifier
        else
            @assert false "fmi2GetModelName(...): FMU does not support ME or CS!"
        end
    elseif type == fmi2TypeCoSimulation
        return md.coSimulation.modelIdentifier
    elseif type == fmi2TypeModelExchange
        return md.modelExchange.modelIdentifier
    end
end

"""
Returns true, if the FMU supports the getting/setting of states

ToDo: update docstring format.
"""
function fmi2CanGetSetState(md::fmi2ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.canGetAndSetFMUstate) || (md.modelExchange != nothing && md.modelExchange.canGetAndSetFMUstate)
end

"""
Returns true, if the FMU state can be serialized

ToDo: update docstring format.
"""
function fmi2CanSerializeFMUstate(md::fmi2ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.canSerializeFMUstate) || (md.modelExchange != nothing && md.modelExchange.canSerializeFMUstate)
end

"""
Returns true, if the FMU provides directional derivatives

ToDo: update docstring format.
"""
function fmi2ProvidesDirectionalDerivative(md::fmi2ModelDescription)
    return (md.coSimulation != nothing && md.coSimulation.providesDirectionalDerivative) || (md.modelExchange != nothing && md.modelExchange.providesDirectionalDerivative)
end

"""
Returns a dictionary `Dict(fmi2ValueReference, Array{String})` of value references and their corresponding names

ToDo: update docstring format.
"""
function fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.valueReferences)
    dict = Dict{fmi2ValueReference, Array{String}}()
    for vr in vrs
        dict[vr] = fmi2ValueReferenceToString(md, vr)
    end
    return dict
end

function fmi2GetValueReferencesAndNames(fmu::FMU2)
    fmi2GetValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns a array of names corresponding to value references `vrs`

If there are multiple names per value reference, availabel modes are `:first` (default, pick only the first one), `:group` (pick all and group them into an array) and `:flat` (pick all, but flat them out into a 1D-array together with all other names)

ToDo: update docstring format.
"""
function fmi2GetNames(md::fmi2ModelDescription; vrs=md.valueReferences, mode=:first)
    names = []
    for vr in vrs
        ns = fmi2ValueReferenceToString(md, vr)

        if mode == :first
            push!(names, ns[1])
        elseif mode == :group
            push!(names, ns)
        elseif mode == :flat
            for n in ns
                push!(names, n)
            end
        else
            @assert false "fmi2GetNames(...) unknown mode `mode`, please choose between `:first`, `:group` and `:flat`."
        end
    end
    return names
end

function fmi2GetNames(fmu::FMU2; kwargs...)
    fmi2GetNames(fmu.modelDescription; kwargs...)
end

"""
Returns a dict with (vrs, names of inputs)

ToDo: update docstring format.
"""
function fmi2GetInputValueReferencesAndNames(md::fmi2ModelDescription)
    fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.inputValueReferences)
end

function fmi2GetInputValueReferencesAndNames(fmu::FMU2)
    fmi2GetInputValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns names of inputs

ToDo: update docstring format.
"""
function fmi2GetInputNames(md::fmi2ModelDescription; kwargs...)
    fmi2GetNames(md; vrs=md.inputValueReferences, kwargs...)
end

function fmi2GetInputNames(fmu::FMU2; kwargs...)
    fmi2GetInputNames(fmu.modelDescription; kwargs...)
end

"""
ToDo: update docstring format.
"""
function fmi2GetOutputValueReferencesAndNames(md::fmi2ModelDescription)
    fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.outputValueReferences)
end

function fmi2GetOutputValueReferencesAndNames(fmu::FMU2)
    fmi2GetOutputValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns names of outputs

ToDo: update docstring format.
"""
function fmi2GetOutputNames(md::fmi2ModelDescription; kwargs...)
    fmi2GetNames(md; vrs=md.outputValueReferences, kwargs...)
end

function fmi2GetOutputNames(fmu::FMU2; kwargs...)
    fmi2GetOutputNames(fmu.modelDescription; kwargs...)
end

"""
ToDo: update docstring format.
"""
function fmi2GetParameterValueReferencesAndNames(md::fmi2ModelDescription)
    fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.parameterValueReferences)
end

function fmi2GetParameterValueReferencesAndNames(fmu::FMU2)
    fmi2GetParameterValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns names of parameters

ToDo: update docstring format.
"""
function fmi2GetParameterNames(md::fmi2ModelDescription; kwargs...)
    fmi2GetNames(md; vrs=md.parameterValueReferences, kwargs...)
end

function fmi2GetParameterNames(fmu::FMU2; kwargs...)
    fmi2GetParameterNames(fmu.modelDescription; kwargs...)
end

"""
Returns dict(vrs, names of states)

ToDo: update docstring format.
"""
function fmi2GetStateValueReferencesAndNames(md::fmi2ModelDescription)
    fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.stateValueReferences)
end

function fmi2GetStateValueReferencesAndNames(fmu::FMU2)
    fmi2GetStateValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns names of states

ToDo: update docstring format.
"""
function fmi2GetStateNames(md::fmi2ModelDescription; kwargs...)
    fmi2GetNames(md; vrs=md.stateValueReferences, kwargs...)
end

function fmi2GetStateNames(fmu::FMU2; kwargs...)
    fmi2GetStateNames(fmu.modelDescription; kwargs...)
end

"""
ToDo: update docstring format.
"""
function fmi2GetDerivateValueReferencesAndNames(md::fmi2ModelDescription)
    fmi2GetValueReferencesAndNames(md::fmi2ModelDescription; vrs=md.derivativeValueReferences)
end

function fmi2GetDerivateValueReferencesAndNames(fmu::FMU2)
    fmi2GetDerivateValueReferencesAndNames(fmu.modelDescription)
end

"""
Returns names of derivatives

ToDo: update docstring format.
"""
function fmi2GetDerivativeNames(md::fmi2ModelDescription; kwargs...)
    fmi2GetNames(md; vrs=md.derivativeValueReferences, kwargs...)
end

function fmi2GetDerivativeNames(fmu::FMU2; kwargs...)
    fmi2GetDerivativeNames(fmu.modelDescription; kwargs...)
end

"""
Returns a dictionary of variables with their descriptions

ToDo: update docstring format.
"""
function fmi2GetNamesAndDescriptions(md::fmi2ModelDescription)
    Dict(md.modelVariables[i].name => md.modelVariables[i].description for i = 1:length(md.modelVariables))
end

function fmi2GetNamesAndDescriptions(fmu::FMU2)
    fmi2GetNamesAndDescriptions(fmu.modelDescription)
end

"""
Returns a dictionary of variables with their units

ToDo: update docstring format.
"""
function fmi2GetNamesAndUnits(md::fmi2ModelDescription)
    Dict(md.modelVariables[i].name => fmi2GetUnit(md.modelVariables[i]) for i = 1:length(md.modelVariables))
end

function fmi2GetNamesAndUnits(fmu::FMU2)
    fmi2GetNamesAndUnits(fmu.modelDescription)
end

"""
Returns a dictionary of variables with their initials

ToDo: update docstring format.
"""
function fmi2GetNamesAndInitials(md::fmi2ModelDescription)
    Dict(md.modelVariables[i].name => fmi2GetInitial(md.modelVariables[i]) for i = 1:length(md.modelVariables))
end

function fmi2GetNamesAndInitials(fmu::FMU2)
    fmi2GetNamesAndInitials(fmu.modelDescription)
end

"""
Returns a dictionary of variables with their starting values

ToDo: update docstring format.
"""
function fmi2GetNamesAndInitials(md::fmi2ModelDescription)
    Dict(md.modelVariables[i].name => md.modelVariables[i].initial for i = 1:length(md.modelVariables))
end

function fmi2GetNamesAndInitials(fmu::FMU2)
    fmi2GetNamesAndInitials(fmu.modelDescription)
end


text:
# Comments:
# ALS without orthog is really slow (probably getUTU needs to be optimized)
# and almost never used unless an initial tt is inputted which is not site-k

# no initial tt, automatically with orthogonalization
function TT_ALS(tensor::Array{Float64},rnks::Vector{Int64})
    N     = ndims(tensor);
    sizes = size(tensor);
    cores = Vector{Array{Float64,3}}(undef,N);
    for i = 1:N-1 # creating site-N canonical initial tensor train
        tmp = qr(rand(rnks[i]*sizes[i], rnks[i+1]));
        cores[i] = reshape(Matrix(tmp.Q),(rnks[i], sizes[i], rnks[i+1]));
    end
    cores[N] = reshape(rand(rnks[N]*sizes[N]),(rnks[N], sizes[N], 1));
    tt0 = MPT(cores,N);
    return TT_ALS(tensor,tt0)
end

# with / without orthogonalization
function TT_ALS(tensor::Array{Float64},tt::MPS)
    maxiter = 10;
    N       = order(tt);
    rnks    = rank(tt);
    sizes   = size(tt);

    for i = 1:maxiter
        for k = 1:2N-2
            if tt.normcore == 0
                swipe = [collect(1:N)..., collect(N-1:-1:2)...];
                n     = swipe[k];
                UTU   = getUTU(tt,n);
                UTy   = getUTy(tt,tensor,n);
                tt[n] = reshape(inv(UTU)*UTy,(rnks[n][1],sizes[n][1],rnks[n][2]));
            else
                swipe = [collect(N:-1:2)..., collect(1:N-1)...];
                Dir   = Int.([-ones(1,N-1)...,ones(1,N-1)...]);
                n     = swipe[k];
                UTy   = getUTy(tt,tensor,n);
                tt[n] = reshape(UTy,(rnks[n][1],sizes[n][1],rnks[n][2]));
                shiftMPTnorm(tt,n,Dir[k])
            end
        end
    end
    return tt
end

# TT-ALS for a vector without initial tt
function TT_ALS(vector::Vector{Float64},middlesizes::Matrix{Int64},rnks::Vector{Int64})
    tensor = reshape(vector,Tuple(middlesizes));
    return TT_ALS(tensor,rnks);
end

# TT-ALS for vector with initial TT
function TT_ALS(vector::Vector{Float64},tt0::MPT)
    tensor = reshape(vector,Tuple([size(tt0)[i][1] for i = 1:order(tt0)]));
    return TT_ALS(tensor,tt0);
end

##########################################
# functions for ALS with orthogonalization
function getUTy(tt::MPS,tensor::Array{Float64},n::Int64)
    N     = order(tt);
    sizes = size(tensor);
    rnks  = rank(tt);
    if n == N
        Gleft    = supercores(tt,N);
        newsizes = (prod(sizes[1:N-1]), sizes[N]);
        UTy      = Gleft*reshape(tensor,Tuple(newsizes));
    elseif n == 1
        Gright   = supercores(tt,1);
        newsizes = (sizes[1], prod(sizes[2:N]));
        UTy      = reshape(tensor,Tuple(newsizes))*Gright;
    else
        Gleft, Gright = supercores(tt,n);
        newsizes1     = (prod(sizes[1:n-1]), prod(sizes[n:N]));
        tmp           = Gleft*reshape(tensor,newsizes1);
        newsizes2     = (rnks[n][1]*sizes[n], prod(sizes[n+1:N]));
        UTy           = reshape(tmp,newsizes2)*Gright;
    end
    return UTy[:]
end

function supercores(tt::MPT{3},n::Int64)
    N     = order(tt);
    sizes = size(tt);
    rnks  = rank(tt);
    if  n == 1
        Gright = tt[2];
        sright = sizes[2][1];
        for i = 3:N
            Gright = Gright*tt[i];
            sright = sright*sizes[i][1];
        end
        return Matrix(reshape(Gright,(rnks[1][2], sright))')
    elseif n == N
        Gleft = tt[1];
        sright = sizes[1][1];
        for i = 2:N-1
            Gleft = Gleft*tt[i];
            sright = sright*sizes[i][1];
        end
        return Matrix(reshape(Gleft,(sright, rnks[N][1]))')
    else
        Gleft = tt[1];
        sleft = sizes[1][1];
        for i = 2:n-1
            Gleft = Gleft*tt[i];
            sleft = sleft*sizes[i][1]
        end
        Gright = tt[n+1];
        sright = sizes[n+1][1];
        for i = n+2:N
            Gright = Gright*tt[i];
            sright = sright*sizes[i][1];
        end
        return Matrix(reshape(Gleft,(sleft, rnks[n][1]))'), Matrix(reshape(Gright,(rnks[n][2],sright))')
    end
end

function shiftMPTnorm(mpt::MPT,n::Int64,dir::Int64)

        if dir == 1
            Gl = unfold(mpt[n],[ndims(mpt[n])],"right");
            ind   = 1;
            F = qr(Gl);
            R = Matrix(F.R); Q = Matrix(F.Q);
        elseif dir == -1
            Gr = unfold(mpt[n],[1],"left");
            ind  = 3;
            F = qr(Gr');
            R = Matrix(F.R); Qt = Matrix(F.Q);
            Q = Qt';
        end

        mpt[n]       = reshape(Q, size(mpt[n]));
        mpt[n+dir]   = nmodeproduct(R,mpt[n+dir],ind);
        mpt.normcore = mpt.normcore + dir;
end

# function for ALS without orthogonalization
function getUTU(tt::MPS,n::Int64)
    N     = order(tt);
    sizes = size(tt);
    rnks  = rank(tt);

    Gleft = [1];
    for i = 1:n-1
        Gleft = Gleft * contractcores(tt[i],tt[i]);
    end
    Gleft = reshape(Gleft,(rnks[n][1],rnks[n][1]));

    Gright = [1];
    for i = N:-1:n+1
        Gright = contractcores(tt[i],tt[i]) * Gright;
    end
    Gright = reshape(Gright,(rnks[n][2],rnks[n][2]));

    return kron(kron(Gright, 1.0*Matrix(I,sizes[n][1],sizes[n][1])), Gleft)
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL; GL = ViewerGL

# input of primitive shapes
W,(WW,EW,FW) = Lar.cuboid([1,1], true)
square = (W,EW)
# hierarchical assembly
model = Lar.Struct([
            Lar.Struct([
                square,
                Lar.t(.175,.175), Lar.s(.25,.25),
                square
            ]),
            Lar.t(.825,.825), Lar.s(-.25,-.25),
            square
        ])
# visualization of generated wire-frame model
V,EV = Lar.struct2lar(model)
VV = [[k] for k=1:size(V,2)]
GL.VIEW(GL.numbering(.4)((V,[VV, EV])));

triangles = Lar.triangulate2d(V, EV)
GL.VIEW([ GL.GLGrid(V,triangles) ]);

# generate edges
ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
# remove duplicated edges from triangulation
nodupev = collect(Set(ev))
GL.VIEW(GL.numbering(.25)((V,[VV, nodupev, triangles]),GL.COLORS[1],0.2));