text:
using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
      Test, BandedMatrices, FillArrays

Random.seed!(1234)

# These functions create analytic solutions to upwind operators. In the function name,
# the first number indicates the derivative order, the second number indicates
# the approximation order,and the last three characters "Pos" or "Neg" indicating
# the winding direction. For instance analyticOneTwoPos corresponds to the positive
# upwind second order approximation of the first derivative.

function analyticOneOnePos()
      A = zeros(5,7)
      for i in 1:5
            A[i,i+1:i+2] = [-1 1]
      end
      return A
end

function analyticOneOneNeg()
      A = zeros(5,7)
      for i in 1:5
            A[i,i:i+1] = [-1 1]
      end
      return A
end

function analyticOneTwoPos()
      A = zeros(5,7)
      for i in 1:4
            A[i,i+1:i+3] = [-3/2 2 -1/2]
      end
      A[5,5:7] = [-1/2 0 1/2]
      return A
end

function analyticOneTwoNeg()
      A = zeros(5,7)
      A[1,1:3] = [-1/2 0 1/2]
      for i in 2:5
            A[i,i-1:i+1] = [1/2 -2 3/2]
      end
      return A
end

function analyticTwoTwoPos()
      A = zeros(5,7)
      for i in 1:3
            A[i,i+1:i+4] = [2 -5 4 -1]
      end
      A[4,4:7] = [1 -2 1 0]
      A[5, 4:7] = [0 1 -2 1]
      return A
end

function analyticTwoTwoNeg()
      A = zeros(5,7)
      A[1,1:4] = [1 -2 1 0]
      A[2,1:4] = [0 1 -2 1]
      for i in 3:5
            A[i,i-2:i+1] = [-1 4 -5 2]
      end
      return A
end

function analyticTwoThreePos()
      A = zeros(7,9)
      for i in 1:4
            A[i,i+1:i+5] = [35/12 -104/12 114/12 -56/12 11/12]
      end
      A[5,5:9] = [11/12 -20/12 6/12 4/12 -1/12]
      A[6,5:9] = [-1/12 16/12 -30/12 16/12 -1/12]
      A[7,5:9] = [-1/12 4/12 6/12 -20/12 11/12]
      return A
end

function analyticTwoThreeNeg()
      A = zeros(7,9)
      A[1,1:5] = [11/12 -20/12 6/12 4/12 -1/12]
      A[2,1:5] = [-1/12 16/12 -30/12 16/12 -1/12]
      A[3,1:5] = [-1/12 4/12 6/12 -20/12 11/12]
      for i in 4:7
            A[i,i-3:i+1] = [11/12 -56/12 114/12 -104/12 35/12]
      end
      return A
end

# analytically derived upwind operators for the irregular grid:
# [0.0, 0.08, 0.1, 0.15, 0.19, 0.26, 0.29], and spacing:
# dx = [0.08, 0.02, 0.05, 0.04, 0.07, 0.03]

function analyticOneOnePosIrr()
      A = zeros(5,7)
      A[1,2:3] = [-50, 50]
      A[2,3:4] = [-20, 20]
      A[3,4:5] = [-25, 25]
      A[4,5:6] = [-100/7, 100/7]
      A[5,6:7] = [-100/3, 100/3]
      return A
end

function analyticOneOneNegIrr()
      A = zeros(5,7)
      A[1,1:2] = [-25/2, 25/2]
      A[2,2:3] = [-50, 50]
      A[3,3:4] = [-20, 20]
      A[4,4:5] = [-25, 25]
      A[5,5:6] = [-100/7, 100/7]
      return A
end

function analyticOneTwoPosIrr()
      A = zeros(5,7)
      A[1,2:4] = [-450/7, 490/7, -40/7]
      A[2,3:5] = [-280/9, 405/9, -125/9]
      A[3,4:6] = [-2625/77, 3025/77, -400/77]
      A[4,5:7] = [-510/21, 1000/21, -490/21]
      A[5,5:7] = [-90/21, -400/21, 490/21]
      return A
end

function analyticOneTwoNegIrr()
      A = zeros(5,7)
      A[1,1:3] = [-5/2, -75/2, 80/2]
      A[2,1:3] = [5/2, -125/2, 120/2]
      A[3, 2:4] = [250/7, -490/7, 240/7]
      A[4, 3:5] = [80/9, -405/9, 325/9]
      A[5, 4:6] = [1225/77, -3025/77, 1800/77]
      return A
end

function analyticTwoTwoPosIrr()
      A = zeros(5,7)
      A[1,2:5] = [200000/77, -308000/77, 143000/77, -35000/77]
      A[2,3:6] = [27500/33, -75000/33, 55000/33, -7500/33]
      A[3,4:7] = [72500/77, -137500/77, 120000/77, -55000/77]
      A[4,4:7] = [42500/77, -71500/77, 40000/77, -11000/77]
      A[5,4:7] = [-10000/77, 44000/77, -100000/77, 66000/77]
      return A
end

function analyticTwoTwoNegIrr()
      A = zeros(5,7)
      A[1,1:4] = [1050/7, -1250/7, -1400/7, 1600/7]
      A[2,1:4] = [350/7, 6250/7, -9800/7, 3200/7]
      A[3,1:4] = [-1400/7, 25000/7, -30800/7, 7200/7]
      A[4,2:5] = [-390000/231, 770000/231, -660000/231, 280000/231]
      A[5,3:6] = [-38500/77, 161000/77, -165000/77, 42500/77]
      return A
end

@testset "Test: Derivative Order = 1, Approx Order = 1, Winding = Positive" begin
      N = 5
      L = UpwindDifference(1,1, 1.0, N, 1.0)
      analyticL = analyticOneOnePos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 1, Approx Order = 1, Winding = Negative" begin
      N = 5
      L = UpwindDifference(1,1, 1.0, N, -1.0)
      analyticL = -1*analyticOneOneNeg()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 1, Approx Order = 2, Winding = Positive" begin
      N = 5
      L = UpwindDifference(1,2, 1.0, N, 1.0)
      analyticL = analyticOneTwoPos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 1, Approx Order = 2, Winding = Negative" begin
      N = 5
      L = UpwindDifference(1,2, 1.0, N, -1.0)
      analyticL = -1*analyticOneTwoNeg()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 2, Approx Order = 2, Winding = Positive" begin
      N = 5
      L = UpwindDifference(2,2, 1.0, N, 1.0)
      analyticL = analyticTwoTwoPos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 2, Approx Order = 2, Winding = Negative" begin
      N = 5
      L = UpwindDifference(2,2, 1.0, N, -1.0)
      analyticL = -1*analyticTwoTwoNeg()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) == analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

# Here the operators are too big for five grid points, so weird corner cases must be accounted for
# We should be able to assume that users will not have cases like this.
@testset "Test: Derivative Order = 2, Approx Order = 3, Winding = Positive" begin
      N = 7
      L = UpwindDifference(2,3, 1.0, N, 1.0)
      analyticL = analyticTwoThreePos()
      x = rand(9)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 2, Approx Order = 3, Winding = Negative" begin
      N = 7
      L = UpwindDifference(2,3, 1.0, N, -1.0)
      analyticL = -1*analyticTwoThreeNeg()
      x = rand(9)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 1, Approx Order = 1, Winding = Positive, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(1,1, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, 1.0)
      analyticL = analyticOneOnePosIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv
end

@testset "Test: Derivative Order = 1, Approx Order = 1, Winding = Negative, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(1,1, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, -1.0)
      analyticL = -1*analyticOneOneNegIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv
end

@testset "Test: Derivative Order = 1, Approx Order = 2, Winding = Positive, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(1,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, 1.0)
      analyticL = analyticOneTwoPosIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 1, Approx Order = 2, Winding = Negative, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(1,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, -1.0)
      analyticL = -1*analyticOneTwoNegIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 2, Approx Order = 2, Winding = Positive, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(2,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, 1.0)
      analyticL = analyticTwoTwoPosIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Derivative Order = 2, Approx Order = 2, Winding = Negative, Grid = Irregular" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(2,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, -1.0)
      analyticL = -1*analyticTwoTwoNegIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Scaling by dx and Derivative Order in Uniform Case" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(1,2, 0.1, N, 1.0)
      analyticL = 10.0*analyticOneTwoPos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(1,2, 0.1, N, -1.0)
      analyticL = -10.0*analyticOneTwoNeg()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(2,2, 0.1, N, 1.0)
      analyticL = 100.0*analyticTwoTwoPos()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(2,2, 0.1, N, -1.0)
      analyticL = -100.0*analyticTwoTwoNeg()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Non-Trivial Coefficient Handling in Uniform Grid Case" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(2,2, 1.0, N, 4.56)
      analyticL = 4.56*analyticTwoTwoPos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(2,2, 1.0, N, -4.56)
      analyticL = -4.56*analyticTwoTwoNeg()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: dx and Derivative Order Scaling and Non-Trivial Coefficient Handling in Uniform Grid Case" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(2,2, 0.1, N, 4.56)
      analyticL = 4.56*100.0*analyticTwoTwoPos()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(2,2, 0.1, N, -4.56)
      analyticL = -4.56*100.0*analyticTwoTwoNeg()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end

@testset "Test: Coefficient Handling in Non-Uniform Grid Case" begin
      N = 5
      # constructor throws an error at the moment
      L = UpwindDifference(2,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, 4.56)
      analyticL = 4.56*analyticTwoTwoPosIrr()
      x = rand(7)

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

      L = UpwindDifference(2,2, [0.08, 0.02, 0.05, 0.04, 0.07, 0.03], N, -4.56)
      analyticL = -4.56*analyticTwoTwoNegIrr()

      # Test that multiplication agrees with analytic multiplication
      @test L*x ≈ analyticL*x

      # Test that concretized multiplication agrees with analytic multiplication
      @test Array(L)*x ≈ analyticL*x

      # Test that matrix-free multiplication agrees with concretized multiplication
      @test L*x ≈ Array(L)*x

      # Test that concretized matrix agrees with analytic matrix
      @test Array(L) ≈ analyticL

      # Test Banded and Sparse concretizations
      @test Array(L) == BandedMatrix(L)
      @test Array(L) == SparseMatrixCSC(L)
      @test typeof(BandedMatrix(L)) <: BandedMatrix
      @test typeof(SparseMatrixCSC(L)) <: SparseMatrixCSC

      # Test convolutions against bpv
      bpv = PeriodicBC(Float64)*x[1:N]
      @test L*bpv ≈ analyticL*bpv

