text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities
import MutableArithmetics
const MA = MutableArithmetics

w = MOI.VariableIndex(0)
fw = MOI.SingleVariable(w)
x = MOI.VariableIndex(1)
fx = MOI.SingleVariable(x)
y = MOI.VariableIndex(2)
fy = MOI.SingleVariable(y)
z = MOI.VariableIndex(3)
fz = MOI.SingleVariable(z)

# Number-like but not subtype of `Number`
struct NonNumber
    value::Int
end
Base.:*(a::NonNumber, b::NonNumber) = NonNumber(a.value * b.value)
Base.:+(a::NonNumber, b::NonNumber) = NonNumber(a.value + b.value)
Base.zero(::Type{NonNumber}) = NonNumber(0)
Base.isapprox(a::NonNumber, b::NonNumber; kws...) = isapprox(a.value, b.value; kws...)
@testset "NonNumber" begin
    two = NonNumber(2)
    three = NonNumber(3)
    six = NonNumber(6)
    three_x = MOIU.operate(*, NonNumber, three, fx)
    six_x = MOIU.operate(*, NonNumber, six, fx)
    @test six_x ≈ two * three_x
    @test six_x ≈ three_x * two
end

@testset "Vectorization" begin
    g = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 1],
                                                       MOI.ScalarAffineTerm.([5, 2],
                                                                             [y, x])),
                                                       [3, 1, 4])
    @testset "vectorize" begin
        g1 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2, x)], 3)
        g2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Int}[], 1)
        g3 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(5, y)], 4)
        @test g ≈ MOIU.vectorize([g1, g2, g3])
        vov = MOIU.vectorize(MOI.SingleVariable[])
        @test MOI.output_dimension(vov) == 0
        @test vov isa MOI.VectorOfVariables
        aff = MOIU.vectorize(MOI.ScalarAffineFunction{Int}[])
        @test MOI.output_dimension(aff) == 0
        @test aff isa MOI.VectorAffineFunction{Int}
        quad = MOIU.vectorize(MOI.ScalarQuadraticFunction{Int}[])
        @test MOI.output_dimension(quad) == 0
        @test quad isa MOI.VectorQuadraticFunction{Int}
    end
    @testset "operate vcat" begin
        v = MOI.VectorOfVariables([y, w])
        wf = MOI.SingleVariable(w)
        xf = MOI.SingleVariable(x)
        @testset "Variable with $T" for T in [Int, Float64]
            @test MOI.VectorOfVariables == MOIU.promote_operation(vcat, T, typeof(wf), typeof(v), typeof(xf))
            vov = MOIU.operate(vcat, T, wf, v, xf)
            @test vov.variables == [w, y, w, x]
            @test MOI.VectorOfVariables == MOIU.promote_operation(vcat, T, typeof(v), typeof(wf), typeof(xf))
            vov = MOIU.operate(vcat, T, v, wf, xf)
            @test vov.variables == [y, w, w, x]
            @test MOI.VectorOfVariables == MOIU.promote_operation(vcat, T, typeof(wf), typeof(xf), typeof(v))
            vov = MOIU.operate(vcat, T, wf, xf, v)
            @test vov.variables == [w, x, y, w]
        end
        f = MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([2, 4], [x, z]), 5)
        g = MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.([3, 1], MOI.ScalarAffineTerm.([5, 2], [y, x])),
            [3, 1, 4])
        @testset "Affine" begin
            @test MOIU.promote_operation(vcat, Int, typeof(wf), typeof(f),
                                         typeof(v), Int, typeof(g), typeof(xf),
                                         Int) == MOI.VectorAffineFunction{Int}
            F = MOIU.operate(vcat, Int, wf, f, v, 3, g, xf, -4)
            expected_terms = MOI.VectorAffineTerm.(
                [1, 2, 2, 3, 4, 8, 6, 9],
                MOI.ScalarAffineTerm.([1, 2, 4, 1, 1, 5, 2, 1],
                                      [w, x, z, y, w, y, x, x]))
            expected_constants = [0, 5, 0, 0, 3, 3, 1, 4, 0, -4]
            F = MOIU.operate(vcat, Int, wf, f, v, 3, g, xf, -4)
            @test F.terms == expected_terms
            @test F.constants == expected_constants
        end
        @testset "Quadratic" begin
            @test MOIU.promote_operation(
                vcat, Int, MOI.VectorQuadraticFunction{Int}, typeof(wf),
                typeof(f), typeof(v), Int, MOI.ScalarQuadraticFunction{Int},
                typeof(g), typeof(xf), Int) == MOI.VectorQuadraticFunction{Int}
            # TODO
        end
    end
end
@testset "MultirowChange construction" begin
    chg1 = MOI.MultirowChange(w, [(Int32(2), 2.0), (Int32(1), 3.0)])
    chg2 = MOI.MultirowChange(w, [(Int64(2), 2.0), (Int64(1), 3.0)])
    @test chg1.variable == chg2.variable
    @test chg1.new_coefficients == chg2.new_coefficients
end
@testset "VectorAffineTerm/VectorQuadraticTerm construction" begin
    scalaraffine = MOI.ScalarAffineTerm(2.0, z)
    @test MOI.VectorAffineTerm(Int32(3), scalaraffine) === MOI.VectorAffineTerm(Int64(3), scalaraffine)
    @test MOI.VectorAffineTerm{Float64}(Int32(3), scalaraffine) === MOI.VectorAffineTerm(Int64(3), scalaraffine)
    scalarquad = MOI.ScalarQuadraticTerm(2.0, y, z)
    @test MOI.VectorQuadraticTerm(Int32(3), scalarquad) === MOI.VectorQuadraticTerm(Int64(3), scalarquad)
    @test MOI.VectorQuadraticTerm{Float64}(Int32(3), scalarquad) === MOI.VectorQuadraticTerm(Int64(3), scalarquad)
end
@testset "eval_variables" begin
    # We do tests twice to make sure the function is not modified
    vals = Dict(w=>0, x=>3, y=>1, z=>5)
    fsv = MOI.SingleVariable(z)
    @test MOI.output_dimension(fsv) == 1
    @test MOIU.eval_variables(vi -> vals[vi], fsv) ≈ 5
    @test MOIU.eval_variables(vi -> vals[vi], fsv) ≈ 5
    fvv = MOI.VectorOfVariables([x, z, y])
    @test MOI.output_dimension(fvv) == 3
    @test MOIU.eval_variables(vi -> vals[vi], fvv) ≈ [3, 5, 1]
    @test MOIU.eval_variables(vi -> vals[vi], fvv) ≈ [3, 5, 1]
    fsa = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(3.0, z), MOI.ScalarAffineTerm(2.0, y)], 2.0)
    @test MOI.output_dimension(fsa) == 1
    @test MOIU.eval_variables(vi -> vals[vi], fsa) ≈ 22
    @test MOIU.eval_variables(vi -> vals[vi], fsa) ≈ 22
    fva = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 1, 2], MOI.ScalarAffineTerm.([1.0, 3.0, 2.0], [x, z, y])), [-3.0, 2.0])
    @test MOI.output_dimension(fva) == 2
    @test MOIU.eval_variables(vi -> vals[vi], fva) ≈ [12, 7]
    @test MOIU.eval_variables(vi -> vals[vi], fva) ≈ [12, 7]
    fsq = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(1.0, [x, y]),
                                      MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]), -3.0)
    @test MOI.output_dimension(fsq) == 1
    @test MOIU.eval_variables(vi -> vals[vi], fsq) ≈ 16
    @test MOIU.eval_variables(vi -> vals[vi], fsq) ≈ 16
    fvq = MOI.VectorQuadraticFunction(MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
                                      MOI.VectorQuadraticTerm.([1, 2, 2], MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y])), [-3.0, -2.0])
    @test MOI.output_dimension(fvq) == 2
    @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
    @test MOIU.eval_variables(vi -> vals[vi], fvq) ≈ [13, 1]
end
@testset "substitute_variables" begin
    # We do tests twice to make sure the function is not modified
    subs = Dict(w => 1.0fy + 1.0fz, x => 2.0fy + 1.0, y => 1.0fy, z => -1.0fw)
    vals = Dict(w => 0.0, x => 3.0, y => 1.0, z => 5.0)
    subs_vals = Dict(w => 6.0, x => 3.0, y => 1.0, z => 0.0)
    fsa = fx + 3.0fz + 2.0fy + 2.0
    subs_sa = -3.0fw + 4.0fy + 3.0
    @test MOIU.eval_variables(vi -> subs_vals[vi], fsa) == MOIU.eval_variables(vi -> vals[vi], subs_sa)
    @test MOIU.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa
    @test MOIU.substitute_variables(vi -> subs[vi], fsa) ≈ subs_sa
    fva = MOIU.operate(vcat, Float64, 3.0fz - 3.0, fx + 2.0fy + 2.0)
    subs_va = MOIU.operate(vcat, Float64, -3.0fw - 3.0, 4.0fy + 3.0)
    @test MOIU.eval_variables(vi -> subs_vals[vi], fva) == MOIU.eval_variables(vi -> vals[vi], subs_va)
    @test MOIU.substitute_variables(vi -> subs[vi], fva) ≈ subs_va
    @test MOIU.substitute_variables(vi -> subs[vi], fva) ≈ subs_va
    fsq = 1.0fx  + 1.0fy + 1.0fx * fz + 1.0fw * fz     + 1.0fw   * fy + 2.0fw     * fw      - 3.0
    #     2y + 1 +     y - 2yw - w    + (y + z) * (-w) + (y + z) * y  + 2 (y + z) * (y + z) - 3
    #     2y     +     y - 2yw - w    - yw - zw        + y^2 + zy     + 2y^2 + 2z^2 + 4yz   - 2
    #     3y - w + 3y^2 + 2z^2 - 3yw - zw + 5yz - 2
    subs_sq = 3.0fy - 1.0fw + 3.0fy * fy + 2.0fz * fz - 3.0fy * fw - 1.0fz * fw + 5.0fy * fz - 2.0
    @test MOIU.eval_variables(vi -> subs_vals[vi], fsq) == MOIU.eval_variables(vi -> vals[vi], subs_sq)
    @test MOIU.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq
    @test MOIU.substitute_variables(vi -> subs[vi], fsq) ≈ subs_sq
    fvq = MOIU.operate(vcat, Float64, 1.0fy + 1.0fx*fz - 3.0, 1.0fx + 1.0fw*fz + 1.0fw*fy - 2.0)
    subs_vq = MOIU.operate(vcat, Float64,
        1.0fy - fw - 2.0fy*fw - 3.0,
        2.0fy + 1.0fy*fy - 1.0fw*fy - 1.0fw*fz + 1.0fy*fz - 1.0)
    @test MOIU.eval_variables(vi -> subs_vals[vi], fvq) == MOIU.eval_variables(vi -> vals[vi], subs_vq)
    @test MOIU.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq
    @test MOIU.substitute_variables(vi -> subs[vi], fvq) ≈ subs_vq

    complex_aff = 2.0im * fy
    # Test that variables can be substituted for `MOI.ScalarAffineFunction{S}`
    # in a `MOI.ScalarAffineFunction{T}` where `S != T`.
    @test MOIU.substitute_variables(vi -> 1.5fx, complex_aff) ≈ 3.0im * fx
    float_aff = 2.0 * fy
    @test MOIU.substitute_variables(vi -> 3fx, float_aff) ≈ 6.0 * fx
    complex_quad = 1.5im * fy * fy
    @test MOIU.substitute_variables(vi -> 2fx, complex_quad) ≈ 6.0im * fx * fx
    int_quad = 3 * fy * fx
    @test MOIU.substitute_variables(vi -> true * fy, int_quad) ≈ 3 * fy * fy
