text:
################################################################################
#
#          NfOrd/Ideal.jl : Ideals of orders of number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export NfOrdGenIdl

export deepcopy, parent, order, basis, basis_mat, basis_mat_inv, minimum, norm,
       ==, in, +, *, intersection, lcm, idempotents, mod, pradical

################################################################################
#
#  Hashing
#
################################################################################

# a (bad) hash function
# - slow (due to basis)
# - unless basis matrix is in HNF it is also non-unique
function hash(A::NfOrdIdl)
  return hash(basis_mat(A))
end

################################################################################
#
#  Deepcopy
#
################################################################################

doc"""
***
    deepcopy(I::NfOrdIdl) -> NfOrdIdl

> Returns a copy of the ideal $I$.
"""
function deepcopy(A::NfOrdIdl)
  B = typeof(A)(order(A))
  for i in fieldnames(A)
    if i == :parent
      continue
    end
    if isdefined(A, i)
      if i == :basis
        setfield!(B, i, NfOrdElem{typeof(parent(A))}[ deepcopy(x) for x in A.basis])
      else
        setfield!(B, i, deepcopy(getfield(A, i)))
      end
    end
  end
  return B
end

################################################################################
#
#  Parent
#
################################################################################

parent(a::NfOrdIdl) = a.parent

################################################################################
#
#  Order
#
################################################################################

order(a::NfOrdIdlSet) = a.order

doc"""
***
    parent(I::NfOrdIdl) -> NfOrd

> Returns the order of $I$.
"""
order(a::NfOrdIdl) = order(parent(a))

################################################################################
#
#  Principal ideal creation
#
################################################################################

doc"""
    *(O::NfOrd, x::NfOrdElem) -> NfOrdIdl
    *(x::NfOrdElem, O::NfOrd) -> NfOrdIdl

> Returns the principal ideal $(x)$ of $\mathcal O$.
"""
function *(O::NfOrd, x::NfOrdElem)
  parent(x) != O && error("Order of element does not coincide with order")
  return ideal(O, x)
end

*(x::NfOrdElem, O::NfOrd) = O*x

################################################################################
#
#  Basis
#
################################################################################

doc"""
***
    basis(I::NfOrdIdl) -> Array{NfOrdElem, 1}

> Returns the $\mathbf Z$-basis of $I$.
"""
function basis(a::NfOrdIdl)
  O = order(a)
  if isdefined(a, :basis)
    return a.basis
  end
  B = Array(elem_type(O), degree(order(a)))
  for i in 1:degree(O)
    t = zero(O)
    v = Array(fmpz, degree(O))
    for j in 1:degree(O)
      t += basis(O)[j]*basis_mat(a)[i,j]
      #v[j] = basis_mat(a)[i,j]
    end
    B[i] = t
    #B[i].elem_in_basis = v
  end
  a.basis = B
  return elem_type(O)[ deepcopy(b) for b in a.basis]
end

################################################################################
#
#  (Inverse) basis matrix
#
################################################################################

doc"""
***
    basis_mat(I::NfOrdIdl) -> fmpz_mat

> Returns the basis matrix of $I$ with respect to the basis of the order.
"""
function basis_mat(a::NfOrdIdl)
  return deepcopy(a.basis_mat)
end

doc"""
***
    basis_mat_inv(I::NfOrdIdl) -> FakeFmpqMat

> Returns the inverse of the basis matrix of $I$.
"""
function basis_mat_inv(a::NfOrdGenIdl)
  if isdefined(a, :basis_mat_inv)
    return deepcopy(a.basis_mat_inv)
  else
    m,d = pseudo_inv(a.basis_mat)
    a.basis_mat_inv = FakeFmpqMat(m,d)
    return deepcopy(a.basis_mat_inv)
  end
end

################################################################################
#
#  Minimum
#
################################################################################

doc"""
***
    minimum(I::NfOrdIdl) -> fmpz

> Returns the minimum of $I$, that is, the minimum of
> $I \cap \mathbf Z_{\geq 0}$, where $\mathcal O$ is the order of $I$.
"""
function minimum(a::NfOrdIdl)
  if isdefined(a, :minimum)
    return deepcopy(a.minimum)
  else
    @hassert :NfOrd 2 isone(basis(order(a)))
    a.minimum = basis_mat(a)[1, 1]
    return deepcopy(a.minimum)
  end
end

################################################################################
#
#  Norm
#
################################################################################

doc"""
***
    norm(I::NfOrdIdl) -> fmpz

> Returns the norm of $I$, that is, the cardinality of $\mathcal O/I$, where
> $\mathcal O$ is the order of $I$.
"""
function norm(a::NfOrdIdl)
  if isdefined(a, :norm)
    return deepcopy(a.norm)
  else
    a.norm = det(basis_mat(a))
    return deepcopy(a.norm)
  end
end

################################################################################
#
#  Equality
#
################################################################################

doc"""
***
    ==(x::NfOrdIdl, y::NfOrdIdl)

> Returns whether $x$ and $y$ are equal.
"""
function ==(x::NfOrdIdl, y::NfOrdIdl)
  return basis_mat(x) == basis_mat(y)
end

################################################################################
#
#  Inclusion of order elements in ideals
#
################################################################################

doc"""
***
    in(x::NfOrdElem, y::NfOrdIdl)

> Returns whether $x$ is contained in $y$.
"""
function in(x::NfOrdElem, y::NfOrdIdl)
  parent(x) != order(y) && error("Order of element and ideal must be equal")
  v = transpose(MatrixSpace(FlintZZ, degree(parent(x)), 1)(elem_in_basis(x)))
  t = FakeFmpqMat(v, fmpz(1))*basis_mat_inv(y)
  return t.den == 1
end

################################################################################
#
#  Binary operations
#
################################################################################

doc"""
***
    +(x::NfOrdIdl, y::NfOrdIdl)

> Returns $x + y$.
"""
function +(x::NfOrdIdl, y::NfOrdIdl)
  d = degree(order(x))
  H = vcat(basis_mat(x), basis_mat(y))
  g = gcd(minimum(x), minimum(y))
  H = sub(_hnf_modular_eldiv(H, g, :lowerleft), (d + 1):2*d, 1:d)
  return ideal(order(x), H)
end

doc"""
    intersection(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

> Returns $x \cap y$.
"""
function intersection(x::NfOrdIdl, y::NfOrdIdl)
  d = degree(order(x))
  H = vcat(basis_mat(x), basis_mat(y))
  K = _kernel(H)
  g = lcm(minimum(x),minimum(y))
  return ideal(order(x), _hnf_modular_eldiv(sub(K, 1:d, 1:d)*basis_mat(x), g, :lowerleft))
end

doc"""
    intersection(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl
    lcm(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdIdl

> Returns $x \cap y$.
"""
lcm(x::NfOrdIdl, y::NfOrdIdl) = intersection(x, y)

doc"""
    *(x::NfOrdIdl, y::NfOrdIdl)

> Returns $x \cdot y$.
"""
function *(x::NfOrdIdl, y::NfOrdIdl)
  return _mul(x, y)
end

function _mul(x::NfOrdIdl, y::NfOrdIdl)
  order(x) != order(y) && error("Not compatible")
  O = order(x)
  d = degree(O)
  l = minimum(x)*minimum(y)
  z = MatrixSpace(FlintZZ, degree(O)*degree(O), degree(O))()
  X = basis(x)
  Y = basis(y)
  for i in 1:d
    for j in 1:d
      t = elem_in_basis(X[i]*Y[j])
      for k in 1:d
        z[(i - 1)*d + j, k] = t[k]
      end
    end
  end
  # This is a d^2 x d matrix
  return ideal(O, sub(_hnf_modular_eldiv(z, l, :lowerleft),
                      (d*(d - 1) + 1):d^2, 1:d))
end

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function *(x::NfOrdElem, y::NfOrdIdl)
  parent(x) != order(y) && error("Orders of element and ideal must be equal")
  return ideal(parent(x), x) * y
end

*(x::NfOrdIdl, y::NfOrdElem) = y * x

function *(x::NfOrdIdl, y::fmpz)
  z = ideal(order(x), basis_mat(x)*y)
  if isdefined(x, :princ_gen)
    z.princ_gen = x.princ_gen * y
  end
  if isdefined(x, :princ_gen_special)
    z.princ_gen_special = (2, 0, x.princ_gen_special[x.princ_gen_special[1] + 1] * y)
  end
  return z
end

*(x::fmpz, y::NfOrdIdl) = y * x

*(x::NfOrdIdl, y::Integer) = x * fmpz(y)

*(x::Integer, y::NfOrdIdl) = y * x

################################################################################
#
#  Idempotents
#
################################################################################

doc"""
    idempotents(x::NfOrdIdl, y::NfOrdIdl) -> NfOrdElem, NfOrdElem

> Returns a tuple `(e, f)` consisting of elements `e in x`, `f in y` such that
> `1 = e + f`.
>
> If the ideals are not coprime, an error is raised.
"""
function idempotents(x::NfOrdIdl, y::NfOrdIdl)
  O = order(x)
  d = degree(O)

  # form the matrix
  #
  # ( 1 |  1  | 0 )
  # ( 0 | M_x | I )
  # ( 0 | M_y | 0 )

  V = MatrixSpace(FlintZZ, 1 + 2*d, 1 + 2*d)()

  u = elem_in_basis(one(O))

  V[1, 1] = 1

  for i in 1:d
    V[1, i + 1] = u[i]
  end

  Hecke._copy_matrix_into_matrix(V, 2, 2, basis_mat(x))
  Hecke._copy_matrix_into_matrix(V, 2 + d, 2, basis_mat(y))

  for i in 1:d
    V[1 + i, d + 1 + i] = 1
  end

  H = hnf(V) # upper right

  for i in 2:(1 + d)
    if H[1, i] != 0
      error("Ideals are not coprime")
    end
  end

  z = basis(x)[1]*H[1, d + 2]

  for i in 2:d
    z = z + basis(x)[i]*H[1, d + 1 + i]
  end

  @hassert :NfOrd 2 -z in x
  @hassert :NfOrd 2 1 + z in y

  return -z, 1 + z
end

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

doc"""
***
    mod(x::NfOrdElem, I::NfOrdIdl)

> Returns the unique element $y$ of the ambient order of $x$ with
> $x \equiv y \bmod I$ and the following property: If
> $a_1,\dotsc,a_d \in \Z_{\geq 1}$ are the diagonal entries of the unique HNF
> basis matrix of $I$ and $(b_1,\dotsc,b_d)$ is the coefficient vector of $y$,
> then $0 \leq b_i < a_i$ for $1 \leq i \leq d$.
"""
function mod(x::NfOrdElem, y::NfOrdIdl)
  parent(x) != order(y) && error("Orders of element and ideal must be equal")
  # this function assumes that HNF is lower left
  # !!! This must be changed as soon as HNF has a different shape

  O = order(y)
  a = elem_in_basis(x)
  #a = deepcopy(b)

  if isdefined(y, :princ_gen_special) && y.princ_gen_special[1] != 0
    for i in 1:length(a)
      a[i] = mod(a[i], y.princ_gen_special[1 + y.princ_gen_special[1]])
    end
    return O(a)
  end

  c = basis_mat(y)
  t = fmpz(0)
  for i in degree(O):-1:1
    t = fdiv(a[i], c[i,i])
    for j in 1:i
      a[j] = a[j] - t*c[i,j]
    end
  end
  z = O(a)
  return z
end

function mod(x::NfOrdElem, y::NfOrdIdl, preinv::Array{fmpz_preinvn_struct, 1})
  parent(x) != order(y) && error("Orders of element and ideal must be equal")
  # this function assumes that HNF is lower left
  # !!! This must be changed as soon as HNF has a different shape

  O = order(y)
  a = elem_in_basis(x) # this is already a copy

  if isdefined(y, :princ_gen_special) && y.princ_gen_special[1] != 0
    for i in 1:length(a)
      a[i] = mod(a[i], y.princ_gen_special[1 + y.princ_gen_special[1]])
    end
    return O(a)
  end

  c = basis_mat(y)
  q = fmpz()
  r = fmpz()
  for i in degree(O):-1:1
    fdiv_qr_with_preinvn!(q, r, a[i], c[i, i], preinv[i])
    for j in 1:i
      submul!(a[j], q, c[i, j])
    end
  end

  z = typeof(x)(O, a)
  return z
end

################################################################################
#
#  p-radical
#
################################################################################

