text:
using ConstructionBase: constructorof
using EquationsOfState: Parameters, EquationOfState
using StructHelpers: @batteries

export Murnaghan,
    Murnaghan1st,
    BirchMurnaghan,
    BirchMurnaghan2nd,
    BirchMurnaghan3rd,
    BirchMurnaghan4th,
    PoirierTarantola,
    PoirierTarantola2nd,
    PoirierTarantola3rd,
    # PoirierTarantola4th,
    Vinet,
    EnergyEquation,
    PressureEquation,
    BulkModulusEquation

abstract type FiniteStrainParameters{N,T} <: Parameters{T} end
abstract type BirchMurnaghan{N,T} <: FiniteStrainParameters{N,T} end
abstract type PoirierTarantola{N,T} <: FiniteStrainParameters{N,T} end
abstract type Murnaghan{T} <: Parameters{T} end
"""
    Murnaghan1st(v0, b0, b′0, e0=zero(v0 * b0))

Create a Murnaghan first order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct Murnaghan1st{T} <: Murnaghan{T}
    v0::T
    b0::T
    b′0::T
    e0::T
    Murnaghan1st{T}(v0, b0, b′0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, b′0, e0)
end
"""
    Murnaghan2nd(v0, b0, b′0, b″0, e0=zero(v0 * b0))

Create a Murnaghan second order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `b″0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct Murnaghan2nd{T} <: Murnaghan{T}
    v0::T
    b0::T
    b′0::T
    b″0::T
    e0::T
    Murnaghan2nd{T}(v0, b0, b′0, b″0, e0 = zero(v0 * b0)) where {T} =
        new(v0, b0, b′0, b″0, e0)
end
"""
    BirchMurnaghan2nd(v0, b0, e0=zero(v0 * b0))

Create a Birch–Murnaghan second order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct BirchMurnaghan2nd{T} <: BirchMurnaghan{2,T}
    v0::T
    b0::T
    e0::T
    BirchMurnaghan2nd{T}(v0, b0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, e0)
end
"""
    BirchMurnaghan3rd(v0, b0, b′0, e0=zero(v0 * b0))

Create a Birch–Murnaghan third order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.

!!! note
    The third-order equation (Equation (22)) becomes identical to the second-order equation
    when ``b′0 = 4`` (not ``0``!).
"""
struct BirchMurnaghan3rd{T} <: BirchMurnaghan{3,T}
    v0::T
    b0::T
    b′0::T
    e0::T
    BirchMurnaghan3rd{T}(v0, b0, b′0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, b′0, e0)
end
"""
    BirchMurnaghan4th(v0, b0, b′0, b″0, e0=zero(v0 * b0))

Create a Birch–Murnaghan fourth order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `b″0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.

!!! note
    The fourth-order equation becomes identical to the third-order equation when
    ```math
    b″0 = -\\frac{ 1 }{ 9b0 } (9b′0^2 - 63b′0 + 143).
    ```
"""
struct BirchMurnaghan4th{T} <: BirchMurnaghan{4,T}
    v0::T
    b0::T
    b′0::T
    b″0::T
    e0::T
    BirchMurnaghan4th{T}(v0, b0, b′0, b″0, e0 = zero(v0 * b0)) where {T} =
        new(v0, b0, b′0, b″0, e0)
end
"""
    PoirierTarantola2nd(v0, b0, e0=zero(v0 * b0))

Create a Poirier–Tarantola second order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct PoirierTarantola2nd{T} <: PoirierTarantola{2,T}
    v0::T
    b0::T
    e0::T
    PoirierTarantola2nd{T}(v0, b0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, e0)
end
"""
    PoirierTarantola3rd(v0, b0, b′0, e0=zero(v0 * b0))

Create a Poirier–Tarantola third order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct PoirierTarantola3rd{T} <: PoirierTarantola{3,T}
    v0::T
    b0::T
    b′0::T
    e0::T
    PoirierTarantola3rd{T}(v0, b0, b′0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, b′0, e0)
end
"""
    PoirierTarantola4th(v0, b0, b′0, b″0, e0=zero(v0 * b0))

Create a Poirier–Tarantola fourth order equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `b″0`: the second-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct PoirierTarantola4th{T} <: PoirierTarantola{4,T}
    v0::T
    b0::T
    b′0::T
    b″0::T
    e0::T
    PoirierTarantola4th{T}(v0, b0, b′0, b″0, e0 = zero(v0 * b0)) where {T} =
        new(v0, b0, b′0, b″0, e0)
end
"""
    Vinet(v0, b0, b′0, e0=zero(v0 * b0))

Create a Vinet equation of state.

This equation of state can have units. The units are specified in
[`Unitful.jl`](https://github.com/PainterQubits/Unitful.jl)'s `@u_str` style.

# Arguments
- `v0`: the volume of solid at zero pressure.
- `b0`: the bulk modulus of solid at zero pressure.
- `b′0`: the first-order pressure-derivative bulk modulus of solid at zero pressure.
- `e0`: the energy of solid at zero pressure.
"""
struct Vinet{T} <: Parameters{T}
    v0::T
    b0::T
    b′0::T
    e0::T
    Vinet{T}(v0, b0, b′0, e0 = zero(v0 * b0)) where {T} = new(v0, b0, b′0, e0)
end
struct AntonSchmidt{T} <: Parameters{T}
    v0::T
    b0::T
    b′0::T
    e∞::T
    AntonSchmidt{T}(v0, b0, b′0, e∞ = zero(v0 * b0)) where {T} = new(v0, b0, b′0, e∞)
end
struct Holzapfel{Z,T} <: Parameters{T}
    v0::T
    b0::T
    b′0::T
    e0::T
    function Holzapfel{Z,T}(v0, b0, b′0, e0 = zero(v0 * b0)) where {Z,T}
        @assert 1 <= Z <= 118 "elements are between 1 and 118!"
        return new(v0, b0, b′0, e0)
    end
end

@batteries Murnaghan1st eq = true hash = true
@batteries Murnaghan2nd eq = true hash = true
@batteries BirchMurnaghan2nd eq = true hash = true
@batteries BirchMurnaghan3rd eq = true hash = true
@batteries BirchMurnaghan4th eq = true hash = true
@batteries PoirierTarantola2nd eq = true hash = true
@batteries PoirierTarantola3rd eq = true hash = true
@batteries PoirierTarantola4th eq = true hash = true
@batteries Vinet eq = true hash = true
@batteries AntonSchmidt eq = true hash = true
@batteries Holzapfel eq = true hash = true

function (::Type{T})(args...) where {T<:Parameters}
    E = Base.promote_typeof(args...)
    return constructorof(T){E}(args...)  # Cannot use `T.(args...)`! For `AbstractQuantity` they will fail!
end
function Murnaghan(args...)
    N = length(args)
    if N == 4
        return Murnaghan1st(args...)
    elseif N == 5
        return Murnaghan2nd(args...)
    else
        throw(ArgumentError("unknown number of arguments $N."))
    end
end
"""
    BirchMurnaghan(args...)

Construct a `BirchMurnaghan` based on the length of arguments, where `e0` must be provided.

See also: [`BirchMurnaghan2nd`](@ref), [`BirchMurnaghan3rd`](@ref), [`BirchMurnaghan4th`](@ref)
"""
function BirchMurnaghan(args...)
    N = length(args)
    if N == 3
        return BirchMurnaghan2nd(args...)
    elseif N == 4
        return BirchMurnaghan3rd(args...)
    elseif N == 5
        return BirchMurnaghan4th(args...)
    else
        throw(ArgumentError("unknown number of arguments $N."))
    end
end
"""
    PoirierTarantola(args...)

Construct a `PoirierTarantola` based on the length of arguments, where `e0` must be provided.

See also: [`PoirierTarantola2nd`](@ref), [`PoirierTarantola3rd`](@ref), [`PoirierTarantola4th`](@ref)
"""
function PoirierTarantola(args...)
    N = length(args)
    if N == 3
        return PoirierTarantola2nd(args...)
    elseif N == 4
        return PoirierTarantola3rd(args...)
    elseif N == 5
        return PoirierTarantola4th(args...)
    else
        throw(ArgumentError("unknown number of arguments $N."))
    end
end

"""
    EquationOfStateOfSolids{T<:Parameters}

Represent an equation of state of solids.
"""
abstract type EquationOfStateOfSolids{T} <: EquationOfState{T} end
"""
    EnergyEquation{T} <: EquationOfStateOfSolids{T}
    EnergyEquation(parameters::Parameters)

Construct an equation of state which evaluates the energy of the given `parameters`.
"""
struct EnergyEquation{T} <: EquationOfStateOfSolids{T}
    param::T