end
@testset "map_indices" begin
    fsq = MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.(1.0, [x, y]),
                                      MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y]), -3.0)
    index_map = Dict(x => y, y => z, w => w, z => x)
    gsq = MOIU.map_indices(index_map, fsq)
    sats = MOI.ScalarAffineTerm.(1.0, [y, z])
    sqts = MOI.ScalarQuadraticTerm.(1.0, [y, w, w], [x, x, z])
    @test gsq.affine_terms == sats
    @test gsq.quadratic_terms == sqts
    @test gsq.constant == -3.
    fvq = MOI.VectorQuadraticFunction(MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.(1.0, [x, y])),
                                      MOI.VectorQuadraticTerm.([1, 2, 2], MOI.ScalarQuadraticTerm.(1.0, [x, w, w], [z, z, y])), [-3.0, -2.0])
    gvq = MOIU.map_indices(index_map, fvq)
    @test gvq.affine_terms == MOI.VectorAffineTerm.([2, 1], sats)
    @test gvq.quadratic_terms == MOI.VectorQuadraticTerm.([1, 2, 2], sqts)
    @test MOIU.constant_vector(gvq) == [-3., -2.]

    @testset "Constants" begin
        @test MOIU.map_indices(index_map, :s) == :s
        @test MOIU.map_indices(index_map, [:a, :b]) == [:a, :b]
        @test MOIU.map_indices(index_map, "s") == "s"
        @test MOIU.map_indices(index_map, ["a", "b"]) == ["a", "b"]
    end
end
@testset "Conversion VectorOfVariables -> VectorAffineFunction" begin
    f = MOI.VectorAffineFunction{Int}(MOI.VectorOfVariables([z, x, y]))
    @test f isa MOI.VectorAffineFunction{Int}
    @test f.terms == MOI.VectorAffineTerm.(1:3, MOI.ScalarAffineTerm.(ones(Int, 3), [z, x, y]))
    @test all(iszero.(MOIU.constant_vector(f)))
    f = MOI.VectorAffineFunction{Float64}(MOI.VectorOfVariables([x, w]))
    @test f isa MOI.VectorAffineFunction{Float64}
    @test f.terms == MOI.VectorAffineTerm.(1:2, MOI.ScalarAffineTerm.(1.0, [x, w]))
    @test all(iszero.(MOIU.constant_vector(f)))
end
@testset "Iteration and indexing on VectorOfVariables" begin
    f = MOI.VectorOfVariables([z, w, x, y])
    it = MOIU.eachscalar(f)
    @test length(it) == 4
    @test eltype(it) == MOI.SingleVariable
    @test collect(it) == [MOI.SingleVariable(z), MOI.SingleVariable(w), MOI.SingleVariable(x), MOI.SingleVariable(y)]
    @test it[2] == MOI.SingleVariable(w)
    @test it[end] == MOI.SingleVariable(y)
end
@testset "Indexing on VectorAffineFunction" begin
    f = MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 1, 3, 2, 2, 1, 3, 1, 2],
                                                       MOI.ScalarAffineTerm.([1, 7, 2, 9, 3, 1, 6, 4, 1],
                                                                             [x, y, z, z, y, z, x, x, y])),
                                 [2, 7, 5])
    it = MOIU.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == MOI.ScalarAffineFunction{Int}
    g = it[2]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([1, 9, 3, 1], [x, z, y, y])
    @test g.constant == 7
    g = it[1]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([7, 1, 4], [y, z, x])
    @test g.constant == 2
    g = it[end]
    @test g isa MOI.ScalarAffineFunction
    @test g.terms == MOI.ScalarAffineTerm.([2, 6], [z, x])
    @test g.constant == 5
    h = it[[3, 1]]
    @test h isa MOI.VectorAffineFunction
    @test h.terms == MOI.VectorAffineTerm.([1, 1, 2, 2, 2], MOI.ScalarAffineTerm.([2, 6, 7, 1, 4], [z, x, y, z, x]))
    @test MOIU.constant_vector(h) == [5, 2]
    F = MOIU.operate(vcat, Int, it[[1, 2]], it[3])
    @test F isa MOI.VectorAffineFunction{Int}
    @test F.terms == MOI.VectorAffineTerm.([1, 1, 1, 2, 2, 2, 2, 3, 3], MOI.ScalarAffineTerm.([7, 1, 4, 1, 9, 3, 1, 2, 6], [y, z, x, x, z, y, y, z, x]))
    @test MOIU.constant_vector(F) == MOIU.constant_vector(f)
end
@testset "Indexing on VectorQuadraticFunction" begin
    f = MOI.VectorQuadraticFunction(MOI.VectorAffineTerm.([2, 1, 3, 2, 2],
                                                          MOI.ScalarAffineTerm.([1, 7, 2, 9, 3],
                                                                                [x, y, z, z, y])),
                                    MOI.VectorQuadraticTerm.([2, 3, 1, 2],
                                                             MOI.ScalarQuadraticTerm.([1, 6, 4, 3],
                                                                                [z, x, x, y],
                                                                                [y, z, z, y])),
                                    [2, 7, 5])
    it = MOIU.eachscalar(f)
    @test length(it) == 3
    @test eltype(it) == MOI.ScalarQuadraticFunction{Int}
    g = it[2]
    @test g isa MOI.ScalarQuadraticFunction
    @test g.affine_terms == MOI.ScalarAffineTerm.([1, 9, 3], [x, z, y])
    @test g.quadratic_terms == MOI.ScalarQuadraticTerm.([1, 3], [z, y], [y, y])
    @test g.constant == 7
    g = it[end]
    @test g isa MOI.ScalarQuadraticFunction
    @test g.affine_terms == MOI.ScalarAffineTerm.([2], [z])
    @test g.quadratic_terms == MOI.ScalarQuadraticTerm.([6], [x], [z])
    @test g.constant == 5
