text:
using MathOptInterface
const MOI = MathOptInterface

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex
const SVF = MOI.SingleVariable

function var_bounds(o::SCIP.Optimizer, vi::VI)
    return MOI.get(o, MOI.ConstraintSet(), CI{SVF,MOI.Interval{Float64}}(vi.value))
end

function chg_bounds(o::SCIP.Optimizer, vi::VI, set::S) where S
    ci = CI{SVF,S}(vi.value)
    MOI.set(o, MOI.ConstraintSet(), ci, set)
    return nothing
end

@testset "Binary variables and explicit bounds" begin
    optimizer = SCIP.Optimizer()

    # Binary variable without explicit bounds.
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test_throws KeyError var_bounds(optimizer, x) == MOI.Interval(0.0, 1.0)

    # Binary variable with [0, 1] bounds.
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.Interval(0.0, 1.0))
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 1.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 1.0)

    # Binary variable with [0, 1] bounds (order should not matter).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b = MOI.add_constraint(optimizer, x, MOI.Interval(0.0, 1.0))
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 1.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 1.0)

    # Binary variable fixed to 0.
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(0.0))
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 0.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 0.0)

    # Binary variable fixed to 0 (different order).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(0.0))
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 0.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(0.0, 0.0)

    # Binary variable fixed to 1.
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(1.0))
    @test var_bounds(optimizer, x) == MOI.Interval(1.0, 1.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(1.0, 1.0)

    # Binary variable fixed to 1 (different order).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(1.0))
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test var_bounds(optimizer, x) == MOI.Interval(1.0, 1.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(1.0, 1.0)

    # Tightened bounds for binary variable
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.Interval(-1.0, 2.0))
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test var_bounds(optimizer, x) == MOI.Interval(-1.0, 2.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(-1.0, 2.0)

    # Tightened bounds for binary variable (different order).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b = MOI.add_constraint(optimizer, x, MOI.Interval(-1.0, 2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(-1.0, 2.0)
    MOI.delete(optimizer, t)
    @test var_bounds(optimizer, x) == MOI.Interval(-1.0, 2.0)

    # No error: binary variable with conflicting bounds (infeasible).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.Interval(2.0, 3.0))
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 3.0)

    # No error: binary variable with conflicting bounds (infeasible, different order).
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b = MOI.add_constraint(optimizer, x, MOI.Interval(2.0, 3.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 3.0)
end

@testset "Bound constraints for a general variable." begin
    optimizer = SCIP.Optimizer()
    inf = SCIP.SCIPinfinity(optimizer)

    # Should work: variable without explicit bounds
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    @test var_bounds(optimizer, x) == MOI.Interval(-inf, inf)

    # Should work: variable with range bounds, but only once!
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.Interval(2.0, 3.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 3.0)
    @test_throws ErrorException MOI.add_constraint(optimizer, x, MOI.Interval(3.0, 4.0))

    # Should work: variable with lower bound, but only once!
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.GreaterThan(2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, inf)
    @test_throws ErrorException MOI.add_constraint(optimizer, x, MOI.GreaterThan(3.0))

    # Should work: variable with lower bound, but only once!
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.LessThan(2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(-inf, 2.0)
    @test_throws ErrorException MOI.add_constraint(optimizer, x, MOI.LessThan(3.0))

    # Should work: fixed variable, but only once!
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 2.0)
    @test_throws ErrorException MOI.add_constraint(optimizer, x, MOI.EqualTo(3.0))

    # Mixed constraint types now allowed (when disjoint)!
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    lb = MOI.add_constraint(optimizer, x, MOI.GreaterThan(2.0))
    ub = MOI.add_constraint(optimizer, x, MOI.LessThan(3.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 3.0)
end

@testset "Changing bounds for variable." begin
    optimizer = SCIP.Optimizer()
    inf = SCIP.SCIPinfinity(optimizer)

    # change interval bounds
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.Interval(2.0, 3.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, 3.0)
    chg_bounds(optimizer, x, MOI.Interval(4.0, 5.0))
    @test var_bounds(optimizer, x) == MOI.Interval(4.0, 5.0)

    # change lower bound
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.GreaterThan(2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, inf)
    chg_bounds(optimizer, x, MOI.GreaterThan(4.0))
    @test var_bounds(optimizer, x) == MOI.Interval(4.0, inf)

    # change upper bound
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.LessThan(3.0))
    @test var_bounds(optimizer, x) == MOI.Interval(-inf, 3.0)
    chg_bounds(optimizer, x, MOI.LessThan(5.0))
    @test var_bounds(optimizer, x) == MOI.Interval(-inf, 5.0)

    # change fixed value
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.EqualTo(2.5))
    @test var_bounds(optimizer, x) == MOI.Interval(2.5, 2.5)
    chg_bounds(optimizer, x, MOI.EqualTo(4.5))
    @test var_bounds(optimizer, x) == MOI.Interval(4.5, 4.5)

    # change mixed
    MOI.empty!(optimizer)
    x = MOI.add_variable(optimizer)
    b = MOI.add_constraint(optimizer, x, MOI.GreaterThan(2.0))
    @test var_bounds(optimizer, x) == MOI.Interval(2.0, inf)
    chg_bounds(optimizer, x, MOI.LessThan(4.0))
    @test var_bounds(optimizer, x) == MOI.Interval(-inf, 4.0)
end

@testset "Second Order Cone Constraint" begin
    # Derived from MOI's problem SOC1
    # max 0x + 1y + 1z
    #  st  x            == 1
    #      x >= ||(y,z)||

    optimizer = SCIP.Optimizer(display_verblevel=0)

    @test MOI.supports_constraint(optimizer, MOI.VectorOfVariables, MOI.SecondOrderCone)

    x, y, z = MOI.add_variables(optimizer, 3)

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [y,z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    ceq = MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.EqualTo(1.0))
    csoc = MOI.add_constraint(optimizer, MOI.VectorOfVariables([x, y, z]),
                              MOI.SecondOrderCone(3))

    @test MOI.get(optimizer, MOI.ConstraintFunction(), csoc) == MOI.VectorOfVariables([x, y, z])
    @test MOI.get(optimizer, MOI.ConstraintSet(), csoc) == MOI.SecondOrderCone(3)

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    atol, rtol = 1e-3, 1e-3
    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ √2 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1/√2 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) ≈ 1/√2 atol=atol rtol=rtol
end

@testset "Second Order Cone Constraint (infeasible)" begin
    # Problem SOC3 - Infeasible
    # min 0
    # s.t. y ≥ 2
    #      x ≤ 1
    #      |y| ≤ x

    optimizer = SCIP.Optimizer(display_verblevel=0)

    x, y = MOI.add_variables(optimizer, 2)

    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.Interval(0.0, 1.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.GreaterThan(2.0))
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([x, y]), MOI.SecondOrderCone(2))

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION
end

@testset "Second Order Cone Constraint (error with unbounded variable)" begin
    optimizer = SCIP.Optimizer()
    x, y = MOI.add_variables(optimizer, 2)
    @test_throws ErrorException MOI.add_constraint(
        optimizer, MOI.VectorOfVariables([x, y]), MOI.SecondOrderCone(2))

    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    MOI.add_constraint(optimizer, MOI.VectorOfVariables([x, y]),
                       MOI.SecondOrderCone(2))
    # no error
end

@testset "SOS1" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)

    x, y, z = MOI.add_variables(optimizer, 3)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.LessThan(1.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.LessThan(1.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.LessThan(1.0))

    c = MOI.add_constraint(optimizer, MOI.VectorOfVariables([x,y,z]), MOI.SOS1([1.0,2.0,3.0]))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,3.0], [x,y,z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(optimizer, MOI.ConstraintFunction(), c) == MOI.VectorOfVariables([x,y,z])
    @test MOI.get(optimizer, MOI.ConstraintSet(), c) == MOI.SOS1([1.0,2.0,3.0])

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    atol, rtol = 1e-6, 1e-6
    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 3.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) ≈ 1.0 atol=atol rtol=rtol
end

@testset "SOS2" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)

    x, y, z = MOI.add_variables(optimizer, 3)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.LessThan(1.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.LessThan(1.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.LessThan(1.0))

    c = MOI.add_constraint(optimizer, MOI.VectorOfVariables([x,y,z]), MOI.SOS2([1.0,2.0,3.0]))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,2.0,3.0], [x,y,z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    @test MOI.get(optimizer, MOI.ConstraintFunction(), c) == MOI.VectorOfVariables([x,y,z])
    @test MOI.get(optimizer, MOI.ConstraintSet(), c) == MOI.SOS2([1.0,2.0,3.0])

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    atol, rtol = 1e-6, 1e-6
    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 5.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) ≈ 1.0 atol=atol rtol=rtol
end

@testset "abspower" begin
    # max.  x1 - x2
    # s.t.  z1 = sign(x1)*abs(x1)^2
    #       z2 = sign(x2)*abs(x2)^3
    #       z1 ≤ 4, z2 ≥ -8

    optimizer = SCIP.Optimizer(display_verblevel=0)

    x1, x2, z1, z2 = MOI.add_variables(optimizer, 4)
    MOI.add_constraint(optimizer, MOI.SingleVariable(z1), MOI.LessThan(4.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(z2), MOI.GreaterThan(-8.0))

    c1 = MOI.add_constraint(optimizer, MOI.VectorOfVariables([x1, z1]),
                            SCIP.AbsolutePowerSet(2.0))
    c2 = MOI.add_constraint(optimizer, MOI.VectorOfVariables([x2, z2]),
                            SCIP.AbsolutePowerSet(3.0))

    @test MOI.get(optimizer, MOI.ConstraintFunction(), c1) == MOI.VectorOfVariables([x1, z1])
    @test MOI.get(optimizer, MOI.ConstraintSet(), c1) == SCIP.AbsolutePowerSet(2.0)
    @test MOI.get(optimizer, MOI.ConstraintFunction(), c2) == MOI.VectorOfVariables([x2, z2])
    @test MOI.get(optimizer, MOI.ConstraintSet(), c2) == SCIP.AbsolutePowerSet(3.0)

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], [x1, x2]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    atol, rtol = 1e-6, 1e-6
    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 4.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x1) ≈ 2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z1) ≈ 4.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x2) ≈ -2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z2) ≈ -8.0 atol=atol rtol=rtol
end

@testset "indicator constraints" begin
    optimizer = SCIP.Optimizer()

    x1 = MOI.add_variable(optimizer)
    x2 = MOI.add_variable(optimizer)
    x3 = MOI.add_variable(optimizer)
    y  = MOI.add_variable(optimizer)
    t = MOI.add_constraint(optimizer, y, MOI.ZeroOne())
    iset = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(1.0))
    ind_func = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x3)),
        ], [0.0, 0.0],
    )

    c = MOI.add_constraint(optimizer, ind_func, iset)
    @test MOI.get(optimizer, MOI.ConstraintFunction(), c) ≈ ind_func
    @test MOI.get(optimizer, MOI.ConstraintSet(), c) == iset

    @test MOI.delete(optimizer, c) === nothing

    # adding incorrect function throws
    ind_func_wrong = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y)),
         MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x1)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x3)),
        ], [0.0, 0.0],
    )
    @test_throws ErrorException MOI.add_constraint(optimizer, ind_func_wrong, iset)
    ind_func_wrong2 = MOI.VectorAffineFunction(
        [MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x2)),
         MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x3)),
        ], [0.0, 0.0],
    )
    @test_throws ErrorException MOI.add_constraint(optimizer, ind_func_wrong2, iset)
end

@testset "deleting variables" begin
    optimizer = SCIP.Optimizer()

    # add variable and constraint
    x = MOI.add_variable(optimizer)
    c = MOI.add_constraint(optimizer, MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0), MOI.EqualTo(0.0))
    @test !MOI.is_empty(optimizer)

    # delete them (constraint first!)
    MOI.delete(optimizer, c)
    MOI.delete(optimizer, x)
    @test MOI.is_empty(optimizer)

    # add variable and constraint (again)
    x = MOI.add_variable(optimizer)
    c = MOI.add_constraint(optimizer, MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0), MOI.EqualTo(0.0))
    @test !MOI.is_empty(optimizer)

    # fail to delete them in wrong order
    @test_throws ErrorException MOI.delete(optimizer, x)
end