doc"""
***
    pradical(O::NfOrd, p::fmpz) -> NfOrdIdl

> Given a prime number $p$, this function returns the $p$-radical
> $\sqrt{p\mathcal O}$ of $\mathcal O$, which is
> just $\{ x \in \mathcal O \mid \exists k \in \mathbf Z_{\geq 0} \colon x^k
> \in p\mathcal O \}$. It is not checked that $p$ is prime.
"""
function pradical(O::NfOrd, p::fmpz)
  j = clog(fmpz(degree(O)),p)

  @assert p^(j-1) < degree(O)
  @assert degree(O) <= p^j

  R = ResidueRing(ZZ,p)
  A = MatrixSpace(R, degree(O), degree(O))()
  for i in 1:degree(O)
    t = powermod(basis(O)[i], p^j, p)
    ar = elem_in_basis(t)
    for k in 1:degree(O)
      A[i,k] = ar[k]
    end
  end
  X = kernel(A)
  Mat = MatrixSpace(ZZ, 1, degree(O))
  MMat = MatrixSpace(R, 1, degree(O))
  if length(X) != 0
    m = lift(MMat(X[1]))
    for x in 2:length(X)
      m = vcat(m,lift(MMat(X[x])))
    end
    m = vcat(m,MatrixSpace(ZZ, degree(O), degree(O))(p))
  else
    m = MatrixSpace(ZZ, degree(O), degree(O))(p)
  end
  r = sub(_hnf(m, :lowerleft), rows(m) - degree(O) + 1:rows(m), 1:degree(O))
  return ideal(O, r)
end

doc"""
***
    pradical(O::NfOrd, p::Integer) -> NfOrdIdl

> Given a prime number $p$, this function returns the $p$-radical
> $\sqrt{p\mathcal O}$ of $\mathcal O$, which is
> just $\{ x \in \mathcal O \mid \exists k \in \mathbf Z_{\geq 0} \colon x^k
> \in p\mathcal O \}$. It is not checked that $p$ is prime.
"""
function pradical(O::NfOrd, p::Integer)
  return pradical(O, fmpz(p))
end

################################################################################
#
#  Ring of multipliers
#
################################################################################

doc"""
***
    ring_of_multipliers(I::NfOrdIdl) -> NfOrdGen

> Computes the order $(I : I)$, which is the set of all $x \in K$
> with $xI \subseteq I$.
"""
function ring_of_multipliers(a::NfOrdIdl)
  B = basis(a)
  O = order(a)
  bmatinv = basis_mat_inv(a)
  #print("First basis element is $(B[1]) \n with representation mat \n")
  #@vprint :NfOrd 1 "$(representation_mat(B[1]))\n"
  #@vprint :NfOrd 1 FakeFmpqMat(representation_mat(B[1]),ZZ(1))*bmatinv
  m = to_fmpz_mat(FakeFmpqMat(representation_mat(B[1]),ZZ(1))*bmatinv)
  for i in 2:degree(O)
    m = hcat(to_fmpz_mat(FakeFmpqMat(representation_mat(B[i]),ZZ(1))*basis_mat_inv(a)),m)
  end
  n = hnf(transpose(m))
  # n is upper right HNF
  n = transpose(sub(n, 1:degree(O), 1:degree(O)))
  b, d = pseudo_inv(n)
  #z = frac_ideal(O, FakeFmpqMat(b, d))
  return Order(nf(O), FakeFmpqMat(b, d)*basis_mat(O))
end


text:
function solve(prob::NonlinearProblem,
                          alg::AbstractNonlinearSolveAlgorithm, args...;
                          kwargs...)
  solver = init(prob, alg, args...; kwargs...)
  sol = solve!(solver)
  return sol
end

function init(prob::NonlinearProblem{uType, iip}, alg::AbstractBracketingAlgorithm, args...;
    alias_u0 = false,
    maxiters = 1000,
    kwargs...
  ) where {uType, iip}

  if !(prob.u0 isa Tuple)
    error("You need to pass a tuple of u0 in bracketing algorithms.")
  end

  if eltype(prob.u0) isa AbstractArray
    error("Bracketing Algorithms work for scalar arguments only")
  end

  if alias_u0
    left, right = prob.u0
  else
    left, right = deepcopy(prob.u0)
  end
  f = prob.f
  p = prob.p
  fl = f(left, p)
  fr = f(right, p)
  cache = alg_cache(alg, left, right,p, Val(iip))
  return BracketingImmutableSolver(1, f, alg, left, right, fl, fr, p, false, maxiters, DEFAULT, cache, iip)
end

function init(prob::NonlinearProblem{uType, iip}, alg::AbstractNewtonAlgorithm, args...;
    alias_u0 = false,
    maxiters = 1000,
    tol = 1e-6,
    internalnorm = DEFAULT_NORM,
    kwargs...
  ) where {uType, iip}

  if alias_u0
    u = prob.u0
  else
    u = deepcopy(prob.u0)
  end
  f = prob.f
  p = prob.p
  if iip
    fu = zero(u)
    f(fu, u, p)
  else
    fu = f(u, p)
  end
  cache = alg_cache(alg, f, u, p, Val(iip))
  return NewtonImmutableSolver(1, f, alg, u, fu, p, false, maxiters, internalnorm, DEFAULT, tol, cache, iip)
end

function solve!(solver::AbstractImmutableNonlinearSolver)
  solver = mic_check(solver)
  while !solver.force_stop && solver.iter < solver.maxiters
    solver = perform_step(solver, solver.alg, Val(solver.iip))
    @set! solver.iter += 1
  end
  if solver.iter == solver.maxiters
    @set! solver.retcode = MAXITERS_EXCEED
  end
  sol = get_solution(solver)
  return sol
end

"""
  mic_check(solver::AbstractImmutableNonlinearSolver)
  mic_check!(solver::AbstractNonlinearSolver)

Checks before running main solving iterations.
"""
function mic_check(solver::BracketingImmutableSolver)
  @unpack f, fl, fr = solver
  flr = fl * fr
  fzero = zero(flr)
  (flr > fzero) && error("Non bracketing interval passed in bracketing method.")
  if fl == fzero
    @set! solver.force_stop = true
    @set! solver.retcode = EXACT_SOLUTION_LEFT
  elseif fr == fzero
    @set! solver.force_stop = true
    @set! solver.retcode = EXACT_SOLUTION_RIGHT
  end
  solver
end

function mic_check(solver::NewtonImmutableSolver)
  solver
end

"""
  get_solution(solver::Union{BracketingImmutableSolver, BracketingSolver})
  get_solution(solver::Union{NewtonImmutableSolver, NewtonSolver})

Form solution object from solver types
"""
function get_solution(solver::BracketingImmutableSolver)
  return BracketingSolution(solver.left, solver.right, solver.retcode)
end

function get_solution(solver::NewtonImmutableSolver)
  return NewtonSolution(solver.u, solver.retcode)
end

"""
  reinit!(solver, prob)

Reinitialize solver to the original starting conditions
"""
function reinit!(solver::NewtonImmutableSolver, prob::NonlinearProblem{uType, true}) where {uType}
  @. solver.u = prob.u0
  @set! solver.iter = 1
  @set! solver.force_stop = false
  return solver
end

function reinit!(solver::NewtonImmutableSolver, prob::NonlinearProblem{uType, false}) where {uType}
  @set! solver.u = prob.u0
  @set! solver.iter = 1
  @set! solver.force_stop = false
  return solver
end

ext: jl


text:
function test_mc(mc, x)
    # Check if loaded/replayed mc matches original
    for f in fieldnames(typeof(mc.p))
        @test getfield(mc.p, f) == getfield(x.p, f)
    end
    @test mc.conf == x.conf
    @test mc.model.L == x.model.L
    @test mc.model.dims == x.model.dims
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    # @test mc.model.neighs == x.model.neighs
    @test mc.model.energy[] == x.model.energy[]
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.thermalization_measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            r = getfield(v, f) == getfield(x.measurements[k], f)
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