end
@testset "Scalar" begin
    @testset "Variable" begin
        f = MOI.SingleVariable(MOI.VariableIndex(0))
        g = MOI.SingleVariable(MOI.VariableIndex(1))
        @testset "one" begin
            @test !isone(f)
            @test !isone(g)
        end
        @testset "zero" begin
            @test !iszero(f)
            @test !iszero(g)
            @test f + 1 ≈ 1 + f
            @test (f + 1.0) - 1.0 ≈ (2.0f) / 2.0
            @test (f - 1.0) + 1.0 ≈ (2.0f) / 2.0
            @test (1.0 + f) - 1.0 ≈ (f * 2.0) / 2.0
            @test 1.0 - (1.0 - f) ≈ (f / 2.0) * 2.0
        end
        @testset "complex operations" begin
            @test real(f) === f
            @test MA.promote_operation(real, typeof(f)) == typeof(f)
            @testset "imag with $T" for T in [Int, Int32, Float64]
                z = MOIU.operate(imag, T, f)
                @test z isa MOI.ScalarAffineFunction{T}
                @test iszero(z)
                @test MOIU.promote_operation(imag, T, typeof(f)) == typeof(z)
            end
            @test conj(f) === f
            @test MA.promote_operation(conj, typeof(f)) == typeof(f)
        end
    end
    @testset "Affine" begin
        @testset "zero" begin
            f = @inferred MOIU.zero(MOI.ScalarAffineFunction{Float64})
            @test iszero(f)
            @test MOIU.isapprox_zero(f, 1e-16)
        end
        @testset "complex operations" begin
            fx = MOI.SingleVariable(MOI.VariableIndex(1))
            fy = MOI.SingleVariable(MOI.VariableIndex(2))
            r = 3fx + 4fy
            c = 2fx + 5fy
            f = (1 + 0im) * r + c * im
            @test real(f) ≈ r
            @test MOIU.operate(imag, Int, f) ≈ c
            @test imag(f) ≈ c
            @test conj(f) ≈ (1 + 0im) * r - c * im
            @test MA.promote_operation(real, typeof(f)) == typeof(r)
            @test MA.promote_operation(imag, typeof(f)) == typeof(c)
            @test MA.promote_operation(conj, typeof(f)) == typeof(f)
        end
        @testset "promote_operation" begin
            @test MOIU.promote_operation(
                -, Int, MOI.SingleVariable
            ) == MOI.ScalarAffineFunction{Int}
            @test MOIU.promote_operation(
                -, Int, MOI.ScalarAffineFunction{Int}
            ) == MOI.ScalarAffineFunction{Int}
            @test MOIU.promote_operation(+, Float64, MOI.SingleVariable,
                                         MOI.SingleVariable) == MOI.ScalarAffineFunction{Float64}
            @test MOIU.promote_operation(+, Float64,
                                         MOI.ScalarAffineFunction{Float64},
                                         Float64) == MOI.ScalarAffineFunction{Float64}
            @test MOIU.promote_operation(+, Int,
                                         MOI.ScalarAffineFunction{Int},
                                         MOI.ScalarAffineFunction{Int}) == MOI.ScalarAffineFunction{Int}
        end
        @testset "Comparison" begin
            @test MOIU.operate(+, Float64, MOI.SingleVariable(x),
                               MOI.SingleVariable(z)) + 1.0 ≈
                  MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 1e-7, 1], [x, y, z]), 1.0) atol=1e-6
            f1 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x), MOI.ScalarAffineTerm(1e-7, y)], 1.0)
            f2 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 1.0)
            @test f1 ≈ f2 atol=1e-6
            fdiff = f1 - f2
            @testset "With iszero" begin
                @test !iszero(fdiff)
                @test iszero(f1 - f1)
                @test iszero(f2 - f2)
            end
            @testset "With tolerance" begin
                MOIU.canonicalize!(fdiff)
                @test !MOIU.isapprox_zero(fdiff, 1e-8)
                @test MOIU.isapprox_zero(fdiff, 1e-6)
            end
        end
        @testset "canonical" begin
            f = MOIU.canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, 1, 3, -2, -3], [y, x, z, x, z]), 5))
            @test MOI.output_dimension(f) == 1
            @test f.terms == MOI.ScalarAffineTerm.([-1, 2], [x, y])
            @test f.constant == 5
            f = MOIU.canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 3, 1, 2, -3, 2],
                                                                              [w, y, w, x,  x, z]), 2) +
                               MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1, -2, -2, 3, 2],
                                                                              [ y,  z,  w, x, y]), 3))
            @test f.terms == MOI.ScalarAffineTerm.([2, 4], [x, y])
            @test f.constant == 5
        end
        f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.5], [x, y]), 0.5)
        @testset "convert" begin
            @test_throws InexactError convert(MOI.SingleVariable, f)
            @test_throws InexactError MOIU.convert_approx(MOI.SingleVariable, f)
            @test MOIU.convert_approx(MOI.SingleVariable, f, tol = 0.5) == MOI.SingleVariable(x)
            @test convert(typeof(f), f) === f
            quad_f = MOI.ScalarQuadraticFunction(f.terms,
                                                 MOI.ScalarQuadraticTerm{Float64}[],
                                                 f.constant)
            @test convert(MOI.ScalarQuadraticFunction{Float64}, f) ≈ quad_f
            for g in [convert(MOI.ScalarAffineFunction{Float64}, MOI.SingleVariable(x)),
                      convert(MOI.ScalarAffineFunction{Float64}, 1MOI.SingleVariable(x))]
                @test g isa MOI.ScalarAffineFunction{Float64}
                @test convert(MOI.SingleVariable, g) == MOI.SingleVariable(x)
                @test MOIU.convert_approx(MOI.SingleVariable, g) == MOI.SingleVariable(x)
            end
        end
        @testset "operate with Float64 coefficient type" begin
            f = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 4.0],
                                                               [x, y]),
                                         5.0)
            @test f ≈ 2.0f / 2.0
        end
        @testset "operate with Int coefficient type" begin
            f = MOIU.canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 3, 1, 2, -3, 2],
                                                                              [w, y, w, x,  x, z]), 2) +
                               MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1, -2, -2, 3, 2],
                                                                              [ y,  z,  w, x, y]), 3))
            @test f === +f
            @test f ≈ MOI.SingleVariable(x) + MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 4],
                                                                                [x, y]), 5)
            @test f ≈ f * 1
            @test f ≈ MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1, 2], [x, y]), 2) * 2 + 1
            @test f ≈ MOI.SingleVariable(x) - MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-1, -4], [x, y]), -5)
            @test f ≈ MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3, 4], [x, y]), 5) - MOI.SingleVariable(x)
        end
        @testset "modification" begin
            f = MOIU.modify_function(f, MOI.ScalarConstantChange(6))
            @test f.constant == 6
            g = deepcopy(f)
            @test g ≈ f
            f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(y, 3))
            @test !(g ≈ f)
            @test g.terms == MOI.ScalarAffineTerm.([2, 4], [x, y])
            @test f.terms == MOI.ScalarAffineTerm.([2, 3], [x, y])
            f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(x, 0))
            @test f.terms == MOI.ScalarAffineTerm.([3], [y])
        end
    end
    @testset "Quadratic" begin
        @testset "zero" begin
            f = @inferred MOIU.zero(MOI.ScalarQuadraticFunction{Float64})
            @test MOIU.isapprox_zero(f, 1e-16)
        end
        @testset "promote_operation" begin
            @test MOIU.promote_operation(
                -, Int, MOI.ScalarQuadraticFunction{Int}
            ) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(+, Int,
                                         MOI.ScalarQuadraticFunction{Int},
                                         MOI.ScalarQuadraticFunction{Int}) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(+, Int,
                                         MOI.ScalarQuadraticFunction{Int},
                                         MOI.ScalarAffineFunction{Int}) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(+, Int,
                                         MOI.ScalarAffineFunction{Int},
                                         MOI.ScalarQuadraticFunction{Int}) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(*, Int,
                                         MOI.SingleVariable,
                                         MOI.SingleVariable) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(*, Float64,
                                         MOI.SingleVariable,
                                         MOI.ScalarAffineFunction{Float64}) == MOI.ScalarQuadraticFunction{Float64}
            @test MOIU.promote_operation(*, Int,
                                         MOI.ScalarAffineFunction{Int},
                                         MOI.SingleVariable) == MOI.ScalarQuadraticFunction{Int}
            @test MOIU.promote_operation(*, Float64,
                                         MOI.ScalarAffineFunction{Float64},
                                         MOI.ScalarAffineFunction{Float64}) == MOI.ScalarQuadraticFunction{Float64}
            @test MOIU.promote_operation(/, Float64,
                                         MOI.ScalarQuadraticFunction{Float64},
                                         Float64) == MOI.ScalarQuadraticFunction{Float64}
        end
        f = 7 + 3fx + 1fx * fx + 2fy * fy + 3fx * fy
        MOIU.canonicalize!(f)
        @test MOI.output_dimension(f) == 1
        @testset "Comparison" begin
            @testset "With iszero" begin
                @test !iszero(f)
                @test iszero(0 * f)
                @test iszero(f - f)
            end
            @testset "With tolerance" begin
                @test !MOIU.isapprox_zero(f, 1e-8)
                # Test isapprox_zero with zero terms
                @test MOIU.isapprox_zero(0 * f, 1e-8)
                g = 1.0fx * fy - (1 + 1e-6) * fy * fx
                MOIU.canonicalize!(g)
                @test MOIU.isapprox_zero(g, 1e-5)
                @test !MOIU.isapprox_zero(g, 1e-7)
            end
        end
        @testset "convert" begin
            @test_throws InexactError convert(MOI.SingleVariable, f)
            @test_throws InexactError convert(MOI.ScalarAffineFunction{Int},
                                              f)
            g = convert(MOI.ScalarQuadraticFunction{Float64}, fx)
            @test convert(MOI.SingleVariable, g) == fx
        end
        @testset "Power" begin
            @testset "Affine" begin
                aff = 1fx + 2 + fy
                @test isone(@inferred aff^0)
                @test convert(typeof(f), aff) ≈ @inferred aff^1
                @test aff * aff ≈ @inferred aff^2
                err = ArgumentError("Cannot take $(typeof(aff)) to the power 3.")
                @test_throws err aff^3
            end
            @testset "Quadratic" begin
                @test isone(@inferred f^0)
                @test f ≈ @inferred f^1
                err = ArgumentError("Cannot take $(typeof(f)) to the power 2.")
                @test_throws err f^2
            end
        end
        @testset "operate" begin
            @testset "No zero affine term" begin
                # Test that no affine 0y term is created when multiplying 1fx by fy
                for fxfy in [1fx * fy, fx * 1fy]
                    @test isempty(fxfy.affine_terms)
                    @test length(fxfy.quadratic_terms) == 1
                    @test fxfy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, x, y) ||
                        fxfy.quadratic_terms[1] == MOI.ScalarQuadraticTerm(1, y, x)
                end
                for fxfx in [1fx * fx, fx * 1fx]
                    @test isempty(fxfx.affine_terms)
                    @test length(fxfx.quadratic_terms) == 1
                    @test fxfx.quadratic_terms[1] == MOI.ScalarQuadraticTerm(2, x, x)
                end
            end
            @testset "operate!" begin
                q = 1.0fx + 1.0fy + (1.0fx) * fz + (1.0fw) * fz
                @test q ≈ 1.0fx + 1.0fy + (1.0fw) * fz + (1.0fx) * fz
                # This calls
                aff = 1.0fx + 1.0fy
                # which tries to mutate `aff`, gets a quadratic expression
                # and mutate it with the remaining term
                @test MOIU.operate!(+, Float64, aff, (1.0fx) * fz, (1.0fw) * fz) ≈ q
            end
            @test f ≈ 7 + (fx + 2fy) * (1fx + fy) + 3fx
            @test f ≈ -(-7 - 3fx) + (fx + 2fy) * (1fx + fy)
            @test f ≈ -((fx + 2fy) * (MOIU.operate(-, Int, fx) - fy)) + 3fx + 7
            @test f ≈ 7 + MOIU.operate(*, Int, fx, fx) + 3fx * (fy + 1) + 2fy * fy
            @test f ≈ (fx + 2) * (fx + 1) + (fy + 1) * (2fy + 3fx) + (5 - 3fx - 2fy)
            @test f ≈ begin
                MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(3, x)],
                                            MOI.ScalarQuadraticTerm.([2], [x], [x]), 4) +
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                        MOI.ScalarQuadraticTerm.([4, 3], [y, x], [y, y]), 3)
            end
            @test f ≈ begin
                MOI.ScalarQuadraticFunction([MOI.ScalarAffineTerm(3, x)],
                                            MOI.ScalarQuadraticTerm.([2], [x], [x]), 10) -
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                            MOI.ScalarQuadraticTerm.([-4, -3], [y, x], [y, y]), 3)
            end
            @test f ≈ begin
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) +
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 2)
            end
            @test f ≈ begin
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5) -
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                            MOI.ScalarQuadraticTerm.([-2, -4, -3], [x, y, x], [x, y, y]), -2)
            end
            @test f ≈ begin
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 2) +
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 5)
            end
            @test f ≈ begin
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Int}[],
                                            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 12) -
                MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(-3, x)], 5)
            end
            @test f ≈ begin
                MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([2], [x]),
                                            MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 7) +
                MOI.SingleVariable(x)
            end
            @test f ≈ MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm.([3], [x]),
                                          MOI.ScalarQuadraticTerm.([2, 4, 3], [x, y, x], [x, y, y]), 10) - 3
            @test f ≈ 2.0 * MOI.ScalarQuadraticFunction(
                MOI.ScalarAffineTerm.([3.0], [x]),
                MOI.ScalarQuadraticTerm.([2.0, 4.0, 3.0], [x, y, x], [x, y, y]), 7.0) / 2.0
        end
        @testset "modification" begin
            f = MOIU.modify_function(f, MOI.ScalarConstantChange(9))
            @test f.constant == 9
            f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(y, 0))
            @test f.affine_terms == MOI.ScalarAffineTerm.([3], [x])
            g = deepcopy(f)
            @test f ≈ g
            f = MOIU.modify_function(f, MOI.ScalarCoefficientChange(y, 2))
            @test !(f ≈ g)
            @test g.affine_terms == MOI.ScalarAffineTerm.([3], [x])
            @test f.affine_terms == MOI.ScalarAffineTerm.([3, 2], [x, y])
        end
    end