@testset "single variable objective" begin
    optimizer = SCIP.Optimizer()

    # Happy Path: add objective and retrieve it.
    x = MOI.add_variable(optimizer)
    obj = MOI.SingleVariable(x)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.SingleVariable}(), obj)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(optimizer, MOI.ObjectiveFunction{MOI.SingleVariable}()) == obj
    @test MOI.get(optimizer, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(optimizer, MOI.ObjectiveFunctionType()) == MOI.ScalarAffineFunction{Float64}

    MOI.empty!(optimizer)

    # Error with type mismatch
    x = MOI.add_variable(optimizer)
    y = MOI.add_variable(optimizer)
    aff_obj = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 3.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), aff_obj)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test_throws ErrorException MOI.get(optimizer, MOI.ObjectiveFunction{MOI.SingleVariable}())
end

@testset "set_parameter" begin
    # bool
    optimizer = SCIP.Optimizer(branching_preferbinary=true)
    @test MOI.get(optimizer, MOI.RawParameter("branching/preferbinary")) == true

    # int
    optimizer = SCIP.Optimizer(conflict_minmaxvars=1)
    @test MOI.get(optimizer, MOI.RawParameter("conflict/minmaxvars")) == 1

    # long int
    optimizer = SCIP.Optimizer(heuristics_alns_maxnodes=2)
    @test MOI.get(optimizer, MOI.RawParameter("heuristics/alns/maxnodes")) == 2

    # real
    optimizer = SCIP.Optimizer(branching_scorefac=0.15)
    @test MOI.get(optimizer, MOI.RawParameter("branching/scorefac")) == 0.15

    # char
    optimizer = SCIP.Optimizer(branching_scorefunc='s')
    @test MOI.get(optimizer, MOI.RawParameter("branching/scorefunc")) == 's'

    # string
    optimizer = SCIP.Optimizer(heuristics_alns_rewardfilename="abc.txt")
    @test MOI.get(optimizer, MOI.RawParameter("heuristics/alns/rewardfilename")) == "abc.txt"

    # invalid
    @test_throws ErrorException SCIP.Optimizer(some_invalid_param_name=true)
end

@testset "use RawParameter" begin
    optimizer = SCIP.Optimizer()

    # bool
    MOI.set(optimizer, MOI.RawParameter("branching/preferbinary"), true)
    @test MOI.get(optimizer, MOI.RawParameter("branching/preferbinary")) == true

    # int
    MOI.set(optimizer, MOI.RawParameter("conflict/minmaxvars"), 1)
    @test MOI.get(optimizer, MOI.RawParameter("conflict/minmaxvars")) == 1

    # long int
    MOI.set(optimizer, MOI.RawParameter("heuristics/alns/maxnodes"), 2)
    @test MOI.get(optimizer, MOI.RawParameter("heuristics/alns/maxnodes")) == 2

    # real
    MOI.set(optimizer, MOI.RawParameter("branching/scorefac"), 0.15)
    @test MOI.get(optimizer, MOI.RawParameter("branching/scorefac")) == 0.15

    # char
    MOI.set(optimizer, MOI.RawParameter("branching/scorefunc"), 's')
    @test MOI.get(optimizer, MOI.RawParameter("branching/scorefunc")) == 's'

    # string
    MOI.set(optimizer, MOI.RawParameter("heuristics/alns/rewardfilename"), "abc.txt")
    @test MOI.get(optimizer, MOI.RawParameter("heuristics/alns/rewardfilename")) == "abc.txt"

    # invalid
    @test_throws ErrorException MOI.set(optimizer, MOI.RawParameter("some/invalid/param/name"), true)
end

@testset "Silent" begin
    optimizer = SCIP.Optimizer()

    @test MOI.supports(optimizer, MOI.Silent())

    # "loud" by default
    @test MOI.get(optimizer, MOI.Silent()) == false
    @test MOI.get(optimizer, MOI.RawParameter("display/verblevel")) == 4

    # make it silent
    MOI.set(optimizer, MOI.Silent(), true)
    @test MOI.get(optimizer, MOI.Silent()) == true
    @test MOI.get(optimizer, MOI.RawParameter("display/verblevel")) == 0

    # but a user can override it
    MOI.set(optimizer, MOI.RawParameter("display/verblevel"), 1)
    @test MOI.get(optimizer, MOI.Silent()) == false
    @test MOI.get(optimizer, MOI.RawParameter("display/verblevel")) == 1
end

@testset "Query results (before/after solve)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.GreaterThan(0.0))

    c = MOI.add_constraint(
        optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
        MOI.LessThan(1.0))

    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # before optimize: can not query most results
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION

    @test_throws ErrorException MOI.get(optimizer, MOI.ObjectiveValue())
    @test_throws ErrorException MOI.get(optimizer, MOI.VariablePrimal(), x)
    @test_throws ErrorException MOI.get(optimizer, MOI.ConstraintPrimal(), c)
    @test_throws ErrorException MOI.get(optimizer, MOI.ObjectiveBound())
    @test_throws ErrorException MOI.get(optimizer, MOI.RelativeGap())
    @test MOI.get(optimizer, MOI.SolveTime()) ≈ 0.0 atol=atol rtol=rtol
    @test_throws ErrorException MOI.get(optimizer, MOI.SimplexIterations())
    @test MOI.get(optimizer, MOI.NodeCount()) == 0

    # after optimize
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.ObjectiveBound()) ≈ 2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.RelativeGap()) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.SolveTime()) < 1.0
    @test MOI.get(optimizer, MOI.SimplexIterations()) <= 3  # conservative
    @test MOI.get(optimizer, MOI.NodeCount()) <= 1          # conservative
end

@testset "Primal start values" begin
    # stop after first feasible solution
    optimizer = SCIP.Optimizer(display_verblevel=0, limits_solutions=1)
    atol, rtol = 1e-6, 1e-6

    # x, y binary
    x, y = MOI.add_variables(optimizer, 2)
    b1 = MOI.add_constraint(optimizer, x, MOI.ZeroOne())
    b2 = MOI.add_constraint(optimizer, y, MOI.ZeroOne())

    # x + y <= 1
    c = MOI.add_constraint(
        optimizer,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 1.0], [x, y]), 0.0),
        MOI.LessThan(1.0))

    # max x + 2y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, 2.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # optimal solution: x = 0, y = 1, value = 2
    # we submit the suboptimal x = 1 as start value
    @test MOI.supports(optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)

    # first set only the value for one variable
    MOI.set(optimizer, MOI.VariablePrimalStart(), x, 1.0)
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), x) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), y) === nothing

    # unset the value
    MOI.set(optimizer, MOI.VariablePrimalStart(), x, nothing)
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), x) === nothing
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), y) === nothing

    # set all values (again)
    MOI.set(optimizer, MOI.VariablePrimalStart(), [x, y], [1.0, nothing])
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), x) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimalStart(), y) === nothing

    # "solve"
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.SOLUTION_LIMIT
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    # we should get back our start values, completed to a solution
    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
end

@testset "No dual solution" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.DualStatus()) == MOI.NO_SOLUTION
end


text:
in = parse.(Int,split(readline("input.txt"), ""))
w = 25
h = 6
l = length(in) ÷ w ÷ h
layers = reshape(in, (w, h, l))

function count_integers(a, c)
    # Alternatively, use list comprehension, as shown in https://github.com/kiranshila/AoC2019/blob/master/8.jl,
    # as [f(x) for x in list] is basically equivalent to map(x -> f(x), list).
    # Also, length(filter()) == count().
    # return [count(x -> x == 0, a[:,:,i]) for i in 1:100]
    return map(layer -> length(filter(x -> x == c, a[:,:,layer])), 1:size(a)[3])
end

zero_counts = count_integers(layers, 0)
layer_with_fewest_zeros = findmin(zero_counts)[2]
ones = length(filter(x -> x == 1, layers[:, :, layer_with_fewest_zeros]))
twos = length(filter(x -> x == 2, layers[:, :, layer_with_fewest_zeros]))
println("Day 8 Part 1: $(ones * twos)")

function color(img, x, y)
    for z in 1:size(img)[3]
        # iterate until we get a color
        if img[x,y,z] == 0
            return ' '
        elseif img[x,y,z] == 1
            return '■'
        end
    end
    println(stderr, "Fully-transparent pixel at ($(x),$(y),$(z))")
    return 0
end

println("Day 8 Part 2: ")
for y in 1:size(layers)[2]
    for x in 1:size(layers)[1]
        print(color(layers, x, y))
    end
    println()
end
println()

# Very satisfying day. This was a much easier puzzle after all that Intcode yesterday!

ext: jl


text:
# HTS221 humidity/temperature sensor.
# See http://www.st.com/resource/en/datasheet/hts221.pdf

const HTS221_t_coef  = Ref((NaN,NaN))
const HTS221_h_coef  = Ref((NaN,NaN))

function HTS221_calibrate()
    setaddr(HTS221_ADDRESS)
    @assert smbus_read(0x0F) == 0xbc

    ## read temp calibration data
    # raw values
    t0_raw = smbus_read_pair(0x3c)
    t1_raw = smbus_read_pair(0x3e)

    # known °C
    u = smbus_read(0x35)
    t0_C = Float64(Int16(u & 0b0011) << 8 | smbus_read(0x32))/8
    t1_C = Float64(Int16(u & 0b1100) << 6 | smbus_read(0x33))/8

    m_t = (t1_C - t0_C)/(Float64(t1_raw) - Float64(t0_raw))
    k_t = t1_C - m_t*t1_raw
    HTS221_t_coef[] = (m_t, k_t)

    ## read humidity calibration data
    h0_raw = smbus_read_pair(0x36)
    h1_raw = smbus_read_pair(0x3a)

    h0_rH = Float64(smbus_read(0x30))/2
    h1_rH = Float64(smbus_read(0x31))/2

    m_h = (h1_rH - h0_rH)/(Float64(h1_raw) - Float64(h0_raw))
    k_h = h1_rH - m_h*h1_raw
    HTS221_h_coef[] = (m_h, k_h)

    return nothing
end

function HTS221_raw_temp()
    setaddr(HTS221_ADDRESS)
    CTRL_REG1 = 0x20
    CTRL_REG2 = 0x21

    smbus_write(CTRL_REG1, 0x00) # power down
    smbus_write(CTRL_REG1, 0x84) # power on, block update
    smbus_write(CTRL_REG2, 0x01) # one-shot aquisition

    while true
        sleep(0.025)
        smbus_read(CTRL_REG2) == 0 && break # check if finished
    end
    t_raw = smbus_read_pair(0x2a)
    smbus_write(CTRL_REG1, 0x00) # power down
    return t_raw
end

function HTS221_raw_humidity()
    setaddr(HTS221_ADDRESS)
    CTRL_REG1 = 0x20
    CTRL_REG2 = 0x21

    smbus_write(CTRL_REG1, 0x00) # power down
    smbus_write(CTRL_REG1, 0x84) # power on, block update
    smbus_write(CTRL_REG2, 0x01) # one-shot aquisition

    while true
        sleep(0.025)
        smbus_read(CTRL_REG2) == 0 && break # check if finished
    end
    h_raw = smbus_read_pair(0x28)
    smbus_write(CTRL_REG1, 0x00) # power down
    return h_raw
end

"""
    HTS221_temperature()

The temperature (in °C) from the HTS221 sensor.
"""
function HTS221_temperature()
    (m_t, k_t) = HTS221_t_coef[]
    t_raw = HTS221_raw_temp()
    m_t*t_raw + k_t
end

"""
    HTS221_humidity()

The relative humidity (as a percentage) from the HTS221 sensor.
"""
function HTS221_humidity()
    (m_h, k_h) = HTS221_h_coef[]
    h_raw = HTS221_raw_humidity()
    m_h*h_raw + k_h
end


text:
module GSMs
export GSM, cascade, cascade!, initialize_gsm_file, append_gsm_data, read_gsm_file,
        Gblock, choose_gblocks, gsm_slab_interface, translate_gsm!, pecgsm, pmcgsm

using LinearAlgebra
using StaticArrays: SA
using OffsetArrays
using Unitful: @u_str, ustrip, unit
using JLD2
using ..Layers: Layer, TEorTM, TE, TM
using ..PGF: mysqrt
using ..Constants: η₀, min_elength
using ..Sheets: Sheet, RWGSheet, find_unique_periods
using FileIO: load