end


text:
type DEOptions{uEltype,uEltypeNoUnits,tTypeNoUnits,tType,F2,F3,F4,F5,F6,tstopsType,ECType}
  maxiters::Int
  timeseries_steps::Int
  save_timeseries::Bool
  adaptive::Bool
  abstol::uEltype
  reltol::uEltypeNoUnits
  gamma::tTypeNoUnits
  qmax::tTypeNoUnits
  qmin::tTypeNoUnits
  dtmax::tType
  dtmin::tType
  internalnorm::F2
  tstops::tstopsType
  saveat::tstopsType
  d_discontinuities::tstopsType
  userdata::ECType
  progress::Bool
  progress_steps::Int
  progress_name::String
  progress_message::F5
  timeseries_errors::Bool
  dense_errors::Bool
  beta1::tTypeNoUnits
  beta2::tTypeNoUnits
  qoldinit::tTypeNoUnits
  dense::Bool
  callback::F3
  isoutofdomain::F4
  unstable_check::F6
  verbose::Bool
  calck::Bool
  force_dtmin::Bool
  advance_to_tstop::Bool
  stop_at_next_tstop::Bool
end

type ODEIntegrator{algType<:OrdinaryDiffEqAlgorithm,uType<:Union{AbstractArray,Number},tType,tTypeNoUnits,tdirType,ksEltype,SolType,rateType,F,ProgressType,CacheType,O} <: AbstractODEIntegrator
  sol::SolType
  u::uType
  k::ksEltype
  t::tType
  dt::tType
  f::F
  uprev::uType
  uprev2::uType
  tprev::tType
  alg::algType
  rate_prototype::rateType
  notsaveat_idxs::Vector{Int}
  dtcache::tType
  dtchangeable::Bool
  dtpropose::tType
  tdir::tdirType
  EEst::tTypeNoUnits
  qold::tTypeNoUnits
  q11::tTypeNoUnits
  iter::Int
  saveiter::Int
  saveiter_dense::Int
  prog::ProgressType
  cache::CacheType
  kshortsize::Int
  just_hit_tstop::Bool
  accept_step::Bool
  isout::Bool
  reeval_fsal::Bool
  u_modified::Bool
  opts::O
  fsalfirst::rateType
  fsallast::rateType

  ODEIntegrator(sol,u,k,t,dt,f,uprev,uprev2,tprev,
      alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,
      EEst,qold,q11,
      iter,saveiter,saveiter_dense,prog,cache,
      kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts) = new(
      sol,u,k,t,dt,f,uprev,uprev2,tprev,
      alg,rate_prototype,notsaveat_idxs,dtcache,dtchangeable,dtpropose,tdir,
      EEst,qold,q11,
      iter,saveiter,saveiter_dense,prog,cache,
      kshortsize,just_hit_tstop,accept_step,isout,reeval_fsal,u_modified,opts) # Leave off fsalfirst and last
end

# When this is changed, DelayDiffEq.jl must be changed as well!

ext: jl


text:
"""
Particle belief over subspaces
"""
mutable struct MCVIBelief{S,A}
    space::MCVISubspace{S,A}
    weights::Vector{Float64}
    id::UInt64
end

MCVIBelief(space::MCVISubspace{S,A}...) where {S,A} = MCVIBelief(space...)

Base.length(b::MCVIBelief) = length(b.weights)

"""
Create initial belief particles
Is supposed to represent the initial state distribution.
Not sure how to propagate num_particles without storing it in pomdp or calling it a different function?
"""
function initial_belief(pomdp::POMDPs.POMDP{S,A,O}, num_particles::Int64, rng::AbstractRNG) where {S,A,O}
    particles = Vector{S}(undef, num_particles)
    wts = zeros(num_particles)
    for i in 1:num_particles
        particles[i] = rand(rng, initialstate(pomdp))
        wts[i]= 1.0/float(num_particles)
    end
    ss = MCVISubspace{S,A}(particles, Vector{Reward}(), Dict{A, MCVISubspace{S,A}}())
    return MCVIBelief(ss, wts, UInt64(0))
end

"""
    rand(rng::AbstractRNG, b::MCVIBelief{S,A})

Samples states from belief
"""
function rand(rng::AbstractRNG, b::MCVIBelief)
    r = Random.rand(rng)
    for (i, p) in enumerate(b.weights)
        r -= p
        if r < 0
            return particle(b.space, i)
        end
    end
    @assert false "Shouldn't reach here"
end

"""
    next(bb::MCVIBelief{S,A}, act::POMDPs.Action{A}, pomdp::POMDPs.POMDP{S,A,O})

Get the next belief according to action.
Actions create a new subspace
"""
function next(bb::MCVIBelief{S,A}, act::A, pomdp::POMDPs.POMDP, rng::AbstractRNG) where {S,A}
    next_ss = next(bb.space, act, pomdp, rng) # Get the next subspace
    return MCVIBelief(next_ss, bb.weights, UInt64(0)) # FIXME update id?
end

"""
    next(bb::MCVIBelief{S,A}, obs::POMDPs.Observation{O}, pomdp::POMDPs.POMDP{S,A,O})

Get the next belief according to observation.
Observations update the weights of the space created by actions.
"""
function next(bb::MCVIBelief{S}, obs::O, pomdp::POMDPs.POMDP) where {S,O}
    wts = weights(bb.space, obs, pomdp) # Get weights for the observation
    belief_after_obs = bb.weights .* wts
    obs_norm = sum(belief_after_obs)
    @assert obs_norm != 0 "Normalizing constant should not be zero"
    belief_after_obs /= obs_norm
    return MCVIBelief(bb.space, belief_after_obs, UInt64(0))
end

# """
#     lower_bound(b::MCVIBelief{S}, pomdp::POMDPs.POMDP, rng::AbstractRNG)
#
# Provides the minimum sum of particle values in a belief
# """
function lower_bound(lb::Any, pomdp::POMDPs.POMDP, b::MCVIBelief)
    sum = 0.0
    for (i, w) in enumerate(b.weights)
        s = particle(b.space, i)
        w = b.weights[i]
        sum += lower_bound(lb, pomdp, s) * w # Weighted lower bound for the state
    end
    return sum
end

# """
#     upper_bound(b::MCVIBelief{S}, pomdp::POMDPs.POMDP, rng::AbstractRNG)
#
# Provides the maximum sum of particle values in a belief
# """
function upper_bound(ub::Any, pomdp::POMDPs.POMDP, b::MCVIBelief)
    sum = 0.0
    for (i, w) in enumerate(b.weights)
        s = particle(b.space, i)
        w = b.weights[i]
        sum += upper_bound(ub, pomdp, s) * w # Weighted upper bound for the state
    end
    return sum
end

"""
    reward(b::MCVIBelief, pomdp::POMDPs.POMDP)

Returns the immediate reward for a given belief
"""
function reward(b::MCVIBelief, pomdp::POMDPs.POMDP)
    return reward(b.space, b.weights)
end

"""
Returns the mean value of the particle belief
"""
function Statistics.mean(b::MCVIBelief)
    μ = 0.0
    for i in 1:length(b.weights)
        μ += b.weights[i] * particle(b.space, i)
    end
    return μ
end

"""
Returns the standard deviation of the particle belief
"""
function Statistics.std(b::MCVIBelief, μ::Float64=mean(b))
    σ = 0.0
    for i in 1:length(b.weights)
        r = particle(b.space, i) - μ
        σ += b.weights[i] * (r * r)
    end
    return σ
end


text:
import CompatHelper
import Dates
import GitHub
import JSON
import Pkg
import Printf
import Test
import TimeZones

const timestamp_regex = r"integration\/(\d\d\d\d-\d\d-\d\d-\d\d-\d\d-\d\d-\d\d\d)\/"
const integration_regex = r"integration\/"
const compathelper_regex = r"compathelper\/"

function close_all_pull_requests(repo::GitHub.Repo;
                                 auth::GitHub.Authorization,
                                 state::String)
    all_pull_requests = CompatHelper.get_all_pull_requests(repo,
                                                           state;
                                                           auth = auth)
    for pr in all_pull_requests
        try
            GitHub.close_pull_request(repo, pr; auth = auth)
        catch
        end
    end
    return nothing
end

function delete_stale_branches(AUTOMERGE_INTEGRATION_TEST_REPO)
    with_cloned_repo(AUTOMERGE_INTEGRATION_TEST_REPO) do git_repo_dir
        cd(git_repo_dir)
        all_origin_branches = list_all_origin_branches(git_repo_dir)::Vector{String}
        for b in all_origin_branches
            if occursin(timestamp_regex, b) || occursin(integration_regex, b) || occursin(compathelper_regex, b)
                try
                    run(`git push origin --delete $(b)`)
                catch
                end
            end
        end
    end
    return nothing
end

function empty_git_repo(git_repo_dir::AbstractString)
    original_working_directory = pwd()
    cd(git_repo_dir)
    for x in readdir(git_repo_dir)
        if x != ".git"
            path = joinpath(git_repo_dir, x)
            rm(path; force = true, recursive = true)
        end
    end
    cd(original_working_directory)
    return nothing
end

function _generate_branch_name(name::AbstractString)
    sleep(0.1)
    _now = now_localzone()
    _now_utc_string = utc_to_string(_now)
    b = "integration/$(_now_utc_string)/$(rand(UInt32))/$(name)"
    sleep(0.1)
    return b