end
@testset "Vector" begin
    @testset "Variable" begin
        f = MOI.VectorOfVariables([MOI.VariableIndex(1), MOI.VariableIndex(2)])
        @testset "complex operations" begin
            @test real(f) === f
            @test MA.promote_operation(real, typeof(f)) == typeof(f)
            @testset "imag with $T" for T in [Int, Int32, Float64]
                z = MOIU.operate(imag, T, f)
                @test z isa MOI.VectorAffineFunction{T}
                @test isempty(z.terms)
                @test all(iszero, z.constants)
                @test MOI.output_dimension(z) == MOI.output_dimension(f)
                @test MOIU.promote_operation(imag, T, typeof(f)) == typeof(z)
            end
            @test conj(f) === f
            @test MA.promote_operation(conj, typeof(f)) == typeof(f)
        end
    end
    @testset "Affine" begin
        @testset "complex operations" begin
            fx = MOI.VectorOfVariables([MOI.VariableIndex(1)])
            fy = MOI.VectorOfVariables([MOI.VariableIndex(2)])
            r = 3fx + 4fy
            c = 2fx + 5fy
            f = (1 + 0im) * r + c * im
            @test real(f) ≈ r
            @test MOIU.operate(imag, Int, f) ≈ c
            @test imag(f) ≈ c
            @test conj(f) ≈ (1 + 0im) * r - c * im
            @test MA.promote_operation(real, typeof(f)) == typeof(r)
            @test MA.promote_operation(imag, typeof(f)) == typeof(c)
            @test MA.promote_operation(conj, typeof(f)) == typeof(f)
        end

        f = MOIU.canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 1, 2,  1,  1,  2, 2,  2, 2, 1, 1,  2, 1,  2],
                                                                          MOI.ScalarAffineTerm.([3, 2, 3, -3, -1, -2, 3, -2, 1, 3, 5, -2, 0, -1],
                                                                                [x, x, z,  y,  y,  x, y,  z, x, y, y,  x, x,  z])), [5, 7]))
        @test MOI.output_dimension(f) == 2
        @test f.terms == MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y]))
        @test MOIU.constant_vector(f) == [5, 7]
        f = MOIU.modify_function(f, MOI.VectorConstantChange([6, 8]))
        @test MOIU.constant_vector(f) == [6, 8]
        g = deepcopy(f)
        @test f ≈ g
        f = MOIU.modify_function(f, MOI.MultirowChange(y, [(2, 9)]))
        @test !(f ≈ g)
        @test f.terms == MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([2, 4, 9], [x, y, y]))
        @test g.terms == MOI.VectorAffineTerm.([1, 1, 2], MOI.ScalarAffineTerm.([2, 4, 3], [x, y, y]))
        f = MOIU.modify_function(f, MOI.MultirowChange(y, [(1, 0)]))
        @test f.terms == MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([2, 9], [x, y]))
    end
    @testset "Quadratic" begin
        f = MOI.VectorQuadraticFunction(MOI.VectorAffineTerm.([1, 2, 2], MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y])), MOI.VectorQuadraticTerm.([1, 1, 2], MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y])), [7, 3, 4])
        @test MOI.output_dimension(f) == 3
        f = MOIU.modify_function(f, MOI.VectorConstantChange([10, 11, 12]))
        @test MOIU.constant_vector(f) == [10, 11, 12]
        f = MOIU.modify_function(f, MOI.MultirowChange(y, [(2, 0), (1, 1)]))
        @test f.affine_terms == MOI.VectorAffineTerm.([1, 2, 1], MOI.ScalarAffineTerm.([3, 1, 1], [x, x, y]))
        g = deepcopy(f)
        f = MOIU.modify_function(f, MOI.MultirowChange(x, [(1, 0), (3, 4)]))
        @test f.affine_terms == MOI.VectorAffineTerm.([2, 1, 3], MOI.ScalarAffineTerm.([1, 1, 4], [x, y, x]))
        @test g.affine_terms == MOI.VectorAffineTerm.([1, 2, 1], MOI.ScalarAffineTerm.([3, 1, 1], [x, x, y]))
    end
end
@testset "Conversion to canonical form" begin
    function isapprox_ordered(f1::T, f2::T) where {T <: Union{MOI.ScalarAffineFunction, MOI.VectorAffineFunction}}
        ((MOI.term_indices.(f1.terms) == MOI.term_indices.(f2.terms)) &&
         (MOI.constant(f1) ≈ MOI.constant(f2)) &&
         (MOI.coefficient.(f1.terms) ≈ MOI.coefficient.(f2.terms)))
    end
    function isapprox_ordered(f1::T, f2::T) where {T <: Union{MOI.ScalarQuadraticFunction, MOI.VectorQuadraticFunction}}
        ((MOI.term_indices.(f1.affine_terms) == MOI.term_indices.(f2.affine_terms)) &&
         (MOI.term_indices.(f1.quadratic_terms) == MOI.term_indices.(f2.quadratic_terms)) &&
         (MOI.constant(f1) ≈ MOI.constant(f2)) &&
         (MOI.coefficient.(f1.affine_terms) ≈ MOI.coefficient.(f2.affine_terms)) &&
         (MOI.coefficient.(f1.quadratic_terms) ≈ MOI.coefficient.(f2.quadratic_terms)))
    end
    function test_canonicalization(f::T, expected::T) where {T <: Union{MOI.ScalarAffineFunction, MOI.VectorAffineFunction, MOI.ScalarQuadraticFunction, MOI.VectorQuadraticFunction}}
        @test MOIU.is_canonical(expected)
        g = @inferred(MOIU.canonical(f))
        @test isapprox_ordered(g, expected)
        @test MOIU.is_canonical(g)
        @test MOIU.is_canonical(expected)
        @test isapprox_ordered(MOIU.canonical(g), g)
        @test MOIU.canonical(g) !== g
        @test @allocated(MOIU.canonicalize!(f)) == 0
    end
    @testset "ScalarAffine" begin
        @test MOIU.is_canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])), 1.0))
        @test !MOIU.is_canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([1])), 1.0))
        @test !MOIU.is_canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([2, 1])), 2.0))
        @test !MOIU.is_canonical(MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.0], MOI.VariableIndex.([1, 2])), 2.0))

        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]), 1.5),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]), 1.5),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0], [MOI.VariableIndex(1)]), -2.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]), -2.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.0], [MOI.VariableIndex(2)]), -2.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex[]), -2.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0], [MOI.VariableIndex(2)]),  1.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0], [MOI.VariableIndex(2)]),  1.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0, 4.0], MOI.VariableIndex.([2, 1])), 0.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 3.0], MOI.VariableIndex.([1, 2])), 0.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.1], MOI.VariableIndex.([1, 1])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.1], MOI.VariableIndex.([1])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.5, 2.0], MOI.VariableIndex.([1, 3, 1])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0, 0.5], MOI.VariableIndex.([1, 3])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 0.5, 2.0], MOI.VariableIndex.([1, 3, 1])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([3.0, 0.5], MOI.VariableIndex.([1, 3])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0, 3.0, 4.0, 5.0], MOI.VariableIndex.([1, 3, 1, 3, 2])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([4.0, 5.0, 6.0], MOI.VariableIndex.([1, 2, 3])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], MOI.VariableIndex.([7, 7])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0, 0.5], MOI.VariableIndex.([7, 7, 2])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.5], MOI.VariableIndex.([2])), 5.0),
            )
        test_canonicalization(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -0.5, 0.5, -0.5], MOI.VariableIndex.([7, 7, 2, 7])), 5.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([0.5], MOI.VariableIndex.([2])), 5.0),
            )
    end
    @testset "VectorAffine" begin
        @test MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.(Int[], MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[]))), Float64[]))
        @test !MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1], MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([1]))), [1.0]))
        @test !MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 2], MOI.ScalarAffineTerm.([0.0, 1.0], MOI.VariableIndex.([1, 2]))), [1.0]))
        @test !MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 1], MOI.ScalarAffineTerm.([1.0, 1.0], MOI.VariableIndex.([1, 2]))), [1.0]))
        @test !MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1], MOI.ScalarAffineTerm.([1.0, -1.0], MOI.VariableIndex.([1, 1]))), [1.0]))
        @test !MOIU.is_canonical(MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1], MOI.ScalarAffineTerm.([1.0, 1.0], MOI.VariableIndex.([1, 1]))), [1.0]))

        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.(Int[], MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[]))), Float64[]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.(Int[], MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[]))), Float64[])
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3], MOI.ScalarAffineTerm.([0.0], MOI.VariableIndex.([5]))), [1.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.(Int[], MOI.ScalarAffineTerm.(Float64[], MOI.VariableIndex.(Int[]))), [1.0])
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 2], MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([5, 6]))), [1.0, 2.0, 3.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 3], MOI.ScalarAffineTerm.([2.0, 1.0], MOI.VariableIndex.([6, 5]))), [1.0, 2.0, 3.0])
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 3], MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([6, 5]))), [1.0, 2.0, 3.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 3], MOI.ScalarAffineTerm.([2.0, 1.0], MOI.VariableIndex.([5, 6]))), [1.0, 2.0, 3.0])
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3, 3], MOI.ScalarAffineTerm.([1.0, 2.0], MOI.VariableIndex.([1, 1]))), [1.0, 2.0, 3.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([3], MOI.ScalarAffineTerm.([3.0], MOI.VariableIndex.([1]))), [1.0, 2.0, 3.0])
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1, 1], MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], MOI.VariableIndex.([1, 2, 1]))), [4.0, 5.0, 6.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1], MOI.ScalarAffineTerm.([4.0, 2.0], MOI.VariableIndex.([1, 2]))), [4.0, 5.0, 6.0]),
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1, 1], MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], MOI.VariableIndex.([2, 1, 1]))), [4.0, 5.0, 6.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1, 1], MOI.ScalarAffineTerm.([5.0, 1.0], MOI.VariableIndex.([1, 2]))), [4.0, 5.0, 6.0]),
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 3, 3], MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], MOI.VariableIndex.([1, 1, 1]))), [4.0, 5.0, 6.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 3], MOI.ScalarAffineTerm.([1.0, 5.0], MOI.VariableIndex.([1, 1]))), [4.0, 5.0, 6.0]),
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 3, 3], MOI.ScalarAffineTerm.([1.0, -3.0, 3.0], MOI.VariableIndex.([1, 1, 1]))), [4.0, 5.0, 6.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2], MOI.ScalarAffineTerm.([1.0], MOI.VariableIndex.([1]))), [4.0, 5.0, 6.0]),
            )
        test_canonicalization(
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2, 3, 3, 3], MOI.ScalarAffineTerm.([1.0, 3.0, -1.0, -2.0], MOI.VariableIndex.([1, 1, 1, 1]))), [4.0, 5.0, 6.0]),
            MOI.VectorAffineFunction(MOI.VectorAffineTerm.([2], MOI.ScalarAffineTerm.([1.0], MOI.VariableIndex.([1]))), [4.0, 5.0, 6.0]),
            )
    end
    @testset "ScalarQuadratic" begin
        x = MOI.SingleVariable(MOI.VariableIndex(1))
        y = MOI.SingleVariable(MOI.VariableIndex(2))
        @test MOIU.is_canonical(convert(MOI.ScalarQuadraticFunction{Float64}, 1.0))
        @test !MOIU.is_canonical(1.0x*y + 2.0x*x + 2.0)
        @test !MOIU.is_canonical(1.0x*x + 0.0x*y + 2.0)

        test_canonicalization(
            convert(MOI.ScalarQuadraticFunction{Float64}, 1.5),
            convert(MOI.ScalarQuadraticFunction{Float64}, 1.5)
        )
        test_canonicalization(
            0.0x*y - 2.0,
            convert(MOI.ScalarQuadraticFunction{Float64}, -2.0)
        )
        test_canonicalization(
            3.0x*y + 4.0x*x + 0.0,
            4.0x*x + 3.0x*y + 0.0
        )
        test_canonicalization(
            1.0x*y + 0.1x*y + 5.0,
            1.1x*y + 5.0
        )
    end
    @testset "VectorQuadratic" begin
        x = MOI.SingleVariable(MOI.VariableIndex(1))
        y = MOI.SingleVariable(MOI.VariableIndex(2))
        @test MOIU.is_canonical(MOIU.operate(vcat, Float64, convert(MOI.ScalarQuadraticFunction{Float64}, 1.0)))
        @test !MOIU.is_canonical(MOIU.operate(vcat, Float64, 1.0x*y + 2.0x*x + 2.0))
        @test !MOIU.is_canonical(MOIU.operate(vcat, Float64, 1.0x*x + 0.0x*y + 2.0))

        test_canonicalization(
            MOIU.operate(
                vcat,
                Float64,
                1.5,
                0.0x*y - 2.0,
                3.0x*y + 4.0x*x + 0.0,
                1.0x*y + 0.1x*y + 5.0
            ),
            MOIU.operate(
                vcat,
                Float64,
                1.5,
                -2.0,
                4.0x*x + 3.0x*y + 0.0,
                1.1x*y + 5.0
            )
        )
    end