@testset "MC" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=0.66, thermalization=33, sweeps=123, recorder=ConfigRecorder)
    run!(mc, verbose=false)
    save("testfile.jld2", mc)
    x = load("testfile.jld2")
    rm("testfile.jld2")
    test_mc(mc, x)

    x.measurements = MonteCarlo.default_measurements(mc, model)
    x.last_sweep = 0
    replay!(x, verbose=false)
    test_mc(mc, x)

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(model, beta=1.0, sweeps=10_000_000, measure_rate=10_000, recorder=ConfigRecorder)
    state = run!(
        mc, verbose = false,
        safe_before = now() + Second(1),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    mc, state = resume!(
        "resumable_testfile.jld",
        verbose = false,
        safe_before = now() + Second(20),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = IsingModel(dims=2, L=10)
    mc = MC(model, beta=1.0, sweeps=10_000length(cs), measure_rate=10_000, recorder=ConfigRecorder)
    state = run!(mc, verbose = false)
    @test mc.configs.configs == cs.configs
    @test mc.configs.rate == cs.rate
    rm("resumable_testfile.jld")
end

function test_dqmc(mc, x)
    for f in fieldnames(typeof(mc.p))
        @test getfield(mc.p, f) == getfield(x.p, f)
    end
    # @test mc.conf == x.conf
    @test mc.model.mu == x.model.mu
    @test mc.model.t == x.model.t
    @test mc.model.U == x.model.U
    for f in fieldnames(typeof(mc.model.l))
        @test getfield(mc.model.l, f) == getfield(x.model.l, f)
    end
    @test mc.model.flv == x.model.flv
    for (k, v) in mc.thermalization_measurements
        for f in fieldnames(typeof(v))
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                getfield(v, f).B == getfield(x.measurements[k], f).B
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    for (k, v) in mc.measurements
        for f in fieldnames(typeof(v))
            v isa MonteCarlo.DQMCMeasurement && f == :output && continue
            r = if getfield(v, f) isa LightObservable
                # TODO
                # implement == for LightObservable in MonteCarloObservable
                # TODO: implement ≈ for LightObservable, LogBinner, etc
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(getfield(v, f).B.compressors)
                    r = r && (a.B.compressors[i].value ≈ b.B.compressors[i].value)
                    r = r && (a.B.compressors[i].switch ≈ b.B.compressors[i].switch)
                end
                r = r && (a.B.x_sum ≈ b.B.x_sum)
                r = r && (a.B.x2_sum ≈ b.B.x2_sum)
                r = r && (a.B.count ≈ b.B.count)
            elseif getfield(v, f) isa LogBinner
                r = true
                a = getfield(v, f)
                b = getfield(x.measurements[k], f)
                for i in eachindex(a.compressors)
                    r = r && (a.compressors[i].value ≈ b.compressors[i].value)
                    r = r && (a.compressors[i].switch ≈ b.compressors[i].switch)
                end
                r = r && (a.x_sum ≈ b.x_sum)
                r = r && (a.x2_sum ≈ b.x2_sum)
                r = r && (a.count ≈ b.count)
            else
                getfield(v, f) == getfield(x.measurements[k], f)
            end
            r != true && @info "Check failed for $k -> $f"
            @test r
        end
    end
    nothing
end

@testset "DQMC" begin
    model = HubbardModelAttractive(4, 2, t = 1.7, U = 5.5)
    mc = DQMC(model, beta=1.0, thermalization=21, sweeps=117, measure_rate = 1)
    mc[:CDC] = charge_density_correlation(mc, model)
    t = time()
    run!(mc, verbose=false)
    t = time() - t

    save("testfile.jld", mc)
    x = load("testfile.jld")
    rm("testfile.jld")
    @test mc.conf == x.conf

    # Repeat these tests once with x being replayed rather than loaded
    test_dqmc(mc, x)

    # Check everything again with x being a replayed simulation
    x[:CDC] = charge_density_correlation(x, model)
    x.last_sweep = 0
    replay!(x, verbose=false)
    test_dqmc(mc, x)

    # Test resume

    # Run for 1s with known RNG
    Random.seed!(123)
    model = HubbardModelAttractive(2, 2, t = 1.7, U = 5.5)
    mc = DQMC(model, beta=1.0, sweeps=10_000_000, measure_rate=100)
    mc[:CDC] = charge_density_correlation(mc, model)

    state = run!(
        mc, verbose = false,
        safe_before = now() + Second(1),
        grace_period = Millisecond(0),
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) > 1 "No measurements have been taken. Test with more time!"
    L = length(cs)

    # Test whether safe file gets overwritten correctly
    mc, state = resume!(
        "resumable_testfile.jld2",
        verbose = false,
        safe_before = now() + Second(10),
        grace_period = Millisecond(0),
        overwrite = true,
        resumable_filename = "resumable_testfile.jld2"
    )

    @test state == false
    cs = deepcopy(mc.configs)
    @assert length(cs) - L > 1 "No new measurements have been taken. Test with more time!"
    @test isfile("resumable_testfile.jld2")

    # Test whether data from resumed simulation is correct
    Random.seed!(123)
    model = HubbardModelAttractive(2, 2, t = 1.7, U = 5.5)
    mc = DQMC(model, beta=1.0, sweeps=100length(cs), measure_rate=100)
    state = run!(mc, verbose = false)
    @test mc.configs.configs == cs.configs
    @test mc.configs.rate == cs.rate
    rm("resumable_testfile.jld2")
end


text:
function initialize!(solver::iLQRSolver2)
    set_verbosity!(solver.opts)
    clear_cache!(solver.opts)

    solver.ρ[1] = solver.opts.bp_reg_initial
    solver.dρ[1] = 0.0

    # Initial rollout
    rollout!(solver)
    cost!(solver.obj, solver.Z)
end

# Generic solve methods
"iLQR solve method (non-allocating)"
function solve!(solver::iLQRSolver{T}) where T<:AbstractFloat
	initialize!(solver)

    solver.stats.iterations = 0
    Z = solver.Z; Z̄ = solver.Z̄;

    n,m,N = size(solver)
    J = Inf
    _J = get_J(solver.obj)
    J_prev = sum(_J)

    for i = 1:solver.opts.iterations
        J = step!(solver, J_prev)

        # check for cost blow up
        if J > solver.opts.max_cost_value
            # @warn "Cost exceeded maximum cost"
            return solver
        end

        copy_trajectories!(solver)

        dJ = abs(J - J_prev)
        J_prev = copy(J)
        gradient_todorov!(solver)

        record_iteration!(solver, J, dJ)
        evaluate_convergence(solver) ? break : nothing
    end
    return solver
end

function step!(solver::iLQRSolver2, J)
    Z = solver.Z
    state_diff_jacobian!(solver.G, solver.model, Z)
    # discrete_jacobian!(solver.∇F, solver.model, Z)
	dynamics_expansion!(solver.D, solver.model, solver.Z)
    cost_expansion!(solver.Q, solver.obj, solver.Z)
	error_expansion!(solver.D, solver.model, solver.G)
	error_expansion!(solver.Q, solver.model, Z, solver.G)
	if solver.opts.static_bp
    	ΔV = static_backwardpass!(solver)
	else
		ΔV = backwardpass!(solver)
	end
    forwardpass!(solver, ΔV, J)
end

"""
$(SIGNATURES)
Simulate the system forward using the optimal feedback gains from the backward pass,
projecting the system on the dynamically feasible subspace. Performs a line search to ensure
adequate progress on the nonlinear problem.
"""
function forwardpass!(solver::iLQRSolver, ΔV, J_prev)
    Z = solver.Z; Z̄ = solver.Z̄
    obj = solver.obj

    _J = get_J(obj)
    J::Float64 = Inf
    α = 1.0
    iter = 0
    z = -1.0
    expected = 0.0
    flag = true

    while (z ≤ solver.opts.line_search_lower_bound || z > solver.opts.line_search_upper_bound) && J >= J_prev

        # Check that maximum number of line search decrements has not occured
        if iter > solver.opts.iterations_linesearch
            for k in eachindex(Z)
                Z̄[k].z = Z[k].z
            end
            cost!(obj, Z̄)
            J = sum(_J)

            z = 0
            α = 0.0
            expected = 0.0

            regularization_update!(solver, :increase)
            solver.ρ[1] += solver.opts.bp_reg_fp
            break
        end

        # Otherwise, rollout a new trajectory for current alpha
        flag = rollout!(solver, α)

        # Check if rollout completed
        if ~flag
            # Reduce step size if rollout returns non-finite values (NaN or Inf)
            # @logmsg InnerIters "Non-finite values in rollout"
            iter += 1
            α /= 2.0
            continue
        end

        # Calcuate cost
        cost!(obj, Z̄)
        J = sum(_J)

        expected::Float64 = -α*(ΔV[1] + α*ΔV[2])
        if expected > 0.0
            z::Float64  = (J_prev - J)/expected
        else
            z = -1.0
        end

        iter += 1
        α /= 2.0
    end

    if J > J_prev
        error("Error: Cost increased during Forward Pass")
    end

    @logmsg InnerLoop :expected value=expected
    @logmsg InnerLoop :z value=z
    @logmsg InnerLoop :α value=2*α
    @logmsg InnerLoop :ρ value=solver.ρ[1]

    return J

end

function copy_trajectories!(solver::iLQRSolver)
    for k = 1:solver.N
        solver.Z[k].z = solver.Z̄[k].z
    end
end

"""
Stash iteration statistics
"""
function record_iteration!(solver::iLQRSolver, J, dJ)
    solver.stats.iterations += 1
    i = solver.stats.iterations::Int
    solver.stats.cost[i] = J
    solver.stats.dJ[i] = dJ
    solver.stats.gradient[i] = mean(solver.grad)

    @logmsg InnerLoop :iter value=i
    @logmsg InnerLoop :cost value=J
    @logmsg InnerLoop :dJ   value=dJ
    @logmsg InnerLoop :grad value=solver.stats.gradient[i]
    # @logmsg InnerLoop :zero_count value=solver.stats[:dJ_zero_counter][end]
    if solver.opts.verbose
        print_level(InnerLoop)
    end
    return nothing
end

"""
$(SIGNATURES)
    Calculate the problem gradient using heuristic from iLQG (Todorov) solver
"""
function gradient_todorov!(solver::iLQRSolver2)
	tmp = solver.S[end].u
    for k in eachindex(solver.d)
		tmp .= abs.(solver.d[k])
		u = abs.(control(solver.Z[k])) .+ 1
		tmp ./= u
		solver.grad[k] = maximum(tmp)
    end
end

"""
$(SIGNATURES)
Check convergence conditions for iLQR
"""
function evaluate_convergence(solver::iLQRSolver)
    # Get current iterations
    i = solver.stats.iterations

    # Check for cost convergence
    # note the  dJ > 0 criteria exists to prevent loop exit when forward pass makes no improvement
    if 0.0 < solver.stats.dJ[i] < solver.opts.cost_tolerance
        return true
    end

    # Check for gradient convergence
    if solver.stats.gradient[i] < solver.opts.gradient_norm_tolerance
        return true
    end

    # Check total iterations
    if i >= solver.opts.iterations
        return true
    end

    # Outer loop update if forward pass is repeatedly unsuccessful
    if solver.stats.dJ_zero_counter > solver.opts.dJ_counter_limit
        return true
    end

    return false
end

"""
$(SIGNATURES)
Update the regularzation for the iLQR backward pass
"""
function regularization_update!(solver::iLQRSolver,status::Symbol=:increase)
    # println("reg $(status)")
    if status == :increase # increase regularization
        # @logmsg InnerLoop "Regularization Increased"
        solver.dρ[1] = max(solver.dρ[1]*solver.opts.bp_reg_increase_factor, solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = max(solver.ρ[1]*solver.dρ[1], solver.opts.bp_reg_min)
        # if solver.ρ[1] > solver.opts.bp_reg_max
        #     @warn "Max regularization exceeded"
        # end
    elseif status == :decrease # decrease regularization
        # TODO: Avoid divides by storing the decrease factor (divides are 10x slower)
        solver.dρ[1] = min(solver.dρ[1]/solver.opts.bp_reg_increase_factor, 1.0/solver.opts.bp_reg_increase_factor)
        solver.ρ[1] = solver.ρ[1]*solver.dρ[1]*(solver.ρ[1]*solver.dρ[1]>solver.opts.bp_reg_min)
    end
end


text:
"""
simulation & machine learning of diffusion advection PDE
"""

using IterTools
using Plots
using Random
using Flux
using LinearAlgebra
Random.seed!(1)
include("../src/operators.jl")
# using EquivariantOperators

# make grid
n = 2
dx = 0.02
cell = dx * Matrix(I, n, n)
rmax = 1.0
grid = Grid(cell, rmax)
sz = size(grid)
▽ = Op(:▽, cell)
blur=Op(:Gaussian,cell;σ=.1)

# diffusion advection with point emitter
p = [0.2,0.5, 0.5,.8]
D,vx,vy,k=p
v = fill([vx,vy],sz)
u0 = zeros(sz)
s = zeros(sz)
put!(s, grid, [0.0, 0.0], 1.)
s=blur(s)
f(u, p, t) = D * (▽ ⋅ ▽(u)) - v .⋅ ▽(u).+s*k

# simulate PDE
using DifferentialEquations
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-3, abstol = 1e-3)

# plot
using Plots
gr()
anim = Animation()

t = 0:0.02:1
for t in t
heatmap(sol(t)[:, :, 1], clim=(0,10))
frame(anim)
end
gif(anim, "f.gif", fps = 10)

##
data = [(sol(t), f(sol(t), 0, 0)) for t in t]
# op = Op(Radfunc(),-1e-6, 2dx, cell)
# ps=Flux.params(op)
p_ = ones(length(p))
ps = Flux.params(p_)

function loss(u, du)
    D,vx,vy,k=p_
    v = fill([vx,vy],sz)
    duhat= D * (▽ ⋅ ▽(u)) - v .⋅ ▽(u).+s*k
    @show l = nae(duhat, du)
end

loss(data[1]...)
opt = ADAM(0.1)

Flux.@epochs 2 Flux.train!(loss, ps, data, opt)

@show p,p_
# (p, p_) = ([0.2, 0.5, 0.5, 0.8], [0.2113613015968995, 0.4977070976688276, 0.4977070976688435, 0.791498663145966])


text:
using LibGit2

# Set to "true" to save the source download. Otherwise, remove the source and related
# object files after a successful build.
keep_source = false

# Download the source code
url = "https://github.com/pmem/pmdk/"
branch = "stable-1.5"
localdir = joinpath(@__DIR__, "pmdk")

# Cleanup leftovers
ispath(localdir) && rm(localdir; force = true, recursive = true)

LibGit2.clone(url, localdir; branch = branch)

# Navigate and build
olddir = pwd()
try
    cd(localdir)

    nprocs = parse(Int, read(`nproc`, String))
    run(`make -j$nprocs CC=clang CXX=clang++`)
    installdir = joinpath(@__DIR__, "usr")
    run(`make install prefix=$installdir`)

    # Delete cloned directory to cut down on accumulated size.
    keep_source || rm(localdir; force = true, recursive = true)
finally
    cd(olddir)
end


text:
########################################################################
# AutoGrad functions extended for CuArray (copied from karray.jl)

# using CUDA: CuArray, CuPtr
# using Knet: c2i, @knet8
using Base.Broadcast: Broadcasted
using AutoGrad: Value, recording

import AutoGrad: Sparse, matches, ungetindex, addto!, addtoindex!, zeroslike
import Base: copyto!
import LinearAlgebra: axpy!

Sparse(a::CuArray{T,N},v,i) where {T,N} = Sparse{T,N}(a,v,i)

axpy!(a, x::Sparse, y::CuArray) = addto!(y, a*x)

function copyto!(a::CuArray, bc::Broadcasted{S,A,F,X}) where
    {S, A, F <: Union{typeof(+),typeof(-)}, X <: Tuple{Any,Sparse}}
    (b,c) = bc.args
    if !(size(a) == size(b) == size(c.container))
        a .= bc.f.(b, full(c))
        return a
    end
    a === b || copyto!(a, b)
    F <: typeof(-) && (c = -c)
    addto!(a, c)
    return a
end

matches(a::CuArray,b::CuArray)=(size(a)==size(b))

function addto!(a::CuArray{T},b::CuArray{T}) where {T}
    if AutoGrad.recording(); a = copy(a); end  # support highorder gradients
    axpy!(1,b,a) # (a+b)
end

function addto!(a::CuArray,b::Sparse)
    @assert size(a) == size(b.container)
    if AutoGrad.recording(); a = copy(a); end  # support highorder gradients
    for (idx,val) in zip(b.indices, b.values)
        addtoindex!(a, val, idx...)
    end
    return a
end

addto!(a::Sparse, b::CuArray) = addto!(b, a)

import Base: +, -
+(a::CuArray, s::Sparse) = addto!(copy(a), s)
+(s::Sparse, a::CuArray) = addto!(copy(a), s)
-(a::CuArray, s::Sparse) = addto!(copy(a), -s)
-(s::Sparse, a::CuArray) = addto!(-a, s)

function ungetindex(x::CuArray{T},dxi,i) where T
    if isbitstype(T)
        if dxi isa Value
            forw(addto!, zeroslike(x), forw(ungetindex, x, dxi, i))
        elseif recording()
            addtoindex!(zero(x), dxi, i...)
        else
            Sparse(x,[dxi],[i])
        end
    else
        # Using addtoindex! instead of setindex! to handle repeated indices
        addtoindex!(Array{Union{T,Nothing}}(nothing, size(x)), dxi, i...)
    end
end

zeroslike(a::CuArray)=zero(a) # Still need this because zero(::Array{!isbits}) is not defined

# This only works when there are no repeated indices. This is true for index types:
# Real, (Real...), CartesianIndex, Colon, AbstractArray{Bool}, Range, EmptyArray
# and pairs of Union{Real,AbstractUnitRange,Colon} and (Colon,Range)
addtoindex!(A::CuArray, X, I::Union{Real,Colon,AbstractRange,AbstractArray{Bool}}...) =
    setindex!(A, addto!(getindex(A,I...), X), I...)
addtoindex!(A::CuArray, X, I...) =
    setindex!(A, getindex(A,I...) .+ X, I...)

# The following index types may have repeated indices:
# AbstractArray{Real}, AbstractArray{CartesianIndex}, (Colon,AbstractVector{Real}), (AbstractVector{Real},Colon)

addtoindex!(A::CuArray, X, I::AbstractArray{<:CartesianIndex})=addtoindex!(A,X,c2i(size(A),I))

for F in (32,64); T=Symbol("Float$F"); @eval begin

    function addtoindex!(A::CuArray{$T}, X, I::AbstractArray{<:Real})
        I = CuArray{Int32}(I)
        X = CuArray{$T}(X)
        @knet8($("addents_$F"),(Cint,CuPtr{Int},CuPtr{$T},CuPtr{$T}),
               length(I), I, A, X)
        return A
    end

    function addtoindex!(A::CuArray{$T}, X, ::Colon, I::AbstractArray{<:Real})
        I = CuArray{Int32}(I)
        X = CuArray{$T}(X)
        @knet8($("addcols_$F"),(Cint,Cint,Cint,CuPtr{Int},CuPtr{$T},CuPtr{$T}),
               size(A,1), size(A,2), length(I), I, A, X)
        return A
    end

    function addtoindex!(A::CuArray{$T}, X, I::AbstractArray{<:Real}, ::Colon)
        I = CuArray{Int32}(I)
        X = CuArray{$T}(X)
        @knet8($("addrows_$F"),(Cint,Cint,Cint,CuPtr{Int},CuPtr{$T},CuPtr{$T}),
               size(A,1), size(A,2), length(I), I, A, X)
        return A
    end

    addtoindex!(A::CuArray{$T}, X, I::AbstractArray{Bool})=addtoindex!(A,X,findall(vec(I)))
    addtoindex!(A::CuArray{$T}, X, c::Colon, I::AbstractArray{Bool})=addtoindex!(A,X,c,findall(vec(I)))
    addtoindex!(A::CuArray{$T}, X, I::AbstractArray{Bool}, c::Colon)=addtoindex!(A,X,findall(vec(I)),c)

    # Fix #578: addtoindex! ambiguity:
    addtoindex!(A::CuArray{$T}, X, I::UnitRange{<:Real})=addtoindex!(A,X,CuArray{Int32}(I))
    addtoindex!(A::CuArray{$T}, X, c::Colon, I::UnitRange{<:Real})=addtoindex!(A,X,c,CuArray{Int32}(I))
    addtoindex!(A::CuArray{$T}, X, I::UnitRange{<:Real}, c::Colon)=addtoindex!(A,X,CuArray{Int32}(I),c)

end; end


text:
module Gameface

using Niancat.Users
using Niancat.Http

"""
Game represents a word game of some kind.

One can send commands to a `Game`, and get `Response`s.
"""
abstract type Game end

"""
Response represents a human-readable response to a player.
"""
abstract type Response end

"""
GameCommand represents any command sent to a game.
"""
abstract type GameCommand end

"""
GameNotification is a message not sent as a response to a game command.

# Example: A new puzzle
A new puzzle can be set directly in the service, by a system user, and not via a chat
channel like Slack or Discord. In this case, the game must notify the main channels for all teams
in that instance, that a new puzzle is set, along with the scoreboard for the previous round.
This is a notification, as it doesn't fit into the command/response flow.
"""
abstract type GameNotification end

struct GameUserEvent
    eventtype::Int
    user::User
    round::String
    data::String
end

const UserEvent_Solution = 1
const GameEvent_Solution = 1

struct GameEvent
    eventtype::Int
    round::String
    data::String
end

"""
    gamecommand(game::Game, user::User, command::GameCommand) :: Response

`gamecommand` sends a command to a given game.
"""
gamecommand(::Game, ::User, ::GameCommand) :: Response = error("Implement this in Game subtypes")

"""
NoResponse is an empty response object, for when a command should not respond
with any message.
"""
struct NoResponse <: Response end

"""
    gameround(game::Game) :: String

The round identifier for the current round of the game instance `game`.
"""
gameround(::Game) :: String = error("Implement this in Game subtypes")

"""
    gameinstanceid(game::Game) :: Int

The database id for this instance of the game.

The default implementation assumes that the game has a field `gameservice`
of type `GameService`, from which it can find the game instance id.
"""
gameinstanceid(g::Game) :: Int = gameinstanceid(g.gameservice)

"""
Score represents a score for a given round of a game.
"""
struct Score
    gameinstanceid::Int
    round::String
    key::String
    value::Float32
end

"""
GameEventPersistence records game events, like scores.
"""
abstract type GameEventPersistence end

score!(::GameEventPersistence, ::User, ::Score) = error("Implement this in GameEventPersistence subtypes")

"""
InstanceInfo encapsulates static, but not dynamic, information on an instance.

# Example
Instance name, and instance database id.

Does not contain dynamic information, such as which teams belong to an instance.

# Limitations
Note that this implementation implies that instances cannot easily be renamed, as you would have
to update the `instanceinfo` member on the `GameService` struct, which is immutable.
The rename would only take full effect after game instances were reloaded.
"""
struct InstanceInfo
    databaseid::Int
    name::String
end

abstract type LanguageDictionary end

"""
GameService is an interface between the games and the service.
Each game instance has its own GameService object.

For example, sending notifications to users from a game is done via the GameService.
Recording user scores is done via the GameService.
"""
abstract type GameService end

score!(::GameService, ::User, ::Score) = error("Implement score! in GameService subtypes")
notify!(::GameService, ::GameNotification) = error("Implement notify! in GameService subtypes")
event!(::GameService, ::GameUserEvent) = error("Implement event! for GameUserEvent in GameService subtypes")
event!(::GameService, ::GameEvent) = error("Implement event! for GameEvent in GameService subtypes")
finddictionary(::GameService, ::AbstractString) :: LanguageDictionary = error("Implement finddictionary in GameService subtypes")

##
## Exports
##

export Game, Response, NoResponse, gamecommand, gameinstanceid, gameround, Score, GameEventPersistence
export GameCommand, GameNotification
export InstanceInfo

# GameService export
export GameService, score!, notify!, finddictionary

# Events
export event!, GameUserEvent, GameEvent

# LanguageDictionary
export LanguageDictionary

end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    basisfun(i::I,u::F,p::I,U::Vector{F}
            )::Vector{F} where {I<:Integer,F<:AbstractFloat}

Calculates the B-Spline basis function at `u`.

# Arguments:
- `i`: knot span from [`findspan`](@ref)
- `u`: parametric point
- `p`: spline degree
- `U`: knot sequence

# Output:
- `N`: basis function vector of lenght `p+1`

# Examples:
```julia
julia> N = basisfun(i,u,p,U)
```

# Reference:
- Algorithm A2.2 from 'The NURBS BOOK' pg70 ([1997Piegl](@cite))
"""
function basisfun(i::I,u::F,p::I,U::Vector{F}
                 )::Vector{F} where {I<:Integer,F<:AbstractFloat}
# TODO: the original C function got U and N passed as reference
# TODO: this function has optimization potential

i = i + 1;

# allocate result and work space
N     = zeros(F,p+1);
left  = zeros(F,p+1);
right = zeros(F,p+1);

N[1] = 1;
for j in 1:p
  left[j+1]  = u - U[i+1-j];
  right[j+1] = U[i+j] - u;
  saved = 0;

  for r in 0:j-1
    temp = N[r+1] / (right[r+2] + left[j-r+1]);
    N[r+1] = saved + right[r+2] * temp;
    saved = left[j-r+1] * temp;
  end # for r

  N[j+1] = saved;
end # for j
return N;
end # basisfun


text:
using ProfileFile
using Base.Test

include("test.jl")


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Regularized decomposition
# ------------------------------------------------------------
@with_kw mutable struct RDData{T <: AbstractFloat}
    Q̃::T = 1e10
    σ::T = 1.0
    incumbent::Int = 1
    major_iterations::Int = 0
    minor_iterations::Int = 0
end

@with_kw mutable struct RDParameters{T <: AbstractFloat}
    τ::T = 1e-6
    γ::T = 0.9
    σ::T = 1.0
    σ̅::T = 4.0
    σ̲::T = 0.5
end

"""
    RegularizedDecomposition

Functor object for using regularized decomposition regularization in an L-shaped algorithm. Create by supplying an [`RD`](@ref) object through `regularize` in `LShaped.Optimizer` or by setting the [`Regularizer`](@ref) attribute.

...
# Parameters
- `σ::AbstractFloat = 1.0`: Initial value of regularization parameter. Controls the relative penalty of the deviation from the current major iterate.
- `σ̅::AbstractFloat = 4.0`: Maximum value of the regularization parameter.
- `σ̲::AbstractFloat = 0.5`: Minimum value of the regularization parameter.
- `log::Bool = true`: Specifices if L-shaped procedure should be logged on standard output or not.
- `penaltyterm::PenaltyTerm = Quadratic`: Specify penaltyterm variant ([`Quadratic`](@ref), [`InfNorm`](@ref), [`ManhattanNorm`][@ref])
...
"""
struct RegularizedDecomposition{T <: AbstractFloat, A <: AbstractVector, PT <: AbstractPenaltyTerm} <: AbstractRegularization
    data::RDData{T}
    parameters::RDParameters{T}

    decisions::DecisionMap
    projection_targets::Vector{MOI.VariableIndex}
    ξ::Vector{Decision{T}}

    Q̃_history::A
    σ_history::A
    incumbents::Vector{Int}

    penaltyterm::PT

    function RegularizedDecomposition(decisions::DecisionMap, ξ₀::AbstractVector, penaltyterm::AbstractPenaltyTerm; kw...)
        T = promote_type(eltype(ξ₀), Float32)
        A = Vector{T}
        ξ = map(ξ₀) do val
            KnownDecision(val, T)
        end
        PT = typeof(penaltyterm)
        return new{T, A, PT}(RDData{T}(),
                             RDParameters{T}(; kw...),
                             decisions,
                             Vector{MOI.VariableIndex}(undef, length(ξ₀)),
                             ξ,
                             A(),
                             A(),
                             Vector{Int}(),
                             penaltyterm)
    end
end

function initialize_regularization!(lshaped::AbstractLShaped, rd::RegularizedDecomposition{T}) where T <: AbstractFloat
    # Add projection targets
    add_projection_targets!(rd, lshaped.master)
    # Sense-correct initial objective
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    rd.data.Q̃ *= coeff
    # Prepare penalty constant
    rd.data.σ = rd.parameters.σ
    push!(rd.σ_history,rd.data.σ)
    @unpack σ = rd.data
    # Initialize penalty
    initialize_penaltyterm!(rd.penaltyterm,
                            lshaped.master,
                            1 / (2 * σ),
                            all_decisions(rd.decisions),
                            rd.projection_targets)
    return nothing
end

function restore_regularized_master!(lshaped::AbstractLShaped, rd::RegularizedDecomposition)
    # Delete penalty-term
    remove_penalty!(rd.penaltyterm, lshaped.master)
    # Delete projection targets
    for var in rd.projection_targets
        MOI.delete(lshaped.master, var)
        StochasticPrograms.remove_decision!(rd.decisions, var)
    end
    empty!(rd.projection_targets)
    return nothing
end

function filter_variables!(rd::RegularizedDecomposition, list::Vector{MOI.VariableIndex})
    # Filter projection targets
    filter!(vi -> !(vi in rd.projection_targets), list)
    # Filter any auxiliary penaltyterm variables
    remove_penalty_variables!(rd.penaltyterm, list)
    return nothing
end

function filter_constraints!(rd::RegularizedDecomposition, list::Vector{<:CI})
    # Filter any auxiliary penaltyterm constraints
    remove_penalty_constraints!(rd.penaltyterm, list)
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, rd::RegularizedDecomposition)
    @unpack Q̃,σ,incumbent = rd.data
    push!(rd.Q̃_history, Q̃)
    push!(rd.σ_history, σ)
    push!(rd.incumbents, incumbent)
    return nothing
end

function log_regularization!(lshaped::AbstractLShaped, t::Integer, rd::RegularizedDecomposition)
    @unpack Q̃,σ,incumbent = rd.data
    rd.Q̃_history[t] = Q̃
    rd.σ_history[t] = σ
    rd.incumbents[t] =  incumbent
    return nothing
end

function take_step!(lshaped::AbstractLShaped, rd::RegularizedDecomposition)
    @unpack Q,θ = lshaped.data
    @unpack τ = lshaped.parameters
    @unpack σ = rd.data
    @unpack γ,σ̅,σ̲ = rd.parameters
    t = timestamp(lshaped)
    σ̃ = incumbent_trustregion(lshaped, t, rd)
    Q̃ = incumbent_objective(lshaped, t, rd)
    sense = MOI.get(lshaped.master, MOI.ObjectiveSense())
    coeff = sense == MOI.MIN_SENSE ? 1.0 : -1.0
    need_update = false
    if abs(θ-Q) <= τ*(1+abs(θ)) || coeff*Q <= coeff*Q̃ + τ || rd.data.major_iterations == 0
        need_update = true
        x = current_decision(lshaped)
        for i in eachindex(rd.ξ)
            rd.ξ[i].value = x[i]
        end
        rd.data.Q̃ = copy(Q)
        rd.data.incumbent = t
        rd.data.major_iterations += 1
    else
        rd.data.minor_iterations += 1
    end
    rd.data.σ = if coeff*Q + τ <= coeff*(1-γ)*Q̃ + coeff*γ*θ
        max(σ, min(σ̅, 2*σ))
    elseif coeff*Q - τ >= coeff*γ*Q̃ + coeff*(1-γ)*θ
        min(σ, max(σ̲, 0.5*σ))
    else
        σ
    end
    need_update |= abs(rd.data.σ - σ) > τ
    if need_update
        @unpack σ = rd.data
        update_penaltyterm!(rd.penaltyterm,
                            lshaped.master,
                            1 / (2 * σ),
                            all_decisions(rd.decisions),
                            rd.projection_targets)
    end
    return nothing
end

# API
# ------------------------------------------------------------
"""
    RD

Factory object for [`RegularizedDecomposition`](@ref). Pass to `regularize` in `LShaped.Optimizer` or set the [`Regularizer`](@ref) attribute. Equivalent factory calls: `RD`, `WithRD`, `RegularizedDecomposition`, `WithRegularizedDecomposition`. See ?RegularizedDecomposition for parameter descriptions.

"""
mutable struct RD <: AbstractRegularizer
    penaltyterm::AbstractPenaltyTerm
    parameters::RDParameters{Float64}
end
RD(; penaltyterm = Quadratic(), kw...) = RD(penaltyterm, RDParameters(; kw...))
WithRD(; penaltyterm = Quadratic(), kw...) = RD(penaltyterm, RDParameters(; kw...))
RegularizedDecomposition(; penaltyterm = Quadratic(), kw...) = RD(penaltyterm, RDParameters(; kw...))
WithRegularizedDecomposition(; penaltyterm = Quadratic(), kw...) = RD(penaltyterm, RDParameters(; kw...))

function MOI.get(rd::RD, ::RegularizationPenaltyTerm)
    return rd.penaltyterm
end

function MOI.set(rd::RD, ::RegularizationPenaltyTerm, penaltyterm::AbstractPenaltyTerm)
    return rd.penaltyterm = penaltyterm
end

function (rd::RD)(decisions::DecisionMap, x::AbstractVector)
    return RegularizedDecomposition(decisions, x, rd.penaltyterm; type2dict(rd.parameters)...)
end

function str(::RD)
    return "L-shaped using regularized decomposition"
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 this file:
# - the `fmi3InstanceState`-enum which mirrors the internal FMU state (state-machine, not the system state)
# - the `FMU3Instance`-struct
# - the `FMU3`-struct
# - string/enum-converters for FMI-attribute-structs (e.g. `fmi3StatusToString`, ...)

# this is a custom type to catch the internal state of the instance

@enum fmi3InstanceState begin
    fmi3InstanceStateInstantiated       # after instantiation
    fmi3InstanceStateInitializationMode # after finishing initialization
    fmi3InstanceStateEventMode
    fmi3InstanceStateStepMode
    fmi3InstanceStateClockActivationMode
    fmi3InstanceStateContinuousTimeMode
    fmi3InstanceStateConfigurationMode
    fmi3InstanceStateReconfigurationMode
    fmi3InstanceStateTerminated
    fmi3InstanceStateError
    fmi3InstanceStateFatal
end

"""
ToDo.
"""
mutable struct FMU3InstanceEnvironment
    logStatusOK::Bool
    logStatusWarning::Bool
    logStatusDiscard::Bool
    logStatusError::Bool
    logStatusFatal::Bool
    logStatusPending::Bool

    function FMU3InstanceEnvironment()
        inst = new()
        inst.logStatusOK = true
        inst.logStatusWarning = true
        inst.logStatusDiscard = true
        inst.logStatusError = true
        inst.logStatusFatal = true
        inst.logStatusPending = true
        return inst
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.2.1. Header Files and Naming of Functions

The mutable struct represents a pointer to an FMU specific data structure that contains the information needed to process the model equations or to process the co-simulation of the model/subsystem represented by the FMU.
"""
mutable struct FMU3Instance
    compAddr::Ptr{Nothing}
    fmu # ::FMU3
    state::fmi3InstanceState
    instanceEnvironment::FMU3InstanceEnvironment

    loggingOn::Bool
    instanceName::String
    continuousStatesChanged::fmi3Boolean

    t::fmi3Float64             # the system time
    t_offset::fmi3Float64      # time offset between simulation environment and FMU
    x::Union{Array{fmi3Float64, 1}, Nothing}   # the system states (or sometimes u)
    ẋ::Union{Array{fmi3Float64, 1}, Nothing}   # the system state derivative (or sometimes u̇)
    ẍ::Union{Array{fmi3Float64, 1}, Nothing}   # the system state second derivative
    #u::Array{fmi3Float64, 1}   # the system inputs
    #y::Array{fmi3Float64, 1}   # the system outputs
    #p::Array{fmi3Float64, 1}   # the system parameters
    z::Array{fmi3Float64, 1}   # the system event indicators
    z_prev::Union{Nothing, Array{fmi3Float64, 1}}   # the last system event indicators

    realValues::Dict

    x_vrs::Array{fmi3ValueReference, 1}   # the system state value references
    ẋ_vrs::Array{fmi3ValueReference, 1}   # the system state derivative value references
    u_vrs::Array{fmi3ValueReference, 1}   # the system input value references
    y_vrs::Array{fmi3ValueReference, 1}   # the system output value references
    p_vrs::Array{fmi3ValueReference, 1}   # the system parameter value references

    # deprecated
    jac_ẋy_x::Matrix{fmi3Float64}
    jac_ẋy_u::Matrix{fmi3Float64}
    jac_x::Array{fmi3Float64}
    jac_u::Union{Array{fmi3Float64}, Nothing}
    jac_t::fmi3Float64

    # linearization jacobians
    A::Union{Matrix{fmi3Float64}, Nothing}
    B::Union{Matrix{fmi3Float64}, Nothing}
    C::Union{Matrix{fmi3Float64}, Nothing}
    D::Union{Matrix{fmi3Float64}, Nothing}

    jacobianUpdate!             # function for a custom jacobian constructor (optimization)
    skipNextDoStep::Bool    # allows skipping the next `fmi3DoStep` like it is not called
    senseFunc::Symbol       # :auto, :full, :sample, :directionalDerivatives, :adjointDerivatives

    # custom

    rootsFound::Array{fmi3Int32}
    stateEvent::fmi3Boolean
    timeEvent::fmi3Boolean
    stepEvent::fmi3Boolean

    # constructor

    function FMU3Instance()
        inst = new()
        inst.state = fmi3InstanceStateInstantiated
        inst.t = -Inf
        inst.t_offset = 0.0

        inst.senseFunc = :auto

        inst.x = nothing
        inst.ẋ = nothing
        inst.ẍ = nothing
        inst.z_prev = nothing

        inst.realValues = Dict()
        inst.x_vrs = Array{fmi3ValueReference, 1}()
        inst.ẋ_vrs = Array{fmi3ValueReference, 1}()
        inst.u_vrs = Array{fmi3ValueReference, 1}()
        inst.y_vrs = Array{fmi3ValueReference, 1}()
        inst.p_vrs = Array{fmi3ValueReference, 1}()

        # initialize further variables
        inst.jac_x = Array{fmi3Float64, 1}()
        inst.jac_u = nothing
        inst.jac_t = -1.0
        inst.jac_ẋy_x = zeros(fmi3Float64, 0, 0)
        inst.jac_ẋy_u = zeros(fmi3Float64, 0, 0)
        inst.skipNextDoStep = false

        inst.A = nothing
        inst.B = nothing
        inst.C = nothing
        inst.D = nothing
        return inst
    end

    function FMU3Instance(compAddr, fmu)
        inst = FMU3Instance()
        inst.compAddr = compAddr
        inst.fmu = fmu
        return inst
    end
end

"""
Overload the Base.show() function for custom printing of the FMU2Component.
"""
Base.show(io::IO, c::FMU3Instance) = print(io,
"FMU:            $(c.fmu.modelDescription.modelName)
InstanceName:   $(isdefined(c, :instanceName) ? c.instanceName : "[not defined]")
Address:        $(c.compAddr)
State:          $(c.state)
Logging:        $(c.loggingOn)
FMU time:       $(c.t)
FMU states:     $(c.x)"
)

"""
=======
A mutable struct representing the excution configuration of a FMU.
For FMUs that have issues with calls like `fmi3Reset` or `fmi3FreeInstance`, this is pretty useful.
"""
mutable struct FMU3ExecutionConfiguration <: FMUExecutionConfig
    terminate::Bool     # call fmi3Terminate before every training step / simulation
    reset::Bool         # call fmi3Reset before every training step / simulation
    setup::Bool         # call setup functions before every training step / simulation
    instantiate::Bool   # call fmi3Instantiate before every training step / simulation
    freeInstance::Bool  # call fmi3FreeInstance after every training step / simulation

    loggingOn::Bool
    externalCallbacks::Bool

    handleStateEvents::Bool                 # handle state events during simulation/training
    handleTimeEvents::Bool                  # handle time events during simulation/training

    assertOnError::Bool                     # wheter an exception is thrown if a fmi3XXX-command fails (>= fmi3StatusError)
    assertOnWarning::Bool                   # wheter an exception is thrown if a fmi3XXX-command warns (>= fmi3StatusWarning)

    autoTimeShift::Bool                     # wheter to shift all time-related functions for simulation intervals not starting at 0.0

    sensealg                                # algorithm for sensitivity estimation over solve call
    useComponentShadow::Bool                # whether FMU outputs/derivatives/jacobians should be cached for frule/rrule (useful for ForwardDiff)
    rootSearchInterpolationPoints::UInt     # number of root search interpolation points
    inPlace::Bool                           # whether faster in-place-fx should be used
    useVectorCallbacks::Bool                # whether to vector (faster) or scalar (slower) callbacks

    maxNewDiscreteStateCalls::UInt          # max calls for fmi3NewDiscreteStates before throwing an exception

    function FMU3ExecutionConfiguration()
        inst = new()

        inst.terminate = true
        inst.reset = true
        inst.setup = true
        inst.instantiate = false
        inst.freeInstance = false

        inst.loggingOn = false
        inst.externalCallbacks = false

        inst.handleStateEvents = true
        inst.handleTimeEvents = true

        inst.assertOnError = false
        inst.assertOnWarning = false

        inst.autoTimeShift = true

        inst.sensealg = nothing # auto
        inst.useComponentShadow = false
        inst.rootSearchInterpolationPoints = 100
        inst.inPlace = true
        inst.useVectorCallbacks = true

        inst.maxNewDiscreteStateCalls = 100

        return inst
    end
end

# default for a "healthy" FMU - this is the fastetst
FMU3_EXECUTION_CONFIGURATION_RESET = FMU3ExecutionConfiguration()
FMU3_EXECUTION_CONFIGURATION_RESET.terminate = true
FMU3_EXECUTION_CONFIGURATION_RESET.reset = true
FMU3_EXECUTION_CONFIGURATION_RESET.instantiate = false
FMU3_EXECUTION_CONFIGURATION_RESET.freeInstance = false

# if your FMU has a problem with "fmi3Reset" - this is default
FMU3_EXECUTION_CONFIGURATION_NO_RESET = FMU3ExecutionConfiguration()
FMU3_EXECUTION_CONFIGURATION_NO_RESET.terminate = false
FMU3_EXECUTION_CONFIGURATION_NO_RESET.reset = false
FMU3_EXECUTION_CONFIGURATION_NO_RESET.instantiate = true
FMU3_EXECUTION_CONFIGURATION_NO_RESET.freeInstance = true

# if your FMU has a problem with "fmi3Reset" and "fmi3FreeInstance" - this is for weak FMUs (but slower)
FMU3_EXECUTION_CONFIGURATION_NO_FREEING = FMU3ExecutionConfiguration()
FMU3_EXECUTION_CONFIGURATION_NO_FREEING.terminate = false
FMU3_EXECUTION_CONFIGURATION_NO_FREEING.reset = false
FMU3_EXECUTION_CONFIGURATION_NO_FREEING.instantiate = true
FMU3_EXECUTION_CONFIGURATION_NO_FREEING.freeInstance = false

"""
ToDo
"""
struct FMU3Event <: FMUEvent
    t::Union{Float32, Float64}
    indicator::UInt

    x_left::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing}
    x_right::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing}

    function FMU3Event(t::Union{Float32, Float64}, indicator::UInt = 0, x_left::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing} = nothing, x_right::Union{Array{Float64, 1}, Array{Float32, 1}, Nothing} = nothing)
        inst = new(t, indicator, x_left, x_right)
        return inst
    end