end
EnergyEquation(eos::EquationOfStateOfSolids) = EnergyEquation(getparam(eos))
"""
    PressureEquation{T} <: EquationOfStateOfSolids{T}
    PressureEquation(parameters::Parameters)

Construct an equation of state which evaluates the pressure of the given `parameters`.
"""
struct PressureEquation{T} <: EquationOfStateOfSolids{T}
    param::T
end
PressureEquation(eos::EquationOfStateOfSolids) = PressureEquation(getparam(eos))
"""
    BulkModulusEquation{T} <: EquationOfStateOfSolids{T}
    BulkModulusEquation(parameters::Parameters)

Construct an equation of state which evaluates the bulk modulus of the given `parameters`.
"""
struct BulkModulusEquation{T} <: EquationOfStateOfSolids{T}
    param::T
end
BulkModulusEquation(eos::EquationOfStateOfSolids) = BulkModulusEquation(getparam(eos))


text:
###############################################################################
#
#   fmpq_rel_series.jl : Power series over flint fmpq rationals
#
###############################################################################

export fmpq_rel_series, FmpqRelSeriesRing

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

function O(a::fmpq_rel_series)
   val = pol_length(a) + valuation(a) - 1
   val < 0 && throw(DomainError())
   z = fmpq_rel_series(Array{fmpq}(0), 0, val, val)
   z.parent = parent(a)
   return z
end

elem_type(::Type{FmpqRelSeriesRing}) = fmpq_rel_series

parent_type(::Type{fmpq_rel_series}) = FmpqRelSeriesRing

base_ring(R::FmpqRelSeriesRing) = R.base_ring

var(a::FmpqRelSeriesRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

max_precision(R::FmpqRelSeriesRing) = R.prec_max

function normalise(a::fmpq_rel_series, len::Int)
   if len > 0
      c = fmpq()
      ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq_rel_series}, Int), c, a, len - 1)
   end
   while len > 0 && iszero(c)
      len -= 1
      if len > 0
         ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
            (Ref{fmpq}, Ref{fmpq_rel_series}, Int), c, a, len - 1)
      end
   end
   return len
end

function pol_length(x::fmpq_rel_series)
   return ccall((:fmpq_poly_length, :libflint), Int, (Ref{fmpq_rel_series},), x)
end

precision(x::fmpq_rel_series) = x.prec

function polcoeff(x::fmpq_rel_series, n::Int)
   if n < 0
      return fmpq(0)
   end
   z = fmpq()
   ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq_rel_series}, Int), z, x, n)
   return z
end

zero(R::FmpqRelSeriesRing) = R(0)

one(R::FmpqRelSeriesRing) = R(1)

function gen(R::FmpqRelSeriesRing)
   z = fmpq_rel_series([fmpq(1)], 1, max_precision(R) + 1, 1)
   z.parent = R
   return z
end

function deepcopy_internal(a::fmpq_rel_series, dict::ObjectIdDict)
   z = fmpq_rel_series(a)
   z.prec = a.prec
   z.val = a.val
   z.parent = parent(a)
   return z
end

function renormalize!(z::fmpq_rel_series)
   i = 0
   zlen = pol_length(z)
   zval = valuation(z)
   zprec = precision(z)
   while i < zlen && iszero(polcoeff(z, i))
      i += 1
   end
   z.prec = zprec
   if i == zlen
      z.val = zprec
   else
      z.val = zval + i
      ccall((:fmpq_poly_shift_right, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int), z, z, i)
   end
   return nothing
end

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, a::FmpqRelSeriesRing)
   print(io, "Univariate power series ring in ", var(a), " over ")
   show(io, base_ring(a))
end

show_minus_one(::Type{fmpq_rel_series}) = show_minus_one(fmpq)

###############################################################################
#
#   Unary operators
#
###############################################################################

function -(x::fmpq_rel_series)
   z = parent(x)()
   ccall((:fmpq_poly_neg, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}),
               z, x)
   z.prec = x.prec
   z.val = x.val
   return z
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::fmpq_rel_series, b::fmpq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function -(a::fmpq_rel_series, b::fmpq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   lenz = max(lena, lenb)
   z = parent(a)()
   if a.val < b.val
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpq_poly_neg, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}), z, z)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, a, max(0, lenz - a.val + b.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, z, a.val - b.val)
      ccall((:fmpq_poly_sub_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpq_poly_sub_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, b, lenz)
   end
   z.prec = prec
   z.val = val
   renormalize!(z)
   return z
end

function *(a::fmpq_rel_series, b::fmpq_rel_series)
   check_parent(a, b)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z = parent(a)()
   z.val = a.val + b.val
   z.prec = prec + z.val
   if lena == 0 || lenb == 0
      return z
   end
   lenz = min(lena + lenb - 1, prec)
   ccall((:fmpq_poly_mullow, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, b, lenz)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpq_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpq_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, y, x)
   return z
end

*(x::fmpq_rel_series, y::Int) = y * x

function *(x::fmpz, y::fmpq_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpq_poly_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpz}),
               z, y, x)
   return z
end

function *(x::fmpq, y::fmpq_rel_series)
   z = parent(y)()
   z.prec = y.prec
   z.val = y.val
   ccall((:fmpq_poly_scalar_mul_fmpq, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq}),
               z, y, x)
   return z
end

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

*(x::fmpq_rel_series, y::fmpq) = y * x

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

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

*(x::fmpq_rel_series, y::Rational) = x * fmpq(y)

*(x::Rational, y::fmpq_rel_series) = fmpq(x) * y

+(x::fmpq_rel_series, y::Integer) = x + fmpz(y)

+(x::Integer, y::fmpq_rel_series) = fmpz(x) + y

+(x::fmpq_rel_series, y::Rational) = x + fmpq(y)

+(x::Rational, y::fmpq_rel_series) = fmpq(x) + y

-(x::fmpq_rel_series, y::Integer) = x - fmpz(y)

-(x::Integer, y::fmpq_rel_series) = fmpz(x) - y

-(x::fmpq_rel_series, y::Rational) = x - fmpq(y)

-(x::Rational, y::fmpq_rel_series) = fmpq(x) - y

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpq_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   z = fmpq_rel_series(x)
   z.prec = x.prec + len
   z.val = x.val + len
   z.parent = parent(x)
   return z
end

function shift_right(x::fmpq_rel_series, len::Int)
   len < 0 && throw(DomainError())
   xlen = pol_length(x)
   xval = valuation(x)
   z = parent(x)()
   if len >= xlen + xval
      z.prec = max(0, x.prec - len)
      z.val = max(0, x.prec - len)
   else
      z.prec = max(0, x.prec - len)
      z.val = max(0, xval - len)
      zlen = min(xlen + xval - len, xlen)
      ccall((:fmpq_poly_shift_right, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, x, xlen - zlen)
      renormalize!(z)
   end
   return z
end

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(x::fmpq_rel_series, prec::Int)
   prec < 0 && throw(DomainError())
   xlen = pol_length(x)
   xprec = precision(x)
   xval = valuation(x)
   if xprec + xval <= prec
      return x
   end
   z = parent(x)()
   z.prec = prec
   if prec <= xval
      z = parent(x)()
      z.val = prec
      z.prec = prec
   else
      z.val = xval
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, x, min(prec - xval, xlen))
   end
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::fmpq_rel_series, b::Int)
   b < 0 && throw(DomainError())
   if isgen(a)
      z = parent(a)()
      z = setcoeff!(z, 0, fmpq(1))
      z.prec = a.prec + b - 1
      z.val = b
   elseif pol_length(a) == 0
      z = parent(a)()
      z.prec = b*valuation(a)
      z.val = b*valuation(a)
   elseif pol_length(a) == 1
      return parent(a)([polcoeff(a, 0)^b], 1,
                           (b - 1)*valuation(a) + precision(a), b*valuation(a))
   elseif b == 0
      return one(parent(a))
   else
      bit = ~((~UInt(0)) >> 1)
      while (UInt(bit) & b) == 0
         bit >>= 1
      end
      z = a
      bit >>= 1
      while bit != 0
         z = z*z
         if (UInt(bit) & b) != 0
            z *= a
         end
         bit >>= 1
      end
   end
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fmpq_rel_series, y::fmpq_rel_series)
   check_parent(x, y)
   prec = min(x.prec, y.prec)
   if prec <= x.val && prec <= y.val
      return true
   end
   if x.val != y.val
      return false
   end
   xlen = normalise(x, min(pol_length(x), prec - x.val))
   ylen = normalise(y, min(pol_length(y), prec - y.val))
   if xlen != ylen
      return false
   end
   return Bool(ccall((:fmpq_poly_equal_trunc, :libflint), Cint,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               x, y, xlen))
end