mutable struct GSM{T1 <: AbstractMatrix, T2 <: AbstractMatrix}
    s11::T1
    s12::T2
    s21::T2
    s22::T1
end

function GSM(n1::Int, n2::Int)
    gsm = GSM(zeros(ComplexF64, n1, n1), zeros(ComplexF64, n1, n2),
              zeros(ComplexF64, n2, n1), zeros(ComplexF64, n2, n2))
    gsm.s12[diagind(gsm.s12)] .= one(eltype(gsm.s12))
    gsm.s21[diagind(gsm.s21)] .= one(eltype(gsm.s21))
    gsm
end

function Base.getindex(gsm::GSM, i, j)
    (i,j) == (1,1) && (return gsm.s11)
    (i,j) == (1,2) && (return gsm.s12)
    (i,j) == (2,1) && (return gsm.s21)
    (i,j) == (2,2) && (return gsm.s22)
    throw(BoundsError(gsm, (i,j)))
end

Base.size(::GSM) = (2,2)
function Base.size(::GSM, k)
    (k == 1 || k  == 2) && (return 2)
    k > 2 && (return 1)
    error("arraysize: dimension out of range")
end

function pecgsm(n1,n2)
    gsm = GSM(n1,n2)
    gsm.s11[diagind(gsm.s11)] .= -one(eltype(gsm.s11))
    gsm.s22[diagind(gsm.s22)] .= -one(eltype(gsm.s22))
    gsm.s12 .= zero(eltype(gsm.s12))
    gsm.s21 .= zero(eltype(gsm.s21))
    return gsm
end

function pmcgsm(n1,n2)
    gsm = GSM(n1,n2)
    gsm.s11[diagind(gsm.s11)] .= one(eltype(gsm.s11))
    gsm.s22[diagind(gsm.s22)] .= one(eltype(gsm.s22))
    gsm.s12 .= zero(eltype(gsm.s12))
    gsm.s21 .= zero(eltype(gsm.s21))
    return gsm
end

"""
    Gblock <: Any

     Gblock(rng::UnitRange{Int}, j::Int)

`Gblock` constructor.
The `Gblock` type is used to represent GSMblocks. A GSMblock is a contiguous portion
of the composite FSS structure for which a single GSM (generalized scattering
matrix) is defined and computed. A GSMblock may consist of a single dielectric
interface plane (with or without FSS sheet present), or it may consist
of multiple, adjacent interface planes and the intervening dielectric
layers.  In the latter case, there must be an FSS sheet present at exactly
one of the interface planes within the GSMblock. Note
that interface plane k is located between layers k and k+1, so that
there are N-1 interface planes in a FSS structure consisting of N
dielectric layers.

# Arguments
- `rng`: Specifies the interface planes contained in the `Gblock`. `rng` must
         not be empty, and its initial and final values must be between `1`
         and `N-1` inclusive, where `N` is the number of dielectric layers.
- `j`:  Interface plane number containing the FSS sheet.  If the Gblock does
        not contain a sheet then `j` should be set to `0`.
"""
mutable struct Gblock
    rng::UnitRange{Int}  # Indices of the interface planes in this Gblock.
    j::Int   # Index of interface plane containing FSS sheet or 0 IF no sheet.

    function Gblock(rng::UnitRange{Int}, j::Int)
        length(rng) > 0 || error("rng=$rng must have positive length")
        if j < 0
            error("negative integer j = $j")
        elseif j == 0
            length(rng) ≠ 1 && error("length(rng) must be 1 when j==0")
        else
            j ∈ rng || error("j=$j is not a member of rng=$rng")
        end
        new(rng, j)
    end # function

end # struct

import Base.==
==(g1::Gblock, g2::Gblock) = all((getfield(g1,f)==getfield(g2,f) for
                                                    f in fieldnames(Gblock)))

"""
    cascade(a::GSM, b::GSM) -> c::GSM

Cascade a pair of generalized scattering matrices (GSMs).

## Input arguments

- `a`, `b`:  variables containing the input GSMs that are to be cascaded.  The matrices
must be conformable, i.e., `n2a == n1b`, where `n2a` is the number of modes in Region
2 for GSM `a`, and `n1b` is the number of modes in Region 1 of GSM `b`.
"""
function cascade(a::GSM, b::GSM)::GSM
    n2a = size(a.s22, 2)
    n1b = size(b.s11, 1)
    n1b ≠ n2a && error("Non-conformable arrays")
    n1 =  size(a.s11, 1)
    n2 =  size(b.s22, 2)
    ninner = n1b

    # Equation (3.35) of the theory documentation:
    gprod1 = (I - a.s22 * b.s11) \ a.s21
    s21 = b.s21 * gprod1
    s11 = a.s11 + (a.s12 * (b.s11 * gprod1))
    gprod2 = (I - b.s11 * a.s22) \ b.s12
    s12 = a.s12 * gprod2
    s22 = b.s22 + (b.s21 * (a.s22 * gprod2))
    GSM(s11, s12, s21, s22)
end

"""
    cascade!(a::GSM, layer::Layer, t=NaN)

Cascade a GSM (generalized scattering matrix) on the left with a
dielectric slab on the right. On exit, `a` is modified. `t` is an
layer thickness in layer units.  If not passed or set to `NaN`,
the thickness given in the `layer` variable is used.  It is assumed that
the modes corresponding to the side `2` ports of the GSM are already properly
defined (i.e. normalized) consistent with the dielectric slab electrical properties.
"""
function cascade!(a::GSM, layer::Layer, t=NaN)
    n1 = size(a.s11, 1)
    n2 = size(a.s22, 1)
    n2 ≠ length(layer.γ) && @error "# modes not consistent" n2 length(layer.γ) exception=ErrorException
    # Obtain layer thickness in meters:
    if !isnan(t)
        units_per_meter = ustrip(Float64, unit(layer.user_width), 1u"m")
        d = t / units_per_meter        # Layer thickness in meters.
    else
        d = layer.width
    end
    d == 0 && return nothing
    #  Loop over each of the modes:
    for i in 1:n2
        p = exp(-layer.γ[i] * d)
        a.s12[:,i] .*= p
        a.s21[i,:] .*= p
        a.s22[:,i] .*= p
        a.s22[i,:] .*= p
    end
    return nothing
end

"z-component of the cross product of two 3-vectors"
zdotcross(a,b) = a[1] * b[2] - a[2] * b[1]

"""
    gsm_electric_gblock(layers::AbstractVector{Layer}, s::Integer, k0::Float64) -> (gsm, tlgfvi, vincs)

Calculate the partial GSM (generalized scattering matrix) due to incident fields of a
`Gblock` containing a single electric-type FSS surface.  Also compute the quantities
from Section 6.2 of the theory documentation needed to compute incident and scattered fields.

## Arguments

- `layers`:  Contains the layer parameters for the cascade structure.  Note that the first
   and last layer's thicknesses are not accounted for in this function.  They are assumed
   to be semi-infinite.
- `s`: Interface number (within layers) at which the FSS is located.
- `k0`: Free-space wavenumber in rad/m.

## Return Values

- `gsm::GSM`: A variable containing the portion of the generalized scattering matrix due
    to incident fields.

- `tlgfvi::Matrix{ComplexF64}`: `tlgfvi` is the voltage transmission line Green's function
    (TLGF) due to a unit current. `tlgfvi[q,1]` is the TLGF for mode `q` of Region 1, with
    source at ``z=z_s`` and observation point at ``z=z_1``. `tlgfvi[q,2]` is the TLGF for
    mode `q` of Region N, with source at ``z=z_s`` and observation point at ``z=z_{N-1}``.
- `vincs::Matrix{ComplexF64}`: `vincs` is the transmission line voltage evaluated at ``z = z_s``
    due to a Thevenin voltage source ``V_g = 2`` at either the left (``z=z_1``) or right (``z=z_{N-1}``)
    ends of the equivalent circuit.   `vincs[q,1]` is the voltage for mode `q` with source at ``z=z_1``,
    `vincs[q,2]` is the voltage for mode `q` with the source at ``z=z_{N-1}``.
"""
function gsm_electric_gblock(layers::AbstractVector{Layer}, s::Integer, k0::Float64)
    N = length(layers)
    cosha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    sinha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    tanha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    Z0 = zeros(ComplexF64, N)
    Zleft = zeros(ComplexF64, N)
    Zright = zeros(ComplexF64, N)

    k0sq = k0 * k0
    n1 = length(layers[1].P)  # Number of modes in Region 1.
    n2 = length(layers[N].P)  # Number of modes in Region N.
    gsm = GSM(n1,n2)
    n12max = max(n1, n2)
    vincs = zeros(ComplexF64, n12max, 2)
    tlgfvi = zeros(ComplexF64, n12max, 2)

    # Compute square root of the ratio of the Region 1 and Region n
    # unit cell areas.  This is needed to account for the fact that the
    # two regions may not employ the same lattice to define their modes,
    # hence the mode normalizations may not be compatible.
    root_area1oN = sqrt(abs(
        zdotcross(layers[N].β₁, layers[N].β₂) / zdotcross(layers[1].β₁, layers[1].β₂)))

    #  Loop over the region 1 source modes:
    for q1 in 1:n1
        Z0[1] = 1 / layers[1].Y[q1]  # Obtain region 1 modal impedance.
        Zleft[1] = Z0[1]  # Eq. (6.12a) of the theory documentation
        # Compute mag squared of transverse wavenumber (rad^2/m^2):
        β² = layers[1].β[q1] ⋅ layers[1].β[q1]
        # Loop over each region:
        for i in 2:N
            # Obtain wavenumber (Region i) squared (rad^2/m^2):
            k² = k0sq * layers[i].μᵣ * layers[i].ϵᵣ
            γ = mysqrt(β² - k²) # Atten. constant (neper/m).

            if layers[1].P[q1] == TE
                Z0[i] = (im * (k0*η₀) * layers[i].μᵣ) / γ # TE
            else
                Z0[i] = γ / (im * (k0/η₀) * layers[i].ϵᵣ) # TM
            end
            if i < N
                # Compute remaining quantities for internal layers:
                sinha[i], cosha[i] = sinhcosh(layers[i].width * γ)
                tanha[i] = sinha[i] / cosha[i]
                Zleft[i] = Z0[i] *
                    (Zleft[i-1] + Z0[i]*tanha[i]) / (Z0[i] + Zleft[i-1] * tanha[i])  # Eq. (6.12b).
            end
        end
        # Loop over regions in backwards order to get Zright:
        Zright[N-1] = Z0[N] # Eq. (6.4a)
        for i in N-2:-1:1
            Zright[i] = Z0[i+1] *
                (Zright[i+1] + Z0[i+1]*tanha[i+1]) /
                (Z0[i+1] + Zright[i+1] * tanha[i+1])  # Eq. (6.4b).
        end
        # Compute voltage and current at z = z_1 using Eq. (6.3):
        current = 2  / (Zright[1] + Z0[1])
        voltage = current * Zright[1]
        # Compute reflection term using Eq. (6.2):
        gsm.s11[q1,q1] = voltage - 1

        # Step voltage and current up to junction s using Eq. (6.5):
        for i in 2:s
            voltage = voltage * cosha[i] - Z0[i] * current * sinha[i]
            current = voltage / Zright[i]
        end
        vincs[q1,1] = voltage  # Save incident voltage

        # Locate identical mode in region n:
        q2 = find_mode_index(layers[1].P[q1], layers[1].M[q1], layers[1].N[q1], layers[N])
        if q2 ≠ 0  # Do we need to compute S21 entry?
            # Step voltage and current up to z = z_{N-1} junction using Eq. (6.5):
            for i in s+1:N-1
                voltage = voltage * cosha[i] - Z0[i] * current * sinha[i]
                current = voltage / Zright[i]
            end
            # Compute transmission S-parameter using Eq. (8):
            gsm.s21[q2,q1] = voltage * layers[1].c[q1] / layers[N].c[q2] * root_area1oN
        end

        # Now compute TLGF's:
        current =  -Zright[s] / (Zleft[s] + Zright[s]) # Eqs. (6.19c) and (6.19a)
        voltage = -Zleft[s] * current                  # (6.19a)
        # Step down to z_1 using (6.13):
        for i in s-1:-1:1
            voltage = voltage * cosha[i+1] + Z0[i+1] * current * sinha[i+1]
            current = -voltage / Zleft[i]
        end
        tlgfvi[q1,1] = voltage # Save quantity used in (18a).

    end # Done with Region 1 modes.

    for q2 in 1:n2  #  Loop over the region N source modes:
        Z0[N] = 1 / layers[N].Y[q2]  # Obtain region N modal impedance.
        Zright[N-1] = Z0[N]  # Eq. (6.4a)
        # Compute mag squared of transverse wavenumber (rad^2/m^2):
        β² = layers[N].β[q2] ⋅ layers[N].β[q2]
        for i in N-1:-1:1 # Loop over each region
            k² = k0sq * layers[i].μᵣ * layers[i].ϵᵣ
            γ = mysqrt(β² - k²) # Atten. constant (neper/m).
            if layers[N].P[q2] == TE
                Z0[i] = (im * (k0*η₀) * layers[i].μᵣ) / γ
            else
                Z0[i] = γ / (im * (k0/η₀) * layers[i].ϵᵣ)
            end
            if i > 1
                # Compute remaining quantities for internal layers:
                sinha[i], cosha[i] = sinhcosh(layers[i].width * γ)
                tanha[i] = sinha[i] / cosha[i]
                Zright[i-1] = Z0[i] *
                    (Zright[i] + Z0[i]*tanha[i]) / (Z0[i] + Zright[i] * tanha[i])  # Eq. (6.4b)
            end
        end
        # Loop over regions in forward order to get Zleft:
        Zleft[1] = Z0[1] # Eq. (6.12a)
        for i in 2:N-1
            Zleft[i] = Z0[i] *
                (Zleft[i-1] + Z0[i]*tanha[i]) / (Z0[i] + Zleft[i-1] * tanha[i]) # Eq. (6.12b)
        end
        # Compute voltage and current at z = z_1 using Eq. (6.11):
        current = -2  / (Zleft[N-1] + Z0[N])
        voltage = -current * Zleft[N-1]
        gsm.s22[q2,q2] = voltage - 1 # Compute reflection term using Eq. (6.10)

        # Step voltage and current down to junction s using
        # Equation (6.13):
        for i in N-2:-1:s
            voltage = voltage * cosha[i+1] + Z0[i+1] * current * sinha[i+1]
            current = -voltage / Zleft[i]
        end
        vincs[q2,2] = voltage  # Save incident voltage at FSS plane.

        # Locate identical mode in region 1:
        q1 = find_mode_index(layers[N].P[q2], layers[N].M[q2], layers[N].N[q2], layers[1])
        if q1 ≠ 0
            # Step voltage and current down to z = z_{1} junction using Eq. (6.13):
            for i in s-1:-1:1
                voltage = voltage * cosha[i+1] + Z0[i+1] * current * sinha[i+1]
                current = -voltage / Zleft[i]
            end
            # Compute transmission S-parameter using Eq. (6.8)
            gsm.s12[q1,q2] = voltage * layers[N].c[q2] / layers[1].c[q1] / root_area1oN
        end

        # Now compute TLGF's:
        current =  Zleft[s] / (Zleft[s] + Zright[s]) # Eq. (6.19b)
        voltage = Zright[s] * current                  # Eq. (6.19a)
        # Step up to z_{N-1} using (6.13):
        for i in s+1:N-1
            voltage = voltage * cosha[i] - Z0[i] * current * sinha[i]
            current = voltage / Zright[i]
        end
        tlgfvi[q2,2] = voltage # Save quantity used in (6.18b)
    end  # Done with Region N modes
    return (gsm, tlgfvi, vincs)