end

"""
ToDo
"""
mutable struct FMU3Solution <: FMUSolution
    fmu                                             # FMU2
    success::Bool

    states                                          # ODESolution

    values
    valueReferences::Union{Array, Nothing}          # Array{fmi3ValueReference}

    events::Array{FMU2Event, 1}

    function FMU3Solution(fmu)
        inst = new()

        inst.fmu = fmu
        inst.success = false
        inst.states = nothing
        inst.values = nothing
        inst.valueReferences = nothing

        inst.events = []

        return inst
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.1. Header Files and Naming of Functions

The mutable struct representing an FMU in the FMI 3.0 Standard.
Also contains the paths to the FMU and ZIP folder as well als all the FMI 3.0 function pointers
"""
mutable struct FMU3 <: FMU
    modelName::String
    instanceName::String
    fmuResourceLocation::String

    modelDescription::fmi3ModelDescription

    type::fmi3Type
    instanceEnvironment::fmi3InstanceEnvironment
    instances::Array # {fmi3Instance}

    # c-functions
    cInstantiateModelExchange::Ptr{Cvoid}
    cInstantiateCoSimulation::Ptr{Cvoid}
    cInstantiateScheduledExecution::Ptr{Cvoid}

    cGetVersion::Ptr{Cvoid}
    cFreeInstance::Ptr{Cvoid}
    cSetDebugLogging::Ptr{Cvoid}
    cEnterConfigurationMode::Ptr{Cvoid}
    cExitConfigurationMode::Ptr{Cvoid}
    cEnterInitializationMode::Ptr{Cvoid}
    cExitInitializationMode::Ptr{Cvoid}
    cTerminate::Ptr{Cvoid}
    cReset::Ptr{Cvoid}
    cGetFloat32::Ptr{Cvoid}
    cSetFloat32::Ptr{Cvoid}
    cGetFloat64::Ptr{Cvoid}
    cSetFloat64::Ptr{Cvoid}
    cGetInt8::Ptr{Cvoid}
    cSetInt8::Ptr{Cvoid}
    cGetUInt8::Ptr{Cvoid}
    cSetUInt8::Ptr{Cvoid}
    cGetInt16::Ptr{Cvoid}
    cSetInt16::Ptr{Cvoid}
    cGetUInt16::Ptr{Cvoid}
    cSetUInt16::Ptr{Cvoid}
    cGetInt32::Ptr{Cvoid}
    cSetInt32::Ptr{Cvoid}
    cGetUInt32::Ptr{Cvoid}
    cSetUInt32::Ptr{Cvoid}
    cGetInt64::Ptr{Cvoid}
    cSetInt64::Ptr{Cvoid}
    cGetUInt64::Ptr{Cvoid}
    cSetUInt64::Ptr{Cvoid}
    cGetBoolean::Ptr{Cvoid}
    cSetBoolean::Ptr{Cvoid}
    cGetString::Ptr{Cvoid}
    cSetString::Ptr{Cvoid}
    cGetBinary::Ptr{Cvoid}
    cSetBinary::Ptr{Cvoid}
    cGetFMUState::Ptr{Cvoid}
    cSetFMUState::Ptr{Cvoid}
    cFreeFMUState::Ptr{Cvoid}
    cSerializedFMUStateSize::Ptr{Cvoid}
    cSerializeFMUState::Ptr{Cvoid}
    cDeSerializeFMUState::Ptr{Cvoid}
    cGetDirectionalDerivative::Ptr{Cvoid}
    cGetAdjointDerivative::Ptr{Cvoid}
    cEvaluateDiscreteStates::Ptr{Cvoid}
    cGetNumberOfVariableDependencies::Ptr{Cvoid}
    cGetVariableDependencies::Ptr{Cvoid}

    # Co Simulation function calls
    cGetOutputDerivatives::Ptr{Cvoid}
    cEnterStepMode::Ptr{Cvoid}
    cDoStep::Ptr{Cvoid}

    # Model Exchange function calls
    cGetNumberOfContinuousStates::Ptr{Cvoid}
    cGetNumberOfEventIndicators::Ptr{Cvoid}
    cGetContinuousStates::Ptr{Cvoid}
    cGetNominalsOfContinuousStates::Ptr{Cvoid}
    cEnterContinuousTimeMode::Ptr{Cvoid}
    cSetTime::Ptr{Cvoid}
    cSetContinuousStates::Ptr{Cvoid}
    cGetContinuousStateDerivatives::Ptr{Cvoid}
    cGetEventIndicators::Ptr{Cvoid}
    cCompletedIntegratorStep::Ptr{Cvoid}
    cEnterEventMode::Ptr{Cvoid}
    cUpdateDiscreteStates::Ptr{Cvoid}

    # Scheduled Execution function calls
    cSetIntervalDecimal::Ptr{Cvoid}
    cSetIntervalFraction::Ptr{Cvoid}
    cGetIntervalDecimal::Ptr{Cvoid}
    cGetIntervalFraction::Ptr{Cvoid}
    cGetShiftDecimal::Ptr{Cvoid}
    cGetShiftFraction::Ptr{Cvoid}
    cActivateModelPartition::Ptr{Cvoid}

    # paths of ziped and unziped FMU folders
    path::String
    binaryPath::String
    zipPath::String

    # execution configuration
    executionConfig::FMU3ExecutionConfiguration
    hasStateEvents::Union{Bool, Nothing}
    hasTimeEvents::Union{Bool, Nothing}

    # c-libraries
    libHandle::Ptr{Nothing}

    # START: experimental section (to FMIFlux.jl)
    dependencies::Matrix
    # linearization jacobians
    A::Matrix{fmi3Float64}
    B::Matrix{fmi3Float64}
    C::Matrix{fmi3Float64}
    D::Matrix{fmi3Float64}

    # END: experimental section

    # Constructor
    function FMU3()
        inst = new()
        inst.instances = []

        inst.hasStateEvents = nothing
        inst.hasTimeEvents = nothing

        inst.executionConfig = FMU3_EXECUTION_CONFIGURATION_NO_RESET
        return inst
    end
end

""" Overload the Base.show() function for custom printing of the FMU3"""
Base.show(io::IO, fmu::FMU3) = print(io,
"Model name:       $(fmu.modelName)
Instance name:     $(fmu.instanceName)
Model description: $(fmu.modelDescription)
Type:              $(fmu.type)
Instances:        $(fmu.instances)"
)

"""
Formats the fmi3Status into a String.
"""
function fmi3StatusToString(status::fmi3Status)
    if status == fmi3StatusOK
        return "OK"
    elseif status == fmi3StatusWarning
        return "Warning"
    elseif status == fmi3StatusDiscard
        return "Discard"
    elseif status == fmi3StatusError
        return "Error"
    elseif status == fmi3StatusFatal
        return "Fatal"
    else
        return "Unknown"
    end
end

"""
Formats the fmi3Status into a Integer.
"""
function fmi3StatusToString(status::Integer)
    if status == fmi3StatusOK
        return "OK"
    elseif status == fmi3StatusWarning
        return "Warning"
    elseif status == fmi3StatusDiscard
        return "Discard"
    elseif status == fmi3StatusError
        return "Error"
    elseif status == fmi3StatusFatal
        return "Fatal"
    else
        return "Unknown"
    end
end

function fmi3CausalityToString(c::fmi3Causality)
    if c == fmi3CausalityParameter
        return "parameter"
    elseif c == fmi3CausalityCalculatedParameter
        return "calculatedParameter"
    elseif c == fmi3CausalityInput
        return "input"
    elseif c == fmi3CausalityOutput
        return "output"
    elseif c == fmi3CausalityLocal
        return "local"
    elseif c == fmi3CausalityIndependent
        return "independent"
    elseif c == fmi3CausalityStructuralParameter
        return "structuralParameter"
    else
        @assert false "fmi3CausalityToString(...): Unknown causality."
    end
end

function fmi3StringToCausality(s::String)
    if s == "parameter"
        return fmi3CausalityParameter
    elseif s == "calculatedParameter"
        return fmi3CausalityCalculatedParameter
    elseif s == "input"
        return fmi3CausalityInput
    elseif s == "output"
        return fmi3CausalityOutput
    elseif s == "local"
        return fmi3CausalityLocal
    elseif s == "independent"
        return fmi3CausalityIndependent
    elseif s == "structuralParameter"
        return fmi3CausalityStructuralParameter
    else
        @assert false "fmi3StringToCausality($(s)): Unknown causality."
    end
end

function fmi3VariabilityToString(c::fmi3Variability)
    if c == fmi3VariabilityConstant
        return "constant"
    elseif c == fmi3VariabilityFixed
        return "fixed"
    elseif c == fmi3VariabilityTunable
        return "tunable"
    elseif c == fmi3VariabilityDiscrete
        return "discrete"
    elseif c == fmi3VariabilityContinuous
        return "continuous"
    else
        @assert false "fmi3VariabilityToString(...): Unknown variability."
    end
end

function fmi3StringToVariability(s::String)
    if s == "constant"
        return fmi3VariabilityConstant
    elseif s == "fixed"
        return fmi3VariabilityFixed
    elseif s == "tunable"
        return fmi3VariabilityTunable
    elseif s == "discrete"
        return fmi3VariabilityDiscrete
    elseif s == "continuous"
        return fmi3VariabilityContinuous
    else
        @assert false "fmi3StringToVariability($(s)): Unknown variability."
    end
end

function fmi3InitialToString(c::fmi3Initial)
    if c == fmi3InitialApprox
        return "approx"
    elseif c == fmi3InitialExact
        return "exact"
    elseif c == fmi3InitialCalculated
        return "calculated"
    else
        @assert false "fmi3InitialToString(...): Unknown initial."
    end
end

function fmi3StringToInitial(s::String)
    if s == "approx"
        return fmi3InitialApprox
    elseif s == "exact"
        return fmi3InitialExact
    elseif s == "calculated"
        return fmi3InitialCalculated
    else
        @assert false "fmi3StringToInitial($(s)): Unknown initial."
    end
end

function fmi3TypeToString(c::fmi3Type)
    if c == fmi3TypeCoSimulation
        return "coSimulation"
    elseif c == fmi3TypeModelExchange
        return "modelExchange"
    elseif c == fmi3TypeScheduledExecution
        return "scheduledExecution"
    else
        @assert false "fmi3TypeToString(...): Unknown type."
    end
end

function fmi3StringToType(s::String)
    if s == "coSimulation"
        return fmi3TypeCoSimulation
    elseif s == "modelExchange"
        return fmi3TypeModelExchange
    elseif s == "scheduledExecution"
        return fmi3TypeScheduledExecution
    else
        @assert false "fmi3StringToInitial($(s)): Unknown type."
    end
end

function fmi3IntervalQualifierToString(c::fmi3IntervalQualifier)
    if c == fmi3IntervalQualifierIntervalNotYetKnown
        return "intervalNotYetKnown"
    elseif c == fmi3IntervalQualifierIntervalUnchanged
        return "intervalUnchanged"
    elseif c == fmi3IntervalQualifierIntervalChanged
        return "intervalChanged"
    else
        @assert false "fmi3IntervalQualifierToString(...): Unknown intervalQualifier."
    end
end

function fmi3StringToIntervalQualifier(s::String)
    if s == "intervalNotYetKnown"
        return fmi3IntervalQualifierIntervalNotYetKnown
    elseif s == "intervalUnchanged"
        return fmi3IntervalQualifierIntervalUnchanged
    elseif s == "intervalChanged"
        return fmi3IntervalQualifierIntervalChanged
    else
        @assert false "fmi3StringToIntervalQualifier($(s)): Unknown intervalQualifier."
    end
end

function fmi3DependencyKindToString(c::fmi3DependencyKind)
    if c == fmi3DependencyKindIndependent
        return "independent"
    elseif c == fmi3DependencyKindConstant
        return "constant"
    elseif c == fmi3DependencyKindFixed
        return "fixed"
    elseif c == fmi3DependencyKindTunable
        return "tunable"
    elseif c == fmi3DependencyKindDiscrete
        return "discrete"
    elseif c == fmi3DependencyKindDependent
        return "dependent"
    else
        @assert false "fmi3DependencyKindToString(...): Unknown dependencyKind."
    end
end

function fmi3StringToDependencyKind(s::AbstractString)
    if s == "independent"
        return fmi3DependencyKindIndependent
    elseif s == "constant"
        return fmi3DependencyKindConstant
    elseif s == "fixed"
        return fmi3DependencyKindFixed
    elseif s == "tunable"
        return fmi3DependencyKindTunable
    elseif s == "discrete"
        return fmi3DependencyKindDiscrete
    elseif s == "dependent"
        return fmi3DependencyKindDependent
    else
        @assert false "fmi3StringToDependencyKind($(s)): Unknown dependencyKind."
    end
end

function fmi3VariableNamingConventionToString(c::fmi3VariableNamingConvention)
    if c == fmi3VariableNamingConventionFlat
        return "flat"
    elseif c == fmi3VariableNamingConventionStructured
        return "structured"
    else
        @assert false "fmi3VariableNamingConventionToString(...): Unknown variableNamingConvention."
    end
end

function fmi3StringToVariableNamingConvention(s::String)
    if s == "flat"
        return fmi3VariableNamingConventionFlat
    elseif s == "structured"
        return fmi3VariableNamingConventionStructured
    else
        @assert false "fmi3StringToVariableNamingConvention($(s)): Unknown variableNamingConvention."
    end
end


text:
"
A truncation region defines a subset of space to which the distribution is truncated. The basic operation supported is `intruncationregion()`, which returns true if a vector is inside the truncation region.
"
abstract type TruncationRegion end

"""
A state abstract object representing computed quantities of a truncated multivariate distribution.