function isequal(x::fmpq_rel_series, y::fmpq_rel_series)
   if parent(x) != parent(y)
      return false
   end
   if x.prec != y.prec || x.val != y.val || pol_length(x) != pol_length(y)
      return false
   end
   return Bool(ccall((:fmpq_poly_equal, :libflint), Cint,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               x, y, pol_length(x)))
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpq_rel_series, y::fmpq)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = fmpq()
         ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
                       (Ref{fmpq}, Ref{fmpq_rel_series}, Int), z, x, 0)
         return ccall((:fmpq_equal, :libflint), Bool,
               (Ref{fmpq}, Ref{fmpq}, Int), z, y, 0)
      else
         return false
      end
   else
      return iszero(y)
   end
end

function ==(x::fmpq_rel_series, y::fmpz)
   if precision(x) == 0
      return true
   elseif pol_length(x) > 1
      return false
   elseif pol_length(x) == 1
      if x.val == 0
         z = fmpq()
         ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
                       (Ref{fmpq}, Ref{fmpq_rel_series}, Int), z, x, 0)
         return isone(denominator(z)) && ccall((:fmpz_equal, :libflint), Bool,
               (Ref{fmpz}, Ref{fmpz}, Int), numerator(z), y, 0)
      else
         return false
      end
   else
      return iszero(y)
   end
end

==(x::fmpz, y::fmpq_rel_series) = y == x

==(x::fmpq, y::fmpq_rel_series) = y == x

==(x::fmpq_rel_series, y::Integer) = x == fmpz(y)

==(x::Integer, y::fmpq_rel_series) = y == x

==(x::fmpq_rel_series, y::Rational{T}) where T <: Union{Int, BigInt} = x == fmpq(y)

==(x::Rational{T}, y::fmpq_rel_series) where T <: Union{Int, BigInt} = y == x

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpq_rel_series, y::fmpq_rel_series)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   yval = valuation(y)
   xval = valuation(x)
   if yval != 0
      if xval >= yval
         x = shift_right(x, yval)
         y = shift_right(y, yval)
      end
   end
   !isunit(y) && error("Unable to invert power series")
   prec = min(x.prec - x.val, y.prec - y.val)
   z = parent(x)()
   z.val = xval - yval
   z.prec = prec + z.val
   if prec != 0
      ccall((:fmpq_poly_div_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, x, y, prec)
   end
   return z
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpq_rel_series, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpq_poly_scalar_divexact_si, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, x, y)
   return z
end

function divexact(x::fmpq_rel_series, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpq_poly_scalar_divexact_fmpz, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpz}),
               z, x, y)
   return z
end

function divexact(x::fmpq_rel_series, y::fmpq)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   z.prec = x.prec
   z.prec = x.prec
   z.val = x.val
   ccall((:fmpq_poly_scalar_div_fmpq, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq}),
               z, x, y)
   return z
end

divexact(x::fmpq_rel_series, y::Integer) = divexact(x, fmpz(y))