end

@testset "Vector operate tests" begin
    w = MOI.VariableIndex(0)
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)

    @testset "Vector terms tests" begin
        at = MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(10.0, x))
        qt = MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(6.0, x, y))
        @test MOIU.operate_term(*, 3.0, at) ==
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(30.0, x))
        @test MOIU.operate_term(*, at, at) ==
            MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(200.0, x, x))
        @test MOIU.operate_term(*, 3.0, qt) ==
            MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(18.0, x, y))
        @test MOIU.operate_term(/, at, 2.0) ==
            MOI.VectorAffineTerm(3, MOI.ScalarAffineTerm(5.0, x))
        @test MOIU.operate_term(/, qt, 3.0) ==
            MOI.VectorQuadraticTerm(3, MOI.ScalarQuadraticTerm(2.0, x, y))
    end

    @testset "Vector promote tests" begin
        T = Float64
        for t1 in [MOI.VectorAffineFunction{T}, MOI.VectorOfVariables, Vector{T}]
            for t2 in [MOI.VectorAffineFunction{T}, MOI.VectorOfVariables]
                @test MOIU.promote_operation(+, T, t1, t2) == MOI.VectorAffineFunction{T}
                @test MOIU.promote_operation(-, T, t1, t2) == MOI.VectorAffineFunction{T}
            end
        end
        for t1 in [MOI.VectorQuadraticFunction{T}, MOI.VectorAffineFunction{T},
                MOI.VectorOfVariables, Vector{T}]
            for t2 in [MOI.VectorQuadraticFunction{T}]
                @test MOIU.promote_operation(+, T, t1, t2) == MOI.VectorQuadraticFunction{T}
                @test MOIU.promote_operation(-, T, t1, t2) == MOI.VectorQuadraticFunction{T}
            end
        end
        for t in [MOI.VectorOfVariables, MOI.VectorAffineFunction{T}]
            @test MOIU.promote_operation(-, T, t) == MOI.VectorAffineFunction{T}
        end
        t = MOI.VectorQuadraticFunction{T}
        @test MOIU.promote_operation(-, T, t) == MOI.VectorQuadraticFunction{T}
    end

    @testset "Vector + and - Vector" begin
        α = [1, 2, 3]
        v = MOI.VectorOfVariables([y, w, y])
        g = MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.([3, 1],
                    MOI.ScalarAffineTerm.([5, 2], [y, x])),
                [3, 1, 4])
        f = MOI.VectorQuadraticFunction(
            MOI.VectorAffineTerm.([1, 2, 2],
                MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y])),
            MOI.VectorQuadraticTerm.([1, 1, 2],
                MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y])),
            [7, 3, 4])
        v_plus_g = MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([3, 1, 1, 2, 3],
                        MOI.ScalarAffineTerm.([5, 2, 1, 1, 1], [y, x, y, w, y])),
                    [3, 1, 4])
        g_plus_α = MOI.VectorAffineFunction(
                    MOI.VectorAffineTerm.([3, 1],
                        MOI.ScalarAffineTerm.([5, 2],  [y, x])),
                    [4, 3, 7])
        α_minus_v = MOI.VectorAffineFunction(
                        MOI.VectorAffineTerm.([1, 2, 3],
                            MOI.ScalarAffineTerm.([-1, -1, -1], [y, w, y])),
                        [1, 2, 3])
        v_minus_v_plus_v = MOI.VectorAffineFunction(
                            MOI.VectorAffineTerm.([1, 2, 3, 1, 2, 3, 1, 2, 3],
                                MOI.ScalarAffineTerm.([1, 1, 1, -1, -1, -1, 1, 1, 1],
                                    [y, w, y, y, w, y, y, w, y])),
                            [0, 0, 0])
        f_plus_α = MOI.VectorQuadraticFunction(
                    MOI.VectorAffineTerm.([1, 2, 2],
                        MOI.ScalarAffineTerm.([3, 1, 2], [x, x, y])),
                    MOI.VectorQuadraticTerm.([1, 1, 2],
                        MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y])),
                    [8, 5, 7])
        f_minus_g = MOI.VectorQuadraticFunction(
                        MOI.VectorAffineTerm.([1, 2, 2, 3, 1],
                            MOI.ScalarAffineTerm.([3, 1, 2, -5, -2], [x, x, y, y, x])),
                        MOI.VectorQuadraticTerm.([1, 1, 2],
                            MOI.ScalarQuadraticTerm.([1, 2, 3], [x, y, x], [x, y, y])),
                        [4, 2, 0])
        @test v + g ≈ v_plus_g
        @test g + α ≈ g_plus_α
        @test α + g ≈ g_plus_α
        @test α - v ≈ α_minus_v
        @test MOIU.operate(+, Int, MOIU.operate(-, Int, v, v), v) ≈ v_minus_v_plus_v
        @test f + α ≈ f_plus_α
        @test f - g ≈ f_minus_g
        @test f - f + f - g ≈ f_minus_g
        @test v + f + α - v ≈ f_plus_α
        @test v - f - α - v ≈ - f_plus_α
        @test MOIU.operate!(-, Int, v, f) - v  ≈ - f
        @test (g + v + g + v + f) - (v + g + v + g) ≈ f
        @test v - α ≈ - α_minus_v
        @test g - f ≈ - f_minus_g
    end
    @testset "Vector * and / constant" begin
        v = MOI.VectorOfVariables([y, w, y])
        v2 = MOIU.operate(
            vcat, Float64,
            2.0fy,
            2.0fw,
            2.0fy)
        f = MOIU.operate(
            vcat, Float64,
            2.0fx + 3.0fy + 1.0,
            3.0fx + 2.0fy + 3.0)
        f2 = MOIU.operate(
            vcat, Float64,
            4.0fx + 6.0fy + 2.0,
            6.0fx + 4.0fy + 6.0)
        g = MOIU.operate(
            vcat, Float64,
            7.0fx * fy + 2.0fx + 3.0fy + 1.0,
            6.0fx * fx + 5.0fy * fy + 3.0fx + 2.0fy + 3.0)
        g2 = MOIU.operate(
            vcat, Float64,
            14.0fx * fy + 4.0fx + 6.0fy + 2.0,
            12.0fx * fx + 10.0fy * fy + 6.0fx + 4.0fy + 6.0)
        @testset "$(typeof(a))" for (a, a2) in [(v, v2), (f, f2), (g, g2)]
            @test a * 2.0 ≈ a2
            @test 2.0 * a ≈ a2
            @test a / 0.5 ≈ a2
        end
    end
end


text:
# SparseArrays.dropstored!

function FoodGen1(mean::Vector{<:Real},
        std::Vector{<:Real};
        min_pos=(1,1),
        max_pos=(64,64))
    FoodGen(min_pos, max_pos, MvNormal(mean, diagm(std)))
end

new_food(fg::FoodGen) = round.(Int, clip(rand(fg.dist), fg.smallest_pos, fg.biggest_pos))

new_food(fg::FoodGen, n::Int) =
    rand(fg.dist, n) |> eachcol .|> Tuple .|>
    x->clip(x, fg.smallest_pos, fg.biggest_pos) .|>
    x->round(Int, x)

ext: jl


text:
################################################################################
### JULIA PACKAGE IMPORTS ###
using Plots
using Colors
using DelimitedFiles

################################################################################
### CONTROL FLOW ###

# Source data, if any, also the save location for the next run
Xfilename = "LastX.txt"
Yfilename = "LastY.txt"
GifOutputName = "Dynamics.gif"

# Is warm-start data available?
WarmStart = true

# Add noise to prey data?
ShakeUpX = false

# Set a grid size for simulation
GridSize = 512

# Choose the time length of the simulation
TMAX = 24

# Choose the time point to begin animating from (i.e., animate TSTART:TMAX; note min is 1, not 0)
TSTART = 1

################################################################################
### PARAMETERS ###

μ = 4
β = 5
D1 = 0.001
D2 = 0.2

δ1x = 0 # prey min initialization
δ2x = 1 # prey max initialization
δ1y = 0.3 # predator min initialization
δ2y = 0.4 # predator max initialization

################################################################################
### SET UP GRID ###

# Actual grid size
ShadowGridSize = GridSize+2 # +2 because julia is 1-indexed; the size of the grid including the PBC
GridEnd = GridSize+1 # +1 because julia is 1-indexed; the grid end point excluding the PBC

# Shadow grid with PBC elements
I1S = collect(1:ShadowGridSize) # Ghost elements for grid PBC
I2S = I1S # alias for second dimension

# True grid without PBC elements
I1 = collect(2:GridEnd) # This is one of those times 1-indexing is annoying; the grid points excluding the PBC
I2 = I1 # alias for second dimension

# Number of time points
T = collect(1:TMAX)
TANIM = collect(TSTART:TMAX)

# time points where a progress message is displayed
TPrint = []
TPrintAnim = []
UpdateInterval = 10
for i = 1:UpdateInterval
    push!(TPrint,Int(floor((i/UpdateInterval)*TMAX)))
    push!(TPrintAnim,Int(floor((i/UpdateInterval)*(TMAX-TSTART)+TSTART)))
end

################################################################################
### VARIABLES ###

X = zeros(ShadowGridSize,ShadowGridSize,TMAX)
Y = zeros(ShadowGridSize,ShadowGridSize,TMAX)
∇X = zeros(ShadowGridSize,ShadowGridSize,TMAX)
∇Y = zeros(ShadowGridSize,ShadowGridSize,TMAX)

################################################################################
### SIMULATION FUNCTIONS ###

function x_update(t)
    """
    Prey dynamics
    """
    for i1 in I1
        for i2 in I2
            X[i1,i2,t] = μ*X[i1,i2,t-1]*(1 - X[i1,i2,t-1])*exp(-β*Y[i1,i2,t-1]) + D1*∇X[i1,i2,t-1]
            if X[i1,i2,t] < 0
                X[i1,i2,t] = 0
            end
        end
    end
end

function y_update(t)
    """
    Predator dynamics
    """
    for i1 in I1
        for i2 in I2
            Y[i1,i2,t] = X[i1,i2,t-1]*(1 - exp(-β*Y[i1,i2,t-1])) + D2*∇Y[i1,i2,t-1]
            if Y[i1,i2,t] < 0
                Y[i1,i2,t] = 0
            end
        end
    end
end

function delx_update(t)
    """
    Spatial coupling
    """
    for i1 in I1
        for i2 in I2
            ∇X[i1,i2,t] = X[i1-1,i2,t] + X[i1+1,i2,t] + X[i1,i2-1,t] + X[i1,i2+1,t] - 4*X[i1,i2,t]
        end
    end
end

function dely_update(t)
    """
    Spatial coupling
    """
    for i1 in I1
        for i2 in I2
            ∇Y[i1,i2,t] = Y[i1-1,i2,t] + Y[i1+1,i2,t] + Y[i1,i2-1,t] + Y[i1,i2+1,t] - 4*Y[i1,i2,t]
        end
    end
end