Every concrete subtype should expose at least the following two fields.

- `n::Int` The length of the distribution
- `tp::Float64` The probability of falling in the truncation region for the non-truncated case.
- `tp_err::Float64` An estimate of the absolute error of the probability `tp`.

Other subtypes may expose.

- `μ::Vector{Float64}` The mean vector.
- `μ_err::Float64` An estimate of the relative error for the mean vector.
- `Σ::PDMat` The covariance matrix.
- `Σ_err::Float64` An estimate of the relative error for the covariance matrix.

Further subtypes may expose.
- `moment_dict::Dict{Vector{Int},Float64}` A dictionary mapping multivariate moment vectors to estimate quantities.
- `prob_dict` A dictionary mapping probability vectors to estimated quantities.
"""
abstract type TruncatedMvDistributionState end

"""
A truncated multi-variate distribution composed of a Multivariate Distribution, Truncation Region and a State object implementing computable state.
"""
struct TruncatedMvDistribution{D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    untruncated::D
    region::R
    state::S
end

"QQQQ"
function TruncatedMvDistribution{D,R,S}(d::D,r::R) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    TruncatedMvDistribution(d,r,S(d))
end

mutable struct TruncatedMvDistributionSecondOrderState <: TruncatedMvDistributionState
    n::Int              # dimension
    tp::Float64         # Truncation probability (probability under the truncated region)
    μ::Vector{Float64}  # Mean vector
    Σ::PDMat            # Covariance matrix
    tp_err::Float64     # Estimate of the error of truncation probability
    μ_err::Float64      # Estimate of the error of the mean vector
    Σ_err::Float64      # Estimate of the error fo the covariance matrix
    TruncatedMvDistributionSecondOrderState(d::MultivariateDistribution) = new(     length(d),
                                                                                NaN,
                                                                                Vector{Float64}(undef,0),
                                                                                PDMat(Array{Float64,2}(I,length(d),length(d))),
                                                                                Inf, Inf, Inf)
end

# function tp(d::TruncatedMvDistribution{D,R,TruncatedMvDistributionSecondOrderState})
#     where {D <: MultivariateDistribution, R <: TruncationRegion}
#     (d.state.tp, d.state.tp_err)
# end

# function mean(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
#     (d.state.μ, d.state.μ_err)
# end

# function cov(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
#     (d.s)
# end

    # function moment(d::TruncatedMvDistribution{D,R,S}, k::Vector{Int}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}

# function Base.show(io::IO, d::BoxTruncatedMvNormalRecursiveMomentsState)
#     println(io, "Box Truncated MvNormal")
#     println(io, "n = $(d.n)")
#     println(io, "μₑ = $(d.μₑ)" )
#     println(io, "Σₑ = $(d.Σₑ)" )
#     println(io, "α = $(alpha(d))")
#     println("Limits:")

#     for i in 1:d.n
#         println(io, "$i:\t ",(d.a[i],d.b[i]))
#     end

#     if d.momentsComputed
#         println("Moments:")
#         for k in keys(d.momentDict)
#             println(io, k, "\t", moment(d,k))
#         end
#         println(io,"mean:", mean(d))
#         println(io,"cov:", cov(d))
#     else
#         println("Moments not computed")
#     end
# end


text:
__precompile__(true)

module PlanarMaps

import Random,
       AsyPlots,
       DataStructures,
       SpecialFunctions

export NeighborCycle, CyclicRange, to,
       PlanarMap,
       Face,
       rotate, over,
       cw, ccw, ccwrange, insertccw, insertcw,
       isedge, add_edge!, delete_edge!,
       edges, face, faces, neighbors,
       interiorfaces, outerface,
       triangulation,
       RootedTree, root,
       isblue, isred, isgreen,
       vertices, interiorvertices,
       WoodedTriangulation, PartialWoodedTriangulation, shell!, shell,
       outerface, blueroot, redroot, greenroot,
       schnyderwood, allwoods,
       schnyder_blue, schnyder_red, schnyder_green, schnyder_coords,
       EmbeddedMap, draw, GraphPenSet,
       UWT

"""
    NeighborCycle(elements::Vector)