divexact(x::fmpq_rel_series, y::Rational{T}) where T <: Union{Int, BigInt} = divexact(x, fmpq(y))

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(a::fmpq_rel_series)
   iszero(a) && throw(DivideError())
   !isunit(a) && error("Unable to invert power series")
   ainv = parent(a)()
   ainv.prec = a.prec
   ainv.val = 0
   ccall((:fmpq_poly_inv_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               ainv, a, a.prec)
   return ainv
end

###############################################################################
#
#   Special functions
#
###############################################################################

function Base.exp(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in exp")
   if pol_length(a) + valuation(a) == 0 || a.prec == 1
      return parent(a)([fmpq(1)], 1, a.prec, 0)
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_exp_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   log(a::fmpq_rel_series)
> Return log$(a)$. Requires the constant term to be one.
"""
function log(a::fmpq_rel_series)
   (a.val != 0 || coeff(a, 0) != 1) && error("Constant term not one in log")
   if pol_length(a) + valuation(a) == 1 || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_log_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.prec)
   renormalize!(z)
   return z
end

doc"""
   tan(a::fmpq_rel_series)
> Return tan$(a)$. Requires a zero constant term.
"""
function tan(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in tan")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_tan_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   tanh(a::fmpq_rel_series)
> Return tanh$(a)$. Requires a zero constant term.
"""
function tanh(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in tanh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_tanh_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   sin(a::fmpq_rel_series)
> Return sin$(a)$. Requires a zero constant term.
"""
function sin(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in sin")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_truncate, :libflint), Void,
                (Ref{fmpq_rel_series}, Int),
               z, a.prec)
   ccall((:fmpq_poly_sin_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   sinh(a::fmpq_rel_series)
> Return sinh$(a)$. Requires a zero constant term.
"""
function sinh(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in sinh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_sinh_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   cos(a::fmpq_rel_series)
> Return cos$(a)$. Requires a zero constant term.
"""
function cos(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in cos")
   if pol_length(a) + valuation(a) == 0 || a.prec == 1
      return one(parent(a))
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_truncate, :libflint), Void,
                (Ref{fmpq_rel_series}, Int),
               z, a.prec)
   ccall((:fmpq_poly_cos_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   cosh(a::fmpq_rel_series)
> Return cosh$(a)$. Requires a zero constant term.
"""
function cosh(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in cosh")
   if pol_length(a) + valuation(a) == 0 || a.prec == 1
      return one(parent(a))
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_cosh_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   asin(a::fmpq_rel_series)
> Return asin$(a)$. Requires a zero constant term.
"""
function asin(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in asin")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_asin_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   asinh(a::fmpq_rel_series)
> Return asinh$(a)$. Requires a zero constant term.
"""
function asinh(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in asinh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_asinh_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   atan(a::fmpq_rel_series)
> Return atan$(a)$. Requires a zero constant term.
"""
function atan(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in atan")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_atan_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   atanh(a::fmpq_rel_series)
> Return atanh$(a)$. Requires a zero constant term.
"""
function atanh(a::fmpq_rel_series)
   (a.val == 0 && pol_length(a) != 0) && error("Constant term not zero in atanh")
   if iszero(a) || a.prec < 2
      return parent(a)()
   end
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_shift_left, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.val)
   ccall((:fmpq_poly_atanh_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, z, a.prec)
   renormalize!(z)
   return z
end

doc"""
   sqrt(a::fmpq_rel_series)
> Return the power series square root of $a$. Requires a constant term equal to
> one.
"""
function sqrt(a::fmpq_rel_series)
   (a.val != 0 || coeff(a, 0) != 1) && error("Constant term not one in sqrt")
   z = parent(a)()
   z.prec = a.prec
   z.val = 0
   ccall((:fmpq_poly_sqrt_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, a.prec)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::fmpq_rel_series)
   ccall((:fmpq_poly_zero, :libflint), Void,
                (Ref{fmpq_rel_series},), z)
   z.prec = parent(z).prec_max
   return z
end

function setcoeff!(z::fmpq_rel_series, n::Int, x::fmpq)
   ccall((:fmpq_poly_set_coeff_fmpq, :libflint), Void,
                (Ref{fmpq_rel_series}, Int, Ref{fmpq}),
               z, n, x)
   return z
end

function mul!(z::fmpq_rel_series, a::fmpq_rel_series, b::fmpq_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   aval = valuation(a)
   bval = valuation(b)
   prec = min(a.prec - aval, b.prec - bval)
   lena = min(lena, prec)
   lenb = min(lenb, prec)
   z.val = a.val + b.val
   z.prec = prec + c.val
   lenz = min(lena + lenb - 1, prec)
   if lena <= 0 || lenb <= 0
      lenz = 0
   end
   ccall((:fmpq_poly_mullow, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               z, a, b, lenz)
   return z
end

function addeq!(a::fmpq_rel_series, b::fmpq_rel_series)
   lena = pol_length(a)
   lenb = pol_length(b)
   prec = min(a.prec, b.prec)
   val = min(a.val, b.val)
   lena = min(lena, prec - a.val)
   lenb = min(lenb, prec - b.val)
   if a.val < b.val
      z = fmpq_rel_series()
      lenz = max(lena, lenb + b.val - a.val)
      ccall((:fmpq_poly_set_trunc, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, b, max(0, lenz - b.val + a.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            z, z, b.val - a.val)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               a, a, z, lenz)
   elseif b.val < a.val
      lenz = max(lena + a.val - b.val, lenb)
      ccall((:fmpq_poly_truncate, :libflint), Void,
            (Ref{fmpq_rel_series}, Int),
            a, max(0, lenz - a.val + b.val))
      ccall((:fmpq_poly_shift_left, :libflint), Void,
            (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
            a, a, a.val - b.val)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               a, a, b, lenz)
   else
      lenz = max(lena, lenb)
      ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               a, a, b, lenz)
   end
   a.prec = prec
   a.val = val
   renormalize!(a)
   return a
end

function add!(c::fmpq_rel_series, a::fmpq_rel_series, b::fmpq_rel_series)
   lena = length(a)
   lenb = length(b)

   prec = min(a.prec, b.prec)

   lena = min(lena, prec)
   lenb = min(lenb, prec)

   lenc = max(lena, lenb)
   c.prec = prec
   ccall((:fmpq_poly_add_series, :libflint), Void,
                (Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Ref{fmpq_rel_series}, Int),
               c, a, b, lenc)
   return c
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{fmpq_rel_series}, ::Type{T}) where {T <: Integer} = fmpq_rel_series

promote_rule(::Type{fmpq_rel_series}, ::Type{Rational{T}}) where T <: Union{Int, BigInt} = fmpq_rel_series

promote_rule(::Type{fmpq_rel_series}, ::Type{fmpz}) = fmpq_rel_series

promote_rule(::Type{fmpq_rel_series}, ::Type{fmpq}) = fmpq_rel_series

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FmpqRelSeriesRing)()
   z = fmpq_rel_series()
   z.prec = a.prec_max
   z.val = a.prec_max
   z.parent = a
   return z
end

function (a::FmpqRelSeriesRing)(b::Integer)
   if b == 0
      z = fmpq_rel_series()
      z.prec = a.prec_max
   else
      z = fmpq_rel_series([fmpq(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpqRelSeriesRing)(b::fmpz)
   if iszero(b)
      z = fmpq_rel_series()
      z.prec = a.prec_max
   else
      z = fmpq_rel_series([fmpq(b)], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

function (a::FmpqRelSeriesRing)(b::fmpq)
   if iszero(b)
      z = fmpq_rel_series()
      z.prec = a.prec_max
   else
      z = fmpq_rel_series([b], 1, a.prec_max, 0)
   end
   z.parent = a
   return z
end

(a::FmpqRelSeriesRing)(b::Rational{T}) where T <: Union{Int, BigInt} = a(fmpq(b))

function (a::FmpqRelSeriesRing)(b::fmpq_rel_series)
   parent(b) != a && error("Unable to coerce power series")
   return b
end

function (a::FmpqRelSeriesRing)(b::Array{fmpq, 1}, len::Int, prec::Int, val::Int)
   z = fmpq_rel_series(b, len, prec, val)
   z.parent = a
   return z
end

(a::FmpqRelSeriesRing)(b::Array{fmpz, 1}, len::Int, prec::Int, val::Int) =
    a(map(fmpq, b), len, prec, val)

(a::FmpqRelSeriesRing)(b::Array{T, 1}, len::Int, prec::Int, val::Int) where {T <: Integer} =
    a(map(fmpq, b), len, prec, val)

(a::FmpqRelSeriesRing)(b::Array{Rational{T}, 1}, len::Int, prec::Int, val::Int) where {T <: Integer} =
    a(map(fmpq, b), len, prec, val)

ext: jl


text:
using OrdinaryDiffEq
using DiffEqFlux, Flux
using Random
using Plots;

Random.seed!(123);

function sir_ode(u,p,t)
    (s,i,r) = u
    (β,γ) = p
    ds = -β*s*i
    di = β*s*i - γ*i
    dr = γ*i
    [ds,di,dr]
end;

p = [0.5,0.25]
u0 = [0.99, 0.01, 0.0]
tspan = (0.0, 40.0)
δt = 1;

solver = Rodas5();

sir_prob = ODEProblem(sir_ode, u0, tspan, p)
sir_sol = solve(sir_prob, solver, saveat = δt);

train_time = 30.0
train_data = Array(sir_sol(0:δt:train_time));

nhidden = 8
sir_node = FastChain(FastDense(3, nhidden, tanh),
                     FastDense(nhidden, nhidden, tanh),
                     FastDense(nhidden, nhidden, tanh),
                     FastDense(nhidden, 3));

p_ = Float64.(initial_params(sir_node));
function dudt_sir_node(u,p,t)
    s,i,r = u
    ds,di,dr = ann_node([s,i,r],p)
    [ds,di,dr]
end
prob_node = ODEProblem(dudt_sir_node, u0, tspan, p_);

tspan_train = (0,train_time)
prob_node = NeuralODE(sir_node,
                      tspan_train,
                      solver,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))
nump = length(prob_node.p)

function loss(p)
    sol = prob_node(u0,p)
    pred = Array(sol)
    sum(abs2, (train_data .- pred)), pred
end;

const losses = []
callback = function (p, l, pred)
    push!(losses, l)
    numloss = length(losses)
    if numloss % 50 == 0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end;

res_node = DiffEqFlux.sciml_train(loss,
                                   prob_node.p,
                                   cb = callback);

plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss", legend = false)

prob_node = NeuralODE(sir_node,
                      tspan_train,
                      solver,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()),
                      p = res_node)
sol_node = prob_node(u0);

scatter(sir_sol, label=["True Susceptible" "True Infected" "True Recovered"])
plot!(sol_node, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

tspan_test = (0.0, 40.0)
prob_node_test = NeuralODE(sir_node,
                      tspan_test,
                      solver,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()),
                      p = res_node)
sol_node_test = prob_node_test(u0);

p_node = scatter(sol_node_test, legend = :topright, label=["True Susceptible" "True Infected" "True Recovered"], title="Neural ODE Extrapolation: training until t=30")
plot!(p_node,sol_node_test, lw=5, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

newu0 = [0.95, 0.05, 0.0]
sir_prob_u0 = remake(sir_prob,u0=newu0)
sir_sol_u0 = solve(sir_prob_u0, solver, saveat =  δt)
node_sol_u0 = prob_node(newu0)
p_node = scatter(sir_sol_u0, legend = :topright, label=["True Susceptible" "True Infected" "True Recovered"], title="Neural ODE with different initial conditions")
plot!(p_node,node_sol_u0, lw=5, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

tspan_train2 = (0.0,20.0)
prob2 = ODEProblem(sir_ode, u0, tspan_train2, p)
sol2 = solve(prob2, solver, saveat = δt)
data2 = Array(sol2)
solver2 = ExplicitRK()
prob_node2 = NeuralODE(sir_node,
                      tspan_train2,
                      solver2,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()))
function loss2(p)
    sol = prob_node2(u0,p)
    pred = Array(sol)
    sum(abs2, (data2 .- sol)), pred
end
const losses2 = []
callback2 = function (p, l, pred)
    push!(losses2, l)
    numloss = length(losses2)
    if numloss % 50 == 0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end
res_node2 = DiffEqFlux.sciml_train(loss2,
                                   prob_node2.p,
                                   cb = callback2);

prob_node2_test = NeuralODE(sir_node,
                      tspan_test,
                      solver,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()),
                      p = res_node2)
sol_node2_test = prob_node2_test(u0);

p_node2 = scatter(sir_sol, legend = :topright, label=["True Susceptible" "True Infected" "True Recovered"], title="Neural ODE Extrapolation: training until t=20")
plot!(p_node2, sol_node2_test, lw=5, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

function loss3(p)
    sol = prob_node(u0,p)
    pred = Array(sol)
    sum(abs2, (train_data[2,:] .- sol[2,:])), pred
end
const losses3 = []
callback3 = function (p, l, pred)
    push!(losses3, l)
    numloss = length(losses3)
    if numloss % 50 ==0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end
res_node3 = DiffEqFlux.sciml_train(loss3,
                                   prob_node.p,
                                   cb = callback3);

prob_node3_test = NeuralODE(sir_node,
                      tspan_test,
                      solver,
                      saveat=δt,
                      sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()),
                      p = res_node3)
sol_node3_test = prob_node3_test(u0);

p_node3 = scatter(sir_sol, legend = :topright, label=["True Susceptible" "True Infected" "True Recovered"], title="Neural ODE Extrapolation: training on I(t) until t=30")
plot!(p_node3, sol_node3_test, lw=5, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])


text:
# Functions for calculating time response of a system

# XXX : `step` is a function in Base, with a different meaning than it has
# here. This shouldn't be an issue, but it might be.
"""`y, t, x = step(sys[, Tf])` or `y, t, x = step(sys[, t])`

Calculate the step response of system `sys`. If the final time `Tf` or time
vector `t` is not provided, one is calculated based on the system pole
locations.

`y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`"""
function Base.step(sys::StateSpace, t::AbstractVector; method=:cont)
    lt = length(t)
    ny, nu = size(sys)
    nx = sys.nx
    u = (x,t)->[one(eltype(t))]
    x0 = zeros(nx)
    if nu == 1
        y, tout, x, _ = lsim(sys, u, t, x0=x0, method=method)
    else
        x = Array{Float64}(undef, lt, nx, nu)
        y = Array{Float64}(undef, lt, ny, nu)
        for i=1:nu
            y[:,:,i], tout, x[:,:,i],_ = lsim(sys[:,i], u, t, x0=x0, method=method)
        end
    end
    return y, t, x
end

Base.step(sys::LTISystem, Tf::Real; kwargs...) = step(sys, _default_time_vector(sys, Tf); kwargs...)
Base.step(sys::LTISystem; kwargs...) = step(sys, _default_time_vector(sys); kwargs...)
Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys), t::AbstractVector; kwargs...)

"""`y, t, x = impulse(sys[, Tf])` or `y, t, x = impulse(sys[, t])`

Calculate the impulse response of system `sys`. If the final time `Tf` or time
vector `t` is not provided, one is calculated based on the system pole
locations.

`y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`"""
function impulse(sys::StateSpace, t::AbstractVector; method=:cont)
    T = promote_type(eltype(sys.A), Float64)
    lt = length(t)
    ny, nu = size(sys)
    nx = sys.nx
    if iscontinuous(sys) #&& method == :cont
        u = (x,i) -> [zero(T)]
        # impulse response equivalent to unforced response of
        # ss(A, 0, C, 0) with x0 = B.
        imp_sys = ss(sys.A, zeros(T, nx, 1), sys.C, zeros(T, ny, 1))
        x0s = sys.B
    else
        u = (x,i) -> i == t[1] ? [one(T)]/sys.Ts : [zero(T)]
        imp_sys = sys
        x0s = zeros(T, nx, nu)
    end
    if nu == 1 # Why two cases # QUESTION: Not type stable?
        y, t, x,_ = lsim(sys, u, t, x0=x0s[:], method=method)
    else
        x = Array{T}(undef, lt, nx, nu)
        y = Array{T}(undef, lt, ny, nu)
        for i=1:nu
            y[:,:,i], t, x[:,:,i],_ = lsim(sys[:,i], u, t, x0=x0s[:,i], method=method)
        end
    end
    return y, t, x
end

impulse(sys::LTISystem, Tf::Real; kwags...) = impulse(sys, _default_time_vector(sys, Tf); kwags...)
impulse(sys::LTISystem; kwags...) = impulse(sys, _default_time_vector(sys); kwags...)
impulse(sys::TransferFunction, t::AbstractVector; kwags...) = impulse(ss(sys), t; kwags...)

"""`y, t, x = lsim(sys, u, t; x0, method])`

`y, t, x, uout = lsim(sys, u::Function, t; x0, method)`

Calculate the time response of system `sys` to input `u`. If `x0` is ommitted,
a zero vector is used.

`y`, `x`, `uout` has time in the first dimension. Initial state `x0` defaults to zero.

Continuous time systems are simulated using an ODE solver if `u` is a function. If `u` is an array, the system is discretized before simulation. For a lower level inteface, see `?Simulator` and `?solve`

`u` can be a function or a matrix/vector of precalculated control signals.
If `u` is a function, then `u(x,i)` (`u(x,t)`) is called to calculate the control signal every iteration (time instance used by solver). This can be used to provide a control law such as state feedback `u(x,t) = -L*x` calculated by `lqr`.
To simulate a unit step, use `(x,i)-> 1`, for a ramp, use `(x,i)-> i*h`, for a step at `t=5`, use (x,i)-> (i*h >= 5) etc.

Usage example:
```julia
using LinearAlgebra # For identity matrix I
A = [0 1; 0 0]
B = [0;1]
C = [1 0]
sys = ss(A,B,C,0)
Q = I
R = I
L = lqr(sys,Q,R)

u(x,t) = -L*x # Form control law,
t=0:0.1:5
x0 = [1,0]
y, t, x, uout = lsim(sys,u,t,x0)
plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")
```
"""
function lsim(sys::StateSpace, u::AbstractVecOrMat, t::AbstractVector;
        x0::VecOrMat=zeros(sys.nx), method::Symbol=_issmooth(u) ? :foh : :zoh)
    ny, nu = size(sys)
    nx = sys.nx

    if length(x0) != nx
        error("size(x0) must match the number of states of sys")
    end
    if !(size(u) in [(length(t), nu) (length(t),)])
        error("u must be of size (length(t), nu)")
    end

    dt = Float64(t[2] - t[1])
    if !iscontinuous(sys) || method == :zoh
        if iscontinuous(sys) # Looks strange to check iscontinuous again
            dsys = c2d(sys, dt, :zoh)[1]
        else
            if sys.Ts != dt
                error("Time vector must match sample time for discrete system")
            end
            dsys = sys
        end
    else
        dsys, x0map = c2d(sys, dt, :foh)
        x0 = x0map*[x0; transpose(u[1:1,:])]
    end
    x = ltitr(dsys.A, dsys.B, Float64.(u), Float64.(x0))
    y = transpose(sys.C*transpose(x) + sys.D*transpose(u))
    return y, t, x
end

@deprecate lsim(sys, u, t, x0) lsim(sys, u, t; x0=x0)
@deprecate lsim(sys, u, t, x0, method) lsim(sys, u, t; x0=x0, method=method)

function lsim(sys::StateSpace, u::Function, t::AbstractVector;
        x0::VecOrMat=zeros(sys.nx), method::Symbol=:cont)
    ny, nu = size(sys)
    nx = sys.nx
    if length(x0) != nx
        error("size(x0) must match the number of states of sys")
    elseif size(u(x0,1)) != (nu,) && size(u(x0,1)) != (nu,1)
        error("return value of u must be of size nu")
    end
    T = promote_type(Float64, eltype(x0))

    dt = T(t[2] - t[1])
    if !iscontinuous(sys) || method == :zoh
        if iscontinuous(sys)
            dsys = c2d(sys, dt, :zoh)[1]
        else
            if sys.Ts != dt
                error("Time vector must match sample time for discrete system")
            end
            dsys = sys
        end
        x,uout = ltitr(dsys.A, dsys.B, u, t, T.(x0))
    else
        s = Simulator(sys, u)
        sol = solve(s, T.(x0), (t[1],t[end]), Tsit5())
        x = sol(t)'
        uout = Array{eltype(x)}(undef, length(t), ninputs(sys))
        for (i,ti) in enumerate(t)
            uout[i,:] = u(x[i,:],ti)'
        end
    end
    y = transpose(sys.C*transpose(x) + sys.D*transpose(uout))
    return y, t, x, uout
end

lsim(sys::TransferFunction, u, t, args...; kwargs...) = lsim(ss(sys), u, t, args...; kwargs...)

"""`ltitr(A, B, u[,x0])`

`ltitr(A, B, u::Function, iters[,x0])`

Simulate the discrete time system `x[k + 1] = A x[k] + B u[k]`, returning `x`.
If `x0` is not provided, a zero-vector is used.

If `u` is a function, then `u(x,i)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u=-Lx` calculated by `lqr`. In this case, an integrer `iters` must be provided that indicates the number of iterations.
"""
function ltitr(A::Matrix{T}, B::Matrix{T}, u::AbstractVecOrMat{T},
        x0::VecOrMat{T}=zeros(T, size(A, 1))) where T
    n = size(u, 1)
    x = Array{T}(undef, size(A, 1), n)
    for i=1:n
        x[:,i] = x0
        x0 = A * x0 + B * u[i,:]
    end
    return transpose(x)
end

function ltitr(A::Matrix{T}, B::Matrix{T}, u::Function, t,
    x0::VecOrMat{T}=zeros(T, size(A, 1))) where T
    iters = length(t)
    x = Array{T}(undef, size(A, 1), iters)
    uout = Array{T}(undef, size(B, 2), iters)

    for i=1:iters
        x[:,i] = x0
        uout[:,i] = u(x0,t[i])
        x0 = A * x0 + B * uout[:,i]
    end
    return transpose(x), transpose(uout)
end

# HELPERS:

# TODO: This is a poor heuristic to estimate a "good" time vector to use for
# simulation, in cases when one isn't provided.
function _default_time_vector(sys::LTISystem, Tf::Real=-1)
    Ts = _default_Ts(sys)
    if Tf == -1
        Tf = 100*Ts
    end
    return 0:Ts:Tf
end

function _default_Ts(sys::LTISystem)
    if !iscontinuous(sys)
        Ts = sys.Ts
    elseif !isstable(sys)
        Ts = 0.05
    else
        ps = pole(sys)
        r = minimum([abs.(real.(ps));0])
        if r == 0.0
            r = 1.0
        end
        Ts = 0.07/r
    end
    return Ts
end

#TODO a reasonable check
_issmooth(u::Function) = false

# Determine if a signal is "smooth"
function _issmooth(u, thresh::AbstractFloat=0.75)
    u = [zeros(1, size(u, 2)); u]       # Start from 0 signal always
    dist = maximum(u) - minimum(u)
    du = abs.(diff(u, dims=1))
    return !isempty(du) && all(maximum(du) <= thresh*dist)
end


text:
using FluxRC, KitBase, Plots, OrdinaryDiffEq
using ProgressMeter: @showprogress

cd(@__DIR__)
ps = UnstructPSpace("../assets/square.msh")

N = deg = 2
Np = (N + 1) * (N + 2) ÷ 2
ncell = size(ps.cellid, 1)
nface = size(ps.faceType, 1)

J = rs_jacobi(ps.cellid, ps.points)

spg = global_sp(ps.points, ps.cellid, N)
fpg = global_fp(ps.points, ps.cellid, N)

pl, wl = tri_quadrature(N)

V = vandermonde_matrix(N, pl[:, 1], pl[:, 2])
Vr, Vs = ∂vandermonde_matrix(N, pl[:, 1], pl[:, 2])
∂l = ∂lagrange(V, Vr, Vs)

ϕ = correction_field(N, V)

pf, wf = triface_quadrature(N)
ψf = zeros(3, N + 1, Np)
for i = 1:3
    ψf[i, :, :] .= vandermonde_matrix(N, pf[i, :, 1], pf[i, :, 2])
end

lf = zeros(3, N + 1, Np)
for i = 1:3, j = 1:N+1
    lf[i, j, :] .= V' \ ψf[i, j, :]
end

a = -1.0
u = zeros(size(ps.cellid, 1), Np)
for i in axes(u, 1), j in axes(u, 2)
    u[i, j] = max(exp(-100 * ((spg[i, j, 1] - 0.5)^2 + (spg[i, j, 2] - 0.5)^2)), 1e-4)
end

function dudt!(du, u, p, t)
    du .= 0.0

    a, deg, nface = p
    ncell = size(u, 1)
    nsp = size(u, 2)

    f = zeros(ncell, nsp, 2)
    for i in axes(f, 1)
        for j in axes(f, 2)
            fg = a * u[i, j]
            gg = a * u[i, j]
            f[i, j, :] .= inv(J[i]) * [fg, gg]
        end
    end

    u_face = zeros(ncell, 3, deg + 1)
    f_face = zeros(ncell, 3, deg + 1, 2)
    for i = 1:ncell, j = 1:3, k = 1:deg+1
        u_face[i, j, k] = sum(u[i, :] .* lf[j, k, :])
        f_face[i, j, k, 1] = sum(f[i, :, 1] .* lf[j, k, :])
        f_face[i, j, k, 2] = sum(f[i, :, 2] .* lf[j, k, :])
    end

    n = [[0.0, -1.0], [1 / √2, 1 / √2], [-1.0, 0.0]]

    fn_face = zeros(ncell, 3, deg + 1)
    for i = 1:ncell, j = 1:3, k = 1:deg+1
        fn_face[i, j, k] = sum(f_face[i, j, k, :] .* n[j])
    end

    f_interaction = zeros(ncell, 3, deg + 1, 2)
    au = zeros(2)
    for i = 1:ncell, j = 1:3, k = 1:deg+1
        fL = J[i] * f_face[i, j, k, :]

        ni, nj, nk = neighbor_fpidx([i, j, k], ps, fpg)

        fR = zeros(2)
        if ni > 0
            fR .= J[ni] * f_face[ni, nj, nk, :]

            @. au = (fL - fR) / (u_face[i, j, k] - u_face[ni, nj, nk] + 1e-6)
            @. f_interaction[i, j, k, :] =
                0.5 * (fL + fR) - 0.5 * abs(au) * (u_face[i, j, k] - u_face[ni, nj, nk])
        else
            @. f_interaction[i, j, k, :] = 0.0
        end

        f_interaction[i, j, k, :] .= inv(J[i]) * f_interaction[i, j, k, :]
    end

    fn_interaction = zeros(ncell, 3, deg + 1)
    for i = 1:ncell
        for j = 1:3, k = 1:deg+1
            fn_interaction[i, j, k] = sum(f_interaction[i, j, k, :] .* n[j])
        end
    end

    rhs1 = zeros(ncell, nsp)
    for i in axes(rhs1, 1), j in axes(rhs1, 2)
        rhs1[i, j] = -sum(f[i, :, 1] .* ∂l[j, :, 1]) - sum(f[i, :, 2] .* ∂l[j, :, 2])
    end

    rhs2 = zero(rhs1)
    for i = 1:ncell
        xr, yr = ps.points[ps.cellid[i, 2], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
        xs, ys = ps.points[ps.cellid[i, 3], 1:2] - ps.points[ps.cellid[i, 1], 1:2]
        _J = xr * ys - xs * yr

        if ps.cellType[i] != 1
            for j = 1:nsp
                #rhs2[i, j] = - sum((fn_interaction[i, :, :] .- fn_face[i, :, :]) .* ϕ[:, :, j]) / _J
                rhs2[i, j] =
                    -sum((fn_interaction[i, :, :] .- fn_face[i, :, :]) .* ϕ[:, :, j])
                #rhs2[i, j] = - sum((fn_interaction[i, :, :] .- fn_face[i, :, :]))
            end
        end
    end

    for i = 1:ncell
        if ps.cellType[i] == 0
            du[i, :] .= rhs1[i, :] .+ rhs2[i, :]
        end
    end

    return nothing
end

tspan = (0.0, 1.0)
p = (a, N, nface)
prob = ODEProblem(dudt!, u, tspan, p)
#sol = solve(prob, Midpoint(), reltol=1e-8, abstol=1e-8, progress=true)

dt = 0.01
itg = init(prob, Euler(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:10
    step!(itg)
end

write_vtk(ps.points, ps.cellid, itg.u[:, 2])

du = zero(u)
dudt!(du, u, p, 1.0)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct FileStatusesResult <: SwaggerModel
    FileStatuses::Any # spec type: Union{ Nothing, FileStatuses } # spec name: FileStatuses

    function FileStatusesResult(;FileStatuses=nothing)
        o = new()
        validate_property(FileStatusesResult, Symbol("FileStatuses"), FileStatuses)
        setfield!(o, Symbol("FileStatuses"), FileStatuses)
        o
    end
end # type FileStatusesResult

const _property_map_FileStatusesResult = Dict{Symbol,Symbol}(Symbol("FileStatuses")=>Symbol("FileStatuses"))
const _property_types_FileStatusesResult = Dict{Symbol,String}(Symbol("FileStatuses")=>"FileStatuses")
Base.propertynames(::Type{ FileStatusesResult }) = collect(keys(_property_map_FileStatusesResult))
Swagger.property_type(::Type{ FileStatusesResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_FileStatusesResult[name]))}
Swagger.field_name(::Type{ FileStatusesResult }, property_name::Symbol) =  _property_map_FileStatusesResult[property_name]

function check_required(o::FileStatusesResult)
    true
end

function validate_property(::Type{ FileStatusesResult }, name::Symbol, val)
end


text:
## arithmetic.jl operators for NamedArray

## (c) 2013--2014 David A. van Leeuwen

## This code is licensed under the MIT License
## See the file LICENSE.md in this distribution

import Base: +, -, *, /, \

-(n::NamedArray) = NamedArray(-n.array, n.dicts, n.dimnames)

## disambiguation magic

# disambiguation (Argh...)
for op in (:+, :-)
    @eval ($op)(x::AbstractRange{T1}, y::NamedVector{T2}) where {T1<:Number,T2<:Number} = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    @eval ($op)(x::NamedVector{T1}, y::AbstractRange{T2}) where {T1<:Number,T2<:Number} = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
end

for op in (:+, :-)
    ## named %op% named
    @eval begin
        function ($op)(x::NamedArray{T1}, y::NamedArray{T2}) where {T1<:Number, T2<:Number}
            if names(x) == names(y) && x.dimnames == y.dimnames
                NamedArray(($op)(x.array, y.array), x.dicts, x.dimnames)
            else
                @warn("Dropping mismatching names")
                ($op)(x.array, y.array)
            end
        end
        ($op)(x::NamedArray{T1,N}, y::AbstractArray{T2,N}) where {T1<:Number,T2<:Number,N} = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
        ($op)(x::AbstractArray{T1,N}, y::NamedArray{T2,N}) where {T1<:Number,T2<:Number,N} = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    end
end

## scalar arithmetic
## disambiguate
for op in (:+, :-)
    @eval begin
        ($op)(x::NamedArray{Bool}, y::Bool) = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
        ($op)(x::Bool, y::NamedArray{Bool}) = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    end
end

## NamedArray, Number

for op in (:+, :-, :*)
    @eval begin
        ($op)(x::NamedArray{T1}, y::T2) where {T1<:Number,T2<:Number} = NamedArray(($op)(x.array, y), x.dicts, x.dimnames)
        ($op)(x::T1, y::NamedArray{T2}) where {T1<:Number,T2<:Number} = NamedArray(($op)(x, y.array), y.dicts, y.dimnames)
    end
end
/(x::NamedArray{T1}, y::T2) where {T1<:Number,T2<:Number} = NamedArray(x.array / y, x.dicts, x.dimnames)
\(x::T1, y::NamedArray{T2}) where {T1<:Number,T2<:Number} = NamedArray(x \ y.array, y.dicts, y.dimnames)

# import LinearAlgebra: A_mul_B!, A_mul_Bc!, A_mul_Bc, A_mul_Bt!, A_mul_Bt, Ac_mul_B, Ac_mul_B!, Ac_mul_Bc, Ac_mul_Bc!, At_mul_B, At_mul_B!, At_mul_Bt, At_mul_Bt!

## Assume dimensions/names are correct
for op in (:A_mul_B!, :A_mul_Bc!, :A_mul_Bt!, :Ac_mul_B!, :Ac_mul_Bc!, :At_mul_B!, :At_mul_Bt!)
    @eval ($op)(C::NamedMatrix, A::AbstractMatrix, B::AbstractMatrix) = ($op)(C.array, A, B)
end

for op in (:A_mul_Bc, :A_mul_Bt)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[1], B.dicts[1]), (A.dimnames[1], B.dimnames[1]))
    for T in [Union{LinearAlgebra.QRCompactWYQ, LinearAlgebra.QRPackedQ}, StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[1], defaultnamesdict(size(B,1))), (A.dimnames[1], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,1)), B.dicts[1]), (:A, B.dimnames[1]))
    end
end
for op in (:Ac_mul_B, :At_mul_B)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[2], B.dicts[2]), (A.dimnames[2], B.dimnames[2]))
    for T in [StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[2], defaultnamesdict(size(B,2))), (A.dimnames[2], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,2)), B.dicts[2]), (:A, B.dimnames[2]))
    end