end

"""
    gsm_magnetic_gblock(layers::Vector{Layer}, s::Integer, k0::Float64) -> (gsm, tlgfiv, iincs)

Calculate the GSM (generalized scattering matrix) of a `GBLOCK` containing a single magnetic-type
FSS surface.  Also compute the quantities from Section 6.3 of the theory documentation needed to
compute incident and scattered fields.

## Arguments

- `layers`:  Contains the layer parameters for the cascade structure.  Note that the first
 and last layer's thicknesses are not accounted for in this function.  They are assumed to be semi-infinite.
- `s`: Interface number (within layers) at which the FSS is located.
- `k0`: Free-space wavenumber in rad/m.

## Return Values

- `gsm::GSM`: A variable containing the generalized scattering matrix initialized by this function.

- `tlgfiv::Matrix{ComplexF64}` `tlgfiv` is the current transmission line Green's function
 (TLGF) due to a unit voltage. `tlgfiv[q,1]` is the TLGF for mode `q` of Region 1, with source
at ``z=z_s`` and observation point at ``z=z_1``. `tlgfvi[q,2]` is the TLGF for mode `q` of
Region N, with source at ``z=z_s`` and observation point at ``z=z_{N-1}``.

- `iincs::Matrix{ComplexF64}`: `iincs` is the transmission line current evaluated at ``z = z_s``
due to a Norton current source ``I_g = 2`` at either the left (``z=z_1``) or right (``z=z_{N-1}``)
ends of the equivalent circuit.   `iincs[q,1]` is the current for mode `q` with source at ``z=z_1``,
`iincs[q,2]` is the voltage for mode `q` with the source at ``z=z_{N-1}``.
"""
function gsm_magnetic_gblock(layers::AbstractVector{Layer}, s::Integer, k0::Float64)
    N = length(layers)
    cosha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    sinha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    tanha = OffsetArray(zeros(ComplexF64, N-2), 2:N-1)
    Z0 = zeros(ComplexF64, N)
    Zleft = zeros(ComplexF64, N)
    Zright = zeros(ComplexF64, N)

    k0sq = k0 * k0
    n1 = length(layers[1].P)  # Number of modes in Region 1.
    n2 = length(layers[N].P)  # Number of modes in Region N.
    gsm = GSM(n1, n2); gsm.s12 .= zero(eltype(gsm.s12)); gsm.s21 .= zero(eltype(gsm.s21))
    n12max = max(n1,n2)
    iincs = zeros(ComplexF64, n12max, 2)
    tlgfiv = zeros(ComplexF64, n12max, 2)
    # Loop over the region 1 source modes:
    for q1 in 1:n1
        Z0[1] = 1 / layers[1].Y[q1]  # Obtain region 1 modal impedance
        Zleft[1] = Z0[1]  # Semi-infinite line.
        # Compute mag squared of transverse wavenumber (rad^2/m^2):
        β² = layers[1].β[q1] ⋅ layers[1].β[q1]
        # Loop over each region:
        for i in 2:s
            # Obtain wavenumber (Region i) squared (rad^2/m^2):
            k² = k0sq * layers[i].μᵣ * layers[i].ϵᵣ
            γ = mysqrt(β² - k²) # Atten. constant (neper/m).
            if layers[1].P[q1] == TE
                Z0[i] = (im * (k0*η₀) * layers[i].μᵣ) / γ # TE
            else
                Z0[i] = γ / (im * (k0/η₀) * layers[i].ϵᵣ) # TM
            end
            # Compute remaining quantities for internal layers:
            sinha[i], cosha[i] = sinhcosh(layers[i].width * γ)
            tanha[i] = sinha[i] / cosha[i]
            Zleft[i] = Z0[i] *
                (Zleft[i-1] + Z0[i]*tanha[i]) / (Z0[i] + Zleft[i-1] * tanha[i])  # Eq. (6.12b).
        end
        # Loop over regions in backwards order to get Zright:
        Zright[s] = zero(eltype(Zright)) # short-circuited
        #Zright[s-1] = Z0[s] * tanha[s] # Eq. (6.23a)
        for i in s-1:-1:1
            Zright[i] = Z0[i+1] * (Zright[i+1] + Z0[i+1]*tanha[i+1]) /
                (Z0[i+1] + Zright[i+1] * tanha[i+1])  # Eq. (6.23b)
        end

        # Compute voltage and current at z = z_1 using Eq. (6.22):
        current = 2 / (1 + Zright[1] / Z0[1])
        voltage = current * Zright[1]
        # Compute reflection term using Eq. (6.27):
        gsm.s11[q1,q1] = 1 - current

        # Step voltage and current up to junction s using Eq. (6.24):
        for i in 2:s
            current = current * cosha[i] - voltage * sinha[i] / Z0[i]
            voltage = current * Zright[i]
        end
        iincs[q1,1] = current  # Save incident current.

        # Now compute TLGF:
        current =  1 / Zleft[s] # Eq. (6.38a)
        voltage = -one(current)     # Eq. (6.38a)
        # Step down to z_1 using (6.32):
        for i in s-1:-1:1
            current = current * cosha[i+1] + voltage * sinha[i+1] / Z0[i+1]
            voltage = -current * Zleft[i]
        end
        tlgfiv[q1,1] = current # Save quantity used in Eq. (6.37a)

    end  # Done with Region 1 modes

    #  Loop over the region N source modes:
    for q2 in 1:n2
        Z0[N] = 1 / layers[N].Y[q2]  # Obtain region N modal impedance
        Zright[N-1] = Z0[N]  # Semi-infinite line
        # Compute mag squared of transverse wavenumber (rad^2/m^2):
        β² = layers[N].β[q2] ⋅ layers[N].β[q2]
        for i in N-1:-1:s+1 # Loop over each region
            # Obtain wavenumber (Region i) squared (rad^2/m^2):
            k² = k0sq * layers[i].μᵣ * layers[i].ϵᵣ
            γ = mysqrt(β² - k²) # Atten. constant (neper/m).
            if layers[N].P[q2] == TE
                Z0[i] =  (im * (k0*η₀) * layers[i].μᵣ) / γ # TE
            else
                Z0[i] = γ / (im * k0/η₀ * layers[i].ϵᵣ)
            end
            # Compute remaining quantities for internal layers:
            sinha[i], cosha[i] = sinhcosh(layers[i].width * γ)
            tanha[i] = sinha[i] / cosha[i]
            Zright[i-1] = Z0[i] * (Zright[i] + Z0[i]*tanha[i]) /
                (Z0[i] + Zright[i] * tanha[i])  # Eq. (6.23b)
        end

        # Loop over regions in forward order to get Zleft:
        Zleft[s] = zero(eltype(Zleft)) # Short-circuited
        #Zleft(s+1) = Z0(s+1) * tanha(s+1) # Eq. (6.31a)
        for i in s+1:N-1
            Zleft[i] = Z0[i] * (Zleft[i-1] + Z0[i]*tanha[i]) /
                 (Z0[i] + Zleft[i-1] * tanha[i])  # Eq. (6.31b)
        end

        # Compute voltage and current at z = z_{N-1} using Eq. (6.30):
        current = -2  / (1 + Zleft[N-1] / Z0[N])
        voltage = -current * Zleft[N-1]
        # Compute reflection term using Eq. (6.35):
        gsm.s22[q2,q2] = 1 + current
        # Step voltage and current down to junction s using Eq. (6.32):
        for i in N-2:-1:s
            current = current * cosha[i+1] + voltage * sinha[i+1] / Z0[i+1]
            voltage = -current * Zleft[i]
        end
        iincs[q2,2] = current  # Save incident voltage at FSS plane

        # Now compute TLGF's:
        current = 1 / Zright[s] # Eq. (6.38b)
        voltage = one(current) # Eq. (6.38b)
        # Step up to z_{N-1} using (6.24):
        for i in s+1:N-1
            current = current * cosha[i] - voltage * sinha[i] / Z0[i]
            voltage = current * Zright[i]
        end
        tlgfiv[q2,2] = current # Save quantity used in (6.18b)
    end # Done with Region N modes
    return (gsm, tlgfiv, iincs)
end