A circular permutation of distinct integers, representing
the (counterclockwise) neighbors of a vertex in a `PlanarMap`
"""
mutable struct NeighborCycle{T<:Integer}
    elements::Vector{T}
    lookup::Dict{T,T}
end

function remap!(N::NeighborCycle)
    N.lookup = Dict(map(reverse,enumerate(N.elements)));
end

function NeighborCycle(elements::Vector)
    NeighborCycle(elements,Dict(map(reverse,enumerate(elements))))
end

import Base.position
function position(N::NeighborCycle{T}, u::T) where T<:Integer
    N.lookup[u]
end

"""
    insertccw(N::NeighborCycle{T},u::T,v::T)

Return the `NeighborCycle` obtained by inserting the vertex
`v` counterclockwise of the vertex `u` in `N`
"""
function insertccw(N::NeighborCycle{T}, u::T, v::T) where T<:Integer
    nbs = copy(N.elements)
    k = position(N, u)
    insert!(nbs, k+1, v)
    NeighborCycle(nbs)
end

function insertcw(N::NeighborCycle{T}, u::T, v::T) where T<:Integer
    nbs = copy(N.elements)
    k = position(N, u)
    insert!(nbs, k, v)
    NeighborCycle(nbs)
end

"""
    PlanarMap(nbs::Vector{NeighborCycle})