end

function generate_branch(name::AbstractString,
                         path_to_content::AbstractString,
                         parent_branch::AbstractString = "master";
                         repo_url)
    original_working_directory = pwd()
    b = _generate_branch_name(name)
    with_cloned_repo(repo_url) do git_repo_dir
        cd(git_repo_dir)
        run(`git checkout $(parent_branch)`)
        run(`git branch $(b)`)
        run(`git checkout $(b)`)
        empty_git_repo(git_repo_dir)
        for x in readdir(path_to_content)
            src = joinpath(path_to_content, x)
            dst = joinpath(git_repo_dir, x)
            rm(dst; force = true, recursive = true)
            cp(src, dst; force = true)
        end
        cd(git_repo_dir)
        try
            run(`git add -A`)
        catch
        end
        try
            run(`git commit -m "Automatic commit - CompatHelper integration tests"`)
        catch
        end
        try
            run(`git push origin $(b)`)
        catch
        end
        cd(original_working_directory)
        rm(git_repo_dir; force = true, recursive = true)
    end
    return b
end

function generate_master_branch(path_to_content::AbstractString,
                                parent_branch::AbstractString = "master";
                                repo_url)
    name = "master"
    b = generate_branch(name, path_to_content, parent_branch; repo_url = repo_url)
    return b
end

function generate_feature_branch(path_to_content::AbstractString,
                                 parent_branch::AbstractString;
                                 repo_url)
    name = "feature"
    b = generate_branch(name,
                        path_to_content,
                        parent_branch;
                        repo_url = repo_url)
    return b
end

function get_git_current_head(dir)
    original_working_directory = pwd()
    cd(dir)
    result = convert(String, strip(read(`git rev-parse HEAD`, String)))::String
    cd(original_working_directory)
    return result
end

function list_all_origin_branches(git_repo_dir)
    result = Vector{String}(undef, 0)
    original_working_directory = pwd()
    cd(git_repo_dir)
    a = try
        read(`git branch -a`, String)
    catch
        ""
    end
    b = split(strip(a), '\n')
    b_length = length(b)
    c = Vector{String}(undef, b_length)
    for i = 1:b_length
        c[i] = strip(strip(strip(b[i]), '*'))
        c[i] = first(split(c[i], "->"))
        c[i] = strip(c[i])
    end
    my_regex = r"^remotes\/origin\/(.*)$"
    for i = 1:b_length
        if occursin(my_regex, c[i])
            m = match(my_regex, c[i])
            if m[1] != "HEAD"
                push!(result, m[1])
            end
        end
    end
    cd(original_working_directory)
    return result
end

@inline now_localzone() = TimeZones.now(TimeZones.localzone())

function templates(parts...)
    this_filename = @__FILE__
    test_directory = dirname(this_filename)
    templates_directory = joinpath(test_directory, "templates")
    result = joinpath(templates_directory, parts...)
    return result
end

function username(auth::GitHub.Authorization)
    user_information = GitHub.gh_get_json(GitHub.DEFAULT_API,
                                          "/user";
                                          auth = auth)
    return user_information["login"]::String
end

@inline function utc_to_string(zdt::TimeZones.ZonedDateTime)
    zdt_as_utc = TimeZones.astimezone(zdt, TimeZones.tz"UTC")
    year = TimeZones.Year(zdt_as_utc.utc_datetime).value
    month = TimeZones.Month(zdt_as_utc.utc_datetime).value
    day = TimeZones.Day(zdt_as_utc.utc_datetime).value
    hour = TimeZones.Hour(zdt_as_utc.utc_datetime).value
    minute = TimeZones.Minute(zdt_as_utc.utc_datetime).value
    second = TimeZones.Second(zdt_as_utc.utc_datetime).value
    millisecond = TimeZones.Millisecond(zdt_as_utc.utc_datetime).value
    result = Printf.@sprintf "%04d-%02d-%02d-%02d-%02d-%02d-%03d" year month day hour minute second millisecond
    return result
end

function with_cloned_repo(f, repo_url)
    original_working_directory = pwd()
    result = with_temp_dir() do dir
        git_repo_dir = joinpath(dir, "REPO")
        cd(dir)
        try
            run(`git clone $(repo_url) REPO`)
        catch
        end
        cd(git_repo_dir)
        return f(git_repo_dir)
    end
    cd(original_working_directory)
    return result
end

function with_feature_branch(f::Function,
                             path_to_content::AbstractString,
                             parent_branch::AbstractString;
                             repo_url)
    b = generate_feature_branch(path_to_content,
                                parent_branch;
                                repo_url = repo_url)
    result = f(b)
    return result
end

function with_master_branch(f::Function,
                            path_to_content::AbstractString,
                            parent_branch::AbstractString;
                            repo_url)
    b = generate_master_branch(path_to_content,
                               parent_branch;
                               repo_url = repo_url)
    result = f(b)
    return result
end

function with_temp_dir(f)
    original_working_directory = pwd()
    tmp_dir = mktempdir()
    atexit(() -> rm(tmp_dir; force = true, recursive = true))
    cd(tmp_dir)
    result = f(tmp_dir)
    cd(original_working_directory)
    rm(tmp_dir; force = true, recursive = true)
    return result
end


text:
"""
    allocate_jacobian(x, y, sparsity=DensePattern())

Allocates a jacobian matrix given the inputs, outputs, and (optionally) the
sparsity structure.
"""
allocate_jacobian

allocate_jacobian(x, y) = allocate_jacobian(x, y, DensePattern())
allocate_jacobian(x, y, ::DensePattern) = NaN .* y .* x'
function allocate_jacobian(x, y, sp::SparsePattern)
    TF = promote_type(Float64, eltype(x), eltype(y))
    return sparse(sp.rows, sp.cols, fill(TF(NaN), length(sp.rows)))
end

"""
    ExplicitComponent([func,] fin, fout, foutin; kwargs...)

Construct a system component defined by the explicit function `func` with inputs
corresponding to `fin`, outputs corresponding to `fout`, and in-place
outputs corresponding to `foutin`.

# Arguments
 - `func`: Output function, of the form `(fout..., foutin...) = func(foutin..., fin...)`
 - `fin`: Tuple of named variables (see [`NamedVar`](@ref)) corresponding to inputs to `func`
 - `fout`: Tuple of named variables (see [`NamedVar`](@ref)) corresponding to outputs from `func`
 - `foutin`: Tuple of named variables (see [`Namedvar`](@ref)) corresponding to in-place function outputs from `func`

# Keyword Arguments
 - `f`: In-place output function `f(y, x)`. `x` is an input vector containing values
    corresponding to all of the variables in `component_inputs`. `y` is an output
    vector containing values corresponding to all of the variables in `component_outputs`.
 - `df`: In-place jacobian function `df(dydx, x)`.
 - `fdf`: In-place combined output and jacobian function `fdf(y, dydx, x)`.
 - `dydx`: Matrix used to define size and type of the jacobian matrix. If omitted,
    the size and type of the jacobian is inferred.
 - `deriv`: Method used to calculate the jacobian if both `df` and `fdf` are not
    provided.
 - `sparsity = DensePattern()`: Defines the sparsity structure of the jacobian
    if `dydx` is not provided
 - `component_inputs`: Names of the inputs from `fin` which are also used as
    inputs to `f`. Defaults to all variables in `inputs`.
 - `component_outputs`: Names of the outputs from `fout` and `foutin` which are
    also used as outputs from `f`. Defaults to all variables in `fout` and `foutin`.
"""
function ExplicitComponent(func, fin, fout, foutin;
    f=nothing, df=nothing, fdf=nothing, dydx = nothing, deriv=ForwardFD(),
    sparsity = DensePattern(),
    component_inputs = name.(fin),
    component_outputs = (name.(fout)..., name.(foutin)...))

    # ensure the output function is defined
    @assert any((!isnothing(func), !isnothing(f), !isnothing(fdf))) "Output function not defined"

    # combine inputs and outputs
    combined_inputs = (fin...,)
    combined_outputs = (fout..., foutin...)

    # get default inputs and outputs
    default_inputs = value.(combined_inputs)
    default_outputs = value.(combined_outputs)

    # check that specified inputs/outputs correspond to actual inputs/outputs
    @assert all(in(name.(combined_inputs)), component_inputs)
    @assert all(in(name.(combined_outputs)), component_outputs)

    # get indices of component inputs (as a tuple)
    component_input_indices = (
        findall(in(component_inputs), name.(fin))...,
        )

    # get indices of component outputs (as a tuple)
    component_unmut_output_indices = (
        findall(in(component_outputs), name.(fout))...,
        )
    component_mut_output_indices = (
        length(fout).+findall(in(component_outputs), name.(foutin))...,
        )
    #alternatively, we could just do:
    # component_unmut_output_indices = 1:length(fout)
    # component_unmut_output_indices = length(fout)+1:length(component_outputs)
    component_output_indices = (component_unmut_output_indices...,component_mut_output_indices...)

    # named variables corresponding to component inputs/outputs
    argin = getindices(combined_inputs, component_input_indices) #---> this is always the same as combined_inputs
    argout = getindices(combined_outputs, component_output_indices) #---> this is always the same as combined_outputs

    # construct input and output vectors
    x0 = combine(argin)
    y0 = combine(argout)

    # ensure sizes of `x`, `y`, and `dydx` are compatibile
    @assert isnothing(dydx) || (length(y0) == size(dydx, 1) && length(x0) == size(dydx, 2))

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* x0
    x_df = NaN .* x0
    y = NaN .* y0
    dydx = isnothing(dydx) ? allocate_jacobian(x0, y0, sparsity) : NaN .* dydx

    # construct output function (if necessary)
    if isnothing(f)
        if !isnothing(func)
            # NOTE: Any changes to this function need to be tested for avoiding allocations
            # use `func` to construct `f`
            f = function(y, x)
                # get new component inputs/outputs
                # Going from input vectors to tuples of arguments. The arguments in the tuples are views to the vector elements.
                new_component_inputs = separate(argin, x)
                new_component_outputs = separate(argout, y)

                # replace default inputs with new component inputs/outputs
                #  Note: Could we get rid of this? Seems like new_inputs = new_component_inputs, new_combined_outputs=new_component_outputs
                new_inputs = setindices(default_inputs, new_component_inputs, component_input_indices)
                new_combined_outputs = setindices(default_outputs, new_component_outputs, component_output_indices)

                # separate not-in-place and in-place outputs
                new_inplace = new_combined_outputs[length(fout) + 1 : end]
                # call function to get new outputs
                new_outputs = Tuple(func(new_inplace..., new_inputs...))
                # update output vector: replace the first elements
                combine!(y, getindices(new_outputs, component_unmut_output_indices))

                # return result
                return y
            end
        else
            # use `fdf` to construct `f`
            let dydx = copy(dydx)
                f = function(y, x)
                    fdf(y, dydx, x)
                    return y
                end
            end
        end
    end

    # construct jacobian functions (if necessary)
    if isnothing(df) && isnothing(fdf)
        # construct df and fdf functions from f
        df, fdf = create_output_jacobian_functions(f, x0, y0, deriv, sparsity)
    elseif isnothing(df)
        # construct df function from fdf
        let y = copy(y)
            df = function(dydx, x)
                fdf(y, dydx, x)
                return dydx
            end
        end
    elseif isnothing(fdf)
        # construct fdf function from f and df
        fdf = function(y, dydx, x)
            f(y, x)
            df(dydx, x)
            return y, dydx
        end
    end

    return ExplicitComponent(f, df, fdf, x_f, x_df, y, dydx, argin, argout)