"""
    gsm_slab_interface(L1::Layer, L2::Layer, k0) -> gsm::GSM

Calculate the GSM (generalized scattering matrix) of a dielectric junction.

## Arguments

- `L1`, `L2`:  Contains the layer parameters layers 1 and 2, resp.
- `k0`: Free-space wavenumber in rad/m.

## Return Values

- `gsm`: A variable containing the generalized scattering matrix computed in accordance
    with Equations (3.23) and (3.24) of the theory documentation.

"""
function gsm_slab_interface(L1::Layer, L2::Layer, k0)::GSM

      n1 = length(L1.P)  # Number of modes in Region 1.
      n2 = length(L2.P)  # Number of modes in Region 2.
      gsm = GSM(n1, n2)
      #  Loop over the region 1 source modes:
      for i1 in 1:n1
          Y1 = L1.Y[i1]  # Obtain source modal admittance.
          # We now require the modal admittance for the corresponding
          # mode in Region 2. This, however, is not guaranteed to exist
          # since we don't necessarily define the same sets of modes in all
          # regions.  Therefore, we search for it on the other side,
          # but if it doesn't exist, we have to compute the modal
          # admittance on the fly.

          # Locate identical mode in region 2:
          i2 = find_mode_index(L1.P[i1], L1.M[i1], L1.N[i1], L2)
          if i2 == 0
              # Obtain wavenumber (Region 2) squared
              k2sq = L2.β[1] ⋅ L2.β[1] - L2.γ[1]^2
              β² = L1.β[i1] ⋅ L1.β[i1] # for desired mode
              γ = mysqrt(β² - k2sq)
              if L1.P[i1] == TE
                  Y2 = γ / (im * (k0*η₀) * L2.μᵣ)
              else
                  Y2 = im * (k0/η₀) * L2.ϵᵣ / γ
              end
          else
              Y2 = L2.Y[i2]
          end

          # Compute reflection and transmission coefficients:
          R = (Y1 - Y2) / (Y2 + Y1) # Eq. (3.24a)
          gsm.s11[i1,i1] = R
          if i2 != 0
              T = 2 * mysqrt(Y2) *  mysqrt(Y1) / (Y2 + Y1) # Eq. (3.24b)
              gsm.s21[i2,i1] = T
          end
      end

      #  Loop over the region 2 source modes:
      for i2 in 1:n2
          Y2 = L2.Y[i2]  # Obtain source modal admittance.
          # We now require the modal admittance for the corresponding
          # mode in Region 1. This, however, is not guaranteed to exist
          # since we don't necessarily define the same sets of modes in all
          # regions.  Therefore, we search for it on the other side,
          # but if it doesn't exist, we have to compute the modal
          # admittance on the fly.

          # Locate identical mode in region 1:
          i1 = find_mode_index(L2.P[i2], L2.M[i2], L2.N[i2], L1)

          if i1 == 0
              # Obtain wavenumber (Region 1) squared
              k1sq = L1.β[1] ⋅ L1.β[1] - L1.γ[1]^2
              β² = L2.β[i2] ⋅ L2.β[i2] # for desired mode
              γ = mysqrt(β² - k1sq)
              if L2.P[i2] == TE
                  Y1 = γ / (im * (k0*η₀) * L1.μᵣ)
              else
                  Y1 = im * (k0/η₀) * L1.ϵᵣ / γ
              end
          else
              Y1 = L1.Y[i1]
          end
          # Compute reflection and transmission coefficients:
          R = (Y1 - Y2) / (Y2 + Y1) # Eq. (3.24a)
          gsm.s22[i2,i2] = -R
          if i1 != 0
              T = 2 * mysqrt(Y2) * mysqrt(Y1) / (Y2 + Y1) # Eq. (3.24b)
              gsm.s12[i1,i2] = T
          end
      end
      return gsm
  end # function

"""
    initialize_gsm_file(fname::AbstractString, layer1::Layer, layer2::Layer, sh1::Sheet, sh2::Sheet)

Initialize a generalized scattering matrix (GSM) file, by appending the GSM-specific
header info, including the unit cell data.  The file is stored in `JLD2` format.

## Arguments

- `fname`: A string containing the file name to open.  The file will be closed after writing.
- `layer1`, `layer2`:  These provide the mode indices for the input and output regions.
- `sh1`, `sh2`: `SHEET` objects from which the unit
         cell info for the input and output regions will be extracted.  These two arguments
    are optional, but if either is present the other must be also.  If they are absent, then
    a square unit cell of dimension 1 meter will be used for both input and output regions.
"""
function initialize_gsm_file(fname::AbstractString, layer1::Layer, layer2::Layer, sh1::Sheet, sh2::Sheet)
    jldopen(fname, "w") do fid
        fid["FileType"] = "GSMFile"
        fid["FileTypeVersion"] = v"0.0.1"
        # Unit cell information:
        fid["s1in"] = sh1.s₁
        fid["s2in"] = sh1.s₂
        fid["unitsin"] = sh1.units
        fid["s1out"] = sh2.s₁
        fid["s2out"] = sh2.s₂
        fid["unitsout"] = sh2.units

        # Modelists:
        fid["ModelistIn"] = (layer1.P, layer1.M, layer1.N)
        fid["ModelistOut"] = (layer2.P, layer2.M, layer2.N)
    end
end # function

function initialize_gsm_file(fname::AbstractString, layer1::Layer, layer2::Layer)
    sh = RWGSheet()
    sh.s₁ = [1.0, 0.0]
    sh.s₂ = [0.0, 1.0]
    sh.units = u"m"
    initialize_gsm_file(layer1, layer2, fname, sh, sh)
    return
end

"""
    append_gsm_data(fname::AbstractString, gname::String, gsm::GSM, l1::Layer, l2::Layer, case::Dict{String,Float64})

Append a block of GSM data to a GSM file for a particular frequency and pair of scan parameters.  It is
assumed that the file has already been initialized by a previous call to `initialize_gsm_file`.

## Arguments

- `fname`: The name of the GSM file to be appended to.  The file should already have been initialized
via a call to `initialize_gsm_file`.
- `gname`: The unique `JLD2` group name to be used in the file for grouping the data associated with this frequency/scan case.
- `gsm`:  The GSM data to be written to the file.
- `l1`, `l2`: `Layer` data for the input and output regions containing the permeability and permittivity information
to be written to the file.
- `case`: A dictionary containing three key-value pairs which together fully define the frequency/scan case.
The keys are "FGHz" and either ("θ" and "ϕ") or ("ψ₁" and "ψ₂").
"""
function append_gsm_data(fname::AbstractString, gname::String, gsm::GSM, l1::Layer, l2::Layer, case::Dict{String,Float64})
    jldopen(fname, "a") do fid
        group = JLD2.Group(fid, gname)
        group["case"] = case
        group["epsmu1"] = (l1.ϵᵣ, l1.μᵣ)
        group["epsmu2"] = (l2.ϵᵣ, l2.μᵣ)
        group["gsm"] = gsm
    end
    return
end

function read_gsm_file(fname::AbstractString)
    dat = load(fname)
    (haskey(dat, "FileType") && dat["FileType"] == "GSMFile") || error("$fname is not a GSM file")
    dat
end

"""
    translate_gsm!(g::GSM, dx::Real, dy::Real, layer1::Layer, layer2::Layer)

Convert `GSM` `g` into one for an identical structure that is translated by
`dx` meters in x and `dy` meters in y. `layer1` and `layer2` define the modes
on each side of the GSM.

## Arguments:

- `g`: On input contains the GSM for a device with zero x and y translation. After
   return, `g` will be modified by applying appropriate phase shifts due to translation,
   according the to the formula presented in John C. Vardaxoglou, "Frequency Selective Surfaces",
   Wiley, 1997, Eq. (7.57).
- `dx`, `dy`:   The translations in the x and y directions in meters.
- `layer1`, `layer2`:  The β components of these variables are used with
    `dx` and `dy` to determine the effect of the translation.
"""
function translate_gsm!(g, dx, dy, layer1, layer2)

    dxy = SA[dx, dy]

    # Check that layers have correct number of modes:
    size(g.s11,1) == length(layer1.P) || error("N1 mismatch")
    size(g.s21,1) == length(layer2.P) || error("N2 mismatch")

    # Compute phase shift constants:
    t1 = [cis(β ⋅ dxy) for β in layer1.β]
    t2 = [cis(β ⋅ dxy) for β in layer2.β]

    # Fix up S11:
    for n in 1:length(layer1.P), m in 1:length(layer1.P)
        g.s11[m,n] *= t1[n] / t1[m]
    end

    # Fix up S12:
    for n in 1:length(layer2.P), m in 1:length(layer1.P)
        g.s12[m,n] *= t2[n] / t1[m]
    end

    # Fix up S21:
    for n in 1:length(layer1.P), m in 1:length(layer2.P)
        g.s21[m,n] *= t1[n] / t2[m]
    end

    # Fix up S22:
    for m in 1:length(layer2.P), n in 1:length(layer2.P)
        g.s22[m,n] *= t2[n] / t2[m]
    end

    return
end

"""
    choose_gblocks(layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int}, k0min::Float64) -> gbl::Vector{Gblock}

Set up gbl, the array of GBLOCKs the defines the basic GSM building
blocks for the FSS structure.

## Input arguments:

- `layers`: A vector of `Layer` objects.
- `sheets`: A vector of `Sheet` objects.
- `junc`: An integer vector containing in `junc[k]` the sheet number present at junction `k` or `0`
  if no sheet is present.
- `k0min`:    The minimum free-space wavenumber (rad/m) for all analysis cases.

## Return value:

- `gbl`:  A vector of `Gblock` elements defining the GBLOCKS used to analyze the complete FSS structure.
"""
function choose_gblocks(layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int},
                                                                     k0min::Float64)::Vector{Gblock}
    nl = length(layers)
    nj = nl-1 # number of dielectric junctions
    ns = length(sheets)

    # Pure radome case is easy:
    if ns == 0 || all(t.style == "NULL" for t in sheets)
        gbl = [junc[i] == 0 ? Gblock(i:i, 0) : Gblock(i:i,i) for i = 1:nj]
        return gbl
    end

    # Compute min. electrical length of each layer:
    elength = [k0min * l.width * sqrt(real(l.ϵᵣ * l.μᵣ)) for l in layers]
    sint = findall(junc .≠ 0) # sint[k] contains dielectric interface number of k'th sheet
    upa = find_unique_periods(junc, sheets)

    # Initialize junction ownerships to self:
    owner = collect(1:nj) # owner[i] is the interface number to which interface i belongs.

    # Assign junction ownership for thin layers adjacent and prior to location of first sheet:
    if sheets[1].style ≠ "NULL"
        for j in sint[1]:-1:2 # layer index
            elength[j] > min_elength && break
            owner[j-1] = sint[1]
        end
    end
    # Assign junction ownership for thin layers adjacent to and after location of last sheet:
    if sheets[end].style ≠ "NULL"
        for j in 1+sint[ns]:nj # layer index
            elength[j] > min_elength && break
            owner[j] = sint[end]
        end
    end

    # Loop over each adjacent pair of junctions stored in sint:
    for i in 1:ns-1
        (j1, j2) = sint[i:i+1]   # Interface locations
        ## to do:  check here if j2-j1==1, and periods are different
        Nwide = count(elength[(1+j1):j2] .> min_elength) # Numb. of wide enough intervening layers
        if Nwide > 0  # step into the region from each end and assign ownership
            for j in j2:-1:j1+1
                (elength[j] > min_elength || sheets[junc[j2]].style == "NULL") && break
                owner[j-1] = j2
            end
            for j in 1+j1:j2
                (elength[j] > min_elength || sheets[junc[j1]].style == "NULL") && break
                owner[j] = j1
            end
        else
            # Assign neighboring junctions for thin layers adjacent to j1 and j2 to j1 and j2
            #
            if sheets[junc[j1]].style == "NULL" && sheets[junc[j2]].style ≠ "NULL"
                owner[j1+1:j2] .= j2
                continue
            elseif sheets[junc[j1]].style ≠ "NULL" && sheets[junc[j2]].style == "NULL"
                owner[j1:j2-1] .= j1
                continue
            elseif sheets[junc[j1]].style == "NULL" && sheets[junc[j2]].style == "NULL"
                continue
            end
            #
            # At this point, we've determined that both sheets are non-null, and closely separated
            # Pick the electrically thickest layer between them as the dividing line.
            # If several layers tie for max thickness, choose the one closest to the middle.
            elmax = maximum(@view elength[j1+1:j2])
            best_layers = findall(elength .== elmax)
            filter!(j -> j1+1 ≤ j ≤ j2, best_layers)
            layermid = (j1 + j2 + 1) ÷ 2
            (_,ibest) = findmin(abs.(best_layers .- layermid))
            layerbest = best_layers[ibest]
            owner[j1:layerbest-1] .= j1
            owner[layerbest:j2] .= j2
        end
    end
    # Count number of GBLOCKS needed:
    ngbl = 1
    for j in 2:nj
        owner[j] ≠ owner[j-1] && (ngbl += 1)
    end
    # and assign them:
    gbl = [Gblock(1:1, 0) for j in 1:ngbl]
    jg = 1 # Initialize Gblock counter
    for j in 1:nj  # Loop over each interface
        if owner[j] == owner[first(gbl[jg].rng)]
            gbl[jg].rng = first(gbl[jg].rng):j # Increment the end pointer.
        else
            jg += 1
            gbl[jg].rng = j:j
            gbl[jg].j = 0
        end
        junc[j] ≠ 0 && (gbl[jg].j = j)
    end

    # Check that every junction occurs exactly once in the list of Gblocks
    vcat((vec(g.rng) for g in gbl)...) == 1:nj || error("Bad Gblock vector: $gbl")
    return gbl