A planar map with vertex set `1:length(nbs)` and
for which the neighbors of vertex `i` listed in
counterclockwise order are `nbs[i]`
"""
struct PlanarMap{T<:Integer}
    nbs::Vector{NeighborCycle{T}}
end

function PlanarMap(a::Vector{<:Vector{<:Integer}})
    return PlanarMap(map(NeighborCycle,a))
end

vertices(P::PlanarMap) = 1:length(P)

function neighbors(P::PlanarMap)
    P.nbs
end

function neighbors(P::PlanarMap,i::Integer)
    P.nbs[i]
end

"""
    edges(P::PlanarMap)

Return a vector containing both `(a,b)` and `(b,a)`
for every edge `(a,b)` in `P`
"""
function edges(P::PlanarMap)
    vcat([[(u,v) for v in N.elements]
                for (u,N) in enumerate(neighbors(P))]...)
end

struct Face{T<:Integer}
    elements::Vector{T}
end

Face(args...) = Face([args...])

function Base.show(io::IO,F::Face)
    print(io,"Face(")
    show(io,F.elements)
    print(io,")")
end

function isedge(P::PlanarMap{T},
                i::T,
                j::T) where T<:Integer
    j in neighbors(P,i)
end

"""
Find a face which has (i,j) as a diagonal
"""
function containing_face(P::PlanarMap{T},
                         i::T,
                         j::T) where T<:Integer
    for v in neighbors(P,i)
        F = Face(P,i,v)
        k = findfirst(isequal(j), F)
        if k ≠ 0
            return F
        end
    end
    error("No face found which has a diagonal ($i,$j)")
end

"""
    add_edge!(P::PlanarMap,i::Integer,j::Integer,F::Face)
    add_edge!(P::PlanarMap,i::Integer,j::Integer)

Add edge from `i` to `j` through the face `F`
in the planar map `P`. If `F` is not given, the
face discovered first in the `NeighborCycle` of `i`
will be used
"""
function add_edge!(P::PlanarMap{T},
                   i::T,
                   j::T) where T<:Integer
    if isedge(P,i,j)
        error("($i,$j) already an edge of $P")
    end
    for v in neighbors(P,i)
        F = face(P,i,v)
        k = findfirst(isequal(j), F)
        if !isa(k, Nothing)
            P.nbs[i] = insertccw(P.nbs[i], F[2], j)
            P.nbs[j] = insertccw(P.nbs[j], F[k+1], i)
            return
        end
    end
end

function add_edge!(P::PlanarMap{T},
                   i::T,
                   j::T,
                   F::Face) where T<:Integer
    if isedge(P,i,j)
        error("($i,$j) already an edge of $P")
    end
    k = findfirst(isequal(i), F)
    l = findfirst(isequal(j), F)
    if k == 0 || l == 0
        error("$i and $j are not both in $F")
    end
    P.nbs[i] = insertccw(P.nbs[i],F[k+1],j)
    P.nbs[j] = insertccw(P.nbs[j],F[l+1],i);
end

import Base: filter, filter!
Base.filter!(f::Function,N::NeighborCycle) = (Base.filter!(f,N.elements); remap!(N))
Base.filter(f::Function,N::NeighborCycle) = NeighborCycle(filter(f,N.elements))
Base.filter!(f::Function,F::Face) = Base.filter!(f,F.elements)
Base.filter(f::Function,F::Face) = Base.filter(f,F.elements)

function delete_edge!(P::PlanarMap{T},
                      i::T,
                      j::T) where T<:Integer
    if ~(i in neighbors(P,j) || j in neighbors(P,i))
        error("No edge to delete")
    end
    filter!(x->x≠j,P.nbs[i])
    filter!(x->x≠i,P.nbs[j])
    foreach(remap!,(P.nbs[i],P.nbs[j]))
end
"""
    OrientedPlanarMap(nbs::Vector{NeighborCycle{T}},
                      outgoing::Vector{Set{T}}) where T<:Integer

An oriented planar map with vertex set `1:length(nbs)` and
for which the neighbors of vertex `i` listed in counterclockwise order
are `nbs[i]` and for which the outneighbors are `outgoing[i]`
"""
struct OrientedPlanarMap{T<:Integer}
    nbs::Vector{NeighborCycle{T}}
    outneighbors::Vector{Set{T}}
end

function neighbors(O::OrientedPlanarMap)
    return O.nbs
end

function outneighbors(O::OrientedPlanarMap)
    return O.outneighbors
end

function PlanarMap(O::OrientedPlanarMap)
    return PlanarMap(neighbors(O))
end

#-----DISPLAY-----------------------------------
import Base.show
function show(io::IO,N::NeighborCycle)
    print(io,"NeighborCycle(")
    print(io,string(N.elements))
    print(io,")")
end

function show(io::IO,M::Union{PlanarMap,OrientedPlanarMap})
    s = isa(M,PlanarMap) ? "" : "Oriented"
    print(io,"$(s)PlanarMap(")
    print(io,string(length(M)))
    print(io,")")
end
#-------------------------------------------------

import Base.length
length(N::NeighborCycle) = length(N.elements)
length(P::PlanarMap) = length(P.nbs)
length(O::OrientedPlanarMap) = length(O.nbs)

#--- CYCLE FUNCTIONS -----------------------------
"""
    over(N::NeighborCycle,u::T,k::T) where T<: Integer

Return the element which is `k` positions
to the right of `u` in the neighbor cycle `N`
"""
function over(N::NeighborCycle,u::T,k::T) where T<: Integer
    return N.elements[mod1(position(N,u)+k,length(N))]
end

"""
    rotate(N::NeighborCycle,k::Integer)

Return a vector containing the elements of `N`
in order beginning with position k
"""
function rotate(N::NeighborCycle{T},u::T) where T
    k = position(N,u)
    return [N.elements[k:end]; N.elements[1:k-1]]
end

prv(N::NeighborCycle,k::Integer) = over(N,k,-1)
nxt(N::NeighborCycle,k::Integer) = over(N,k,1)

import Base.==
function ==(N::NeighborCycle,M::NeighborCycle)
    if length(N) == 0 || length(M) == 0
        return length(N) == 0 && length(M) == 0
    end
    a = first(N.elements)
    N.elements == rotate(M,a)
end

function pairs(N::NeighborCycle)
    return zip(N.elements,rotate(N,2))
end

cw(P::PlanarMap{T},u::T,v::T) where T<:Integer =
                        over(neighbors(P,u),v,-1)
ccw(P::PlanarMap{T},u::T,v::T) where T<:Integer =
                        over(neighbors(P,u),v,1)

struct CyclicRange
    start::Int64
    done::Int64
end

to(a,b) = CyclicRange(a,b)

length(F::Face) = length(F.elements)
import Base: getindex, iterate, lastindex
Base.getindex(F::Face,i::Integer) = Base.getindex(F.elements,mod1(i,length(F)))
Base.getindex(F::Face,r::UnitRange) = Base.getindex(F.elements,r)
function Base.getindex(F::Face,C::CyclicRange)
    a,b = mod1.((C.start,C.done),length(F))
    if a ≤ b
        F[a:b]
    else
        vcat(F[a:end], F[1:b])
    end
end
function Base.getindex(N::NeighborCycle,C::CyclicRange)
    a,b = map(mod1,(C.start, C.done),(length(N),length(N)))
    if a ≤ b
        N[a:b]
    else
        vcat(N[a:end], N[1:b])
    end
end
Base.getindex(N::PlanarMaps.NeighborCycle{Int64},
              r::UnitRange{Int64}) = Base.getindex(N.elements,r)
Base.iterate(F::Face) = Base.iterate(F.elements)
Base.iterate(F::Face,i::Integer) = Base.iterate(F.elements,i)
Base.lastindex(F::Face) = Base.lastindex(F.elements)
Base.getindex(N::NeighborCycle,i::Integer) = Base.getindex(N.elements,mod1(i,length(N)))
Base.iterate(N::NeighborCycle) = Base.iterate(N.elements)
Base.iterate(N::NeighborCycle,i::Integer) = Base.iterate(N.elements, i)
Base.lastindex(N::NeighborCycle) = Base.lastindex(N.elements)

"""
    face(P::PlanarMap,u::Int64,v::Int64,rotate=cw)