end
for op in (:Ac_mul_Bc, :At_mul_Bt)
    @eval ($op)(A::NamedMatrix, B::NamedMatrix) = NamedArray(($op)(A.array, B.array), (A.dicts[2], B.dicts[1]), (A.dimnames[2], B.dimnames[1]))
    for T in [StridedMatrix, AbstractMatrix] ## v0.4 ambiguity-hell
        @eval ($op)(A::NamedMatrix, B::$T) = NamedArray(($op)(A.array, B), (A.dicts[2], defaultnamesdict(size(B,1))), (A.dimnames[2], :B))
        @eval ($op)(A::$T, B::NamedMatrix) = NamedArray(($op)(A, B.array), (defaultnamesdict(size(A,2)), B.dicts[1]), (:A, B.dimnames[1]))
    end
end

import LinearAlgebra: Givens, LinearAlgebra.BlasFloat, lu!, lu, LU, ipiv2perm, qr!, qr, eigen!, eigen, eigvals!,
    eigvals, hessenberg, hessenberg!, schur!, schur, svd!, svd, svdvals!, svdvals, svd, diag,
    cond, kron, lyap, sylvester, isposdef

## matmul
## ambiguity, this can somtimes be a pain to resolve...
*(x::SparseMatrixCSC{Tx,TiA},y::NamedMatrix{Ty}) where {Tx,TiA,Ty} = x*y.array
*(x::SparseMatrixCSC{Tx,S},y::NamedVector{Ty}) where {Tx,S,Ty} = x*y.array
for t in (:Tridiagonal, :(LinearAlgebra.AbstractTriangular), :Givens, :Bidiagonal)
    @eval *(x::$t, y::NamedMatrix) = NamedArray(x*y.array, ([string(i) for i in 1:size(x,1)],names(y, 2)), y.dimnames)
    @eval *(x::$t, y::NamedVector) = x*y.array