end

"""
    find_mode_index(p::TEorTM, m::Int, n::Int, layer::Layer)

Find the index of the mode (p,m,n).

## Arguments

- `p`, `m`, `n`: Mode indices.
- `layer`: An instance of Layer with modes initialized.

## Return Value

The positive integer index, if found.  Otherwise, zero.
"""
function find_mode_index(p, m::Int, n::Int, layer::Layer)
    for i in 1:length(layer.P)
        if layer.P[i] == p  && layer.M[i] == m && layer.N[i] == n
            return i
        end
    end
    return 0
end

@inline function sinhcosh(x)
    @fastmath begin
        expx = exp(x)
        expmx = 1/expx
        sinhx = 0.5*(expx - expmx)
        coshx = 0.5*(expx + expmx)
    end
    return (sinhx, coshx)
end

end # module


text:
"""
    `CalcVSCoefs(boundary, etas; U_slip, sigma)`

Calculates the vortex sheet RBF coefficients.

# ARGUMENTS
* `boundary::Boundary`                       : Boundary object
* `etas::Array{T}`                           : RBF coefficients for rho
* `U_slip::Union{Nothing, Function}=nothing` : Slip velocity field function: U_slip(X)
* `sigma = 0.2`                              : Gaussian width
where {T<:Real}

# RETURNS
* `alphas::Array{Float64}`                   : Array of RBF coefficients
"""
function CalcVSCoefs(boundary::Boundary,
                     etas::Array{T};
                     U_slip::Union{Nothing, Function}=nothing,
                     sigma=0.2) where {T<:Real}

    # === Extract geometry from boundary ===
    X_body = [[point[1], point[2], 0.0] for point in boundary.bodyPTS]
    n_hats = [[n_hat[1], n_hat[2], 0.0] for n_hat in boundary.nHats]
    t_hats = [[t_hat[1], t_hat[2], 0.0] for t_hat in boundary.tHats]
    NPTS = boundary.NPTS_BODY

    # === Function definitions ===

    # i -> RBF index
    # j -> Integration index
    # k -> Evaluation index

    X_i(i) = X_body[i]                                  # Set_i of points on S
    X_j(j) = X_body[j]                                  # Set_j of points on S
    X_jp1(j) = (j != NPTS) ? X_body[j+1] : X_body[1]    # Point j + 1
    X_k(k) = X_body[k]                                  # Set_k of points on S

    R_kj(k,j) = X_k(k) - X_j(j)                         # X_k - X_j
    R_ji(j,i) = X_j(j) - X_i(i)                         # X_j - X_i
    R_ki(k,i) = X_k(k) - X_i(i)                         # X_k - X_i
    r_kj(k,j) = norm(R_kj(k,j))                         # | X_k - X_j |
    r_ji(j,i) = norm(R_ji(j,i))                         # | X_j - X_i |
    r_ki(k,i) = norm(R_ki(k,i))                         # | X_k - X_i |

    del_S(j) = norm(X_j(j) - X_jp1(j))                  # ΔS_j
    L = sum( [del_S(j) for j in 1:NPTS] )               # Length of surface

    #etas = CalcRhoCoefs(boundary)                       # Coefficients for rho RBF
    rho_1(k) = CalcRho(boundary, etas, X_k(k))          # ρ_1 of Eigendecomposition

    # === Coefficent functions ===
    function Theta_ki(k,i)
        res = 0
        for j in 1:NPTS
            if j != k
                res = res + (1/(2*pi)) * (dot(R_kj(k,j), n_hats[j]) * RBF_gauss(r_ji(j,i), sigma=sigma) * del_S(j))/(r_kj(k,j)^2)
            end
        end
        return res
    end
    #Theta_ki(k,i) = sum([((k != j) ? (1/pi) * (dot(R_kj(k,j), n_hats[j]) * RBF_gauss(r_ji(j,i)) * del_S(j))/(r_kj(k,j)^2) : 0) for j in 1:NPTS] )
    Lambda_ki(k,i) = sum( [RBF_gauss(r_ji(j,i), sigma=sigma) * (rho_1(k))/(L) for j in 1:NPTS] )
    phi_ki(k,i) = RBF_gauss(r_ki(k,i), sigma=sigma)

    # === Calculate coefficient matrix ===
    A = zeros(NPTS,NPTS)
    for k in 1:NPTS
        for i in 1:NPTS
            A[k,i] = phi_ki(k,i) - Theta_ki(k,i) + Lambda_ki(k,i)
        end
    end
    #A = [ [phi_ki(k,i) - Theta_ki(k,i) + Lambda_ki(k,i) for k in 1:NPTS] for i in 1:NPTS]
    #println("A: ",A)

    # === Calculate RHS vector ===
    b = zeros(NPTS)
    for i in 1:NPTS
        b[i] = dot(U_slip(X_i(i)), t_hats[i])
    end
    #b = [dot(U_slip(X_i(i)), t_hats[i]) for i in 1:NPTS]
    #println("b: ",b)

    # === Calculate alpha coefficients ===
    alphas = A\b
    return alphas

end

"""
    `CalcVS(boundary, alpha, X_eval)`

Evaluates the RBF interpolation function: gamma(X) = sum(alpha_i * phi(X - X_i)).

# ARGUMENTS
* `boundary::Boundary` : Boundary object
* `alpha::Array{T}`    : Array containing RBF coefficients
* `X_eval::Array{T}`   : Point to evaluate function at
* `sigma = 0.2`        : Gaussian width
where {T<:Real}

# RETURNS
* `gamma::Float64`     : Evaluated function value
"""
function CalcVS(boundary::Boundary,
                alpha::Array{T},
                X_eval::Array{T};
                sigma=0.2) where {T<:Real}

    # === Unpack geometry ===
    XP = [[point[1], point[2], 0.0] for point in boundary.bodyPTS]
    NPTS = boundary.NPTS_BODY

    # === Summation over all points for gamma ===
    gamma = -sum( [alpha[i] * RBF_gauss(norm(X_eval - XP[i]), sigma=sigma) for i in 1:NPTS ] )
    return gamma

end

# """
#
# """
# function CalcVSVelocityCirlce(boundary::Boundary,
#                               alpha::Array{T},
#                               X_eval::Array{T};
#                               radius=1.0) where {T<:Real}
#
#     # === Unpack geometry ===
#     X_body = [[point[1], point[2], 0.0] for point in boundary.bodyPTS]
#     NPTS = boundary.NPTS_BODY
#
#     # === Function definitions ===
#     X_j(j) = X_body[j]                                       # Set_j of points on S
#     X_jp1(j) = (j != NPTS) ? X_body[j+1] : X_body[1]         # Point j + 1
#     R_xj(j) = X_eval - X_j(j)                                # X_eval - X_j
#     r_xj(j) = norm(R_xj(j))                                  # | X_eval - X_j |
#
#     gamma(j) = CalcVS(boundary, alpha, X_j(j)) .* [0, 0, 1]  # Vortex sheet strength
#     del_S(j) = norm(X_j(j) - X_jp1(j))                       # ΔS_j
#
#     num(j) = cross(X_eval - X_j(j), gamma(j)) * del_S(j)     # Numerator
#     den(j) = 2*pi * r_xj(j)^2                                # Denominator
#
#     # === Calculate velocity ===
#     U = sum( [((r_xj(j) < 1e-12) ? zeros(3) : num(j) ./ den(j)) for j in 1:NPTS] )
#     return U
#
# end


text:
using Test
using Agents

include("forest_fire_defs.jl")
include("schelling_defs.jl")
include("space_test.jl")
include("data_collector_test.jl")
include("schelling_test.jl")


text:
using LightsheetImageTools
using Test

@testset "LightsheetImageTools.jl" begin
    # Write your tests here.
end


text:
using CUDA

export NodeIds, ⋁NodeIds, ⋀NodeIds, BitCircuit,
    num_nodes, num_elements, num_leafs, num_features,
    nodes, elements,
    has_single_child, sibling

#####################
# Bit Circuits
#####################

"Integer identifier for a circuit node"
const NodeId = Int32 # somehow faster than UInt32?

# In a bits circuit, id
# 1 is true, 2 is false
const TRUE_BITS = NodeId(1)
const FALSE_BITS = NodeId(2)
# 3:nf+2 are nf positive literals
# nf+3:2nf+2 are nf negative literals
# 2nf+3:end are inner decision nodes

"The BitCircuit ids associated with a node"
abstract type NodeIds end
mutable struct ⋁NodeIds <: NodeIds #somehow mutable is faster
    layer_id::NodeId
    node_id::NodeId
end
mutable struct ⋀NodeIds <: NodeIds #somehow mutable is faster
    layer_id::NodeId
    prime_id::NodeId
    sub_id::NodeId
    ⋀NodeIds(p, s) = begin
        l = max(p.layer_id, s.layer_id)
        new(l, p.node_id, s.node_id)
    end
end

"""
A bit circuit is a low-level representation of a logical circuit structure.

They are a "flat" representation of a circuit, essentially a bit string,
that can be processed by lower level code (i.e., GPU kernels)

The wiring of the circuit is captured by two matrices: nodes and elements.
  * Nodes are either leafs or decision (disjunction) nodes in the circuit.
  * Elements are conjunction nodes in the circuit.
  * In addition, there is a vector of layers, where each layer is a list of node ids.
    Layer 1 is the leaf/input layer. Layer end is the circuit root.
  * And there is a vector of parents, pointing to element id parents of decision nodes.

Nodes are represented as a 4xN matrix where
  * nodes[1,:] is the first element id belonging to this decision
  * nodes[2,:] is the last element id belonging to this decision
  * nodes[3,:] is the first parent index belonging to this decision
  * nodes[4,:] is the last parent index belonging to this decision

  Elements belonging to node `i` are `elements[:, nodes[1,i]:nodes[2,i]]`
  Parents belonging to node `i` are `parents[nodes[3,i]:nodes[4,i]]`

Elements are represented by a 3xE matrix, where
  * elements[1,:] is the decision node id (parents of the element),
  * elements[2,:] is the prime node id (child of the element)
  * elements[3,:] is the sub node id (child of the element)
"""
struct BitCircuit{V,M}
    layers::Vector{V}
    nodes::M
    elements::M
    parents::V
end

function BitCircuit(circuit::LogicCircuit, data; reset=true, on_decision=noop)
    BitCircuit(circuit, num_features(data); reset, on_decision)
end

