text:
using Test
using BSeries

using BSeries.Latexify: latexify

using StaticArrays: @SArray, @SMatrix, @SVector

using SymEngine: SymEngine
using SymPy: SymPy
using Symbolics: Symbolics

@testset "BSeries" begin

@testset "lazy representation of exact ODE solution" begin
  exact = ExactSolution{Rational{Int}}()
  terms = collect(Iterators.take(exact, 4))
  @test terms == [1//1, 1//2, 1//6, 1//3]
  @test exact == ExactSolution(exact)
end

@testset "latexify" begin
  # explicit midpoint method
  A = @SArray [0 0; 1//2 0]; b = @SArray [0, 1//1]; c = @SArray [0, 1//2];

  series_integrator = bseries(A, b, c, 2)
  @test_nowarn latexify(series_integrator)
  @test_nowarn latexify(series_integrator, cdot=false)
  @test_nowarn latexify(series_integrator, dt=SymEngine.symbols("h"))

  @testset "SymEngine.jl" begin
    α = SymEngine.symbols("α")
    A = [0 0; 1/(2*α) 0]; b = [1-α, α]; c = [0, 1/(2*α)]
    series_integrator = bseries(A, b, c, 3)
    @test_nowarn latexify(series_integrator)
  end

  @testset "SymPy.jl" begin
    α = SymPy.symbols("α", real=true)
    A = [0 0; 1/(2*α) 0]; b = [1-α, α]; c = [0, 1/(2*α)]
    series_integrator = bseries(A, b, c, 3)
    @test_nowarn latexify(series_integrator)
  end

  @testset "Symbolics.jl" begin
    Symbolics.@variables α
    A = [0 0; 1/(2*α) 0]; b = [1-α, α]; c = [0, 1/(2*α)]
    series_integrator = bseries(A, b, c, 3)
    @test_nowarn latexify(series_integrator)
  end
end

@testset "AbstractDict interface" begin
  # explicit midpoint method
  A = @SArray [0 0; 1//2 0]; b = @SArray [0, 1//1]; c = @SArray [0, 1//2];
  series_integrator = bseries(A, b, c, 2)

  # These are mostly simple smoke tests
  @test_nowarn begin
    @test isempty(empty(series_integrator, RootedTrees.RootedTree{Int32, Vector{Int32}}, Rational{Int32}))
  end
  empty!(series_integrator)
  @test isempty(series_integrator)
  @test isempty(keys(series_integrator))
  @test isempty(values(series_integrator))
  @test length(series_integrator) == 0
  @test get(series_integrator, rootedtree([1]), :default) == :default
  @test get(series_integrator, rootedtree([1])) do
    :default
  end == :default
  @test getkey(series_integrator, rootedtree([1]), :default) == :default
  @test_nowarn sizehint!(series_integrator, 10)

  @test get!(series_integrator, rootedtree([1]), 2) == 2
  @test pop!(series_integrator, rootedtree([1])) == 2
  @test pop!(series_integrator, rootedtree([1]), :default) == :default
  series_integrator[rootedtree([1])] = 3
  @test isempty(delete!(series_integrator, rootedtree([1])))
end

@testset "substitute" begin
  Symbolics.@variables a1 a2 a31 a32
  a = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Symbolics.Num}(
        rootedtree(Int[])     => 1,
        rootedtree([1])       => a1,
        rootedtree([1, 2])    => a2,
        rootedtree([1, 2, 2]) => a31,
        rootedtree([1, 2, 3]) => a32)

  Symbolics.@variables b1 b2 b31 b32
  b = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Symbolics.Num}(
        rootedtree(Int[])     => 0,
        rootedtree([1])       => b1,
        rootedtree([1, 2])    => b2,
        rootedtree([1, 2, 2]) => b31,
        rootedtree([1, 2, 3]) => b32)

  # See equation (6) of
  # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
  #   Numerical integrators based on modified differential equations
  #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)

  ba = @inferred substitute(b, a)
  t = rootedtree(Int[])
  @test isequal(ba[t], a[t])

  t = rootedtree([1])
  coef = a1 * b1
  @inferred substitute(b, a, t)
  @test isequal(substitute(b, a, t), coef)
  @test isequal(ba[t], coef)

  t = rootedtree([1, 2])
  coef = a1 * b2 + a2 * b1^2
  @inferred substitute(b, a, t)
  @test isequal(substitute(b, a, t), coef)
  @test isequal(ba[t], coef)

  t = rootedtree([1, 2, 2])
  coef = a1 * b31 + 2 * a2 * b1 * b2 + a31 * b1^3
  @inferred substitute(b, a, t)
  @test isequal(substitute(b, a, t), coef)
  @test isequal(ba[t], coef)

  t = rootedtree([1, 2, 3])
  coef = a1 * b32 + 2 * a2 * b1 * b2 + a32 * b1^3
  @inferred substitute(b, a, t)
  @test isequal(substitute(b, a, t), coef)
  @test isequal(ba[t], coef)
end

@testset "compose" begin
  Symbolics.@variables a0 a1 a2 a31 a32
  a = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Symbolics.Num}(
        rootedtree(Int[])     => a0,
        rootedtree([1])       => a1,
        rootedtree([1, 2])    => a2,
        rootedtree([1, 2, 2]) => a31,
        rootedtree([1, 2, 3]) => a32)

  Symbolics.@variables b1 b2 b31 b32
  b = OrderedDict{RootedTrees.RootedTree{Int, Vector{Int}}, Symbolics.Num}(
        rootedtree(Int[])     => 0,
        rootedtree([1])       => b1,
        rootedtree([1, 2])    => b2,
        rootedtree([1, 2, 2]) => b31,
        rootedtree([1, 2, 3]) => b32)

  # See equation (5) of
  # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
  #   Numerical integrators based on modified differential equations
  #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)
  t = rootedtree([1])
  @inferred compose(b, a, t)
  @test isequal(compose(b, a, t), a0 * b1 + a1)

  t = rootedtree([1, 2])
  @inferred compose(b, a, t)
  @test isequal(compose(b, a, t), a0 * b2 + a1 * b1 + a2)

  t = rootedtree([1, 2, 2])
  @inferred compose(b, a, t)
  @test isequal(compose(b, a, t), a0 * b31 + a1 * b1^2 + 2 * a2 * b1 + a31)

  t = rootedtree([1, 2, 3])
  @inferred compose(b, a, t)
  @test isequal(compose(b, a, t), a0 * b32 + a1 * b2 + a2 * b1 + a32)

  @testset "Composing RK4 with itself" begin
    # classical fourth-order Runge-Kutta method
    A = @SArray [0 0 0 0;
                 1//2 0 0 0;
                 0 1//2 0 0;
                 0 0 1 0]
    b = @SArray [1//6, 1//3, 1//3, 1//6]
    c = @SArray [0, 1//2, 1//2, 1]

    series_rk4 = bseries(A, b, c, 8)
    series_rk4_composed = @inferred compose(series_rk4, series_rk4,
                                            normalize_stepsize=true)

    # Butcher coefficients of the composition (with normalized step size)
    A = @SArray [0     0     0     0     0     0     0     0;
                 1//4  0     0     0     0     0     0     0;
                 0     1//4  0     0     0     0     0     0;
                 0     0     1//2  0     0     0     0     0;
                 1//12 1//6  1//6  1//12 0     0     0     0;
                 1//12 1//6  1//6  1//12 1//4  0     0     0;
                 1//12 1//6  1//6  1//12 0     1//4  0     0;
                 1//12 1//6  1//6  1//12 0     0     1//2  0]
    b = @SArray [1//12, 1//6, 1//6, 1//12, 1//12, 1//6, 1//6, 1//12]
    c = @SArray [0, 1//4, 1//4, 1//2, 1//2, 3//4, 3//4, 1//1]

    series_2rk4 = bseries(A, b, c, order(series_rk4))

    @test series_rk4_composed == series_2rk4

    # Ensure that the normalization uses the correct factor
    unnormalized_series = @inferred compose(series_rk4, series_rk4,
                                            normalize_stepsize=false)
    for t in keys(unnormalized_series)
      @test unnormalized_series[t] == series_rk4_composed[t] * 2^order(t)
    end
  end
end

@testset "modified_equation" begin
  t1  = rootedtree([1])
  t2  = rootedtree([1, 2])
  t31 = rootedtree([1, 2, 2])
  t32 = rootedtree([1, 2, 3])
  t41 = rootedtree([1, 2, 2, 2])
  t42 = rootedtree([1, 2, 2, 3])
  t43 = rootedtree([1, 2, 3, 3])
  t44 = rootedtree([1, 2, 3, 4])
  t51 = rootedtree([1, 2, 2, 2, 2])
  t52 = rootedtree([1, 2, 2, 2, 3])
  t53 = rootedtree([1, 2, 2, 3, 3])
  t54 = rootedtree([1, 2, 2, 3, 4])
  t55 = rootedtree([1, 2, 3, 2, 3])
  t56 = rootedtree([1, 2, 3, 3, 3])
  t57 = rootedtree([1, 2, 3, 3, 4])
  t58 = rootedtree([1, 2, 3, 4, 4])
  t59 = rootedtree([1, 2, 3, 4, 5])

  @testset "SSPRK(2, 2) aka Heun's method" begin
    A = @SArray [0.0 0.0;
                1.0 0.0]
    b = @SArray [1/2, 1/2]
    c = @SArray [0.0, 1.0]
    order = 5
    series = modified_equation(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] ≈ 1.0                   atol=10*eps()
    @test series[t2 ] ≈ 0                     atol=10*eps()
    @test series[t31] ≈ 0.166666666666667     atol=10*eps()
    @test series[t32] ≈ -0.166666666666667    atol=10*eps()
    @test series[t41] ≈ -1.11022302462516e-16 atol=10*eps()
    @test series[t42] ≈ -0.125000000000000    atol=10*eps()
    @test series[t43] ≈ -2.77555756156289e-17 atol=10*eps()
    @test series[t44] ≈ 0.125000000000000     atol=10*eps()
    @test series[t51] ≈ -0.0333333333333332   atol=10*eps()
    @test series[t52] ≈ -0.0583333333333333   atol=10*eps()
    @test series[t55] ≈ 0.0750000000000000    atol=10*eps()
    @test series[t53] ≈ 0.0583333333333333    atol=10*eps()
    @test series[t56] ≈ 0.0333333333333334    atol=10*eps()
    @test series[t54] ≈ 0.0500000000000000    atol=10*eps()
    @test series[t57] ≈ 0.0583333333333333    atol=10*eps()
    @test series[t58] ≈ -0.0583333333333333   atol=10*eps()
    @test series[t59] ≈ -0.0500000000000000   atol=10*eps()
  end

  @testset "SSPRK(3, 3)" begin
    A = @SArray [0.0 0.0 0.0;
                1.0 0.0 0.0;
                1/4 1/4 0.0]
    b = @SArray [1/6, 1/6, 2/3]
    c = @SArray [0.0, 1.0, 1/2]
    order = 5
    series = modified_equation(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] ≈ 1.0                  atol=10*eps()
    @test series[t2 ] ≈ 0                    atol=10*eps()
    @test series[t31] ≈ 0                    atol=10*eps()
    @test series[t32] ≈ 0                    atol=10*eps()
    @test series[t41] ≈ 0                    atol=10*eps()
    @test series[t42] ≈ -0.0416666666666667  atol=10*eps()
    @test series[t43] ≈ 0.0833333333333333   atol=10*eps()
    @test series[t44] ≈ -0.0416666666666667  atol=10*eps()
    @test series[t51] ≈ 0.00833333333333330  atol=10*eps()
    @test series[t52] ≈ -0.0166666666666667  atol=10*eps()
    @test series[t55] ≈ 0.0333333333333333   atol=10*eps()
    @test series[t53] ≈ 0.0166666666666667   atol=10*eps()
    @test series[t56] ≈ -0.00833333333333333 atol=10*eps()
    @test series[t54] ≈ 0.00833333333333333  atol=10*eps()
    @test series[t57] ≈ -0.0250000000000000  atol=10*eps()
    @test series[t58] ≈ -0.0166666666666667  atol=10*eps()
    @test series[t59] ≈ 0.0333333333333333   atol=10*eps()
  end

  @testset "Explicit midpoint method (Runge's method)" begin
    A = @SArray [0 0; 1//2 0]
    b = @SArray [0, 1//1]
    c = @SArray [0, 1//2]
    order = 5
    series = modified_equation(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] == 1
    @test series[t2 ] == 0
    @test series[t31] == -1//12
    @test series[t32] == -1//6
    @test series[t41] == 0
    @test series[t42] == 0
    @test series[t43] == 1//8
    @test series[t44] == 1//8
    @test series[t51] == 7//240
    @test series[t52] == 1//40
    @test series[t55] == 1//30
    @test series[t53] == 3//80
    @test series[t56] == -7//240
    @test series[t54] == 7//240
    @test series[t57] == -1//40
    @test series[t58] == -19//240
    @test series[t59] == -1//20
  end
end

@testset "elementary differentials" begin
  @testset "Bicolored trees" begin
    @testset "Lotka-Volterra" begin
      # Verified with Mathematica
      # u = {q, p};
      # f1 = {q * (p - 1), 0};
      # f2 = {0, p * (2 - q)};
      #
      # f1p = Simplify@D[f1, {u}];
      # f2p = Simplify@D[f2, {u}];
      #
      # f1pp = Simplify@D[f1, {u, 2}];
      # f2pp = Simplify@D[f2, {u, 2}];
      #
      # (* Trees of order 2 *)
      # Map[InputForm, f1p.f1]
      # Map[InputForm, f2p.f1]
      # Map[InputForm, f1p.f2]
      # Map[InputForm, f2p.f2]
      #
      # (* Tall trees of order 3 *)
      # Map[InputForm, f1p.(f1p.f1)]
      # Map[InputForm, f2p.(f1p.f1)]
      # Map[InputForm, f1p.(f2p.f1)]
      # Map[InputForm, f2p.(f2p.f1)]
      # Map[InputForm, f1p.(f1p.f2)]
      # Map[InputForm, f2p.(f1p.f2)]
      # Map[InputForm, f1p.(f2p.f2)]
      # Map[InputForm, f2p.(f2p.f2)]
      #
      # (* Bushy trees of order 3 *)
      # Map[InputForm, f1pp.f1.f1]
      # Map[InputForm, f2pp.f1.f1]
      # Map[InputForm, f1pp.f2.f1]
      # Map[InputForm, f2pp.f2.f1]
      # (* The other choice of colors of the leaves is redundant *)
      # Map[InputForm, f1pp.f2.f2]
      # Map[InputForm, f2pp.f2.f2]
      @testset "SymEngine" begin
        q, p = u = SymEngine.symbols("q, p")
        f = ([q * (p - 1), 0], [0, p * (2 - q)])
        differentials = elementary_differentials(f, u, 3)
        @test differentials[rootedtree(Int64[], Bool[])         ] == [1, 1]
        @test differentials[rootedtree([1], Bool[0])            ] == [q * (p - 1), 0]
        @test differentials[rootedtree([1], Bool[1])            ] == [0, p * (2 - q)]
        @test differentials[rootedtree([1, 2], Bool[0, 0])      ] == [(-1 + p)^2*q,0]
        @test differentials[rootedtree([1, 2], Bool[1, 0])      ] == [0,-((-1 + p)*p*q)]
        @test differentials[rootedtree([1, 2], Bool[0, 1])      ] == [p*(2 - q)*q,0]
        @test differentials[rootedtree([1, 2], Bool[1, 1])      ] == [0,p*(2 - q)^2]
        @test differentials[rootedtree([1, 2, 3], Bool[0, 0, 0])] == [(-1 + p)^3*q,0]
        @test differentials[rootedtree([1, 2, 3], Bool[1, 0, 0])] == [0,-((-1 + p)^2*p*q)]
        @test differentials[rootedtree([1, 2, 3], Bool[0, 1, 0])] == [-((-1 + p)*p*q^2),0]
        @test differentials[rootedtree([1, 2, 3], Bool[1, 1, 0])] == [0,-((-1 + p)*p*(2 - q)*q)]
        @test differentials[rootedtree([1, 2, 3], Bool[0, 0, 1])] == [(-1 + p)*p*(2 - q)*q,0]
        @test differentials[rootedtree([1, 2, 3], Bool[1, 0, 1])] == [0,-(p^2*(2 - q)*q)]
        @test differentials[rootedtree([1, 2, 3], Bool[0, 1, 1])] == [p*(2 - q)^2*q,0]
        @test differentials[rootedtree([1, 2, 3], Bool[1, 1, 1])] == [0,p*(2 - q)^3]
        @test differentials[rootedtree([1, 2, 2], Bool[0, 0, 0])] == [0, 0]
        @test differentials[rootedtree([1, 2, 2], Bool[1, 0, 0])] == [0, 0]
        @test differentials[rootedtree([1, 2, 2], Bool[0, 1, 0])] == [(-1 + p)*p*(2 - q)*q,0]
        @test differentials[rootedtree([1, 2, 2], Bool[1, 1, 0])] == [0,-((-1 + p)*p*(2 - q)*q)]
        @test differentials[rootedtree([1, 2, 2], Bool[0, 1, 1])] == [0, 0]
        @test differentials[rootedtree([1, 2, 2], Bool[1, 1, 1])] == [0, 0]
      end
    end
  end
end

@testset "modified_equation with elementary differentials" begin
  @testset "SymEngine.jl" begin
    @testset "Explicit Euler" begin
      # Lotka-Volterra model
      dt = SymEngine.symbols("dt")
      p, q = u = SymEngine.symbols("p, q")
      f = [p * (2 - q), q * (p - 1)]

      # Explicit Euler method
      A = @SArray [0//1;]
      b = @SArray [1//1]
      c = @SArray [0//1]

      # tested with the Python package BSeries
      series2 = modified_equation(f, u, dt, A, b, c, 2)
      series2_reference = [
        -dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
        -dt*( p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
      ]
      @test mapreduce(isequal, &, series2, series2_reference)

      # tested with the Python package BSeries
      series3 = modified_equation(f, u, dt, A, b, c, 3)
      series3_reference = [
        -dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
        dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
      ]
      @test mapreduce(iszero ∘ SymEngine.expand, &, series3 - series3_reference)

      # tested with the Python package BSeries
      series4 = modified_equation(f, u, dt, A, b, c, 4)
      series4_reference = [
        -dt^3*(-2*p^2*q*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(-p^2*q*(2 - q)^2 + p*q^2*(p - 1)^2 + p*q*(1 - p)*(2 - q)*(p - 1) + p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(p^2*q^2*(p - 1) - 2*p^2*q*(2 - q)^2 - p^2*q*(2 - q)*(p - 1) - p*q*(2 - q)^2*(p - 1) - p*q*(2 - q)*(p - 1)^2 - p*q*(p - 1)^3 + p*(2 - q)^4)/4 - dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
          -dt^3*(-2*p*q^2*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(p^2*q*(2 - q)^2 - p*q^2*(p - 1)^2 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(-p^2*q^2*(2 - q) - p*q^2*(2 - q)*(p - 1) - 2*p*q^2*(p - 1)^2 + p*q*(2 - q)^3 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2 + q*(p - 1)^4)/4 + dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
      ]
      @test mapreduce(iszero ∘ SymEngine.expand, &, series4 - series4_reference)
    end

    @testset "IMEX Euler (partitioned Euler, symplectic Euler)" begin
      # Lotka-Volterra model
      dt = SymEngine.symbols("dt")
      q, p = u = SymEngine.symbols("q, p")
      f = ([0, p * (2 - q)], [q * (p - 1), 0])

      # Symplectic Euler method
      ex_euler = RungeKuttaMethod(
        @SMatrix([0]), @SVector [1]
      )
      im_euler = RungeKuttaMethod(
        @SMatrix([1]), @SVector [1]
      )
      ark = AdditiveRungeKuttaMethod([im_euler, ex_euler])

      # Hairer, Lubich, Wanner (2006) Geometric numerical integration
      # Example IX.1.3(b)
      series_integrator = bseries(ark, 2)
      series = modified_equation(f, u, dt, series_integrator)
      series_reference = [
        q*(p - 1) - dt / 2 * q * (p^2 + p * q - 4 * p + 1),
        -p * (q - 2) + dt / 2 * p * (q^2 + p * q - 5 * q + 4)
      ]
      @test mapreduce(iszero ∘ SymEngine.expand, &, series - series_reference)
    end
  end

  @testset "SymPy.jl" begin
    # Lotka-Volterra model
    dt = SymPy.symbols("dt")
    p, q = u = SymPy.symbols("p, q")
    f = [p * (2 - q), q * (p - 1)]

    # Explicit Euler method
    A = @SArray [0//1;]
    b = @SArray [1//1]
    c = @SArray [0//1]

    # tested with the Python package BSeries
    series2 = modified_equation(f, u, dt, A, b, c, 2)
    series2_reference = [
      -dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
      -dt*( p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(isequal, &, series2, series2_reference)

    # tested with the Python package BSeries
    series3 = modified_equation(f, u, dt, A, b, c, 3)
    series3_reference = [
      -dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
      dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(iszero ∘ SymPy.expand, &, series3 - series3_reference)

    # tested with the Python package BSeries
    series4 = modified_equation(f, u, dt, A, b, c, 4)
    series4_reference = [
      -dt^3*(-2*p^2*q*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(-p^2*q*(2 - q)^2 + p*q^2*(p - 1)^2 + p*q*(1 - p)*(2 - q)*(p - 1) + p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(p^2*q^2*(p - 1) - 2*p^2*q*(2 - q)^2 - p^2*q*(2 - q)*(p - 1) - p*q*(2 - q)^2*(p - 1) - p*q*(2 - q)*(p - 1)^2 - p*q*(p - 1)^3 + p*(2 - q)^4)/4 - dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
        -dt^3*(-2*p*q^2*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(p^2*q*(2 - q)^2 - p*q^2*(p - 1)^2 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(-p^2*q^2*(2 - q) - p*q^2*(2 - q)*(p - 1) - 2*p*q^2*(p - 1)^2 + p*q*(2 - q)^3 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2 + q*(p - 1)^4)/4 + dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(iszero ∘ SymPy.expand, &, series4 - series4_reference)
  end

  @testset "Symbolics.jl" begin
    # Lotka-Volterra model
    Symbolics.@variables dt
    u = Symbolics.@variables p q
    f = [p * (2 - q), q * (p - 1)]

    # Explicit Euler method
    A = @SArray [0//1;]
    b = @SArray [1//1]
    c = @SArray [0//1]

    # tested with the Python package BSeries
    series2 = modified_equation(f, u, dt, A, b, c, 2)
    series2_reference = [
      -dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
      -dt*( p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(isequal, &, series2, series2_reference)

    # tested with the Python package BSeries
    series3 = modified_equation(f, u, dt, A, b, c, 3)
    series3_reference = [
      -dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
      dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(iszero ∘ Symbolics.expand, &, series3 - series3_reference)

    # tested with the Python package BSeries
    series4 = modified_equation(f, u, dt, A, b, c, 4)
    series4_reference = [
      -dt^3*(-2*p^2*q*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(-p^2*q*(2 - q)^2 + p*q^2*(p - 1)^2 + p*q*(1 - p)*(2 - q)*(p - 1) + p*q*(2 - q)*(p - 1)*(q - 2))/12 - dt^3*(p^2*q^2*(p - 1) - 2*p^2*q*(2 - q)^2 - p^2*q*(2 - q)*(p - 1) - p*q*(2 - q)^2*(p - 1) - p*q*(2 - q)*(p - 1)^2 - p*q*(p - 1)^3 + p*(2 - q)^4)/4 - dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p^2*q*(2 - q) - p*q*(2 - q)*(p - 1) - p*q*(p - 1)^2 + p*(2 - q)^3)/3 - dt*(-p*q*(p - 1) + p*(2 - q)^2)/2 + p*(2 - q),
        -dt^3*(-2*p*q^2*(2 - q)*(p - 1) + 2*p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(p^2*q*(2 - q)^2 - p*q^2*(p - 1)^2 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2)/12 - dt^3*(-p^2*q^2*(2 - q) - p*q^2*(2 - q)*(p - 1) - 2*p*q^2*(p - 1)^2 + p*q*(2 - q)^3 + p*q*(2 - q)^2*(p - 1) + p*q*(2 - q)*(p - 1)^2 + q*(p - 1)^4)/4 + dt^2*p*q*(2 - q)*(p - 1)/6 + dt^2*(-p*q^2*(p - 1) + p*q*(2 - q)^2 + p*q*(2 - q)*(p - 1) + q*(p - 1)^3)/3 - dt*(p*q*(2 - q) + q*(p - 1)^2)/2 + q*(p - 1)
    ]
    @test mapreduce(iszero ∘ Symbolics.expand, &, series4 - series4_reference)
  end
end

@testset "modifying_integrator" begin
  t1  = rootedtree([1])
  t2  = rootedtree([1, 2])
  t31 = rootedtree([1, 2, 2])
  t32 = rootedtree([1, 2, 3])
  t41 = rootedtree([1, 2, 2, 2])
  t42 = rootedtree([1, 2, 2, 3])
  t43 = rootedtree([1, 2, 3, 3])
  t44 = rootedtree([1, 2, 3, 4])
  t51 = rootedtree([1, 2, 2, 2, 2])
  t52 = rootedtree([1, 2, 2, 2, 3])
  t53 = rootedtree([1, 2, 2, 3, 3])
  t54 = rootedtree([1, 2, 2, 3, 4])
  t55 = rootedtree([1, 2, 3, 2, 3])
  t56 = rootedtree([1, 2, 3, 3, 3])
  t57 = rootedtree([1, 2, 3, 3, 4])
  t58 = rootedtree([1, 2, 3, 4, 4])
  t59 = rootedtree([1, 2, 3, 4, 5])

  @testset "SSPRK(2, 2) aka Heun's method" begin
    A = @SArray [0.0 0.0;
                1.0 0.0]
    b = @SArray [1/2, 1/2]
    c = @SArray [0.0, 1.0]
    order = 5
    series = modifying_integrator(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] ≈ 1                    atol=10*eps()
    @test series[t2 ] ≈ 0                    atol=10*eps()
    @test series[t31] ≈ -0.166666666666667   atol=10*eps()
    @test series[t32] ≈ 1/6                  atol=10*eps()
    @test series[t41] ≈ 8.32667268468867e-17 atol=10*eps()
    @test series[t42] ≈ 0.125000000000000    atol=10*eps()
    @test series[t43] ≈ 1.38777878078145e-17 atol=10*eps()
    @test series[t44] ≈ -0.125000000000000   atol=10*eps()
    @test series[t51] ≈ 0.200000000000000    atol=10*eps()
    @test series[t52] ≈ 0.0583333333333333   atol=10*eps()
    @test series[t55] ≈ 0.00833333333333335  atol=10*eps()
    @test series[t53] ≈ -0.0583333333333333  atol=10*eps()
    @test series[t56] ≈ -0.200000000000000   atol=10*eps()
    @test series[t54] ≈ -0.133333333333333   atol=10*eps()
    @test series[t57] ≈ -0.0583333333333333  atol=10*eps()
    @test series[t58] ≈ 0.0583333333333333   atol=10*eps()
    @test series[t59] ≈ 0.133333333333333    atol=10*eps()
  end

  @testset "SSPRK(3, 3)" begin
    A = @SArray [0.0 0.0 0.0;
                1.0 0.0 0.0;
                1/4 1/4 0.0]
    b = @SArray [1/6, 1/6, 2/3]
    c = @SArray [0.0, 1.0, 1/2]
    order = 5
    series = modifying_integrator(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] ≈ 1                    atol=10*eps()
    @test series[t2 ] ≈ 0                    atol=10*eps()
    @test series[t31] ≈ 0                    atol=10*eps()
    @test series[t32] ≈ 0                    atol=10*eps()
    @test series[t41] ≈ 0                    atol=10*eps()
    @test series[t42] ≈ 0.0416666666666667   atol=10*eps()
    @test series[t43] ≈ -0.0833333333333333  atol=10*eps()
    @test series[t44] ≈ 1/24                 atol=10*eps()
    @test series[t51] ≈ -0.00833333333333330 atol=10*eps()
    @test series[t52] ≈ 0.0166666666666667   atol=10*eps()
    @test series[t55] ≈ -0.0333333333333333  atol=10*eps()
    @test series[t53] ≈ -0.0166666666666667  atol=10*eps()
    @test series[t56] ≈ 0.00833333333333332  atol=10*eps()
    @test series[t54] ≈ -0.00833333333333334 atol=10*eps()
    @test series[t57] ≈ 0.0250000000000000   atol=10*eps()
    @test series[t58] ≈ 1/60                 atol=10*eps()
    @test series[t59] ≈ -0.0333333333333333  atol=10*eps()
  end

  @testset "Explicit midpoint method (Runge's method)" begin
    A = @SArray [0 0; 1//2 0]
    b = @SArray [0, 1//1]
    c = @SArray [0, 1//2]
    order = 5
    series = modifying_integrator(A, b, c, order)

    # tested with the Python package BSeries
    @test series[t1 ] == 1
    @test series[t2 ] == 0
    @test series[t31] == 1//12
    @test series[t32] == 1//6
    @test series[t41] == 0
    @test series[t42] == 0
    @test series[t43] == -1//8
    @test series[t44] == -1//8
    @test series[t51] == 1//80
    @test series[t52] == 1//60
    @test series[t55] == 7//240
    @test series[t53] == 1//240
    @test series[t56] == 9//80
    @test series[t54] == 1//80
    @test series[t57] == 13//120
    @test series[t58] == 13//80
    @test series[t59] == 2//15
  end
end

@testset "modifying_integrator with elementary differentials" begin
  @testset "SymEngine.jl" begin
    dt = SymEngine.symbols("dt")
    α, β, γ = SymEngine.symbols("alpha, beta, gamma")
    u1, u2, u3 = u = SymEngine.symbols("u1 u2 u3")
    f = [α*u[2]*u[3], β*u[3]*u[1], γ*u[1]*u[2]]

    # Implicit midpoint method
    A = @SArray [1//2;]
    b = @SArray [1//1]
    c = @SArray [1//2]

    # See eq. (12) of
    # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
    #   Numerical integrators based on modified differential equations
    #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)

    series = modifying_integrator(f, u, dt, A, b, c, 5)
    # Dirty workaround since there is no way to get the polynomial coefficients
    # at the moment.
    # We differentiate the expression and set `dt` to zero to get the corresponding
    # coefficient divided by factorial(how often we needed to differentiate).
    s3_reference = -(α * β * u3^2 + α * γ * u2^2 + β * γ * u1^2) / 12
    for i in eachindex(f)
      s3 = SymEngine.subs(
        1//2 * SymEngine.diff(SymEngine.diff((series[i] - f[i]) / f[i], dt), dt),
        dt => 0)
      @test iszero(SymEngine.expand(s3 - s3_reference))
    end

    s5_reference = 6//5 * s3_reference^2 + 1//60 * α * β * γ * (
      β * u1^2 * u3^2 + γ * u2^2 * u1^2 + α * u3^2 * u2^2
    )
    for i in eachindex(f)
      s5 = SymEngine.subs(
        1//24 * SymEngine.diff(SymEngine.diff(SymEngine.diff(SymEngine.diff(
          (series[i] - f[i]) / f[i], dt), dt), dt), dt),
        dt => 0)
      @test iszero(SymEngine.expand(s5 - s5_reference))
    end
  end

  @testset "SymPy.jl" begin
    dt = SymPy.symbols("dt")
    α, β, γ = SymPy.symbols("alpha, beta, gamma")
    u1, u2, u3 = u = SymPy.symbols("u1 u2 u3")
    f = [α*u[2]*u[3], β*u[3]*u[1], γ*u[1]*u[2]]

    # Implicit midpoint method
    A = @SArray [1//2;]
    b = @SArray [1//1]
    c = @SArray [1//2]

    # See eq. (12) of
    # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
    #   Numerical integrators based on modified differential equations
    #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)

    series = modifying_integrator(f, u, dt, A, b, c, 5)
    # Dirty workaround used also for the other symbolic setups - just make
    # it consistent here, although we could use another approach with SymPy.jl.
    # We differentiate the expression and set `dt` to zero to get the corresponding
    # coefficient divided by factorial(how often we needed to differentiate).
    s3_reference = -(α * β * u3^2 + α * γ * u2^2 + β * γ * u1^2) / 12
    for i in eachindex(f)
      s3 = SymPy.subs(
        1//2 * SymPy.diff(SymPy.diff((series[i] - f[i]) / f[i], dt), dt),
        dt => 0)
      @test iszero(SymPy.expand(s3 - s3_reference))
    end

    s5_reference = 6//5 * s3_reference^2 + 1//60 * α * β * γ * (
      β * u1^2 * u3^2 + γ * u2^2 * u1^2 + α * u3^2 * u2^2
    )
    for i in eachindex(f)
      s5 = SymPy.subs(
        1//24 * SymPy.diff(SymPy.diff(SymPy.diff(SymPy.diff(
          (series[i] - f[i]) / f[i], dt), dt), dt), dt),
        dt => 0)
      @test iszero(SymPy.expand(s5 - s5_reference))
    end
  end

  @testset "Symbolics.jl" begin
    Symbolics.@variables dt
    Symbolics.@variables α β γ
    u = Symbolics.@variables u1 u2 u3
    f = [α*u[2]*u[3], β*u[3]*u[1], γ*u[1]*u[2]]

    # Implicit midpoint method
    A = @SArray [1//2;]
    b = @SArray [1//1]
    c = @SArray [1//2]

    # See eq. (12) of
    # - Philippe Chartier, Ernst Hairer and Gilles Vilmart (2007)
    #   Numerical integrators based on modified differential equations
    #   [DOI: 10.1090/S0025-5718-07-01967-9](https://doi.org/10.1090/S0025-5718-07-01967-9)

    series = modifying_integrator(f, u, dt, A, b, c, 5)
    # Dirty workaround since there is no way to get the polynomial coefficients
    # in Symbolics.jl at the moment, see
    # https://github.com/JuliaSymbolics/Symbolics.jl/issues/216
    # We differentiate the expression and set `dt` to zero to get the corresponding
    # coefficient divided by factorial(how often we needed to differentiate).
    d = Symbolics.Differential(dt)
    s3_reference = -(α * β * u3^2 + α * γ * u2^2 + β * γ * u1^2) / 12
    for i in eachindex(f)
      s3 = Symbolics.simplify_fractions(Symbolics.substitute(
        Symbolics.expand_derivatives(1//2 * d(d((series[i] - f[i]) / f[i]))),
        dt => 0))
      @test isequal(s3, s3_reference)
    end

    s5_reference = 6//5 * s3_reference^2 + 1//60 * α * β * γ * (
      β * u1^2 * u3^2 + γ * u2^2 * u1^2 + α * u3^2 * u2^2
    )
    for i in eachindex(f)
      s5 = Symbolics.simplify_fractions(Symbolics.substitute(
        Symbolics.expand_derivatives(1//24 * d(d(d(d((series[i] - f[i]) / f[i]))))),
        dt => 0))
      @test iszero(Symbolics.expand(s5 - s5_reference))
    end
  end
end

@testset "nonlinear oscillator" begin
  @testset "SymEngine.jl" begin
    dt = SymEngine.symbols("dt")
    u = SymEngine.symbols("u_1, u_2")
    f = [-u[2], u[1]] / (u[1]^2 + u[2]^2)

    # explicit midpoint method
    A = @SArray [0 0; 1//2 0]
    b = @SArray [0, 1//1]
    c = @SArray [0, 1//2];

    # tested with Mathematica using
    # ```
    # ClearAll["Global`*"];
    #
    # solution[t_,
    #   u0_] := {u0[[1]]*Cos[t/(u0[[1]]^2 + u0[[2]]^2)] -
    #   u0[[2]]*Sin[t/(u0[[1]]^2 + u0[[2]]^2)],
    #   u0[[2]]*Cos[t/(u0[[1]]^2 + u0[[2]]^2)] +
    #   u0[[1]]*Sin[t/(u0[[1]]^2 + u0[[2]]^2)]}
    #
    # factorF2[u_] := 1/(12*(u[[1]]^2 + u[[2]]^2)^2)
    # (*factorF2[u_]:=ff2*)
    # factorF4[u_] := 1/(20*(u[[1]]^2 + u[[2]]^2)^4)
    # (*factorF4[u_]:=ff4*)
    # factorF6[u_] := 127/(2016*(u[[1]]^2 + u[[2]]^2)^6)
    # (*factorF6[u_]:=ff6*)
    # factorF8[u_] := 1513/(12960*(u[[1]]^2 + u[[2]]^2)^8)
    # (*factorF8[u_]:=ff8*)
    #
    # factorU3[u_] := 0
    # (*factorU3[u_]:=fu3*)
    # factorU5[u_] := 1/(48*(u[[1]]^2 + u[[2]]^2)^6)
    # (*factorU5[u_]:=fu5*)
    # factorU7[u_] := 31/(640*(u[[1]]^2 + u[[2]]^2)^8)
    # (*factorU7[u_]:=fu7*)
    # factorU9[u_] := 8969/(80640*(u[[1]]^2 + u[[2]]^2)^10)
    # (*factorU9[u_]:=fu9*)
    #
    # newModifiedF[
    #   u_] :=
    # (1 + factorF2[u]*h^2 + factorF4[u]*h^4 + factorF6[u]*h^6)*
    #   originalF[
    #     u] + (factorU3[u]*h^3 + factorU5[u]*h^5 + factorU7[u]*h^7)*u
    #
    # u0 = {u01, u02};
    # f = newModifiedF;
    # y2 = u0 + h/2*f[u0];
    # unew = u0 + h*f[y2];
    # difference = Simplify@Series[unew - solution[h, u0], {h, 0, 10}];
    #
    # Simplify[difference[[1]]]
    # Simplify[difference[[2]]]
    # ```
    series = modifying_integrator(f, u, dt, A, b, c, 10)
    terms = SymEngine.subs.(series, (Dict(u[1] => 1//1, u[2] => 0//1), ))

    @test isequal(terms[1], 1//48 * dt^5 + 31//640 * dt^7 + 8969//80640 * dt^9)
    @test isequal(terms[2], 1 + 1//12 * dt^2 + 1//20 * dt^4 + 127//2016 * dt^6 + 1513//12960 * dt^8)
  end
end

@testset "integer coefficients" begin
  # reported by David Ketcheson on 2021-12-09
  A = reshape([1], 1, 1); b = [1]; c = [1]
  @inferred modified_equation(A, b, c, 4)
end

@testset "Runge-Kutta methods interface" begin
  # classical RK4
  A = [0 0 0 0
       1//2 0 0 0
       0 1//2 0 0
       0 0 1 0]
  b = [1//6, 1//3, 1//3, 1//6]
  rk = RungeKuttaMethod(A, b)

  # fourth-order accurate
  series_integrator = @inferred bseries(rk, 4)
  series_exact = @inferred ExactSolution(series_integrator)
  @test mapreduce(==, &, values(series_integrator), values(series_exact))

  # not fifth-order accurate
  series_integrator = @inferred bseries(rk, 5)
  series_exact = @inferred ExactSolution(series_integrator)
  @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false
end

@testset "additive Runge-Kutta methods interface" begin
  @testset "IMEX Euler (partitioned Euler, symplectic Euler)" begin
    ex_euler = RungeKuttaMethod(
      @SMatrix([0]), @SVector [1]
    )
    im_euler = RungeKuttaMethod(
      @SMatrix([1]), @SVector [1]
    )
    ark = AdditiveRungeKuttaMethod([im_euler, ex_euler])

    # first-order accurate
    series_integrator = @inferred bseries(ark, 1)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact))

    # not second-order accurate
    series_integrator = @inferred bseries(ark, 2)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false

    # modified equations and modifying integrators
    let series_integrator = @inferred bseries(ark, 3)
      @testset "modified_equation" begin
        mod_eq = @inferred modified_equation(series_integrator)

        # Hairer, Lubich, Wanner (2006) Geometric numerical integration
        # Table IX.10.1, p. 383
        # Black nodes are `0`, white nodes are `1`
        mod_eq_reference = Dict(
          rootedtree(Int[], Bool[])            => 0//1,
          rootedtree([1], Bool[0])             => 1//1,
          rootedtree([1], Bool[1])             => 1//1,
          rootedtree([1, 2], Bool[0, 0])       => 1//2,
          rootedtree([1, 2], Bool[1, 0])       => 1//2,
          rootedtree([1, 2], Bool[0, 1])       => -1//2,
          rootedtree([1, 2], Bool[1, 1])       => -1//2,
          rootedtree([1, 2, 3], Bool[0, 0, 0]) => 1//3,
          rootedtree([1, 2, 3], Bool[1, 0, 0]) => 1//3,
          rootedtree([1, 2, 3], Bool[0, 1, 0]) => -1//6,
          rootedtree([1, 2, 3], Bool[1, 1, 0]) => -1//6,
          rootedtree([1, 2, 3], Bool[0, 0, 1]) => -1//6,
          rootedtree([1, 2, 3], Bool[1, 0, 1]) => -1//6,
          rootedtree([1, 2, 3], Bool[0, 1, 1]) => 1//3,
          rootedtree([1, 2, 3], Bool[1, 1, 1]) => 1//3,
          rootedtree([1, 2, 2], Bool[0, 0, 0]) => 1//6,
          rootedtree([1, 2, 2], Bool[1, 0, 0]) => 1//6,
          rootedtree([1, 2, 2], Bool[0, 1, 0]) => -1//3,
          rootedtree([1, 2, 2], Bool[1, 1, 0]) => -1//3,
          rootedtree([1, 2, 2], Bool[0, 1, 1]) => 1//6,
          rootedtree([1, 2, 2], Bool[1, 1, 1]) => 1//6,
        )
        for (key, val) in mod_eq
          @test val == mod_eq_reference[key]
        end
      end

      @testset "modifying_integrator" begin
        mod_int = @inferred modifying_integrator(series_integrator)
      end
    end
  end

  @testset "Störmer-Verlet" begin
    # Hairer, Lubich, Wanner (2002)
    # Geometric numerical integration
    # Table II.2.1
    As = [
      [0 0; 1//2 1//2],
      [1//2 0; 1//2 0]
    ]
    bs = [
      [1//2, 1//2],
      [1//2, 1//2]
    ]
    ark = AdditiveRungeKuttaMethod(As, bs)

    # second-order accurate
    series_integrator = @inferred bseries(ark, 2)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact))

    # not third-order accurate
    series_integrator = @inferred bseries(ark, 3)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false

    # modified equations and modifying integrators
    let series_integrator = @inferred bseries(ark, 4)
      @testset "modified_equation" begin
        mod_eq = @inferred modified_equation(series_integrator)

        # second-order accurate
        for t in BicoloredRootedTreeIterator(2)
          @test iszero(mod_eq[t])
        end

        # Hairer, Lubich, Wanner (2006) Geometric numerical integration
        # Table IX.10.2, p. 386
        # Black nodes are `0`, white nodes are `1`
        series_reference = Dict(
          rootedtree(Int[], Bool[])            => 1//1,
          rootedtree([1], Bool[0])             => 1//1,
          rootedtree([1], Bool[1])             => 1//1,
          rootedtree([1, 2], Bool[0, 0])       => 1//2,
          rootedtree([1, 2], Bool[1, 0])       => 1//2,
          rootedtree([1, 2], Bool[0, 1])       => 1//2,
          rootedtree([1, 2], Bool[1, 1])       => 1//2,
          rootedtree([1, 2, 3], Bool[0, 0, 0]) => 1//4,
          rootedtree([1, 2, 3], Bool[1, 0, 0]) => 1//4,
          rootedtree([1, 2, 3], Bool[0, 1, 0]) => 0//1,
          rootedtree([1, 2, 3], Bool[1, 1, 0]) => 0//1,
          rootedtree([1, 2, 3], Bool[0, 0, 1]) => 1//4,
          rootedtree([1, 2, 3], Bool[1, 0, 1]) => 1//4,
          rootedtree([1, 2, 3], Bool[0, 1, 1]) => 1//4,
          rootedtree([1, 2, 3], Bool[1, 1, 1]) => 1//4,
          rootedtree([1, 2, 2], Bool[0, 0, 0]) => 1//2,
          rootedtree([1, 2, 2], Bool[1, 0, 0]) => 1//2,
          rootedtree([1, 2, 2], Bool[0, 1, 0]) => 1//4,
          rootedtree([1, 2, 2], Bool[1, 1, 0]) => 1//4,
          rootedtree([1, 2, 2], Bool[0, 1, 1]) => 1//4,
          rootedtree([1, 2, 2], Bool[1, 1, 1]) => 1//4,
        )
        for t in keys(series_reference)
          @test series_reference[t] == series_integrator[t]
        end
        mod_eq_reference = Dict(
          rootedtree(Int[], Bool[])            => 0//1,
          rootedtree([1], Bool[0])             => 1//1,
          rootedtree([1], Bool[1])             => 1//1,
          rootedtree([1, 2], Bool[0, 0])       => 0//1,
          rootedtree([1, 2], Bool[1, 0])       => 0//1,
          rootedtree([1, 2], Bool[0, 1])       => 0//1,
          rootedtree([1, 2], Bool[1, 1])       => 0//1,
          rootedtree([1, 2, 3], Bool[0, 0, 0]) => 1//12,
          rootedtree([1, 2, 3], Bool[1, 0, 0]) => 1//12,
          rootedtree([1, 2, 3], Bool[0, 1, 0]) => -1//6,
          rootedtree([1, 2, 3], Bool[1, 1, 0]) => -1//6,
          rootedtree([1, 2, 3], Bool[0, 0, 1]) => 1//12,
          rootedtree([1, 2, 3], Bool[1, 0, 1]) => 1//12,
          rootedtree([1, 2, 3], Bool[0, 1, 1]) => 1//12,
          rootedtree([1, 2, 3], Bool[1, 1, 1]) => 1//12,
          rootedtree([1, 2, 2], Bool[0, 0, 0]) => 1//6,
          rootedtree([1, 2, 2], Bool[1, 0, 0]) => 1//6,
          rootedtree([1, 2, 2], Bool[0, 1, 0]) => -1//12,
          rootedtree([1, 2, 2], Bool[1, 1, 0]) => -1//12,
          rootedtree([1, 2, 2], Bool[0, 1, 1]) => -1//12,
          rootedtree([1, 2, 2], Bool[1, 1, 1]) => -1//12,
        )
        for t in keys(mod_eq_reference)
          @test mod_eq_reference[t] == mod_eq[t]
        end

        # Hairer, Lubich, Wanner (2003)
        # Geometric numerical integration illustrated by the Störmer-Verlet method
        # https://doi.org/10.1017/S0962492902000144
        # equation (4.8)
        @testset "Pendulum, SymEngine" begin
          dt = SymEngine.symbols("dt")
          q, v = u = SymEngine.symbols("q, v")
          f = ([v, 0], [0, -sin(q)])

          series = modified_equation(f, u, dt, series_integrator)
          series_reference = f[1] + f[2] + 1//12 * dt^2 * [2 * cos(q) * v, sin(q) * cos(q) + sin(q) * v^2]
          @test mapreduce(iszero ∘ SymEngine.expand, &, series - series_reference)
        end
      end # @testset "modified_equation"

      @testset "modifying_integrator" begin
        mod_int = @inferred modifying_integrator(series_integrator)

        for t in BicoloredRootedTreeIterator(2)
          @test iszero(mod_int[t])
        end
      end
    end
  end

  @testset "SSPRK33 two times" begin
    # Using the same method multiple times is equivalent to using a plain RK
    # method without any splitting/partitioning/decomposition
    A = @SArray [0 0 0; 1 0 0; 1//4 1//4 0]
    b = @SArray [1//6, 1//6, 2//3]
    rk = RungeKuttaMethod(A, b)
    ark = AdditiveRungeKuttaMethod([rk, rk])

    # third-order accurate
    series_integrator = @inferred bseries(ark, 3)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact))

    # not fourth-order accurate
    series_integrator = @inferred bseries(ark, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false

    # modified equations and modifying integrators
    let series_rk = bseries(rk, 5)
      series_ark = bseries(ark, order(series_rk))

      for colored_tree in keys(series_ark)
        tree = rootedtree(colored_tree.level_sequence)
        @test series_rk[tree] == series_ark[colored_tree]
      end

      @testset "modified_equation" begin
        mod_eq_rk  = modified_equation(series_rk)
        mod_eq_ark = modified_equation(series_ark)
        for colored_tree in keys(mod_eq_ark)
          tree = rootedtree(colored_tree.level_sequence)
          @test mod_eq_rk[tree] == mod_eq_ark[colored_tree]
        end
      end

      @testset "modifying_integrator" begin
        mod_int_rk  = modifying_integrator(series_rk)
        mod_int_ark = modifying_integrator(series_ark)
        for colored_tree in keys(mod_int_ark)
          tree = rootedtree(colored_tree.level_sequence)
          @test mod_int_rk[tree] == mod_int_ark[colored_tree]
        end
      end
    end
  end
end

@testset "multirate infinitesimal split methods interface" begin
  # NOTE: This is still considered experimental and might change at any time!

  @testset "KW43" begin
    # Oswald Knoth, Ralf Wolke
    # Implicit-explicit Runge-Kutta methods for computing atmospheric reactive flows
    # Applied Numerical Mathematics, Volume 28, (1998) Issue 2-4
    # https://doi.org/10.1016/S0168-9274(98)00051-8
    ns = 4
    A = zeros(Polynomial{Rational{Int}, :x}, ns+1, ns)
    D = zeros(Int, ns+1, ns)
    G = zeros(Rational{Int}, ns+1, ns)
    c = zeros(Rational{Int}, ns+1)
    dts = zeros(Rational{Int}, ns+1)

    A[2, 1] = 1 // 2
    A[3, 1] = -1 // 6
    A[3, 2] = 2 // 3
    A[4, 1] = 1 // 3
    A[4, 2] = -1 // 3
    A[4, 3] = 1
    A[5, 1] = 1 // 6
    A[5, 2] = 1 // 3
    A[5, 3] = 1 // 3
    A[5, 4] = 1 // 6
    A[5, :] = A[5, :] - A[4, :]
    A[4, :] = A[4, :] - A[3, :]
    A[3, :] = A[3, :] - A[2, :]

    c[2] = 1 // 2
    c[3] = 1 // 2
    c[4] = 1
    c[5] = 1

    D[2,1] = 1
    D[3,2] = 1
    D[4,3] = 1
    D[5,4] = 1

    dts[1] = 0
    dts[2] = 1 // 2
    dts[3] = 0
    dts[4] = 1 // 2
    dts[5] = 0

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # third-order accurate
    series_integrator = @inferred bseries(mis, 3)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact))

    # not fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(==, &, values(series_integrator), values(series_exact)) == false
  end

  @testset "MRI-GARK-ERK33" begin
    # Sandu, Adrian.
    # "A class of multirate infinitesimal GARK methods".
    # SIAM Journal on Numerical Analysis 57, no. 5 (2019): 2300-2327.
    # https://doi.org/10.1137/18M1205492
    ns = 3
    A = zeros(Polynomial{Float64, :x}, ns+1, ns)
    D = zeros(ns+1, ns)
    G = zeros(ns+1, ns)
    c = zeros(ns+1)
    dts = zeros(ns+1)

    delta = -0.5

    # Γ^i in Sandu (2019) is the i-th term
    A[2, 1] = Polynomial([1 / 3])
    A[3, 1] = Polynomial([(-6 * delta - 7) / 12, (2 * delta + 1) / 2])
    A[3, 2] = Polynomial([(6 * delta + 11) / 12, -(2 * delta + 1) / 2])
    A[4, 1] = Polynomial([0.0, 0.5])
    A[4, 2] = Polynomial([(6 * delta - 5) / 12, -(2 * delta + 1) / 2])
    A[4, 3] = Polynomial([(3 - 2 * delta) / 4, delta])

    # These should be the usual `c` coefficients → c[3] = 2 / 3?
    c[2] = 1 / 3
    c[3] = 2 / 3
    c[4] = 1

    # All one???
    D[2, 1] = 1
    D[3, 2] = 1
    D[4, 3] = 1

    # Differences of `c`s?
    dts[1] = 0
    dts[2] = 1 / 3
    dts[3] = 1 / 3
    dts[4] = 1 / 3

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # third-order accurate
    series_integrator = @inferred bseries(mis, 3)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact))

    # not fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false
  end

  @testset "MRI-GARK-ERK45a" begin
    # Sandu, Adrian.
    # "A class of multirate infinitesimal GARK methods".
    # SIAM Journal on Numerical Analysis 57, no. 5 (2019): 2300-2327.
    # https://doi.org/10.1137/18M1205492
    ns = 5
    A = zeros(Polynomial{Float64, :x}, ns+1, ns)
    D = zeros(ns+1, ns)
    G = zeros(ns+1, ns)
    c = zeros(ns+1)
    dts = zeros(ns+1)

    # Γ^i in Sandu (2019) is the i-th term
    A[2,1] = Polynomial([1 / 5])

    A[3,1] = Polynomial([-53 / 16, 503 / 80])
    A[3,2] = Polynomial([281 / 80, -503 / 80])

    A[4,1] = Polynomial([-36562993.0/71394880.0, -1365537.0/35697440.0])
    A[4,2] = Polynomial([34903117.0/17848720.0, 4963773.0/7139488.0])
    A[4,3] = Polynomial([-88770499.0/71394880.0, -1465833.0/2231090.0])

    A[5,1] = Polynomial([-7631593.0/71394880.0, 66974357.0/35697440.0])
    A[5,2] = Polynomial([-166232021.0/35697440.0, 21445367.0/7139488.0])
    A[5,3] = Polynomial([6068517.0/1519040.0, -3.0])
    A[5,4] = Polynomial([8644289.0/8924360.0, -8388609.0/4462180.0])

    A[6,1] = Polynomial([277061.0/303808.0, -18227.0/7520.0])
    A[6,2] = Polynomial([-209323.0/1139280.0, 2.0])
    A[6,3] = Polynomial([-1360217.0/1139280.0, 1.0])
    A[6,4] = Polynomial([-148789.0/56964.0, 5.0])
    A[6,5] = Polynomial([147889.0/45120.0, -41933.0/7520.0])

    # These should be the usual `c` coefficients
    c[2] = 1 / 5
    c[3] = 2 / 5
    c[4] = 3 / 5
    c[5] = 4 / 5
    c[6] = 1

    D[2,1] = 0
    D[3,2] = 1
    D[4,3] = 1
    D[5,4] = 1
    D[6,5] = 1

    # Differences of `c`s?
    dts[1] = 0.0
    dts[2] = 1 / 5
    dts[3] = 1 / 5
    dts[4] = 1 / 5
    dts[5] = 1 / 5
    dts[6] = 1 / 5

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact))

    # not fifth-order accurate
    series_integrator = @inferred bseries(mis, 5)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false
  end

  @testset "EB4" begin
    # Krogstad, Stein
    # "Generalized integrating factor methods for stiff PDEs".
    # Journal of Computational Physics 203, no. 1 (2005): 72-88.
    # https://doi.org/10.1016/j.jcp.2004.08.006
    ns = 4
    A = zeros(Polynomial{Float64, :x}, ns+1, ns)
    D = zeros(ns+1, ns)
    G = zeros(ns+1, ns)
    c = zeros(ns+1)
    dts = zeros(ns+1)

    fac1 = 1.0
    fac2 = 2.0

    A[2,1] = Polynomial([0.5])

    A[3,1] = Polynomial([0.5, -1 / fac1])
    A[3,2] = Polynomial([0.0, 1 / fac1])

    A[4,1] = Polynomial([1.0, -2 / fac1])
    A[4,3] = Polynomial([0.0, 2 / fac1])

    A[5,1] = Polynomial([1.0, -3 / fac1, 4 / fac2])
    A[5,2] = Polynomial([0.0, 2 / fac1, -4 / fac2])
    A[5,3] = Polynomial([0.0, 2 / fac1, -4 / fac2])
    A[5,4] = Polynomial([0.0, -1 / fac1, 4 / fac2])

    c[1] = 0
    c[2] = 0.5
    c[3] = 0.5
    c[4] = 1
    c[5] = 1

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact))

    # not fifth-order accurate
    series_integrator = @inferred bseries(mis, 5)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false
  end

  @testset "OwrenRK4" begin
    # Owren, Brynjulf.
    # "Order conditions for commutator-free Lie group methods".
    # Journal of Physics A: Mathematical and General 39, no. 19 (2006): 5585.
    # https://doi.org/10.1088/0305-4470/39/19/S15
    ns = 5
    A = zeros(Polynomial{Float64, :x}, ns+1, ns)
    D = zeros(ns+1, ns)
    G = zeros(ns+1, ns)
    c = zeros(ns+1)
    dts = zeros(ns+1)

    A[2,1] = Polynomial([1 / 2])

    A[3,2] = Polynomial([1 / 2])

    A[4,1] = Polynomial([-1 / 2])
    A[4,3] = Polynomial([1.0])

    A[5,1] = Polynomial([1 / 4])
    A[5,2] = Polynomial([1 / 6])
    A[5,3] = Polynomial([1 / 6])
    A[5,4] = Polynomial([-1 / 12])

    A[6,1] = Polynomial([-1 / 12])
    A[6,2] = Polynomial([1 / 6])
    A[6,3] = Polynomial([1 / 6])
    A[6,4] = Polynomial([1 / 4])

    D[4,2] = 1
    D[6,5] = 1

    dts[1] = 0
    dts[2] = 1 / 2
    dts[3] = 1 / 2
    dts[4] = 1 / 2
    dts[5] = 1 / 2
    dts[6] = 1 / 2

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact))

    # not fifth-order accurate
    series_integrator = @inferred bseries(mis, 5)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false
  end

  @testset "MISBI4_4" begin
    # Knoth, Oswald, and Joerg Wensch.
    # "Generalized split-explicit Runge–Kutta methods for the compressible Euler equations".
    # Monthly Weather Review 142, no. 5 (2014): 2067-2081.
    # https://doi.org/10.1175/MWR-D-13-00068.1
    # Table 4, method MIS4a
    ns = 4
    A = zeros(Polynomial{Float64, :x}, ns+1, ns)
    D = zeros(ns+1, ns)
    G = zeros(ns+1, ns)
    c = zeros(ns+1)
    dts = zeros(ns+1)

    # β in Knoth, Wensch (2014)
    A[2,1] = Polynomial([1 / 2])
    A[2,1] = Polynomial([0.38758444641450318])
    A[3,1] = Polynomial([-2.5318448354142823e-002])
    A[3,2] = Polynomial([0.38668943087310403])
    A[4,1] = Polynomial([0.20899983523553325])
    A[4,2] = Polynomial([-0.45856648476371231])
    A[4,3] = Polynomial([0.43423187573425748])
    A[5,1] = Polynomial([-0.10048822195663100])
    A[5,2] = Polynomial([-0.46186171956333327])
    A[5,3] = Polynomial([0.83045062122462809])
    A[5,4] = Polynomial([0.27014914900250392])

    # seems to be unused?
    c[2] = 0.38758444641450318
    c[3] = 0.61521685655017821
    c[4] = 0.23254717315441453
    c[5] = 1.0000000000000002

    # α in Knoth, Wensch (2014)
    D[3,2] =  0.52349249922385610
    D[4,2] =  1.1683374366893629
    D[4,3] = -0.75762080241712637
    D[5,2] = -3.6477233846797109e-002
    D[5,3] =  0.56936148730740477
    D[5,4] =  0.47746263002599681

    # γ in Knoth, Wensch (2014)
    G[3,2] =  0.13145089796226542
    G[4,2] = -0.36855857648747881
    G[4,3] =  0.33159232636600550
    G[5,2] = -6.5767130537473045e-002
    G[5,3] =  4.0591093109036858e-002
    G[5,4] =  6.4902111640806712e-002

    for i in 1:ns+1
      for j in 1:i-1
        dts[i] = dts[i] + A[i,j][0]
      end
    end

    mis = @inferred MultirateInfinitesimalSplitMethod(A, D, G, c)

    # third-order accurate
    series_integrator = @inferred bseries(mis, 3)
    series_exact = @inferred ExactSolution(series_integrator)
    @test_broken mapreduce(≈, &, values(series_integrator), values(series_exact))

    # not fourth-order accurate
    series_integrator = @inferred bseries(mis, 4)
    series_exact = @inferred ExactSolution(series_integrator)
    @test mapreduce(≈, &, values(series_integrator), values(series_exact)) == false
  end
end

end # @testset "BSeries"


text:
using JuMP
using FlexJuMP
using Gurobi

# Setup the uncertainty set parameters
means = [620; 388; 583; 313]
a = 5
covar = [11.11 a a a; a 11.11 a a; a a 11.11 a; a a a 11.11]
box_dev = ones(4) * 10

# Setup the model
m = FlexibilityModel(solver = GurobiSolver(OutputFlag = 0))

# Define variables
@randomvariable(m, T[i = 1:4], mean = means[i])
@recoursevariable(m, Qc)
@variable(m, x)

# Define the constraints
@constraint(m, -100 - 0.67Qc + 2T[2] + x <= 0.0)
@constraint(m, -250 - T[2] == x)
@constraint(m, 1388.5 + 0.5Qc - 0.75T[1] - T[2] - T[3] <= 0.0)
@constraint(m, 2044 + Qc - 1.5T[1] - 2T[2] - T[3] <= 0.0)
@constraint(m, 2830 + Qc - 1.5T[1] - 2T[2] - T[3] - 2T[4] <= 0.0)
@constraint(m, -3153 - Qc + 1.5T[1] + 2T[2] + T[3] + 3T[4] <= 0.0)

# Define the uncertainty set
setcovariance(m, covar)
# setuncertaintyset(m, :Ellipsoid)
setuncertaintyset(m, :Hyperbox, [[box_dev]; [box_dev]])
# setuncertaintyset(m, :PNorm, 1)
# setuncertaintyset(m, :PNorm, 2)
# setuncertaintyset(m, :PNorm, Inf)

# Ensure the mean is valid
if ismeanfeasible(m)

    # Can estimate the SF without first solving the flexibility model
    SF = findstochasticflexibility(m, num_pts = 2000, use_vulnerability_model = true)

    # Solve
    status = solve(m, active_constr = true)

    if status == :Optimal
        # Retrieve optimized data
        flexibility_index = getflexibilityindex(m)
        temperatures = getvalue(T)
        cooling = getvalue(Qc)
        state = getvalue(x)
        actives = getactiveconstraints(m)

        # Parse data
        data = getflexibilitydata(m)

        # Get the confidence level
        if data.uncertainty_set.name == :Ellipsoid
            conf_lvl = getconfidencelevel(m)
        end

        # Rank the inequality constraints
        rank_data = rankinequalities(m, max_ranks = 3, active_constr = true)
    end
end

ext: jl


text:
function loadCommunityFromCSV(abm::Model, file::String)

    #Global Array
    d = DataFrame(CSV.File(string(file,"_global.csv")))
    com = Community(abm,N=d.N[1]) #Create Community
    com.t = d.t[1]
    for i in com.declaredSymbols_["Global"]
        setproperty!(com, i, d[!,i][1])
    end
    for i in com.declaredSymbols_["GlobalArray"]
        values = sort([j for j in names(d) if occursin(string(i,"__"),j)])
        indexes = [Meta.parse.(split(split(j,"__")[2],"_")) for j in values]
        m = zeros(indexes[end]...)
        for (v,ind) in zip(values,indexes)
            m[ind...] = d[1,v]
        end
        setproperty!(com, i, m)
    end

    #Local and Identity load
    d = DataFrame(CSV.File(string(file,"_local.csv")))
    for i in com.declaredSymbols_["Local"]
        setproperty!(com, i, d[!,i])
    end
    for i in com.declaredSymbols_["Identity"]
        setproperty!(com, i, Int.(d[!,i]))
    end

    return com
end

function loadCommunityInTimeFromCSV(abm::Model, file)

    out = false
    count = 1

    com = CommunityInTime()
    while !out
        fileName = string(file,"_",count,".csv")
        if isfile(fileName)
            push!(com,loadCommunityFromCSV(abm,fileName))
            count += 1
        else
            out = true
        end
    end

    return com
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
struct Report
	text::Any
	values::Dict
	table::Any
end

function Report(; text="Nothing to report", values=Dict(), table=Dict())
    Report(text, values, table)
end

function Base.show(io::IO, report::Report)
	println(io, report.text)
end

function table(report::Report)
	return report.table
end

function Base.values(report::Report)
	return report.values
end

function format_text(text)
	text = replace(text, " -0.0 " => " 0.0 ")
	text = replace(text, " 0.0 " => " 0 ")
	text = replace(text, " -0.0," => " 0.0,")
	text = replace(text, " 0.0," => " 0,")
	text = replace(text, " -0.0)" => " 0.0)")
	text = replace(text, " 0.0)" => " 0)")
end


text:
module PSSR
# External package
using Revise
using ProgressMeter: showprogress
using ProgressMeter
using Flux
using Flux: train!
using Flux: throttle
using IterTools
using Statistics: mean

# Program packages
export get_protein_data
include("./get_protein_data.jl")

export check_data
export count_size
include("./test.jl")

export prepare_data
export build_matrix_X
export build_matrix_Y
include("./prepare_data.jl")

export nettrain
export accuracy
export onecoldencoding
include("./Neural_network.jl")

export write_dataset
include("./write_dataset.jl")
# Write your package code here.

end


text:
__precompile__()

module GradientNoise

    include("perlin.jl")
    include("opensimplex.jl")

    export perlin,
        octaveperlin
        smoothstep,
        smootherstep
end


text:
"""
This package provides Julia bindings for the
[extxyz](https://github.com/libAtoms/extxyz) C library which implements a parser
and writer for the extended XYZ file format used in materials and molecular
modelling, following the
[specification](https://github.com/libAtoms/extxyz#extended-xyz-specification-and-parsing-tools)
set out in the extxyz repo.

## Basic Usage

Four key functions are exported: `read_frame()` and `write_frame()` for reading
and writing single configurations (snapshots), respectively, and `read_frames()`
and `write_frames()` for reading and writing trajectories. All functions can
work with string filenames, an open `Base.IO` instance or (intended primarily
for internal use) a C `FILE*` pointer, stored as a `Ptr{Cvoid}` type.

```julia
using ExtXYZ

frame = read_frame("input.xyz")  # single atomic configuration, represented as a Dict{String}{Any}
write_frame("output.xyz", frame) # write a single frame to `output.xyz`.

frame10 = read_frame("input.xyz", 10) # read a specific frame, counting from 1 for first frame in file

all_frames = read_frames("seq.xyz")  # read all frames, returns Vector{Dict{String}{Any}}
frames = read_frames("seq.xyz", 1:4) # specific range of frames

write_frames("output.xyz", frames, append=true) # append four frames to output
```
"""
module ExtXYZ

using extxyz_jll

cfopen(filename::String, mode::String) = ccall(:fopen,
                                               Ptr{Cvoid},
                                               (Cstring, Cstring),
                                               filename, mode)

cfclose(fp::Ptr{Cvoid}) = ccall(:fclose,
                                Cint,
                                (Ptr{Cvoid},),
                                fp)

function cfopen(f::Function, iostream::IOStream, mode::String="r")
    newfd = Libc.dup(RawFD(fd(iostream)))
    fp = ccall(@static(Sys.iswindows() ? :_fdopen : :fdopen),
               Ptr{Cvoid}, (Cint, Cstring), newfd, mode)
    try
        f(fp)
    finally
        cfclose(fp)
    end
end

function cfopen(f::Function, filename::String, mode::String="r")
    fp = cfopen(filename, mode)
    try
        f(fp)
    finally
        cfclose(fp)
    end
end

struct DictEntry
    key::Ptr{Cchar}
    data::Ptr{Cvoid}
    data_t::Cint
    nrows::Cint
    ncols::Cint
    next::Ptr{DictEntry}
    first_data_ll::Ptr{Cvoid}
    last_data_ll::Ptr{Cvoid}
    n_in_row::Cint
end

const _kv_grammar = Ref{Ptr{Cvoid}}(0)

function __init__()
    _kv_grammar[] = ccall((:compile_extxyz_kv_grammar, libextxyz),
                           Ptr{Cvoid}, ())
    nothing
end

cfree_dict(dict::Ptr{DictEntry}) = ccall((:free_dict, libextxyz),
                                        Cvoid,
                                        (Ptr{DictEntry},),
                                        dict)

cprint_dict(dict::Ptr{DictEntry}) = ccall((:print_dict, libextxyz),
                                        Cvoid,
                                        (Ptr{DictEntry},),
                                        dict)

const DATA_I = 1
const DATA_F = 2
const DATA_B = 3
const DATA_S = 4

const TYPE_MAP = Dict(DATA_I => Cint,
                      DATA_F => Cdouble,
                      DATA_B => Cint,
                      DATA_S => Cstring)

import Base: convert

function convert(::Type{Dict{String,Any}}, c_dict::Ptr{DictEntry}; transpose_arrays=false)
    result = Dict{String,Any}()
    node_ptr = c_dict
    while node_ptr != C_NULL
        node = unsafe_load(node_ptr)
        data_ptr = reinterpret(Ptr{TYPE_MAP[node.data_t]}, node.data)

        if node.nrows == 0 && node.ncols == 0
            # scalar
            value = unsafe_load(data_ptr)
            # convert to primitive types
            if node.data_t == DATA_S
                value = unsafe_string(value)
            elseif node.data_t == DATA_B
                value = convert(Bool, value)
            end
        else
            # array, either 1D or 2D
            if node.nrows == 0
                # vector (1D array)
                dims = (node.ncols, )
            else
                # matrix (2D array)
                dims = (node.nrows, node.ncols)
            end

            value = unsafe_wrap(Array,
                                reinterpret(Ptr{TYPE_MAP[node.data_t]}, node.data),
                                dims)

            if node.data_t == DATA_S
                value = unsafe_string.(value)
            elseif node.data_t == DATA_B
                value = convert(Array{Bool}, value)
            else
                value = copy(value)
            end

            if node.nrows != 0 && node.ncols != 0
                value = reshape(value, size(value, 2), size(value, 1))
                transpose_arrays && (value = permutedims(value, (2, 1)))
            end
        end

        key = unsafe_string(node.key)
        result[key] = value
        node_ptr = node.next
    end
    return result
end

# utility functions for converting from Julia types to corresponding C type and value

Ctype(::Type{S}) where S <: AbstractArray{T,N} where {T,N} = Ctype(T)
Ctype(::Type{<:Integer}) = (Cint, DATA_I)
Ctype(::Type{Bool}) = (Cint, DATA_B)
Ctype(::Type{<:Real}) = (Cdouble, DATA_F)
Ctype(::Type{String}) = (Cstring, DATA_S)

Cvalue(value::T; transpose_arrays=nothing) where {T<:Union{Integer,Bool,Real}} = convert(Ctype(T)[1], value)

function Cvalue(value::String; transpose_arrays=nothing)
    ptr = pointer(Base.unsafe_convert(Cstring, Base.cconvert(Cstring, value)))
    new = Ptr{Cchar}(Libc.malloc(sizeof(value)+1))
    unsafe_copyto!(new, ptr, sizeof(value))
    unsafe_store!(new, C_NULL, sizeof(value)+1)
    return new
end

function Cvalue(value::Array{T,N}; transpose_arrays=false) where {T<:Union{Integer,Bool,Real},N}
    if ndims(value) == 2
        value = reshape(value, size(value, 2), size(value, 1))
    end
    value = convert(Array{Ctype(T)[1]}, value)
    if ndims(value) == 2 && transpose_arrays
        value = permutedims(value, (2, 1))
    end
    return value
end

function Cvalue(value::Array{String,1}; transpose_arrays=nothing)
    result = Array{Ptr{Cchar}}(undef, length(value))
    result .= Cvalue.(value)
    return result
end

# dims(value) -> (nrows, ncols) for scalar, vector and matrix value
dims(value) = (0, 0)
dims(value::AbstractVector) = (0, size(value, 1))
dims(value::AbstractMatrix) = (size(value, 1), size(value, 2))

function convert(::Type{Ptr{DictEntry}}, dict::Dict{String}{Any}; ordered_keys=nothing, transpose_arrays=false)
    c_dict_ptr = Ptr{DictEntry}(Libc.malloc(sizeof(DictEntry)))
    node_ptr = c_dict_ptr

    ordered_keys === nothing && (ordered_keys = keys(dict))
    for (idx, key) in enumerate(ordered_keys)
        value = dict[key]
        ckey = Cvalue(key)
        cvalue = Cvalue(value; transpose_arrays=transpose_arrays)
        nrow, ncol = dims(cvalue)
        type, data_t = Ctype(typeof(value))
        data = Ptr{type}(Libc.malloc(sizeof(cvalue)))
        if nrow == 0 && ncol == 0
            unsafe_store!(data, cvalue)
        else
            ptr = pointer(cvalue)
            data_t == DATA_S && (ptr = reinterpret(Ptr{type}, ptr)) # char** -> char*
            unsafe_copyto!(data, ptr, length(cvalue))
        end

        # allocate another DictEntry struct unless we're on the last one already
        next_ptr = C_NULL
        idx != length(dict) && (next_ptr = Ptr{DictEntry}(Libc.malloc(sizeof(DictEntry))))

        node = DictEntry(ckey, data, data_t, nrow, ncol,
                         next_ptr, C_NULL, C_NULL, 0)
        unsafe_store!(node_ptr, node) # in place mutation of node_ptr
        node_ptr = next_ptr
    end
    return c_dict_ptr
end

function read_frame_dicts(fp::Ptr{Cvoid}; verbose=false)
    nat = Ref{Cint}(0)
    info = Ref{Ptr{DictEntry}}()
    arrays = Ref{Ptr{DictEntry}}()
    eof = false
    try
        res =  ccall((:extxyz_read_ll, libextxyz),
                      Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Cint}, Ptr{Ptr{DictEntry}}, Ptr{Ptr{DictEntry}}),
                      _kv_grammar[], fp, nat, info, arrays)
        if res != 1
            eof = true
            throw(EOFError())
        end

        if verbose
            cprint_dict(info[])
            cprint_dict(arrays[])
        end

        pinfo = reinterpret(Ptr{DictEntry}, info[])
        parrays = reinterpret(Ptr{DictEntry}, arrays[])
        jinfo = convert(Dict{String,Any}, pinfo, transpose_arrays=true)
        jarrays = convert(Dict{String,Any}, parrays, transpose_arrays=false)
        return nat[], jinfo, jarrays

    finally
        if !eof
            cfree_dict(info[])
            cfree_dict(arrays[])
        end
    end
end

"""
extract "Lattice" entry and apply semantic conversions
"""
function extract_lattice!(result_dict)
    "Lattice" in keys(result_dict) || return nothing
    lattice = pop!(result_dict, "Lattice")
    if size(lattice) == (3, 3)
        lattice = convert(Array{Float64}, lattice)
    elseif size(lattice) == (3,)
        lattice = convert(Array{Float64}, diagm(lattice))
    elseif lattice.shape == (9,)
        lattice = convert(Array{Float64}, reshape(lattice, 3, 3))
    else
        error("Lattice has wrong shape!")
    end
    return lattice
end

"""
    read_frame(file)

Read a single frame from the ExtXYZ file `file`, which can be a file pointer,
open IO stream or a string filename.
"""
function read_frame(fp::Ptr{Cvoid}; verbose=false)
    nat, info, arrays = try
        read_frame_dicts(fp; verbose=verbose)
    catch err
        if isa(err, EOFError)
            return nothing
        end
        rethrow()
    end

    dict = Dict{String, Any}()
    dict["N_atoms"] = nat # number of atoms

    # periodic boundary conditions
    if "pbc" in keys(info)
        dict["pbc"] = pop!(info, "pbc")
    end

    # cell is transpose of the stored lattice
    lattice = extract_lattice!(info)
    dict["cell"] = permutedims(lattice, (2, 1))

    delete!(info, "Properties")

    # everything else stays in info and arrays
    dict["info"] = info
    dict["arrays"] = arrays

    return dict
end

read_frame(file::Union{String,IOStream}, index; kwargs...) = read_frames(file, index; kwargs...)[1]
read_frame(file::Union{String,IOStream}; kwargs...) = read_frame(file, 1; kwargs...)

"""
    iread_frames(file[, range])

Return a Channel for reading from an ExtXYZ file. Frames are yielded one by one.
`range` can be a single integer, range object or integer array of frame indices.

Example usage:

```julia
ch = iread_frames("file.xyz")
for frame in ch
    process(frame)
end
```
"""
function iread_frames(fp::Ptr{Cvoid}, range; close_fp=false, kwargs...)
    Channel() do channel
        for frame in 1:first(range)-1
            atoms = read_frame(fp; kwargs...)
            atoms === nothing && break
        end
        for frame in range
            atoms = read_frame(fp; kwargs...)
            atoms === nothing && break
            put!(channel, atoms)
        end
        if (close_fp)
            cfclose(fp)
        end
        channel
    end
end

function iread_frames(file::Union{String,IOStream}, range; kwargs...)
    fp = cfopen(file, "r")
    fp == C_NULL && error("file $file cannot be opened for reading")
    iread_frames(fp, range; close_fp=true, kwargs...)
end

iread_frames(file::Union{String,IOStream}, index::Int; kwargs...) = iread_frames(file, [index]; kwargs...)
iread_frames(file::Union{String,IOStream}; kwargs...) = iread_frames(file, Iterators.countfrom(1); kwargs...)

"""
    read_frames(file[, range])

Read a sequence of frames from the ExtXYZ `file`, which can be specified by a file pointer, filename or IOStream.

`range` can be a single integer, range object or integer array of frame indices.
"""
read_frames(fp::Ptr{Cvoid}, range; kwargs...) = collect(iread_frames(fp, range; kwargs...))

function read_frames(file::Union{String,IOStream}, range; kwargs...)
    cfopen(file) do fp
        fp == C_NULL && error("file $file cannot be opened for reading")
        read_frames(fp, range; kwargs...)
    end
end

read_frames(file::Union{String,IOStream}, index::Int; kwargs...) = read_frames(file, [index]; kwargs...)
read_frames(file::Union{String,IOStream}; kwargs...) = read_frames(file, Iterators.countfrom(1); kwargs...)

function write_frame_dicts(fp::Ptr{Cvoid}, nat, info, arrays; verbose=false)
    nat = Cint(nat)
    cinfo = convert(Ptr{DictEntry}, info; transpose_arrays=true)

    # ensure "species" goes in column 1 and "pos" goes in column 2
    ordered_keys = collect(keys(arrays))
    species_idx = findfirst(isequal("species"), ordered_keys)
    ordered_keys[1], ordered_keys[species_idx] = ordered_keys[species_idx], ordered_keys[1]
    pos_idx = findfirst(isequal("pos"), ordered_keys)
    ordered_keys[2], ordered_keys[pos_idx] = ordered_keys[pos_idx], ordered_keys[2]
    carrays = convert(Ptr{DictEntry}, arrays; ordered_keys=ordered_keys)

    if verbose
        cprint_dict(cinfo)
        cprint_dict(carrays)
    end
    try
        res = ccall((:extxyz_write_ll, libextxyz),
                     Cint, (Ptr{Cvoid}, Cint, Ptr{DictEntry}, Ptr{DictEntry}), fp, nat, cinfo, carrays)
        res != 0 && error("error writing to file")
    finally
        cfree_dict(cinfo)
        cfree_dict(carrays)
    end
end

"""
    write_frame(file, dict)

Write a single atomic configuration represented by `dict` to `file`, which
can be a file pointer, open IO stream or string filename.
"""
function write_frame(fp::Ptr{Cvoid}, dict; verbose=false)
    nat = dict["N_atoms"]
    info = copy(dict["info"])
    info["Lattice"] = permutedims(dict["cell"], (2, 1))
    info["pbc"] = get(dict, "pbc", [true, true, true])

    write_frame_dicts(fp, nat, info, dict["arrays"]; verbose=verbose)
end

"""
    write_frames(file, dicts)

Write a sequence of atomic configurations to `file`. Can also be used asynchronously
by passing a Channel in place of `dicts`, e.g.

```julia
Channel() do ch
    @async write_frames(outfile, ch)

    for frame in frames
        put!(ch, frame)
    end
end
```
"""
write_frames(fp::Ptr{Cvoid}, dicts; kwargs...) = write_frame.(dicts)

function write_frames(file::Union{String,IOStream}, dicts; append=false, kwargs...)
    mode = append ? "a" : "w"
    cfopen(file, mode) do fp
        fp == C_NULL && error("file $file cannot be opened for writing")
        for dict in dicts
            write_frame(fp, dict; kwargs...)
        end
    end
end

write_frame(file::Union{String,IOStream}, dict; kwargs...) = write_frames(file, [dict]; kwargs...)

export read_frame, read_frames, iread_frames, write_frame, write_frames

end


text:
module Input

#=
This function reads the data from a txt or csv file.
:param infile: the path of the file
:return: an array with the read data
=#
function read_txt(infile::String)
    a=read(open(infile,"r"),String)
    return parse.(Float64,filter!(e->e!="",split(a,r"[^\d.]")))
end # infile

#=
This function reads the data from a wfdb file.
:param record: the record name
:param annotator: the annotator of the record
:return: an array with the read data
=#
function read_wfdb(record::String,annotator::String)
    temp=string(record,"_temp.txt")
    run(pipeline(`ann2rr -r "$record" -a "$annotator" -i s -c`,stdout="$temp"))
    a=read_txt("$temp")
    run(`rm "$temp"`)
    return a*1000
end # infile

end # module


text:
##############################   Born-Mayer  ###################################
# mutable struct BornMayer <: EmpiricalPotential
#     A::Float64
#     ρ::Float64
# end

# function BornMayer()
#     #ToDO
#     return BornMayer(1.0, 1.0)
# end

# function get_trainable_params(bm::BornMayer)
#     return (A = bm.A, ρ = bm.ρ)
# end

# function get_nontrainable_params(lj::BornMayer)
#     p = Parameter{}(())
#     return p
# end

# ##############################   Energy  ###################################

# function potential_energy(r::Vector{<:Real}, p::BornMayer)
#     return p.A * exp(-norm(r) / p.ρ)
# end

# ##############################   Force   ###################################

# function force(r::Vector{<:Real}, p::BornMayer)
#     d = norm(r)
#     return p.A /p.ρ * exp(-d / p.ρ ) .* r ./ d
# end

# ##############################   Gradients  ###################################

# function grad_potential_energy(r::Vector{<:Real}, p::BornMayer)
#     d = norm(r)
#     return (dpdA = exp(-d / p.ρ),
#             dpdρ = p.A * d * exp(-d/p.ρ) / p.ρ^2 )
# end

# function grad_force(r::Vector{<:Real}, p::BornMayer)
#     d = norm(r)
#     return (dfdA = 1.0 /p.ρ * exp(-d / p.ρ ) .* r ./ d,
#             dfdρ = p.A / p.ρ^3 * exp(-d / p.ρ )*(d - p.ρ)  .* r ./ d)
# end

# function grad_virial(r::Vector{<:Real}, p::BornMayer)
#     df = grad_force(r, p)
#     dfdA = df[:dfdA]
#     dfdrho = df[:dfdρ]
#     return (dvdA = dfdA[1]*r[1] + dfdA[2]*r[2] + dfdA[3]*r[3],
#             dvdρ = dfdrho[1]*r[1] + dfdrho[2]*r[2] + dfdrho[3]*r[3])
# end


text:
#
# Part of the Manager module for working with backups. It may
# save and load application state into/from the file.
#
# TODO: describe index based saving format
#
# @author DeadbraiN
#
import Backup
import Creature
import Server
import Client
import Creature
import ManagerTypes
import Helper

export recover
export backup
#
# This function is used for recovering a manager's data from
# backup file. It means that an application was crashed before
# and now we have to recover it with last correct backup. Works
# in pair with backup() function.
# @param man Manager data type
# @return {Bool} recover status
#
function recover(man::ManagerTypes.ManagerData)
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info(string("Recovering from backup: ", Backup.lastFile())) end
  local data = Backup.load()
  local task::ManagerTypes.OrganismTask
  #
  # In case if we haven't recovered from backup file,
  # we have to remove it, because it's broken
  #
  if data === nothing
    if !_removeNew(man) return false end
    return Backup.lastFile() !== "" ? recover(man) : false
  end
  #
  # Main task reference should be set before the loop, because
  # it may be used inside organism tasks
  #
  man.task = current_task()
  for task in man.tasks
    task.organism.codeFn = Creature.eval(task.organism.code)
    #
    # This line prevents observers to be shared between organisms
    # Every organism should have it's own observer's instance
    #
    task.organism.observer = Event.create()
    bindEvents(man, task.organism)
    #
    # We must call yieldto(), because born() function parameters
    # will not be passed into the task and only last organism will
    # be run.
    #
    Manager._updateOrgTask(man, task)
  end
  #
  # We don't need to load\create servers and clients every
  # time on load backup. They should be created every time
  # on Manager application starts. So we have to skip loading
  # of man.cons property.
  #
  man.cfg            = data.cfg
  man.world          = data.world
  man.positions      = data.positions
  man.organisms      = data.organisms
  man.tasks          = data.tasks
  man.params         = data.params
  man.organismId     = data.organismId
  man.totalOrganisms = data.totalOrganisms
  man.maxEnergy      = data.maxEnergy
  man.plugins        = data.plugins
  man.killed         = data.killed
  #
  # We don't need to store amount of eval() calls
  # between app runs
  #
  data.cfg.orgEvals = 0
  #
  # We have to remove all event handlers from observers
  # after backup loading, because they may be multiplied
  # on every app running.
  #
  man.world.obs = Event.create()

  true
end
#
# Makes a dump of Manager data and saves it into the backup file.
# Works in pair with recover().
# @return {Bool} Backup status
#
function backup(man::ManagerTypes.ManagerData)
  local manCopy::ManagerTypes.ManagerData = deepcopy(man)
  local tmpTask::Task = Task(Helper.emptyFn)
  local task::ManagerTypes.OrganismTask
  local ret::Bool
  local org::Creature.Organism
  local posPair::Pair{Int, Creature.Organism}
  local orgPair::Pair{UInt, Creature.Organism}
  #
  # We have to stop the task before it will be saved into the backup file.
  # yield() call is needed, because Julia has strange issue with yieldto()
  # or i don't understand it's logic. Anyway, it stucks without this call.
  # TODO: do we need this?
  # TODO: move this yieldto() and fire() in separate function
  # TODO: reuse this approach in every file...
  yield()
  Event.fire(man.obs, "yield", man)
  consume(tmpTask)
  #
  # This is a small trick. We have to set all tasks in "done"
  # state for serializing into the file. Julia can't save active
  # tasks. After recover we will rerun them. It also can't save
  # anonymous functions, so we have to have some empty function
  # as a part of current Manager module (Helper.emptyFn()).
  #
  for task in manCopy.tasks
    org          = task.organism
    task.task    = tmpTask
    org.codeFn   = Helper.emptyFn
    org.observer = Event.create()
  end
  for posPair in manCopy.positions
    org          = posPair[2]
    org.codeFn   = Helper.emptyFn
    org.observer = Event.create()
  end
  for orgPair in manCopy.organisms
    org          = orgPair[2]
    org.codeFn   = Helper.emptyFn
    org.observer = Event.create()
  end
  manCopy.world.obs    = Event.create()
  manCopy.dotCallback  = Helper.emptyFn
  manCopy.moveCallback = Helper.emptyFn
  manCopy.task         = tmpTask
  manCopy.cons         = ManagerTypes.Connections()
  #
  # Before saving backup file, we need to make sure that backup folder is in place.
  # Otherwise - create one.
  #
  Backup.folder()
  #
  # These code lines create new backup file and remove last one
  #
  if (ret = Backup.save(manCopy))
    _removeOld(man)
    if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info(string("Backup has created: ", Backup.lastFile())) end
  end
  Event.fire(man.obs, "backup", man)

  ret
end
#
# Removes first last backup by date. Only last backupAmount files are
# available at the moment.
# @param man Reference to manager data object
#
function _removeOld(man::ManagerTypes.ManagerData)
  local files::Array{String, 1} = Backup.getFiles(Backup.FOLDER_NAME)
  local len::Int = length(files)
  local i::Int

  if len <= man.cfg.backupAmount return true end
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Removing old backup files...") end
  for i = 1:(len - man.cfg.backupAmount)
    rm(Backup.FOLDER_NAME * "/" * files[i])
  end

  true
end
#
# Removes first new backup by date. It's needed for removing broken
# backup files.
# @param man Reference to manager data object
#
function _removeNew(man::ManagerTypes.ManagerData)
  local files::Array{String, 1} = Backup.getFiles(Backup.FOLDER_NAME)
  local len::Int = length(files)
  local file::String = len > 0 ? files[len] : ""

  if len < 1 return false end
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info(string("Backup has removed: ", file)) end
  rm(Backup.FOLDER_NAME * "/" * file)

  true
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
    full_introspection!(client::Client)

Performs a full instrospection of the GraphQL schema of the `client`. The results are stored
in the client struct.

See also: [`Client`](@ref)
"""
function full_introspection!(client::Client)
    query_str = """
    query IntrospectionQuery {
        __schema {
            queryType {
                name
            }
            mutationType {
                name
            }
            subscriptionType {
                name
            }
            types {
                ...FullType
            }
            directives {
                name
                description
                locations
                args {
                    ...InputValue
                }
            }
        }
    }

    fragment FullType on __Type {
        kind
        name
        description
        fields(includeDeprecated: true) {
            name
            description
            args {
                ...InputValue
            }
            type {
                ...TypeRef
            }
            isDeprecated
            deprecationReason
        }
        inputFields {
            ...InputValue
        }
        interfaces {
            ...TypeRef
        }
        enumValues(includeDeprecated: true) {
            name
            description
            isDeprecated
            deprecationReason
        }
        possibleTypes {
            ...TypeRef
        }
    }
    fragment InputValue on __InputValue {
        name
        description
        type {
            ...TypeRef
        }
        defaultValue
    }
    fragment TypeRef on __Type {
        kind
        name
        ofType {
            kind
            name
            ofType {
                kind
                name
                ofType {
                    kind
                    name
                    ofType {
                        kind
                        name
                        ofType {
                            kind
                            name
                            ofType {
                                kind
                                name
                                ofType {
                                    kind
                                    name
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    """

    # Execute introspection query and extract schema
    body = execute(client, query_str)
    schema = body.data["__schema"]

    keys_to_include = ["queryType"]
    query_type = schema["queryType"]["name"]
    server_has_mutations = !isnothing(schema["mutationType"])
    server_has_subs = !isnothing(schema["subscriptionType"])
    if server_has_mutations
        mutation_type = schema["mutationType"]["name"]
        push!(keys_to_include, "mutationType")
    end
    if server_has_subs
        subscription_type = schema["subscriptionType"]["name"]
        push!(keys_to_include, "subscriptionType")
    end
    gql_types = [schema[k]["name"] for k in keys_to_include]

    # Populate input_object_fields_to_type_map
    empty!(client.input_object_fields_to_type_map) # Reset
    for gql_type in schema["types"]
        if gql_type["kind"] == "INPUT_OBJECT"
            field_to_type_map = Dict()
            for field in gql_type["inputFields"]
                field_to_type_map[field["name"]] = get_field_type_string(field)
            end
            client.input_object_fields_to_type_map[gql_type["name"]] = field_to_type_map
        end
    end

    # Populate query_to_type_map
    empty!(client.query_to_type_map) # Reset
    for gql_type in schema["types"]
        if gql_type["name"] in gql_types
            for field in gql_type["fields"]
                field_type = _recursive_get_value(
                    field["type"],
                    "ofType"
                )["name"]
                client.query_to_type_map[field["name"]] = field_type
                gql_type["name"] == query_type && push!(client.queries, field["name"])
                server_has_mutations && gql_type["name"] == mutation_type && push!(client.mutations, field["name"])
                server_has_subs && gql_type["name"] == subscription_type && push!(client.subscriptions, field["name"])
            end
        end
    end

    # Populate query_to_args_map
    empty!(client.query_to_args_map) # Reset
    for gql_type in schema["types"]
        if gql_type["name"] in gql_types
            for field in gql_type["fields"]
                for arg in field["args"]
                    get!(client.query_to_args_map, field["name"], Dict{String, String}()) # Add if it doesn't exist
                    client.query_to_args_map[field["name"]][arg["name"]] = get_field_type_string(arg)
                end
            end
        end
    end

    # Populate type_to_fields_map
    empty!(client.type_to_fields_map) # Reset
    for gql_type in schema["types"]
        if get(gql_type, "fields", nothing) !== nothing
            client.type_to_fields_map[gql_type["name"]] = Dict(f["name"] => f for f in gql_type["fields"])
        elseif get(gql_type, "inputFields", nothing) !== nothing
            client.type_to_fields_map[gql_type["name"]] = Dict(f["name"] => f for f in gql_type["inputFields"])# [f["name"] for f in gql_type["inputFields"]]
        end
    end

    client.schema = schema
    client.introspection_complete = true

    return client
end

"""
    introspect_node(node)

Introspect single node in GQL schema.
"""
function introspect_node(client::Client, node)
    query_str = "{
            __type(name:\"$node\") {
                name
                fields {
                    description
                    type {
                        name
                        kind
                        ofType {
                            name
                        }
                    }
                }
            }
        }"
    body = execute(client, query_str)
    fields = [field["description"] for field in body["data"]["__type"]["fields"]]
    return fields
end

"""
    _recursive_get_value(dict, key)

Recursively iterate through a dict until value of key is nothing.
"""
function _recursive_get_value(dict, key)
    if dict[key] !== nothing
        return _recursive_get_value(dict[key], key)
    else
        return dict
    end
end


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

rpsl = "redcsl"

ρ = (!Sys.iswindows() && occursin("32",read(`uname -m`,String))) ?  3 : 4

date = Dict(
  0 =>  "2017-05-16",
  1 =>  "2018-01-17",
  2 =>  "2018-04-06",
  3 =>  "2019-04-13",
  4 =>  "2019-09-15")

rsvn = Dict(
  0 =>  "4052",
  1 =>  "4372",
  2 =>  "4567",
  3 =>  "4961",
  4 =>  "5129")

function _spawn(cmd, input=devnull, output=devnull)
    run(pipeline(cmd, stdin=input, stdout=output, stderr=stderr), wait=false)
end


text:
mutable struct State{T, H, P}
        Qu :: T
        Qv :: T
        Qθ :: T
        Qs :: T
        Qb :: T
         h :: T
    h_crit :: H
     plume :: P
end

h_criterion(args...) = nothing

function State(diffusivity, nonlocalflux, mixingdepth, grid, T=Float64)
    plume = instantiate_plume(nonlocalflux, grid)
    h_crit = h_criterion(mixingdepth, grid)
    return State(zero(T), zero(T), zero(T), zero(T), zero(T), zero(T),
                 h_crit, plume)
end

"""
    update_state!(model)

Update the top flux conditions and mixing depth for `model`
and store in `model.state`.
"""
function update_state!(m)

    m.state.Qu = getbc(m, m.bcs.U.top)
    m.state.Qv = getbc(m, m.bcs.V.top)
    m.state.Qθ = getbc(m, m.bcs.T.top)
    m.state.Qs = getbc(m, m.bcs.S.top)
    m.state.Qb = m.constants.g * (m.constants.α * m.state.Qθ - m.constants.β * m.state.Qs)

    update_mixing_depth!(m)

    update_nonlocal_flux!(m)

    return nothing
end


text:
module Indexes
using ..BigArrays

export colon2unit_range, chunkid2global_range, index2chunkid, index2unit_range
export global_range2buffer_range, global_range2chunk_range
export cartesian_range2unit_range, unit_range2string, cartesian_range2string

# make Array accept cartetian range as index
function cartesian_range2unit_range(r::CartesianIndices{N}) where N
    map((x,y)->x:y, first(r).I, last(r).I)
end

"""
    global_range2buffer_range(globalRange::CartesianRange, bufferGlobalRange::CartesianRange)

Transform a global range to a range inside buffer.
"""
function global_range2buffer_range(globalRange::CartesianIndices{N},
                                 bufferGlobalRange::CartesianIndices{N}) where N
    start = first(globalRange) - first(bufferGlobalRange) + one(first(globalRange))
    stop  = last(globalRange)  - first(bufferGlobalRange) + one(last(globalRange))
    return CartesianIndices( map((x,y)->x:y, start.I, stop.I) )
end

"""
    global_range2chunk_range(globalRange::CartesianRange, bufferGlobalRange::CartesianRange)

Transform a global range to a range inside chunk.
"""
function global_range2chunk_range(globalRange::CartesianIndices{N},
                                 chunkSize::NTuple{N};
                                 offset::CartesianIndex{N} = CartesianIndex{N}()-1) where N
    chunkID = index2chunkid(first(globalRange), chunkSize; offset=offset)
    start = index2cartesian_index( map((s,i,sz,o)->s-(i-1)*sz-o, first(globalRange).I,
                                                            chunkID, chunkSize, offset.I))
    stop  = index2cartesian_index( map((s,i,sz,o)->s-(i-1)*sz-o, last(globalRange).I,
                                                            chunkID, chunkSize, offset.I))
    return CartesianIndices(map((x,y)->x:y, start.I, stop.I))
end

function index2chunkid(idx::CartesianIndex{N}, chunkSize::NTuple{N};
                       offset::CartesianIndex{N} = CartesianIndex{N}()-1) where N
    # the offset is actually start of the real data, it could be not aligned to 0
    ( map((x,y,o)->fld(x-1-o, y)+1, idx.I, chunkSize, offset.I) ... ,)
end

function chunkid2global_range(chunkID::CartesianIndex{N}, chunkSize::NTuple{N};
                              offset::CartesianIndex{N} = CartesianIndex{N}()-1) where N
    start = index2cartesian_index( map((x,y)->(x-1)*y+1, chunkID.I, chunkSize) ) + offset
    stop  = index2cartesian_index( map((x,y)->x*y,       chunkID.I, chunkSize) ) + offset
    range = map((x,y)->x:y, start.I, stop.I)
    return CartesianIndices(range)
end

"""
replace Colon of indexes by UnitRange
"""
function colon2unit_range(buf::AbstractArray, indexes::Tuple)
    colon2unit_range(size(buf), indexes)
end

function colon2unit_range(sz::NTuple{N}, indexes::Tuple) where N
    map((x,y)-> x==Colon() ? UnitRange(1:y) : x, indexes, sz)
end

function index2cartesian_index( idx::Union{Tuple, Vector} )
    CartesianIndex((idx...,))
end

# function Base.string{N}(r::CartesianRange{CartesianIndex{N}})
#     ret = ""
#     for d in 1:N
#         s = r.start.I[d]
#         e = r.stop.I[d]
#         ret *= "$s:$e_"
#     end
#     return ret[1:end-1]
# end

function unit_range2string(idxes::Union{Tuple,Vector})
    ret = ""
    for idx in idxes
        ret *= "$(first(idx)-1)-$(last(idx))_"
    end
    return ret[1:end-1]
end

function cartesian_range2string(r::CartesianIndices{N}) where N
    ret = ""
    for i in 1:3
        ret *= "$(first(r)[i]-1)-$(last(r)[i])_"
    end
    return ret[1:end-1]
end

"""
transform string to UnitRange
the format of string should look like this:
2968-3480_1776-2288_16912-17424
-1024--896_-1024--896_1428-1429
"""
function string2unit_range( str::AbstractString )
    groups = match(r"(-?\d+)-(-?\d+)_(-?\d+)-(-?\d+)_(-?\d+)-(-?\d+)(?:\.gz)?$", str)
    idxes = map(Meta.parse, groups.captures)
    [idxes[1]+1:idxes[2], idxes[3]+1:idxes[4], idxes[5]+1:idxes[6]]
end

"""
    adjust bounding box range when fitting in new subarray
"""
function union(globalRange::CartesianIndices, idxes::CartesianIndices)
    start = map(min, first(globalRange).I, first(idxes).I)
    stop  = map(max, last(globalRange).I,  last(idxes).I)
    return CartesianIndices( map((x,y)->x:y, start.I, stop.I) )
end
function union!(r1::CartesianIndices, r2::CartesianIndices)
    r1 = union(r1, r2)
end

"""
    transform x1:x2_y1:y2_z1:z2 style string to CartesianRange
"""
function string2cartesian_range( s::String )
    secs = split(s, "_")
    starts = map( x->Meta.parse(split(x,"-")[1])+1, secs )
    stops  = map( x->Meta.parse(split(x,"-")[2]), secs )
    CartesianIndices( map((x,y)->x:y, starts, stops) )
end

function index2unit_range(x::UnitRange)
    x
end
function index2unit_range(x::Int)
    x:x
end

end # end of module


text:
using JLD, Plots, LaTeXStrings
font = Plots.font("DejaVu Sans", 24)
pyplot(size=(1200,800),linewidth=4,xscale=:log10,yscale=:log10,
       xtickfont=font,ytickfont=font,markersize=8,
       legendfont=font,guidefont=font,titlefont=font,
       xlabel="Δx",ylabel="||error||",
       title="Poisson Equation Convergence",grid=false)

function main()
    # import jld files
    jldfile = "data/TEMPpoisson_validation.jld"

    o1L1arr   = load(jldfile,"o1L1arr")
    o1L2arr   = load(jldfile,"o1L2arr")
    o1LInfarr = load(jldfile,"o1LInfarr")
    o2L1arr   = load(jldfile,"o2L1arr")
    o2L2arr   = load(jldfile,"o2L2arr")
    o2LInfarr = load(jldfile,"o2LInfarr")
    harr      = load(jldfile,"harr")
    κarr      = load(jldfile,"κarr")
    o1timearr = load(jldfile,"o1timearr")
    o2timearr = load(jldfile,"o2timearr")

    # plot files
    MS = 10         # marker size
    LW = 3          # line width
    FS = 15         # font size
    FStitle = FS+5

    # create plots
    # error plots
    p1 = plot(harr,o1L1arr,color=:red,linestyle=:dash,label="L1 (linear ϕ)")
    plot!(p1,harr,o1L2arr,color=:blue,linestyle=:dash,label="L2")
    plot!(p1,harr,o1LInfarr,color=:green,linestyle=:dash,label="Linf")
    scatter!(p1,harr,o1L1arr,label="",color=:black)
    scatter!(p1,harr,o1L2arr,label="",color=:black)
    scatter!(p1,harr,o1LInfarr,label="",color=:black)
    plot!(p1,harr,harr.^2,color=:black,linestyle=:dash,label="O(Δx^2)")

    plot!(p1,harr,o2L1arr,color=:red,label="L1 (quadratic ϕ)")
    plot!(p1,harr,o2L2arr,color=:blue,label="L2")
    plot!(p1,harr,o2LInfarr,color=:green,label="Linf")
    scatter!(p1,harr,o2L1arr,color=:black,label="")
    scatter!(p1,harr,o2L2arr,color=:black,label="")
    scatter!(p1,harr,o2LInfarr,color=:black,label="")
    plot!(p1,harr,0.01*harr.^4,color=:black,label="O(Δx^4)")

    # time plots
    p2 = plot(title="Poisson Time-Accuracy",xlabel="||error||",ylabel="elapsed time (s)")
    colorArr = [:red,:green,:blue,:orange,:yellow,:purple]
    for i=1:length(o1timearr)
        scatter!([o1LInfarr[i]],[o1timearr[i]],
                 color=colorArr[i],label="",markershape=:star,
                 markersize=16)
        scatter!([o2LInfarr[i]],[o2timearr[i]],
                 color=colorArr[i],label="",markershape=:circle,
                 markersize=12)
    end

    # condition number plot
    p3 = plot(harr,κarr,label="",linewidth=LW,color=:blue)
    scatter!(p3,harr,κarr,label="",ylabel=L"\kappa",
            title="Poisson condition number",markersize=MS,
            color=:blue)

    # export plots
    run(`mkdir -p figures`)
    plot(p1)
    savefig(p1,"figures/TEMPpoisson_validation.png")
    plot(p2)
    savefig(p2,"figures/TEMPpoisson_time_accuracy.png")
    plot(p3)
    savefig(p3,"figures/TEMP_poisson_cond.png")

    nothing
end

main()

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
"""
    abstract type StringWrapper <: AbstractString end

Provides default `sw::StringWrapper<:AbstractString` API methods defering to `representation(sw)` call, by default `sw.x`.
"""
abstract type StringWrapper <: AbstractString end
@inline representation(sw::StringWrapper) = sw.x

@inline Base.@propagate_inbounds Base.firstindex(sw::StringWrapper) =
    firstindex(representation(sw))
@inline Base.@propagate_inbounds Base.lastindex(sw::StringWrapper) =
    lastindex(representation(sw))
@inline Base.@propagate_inbounds Base.length(sw::StringWrapper) =
    length(representation(sw))
@inline Base.@propagate_inbounds Base.ncodeunits(sw::StringWrapper) =
    ncodeunits(representation(sw))
@inline Base.@propagate_inbounds Base.sizeof(sw::StringWrapper) =
    sizeof(representation(sw))

@inline Base.@propagate_inbounds Base.getindex(sw::StringWrapper,i::Integer) =
    getindex(representation(sw),i)
@inline Base.@propagate_inbounds Base.codeunit(s::StringWrapper, i::Integer) =
    codeunit(representation(sw), i)
@inline Base.@propagate_inbounds Base.isvalid(sw::StringWrapper,i::Int) =
    isvalid(representation(sw), i)

@inline Base.@propagate_inbounds Base.prevind(sw::StringWrapper,i::Int,n::Int) =
    prevind(representation(sw),i,n)
@inline Base.@propagate_inbounds Base.nextind(sw::StringWrapper,i::Int,n::Int) =
    nextind(representation(sw),i,n)
@inline Base.@propagate_inbounds Base.prevind(sw::StringWrapper,i::Int) =
    prevind(representation(sw),i)
@inline Base.@propagate_inbounds Base.nextind(sw::StringWrapper,i::Int) =
    nextind(representation(sw),i)
@inline Base.@propagate_inbounds Base.thisind(sw::StringWrapper,i::Int) =
    thisind(representation(sw), i)

@inline Base.@propagate_inbounds Base.iterate(sw::StringWrapper) =
    iterate(representation(sw))
@inline Base.@propagate_inbounds Base.iterate(sw::StringWrapper,i::Integer) =
    iterate(representation(sw),i)
@inline Base.@propagate_inbounds Base.SubString(sw::StringWrapper,start::Int,stop::Int) =
    SubString(representation(sw),start,stop)