end

## There is no such thing as a A_mul_B
## Named * Named
*(A::NamedMatrix, B::NamedMatrix) = NamedArray(A.array * B.array, (A.dicts[1], B.dicts[2]), (A.dimnames[1], B.dimnames[2]))
function *(A::NamedMatrix, B::NamedVector)
    result = A.array * B.array
    if isa(result, AbstractArray)
        return NamedArray(result, (A.dicts[1],), (B.dimnames[1],))
    else
        return result
    end
end

if  @isdefined RowVector
    *(A::NamedRowVector, B::NamedVector) = A.array * B.array
end
## Named * Abstract
*(A::NamedMatrix, B::AbstractMatrix) = NamedArray(A.array * B, (A.dicts[1], defaultnamesdict(size(B,2))), A.dimnames)
*(A::AbstractMatrix, B::NamedMatrix) = NamedArray(A * B.array, (defaultnamesdict(size(A,1)), B.dicts[2]), B.dimnames)
*(A::NamedMatrix, B::AbstractVector) = NamedArray(A.array * B, (A.dicts[1],), (A.dimnames[1],))
*(A::AbstractMatrix, B::NamedVector) = A * B.array
if @isdefined RowVector
    *(A::NamedRowVector, B::AbstractVector) = A.array * B
end
## \ --- or should we overload A_div_B?
## Named \ Named
\(x::NamedVector, y::NamedVector) = x.array \ y.array
\(x::NamedMatrix, y::NamedVector) = NamedArray(x.array\y.array, (names(x, 2),), (x.dimnames[2],))
\(x::NamedVector, y::NamedMatrix) = NamedArray(x.array\y.array, (["1"],names(y, 2)), (:A, y.dimnames[2]))
\(x::NamedMatrix, y::NamedMatrix) = NamedArray(x.array\y.array, (names(x, 2),names(y, 2)), (x.dimnames[2], y.dimnames[2]))