end

ExplicitComponent(fin, fout, foutin; kwargs...) = ExplicitComponent(nothing, fin, fout, foutin; kwargs...)

"""
    ExplicitComponent(component::AbstractImplicitComponent [, fout]; kwargs...)

Couple an implicit component with a solver to construct an explicit component.

 # Arguments:
 - `component`: Implicit component or system to convert to an explicit component
 - `fout`: (optional) Tuple of named variables (see [`NamedVar`](@ref))
    corresponding to outputs from the resulting explicit component

 # Keyword Arguments
 - `solver = Newton()`: Solver, either a function of the form `y = f(component, x, y0)`
        or an object of type [`AbstractSolver`](@ref)
 - `dydx`: Provides size and type of jacobian output, otherwise it will be
        allocated based on the inputs, outputs, and the sparsity structure in `sparsity`
 - `sparsity = DensePattern()`: Defines the sparsity structure of the jacobian if `dydx` is not provided
 - `mode`: Mode in which to compute the analytic sensitivity equations.  May be
    either [`Direct()`](@ref) or [`Adjoint()`](@ref).  Defaults to `Direct()` if
    the number of inputs is less than the number of outputs and `Adjoint()`
    otherwise.
"""
function ExplicitComponent(component::AbstractImplicitComponent, argout=component.argout; solver=Newton(),
    dydx = nothing, sparsity=DensePattern(), mode=nothing)

    # input variables
    argin = component.argin

    # state variables
    argstate = component.argout

    # initial guess comes from default values for state variables
    u0 = combine(argstate)

    # construct dfdx and dfdu
    dfdx, dfdu = output_mapping_matrices(argin, argstate, argout)

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* inputs(component) # inputs are the same as component
    x_df = NaN .* inputs(component) # inputs are the same as component
    y = NaN .* combine(argout) # outputs correspond to specified outputs
    dydx = isnothing(dydx) ? allocate_jacobian(x_f, y, sparsity) : NaN .* dydx

    # choose mode in which to apply analytic sensitivity equations (if not specified)
    if isnothing(mode)
        mode = ifelse(length(x_f) < length(y), Direct(), Adjoint())
    end

    # NOTE: xcache stores the current inputs to the solver
    # As xcache is updated, the function upon which the solver operates is
    # updated as well.
    xcache = similar(x_f)

    # we distinguish between mutating and non-mutating cache so that we can use
    # automatic differentiation on the non-mutating cache.
    cache = create_mutating_solver_cache(component, xcache, solver)

    # output function
    f = function(y, x)
        # determine whether to use internal cache
        if eltype(x) <: eltype(xcache)
            # use existing cache since we can store the inputs
            copyto!(xcache, x)
            # call solver to get state variables
            u = solver(cache, x, u0)
        else
            # we can't use the existing cache, create a new cache
            new_cache = create_nonmutating_solver_cache(component, x, solver)
            # call solver to get state variables
            u = solver(new_cache, x, eltype(y).(u0))
        end
        # map inputs and state variables to outputs: y = dfdx*x + dfdu*u
        y = mul!(mul!(y, dfdx, x), dfdu, u, 1, 1)
        return y
    end

    # function and partials
    fdf = function(y, dydx, x)
        # determine whether to use internal cache
        if eltype(x) <: eltype(xcache)
            # use existing cache if we can store the inputs
            copyto!(xcache, x)
            # call solver to get state variables
            u = solver(cache, x, u0)
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian!(component, x, u)
            # get jacobian of the residual function with respect to the outputs
            drdu = residual_output_jacobian!(component, x, u)
            # map inputs and state variables to outputs: y = dfdx*x + dfdu*u
            y = mul!(mul!(y, dfdx, x), dfdu, u, 1, 1)
        else
            # we can't use the existing cache, create a new cache
            new_cache = create_nonmutating_solver_cache(component, x, solver)
            # call solver to get state variables
            u = solver(new_cache, x, eltype(y).(u0))
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian(component, x, u)
            # get jacobian of the residual function with respect to the outputs
            drdu = residual_output_jacobian(component, x, u)
            # map inputs and state variables to outputs: y = dfdx*x + dfdu*u
            y = mul!(mul!(y, dfdx, x), dfdu, u, 1, 1)
        end
        # get analytic sensitivities
        copyto!(dydx, analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, mode))
        # return outputs
        return y, dydx
    end

    # partials only
    df = function(dydx, x)
        y, dydx = fdf(y, dydx, x)
        return dydx
    end

    return ExplicitComponent(f, df, fdf, x_f, x_df, y, dydx, argin, argout)
end

"""
    ExplicitComponent(component::AbstractImplicitComponent, output_component; kwargs...)

Couple an implicit component with a solver to construct an explicit component.

Reduce the dimensionality of the outputs by incorporating the output function
contained in `output_component`.  The vector of inputs to `output_component` are
assumed to contain the input variables to `component` followed by the
output/state variables of `component` .

 # Keyword Arguments
 - `solver = Newton()`: Solver, either a function of the form `y = f(component, x, y0)`
        or an object of type [`AbstractSolver`](@ref)
 - `dydx`: Provides size and type of jacobian output, otherwise it will be
        allocated based on the inputs, outputs, and the sparsity structure in `sparsity`
 - `sparsity = DensePattern()`: Defines the sparsity structure of the jacobian if `dydx` is not provided
 - `mode`: Mode in which to compute the analytic sensitivity equations.  May be
    either [`Direct()`](@ref) or [`Adjoint()`](@ref).  Defaults to `Direct()` if
    the number of inputs is less than the number of outputs and `Adjoint()`
    otherwise.
"""
function ExplicitComponent(component::AbstractImplicitComponent, output_component::AbstractExplicitComponent;
    solver=Newton(), u0=rand(length(component.y_f)), dydx = nothing,
    sparsity=DensePattern(), mode = nothing)

    # input variables
    argin = component.argin

    # state variables
    argstate = component.argout

    # output variables
    argout = output_component.argout

    # initial guess comes from default values for state variables
    u0 = combine(argstate)

    # pre-allocated storage vector for inputs to output function
    v = combine(output_component.argin)

    # pre-allocated storage matrices for jacobians of output function
    dfdx = outputs(output_component) * inputs(component)' .* NaN
    dfdu = outputs(output_component) * outputs(component)' .* NaN

    # construct dvdx and dvdu where `v` is a vector of inputs to `output_component`
    dvdx, dvdu = output_mapping_matrices(argin, argstate, output_component.argin)

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* inputs(component)
    x_df = NaN .* inputs(component)
    y = NaN .* outputs(output_component)
    dydx = isnothing(dydx) ? allocate_jacobian(x_f, y, sparsity) : NaN .* dydx

    # choose mode in which to apply analytic sensitivity equations (if not specified)
    if isnothing(mode)
        mode = ifelse(length(x_f) < length(y), Direct(), Adjoint())
    end

    # input and state variable dimensions
    nx = length(component.x_f)
    nu = length(component.y_f)

    # NOTE: xcache stores the current inputs to the solver
    # As xcache is updated, the function upon which the solver operates is
    # updated as well.
    xcache = similar(x_f)

    # we distinguish between mutating and non-mutating cache so that we can
    # use different types when using the non-mutating cache
    cache = create_mutating_solver_cache(component, xcache, solver)

    # output function
    f = function(y, x)
        if eltype(x) <: eltype(xcache)
            # use existing cache since we can store the inputs
            copyto!(xcache, x)
            # call solver to get state variables
            u = solver(cache, x, u0)
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = mul!(mul!(v, dvdx, x), dvdu, u, 1, 1)
            # get outputs from output function
            copyto!(y, outputs!(output_component, v))
        else
            # we can't use the existing cache, create a new cache
            new_cache = create_nonmutating_solver_cache(component, x, solver)
            # call solver to get state variables
            u = solver(new_cache, x, eltype(y).(u0))
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = dvdx*x + dvdu*u
            # get outputs from output function
            copyto!(y, outputs(output_component, v))
        end
        # return result
        return y
    end

    # partials
    df = function(dydx, x)
        if eltype(x) <: eltype(xcache)
            # use existing cache since we can store the inputs
            copyto!(xcache, x)
            # call solver to get state variables
            u = solver(cache, x, u0)
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = mul!(mul!(v, dvdx, x), dvdu, u, 1, 1)
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian!(component, x, u)
            # get jacobian of the residual function with respect to the state variables
            drdu = residual_output_jacobian!(component, x, u)
            # get jacobian of the output component with respect to its inputs
            dfdv = jacobian!(output_component, v)
            # extract output jacobian corresponding to inputs
            dfdx = mul!(dfdx, dfdv, dvdx)
            # extract output jacobian corresponding to state variables
            dfdu = mul!(dfdu, dfdv, dvdu)
        else
            # we can't use the existing cache, create a new cache
            new_cache = create_nonmutating_solver_cache(component, x, solver)
            # u0 might not have the right type so we convert it
            u = solver(new_cache, x, eltype(x).(u0))
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian(component, x, u)
            # get jacobian of the residual function with respect to the state variables
            drdu = residual_output_jacobian(component, x, u)
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = dvdx*x + dvdu*u
            # get jacobian of the output component with respect to its inputs
            dfdv = jacobian(output_component, v)
            # extract output jacobian corresponding to inputs
            dfdx = dfdv*dvdx
            # extract output jacobian corresponding to state variables
            dfdu = dfdv*dvdu
        end
        # apply analytic sensitivity equation in direct or adjoint mode
        copyto!(dydx, analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, mode))
        return dydx
    end

    fdf = function(y, dydx, x)
        if eltype(x) <: eltype(xcache)
            # use existing cache since we can store the inputs
            copyto!(xcache, x)
            # call solver to get state variables
            u = solver(cache, x, u0)
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = mul!(mul!(v, dvdx, x), dvdu, u, 1, 1)
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian!(component, x, u)
            # get jacobian of the residual function with respect to the state variables
            drdu = residual_output_jacobian!(component, x, u)
            # get jacobian of the output component with respect to its inputs
            y[:], dfdv = outputs_and_jacobian!(output_component, v)
            # extract output jacobian corresponding to inputs
            dfdx = mul!(dfdx, dfdv, dvdx)
            # extract output jacobian corresponding to state variables
            dfdu = mul!(dfdu, dfdv, dvdu)
        else
            # we can't use the existing cache, create a new cache
            new_cache = create_nonmutating_solver_cache(component, x, solver)
            # u0 might not have the right type so we convert it
            u = solver(new_cache, x, eltype(x).(u0))
            # get jacobian of the residual function with respect to the inputs
            drdx = residual_input_jacobian(component, x, u)
            # get jacobian of the residual function with respect to the state variables
            drdu = residual_output_jacobian(component, x, u)
            # map inputs and state variables to output function inputs: v = dvdx*x + dvdu*u
            v = dvdx*x + dvdu*u
            # get jacobian of the output component with respect to its inputs
            y[:], dfdv = outputs_and_jacobian(output_component, v)
            # extract output jacobian corresponding to inputs
            dfdx = dfdv*dvdx
            # extract output jacobian corresponding to state variables
            dfdu = dfdv*dvdu
        end
        # apply analytic sensitivity equation in direct or adjoint mode
        copyto!(dydx, analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, mode))
        return y, dydx
    end

    return ExplicitComponent(f, df, fdf, x_f, x_df, y, dydx, argin, argout)