function PBC(t)
    """
    Period Boundary Condition
    """
    for i1 in I1
        X[i1,1,t] = X[i1,GridEnd,t]
        X[i1,ShadowGridSize,t] = X[i1,2,t]
        Y[i1,1,t] = Y[i1,GridEnd,t]
        Y[i1,ShadowGridSize,t] = Y[i1,2,t]
    end
    for i2 in I2
        X[1,i2,t] = X[GridEnd,i2,t]
        X[ShadowGridSize,i2,t] = X[2,i2,t]
        Y[1,i2,t] = Y[GridEnd,i2,t]
        Y[ShadowGridSize,i2,t] = Y[2,i2,t]
    end
end

function XAddNoise(t)
    """
    Add noise to data; i.e., if importing results that are steady-state.
    Noise centred at 0, uniform in +/- δ2x (upper initialization limit)
    After noise addition, values >1 -> 1, value <0 -> 0
    """
    X[:,:,t] += 2*(rand(ShadowGridSize,ShadowGridSize)-0.5*ones(ShadowGridSize,ShadowGridSize))*δ2x
    for i1 in I1
        for i2 in I2
            if X[i1,i2,t] > 1
                X[i1,i2,t] = 1
            elseif X[i1,i2,t] < 0
                X[i1,i2,t] = 0
            end
        end
    end
end

################################################################################
### INITIALIZE FOR t=1 ###
if WarmStart
    # Assign warm start data as initial value
    X[:,:,1] = readdlm(Xfilename, '|', Any)
    Y[:,:,1] = readdlm(Yfilename, '|', Any)
else
    # Assign random values
    X[:,:,1] = rand(ShadowGridSize,ShadowGridSize)*(δ2x-δ1x) + ones(ShadowGridSize,ShadowGridSize)*δ1x
    Y[:,:,1] = rand(ShadowGridSize,ShadowGridSize)*(δ2y-δ1y) + ones(ShadowGridSize,ShadowGridSize)*δ1y
end

# Add noise to prey data; useful to add disturbance to steady-state warm-start results
if ShakeUpX
    XAddNoise(1)
end

################################################################################
### RUN SIMULATION ###
println("Starting simulation...")

# Iterate over time horizon starting from t=1
for t = 2:TMAX
    # Apply boundary conditions to previous time step
    PBC(t-1)

    # Calculate derivative analogues (using PBC)
    delx_update(t-1)
    dely_update(t-1)

    # Update X and Y for period t
    x_update(t)
    y_update(t)

    if t in TPrint
        println("Simulation progress: "*string(Int(floor(100*t/TMAX)))*"%")
    end
end
println("Done simulation!")

# Save last time point for warm start
file = open(Xfilename,"w")
    writedlm(file, X[:,:,TMAX], '|')
close(file)
file = open(Yfilename,"w")
    writedlm(file, Y[:,:,TMAX], '|')
close(file)

################################################################################
### ANIMATE ###
# Start prey animation loop
println("Animating dynamics...")
anim = @animate for t in TANIM
    # Takes the t-th slice and plot
    # for a different colorscheme: color=:halide or, color=:magma
    h1 = heatmap(X[:,:,t], clims=(0,1), xlabel="Prey")
    h2 = heatmap(Y[:,:,t], clims=(0,1), xlabel="Predators")
    plot(h1,h2, layout=(1,2), legend=false, wsize=(1200,600))

    # Add the number to the heatmap
    annotate!(0,0,string(lpad(t - (TSTART-1),ndigits(TMAX),"0")))

    # Progress message:
    if t in TPrintAnim
        println("Animation progress: "*string(Int(floor(100*(t-TSTART)/(TMAX-(TSTART-1)))))*"%")
    end
end
gif(anim, GifOutputName, fps = 7)
println("Done animation!")

max_stars_repo_name: Tominapa/WeekendProjects

max_issues_repo_name: Tominapa/WeekendProjects

max_forks_repo_name: Tominapa/WeekendProjects


text:
module typedb
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  module protocol
    const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
    include("concept_pb.jl")
    include("answer_pb.jl")
    include("core_database_pb.jl")
    include("options_pb.jl")
    include("session_pb.jl")
    include("logic_pb.jl")
    include("query_pb.jl")
    include("transaction_pb.jl")
    include("core_service_pb.jl")
  end
end


text:
using Documenter, HybridSystems

makedocs(
    strict = true,
    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
    html_prettyurls = get(ENV, "CI", nothing) == "true",
    modules = [HybridSystems],
    sitename = "HybridSystems.jl",
    pages = [
        "Home" => "index.md",
        "Library" => Any[
        "Types" => "lib/types.md",
        "Methods" => "lib/methods.md"]
    ]
)

deploydocs(
    repo = "github.com/blegat/HybridSystems.jl.git",
)


text:
#=
using Pkg

tmp_packages = ["Kronecker", "Documenter"]

push!(LOAD_PATH,"../src/")

Pkg.activate(".")

Pkg.add.(tmp_packages) # IMPORTANT
=#

using Documenter, Kronecker, LinearAlgebra

makedocs(
    sitename = "Kronecker.jl",
    authors = "Michiel Stock",
    format = Documenter.HTML(),
    modules = [Kronecker],
    pages = Any[
        "Basic use" => "man/basic.md",
        "Types" => "man/types.md",
        "Multiplication" => "man/multiplication.md",
        "Factorization methods" => "man/factorization.md",
        "Indexed Kronecker products" => "man/indexed.md",
        "Kronecker sums" => "man/kroneckersums.md",
        "Kronecker powers and graphs" => "man/kroneckerpowers.md"
    ]
)

deploydocs(
    repo = "github.com/MichielStock/Kronecker.jl.git",
)


text:
using PowerDynamics
using NetworkDynamics
using LaTeXStrings
include("plotting.jl")
include("ThirdOrderEq.jl")
include("ThirdOrderEqGovernorIEEEG1.jl")
include("ThirdOrderEqExciterIEEEDC1A.jl")

line_list = []
append!(line_list, [StaticLine(from = 1, to = 2, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 3, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 3, to = 4, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 2, to = 4, Y = -1im / 0.02)])
append!(line_list, [StaticLine(from = 1, to = 4, Y = -1im / 0.02)])