## Named \ Abstract
\(x::NamedVector, y::AbstractVecOrMat) = x.array \ y
\(x::NamedMatrix, y::AbstractVector) = NamedArray(x.array \ y, (x.dicts[2],), (x.dimnames[2],))
\(x::NamedMatrix, y::AbstractMatrix) = NamedArray(x.array \ y, (names(x, 2),[string(i) for i in 1:size(y,2)]), (x.dimnames[2],:B))
## Abstract \ Named
## ambiguity
\(x::Diagonal{Tx}, y::NamedVector{Ty}) where {Tx<:Number,Ty<:Number} = x \ y.array
\(x::Union{Bidiagonal{Tx},LinearAlgebra.AbstractTriangular{Tx}}, y::NamedVector{Ty}) where {Tx<:Number,Ty<:Number} = x \ y.array
\(x::Union{Bidiagonal{Tx},LinearAlgebra.AbstractTriangular{Tx}}, y::NamedMatrix{Ty}) where {Tx<:Number,Ty<:Number} = NamedArray(x \ y.array, (defaultnamesdict(size(x,1)), y.dicts[2]), (:A, y.dimnames[2]))

\(x::Bidiagonal,y::NamedVector) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y, 2)), (:A, y.dimnames[2]))
\(x::Bidiagonal,y::NamedMatrix) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y, 2)), (:A, y.dimnames[2]))

## AbstractVectorOrMat gives us more ambiguities than separate entries...
\(x::AbstractVector, y::NamedVector) = x \ y.array
\(x::AbstractMatrix, y::NamedVector) = x \ y.array
\(x::AbstractVector, y::NamedMatrix) = NamedArray(x \ y.array, (["1"],names(y, 2)), (:A, y.dimnames[2]))
\(x::AbstractMatrix, y::NamedMatrix) = NamedArray(x \ y.array, ([string(i) for i in 1:size(x,2)], names(y, 2)), (:A, y.dimnames[2]))

## keeping names for some matrix routines
for f in (:inv, :sqrt, :pinv, :exp) # :chol
    @eval (LinearAlgebra.$f)(n::NamedArray) = NamedArray(($f)(n.array), n.dicts, n.dimnames)
end
# TODO :  `chol(A::AbstractMatrix)` is deprecated, use `(cholesky(A)).U` instead.

## tril, triu
LinearAlgebra.tril!(n::NamedMatrix, k::Integer) = (tril!(n.array, k); n)
LinearAlgebra.triu!(n::NamedMatrix, k::Integer) = (triu!(n.array, k); n)

## LU factorization
function LinearAlgebra.lu!(n::NamedArray{T, N, AT, DT}, pivot = Val(true); kargs...) where {T, N, AT, DT}
    luf = lu!(n.array, pivot; kargs...)
    LU(n, luf.ipiv, luf.info)
end

# From /stdlib/LinearAlgebra/src/lu.jl
function Base.getproperty(A::LU{T,NamedArray{T, N, AT, DT}}, d::Symbol) where {T, N, AT, DT}
    m, n = size(A)
    if d == :L
        L = tril!(getfield(A, :factors)[1:m, 1:min(m,n)])
        for i = 1:min(m,n); L[i,i] = one(T); end
        setnames!(L, defaultnames(L,2), 2)
        setdimnames!(L, :LU, 2)
        return L
    elseif d == :U
        U = triu!(getfield(A, :factors)[1:min(m,n), 1:n])
        setnames!(U, defaultnames(U,1), 1)
        setdimnames!(U, :LU, 1)
        return U
    elseif d == :p
        return ipiv2perm(getfield(A, :ipiv), m)
    elseif d == :P
        return Matrix{T}(I, m, m)[:,invperm(A.p)]
    else
        getfield(F, d)
    end