end

"""
    analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, mode)

Compute the sensitivities of an implicit system using the analytic sensitivity
equation in `Direct()` or `Adjoint()` mode.
"""
analytic_sensitivity_equation
analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, ::Direct) = dfdx - dfdu*(drdu\drdx)
analytic_sensitivity_equation(dfdx, dfdu, drdx, drdu, ::Adjoint) = dfdx - transpose(transpose(drdu)\transpose(dfdu))*drdx

"""
    ExplicitSystem(components, argin, argout; kwargs...)

Construct an explicit system component from a collection of explicit system
components.

# Arguments
- `components`: Collection of components
- `argin`: Tuple of system variables (see [`NamedVar`](@ref)) corresponding to system inputs
- `argout`: Tuple of system variables (see [`NamedVar`](@ref)) corresponding to system outputs

# Keyword Arguments
 - `dydx`: Matrix used to define size and type of the jacobian matrix. If omitted,
    the jacobian matrix size and type will be infered from the inputs, outputs, and
    the sparsity structure.
 - `sparsity`: Sparsity structure of the system jacobian matrix
 - `mode`: Mode used to calculate system derivatives using the chain rule.  May
    be either [`Forward()`](@ref) or [`Reverse()`](@ref), defaults to
    [`Forward`](@ref) if the number of outputs exceeds the number of inputs and
    [`Reverse()`](@ref) otherwise.
"""
function ExplicitSystem(components, argin, argout;
    dydx = nothing, sparsity=DensePattern(), mode=nothing)

    # number of components
    nc = length(components)

    # construct input and output vectors
    x0 = combine(argin)
    y0 = combine(argout)

    # choose mode in which to apply analytic sensitivity equations (if not specified)
    if isnothing(mode)
        mode = ifelse(length(x0) < length(y0), Forward(), Reverse())
    end

    # construct system mapping
    component_input_mapping = system_component_mapping(components, argin)

    # construct output mapping
    output_mapping = system_output_mapping(components, argin, argout)

    # ensure all components are explicit
    @assert all(isa.(components, AbstractExplicitComponent)) "All components of an explicit system must be explicit"

    # ensure sizes of `x`, `y`, and `dydx` are compatibile
    @assert isnothing(dydx) || (length(y0) == size(dydx, 1) && length(x0) == size(dydx, 2))

    # check that there is an input mapping for each component
    @assert length(components) == length(component_input_mapping) "There must be an input mapping for each component"
    for ic = 1:nc
        # check that there is an input mapping for each input to the component
        @assert length(components[ic].x_f) == length(component_input_mapping[ic]) "Input mapping for component $ic must have the same length as the number of inputs to component $ic"
        for ix = 1:length(component_input_mapping[ic])
            (jc, jy) = component_input_mapping[ic][ix]
            if iszero(jc)
                # check that component input uses available system input
                @assert jy <= length(x0) "Input $ix of "*
                    "component $ic corresponds to non-existant system input $jy"
            else
                # check that each component input only uses available component outputs
                @assert jc <= nc "Input $ix of component $ic comes from "*
                    "non-existant component $jc"
                @assert jc < ic "Implicit dependency in explicit system: Input $ix "*
                    "of component $ic depends on the output of component $jc"
                @assert jy <= length(outputs(components[jc])) "Input $ix of component $ic "*
                    "corresponds to non-existant output $jx of component $jc"
            end
        end
        for iy = 1:length(output_mapping)
            (jc, jy) = output_mapping[iy]
            if iszero(jc)
                # check that system output uses available system inputs
                @assert jy <= length(x0) "System output $iy "*
                    "corresponds to non-existant system input $jy"
            else
                # check that each component input only uses available component outputs
                @assert jc <= nc "System output $iy comes from non-existant "
                    "component $jc"
                @assert jy <= length(outputs(components[jc])) "System output $iy "*
                    "corresponds to non-existant output $jy of component $jc"
            end
        end
    end

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* x0
    x_df = NaN .* x0
    y = NaN .* y0
    dydx = isnothing(dydx) ? allocate_jacobian(x0, y0, sparsity) : NaN .* dydx

    # construct forward mode mapping
    input_mapping, component_output_mapping = forward_mode_mapping(x0,
        components, component_input_mapping, output_mapping)

    return ExplicitSystem(components, input_mapping, component_output_mapping,
        component_input_mapping, output_mapping, x_f, x_df, y, dydx, argin, argout, mode)
end

"""
    forward_mode_mapping(inputs, components, component_input_mapping, output_mapping)

Returns (`input_mapping, component_output_mapping`) where `input_mapping` contains
a mapping from each system input to component inputs and/or system outputs and
`reversed_component_mapping` contains a mapping from each component output to
corresponding component inputs and/or system outputs.
"""
function forward_mode_mapping(inputs, components, component_input_mapping, output_mapping)

    # get number of inputs and components
    nx = length(inputs)
    nc = length(components)

    # initialize input mapping
    input_mapping = Vector{NTuple{2, Vector{Int}}}(undef, nx)
    # populate with empty arrays
    for ix = 1:nx
        input_mapping[ix] = (Int[], Int[])
    end

    # initialize reversed component mapping (outputs to inputs)
    component_output_mapping = Vector{Vector{NTuple{2, Vector{Int}}}}(undef, nc)
    for ic = 1:nc
        # number of outputs from this component
        ncy = length(outputs(components[ic]))
        # initialize output to input mapping for this component
        component_output_mapping[ic] = Vector{NTuple{2, Vector{Int}}}(undef, ncy)
        # populate with empty arrays
        for icy = 1:ncy
            component_output_mapping[ic][icy] = (Int[], Int[])
        end
    end

    # add mappings from component_mapping
    for jc = 1:nc
        # number of inputs from this component
        ncx = length(component_input_mapping[jc])
        for jcx = 1:ncx
            ic, icy = component_input_mapping[jc][jcx]
            if iszero(ic)
                ix = icy
                push!(input_mapping[ix][1], jc)
                push!(input_mapping[ix][2], jcx)
            else
                push!(component_output_mapping[ic][icy][1], jc)
                push!(component_output_mapping[ic][icy][2], jcx)
            end
        end
    end

    # add mappings from output_mapping
    ny = length(output_mapping)
    for iy = 1:ny
        ic, icy = output_mapping[iy]
        if iszero(ic)
            ix = icy
            push!(input_mapping[ix][1], nc+1)
            push!(input_mapping[ix][2], iy)
        else
            push!(component_output_mapping[ic][icy][1], nc+1)
            push!(component_output_mapping[ic][icy][2], iy)
        end
    end

    return input_mapping, component_output_mapping
end