"construct a new `BitCircuit` accomodating the given number of features"
function BitCircuit(circuit::LogicCircuit, num_features::Int; reset=true, on_decision=noop)
    #TODO: consider not using foldup_aggregate and instead calling twice to ensure order but save allocations
    #TODO add inbounds annotations

    f_con(n) = ⋁NodeIds(one(NodeId), istrue(n) ? TRUE_BITS : FALSE_BITS)

    f_lit(n) = ⋁NodeIds(one(NodeId),
        ispositive(n) ? NodeId(2+variable(n)) : NodeId(2+num_features+variable(n)))

    # store data in vectors to facilitate push!
    num_leafs = 2+2*num_features
    layers::Vector{Vector{NodeId}} = Vector{NodeId}[collect(1:num_leafs)]
    nodes::Vector{NodeId} = zeros(NodeId, 4*num_leafs)
    elements::Vector{NodeId} = NodeId[]
    parents::Vector{Vector{NodeId}} = Vector{NodeId}[NodeId[] for i = 1:num_leafs]
    last_dec_id::NodeId = 2*num_features+2
    last_el_id::NodeId = zero(NodeId)

    to⋁NodeIds(c::⋁NodeIds) = c
    to⋁NodeIds(c::⋀NodeIds) = begin
            # need to add a dummy decision node in between AND nodes
            last_dec_id += one(NodeId)
            last_el_id += one(NodeId)
            push!(elements, last_dec_id, c.prime_id, c.sub_id)
            push!(parents[c.prime_id], last_el_id)
            push!(parents[c.sub_id], last_el_id)
            layer_id = c.layer_id + one(NodeId)
            push!(nodes, last_el_id, last_el_id, zero(NodeId), zero(NodeId))
            push!(parents, NodeId[])
            length(layers) < layer_id && push!(layers, NodeId[])
            push!(layers[layer_id], last_dec_id)
            on_decision(nothing, c, layer_id, last_dec_id, last_el_id, last_el_id)
            ⋁NodeIds(layer_id, last_dec_id)
    end

    f_and(n, cs) = begin
        @assert length(cs) > 1 "BitCircuits only support AND gates with at least two children"
        a12 = ⋀NodeIds(to⋁NodeIds(cs[1]), to⋁NodeIds(cs[2]))
        if length(cs) == 2
            a12
        elseif length(cs) > 1000 # Avoid using recursion here since it will cause StackOverflowError
            ⋀node_idx = cs[end]
            for idx = length(cs) - 1 : -1 : 1
                ⋀node_idx = ⋀NodeIds(to⋁NodeIds(cs[idx]), to⋁NodeIds(⋀node_idx))
            end
            ⋀node_idx
        else
            f_and(n, [a12, cs[3:end]...])
        end
    end

    f_or(n, cs) = begin
        first_el_id::NodeId = last_el_id + one(NodeId)
        layer_id::NodeId = zero(NodeId)
        last_dec_id::NodeId += one(NodeId)

        f_or_child(c::⋀NodeIds) = begin
            layer_id = max(layer_id, c.layer_id)
            last_el_id += one(NodeId)
            push!(elements, last_dec_id, c.prime_id, c.sub_id)
            @inbounds push!(parents[c.prime_id], last_el_id)
            @inbounds push!(parents[c.sub_id], last_el_id)
        end

        f_or_child(c::⋁NodeIds) = begin
            layer_id = max(layer_id, c.layer_id)
            last_el_id += one(NodeId)
            push!(elements, last_dec_id, c.node_id, TRUE_BITS)
            @inbounds push!(parents[c.node_id], last_el_id)
            @inbounds push!(parents[TRUE_BITS], last_el_id)
        end

        foreach(f_or_child, cs)

        layer_id += one(NodeId)
        length(layers) < layer_id && push!(layers, NodeId[])
        push!(nodes, first_el_id, last_el_id, zero(NodeId), zero(NodeId))
        push!(parents, NodeId[])
        push!(layers[layer_id], last_dec_id)
        on_decision(n, cs, layer_id, last_dec_id, first_el_id, last_el_id)
        ⋁NodeIds(layer_id, last_dec_id)
    end

    r = foldup_aggregate(circuit, f_con, f_lit, f_and, f_or, NodeIds; reset)
    to⋁NodeIds(r)

    nodes_m = reshape(nodes, 4, :)
    elements_m = reshape(elements, 3, :)
    parents_m = Vector{NodeId}(undef, size(elements_m,2)*2)
    last_parent = zero(NodeId)
    @assert last_dec_id == size(nodes_m,2) == size(parents,1)
    @assert sum(length, parents) == length(parents_m)
    for i in 1:last_dec_id-1
        if !isempty(parents[i])
            nodes_m[3,i] = last_parent + one(NodeId)
            parents_m[last_parent + one(NodeId):last_parent + length(parents[i])] .= parents[i]
            last_parent += length(parents[i])
            nodes_m[4,i] = last_parent
        else
            @assert i <= num_leafs "Only root and leaf nodes can have no parents: $i"
        end
    end

    return BitCircuit(layers, nodes_m, elements_m, parents_m)
end

import .Utils: num_nodes # extend

nodes(c::BitCircuit) = c.nodes
elements(c::BitCircuit) = c.elements

"How many nodes are indexed by a given bit circuit?"
num_nodes(c::BitCircuit) = size(c.nodes, 2)

"Number of elements (conjunctions) in layer or bit circuit"
num_elements(c::BitCircuit) = size(c.elements, 2)
num_elements(nodes, id) =
    @inbounds (nodes[2,id] - nodes[1,id] + 1)

import .Utils: num_features #extend

num_leafs(c::BitCircuit) =
    num_leafs(length(c.layers[1]))
num_leafs(nl::Int) = nl

num_features(c::BitCircuit) = num_features(num_leafs(c))
num_features(nl::Int) = (nl-2) ÷ 2

isliteralgate(c::BitCircuit, nodeid) =
    isliteralgate(num_leafs(c), nodeid)
isliteralgate(nl::Int, nodeid) =
    2 < nodeid <= nl

isleafgate(c::BitCircuit, nodeid) =
    isleafgate(num_leafs(c), nodeid)
isleafgate(nl::Int, nodeid) =
    nodeid <= nl

function literal(c, nodeid)
    @assert 2 < nodeid <= num_leafs(c)
    if nodeid <= 2 + num_features(c)
        Lit(nodeid-2)
    else
        -Lit(nodeid-2-num_features(c))
    end
end

"Does the bitcircuit node have a single child?"
@inline has_single_child(nodes, id) =
    @inbounds (nodes[1,id] == nodes[2,id])

"Get the other child of a parent element"
@inline sibling(els, par, me) =
    @inbounds ((els[2,par] == me) ? els[3,par] : els[2,par])

import .Utils: to_gpu, to_cpu, isgpu #extend

to_gpu(c::BitCircuit) =
    BitCircuit(map(to_gpu, c.layers), to_gpu(c.nodes), to_gpu(c.elements), to_gpu(c.parents))

to_cpu(c::BitCircuit) =
    BitCircuit(map(to_cpu, c.layers), to_cpu(c.nodes), to_cpu(c.elements), to_cpu(c.parents))

isgpu(c::BitCircuit{<:CuArray,<:CuArray}) = true
isgpu(c::BitCircuit{<:Array,<:Array}) = false


text:
import FFTW
import Luna: RK45
import Test: @test

function testinit()
    trange = 32
    ωrange = 160
    samples = 2^(ceil(Int, log2(trange*ωrange/2π)))
    δt = trange/samples
    δω = 2π/trange
    n = collect(range(0, length=samples))
    t = @. (n-samples/2)*δt
    ω = @. (n-samples/2)*δω
    N = 5
    z0 = π/2
    zmax = z0*2

    At = @. N*sech(t)
    Aω = FFTW.fftshift(FFTW.fft(At))

    At = zero(Aω)
    nl = zero(Aω)
    nlω = zero(Aω)

    FT = FFTW.plan_fft!(Aω)
    IFT = FFTW.plan_ifft!(Aω)
    f! = let FT=FT, IFT=IFT, ω=ω
        function f!(out, Aω, z)
                out .= Aω
                IFT*out
                out .= abs2.(out).*out
                FT*out
                @. out = im*out - im/2*ω^2*Aω
                return out
        end
    end

    Lin = @. -im/2*ω^2

    linout = similar(Aω)
    Linfunc = let ω=ω
        function Linfunc(out, z)
            @. out = -im/2*ω^2
        end
    end

    fnl! = let FT=FT, IFT=IFT
        function fnl!(out, Aω, z)
                out .= Aω
                IFT*out
                out .= abs2.(out).*out
                FT*out
                @. out = im*out
                return out
        end
    end

    return t, ω, zmax, Aω, f!, Lin, fnl!, Linfunc
end

function test_precon(plot=false)
    t, ω, zmax, Aω, f!, Lin, fnl!, Linfunc = testinit()
    z = 0
    dz = 1e-3

    saveN = 501

    zarr, Aarr = RK45.solve_precon(fnl!, Lin, Aω, z, dz, zmax, saveN, rtol=1e-6, status_period=2)

    if plot
        Atarr = FFTW.ifft(FFTW.ifftshift(Aarr, 1), 1)
        energy = dropdims(sum(abs2.(Atarr), dims=1), dims=1)
        pygui(true)
        plt.figure()
        plt.pcolormesh(t, zarr, abs2.(transpose(Atarr)))
        plt.colorbar()
        plt.figure()
        plt.pcolormesh(ω, zarr, abs2.(transpose(Aarr)))
        plt.figure()
        plt.plot(zarr, 1 .- energy/energy[1])
        plt.figure()
        plt.plot(t, abs2.(Atarr[:, 1]))
        plt.plot(t, abs2.(Atarr[:, end]))
    end

    return zarr, Aarr
end

function test_noprecon(plot=false)
    t, ω, zmax, Aω, f!, Lin, fnl!, Linfunc = testinit()
    z = 0
    dz = 1e-3

    saveN = 201

    zarr, Aarr = RK45.solve(f!, Aω, z, dz, zmax, saveN)

    if plot
        Atarr = FFTW.ifft(FFTW.ifftshift(Aarr, 1), 1)
        energy = dropdims(sum(abs2.(Atarr), dims=1), dims=1)
        pygui(true)
        plt.figure()
        plt.pcolormesh(t, zarr, abs2.(transpose(Atarr)))
        plt.colorbar()
        plt.figure()
        plt.pcolormesh(ω, zarr, abs2.(transpose(Aarr)))
        plt.figure()
        plt.plot(zarr, 1 .- energy/energy[1])
        plt.figure()
        plt.plot(t, abs2.(Atarr[:, 1]))
        plt.plot(t, abs2.(Atarr[:, end]))
    end

    return zarr, Aarr
end

t, ω, zmax, Aω, f!, Lin, fnl!, Linfunc = testinit()
z = 0
dz = 1e-3

zarr, Aarr = RK45.solve(f!, copy(Aω), z, dz, zmax, rtol=1e-8, output=true, outputN=501)
zarrp, Aarrp = RK45.solve_precon(fnl!, Lin, copy(Aω), z, dz, zmax,
                                 rtol=1e-8, output=true, outputN=501)
zarrpf, Aarrpf = RK45.solve_precon(fnl!, Linfunc, copy(Aω), z, dz, zmax,
                                   rtol=1e-8, output=true, outputN=501)
# Is the initial spectrum restored after 2 soliton periods? (without preconditioner)
@test isapprox(abs2.(Aarr[:, 1]), abs2.(Aarr[:, end]), rtol=1e-4)
# Is the initial spectrum restored after 2 soliton periods? (with preconditioner)
@test isapprox(abs2.(Aarrp[:, 1]), abs2.(Aarrp[:, end]), rtol=1e-3)
# Is the initial spectrum restored after 2 soliton periods?
# (with preconditioner and z-dependent linear part)
@test isapprox(abs2.(Aarrpf[:, 1]), abs2.(Aarrpf[:, end]), rtol=1e-3)
# Is there a difference if the linear part is a function (but constant)?
@test all(abs2.(Aarrp) .== abs2.(Aarrpf))


text:
# Cross Entropy loss function

export CrossEntropy

"""
**Cross Entropy loss**

    CrossEntropy(b)

Returns the function
```math
f(x) = -\\frac{1}{N} \\sum_{i = 1}^{N} b_i \\log (x_i)+(1-b_i) \\log (1-x_i),
```
where `b` is an array such that `0 ≤ b ≤ 1` component-wise.
"""
struct CrossEntropy{R <: Real, T <: AbstractArray{R}} <: ProximableFunction
    b::T
    function CrossEntropy{R, T}(b::T) where {R <: Real, T <: AbstractArray{R}}
        if !(all(0 .<= b .<= 1. ))
            error("b must be 0 ≤ b ≤ 1 ")
        else
            new(b)
        end
    end
end

is_convex(f::CrossEntropy) = true
is_smooth(f::CrossEntropy) = true

CrossEntropy(b::T) where {R <: Real, T <: AbstractArray{R}} = CrossEntropy{R, T}(b)

function (f::CrossEntropy{R})(x::AbstractArray{R}) where {R <: Real}
    sum = R(0)
    for i in eachindex(f.b)
        sum += f.b[i]*log(x[i])+(1-f.b[i])*log(1-x[i])
    end
    return -sum/length(f.b)
end