end

# TODO: function cholfact!(n::NamedArray{T}, uplo::Symbol=:U) where T<:LinearAlgebra.BlasFloat
# TODO:     ishermitian(n) || LinearAlgebra.non_hermitian_error("cholfact!")
# TODO:     return cholfact!(Hermitian(n, uplo))
# TODO: end

# TODO: cholfact(n::NamedArray{T}, uplo::Symbol=:U) where {T<:LinearAlgebra.BlasFloat} = cholfact!(copy(n), uplo)

## ldlt skipped

## from factorization
function qr!(n::NamedMatrix, pivot)
    qr = qr!(n.array, pivot)
    LinearAlgebra.QRCompactWY(NamedArray(qr.factors, n.dicts, n.dimnames), qr.T)
end
qr!(n::NamedArray) = qr!(n, Val(false))

LAPACK.gemqrt!(side::Char, trans::Char, V::NamedArray{BF}, T::StridedMatrix{BF}, C::StridedVecOrMat{BF}) where {BF<:LinearAlgebra.BlasFloat} = LAPACK.gemqrt!(side, trans, V.array, T, C)

qr(n::NamedMatrix{T}, pivot = Val(false)) where {T<:LinearAlgebra.BlasFloat} = qr!(copy(n), pivot)

eigen!(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigen!(n.array, permute=permute, scale=scale)
eigen(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigen!(copy(n.array), permute=permute, scale=scale)

eigvals!(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigvals!(n.array, permute=permute, scale=scale)
eigvals(n::NamedMatrix; permute::Bool=true, scale::Bool=true) = eigvals!(copy(n.array), permute=permute, scale=scale)

hessenberg!(n::NamedMatrix) = hessenberg!(n.array)
hessenberg(n::NamedMatrix) = hessenberg(copy(n.array))

schur!(n::NamedMatrix) = schur!(n.array)
schur(n::NamedMatrix) = schur!(copy(n.array))
schur(A::NamedMatrix, B::AbstractMatrix) = schur(A.array, B)
schur!(A::NamedMatrix, B::AbstractMatrix) = schur!(A.array, B)

svd!(n::NamedMatrix; full::Bool=false) = svd!(n.array; full=full)
svd(A::NamedMatrix{T}; full=false) where {T<:LinearAlgebra.BlasFloat} = svd!(copy(A), full=full)

svdvals!(n::NamedArray) = svdvals!(n.array)
svdvals(n::NamedArray) = svdvals(copy(n.array))

# TODO : diag(n::NamedMatrix) = NamedArray(diag(n.array), n.dicts[1:1], n.dimnames[1:1])

# TODO : diagm(n::NamedVector) = NamedArray(diagm(n.array), n.dicts[[1,1]], n.dimnames[[1,1]])

# rank, vecnorm, norm, condskeel, trace, det, logdet OK
cond(n::NamedArray) = cond(n.array)

# null(n::NamedArray) = null(n.array)

function kron(a::NamedArray, b::NamedArray)
    n = Array{typeof(AbstractString[])}(undef, 2)
    dn = AbstractString[]
    for dim in 1:2
        n[dim] = AbstractString[]
        for i in names(a, dim)
            for j in names(b, dim)
                push!(n[dim], string(i, "×", j))
            end
        end
        push!(dn, string(dimnames(a,dim), "×", dimnames(b,dim)))
    end
    NamedArray(kron(a.array, b.array), tuple(n...), tuple(dn...))
end

# linreg(x::NamedVector, y::AbstractVector) = linreg(x.array, y)

lyap(A::NamedMatrix, C::AbstractMatrix) = NamedArray(lyap(A.array,C), A.dicts, A.dimnames)

sylvester(A::NamedMatrix, B::AbstractMatrix, C::AbstractMatrix) = NamedArray(sylvester(A.array, B, C), A.dicts, A.dimnames)

## issym, istriu, istril OK
isposdef(n::NamedArray) = isposdef(n.array)

## eigs OK


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct UpgradeOperationHistoryStatus <: SwaggerModel
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    startTime::Any # spec type: Union{ Nothing, DateTime } # spec name: startTime
    endTime::Any # spec type: Union{ Nothing, DateTime } # spec name: endTime

    function UpgradeOperationHistoryStatus(;code=nothing, startTime=nothing, endTime=nothing)
        o = new()
        validate_property(UpgradeOperationHistoryStatus, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(UpgradeOperationHistoryStatus, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        validate_property(UpgradeOperationHistoryStatus, Symbol("endTime"), endTime)
        setfield!(o, Symbol("endTime"), endTime)
        o
    end
end # type UpgradeOperationHistoryStatus

const _property_map_UpgradeOperationHistoryStatus = Dict{Symbol,Symbol}(Symbol("code")=>Symbol("code"), Symbol("startTime")=>Symbol("startTime"), Symbol("endTime")=>Symbol("endTime"))
const _property_types_UpgradeOperationHistoryStatus = Dict{Symbol,String}(Symbol("code")=>"String", Symbol("startTime")=>"DateTime", Symbol("endTime")=>"DateTime")
Base.propertynames(::Type{ UpgradeOperationHistoryStatus }) = collect(keys(_property_map_UpgradeOperationHistoryStatus))
Swagger.property_type(::Type{ UpgradeOperationHistoryStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpgradeOperationHistoryStatus[name]))}
Swagger.field_name(::Type{ UpgradeOperationHistoryStatus }, property_name::Symbol) =  _property_map_UpgradeOperationHistoryStatus[property_name]

const _allowed_UpgradeOperationHistoryStatus_code = ["RollingForward", "Cancelled", "Completed", "Faulted"]

function check_required(o::UpgradeOperationHistoryStatus)
    true
end

function validate_property(::Type{ UpgradeOperationHistoryStatus }, name::Symbol, val)
    if name === Symbol("code")
        Swagger.validate_param(name, "UpgradeOperationHistoryStatus", :enum, val, _allowed_UpgradeOperationHistoryStatus_code)
    end
end


text:
function update_weights!(
    new_weights::Vector,
    weights::Vector,
    training_proportion::Real,
)
    n = round(Int, training_proportion * sum(weights))
    π = weights ./ sum(weights)
    d = Distributions.Multinomial(n, π)
    rand!(d, new_weights)
    return
end


text:
# Functions for generating gridworld PDDL problems
using PDDL

"Converts ASCII gridworlds to PDDL problem."
function ascii_to_pddl(str::String, name="doors-keys-gems-problem")
    rows = split(str, "\n", keepempty=false)
    width, height = maximum(length.(strip.(rows))), length(rows)
    doors, keys, gems = Const[], Const[], Const[]
    walls = parse_pddl("(= walls (new-bit-matrix false $width $height))")
    init = Term[walls]
    start, goal = Term[], pddl"(true)"
    for (y, row) in enumerate(rows)
        for (x, char) in enumerate(strip(row))
            if char == '.' # Unoccupied
                continue
            elseif char == 'W' # Wall
                wall = parse_pddl("(= walls (set-index walls true $y $x))")
                push!(init, wall)
            elseif char == 'D' # Door
                d = Const(Symbol("door$(length(doors)+1)"))
                push!(doors, d)
                append!(init, parse_pddl("(= (xloc $d) $x)", "(= (yloc $d) $y)"))
                push!(init, parse_pddl("(locked $d)"))
            elseif char == 'k' # Key
                k = Const(Symbol("key$(length(keys)+1)"))
                push!(keys, k)
                append!(init, parse_pddl("(= (xloc $k) $x)", "(= (yloc $k) $y)"))
            elseif char == 'g' || char == 'G' # Gem
                g = Const(Symbol("gem$(length(gems)+1)"))
                push!(gems, g)
                append!(init, parse_pddl("(= (xloc $g) $x)", "(= (yloc $g) $y)"))
                if char == 'G' goal = parse_pddl("(has $g)") end
            elseif char == 's' # Start position
                start = parse_pddl("(= xpos $x)", "(= ypos $y)")
            end
        end
    end
    append!(init, start)
    objtypes = merge(Dict(d => :door for d in doors),
                     Dict(k => :key for k in keys),
                     Dict(g => :gem for g in gems))
    problem = GenericProblem(Symbol(name), Symbol("doors-keys-gems"),
                             [doors; keys; gems], objtypes, init, goal, nothing)
    return problem
end

function load_ascii_problem(path::AbstractString)
    str = open(f->read(f, String), path)
    return ascii_to_pddl(str)
end

function convert_ascii_problem(path::String)
    str = open(f->read(f, String), path)
    str = ascii_to_pddl(str)
    new_path = join(split(path, ".")[1:end-1]) * ".pddl"
    write(new_path, write_problem(str))
    return new_path
end