"""
    ImplicitComponent([func,] fin, fout [, r0]; kwargs...)

Construct a system component defined by the in-place implicit residual function
`func` with inputs corresponding to `fin` and outputs corresponding to `fout`.
The output variables are identical to the state variables of the component.

# Arguments
 - `func`: In-place residual function, of the form `residuals = func(r, fin..., fout...)`
 - `fin`: Tuple of named variables (see [`NamedVar`](@ref)) corresponding to function inputs
 - `fout`: Tuple of named variables (see [`NamedVar`](@ref)) corresponding to function outputs
 - `r0`: Vector which defines the size and type of the component's residuals.
    Defaults to the size and type of the output vector. It's length must
    correspond to the total number of component outputs/states.

# Keyword Arguments
 - `f`: In-place residual function `f(r, x, y)`. `x` is an input vector containing values
    corresponding to all of the variables in `component_inputs`. `y` is an output
    vector containing values corresponding to all of the variables in
    `component_outputs`.
 - `dfdx`: In-place residual jacobian function with respect to the inputs `dfdx(drdx, x, y)`
 - `dfdy`: In-place residual jacobian function with respect to the outputs `dfdy(drdy, x, y)`
 - `fdfdx`: In-place combined residual and jacobian with respect to the inputs function `fdfdx(r, drdx, x, y)`.
 - `fdfdy`: In-place combined residual and jacobian with respect to the outputs function `fdfdy(r, drdy, x, y)`.
 - `fdf`: In-place combined residual and jacobians function `fdf(r, drdx, drdy, x, y)`.
 - `drdx`: Matrix used to define size and type of the jacobian of the residual
    with respect to the inputs. If omitted, a dense jacobian is used.
 - `drdy`: Matrix used to define size and type of the jacobian of the residual
     with respect to the outputs. If omitted, a dense jacobian is used.
 - `xderiv`: Method used to calculate the residual jacobian with respect to the
    inputs if it is not provided. Defaults to [`ForwardFD`](@ref).
 - `yderiv`: Method used to calculate the residual jacobian with respect to the
    outputs if it is not provided. Defaults to [`ForwardFD`](@ref).
 - `xsparsity = DensePattern()`: Defines the sparsity structure of the jacobian
    with respect to the inputs if `drdx` is not provided
 - `ysparsity = DensePattern()`: Defines the sparsity structure of the jacobian
    with respect to the outputs if `drdy` is not provided
 - `component_inputs`: Names of the inputs from `fin` which are also used as
    inputs to `f`. Defaults to all variables in `fin`.
 - `component_outputs`: Names of the outputs from `fout` which are also used
    as outputs from `f`. Defaults to all variables in `fout`.
"""
function ImplicitComponent(func, fin, fout, r0=nothing;
    f=nothing, dfdx=nothing, dfdy=nothing, fdfdx=nothing, fdfdy=nothing, fdf=nothing,
    drdx=nothing, drdy=nothing, xderiv=ForwardFD(), yderiv=ForwardFD(),
    xsparsity = DensePattern(), ysparsity = DensePattern(),
    component_inputs = name.(fin), component_outputs = name.(fout))

    # ensure the residual function is defined
    @assert any((!isnothing(func), !isnothing(f), !isnothing(fdfdx), !isnothing(fdfdy), !isnothing(fdf))) "Residual function not defined"

    # get default inputs and outputs
    default_inputs = value.(fin)
    default_outputs = value.(fout)

    # check that specified argin/argout correspond to actual argin/argout
    @assert all(in(name.(fin)), component_inputs)
    @assert all(in(name.(fout)), component_outputs)

    # get indices of component inputs (as a tuple)
    component_input_indices = (
        findall(in(component_inputs), name.(fin))...,
        )

    # get indices of component outputs (as a tuple)
    component_output_indices = (
        findall(in(component_outputs), name.(fout))...,
        )

    # named variables corresponding to component inputs/outputs
    argin = getindices(fin, component_input_indices)
    argout = getindices(fout, component_output_indices)

    # construct input and output vectors
    x0 = combine(argin)
    y0 = combine(argout)
    r0 = isnothing(r0) ? similar(y0, promote_type(eltype(x0), eltype(y0))) : r0

    # ensure sizes of `r` and `y` are compatible
    @assert length(r0) == length(y0) "The length of the residual vector must match the number of outputs/states"

    # ensure sizes of `x`, `r`, and `drdx` are compatibile
    @assert isnothing(drdx) || (length(r0) == size(drdx, 1) && length(x0) == size(drdx, 2))

    # ensure sizes of `y`, `r`, and `drdy` are compatibile
    @assert isnothing(drdy) || (length(r0) == size(drdy, 1) && length(y0) == size(drdy, 2))

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* x0
    y_f = NaN .* y0
    x_dfdx = NaN .* x0
    y_dfdx = NaN .* y0
    x_dfdy = NaN .* x0
    y_dfdy = NaN .* y0
    r = NaN .* r0
    drdx = isnothing(drdx) ? allocate_jacobian(x0, r0) : drdx .* NaN
    drdy = isnothing(drdy) ? allocate_jacobian(y0, r0) : drdy .* NaN

    # construct residual function if necessary
    if isnothing(f)
        if !isnothing(func)
            # construct residual function from `func`
            f = function(r, x, y)
                # get new component inputs/outputs
                new_component_inputs = separate(argin, x)
                new_component_outputs = separate(argout, y)
                # replace default inputs with new component inputs/outputs
                new_inputs = setindices(default_inputs, new_component_inputs, component_input_indices)
                new_outputs = setindices(default_outputs, new_component_outputs, component_output_indices)
                # call function to update residual values
                func(r, new_inputs..., new_outputs...)
                # return result
                return r
            end
        elseif !isnothing(fdfdx)
            # construct residual function from fdfdx
            let drdx = copy(drdx)
                f = function(r, x, y)
                    fdfdx(r, drdx, x, y)
                    return y
                end
            end
        elseif !isnothing(fdfdy)
            # construct residual function from fdfdy
            let drdy = copy(drdy)
                f = function(r, x, y)
                    fdf(r, drdy, x, y)
                    return y
                end
            end
        else
            # construct residual function from fdf
            let drdx = copy(drdx), drdy = copy(drdy)
                f = function(r, x, y)
                    fdf(r, drdx, drdy, x, y)
                    return y
                end
            end
        end
    end

    # construct jacobian functions with respect to the inputs
    if isnothing(dfdx) && isnothing(fdfdx) && isnothing(fdf)
        # construct dfdx and fdfdx from f
        dfdx, fdfdx = create_residual_input_jacobian_functions(f, x_f, y_f, r, xderiv, xsparsity)
    elseif isnothing(dfdx) && isnothing(fdfdx)
        # construct dfdx from fdf
        let r = copy(r), drdy = copy(drdy)
            dfdx = function(drdx, x, y)
                fdf(r, drdx, drdy, x, y)
                return drdx
            end
        end

        # construct fdfdx from fdf
        let drdy = copy(drdy)
            fdfdx = function(r, drdx, x, y)
                fdf(r, drdx, drdy, x, y)
                return r, drdx
            end
        end
    elseif isnothing(dfdx)
        # construct dfdx from fdfdx
        let r = copy(r)
            dfdx = function(drdx, x, y)
                fdfdx(r, drdx, x, y)
                return drdx
            end
        end
    elseif isnothing(fdfdx)
        # construct fdfdx from f and dfdx
        fdfdx = function(r, drdx, x, y)
            f(r, x, y)
            dfdx(drdx, x, y)
            return r, drdx
        end
    end

    # construct jacobian functions with respect to the outputs
    if isnothing(dfdy) && isnothing(fdfdy) && isnothing(fdf)
        # construct dfdy and fdfdy from f
        dfdy, fdfdy = create_residual_output_jacobian_functions(f, x_f, y_f, r, yderiv, ysparsity)
    elseif isnothing(dfdy) && isnothing(fdfdy)
        # construct dfdy from fdf
        let r = copy(r), drdx = copy(drdx)
            dfdy = function(drdy, x, y)
                fdf(r, drdx, drdy, x, y)
                return drdy
            end
        end

        # construct fdfdy from fdf
        let drdx = copy(drdx)
            fdfdy = function(r, drdy, x, y)
                fdf(r, drdx, drdy, x, y)
                return r, drdy
            end
        end
    elseif isnothing(dfdy)
        # construct dfdy from fdfdy
        let r = copy(r)
            dfdy = function(drdy, x, y)
                fdfdy(r, drdy, x, y)
                return drdy
            end
        end
    elseif isnothing(fdfdy)
        # construct fdfdy from f and dfdy
        fdfdy = function(r, drdy, x, y)
            f(r, x, y)
            dfdy(drdy, x, y)
            return r, drdy
        end
    end

    # construct combined output and jacobians function if necessary
    if isnothing(fdf)
        # construct fdf from dfdx and fdfdy
        fdf = function(r, drdx, drdy, x, y)
            dfdx(drdx, x, y)
            fdfdy(r, drdy, x, y)
            return r, drdx, drdy
        end
    end

    return ImplicitComponent(f, dfdx, dfdy, fdfdx, fdfdy, fdf, x_f, y_f, x_dfdx,
        y_dfdx, x_dfdy, y_dfdy, r, drdx, drdy, argin, argout)
end

ImplicitComponent(fin, fout, args...; kwargs...) = ImplicitComponent(nothing, fin, fout, args...; kwargs...)

make_implicit(component::ExplicitComponent) = ImplicitComponent(component)
make_implicit(component::ExplicitSystem) = ImplicitComponent(component)
make_implicit(component::ImplicitComponent) = component
make_implicit(component::ImplicitSystem) = component