node_list = []
append!(node_list, [SlackAlgebraic(U = 1, Y_n = 1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = 0.6, D = 0.01, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
append!(node_list, [ThirdOrderEq(H = 5.0, P = -0.6, D = 0.01, Ω = 50, E_f = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1)])
append!(node_list, [ThirdOrderEqExciterIEEED1A(S_E_max = 0.86, S_E_tq = 0.5, V_R_max = 7.6,H = 5.0, P = 1.0, D = 0.01,  Ω = 50,U = 0.2, X_d_dash = 0.111, T_d_dash = 0.1, X_q_dash = 0.103, X_d = 0.1, K_e = 0.2 , K_f = 0.4, K_a = 0.5, U_ref = 1.2, U_ref2 = 0.9, U_rmax = 1.75, U_rmin = 0.0, T_e = 0.8, T_f = 0.3, T_a = 0.3)])

powergrid = PowerGrid(node_list, line_list)

operationpoint = find_operationpoint(powergrid)

result = simulate(
    Perturbation(3, :ω, Inc(0.00001)),
    powergrid,
    operationpoint,
    timespan = (0.0, 50),
)

plot_res(result, powergrid, 3)

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
function GridAxis(arg0::Vector{jdouble}, arg1::jint)
    return GridAxis((Vector{jdouble}, jint), arg0, arg1)
end

function get_n(obj::GridAxis)
    return jcall(obj, "getN", jint, ())
end

function interpolation_index(obj::GridAxis, arg0::jdouble)
    return jcall(obj, "interpolationIndex", jint, (jdouble,), arg0)
end

function node(obj::GridAxis, arg0::jint)
    return jcall(obj, "node", jdouble, (jint,), arg0)
end

function size(obj::GridAxis)
    return jcall(obj, "size", jint, ())
end


text:
#
# implementation greatly inspired from: https://github.com/geodynamics/seismic_cpml/blob/master

# # Credits
# Author: Pawan Bharadwaj
#         (bharadwaj.pawan@gmail.com)
#
# * original code in FORTRAN90: March 2013
# * modified: 11 Sept 2013
# * major update: 25 July 2014
# * code optimization with help from Jan Thorbecke: Dec 2015
# * rewritten in Julia: June 2017
# * added parrallelization over supersources in Julia: July 2017
# * efficient parrallelization using distributed arrays: Sept 2017
# * optimized memory allocation: Oct 2017

#

# union of old stuff without using ParallelStencils (should eventually remove this)
FdtdOld = Union{FdtdAcoustic,FdtdAcousticVisco,FdtdAcousticBorn}

# global const nlayer_rand = 0

include("types.jl")
include("attenuation.jl")

#
#As forward modeling method, the
#finite-difference method is employed.
#It uses a discrete version of the two-dimensional isotropic acoustic wave equation.
#
#```math
#\pp[\tzero] - \pp[\tmo] = \dt \mB \left({\partial_x\vx}[\tmh]
# + \partial_z \vz[\tmh]  + \dt\sum_{0}^{\tmo}\sfo\right)
# ```
# ```math
#\pp[\tpo] - \pp[\tzero] = \dt \mB \left(\partial_x \vx[\tph]
# + {\partial_z \vz}[\tph]  + \dt\sum_{0}^{\tzero}\sfo\right)
# ```

#attenuation related
#Rmemory::StackedArray2DVector{Array{Float64,3}} # memory variable for attenuation, see Carcione et al (1988), GJ
#Rmemoryp::StackedArray2DVector{Array{Float64,3}} # at previous time step

"""
```julia
pa = SeisForwExpt(attrib_mod; ageom, srcwav, medium, tgrid);
```
Method to create an instance of `SeisForwExpt`.
The output of this method can be used as an input to the in-place method `update!`, to actually perform a
finite-difference modeling.

# Arguments
* `attrib_mod` : attribute to choose the type of modeling. Choose from
  * `=FdtdAcoustic()` for acoustic modeling
  * `=FdtdElastic()` for solving the isotropic elastic wave equation
  * `=FdtdAcousticBorn()` for Born modeling

# Keyword Arguments
* `medium` : medium parameters
* `tgrid` : modeling time grid, maximum time in `tgrid` should be greater than or equal to maximum source time, same sampling interval as in `srcwav`
* `ageom` :  acquisition
* `srcwav` : source wavelets

# Optional Keyword Arguments

* `sflags=2` : source related flags
  * `=0` inactive sources
  * `=1` sources with injection rate
  * `=2` volume injection sources
  * `=3` sources input after time reversal (use only during backpropagation)
* `rflags=1` : receiver related flags
  * `=0` receivers do not record (or) inactive receivers
  * `=1` receivers are active only for the second propagating wavefield
* `rfields=[:p]` : multi-component receiver flag. Choose `Vector{Symbol}`
  * `=[:p]` record pressure
  * `=[:vx]` record horizontal component of particle velocity
  * `=[:vz]` record vertical component of particle velocity
  * `=[:p, :vx]` record both pressure and velocity
* `tsnaps` : store snapshots at these modeling times (defaults to half time)
  * `=[0.1,0.2,0.3]` record at these instances of tgrid
* `snaps_field::Symbol=:tauxx` : which field to record?
* `zfree` : z-coordinate of a stress-free horizontal layer
* `verbose::Bool=false` : verbose flag
"""
SeisForwExpt(
    attrib_mod::Union{FdtdAcoustic,FdtdElastic,FdtdAcousticBorn,FdtdAcousticVisco},
    args1...;
    args2...,
) = PFdtd(attrib_mod, args1...; args2...)

"""
Primary method to generate Expt variable when FdtdAcoustic() and FdtdAcousticBorn() are used.

# Some internal arguments
* `jobname::String` : name
* `npw::Int64=1` : number of independently propagating wavefields in `medium`
* `backprop_flag::Bool=Int64` : save final state variables and the boundary conditions for later use
  * `=1` save boundary and final values in `boundary`
  * `=-1` use stored values in `boundary` for back propagation
* `gmodel_flag=false` : flag that is used to output gradient; there should be atleast two propagating wavefields in order to do so: 1) forward wavefield and 2) adjoint wavefield
* `illum_flag::Bool=false` : flag to output wavefield energy or source illumination; it can be used as preconditioner during inversion
* `pml_edges::Vector{Symbol}=[:zmin, :zmax, :xmax, :xmin]` : use absorbing PML boundary conditions or not
  * `=[:zmin, :zmax]` apply PML conditions only at the top and bottom of the medium
  * `=[:zmax, :xmax, :xmin]` top is reflecting
"""
function PFdtd(
    attrib_mod;
    jobname::Symbol = :forward_propagation,
    npw::Int64 = 1,
    medium::Medium = nothing,
    pml_edges::Vector{Symbol} = [:zmin, :zmax, :ymin, :ymax, :xmax, :xmin],
    tgrid::StepRangeLen = nothing,
    ageom::Union{AGeom,Vector{AGeom}} = nothing,
    srcwav::Union{SrcWav,Vector{SrcWav}} = nothing,
    sflags::Union{Int,Vector{Int}} = fill(2, npw),
    rflags::Union{Int,Vector{Int}} = fill(1, npw),
    rfields::Vector{Symbol} = [:vz],
    zfree = nothing,
    backprop_flag::Int64 = 0,
    gmodel_flag::Bool = false,
    illum_flag::Bool = false,
    tsnaps::AbstractVector{Float64} = fill(0.5 * (tgrid[end] + tgrid[1]), 1),
    snaps_field = nothing,
    verbose::Bool = false,
    nworker = nothing,
)
    N = ndims(medium)

    # convert to vectors
    if (typeof(ageom) == AGeom)
        ageom = [ageom]
    end
    if (typeof(srcwav) == SrcWav)
        srcwav = [srcwav]
    end
    if (typeof(sflags) == Int)
        sflags = [sflags]
    end
    if (typeof(rflags) == Int)
        rflags = [rflags]
    end

    # check if fields input are meaningful for given attrib_mod and ndims
    @assert N == _fd.ndims "Cannot initiate SeisForwExpt due to ndims inconsistency with @init_parallel_stencil"
    !(snaps_field === nothing) && @assert snaps_field ∈ Fields(attrib_mod, ndims = N)
    foreach(rf -> @assert(rf ∈ Fields(attrib_mod; ndims = N)), rfields)
    foreach(
        s -> foreach(
            ss -> foreach(f -> @assert(f ∈ Fields(attrib_mod, ndims = N)), names(ss.d)[1]),
            s,
        ),
        srcwav,
    )

    # check sizes and errors based on input
    #(length(TDout) ≠ length(findn(rflags))) && error("TDout dimension")
    (length(ageom) ≠ npw) && error("ageom dimension")
    (length(srcwav) ≠ npw) && error("srcwav dimension")
    (length(sflags) ≠ npw) && error("sflags dimension")
    (length(rflags) ≠ npw) && error("rflags dimension")
    (maximum(tgrid) < maximum(srcwav[1][1].grid)) &&
        error("modeling time is less than source time")
    #(any([getfield(TDout[ip],:tgrid).δx < tgridmod.δx for ip=1:length(TDout)])) && error("output time grid sampling finer than modeling")
    #any([maximum(getfield(TDout[ip],:tgrid).x) > maximum(tgridmod) for ip=1:length(TDout)]) && error("output time > modeling time")

    # necessary that nss and fields should be same for all nprop
    nss = length(ageom[1])
    sfields = [names(srcwav[ipw][1].d)[1] for ipw = 1:npw]

    # all the propagating wavefields should have same supersources? check that?
    (fill(nss, npw) != [length(ageom[ip]) for ip = 1:npw]) && error("different supersources")

    # check if all sources are receivers are inside medium
    any(.![(ageom[ip] ∈ medium.mgrid) for ip = 1:npw]) && error("sources or receivers not inside medium")

    all([issimilar(ageom[ip], srcwav[ip]) for ip = 1:npw]) || error("ageom and srcwav mismatch")

    #(verbose) &&	println(string("\t> number of super sources:\t",nss))

    # indices where the fields are stored for backpropagation (used for RTM and FWI)
    bindices = get_boundary_indices(medium.mgrid, attrib_mod)

    # extend mediums in the PML layers
    exmedium = padarray(medium, _fd.npextend, pml_edges)
    if (typeof(attrib_mod) == FdtdAcousticBorn)
        (npw ≠ 2) && error("born modeling needs npw=2")
    end

    mod = NamedArray(
        [
            Data.Array(zeros(length.(exmedium.mgrid)...)) for
            name in get_medium_names(attrib_mod)
        ],
        get_medium_names(attrib_mod),
    )
    δmod = NamedArray(
        [zeros(length.(exmedium.mgrid)...) for name in get_medium_names(attrib_mod)],
        get_medium_names(attrib_mod),
    )

    δmodall = zeros(2 * prod(length.(medium.mgrid)))

    # PML
    pml = get_pml(attrib_mod, exmedium.mgrid)

    # initialize gradient arrays
    gradient = zeros(2 * prod(length.(medium.mgrid)))

    # shared arrays required to reduce all the gradient from individual workers
    grad_mod = NamedArray(
        [
            SharedArray{Float64}(zeros(length.(exmedium.mgrid)...)) for
            name in get_medium_names(attrib_mod)
        ],
        get_medium_names(attrib_mod),
    )

    illum_stack = SharedArray{Float64}(zeros(length.(medium.mgrid)...))

    if (!(snaps_field === nothing))
        @assert eval(snaps_field) <: Fields "invalid snaps field"
        itsnaps = NamedArray(
            [argmin(abs.(tgrid .- tsnaps[i])) for i = 1:length(tsnaps)],
            string.(tsnaps),
        )
    else
        snaps_field = :Nothing
        itsnaps = NamedArray([0], ["0"])
    end

    nrmat = [ageom[ipw][iss].nr for ipw = 1:npw, iss = 1:nss]
    datamat = SharedArray{Float64}(length(tgrid), maximum(nrmat), nss)
    data = [Records(tgrid, ageom[ip], rfields) for ip = 1:length(findall(!iszero, rflags))]

    # default is all prpagating wavefields are active
    activepw = [ipw for ipw = 1:npw]

    # dont need visco parameters, initialize with proper sizes later
    nsls = Int32(0)
    memcoeff1 = zeros(1, 1, 1)
    memcoeff2 = zeros(1, 1, 1)

    fc = get_fc(exmedium, tgrid)
    ic = get_ic(exmedium, tgrid, nsls, npw)

    # check stress-free surface condition
    # index of the free surface on the exmedium
    izfree = (zfree === nothing) ? 0 : argmin(abs.(exmedium.mgrid[1] .- zfree))
    ic = vcat(ic, NamedArray([izfree], ([:izfree],)))

    pac = P_common(
        jobname,
        attrib_mod,
        activepw,
        exmedium,
        medium,
        ageom,
        srcwav,
        pml_edges,
        sfields,
        sflags,
        rfields,
        rflags,
        fc,
        ic,
        pml,
        mod,
        NamedArray([memcoeff1, memcoeff2], ([:memcoeff1, :memcoeff2],)),
        δmod,
        δmodall,
        gradient,
        grad_mod,
        illum_flag,
        illum_stack,
        backprop_flag,
        snaps_field,
        itsnaps,
        gmodel_flag,
        bindices,
        datamat,
        data,
        verbose,
    )

    # update mod arrays in pac
    update!(pac, medium)

    # dividing the supersources to workers
    if (nworker === nothing)
        nworker = min(nss, Distributed.nworkers())
    end
    work = Distributed.workers()[1:nworker]
    ssi = [round(Int, s) for s in range(0, stop = nss, length = nworker + 1)]
    sschunks = Array{UnitRange{Int64}}(undef, nworker)
    for ib = 1:nworker
        sschunks[ib] = ssi[ib]+1:ssi[ib+1]
    end

    # a distributed array of P_x_worker --- note that the parameters for each super source are efficiently distributed here
    papa = ddata(
        T = _fd.use_gpu ?
            Vector{P_x_worker_x_pw{N,CUDA.CuArray{_fd.datatype,N,CUDA.Mem.DeviceBuffer}}} : Vector{P_x_worker_x_pw{N,Array{_fd.datatype,N}}},
        init = I -> Vector{P_x_worker_x_pw}(sschunks[I...][1], pac),
        pids = work,
    )

    pa = PFdtd(sschunks, papa, pac)

    # update source wavelets
    update!(pa, pa.c.srcwav, pa.c.sflags, verbose = true)

    # update pml coefficient, now that freqpeak is scanned from the source wavelets
    update_pml!(pa.c)

    check_stability(pa, verbose, H = div(20, _fd.order))

    # update viscoelastic/ viscoacoustic parameters here
    # if(typeof(attrib_mod)==FdtdAcousticVisco)
    # 	nsls=Int32(3)
    # 	exmedium.ic=vcat(exmedium.ic,NamedArray([nsls], ([:nsls],)))
    # 	exmedium.fc=vcat(exmedium.fc,NamedArray(2*pi .* [pa.c.fc[:freqmin], pa.c.fc[:freqmax]], ([:freqmin,:freqmax],)))
    # 	memcoeff1, memcoeff2=get_memcoeff(exmedium)
    # else

    return pa
end

"""
Get some integer constants to store them in pac, e.g., model sizes. ``
The idea is to use them later for a cleaner code.
"""
function get_ic(medium, tgrid, nsls, npw)
    N = ndims(medium)
    return NamedArray(
        vcat(length.(medium.mgrid), [length(tgrid), nsls, npw]),
        vcat(dim_names(N, "n"), [:nt, :nsls, :npw]),
    )
end

"""
Get some float constants, and store then in pac, e.g, spatial sampling.
The idea is to use them later inside the loops for faster modelling.
"""
function get_fc(medium, tgrid)
    N = ndims(medium)
    ds = step.(medium.mgrid)
    # denominator depending on _fd.order
    dsI = inv.(ds)
    if (_fd.order == 4)
        dsI = inv.(ds .* 24.0)
    elseif (_fd.order == 6)
        dsI = inv.(ds .* 1920.0)
    elseif (_fd.order == 8)
        dsI = inv.(ds .* 107520.0)
    end
    dt = step(tgrid)
    dtI = inv(dt)

    return NamedArray(
        vcat([dt, dtI], ds, dsI),
        vcat([:dt, :dtI], dim_names(N, "d"), dim_names(N, "d", "I")),
    )
end

function get_medium_names(::FdtdAcoustic)
    "bulk modulus, and inverse density values"
    return [:KI, :K, :rhoI]
end

function get_medium_names(::FdtdElastic)
    "bulk modulus, and density values on vx and vz stagerred grids"
    return [:lambda, :M, :mu, :rho]
end

# 	# memory fields for all derivatives
# 	memory_pml=NamedArray([NamedArray(fill(zeros(n...), length(fields[d])), fields[d]) for d in dnames], dnames)

"""
Create field arrays for each worker.
Each worker performs the modeling of supersources in `sschunks`.
The parameters common to all workers are stored in `pac`.
"""
function P_x_worker_x_pw(ipw, sschunks::UnitRange{Int64}, pac::P_common{T,N}) where {T,N}
    n = length.(pac.exmedium.mgrid)

    born_svalue_stack = zeros(n...)

    fields = Fields(pac.attrib_mod)
    w1 = NamedArray(
        [NamedArray([zeros(eval(f)(), T(), n...) for f in fields], Symbol.(fields))],
        [:t],
    )

    # temp field arrays copied on the CPU before performing scalar operations
    wr = NamedArray(
        [
            _fd.use_gpu ? Array(zeros(eval(f)(), T(), n...)) :
            zeros(Data.Number, fill(1, N)...) for f in pac.rfields
        ],
        Symbol.(pac.rfields),
    )

    # dummy (use for viscoelastic modeling later)
    w2 = NamedArray(
        [NamedArray([@zeros(fill(1, N)...) for i in [:r]], ([:r],)) for i = 1:2],
        ([:t, :tp],),
    )

    # memory fields for all derivatives
    dfields = Fields(pac.attrib_mod, "d")
    memory_pml = NamedArray(
        [zeros(eval(f)(), T(), n..., pml = true) for f in dfields],
        Symbol.(dfields),
    )

    ss = [P_x_worker_x_pw_x_ss(ipw, iss, pac) for (issp, iss) in enumerate(sschunks)]

    return P_x_worker_x_pw(ss, w1, wr, w2, memory_pml, born_svalue_stack)
end

function Vector{P_x_worker_x_pw}(sschunks::UnitRange{Int64}, pac::P_common)
    return [P_x_worker_x_pw(ipw, sschunks, pac) for ipw = 1:pac.ic[:npw]]
end

# function P_x_worker_x_pw_x_ss(ipw, iss::Int64, pac::T) where T<: P_common{<:FdtdOld}

# 	rfields=pac.rfields
# 	sfields=pac.sfields
# 	nt=pac.ic[:nt]
# 	nx=pac.ic[:nx]; nz=pac.ic[:nz]
# 	nzd,nxd=length.(pac.medium.mgrid)
# 	ageom=pac.ageom
# 	srcwav=pac.srcwav
# 	sflags=pac.sflags
# 	mesh_x, mesh_z = pac.exmedium.mgrid[2], pac.exmedium.mgrid[1]

# 	# records_output, distributed array among different procs
# 	records = NamedArray([zeros(nt,pac.ageom[ipw][iss].nr) for i in 1:length(rfields)], (rfields,))

# 	# gradient outputs
# 	grad_modKI = zeros(nz, nx)
# 	grad_modrhovxI = zeros(nz, nx)
# 	grad_modrhovzI = zeros(nz, nx)

# 	# saving illum
# 	illum =  (pac.illum_flag) ? zeros(nz, nx) : zeros(1,1)

# 	snaps = (pac.snaps_flag) ? zeros(nzd,nxd,length(pac.itsnaps)) : zeros(1,1,1)

# 	# source wavelets
# 	wavelets = [NamedArray([zeros(pac.ageom[ipw][iss].ns) for i in 1:length(sfields[ipw])],(sfields[ipw],)) for it in 1:nt]
# 	fill_wavelets!(ipw, iss, wavelets, srcwav, sflags)

# 	# storing boundary values for back propagation
# 	nz1, nx1=length.(pac.medium.mgrid)
# 	if(pac.backprop_flag ≠ 0)
# 		boundary=[zeros(3,nx1+6,nt),
# 		  zeros(nz1+6,3,nt),
# 		  zeros(3,nx1+6,nt),
# 		  zeros(nz1+6,3,nt),
# 		  zeros(nz1+2*_fd.npml,nx1+2*_fd.npml,3)]
# 	else
# 		boundary=[zeros(1,1,1) for ii in 1:5]
# 	end

# 	coords=[:x1,:x2,:z1,:z2]

# 	is=NamedArray([zeros(Int64,ageom[ipw][iss].ns) for i in coords], (coords,))
# 	# source_spray_weights per supersource
# 	ssprayw = zeros(4,ageom[ipw][iss].ns)
# 	# denomsI = zeros(ageom[ipw][iss].ns)

# 	# receiver interpolation weights per sequential source
# 	rinterpolatew = zeros(4,ageom[ipw][iss].nr)
# 	# denomrI = zeros(ageom[ipw][iss].nr)

# 	pass=P_x_worker_x_pw_x_ss(iss,wavelets,ssprayw,records,rinterpolatew,
# 			  boundary,snaps,illum,# grad_modKI,grad_modrhovxI,grad_modrhovzI)
# 			   NamedArray([grad_modKI,grad_modrhovxI,grad_modrhovzI],([:KI,:rhovxI,:rhovzI],)))

# 	# update acquisition
# 	update!(pass,ipw,iss,ageom[ipw][iss],pac)
# 	return pass
# end

"""
Create modeling parameters for each supersource.
Every worker mediums one or more supersources.
"""
function P_x_worker_x_pw_x_ss(ipw, iss::Int64, pac::P_common{T,N}) where {T,N}
    rfields = pac.rfields
    sfields = pac.sfields
    nt = pac.ic[:nt]
    n = length.(pac.exmedium.mgrid)
    ageom = pac.ageom
    srcwav = pac.srcwav
    sflags = pac.sflags

    # records_output, distributed array among different procs
    records = NamedArray(
        [[Data.Array(zeros(ageom[ipw][iss].nr)) for it = 1:nt] for i = 1:length(rfields)],
        (rfields,),
    )

    # gradient outputs
    grad_modlambda = zeros(1, 1)

    # saving illum
    # illum =  (pac.illum_flag) ? zeros(nz, nx) : zeros(1,1)
    illum = zeros(1, 1)

    snaps = NamedArray(
        [
            (eval(pac.snaps_field) <: Fields) ?
            Array(zeros(eval(pac.snaps_field)(), T(), n...)) :
            zeros(Data.Number, fill(1, N)...) for i = 1:length(pac.itsnaps)
        ],
        names(pac.itsnaps)[1],
    )
    # source wavelets
    wavelets = [
        NamedArray(
            [Data.Array(zeros(ageom[ipw][iss].ns)) for i = 1:length(sfields[ipw])],
            (sfields[ipw],),
        ) for it = 1:nt
    ]
    # fill_wavelets!(ipw, iss, wavelets, srcwav, sflags)
    # storing boundary values for back propagation
    # nz1, nx1=length.(pac.medium.mgrid)
    # if(pac.backprop_flag ≠ 0)
    # 	boundary=[zeros(3,nx1+6,nt),
    # 	  zeros(nz1+6,3,nt),
    # 	  zeros(3,nx1+6,nt),
    # 	  zeros(nz1+6,3,nt),
    # 	  zeros(nz1+2*_fd.npml,nx1+2*_fd.npml,3)]
    # else
    boundary = [zeros(1, 1, 1) for ii = 1:5]
    # end

    # initialize source_spray_weights per supersource and receiver interpolation weights per sequential source
    ssprayw = NamedArray(
        [
            spzeros(
                Data.Number,
                prod(length.(get_mgrid(eval(sf)(), T(), pac.exmedium.mgrid...))),
                ageom[ipw][iss].ns,
            ) for sf in sfields[ipw]
        ],
        sfields[ipw],
    )
    rinterpolatew = NamedArray(
        [
            spzeros(
                Data.Number,
                ageom[ipw][iss].nr,
                prod(length.(get_mgrid(eval(rf)(), T(), pac.exmedium.mgrid...))),
            ) for rf in rfields
        ],
        rfields,
    )

    # transfer to GPUs if
    if (_fd.use_gpu)
        ssprayw=map(x -> CuSparseMatrixCSC(x), ssprayw)
        rinterpolatew=map(x -> CuSparseMatrixCSC(x), rinterpolatew)
    end

    pass = P_x_worker_x_pw_x_ss(
        iss,
        wavelets,
        ssprayw,
        records,
        rinterpolatew,
        boundary,
        snaps,
        illum,# grad_modKI,grad_modrhovxI,grad_modrhovzI)
        NamedArray([grad_modlambda], ([:lambda],)),
    )

    # update acquisition
    update!(pass, ipw, iss, ageom[ipw][iss], pac)
    return pass
end

include("proj_mat.jl")
include("source.jl")
include("receiver.jl")
include("dirichlet.jl")
include("cpml.jl")
include("advance_acou.jl")
include("advance_elastic.jl")
include("rho_projection.jl")
include("gradient.jl")
include("born.jl")
include("boundary.jl")
include("gallery.jl")

# update TDout after forming a vector and resampling
#	ipropout=0;
#	for iprop in 1:pac.ic[:npw]
#		if(pac.rflags[iprop] ≠ 0)
#			ipropout += 1
##			Records.TD_resamp!(pac.data[ipropout], Records.TD_urpos((Array(records[:,:,iprop,:,:])), rfields, tgridmod, ageom[iprop],
##				ageom_urpos[1].nr[1],
##				(ageom_urpos[1].r[:z][1], ageom_urpos[1].r[:x][1])
##				))
#		end
#	end
# return without resampling for testing
#return [Records.TD(reshape(records[1+(iprop-1)*nd : iprop*nd],tgridmod.nx,recv_n,nss),
#		       tgridmod, ageom[1]) for iprop in 1:npw]

function stack_illums!(pac::P_common, pap::Vector{P_x_worker_x_pw{N}}) where {N}
    nx, nz = pac.ic[:nx], pac.ic[:nz]
    illums = pac.illum_stack
    pass = pap[1].ss
    for issp = 1:length(pass)
        gs = pass[issp].illum
        gss = view(gs, _fd.npml+1:nz-_fd.npml, _fd.npml+1:nx-_fd.npml)
        @. illums += gss
    end
end

# Need illumination to estimate the approximate diagonal of Hessian
@inbounds @fastmath function compute_illum!(
    issp::Int64,
    pap::Vector{P_x_worker_x_pw{N}},
) where {N}
    # saving illumination to be used as preconditioner
    p = pap[1].w1[:t][:p]
    illum = pap[1].ss[issp].illum
    for i in eachindex(illum)
        illum[i] += abs2(p[i])
    end
end

"""
Save snapshots for ipw=1
"""
function snaps_save!(itsnap, issp::Int64, pac::P_common, pap)
    p = pap[1].w1[:t][pac.snaps_field]
    snaps = pap[1].ss[issp].snaps[itsnap]
    copyto!(snaps, p)
end

include("stability.jl")
include("updates.jl")
include("getprop.jl")


text:
include("helper.jl")
include("minitest.jl")
Minitest.do_tests() do
  Minitest.include("test_entropy.jl")
end