function (f::CrossEntropy{B})(x::AbstractArray{R}) where {B <: Bool, R <: Real}
    sum = R(0)
    for i in eachindex(f.b)
        sum += f.b[i] ? log(x[i]) : (1-f.b[i])*log(1-x[i])
    end
    return -sum/length(f.b)
end

function gradient!(y::AbstractArray{R}, f::CrossEntropy{R}, x::AbstractArray{R}) where {R <: Real}
    sum = R(0)
    for i in eachindex(x)
        y[i] = 1/length(f.b)*( - f.b[i]/x[i] + (1-f.b[i])/(1-x[i]) )
        sum += f.b[i]*log(x[i])+(1-f.b[i])*log(1-x[i])
    end
    return -1/length(f.b)*sum
end

function gradient!(y::AbstractArray{R}, f::CrossEntropy{B}, x::AbstractArray{R}) where {R <: Real, B <: Bool}
    sum = R(0)
    for i in eachindex(x)
        y[i] = f.b[i] ? - 1/x[i] : 1/(1-x[i])
        y[i] *= 1/length(f.b)
        sum += f.b[i] ? log(x[i]) : log(1-x[i])
    end
    return -1/length(f.b)*sum
end

function prox!(y::AbstractArray{R}, f::CrossEntropy, x::AbstractArray{R}, gamma::R=R(1)) where {R}
    # TODO: fill-in here
    error("not implemented")
end

fun_name(f::CrossEntropy) = "cross entropy loss"
fun_dom(f::CrossEntropy) = "AbstractArray{Real}"
fun_expr(f::CrossEntropy) = "x ↦ 1/N sum( b*log(x)+(1-b)*log(1-x), i=1,...,N )"
fun_params(f::CrossEntropy) = "b = $(typeof(f.b))"


text:
export VCAT

"""
`VCAT(A::AbstractOperator...)`

Shorthand constructors:

`[A1; A2 ...]`

`vcat(A...)`

Vertically concatenate `AbstractOperator`s. Notice that all the operators must share the same domain dimensions and type, e.g. `size(A1,2) == size(A2,2)` and `domainType(A1) == domainType(A2)`.

```julia
julia> VCAT(DFT(4,4),Variation((4,4)))
[ℱ;Ʋ]  ℝ^(4, 4) -> ℂ^(4, 4)  ℝ^(16, 2)

julia> V = [Eye(3); DiagOp(2*ones(3))]
[I;╲]  ℝ^3 -> ℝ^3  ℝ^3

julia> vcat(V,FiniteDiff((3,)))
VCAT  ℝ^3 -> ℝ^3  ℝ^3  ℝ^2
```

When multiplying a `VCAT` with an array of the proper size, the result will be a `Tuple` containing arrays with the `VCAT`'s codomain type and size.

```julia
julia> V*ones(3)
([1.0, 1.0, 1.0], [2.0, 2.0, 2.0])

```

"""
struct VCAT{N, # number of AbstractOperator
            L <: NTuple{N,AbstractOperator},
            P <: Tuple,
            C <: AbstractArray,
           } <: AbstractOperator
	A::L     # tuple of AbstractOperators
	idxs::P  # indices
	         # H = VCAT(Eye(n),VCAT(Eye(n),Eye(n))) has H.idxs = (1,2,3)
           # `AbstractOperators` are flatten
           # H = VCAT(Eye(n),Compose(MatrixOp(randn(n,n)),VCAT(Eye(n),Eye(n))))
           # has H.idxs = (1,(2,3))
           # `AbstractOperators` are stack
	buf::C   # buffer memory
  function VCAT(A::L, idxs::P, buf::C) where {N,
                  L <: NTuple{N,AbstractOperator},
                  P <: Tuple,
                  C <: AbstractArray,
                  }
    if any([size(A[1],2) != size(a,2) for a in A])
      throw(DimensionMismatch("operators must have the same domain dimension!"))
    end
    if any([domainType(A[1]) != domainType(a) for a in A])
      throw(error("operators must all share the same domainType!"))
    end
    new{N,L,P,C}(A, idxs, buf)
  end
end

function VCAT(A::Vararg{AbstractOperator})
	if any((<:).(typeof.(A),VCAT)) #there are VCATs in A
		AA = ()
		for a in A
			if typeof(a) <: VCAT # flatten
				AA = (AA...,a.A...)
			else                 # stack
				AA = (AA...,a)
			end
		end
		# use buffer from VCAT in A
		buf = A[findfirst( (<:).(typeof.(A),VCAT) ) ].buf
	else
		AA = A
		s = size(AA[1],2)
		t = domainType(AA[1])
		# generate buffer
    buf = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s))
	end

	return VCAT(AA, buf)
end

function VCAT(AA::NTuple{N,AbstractOperator}, buf::C) where {N,C}
	if N == 1
		return AA[1]
	else
		# build H.idxs
		K = 0
		idxs = []
		for i in eachindex(ndoms.(AA,1))
			if ndoms(AA[i],1) == 1 # flatten operator
				K += 1
				push!(idxs,K)
			else                   # stacked operator
				idxs = push!(idxs,(collect(K+1:K+ndoms(AA[i],2))...,) )
				for ii = 1:ndoms(AA[i],1)
					K += 1
				end
			end
		end
		return VCAT(AA, (idxs...,), buf)
	end
end

VCAT(A::AbstractOperator) = A

# Mappings
@generated function mul!(y::C, A::AdjointOperator{VCAT{N,L,P,C}}, b::DD) where {N,L,P,C,DD <: ArrayPartition}
  ex = :(H = A.A)

  if fieldtype(P,1) <: Int
    # flatten operator
    # build mul!(y, H.A[1]', b.x[H.idxs[1]])
    bb = :(b.x[H.idxs[1]])
	else
    # stacked operator
    # build mul!(y, H.A[1]',ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ...  ))
    bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))]
    bb = :( ArrayPartition($(bb...)) )
	end
	ex = :($ex; mul!(y,H.A[1]',$bb)) # write on y

  for i = 2:N
    if fieldtype(P,i) <: Int
      # flatten operator
      # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]])
      bb = :(b.x[H.idxs[$i]])
		else
      # stacked operator
      # build mul!(H.buf, H.A[i]',( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ...  ))
      bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))]
      bb = :( ArrayPartition( $(bb...) ) )
		end
    ex = :($ex; mul!(H.buf,H.A[$i]',$bb)) # write on H.buf
    # sum H.buf with y
    ex = :($ex; y .+= H.buf)
	end
	ex = :($ex; return y)
	return ex
end

@generated function mul!(y::DD, H::VCAT{N,L,P,C}, b::C) where {N,L,P,C,DD <: ArrayPartition}
  ex = :()
  for i = 1:N
    if fieldtype(P,i) <: Int
      # flatten operator
      # build mul!(y.x[H.idxs[i]], H.A[i], b)
      yy = :(y.x[H.idxs[$i]])
    else
      # stacked operator
      # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ...  ), H.A[i], b)
      yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))]
      yy = :(ArrayPartition( $(yy...) ) )
    end
    ex = :($ex; mul!($yy,H.A[$i],b))
	end
	ex = :($ex; return y)
	return ex
end

## same as mul! but skips `Zeros`
@generated function mul_skipZeros!(y::C, A::AdjointOperator{VCAT{N,L,P,C}}, b::DD) where {N,L,P,C,DD <: ArrayPartition}
  ex = :(H = A.A)

  if fieldtype(P,1) <: Int
    # flatten operator
    # build mul!(y, H.A[1]', b.x[H.idxs[1]])
    bb = :(b.x[H.idxs[1]])
	else
    # stacked operator
    # build mul!(y, H.A[1]',ArrayPartition( b.x[H.idxs[1][1]], b.x[H.idxs[1][2]] ...  ))
    bb = [ :(b.x[H.idxs[1][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,1)))]
    bb = :( ArrayPartition($(bb...)) )
	end
  ex = :($ex; mul!(y,H.A[1]',$bb)) # write on y

  for i = 2:N
    if !(fieldtype(L,i) <: Zeros)
      if fieldtype(P,i) <: Int
        # flatten operator
        # build mul!(H.buf, H.A[i]', b.x[H.idxs[i]])
        bb = :(b.x[H.idxs[$i]])
      else
        # stacked operator
        # build mul!(H.buf, H.A[i]',( b.x[H.idxs[i][1]], b.x[H.idxs[i][2]] ...  ))
        bb = [ :( b.x[H.idxs[$i][$ii]] ) for ii in eachindex(fieldnames(fieldtype(P,i)))]
        bb = :( ArrayPartition( $(bb...) ) )
      end
      ex = :($ex; mul!(H.buf,H.A[$i]',$bb)) # write on H.buf
      # sum H.buf with y
      ex = :($ex; y .+= H.buf)
    end
	end
	ex = :($ex; return y)
	return ex
end

@generated function mul_skipZeros!(y::DD, H::VCAT{N,L,P,C}, b::C) where {N,L,P,C,DD <: ArrayPartition}
  ex = :()
  for i = 1:N
    if !(fieldtype(L,i) <: Zeros)
      if fieldtype(P,i) <: Int
        # flatten operator
        # build mul!(y.x[H.idxs[i]], H.A[i], b)
        yy = :(y.x[H.idxs[$i]])
      else
        # stacked operator
        # build mul!(ArrayPartition( y[.xH.idxs[i][1]], y.x[H.idxs[i][2]] ...  ), H.A[i], b)
        yy = [ :(y.x[H.idxs[$i][$ii]]) for ii in eachindex(fieldnames(fieldtype(P,i)))]
        yy = :(ArrayPartition( $(yy...) ) )
      end
      ex = :($ex; mul!($yy,H.A[$i],b))
    end
	end
	ex = :($ex; return y)
	return ex
end

# Properties

function size(H::VCAT)
	size_out = []
	for s in size.(H.A,1)
		eltype(s) <: Int ? push!(size_out,s) : push!(size_out,s...)
	end
	p = vcat([[idx... ] for idx in H.idxs]...)
	invpermute!(size_out,p)

	(size_out...,), size(H.A[1],2)
end

fun_name(L::VCAT) = length(L.A) == 2 ? "["*fun_name(L.A[1])*";"*fun_name(L.A[2])*"]" : "VCAT"

function codomainType(H::VCAT)
	codomain = vcat([typeof(d)<:Tuple ? [d...] : d  for d in codomainType.(H.A)]...)
	p = vcat([[idx... ] for idx in H.idxs]...)
	invpermute!(codomain,p)
	return (codomain...,)
end
domainType(L::VCAT) = domainType.(Ref(L.A[1]))

is_linear(L::VCAT) = all(is_linear.(L.A))
is_AcA_diagonal(L::VCAT) = all(is_AcA_diagonal.(L.A))
is_full_column_rank(L::VCAT) = any(is_full_column_rank.(L.A))

diag_AcA(L::VCAT) = (+).(diag_AcA.(L.A)...,)

# utils
function permute(H::VCAT{N,L,P,C}, p::AbstractVector{Int}) where {N,L,P,C}

	unfolded = vcat([[idx... ] for idx in H.idxs]...)
	invpermute!(unfolded,p)

	new_part = ()
	cnt = 0
	for z in length.(H.idxs)
		new_part = (new_part..., z == 1 ? unfolded[cnt+1] : (unfolded[cnt+1:z+cnt]...,))
		cnt += z
	end

	VCAT{M,N,L,P,C}(H.A,new_part,H.buf)
end

remove_displacement(V::VCAT) = VCAT(remove_displacement.(V.A), V.idxs, V.buf)


text:
@testset "StrainTrait" begin

    mat = NeoHook(λ = 1.0, μ = 1.0)
    state = initial_material_state(mat)
    F = Tensor{2,3}((2.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0))
    Δt = 0.0
    I = one(SymmetricTensor{2,3,Float64})

    S, dSdC, state = material_response(∂S∂C(), mat, F, state, Δt)

    Pᵀ, dPᵀdF, state = material_response(∂Pᵀ∂F(), mat, F, state, Δt)
    @test Pᵀ == S⋅F'
    @test dPᵀdF == otimesu(F,I) ⊡ dSdC ⊡ otimesu(F',I) + otimesu(S,I)

    S_E, dSdE, state = material_response(∂S∂E(), mat, F, state, Δt)
    @test S_E == S
    @test 2dSdC == dSdE

end