"""
    ImplicitComponent(component::ExplicitComponent)

Constructs an implicit system component from an explicit system component.
"""
function ImplicitComponent(component::ExplicitComponent)

    # allocate storage and initialize with NaNs (since values are as of yet undefined)
    x_f = NaN .* component.x_f
    y_f = NaN .* component.y
    x_dfdx = NaN .* component.x_f
    y_dfdx = NaN .* component.y
    x_dfdy = NaN .* component.x_f
    y_dfdy = NaN .* component.y
    r = NaN .* component.y
    drdx = NaN .* component.dydx
    drdy = Matrix(-I, length(r), length(y_f))

    f = function(r, x, y)
        # update outputs, store in residual
        outputs!(component, r, x)
        # subtract provided outputs
        r .-= y
        # return residual
        return r
    end

    dfdx = (drdx, x, y) -> jacobian!(component, drdx, x)

    fdfdx = function(r, drdx, x, y)
        # update outputs, store in residual
        outputs_and_jacobian!(component, r, drdx, x)
        # subtract provided outputs
        r .-= y
        # return residual
        return r
    end

    dfdy = function(drdy, x, y)
        # negative identity matrix
        drdy .= 0
        for i = 1:length(y)
            drdy[i,i] = -1
        end
    end

    fdfdy = function(r, drdy, x, y)
        f(r, x, y)
        dfdy(drdy, x, y)
        return r, drdy
    end

    fdf = function(r, drdx, drdy, x, y)
        fdfdx(r, drdx, x, y)
        dfdy(drdy, x, y)
        return r, drdx, drdy
    end

    # input/output arguments are preserved
    argin = component.argin
    argout = component.argout

    return ImplicitComponent(f, dfdx, dfdy, fdfdx, fdfdy, fdf, x_f, y_f, x_dfdx,
        y_dfdx, x_dfdy, y_dfdy, r, drdx, drdy, argin, argout)
end

"""
    ImplicitSystem(components, argin; kwargs...)

Constructs an implicit system component from a collection of implicit components.

# Arguments
 - `components`: Collection of components
 - `argin`: Tuple of system variables (see [`NamedVar`](@ref)) corresponding to system inputs

# Keyword Arguments
- `drdx`: Matrix used to define size and type of the jacobian of the residual
   with respect to the inputs.  If omitted, the jacobian matrix size and type
   will be inferred from the inputs, outputs, and the sparsity structure.
- `drdy`: Matrix used to define size and type of the jacobian of the residual
    with respect to the outputs.  If omitted, the jacobian matrix size and type
    will be infered from the inputs, outputs, and the sparsity structure.
- `xsparsity`: Sparsity structure of the jacobian of the residual with respect
    to the inputs. Defaults to [`DensePattern()`](@ref).
- `ysparsity`: Sparsity struction of the jacobian of the residual with respect
    to the outputs.  Defaults to [`DensePattern()`](@ref).
"""
function ImplicitSystem(components, argin;
    drdx=nothing, drdy=nothing, xsparsity=DensePattern(), ysparsity=DensePattern())

    # number of components
    nc = length(components)

    # convert all components to implicit components
    components = make_implicit.(components)

    # construct input vector
    x0 = combine(argin)

    # construct system mapping
    component_input_mapping = system_component_mapping(components, argin)

    # check that there is an input mapping for each component
    @assert length(components) == length(component_input_mapping) "There must be an input mapping for each component"
    for ic = 1:nc
        # check that there is an input mapping for each input to the component
        @assert length(components[ic].x_f) == length(component_input_mapping[ic]) "Input mapping for component $ic must have the same length as the number of inputs to component $ic"
        for ix = 1:length(component_input_mapping[ic])
            (jc, jy) = component_input_mapping[ic][ix]
            if iszero(jc)
                @assert jy <= length(x0) "Input $ix of "*
                    "component $ic corresponds to non-existant system input $jy"
            else
                # check that each component input only uses available component outputs
                @assert jc <= nc "Input $ix of component $ic comes from "*
                    "non-existant component $jc"
                @assert jy <= length(outputs(components[jc])) "Input $ix of component $ic "*
                    "corresponds to non-existant output $jx of component $jc"
            end
        end
    end

    # allocate storage and initialize with NaNs (since values are not yet defined)
    x_f = x0 .* NaN
    y_f = vcat(outputs.(components)...) .* NaN
    x_dfdx = copy(x_f)
    y_dfdx = copy(y_f)
    x_dfdy = copy(x_f)
    y_dfdy = copy(y_f)
    r = copy(y_f)
    drdx = isnothing(drdx) ? allocate_jacobian(x_f, r, xsparsity) : drdx .* NaN
    drdy = isnothing(drdy) ? allocate_jacobian(y_f, r, ysparsity) : drdy .* NaN

    # ensure jacobian sizes are compatibile
    @assert (length(r) == size(drdx, 1) && length(x_f) == size(drdx, 2))
    @assert (length(r) == size(drdy, 1) && length(y_f) == size(drdy, 2))

    # system outputs are all of the component outputs concatenated
    argout = (Iterators.flatten([comp.argout for comp in components])...,)

    # construct index for accessing outputs/residuals for each component
    idx = vcat(0, cumsum(length.(residuals.(components)))...)

    return ImplicitSystem(Tuple(components), collect(component_input_mapping),
        x_f, y_f, x_dfdx, y_dfdx, x_dfdy, y_dfdy, r, drdx, drdy, argin, argout, idx)
end


text:
using Sigma

println("running examples")
include("../examples/runexamples.jl")

# tests = ["controlflow",
#          "hyperbox",
#          "dreal",
#          "distributions"]
#
# println("Running tests:")
#
# for t in tests
#     test_fn = "$t.jl"
#     println(" * $test_fn")
#     include(test_fn)
# end


text:
import UUIDs: UUID, uuid1

###
# CLIENT
###

mutable struct ClientSession
    id::Symbol
    stream::Any
    connected_notebook::Union{Notebook,Nothing}
    pendingupdates::Channel
end

ClientSession(id::Symbol, stream) = let
    ClientSession(id, stream, nothing, Channel(1024))
end

"A combination of _client ID_ and a _request ID_. The front-end generates a unqique ID for every request that it sends. The back-end (the stuff you are currently reading) can respond to a specific request. In that case, the response does not go through the normal message handlers in the front-end, but it flies directly to the place where the message was sent. (It resolves the promise returned by `send(...)`.)"
struct Initiator
    client_id::Symbol
    request_id::Symbol
end

###
# SERVER
###

"""
The `ServerSession` keeps track of:

- `connected_clients`: connected (web) clients
- `notebooks`: running notebooks
- `secret`: the web acces token
- `default_environment_path`: specifiy the project/package environment to start new notebooks in. Set to `nothing` to default to the global environment.
"""
Base.@kwdef struct ServerSession
    connected_clients::Dict{Symbol,ClientSession} = Dict{Symbol,ClientSession}()
    notebooks::Dict{UUID,Notebook} = Dict{UUID,Notebook}()
    secret::UUID = uuid1()
    default_environment_path::Union{Nothing,AbstractString} = nothing
end

###
# UPDATE MESSAGE
###

struct UpdateMessage
    type::Symbol
    message::Any
    notebook::Union{Notebook,Nothing}
    cell::Union{Cell,Nothing}
    initiator::Union{Initiator,Missing}
end

UpdateMessage(type::Symbol, message::Any) = UpdateMessage(type, message, nothing, nothing, missing)
UpdateMessage(type::Symbol, message::Any, notebook::Notebook) = UpdateMessage(type, message, notebook, nothing, missing)

function clientupdate_cell_output(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_output,
        Dict(
            :queued => cell.queued,
            :running => cell.running,
            :runtime => cell.runtime,
            :errored => cell.errored,
            :output => Dict(
                :mime => cell.repr_mime,
                :body => cell.output_repr,
                :rootassignee => cell.rootassignee,
            )
        ),
        notebook, cell, initiator)
end

function clientupdate_cell_input(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_input,
        Dict(:code => cell.code, :folded => cell.code_folded), notebook, cell, initiator)
end

function clientupdate_cell_added(notebook::Notebook, cell::Cell, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_added,
        Dict(
            :index => new_index - 1, # 1-based index (julia) to 0-based index (js)
        ), notebook, cell, initiator)
end