Return the face of the planar map P which includes the
edge (u,v) and finds each new edge by rotating in the
specified direction about the current vertex. Note
that this traces out the face in the direction *opposite*
to `rotate`
"""
function face(P::PlanarMap,u::Int64,v::Int64;rotate=cw)
    f = [u,v]
    while length(f) == 2 || (f[end-1],f[end]) ≠ (u,v)
        push!(f,rotate(P,f[end],f[end-1]))
    end
    return Face(f[1:end-2])
end

function edges(F::Face)
    ((F[i],F[mod1(i+1,length(F))]) for i=1:length(F))
end

function rotate(F::Face,k::Integer)
    return circshift(F.elements,-(k-1))
end

import Base.keys
keys(F::Face) = keys(F.elements)

function ==(F::Face,G::Face)
    a = first(F)
    places = findall(isequal(a), G)
    any(rotate(G,p) == F.elements for p in places)
end

struct FaceList
    interiorfaces::Vector{Face}
    outerface::Face
end

"""
    faces(P::PlanarMap)

Return a vector of faces of `P`
"""
function faces(RNG::Random.AbstractRNG,
               P::PlanarMap;
               outeredge=(1,neighbors(P,1)[1]),
               rotate=cw)
    rem_edges = Set(edges(P))
    int_faces = Face[]
    outerface = face(P,pop!(rem_edges,outeredge)...;rotate=rotate)
    for e in edges(outerface)
        delete!(rem_edges,e)
    end
    while length(rem_edges) > 0
        e = rand(RNG,rem_edges)
        f = face(P,pop!(rem_edges,e)...;rotate=rotate)
        for e in edges(f)
	        delete!(rem_edges,e)
        end
        push!(int_faces,f)
    end
    return FaceList(int_faces,outerface)
end

faces(P::PlanarMap;kwargs...) = faces(Random.MersenneTwister(0),P;kwargs...)

function allfaces(RNG::Random.AbstractRNG,
                  P::PlanarMap;
                  outeredge=(1,neighbors(P,1)[1]),
                  rotate=cw)
    L = faces(RNG,P;outeredge=outeredge,rotate=rotate)
    return vcat(interiorfaces(L),[outerface(L)])
end

allfaces(P::PlanarMap;kwargs...) =
            allfaces(Random.MersenneTwister(0),P;kwargs...)

function interiorfaces(L::FaceList)
    L.interiorfaces
end

function outerface(L::FaceList)
    L.outerface
end

function find_drawable_edge(P::PlanarMap,F::Face)
    if length(F) == 3
        error("No drawable edge; $F is a triangle")
    end
    for (i,v) in enumerate(F)
        j = findfirst(x-> x≠v && x∉neighbors(P,v), F)
        !isa(j,Nothing) && return (i,j)
    end
    error("No drawable edge found")
end

function poprand!(G::Random.AbstractRNG,collection)
    splice!(collection,rand(G,1:length(collection)))
end

function facewithcorners(P::PlanarMap{T},blue::T,green::T,red::T) where T<:Integer
    for v in neighbors(P,blue)
        F = face(P,blue,v)
        k = findfirst(isequal(green), F.elements)
        # look for a red *after* the green already found:
        l = k+findfirst(isequal(red), F.elements[k+1:end])
        if 1 < k < l
            return F
        end
    end
    error("No face containing $blue, $green, $red")
end

poprand!(collection) = poprand(Random.GLOBAL_RNG,collection)

"""
    triangulation(P::PlanarMap)

Find a simple plane triangulation obtained from P
by adding edges
"""
function triangulation(RNG::Random.AbstractRNG,
                       P::PlanarMap;
                       outeredge=(1,neighbors(P,1)[1]),
                       corners=nothing)
    T = deepcopy(P)
    if corners ≠ nothing
        F = PlanarMaps.facewithcorners(T,corners...)
        for (a,b) in edges(Face(corners...))
            if !isedge(T,a,b)
                add_edge!(T,a,b,F)
            end
        end
        outeredge=corners[1:2]
    end
    facelist = [F for F in allfaces(T;outeredge=outeredge) if length(F) > 3]
    while !isempty(facelist)
        F = poprand!(RNG,facelist)
        (i,j) = find_drawable_edge(T,F)
        T.nbs[F[i]] = insertccw(T.nbs[F[i]],F[i+1],F[j])
        T.nbs[F[j]] = insertccw(T.nbs[F[j]],F[j+1],F[i])
        newfaces = [f for f in map(Face,(F[to(i,j)],F[to(j,i)])) if length(f) > 3]
        append!(facelist,newfaces)
    end
    return T
end

function triangulation(P::PlanarMap;kwargs...)
    triangulation(Random.MersenneTwister(0),P;kwargs...)
end

import Base.in
Base.in(a,N::NeighborCycle) = Base.in(a,keys(N.lookup))

pairs(F::Face) = zip(F.elements,circshift(F.elements,-1))

include("woodedtriangulations.jl")
include("uniformwoodedtriangulations.jl")
include("graphing.jl")

end # module


text:
# ####################### Parameters #############################################
# # Abstract type
# import Base.NamedTuple as Parameter

# # Operators
# ⊕(p::Parameter, q::Parameter) = merge(p, q)
# ≅(p::Tuple, q::Tuple) = length(p) == length(q) && all([size(pp) == size(qq) for (pp, qq) in zip(p, q)])
# ≅(x::NamedTuple{N,T}, y::NamedTuple{N2,T2}) where {N,T,N2,T2} =
#   all( [key1 == key2 for (key1, key2) in zip(keys(x), keys(y))] ) && (values(x) ≅ values(y))

# +(p::Tuple, q::Tuple) = (length(p) == length(q)) ? [v1 + v2 for (v1, v2) in zip(values(p), values(q))] : AssertionError("p and q must have the same length and sizes.")
# +(p::Parameter, q::Parameter) = ( p ≅ q ) ? Parameter{Tuple(keys(p))}( values(p) + values(q) ) : AssertionError("p and q must have the same keys and sizes.")

# -(p::Tuple, q::Tuple) = (length(p) == length(q)) ? [v1 - v2 for (v1, v2) in zip(values(p), values(q))] : AssertionError("p and q must have the same length and sizes.")
# -(p::Parameter, q::Parameter) = ( p ≅ q ) ? Parameter{Tuple(keys(p))}( values(p) - values(q) ) : AssertionError("p and q must have the same keys and sizes.")

# *(a::Real, p::Parameter) = Parameter{Tuple(keys(p))}( a .* values(p) )

# ## Convert parameter (named tuple) to vector for use in fitting
# function parameter_to_vec( p :: Parameter )
#     l = length(p)
#     v = Vector{Float64}()
#     for (k, vals) in zip(keys(p), p)
#         if typeof(vals) <: AbstractFloat
#             append!(v, vals)
#         else
#             append!(v, vec(vals))
#         end
#     end
#     return v
# end

# ## Convert vector back to parameter for use after fitting
# function vec_to_parameter!(p::Parameter, v::Vector{Float64})
#     i = 1
#     j = 1
#     t = Vector{Any}(undef, length(p))
#     for val in p
#         l = length(val)
#         if l == 1
#             t[j] =  v[i]
#             i += 1
#         else
#             s = size(val)
#             t[j] = reshape( v[i:(i+l-1)], s)
#             i += l
#         end
#         j += 1
#     end
#     p = Parameter{keys(p)}(t)
#     return p
# end


text:
function surf_exch(ebm::EBM, cn::Constants, Ta, Ua)

    zTs = ebm.zT
    if (ebm.zvar)
        zTs = ebm.zT - sum(ebm.Ds)
        zTs = max(zTs, 1.0)
    end

    # Neutral exchange coefficients
    z0h = 0.1 * ebm.z0
    CD = (cn.vkman / log(ebm.zU / ebm.z0))^2.0
    ebm.CH = cn.vkman^2.0 / (log(ebm.zU / ebm.z0) * log(zTs / z0h))

    # Stability correction (Louis et al. 1982, quoted by Beljaars 1992)
    if (ebm.em == 1)
        RiB = cn.g * (Ta - ebm.Tsurf) * ebm.zU^2.0 / (zTs * Ta * Ua^2.0)
        if (RiB > 0.0)
            fh = 1.0 / (1.0 + 3.0 * ebm.bstb * RiB * sqrt(1.0 + ebm.bstb * RiB))
        else
            fh = 1.0 - 3.0 * ebm.bstb * RiB / (1.0 + 3.0 * ebm.bstb^2.0 * CD * sqrt(-RiB * ebm.zU / ebm.z0))
        end
        ebm.CH = fh * ebm.CH
    end

    return nothing

end


text:
#!/usr/bin/env julia
using RobotOS
@rosimport mavs_control.msg: Control
rostypegen()
using mavs_control.msg

#http://wiki.ros.org/message_filters#Time_Synchronizer

# Set up service to get Gazebo link state
const get_link_state = ServiceProxy("/gazebo/get_link_state",GetLinkState)
println("Waiting for 'gazebo/get_link_state' service...")
wait_for_service("gazebo/get_link_state")

function main()
  init_node("ros_ex")
  const pub = Publisher{Control}("/mavs/optimal_control", queue_size=10)
  loop_rate = Rate(5.0)
  while !is_shutdown()
    msg = Control()
    msg.sa = [1.0, 2.0]
    msg.vx = [ 1.9, 3.0]
    msg.x = [ 1.9, 3.0]
    msg.y = [ 1.9, 3.0]
    msg.psi = [ 1.9, 3.0]
    publish(pub,msg)
    rossleep(loop_rate)
  end
end

main()


text:
# imports
using FMI
using FMIFlux
using FMIZoo
using Flux
using DifferentialEquations: Tsit5
import Plots

# set seed
import Random
Random.seed!(42);

tStart = 0.0
tStep = 0.01
tStop = 5.0
tSave = collect(tStart:tStep:tStop)

realFMU = fmiLoad("SpringFrictionPendulum1D", "Dymola", "2022x")
fmiInfo(realFMU)

initStates = ["s0", "v0"]
x₀ = [0.5, 0.0]
params = Dict(zip(initStates, x₀))
vrs = ["mass.s", "mass.v", "mass.a", "mass.f"]

realSimData = fmiSimulate(realFMU, tStart, tStop; parameters=params, recordValues=vrs, saveat=tSave)
fmiPlot(realSimData)

velReal = fmi2GetSolutionValue(realSimData, "mass.v")
posReal = fmi2GetSolutionValue(realSimData, "mass.s")
fmiUnload(realFMU)
simpleFMU = fmiLoad("SpringPendulum1D", "Dymola", "2022x")
fmiInfo(simpleFMU)

vrs = ["mass.s", "mass.v", "mass.a"]
simpleSimData = fmiSimulate(simpleFMU, tStart, tStop; recordValues=vrs, saveat=tSave, reset=false)
fmiPlot(simpleSimData)

velSimple = fmi2GetSolutionValue(simpleSimData, "mass.v")
posSimple = fmi2GetSolutionValue(simpleSimData, "mass.s")

# loss function for training
function lossSum()
    global posReal
    solution = neuralFMU(x₀, tStart)

    posNet = fmi2GetSolutionState(solution, 1; isIndex=true)
    # velNet = fmi2GetSolutionState(solution, 2; isIndex=true)

    Flux.Losses.mse(posReal, posNet) #+ Flux.Losses.mse(velReal, velNet)
end

# callback function for training
global counter = 0
function callb()
    global counter += 1
    if counter % 20 == 1
        avgLoss = lossSum()
        @info "Loss [$counter]: $(round(avgLoss, digits=5))   Avg displacement in data: $(round(sqrt(avgLoss), digits=5))"
    end
end

# NeuralFMU setup
numStates = fmiGetNumberOfStates(simpleFMU)

net = Chain(inputs -> fmiEvaluateME(simpleFMU, inputs),
            Dense(numStates, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))

neuralFMU = ME_NeuralFMU(simpleFMU, net, (tStart, tStop), Tsit5(); saveat=tSave);

solutionBefore = neuralFMU(x₀, tStart)
fmiPlot(solutionBefore)

# train
paramsNet = Flux.params(neuralFMU)

optim = ADAM()
Flux.train!(lossSum, paramsNet, Iterators.repeated((), 300), optim; cb=callb)

# plot results mass.s
solutionAfter = neuralFMU(x₀, tStart)

fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
                 xtickfontsize=12, ytickfontsize=12,
                 xguidefontsize=12, yguidefontsize=12,
                 legendfontsize=8, legend=:topright)

posNeuralFMU = fmi2GetSolutionState(solutionAfter, 1; isIndex=true)

Plots.plot!(fig, tSave, posSimple, label="SimpleFMU", linewidth=2)
Plots.plot!(fig, tSave, posReal, label="RealFMU", linewidth=2)
Plots.plot!(fig, tSave, posNeuralFMU, label="NeuralFMU (300 epochs)", linewidth=2)
fig

Flux.train!(lossSum, paramsNet, Iterators.repeated((), 700), optim; cb=callb)
# plot results mass.s
solutionAfter = neuralFMU(x₀, tStart)
posNeuralFMU = fmi2GetSolutionState(solutionAfter, 1; isIndex=true)
Plots.plot!(fig, tSave, posNeuralFMU, label="NeuralFMU (1000 epochs)", linewidth=2)
fig

fmiUnload(simpleFMU)