function clientupdate_cell_deleted(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_deleted,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cells_moved(notebook::Notebook, cells::Vector{Cell}, new_index::Integer; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cells_moved,
        Dict(
            :cells => [cell.cell_id for cell in cells],
            :index => new_index - 1, # 1-based index (julia) to 0-based index (js)
        ), notebook, nothing, initiator)
end

function clientupdate_cell_queued(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_queued,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cell_running(notebook::Notebook, cell::Cell; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_running,
        Dict(), notebook, cell, initiator)
end

function clientupdate_cell_folded(notebook::Notebook, cell::Cell, folded::Bool; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:cell_folded,
        Dict(:folded => folded), notebook, cell, initiator)
end

function clientupdate_notebook_list(notebooks; initiator::Union{Initiator,Missing}=missing)
    UpdateMessage(:notebook_list,
        Dict(
            :notebooks => [
                Dict(
                    :notebook_id => notebook.notebook_id,
                    :path => notebook.path,
                    :in_temp_dir => startswith(notebook.path, tempdir()),
                    :shortpath => basename(notebook.path)
                ) for notebook in values(notebooks)
            ]
        ), nothing, nothing, initiator)
end


text: ### Placeholder for merging chains


text:
######################################################################
# CancellableTasks UTs
# -----
# Licensed under MIT License
module TestCancellableTasks
using Test
using ExtraFun.CancellableTasks
import ExtraFun: CancellationError, TimeoutError, cancel

function task_fails(cb, Error::Type{<:Exception})
    try
        cb()
        return false
    catch ex
        @assert ex isa TaskFailedException
        return ex.task.result isa Error
    end
end
function task_fails(cb, err)
    try
        cb()
        return false
    catch ex
        @assert ex isa TaskFailedException
        return ex.task.result === err
    end
end

@testset "Cancel Tasks" begin
    @testset "Cancellable" begin
        @test fetch(with_cancel(() -> 42)) == 42

        @test task_fails(CancellationError) do
            task = with_cancel(() -> sleep(1))
            cancel(task)
            wait(task)
        end

        @test task_fails("foobar") do
            wait(with_cancel() do
                throw("foobar")
            end)
        end

        @testset "late cancel" begin
            let task = with_cancel(() -> sleep(1))
                @assert cancel(task) === task
                @test task_fails(CancellationError) do; wait(task); end
                @test cancel(task) === task
                @test task_fails(CancellationError) do; wait(task); end
            end

            let task = with_cancel(() -> 42)
                @assert fetch(task) == 42
                @test cancel(task) === task
                @test fetch(task) == 42
            end
        end
    end

    @testset "Timeouts" begin
        @test fetch(with_timeout(() -> 42, 1)) == 42

        @test task_fails(TimeoutError) do
            wait(with_timeout(.5) do; sleep(1) end)
        end

        @test task_fails("foobar") do
            wait(with_timeout(1) do
                throw("foobar")
            end)
        end

        @testset "late cancel" begin
            let task = with_timeout(() -> sleep(1), .5)
                @assert task_fails(TimeoutError) do; wait(task); end
                @test cancel(task) === task
                @test task_fails(TimeoutError) do; wait(task); end
            end

            let task = with_timeout(() -> 42, .5)
                @assert fetch(task) == 42
                @test cancel(task) === task
                @test fetch(task) == 42
            end
        end
    end
end

end # module TestCancellableTasks

max_stars_repo_name: Kirusifix/ExtraFun

max_issues_repo_name: Kirusifix/ExtraFun

max_forks_repo_name: Kirusifix/ExtraFun


text:
"""
DuckDB value
"""
mutable struct Value
    handle::duckdb_value

    function Value(handle::duckdb_value)
        result = new(handle)
        finalizer(_destroy_value, result)
        return result
    end
end

function _destroy_value(val::Value)
    if val.handle != C_NULL
        duckdb_destroy_value(val.handle)
    end
    val.handle = C_NULL
    return
end

getvalue(val::Value, ::Type{T}) where {T <: Int64} = duckdb_get_int64(val.handle)
function getvalue(val::Value, ::Type{T}) where {T <: String}
    ptr = duckdb_get_varchar(val.handle)
    result = unsafe_string(ptr)
    duckdb_free(ptr)
    return result
end
function getvalue(val::Value, ::Type{T}) where {T}
    throw(NotImplementedException("Unsupported type for getvalue"))
end

create_value(val::T) where {T <: Int64} = Value(duckdb_create_int64(val))
create_value(val::T) where {T <: AbstractString} = Value(duckdb_create_varchar_length(val, length(val)))
function create_value(val::T) where {T}
    throw(NotImplementedException("Unsupported type for getvalue"))
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
using Pkg;
current_env = Base.load_path()[1]
Pkg.activate(@__DIR__)

using RadialBasisFunctionModels
using Documenter

include("make_literate.jl")

project_path = joinpath(@__DIR__, "..", "Project.toml")

test_path = joinpath(@__DIR__, "..", "test", "Project.toml")

if !(project_path ∈ Base.load_path())
    push!(LOAD_PATH, project_path)
end
if !(test_path ∈ Base.load_path())
    push!(LOAD_PATH, test_path)
end
#%%

DocMeta.setdocmeta!(RadialBasisFunctionModels, :DocTestSetup, :(using RadialBasisFunctionModels); recursive=true)

makedocs(;
    modules=[RadialBasisFunctionModels],
    authors="Manuel Berkemeier",
    repo="https://github.com/manuelbb-upb/RadialBasisFunctionModels.jl/blob/{commit}{path}#{line}",
    sitename="RadialBasisFunctionModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://manuelbb-upb.github.io/RadialBasisFunctionModels.jl",
        assets=String[],
        mathengine = Documenter.MathJax2()
    ),
    pages=[
        "Home" => "index.md",
        "Main Module" => "RadialBasisFunctionModels.md"
    ],
)

deploydocs(;
    repo="github.com/manuelbb-upb/RadialBasisFunctionModels.jl",
)

Pkg.activate(current_env)


text:
#https://arxiv.org/pdf/1901.05824.pdf

struct WI <: Algorithm
end

@with_kw struct WII <: Algorithm
    tol::Float64 = Defaults.tol
    maxiter::Int = Defaults.maxiter
end

function make_time_mpo(ham::MPOHamiltonian{S,T},dt,alg::WI) where {S,T}
    δ = dt*(-1im);
    W1 = PeriodicArray{T,3}(undef,ham.period,ham.odim-1,ham.odim-1);
    for loc in 1:size(ham,1)
        W1[loc,1,1] = ham[loc,1,1] + δ * ham[loc,1,ham.odim];

        for j in 2:size(ham,2)-1,
            k in 2:size(ham,3)-1
            W1[loc,j,k] = ham[loc,j,k]

            W1[loc,1,k] = ham[loc,1,k]*sqrt(δ)
            W1[loc,j,1] = ham[loc,j,ham.odim]*sqrt(δ)
        end
    end

    _makedense(W1)
end

function make_time_mpo(ham::MPOHamiltonian{S,T},dt,alg::WII) where {S,T}
    WA = PeriodicArray{T,3}(undef,ham.period,ham.odim-2,ham.odim-2);
    WB = PeriodicArray{T,2}(undef,ham.period,ham.odim-2);
    WC = PeriodicArray{T,2}(undef,ham.period,ham.odim-2);
    WD = PeriodicArray{T,1}(undef,ham.period);

    δ = dt*(-1im);

    for i in 1:ham.period,
        j in 2:ham.odim-1,
        k in 2:ham.odim-1

        init = [one(ham[i,1,ham.odim]),zero(ham[i,1,k]),zero(ham[i,j,ham.odim]),zero(ham[i,j,k])]

        (y,convhist) = exponentiate(1.0,RecursiveVec(init),Lanczos(tol = alg.tol,maxiter = alg.maxiter)) do x
            out = similar(x.vecs);

            @tensor out[1][-1 -2;-3 -4] := δ*x[1][-1,1,-3,-4]*ham[i,1,ham.odim][2,-2,2,1]

            @tensor out[2][-1 -2;-3 -4] := δ*x[2][-1,1,-3,-4]*ham[i,1,ham.odim][2,-2,2,1]
            @tensor out[2][-1 -2;-3 -4] += sqrt(δ)*x[1][1,2,1,-4]*ham[i,1,k][-1,-2,-3,2]

            @tensor out[3][-1 -2;-3 -4] := δ*x[3][-1,1,-3,-4]*ham[i,1,ham.odim][2,-2,2,1]
            @tensor out[3][-1 -2;-3 -4] += sqrt(δ)*x[1][1,2,1,-4]*ham[i,j,ham.odim][-1,-2,-3,2]

            @tensor out[4][-1 -2;-3 -4] := δ*x[4][-1,1,-3,-4]*ham[i,1,ham.odim][2,-2,2,1]
            @tensor out[4][-1 -2;-3 -4] += x[1][1,2,1,-4]*ham[i,j,k][-1,-2,-3,2]
            @tensor out[4][-1 -2;-3 -4] += sqrt(δ)*x[2][1,2,-3,-4]*ham[i,j,ham.odim][-1,-2,1,2]
            @tensor out[4][-1 -2;-3 -4] += sqrt(δ)*x[3][-1,1,2,-4]*ham[i,1,k][2,-2,-3,1]

            RecursiveVec(out)
        end
        convhist.converged == 0 && @warn "failed to exponentiate $(convhist.normres)"

        WA[i,j-1,k-1] = y[4];
        WB[i,j-1] = y[3];
        WC[i,k-1] = y[2];
        WD[i] = y[1];
    end

    W2 = PeriodicArray{T,3}(undef,ham.period,ham.odim-1,ham.odim-1);
    W2[:,2:end,2:end] = WA
    W2[:,2:end,1] = WB
    W2[:,1,2:end] = WC
    W2[:,1,1] = WD

    for i in 1:ham.odim-1,
        j in 1:ham.odim-1,
        k in 1:ham.odim-1

        @assert space(W2[1,i,j],3) == space(W2[1,j,k],1)'
    end

    _makedense(W2)
end

function _makedense(ham)
    @assert length(size(ham)) == 3

    domspaces = map(1:size(ham,1)) do loc
        [space(h,1) for h in ham[loc,:,1]]
    end

    embeds = PeriodicArray(_embedders.(domspaces));

    data = PeriodicArray(map(1:size(ham,1)) do loc
        reduce(+,map(Iterators.product(1:size(ham,2),1:size(ham,3))) do (i,j)
            @tensor temp[-1 -2;-3 -4]:=embeds[loc][i][-1,1]*ham[loc,i,j][1,-2,2,-4]*conj(embeds[loc+1][j][-3,2])
        end)
    end)

    InfiniteMPO(data)
end


text:
"""
    solve_multigrad(problem, gradients, args...; kwargs...)

Solve `problem` for multiple magnetic field gradients.
"""
function solve_multigrad(problem, gradients, args...; kwargs...)
    itertimes = zeros(0)
    results = Vector{output_type(problem)}()
    @info "Solving" typeof(problem)
    for grad ∈ gradients
        @info grad
        starttime = time()
        push!(results, solve(problem, grad, args...; kwargs...))
        push!(itertimes, time() - starttime)
    end
    results, itertimes
end


text:
using Gadfly, RDatasets, DataFrames

set_default_plot_size(6inch, 3inch)

iris = dataset("datasets", "iris")

plots = Plot[]
for species in ["setosa", "versicolor", "virginica"]
    iris_subset = iris[iris[:Species] .== species, :]
    push!(plots, plot(iris_subset, x=:SepalLength, y=:SepalWidth, Geom.point))
end
hstack(plots)


text:
# standardize and normalize
function stnorm(X, mX=0.0, sX=1.0)
    if mX == 0.0
        mX = mean(X,1)
        sX = std(X,1)
    end
    X = (X .- mX) ./ sX
    return X, mX, sX
end