text:
using LinearAlgebra, SparseArrays, InfiniteArrays, Infinities, FillArrays, LazyArrays, Statistics, DSP, BandedMatrices, LazyBandedMatrices, Test, Base64
import InfiniteArrays: InfUnitRange, InfStepRange, OneToInf, NotANumber, oneto, unitrange
import LazyArrays: CachedArray, MemoryLayout, LazyLayout, DiagonalLayout, LazyArrayStyle, colsupport, DualLayout
import BandedMatrices: _BandedMatrix, BandedColumns
import Base.Broadcast: broadcasted, Broadcasted, instantiate

@testset "construction" begin
    @testset "Array constructor errors" begin
        @test_throws ArgumentError Array{Float64}(undef, ∞)
        @test_throws ArgumentError Array{Float64}(undef, ∞, ∞)
        @test_throws ArgumentError Array{Float64}(undef, 1, ∞)
        @test_throws ArgumentError Array{Float64}(undef, ∞, 1)

        @test_throws ArgumentError Vector{Float64}(undef, ∞)
        @test_throws ArgumentError Matrix{Float64}(undef, ∞, ∞)
        @test_throws ArgumentError Matrix{Float64}(undef, 1, ∞)
        @test_throws ArgumentError Matrix{Float64}(undef, ∞, 1)

        @test_throws ArgumentError Array{Float64}(undef, (∞,))
        @test_throws ArgumentError Array{Float64}(undef, (∞, ∞))
        @test_throws ArgumentError Array{Float64}(undef, (1, ∞))
        @test_throws ArgumentError Array{Float64}(undef, (∞, 1))

        @test_throws ArgumentError Vector{Float64}(undef, (∞,))
        @test_throws ArgumentError Matrix{Float64}(undef, (∞, ∞))
        @test_throws ArgumentError Matrix{Float64}(undef, (1, ∞))
        @test_throws ArgumentError Matrix{Float64}(undef, (∞, 1))

        @test Array{Float64}(undef, ()) isa Array{Float64,0}
        @test Array{Float64,0}(undef, ()) isa Array{Float64,0}
    end

    @testset "similar" begin
        a = 1:∞
        @test similar(a) isa CachedArray{Int}
        @test similar(a, Float64) isa CachedArray{Float64}
        @test similar(a, 5) isa Vector{Int}
        @test similar(a, (6,)) isa Vector{Int}
        @test similar(a, Float64, 5) isa Vector{Float64}
        @test similar(a, Float64, (6,)) isa Vector{Float64}
        @test similar(a, Float64, Base.OneTo(5)) isa Vector{Float64}
        @test similar(a, Float64, (Base.OneTo(5),)) isa Vector{Float64}
        @test similar(a, ∞) isa CachedArray{Int}
        @test similar(a, (∞,)) isa CachedArray{Int}
        @test similar(a, Float64, ∞) isa CachedArray{Float64}
        @test similar(a, Float64, (∞,)) isa CachedArray{Float64}
        @test similar(a, Float64, (∞,∞)) isa CachedArray{Float64}
        @test similar(a, Float64, oneto(∞)) isa CachedArray{Float64}
        @test similar(a, Float64, (oneto(∞),)) isa CachedArray{Float64}
        @test similar(a, Float64, (oneto(∞),oneto(∞))) isa CachedArray{Float64}

        @test similar([1,2,3],Float64,()) isa Array{Float64,0}

        @test similar(a, Float64, (2,∞)) isa CachedArray{Float64}
        @test similar(a, Float64, (∞,2)) isa CachedArray{Float64}

        @test similar(Array{Float64}, (∞,)) isa CachedArray{Float64}
        @test similar(Array{Float64}, (∞,∞)) isa CachedArray{Float64}
        @test similar(Array{Float64}, (2,∞)) isa CachedArray{Float64}
        @test similar(Array{Float64}, (∞,2)) isa CachedArray{Float64}

        @test similar(Array{Float64}, (oneto(∞),)) isa CachedArray{Float64}
        @test similar(Array{Float64}, (oneto(∞),oneto(∞))) isa CachedArray{Float64}
        @test similar(Array{Float64}, (oneto(∞),oneto(5))) isa CachedArray{Float64}
        @test similar(Array{Float64}, (oneto(5),oneto(∞))) isa CachedArray{Float64}
    end

    @testset "zeros/fill/ones" begin
        a = zeros(1,∞)
        @test length(a) === ℵ₀
        @test size(a) === (1,ℵ₀)
        @test a isa CachedArray{Float64}
        @test all(iszero,a[1,1:100])
        a[5] = 1
        @test a[1,1:100] == [zeros(4); 1; zeros(95)]

        a = fill(1,∞)
        @test length(a) === ℵ₀
        @test size(a) === (ℵ₀,)
        @test a isa CachedArray{Int}
        @test all(x -> x===1,a[1:100])
        a[5] = 2
        @test a[1:100] == [fill(1,4); 2; fill(1,95)]

        a = ones(∞)
        @test a isa CachedArray{Float64}
        a[5] = 2
        @test a[1:100] == [fill(1,4); 2; fill(1,95)]

        @test ones(5,∞)[:,1:10] == ones(5,10)
        @test ones(∞,5)[1:10,:] == ones(10,5)
        @test ones(∞,∞)[1:5,1:5] == ones(5,5)

        @test zeros(∞, 5)[1:10,:] == zeros(10,5)
        @test zeros(Int, ∞, 5)[1:10,:] == zeros(10,5)
        @test zeros(Int, 5, ∞)[:,1:10] == zeros(5, 10)
    end
end

@testset "ranges" begin
    @test size(10:1:∞) == (ℵ₀,)
    @testset "colon" begin
        @test @inferred(10:1:∞) === @inferred(range(10; step=1, length=ℵ₀))
        @inferred(1:0.2:∞)  === @inferred(range(1; step=0.2, length=ℵ₀))
        @inferred(1.0:0.2:∞)  === @inferred(range(1.0; step=0.2, length=ℵ₀))
        @inferred(1:∞) === @inferred(range(1; length=ℵ₀))
    end
    @test_throws ArgumentError 2:-.2:∞
    @test_throws ArgumentError 0.0:-∞
    @test_throws ArgumentError ∞:-1:1

    @test_throws ArgumentError (2:.2:-∞)
    @test_throws ArgumentError (2.:.2:-∞)
    @test_throws ArgumentError (1:-∞)

    @test ∞:1 ≡ 1:0

    @testset "indexing" begin
        L32 = @inferred(Int32(1):∞)
        L64 = @inferred(Int64(1):∞)
        @test @inferred(L32[1]) === Int32(1) && @inferred(L64[1]) === Int64(1)
        @test L32[2] == 2 && L64[2] == 2
        @test L32[3] == 3 && L64[3] == 3
        @test L32[4] == 4 && L64[4] == 4

        @test @inferred((1.0:∞)[1]) === 1.0
        @test @inferred((1.0f0:∞)[1]) === 1.0f0
        @test @inferred((Float16(1.0):∞)[1]) === Float16(1.0)

        @test @inferred((0.1:0.1:∞)[2]) === 0.2
        @test @inferred((0.1f0:0.1f0:∞)[2]) === 0.2f0
        @test @inferred((1:∞)[1:4]) === 1:4
        @test @inferred((1.0:∞)[1:4]) === 1.0:4
        @test (2:∞)[1:4] == 2:5
        @test (1:∞)[2:5] === 2:5
        @test (1:∞)[2:2:5] === 2:2:4
        @test (1:2:∞)[2:6] === 3:2:11
        @test (1:2:∞)[2:3:7] === 3:6:13

        @test isempty((1:∞)[5:4])
        @test_throws BoundsError (1:∞)[8:-1:-2]
    end
    @testset "length" begin
        @test length(.1:.1:∞) == ℵ₀
        @test length(1.1:1.1:∞) == ℵ₀
        @test length(1.1:1.3:∞) == ℵ₀
        @test length(1:1:∞) == ℵ₀
        @test length(1:.2:∞) == ℵ₀
        @test length(1.:.2:∞) == ℵ₀
        @test length(2:-.2:-∞) == ℵ₀
        @test length(2.:-.2:-∞) == ℵ₀
    end

    @testset "intersect" begin
        @test intersect(1:∞, 2:3) == 2:3
        @test intersect(2:3, 1:∞) == 2:3
        @test intersect(1:∞, 2:∞) == 2:∞
        @test intersect(2:∞, 1:∞) == 2:∞
        @test intersect(-3:∞, 2:8) == 2:8
        @test intersect(1:∞, -2:3:15) == 1:3:15
        @test intersect(1:∞, -2:3:∞) == 1:3:∞
        @test intersect(1:11, -2:2:∞) == intersect(-2:2:∞,1:11) == 2:2:10
        @test intersect(1:∞, -2:1:15) == 1:15
        @test intersect(1:∞, 15:-1:-2) == 1:15
        @test intersect(1:∞, 15:-4:-2) == 3:4:15
        @test intersect(-20:∞, -10:3:-2) == -10:3:-2
        @test isempty(intersect(-5:5, -6:13:∞))
        @test isempty(intersect(1:∞, 15:4:-2))

        @test @inferred(intersect(0:3:∞, 0:4:∞)) == intersect(0:4:∞, 0:3:∞) == 0:12:∞

        @test intersect(24:-3:0, 0:4:∞) == 0:12:24
        @test_throws ArgumentError intersect(1:6:∞, 0:4:∞) # supporting empty would break type inferrence

        @test intersect(1:∞,3) == 3:3
        @test intersect(1:∞, 2:∞, UnitRange(3,7), UnitRange(4,6)) == UnitRange(4,6)

        @test intersect(1:∞, 2) === intersect(2, 1:∞) === 2:2
        @test_skip intersect(1.0:∞, 2) == intersect(2, 1.0:∞) == [2.0] # gives infinite loop
    end

    @testset "sort/sort!/partialsort" begin
        @test sort(1:∞) ≡ sort!(1:∞) ≡ 1:∞
        @test sort(2:2:∞) ≡ sort!(2:2:∞) ≡ 2:2:∞
        @test_throws ArgumentError sort(2:-2:-∞)
        @test_throws ArgumentError sort!(2:-2:-∞)
    end
    @testset "in" begin
        @test 0 in UInt(0):100:∞
        @test 0 in 0:-100:-∞
        @test ∞ ∉ UInt(0):100:∞

        @test !(3.5 in 1:∞)
        @test (3 in 1:∞)
        @test (3 in 5:-1:-∞)

        let r = 0.0:0.01:∞
            @test (r[30] in r)
        end
        let r = (-4*Int64(maxintfloat(Int === Int32 ? Float32 : Float64))):∞
            @test (3 in r)
            @test (3.0 in r)
        end
    end
    @testset "in() works across types, including non-numeric types (#21728)" begin
        @test 1//1 in 1:∞
        @test 1//1 in 1.0:∞
        @test !(5//1 in 6:∞)
        @test !(5//1 in 6.0:∞)
        @test Complex(1, 0) in 1:∞
        @test Complex(1, 0) in 1.0:∞
        @test Complex(1.0, 0.0) in 1:∞
        @test Complex(1.0, 0.0) in 1.0:∞
        @test_skip !(Complex(1, 1) in 1:∞)  # this is an infinite-loop at the moment
        @test_skip !(Complex(1, 1) in 1.0:∞)
        @test_skip !(Complex(1.0, 1.0) in 1:∞)
        @test_skip !(Complex(1.0, 1.0) in 1.0:∞)
        @test !(π in 1:∞)
    end

    @testset "indexing range with empty range (#4309)" begin
        @test (3:∞)[5:4] == 7:6
        @test (0:2:∞)[7:6] == 12:2:10
    end

    @testset "indexing with negative ranges (#8351)" begin
        for a=[3:∞, 0:2:∞], b=[0:1, 2:-1:0]
            @test_throws BoundsError a[b]
        end
    end

    @testset "sums of ranges" begin
        @test sum(1:∞) ≡ mean(1:∞) ≡ median(1:∞) ≡ ℵ₀
        @test sum(0:∞) ≡ mean(1:∞) ≡ median(1:∞) ≡ ℵ₀
        @test sum(0:2:∞) ≡ mean(0:2:∞) ≡ median(0:2:∞) ≡ RealInfinity(∞)
        @test sum(0:-2:-∞) ≡ mean(0:-2:-∞) ≡ median(0:-2:-∞) ≡ -∞
    end

    @testset "broadcasted operations with scalars" begin
        @test Base.BroadcastStyle(typeof(1:∞)) isa LazyArrayStyle{1}
        @test Base.BroadcastStyle(typeof(Base.Slice(1:∞))) isa LazyArrayStyle{1}

        @test broadcast(-, 1:∞, 2) ≡ -1:∞
        @test broadcast(-, 1:∞, 0.25) ≡ 1-0.25:∞
        @test broadcast(+, 1:∞, 2) ≡ 3:∞
        @test broadcast(+, 1:∞, 0.25) ≡ 1+0.25:∞
        @test broadcast(+, 1:2:∞, 1) ≡ 2:2:∞
        @test broadcast(+, 1:2:∞, 0.3) ≡ 1+0.3:2:∞
        @test broadcast(-, 1:2:∞, 1) ≡ 0:2:∞
        @test broadcast(-, 1:2:∞, 0.3) ≡ 1-0.3:2:∞
        @test broadcast(-, 2, 1:∞) ≡ 1:-1:-∞
        @test exp.((1:∞)') ≡ broadcast(exp, (1:∞)') ≡ exp.(1:∞)'
        @test exp.(transpose(1:∞)) ≡ broadcast(exp, transpose(1:∞)) ≡ transpose(exp.(1:∞))
        @test 1 .+ (1:∞)' ≡ broadcast(+, 1, (1:∞)') ≡ (2:∞)'
        @test 1 .+ transpose(1:∞) ≡ broadcast(+, 1, transpose(1:∞)) ≡ transpose(2:∞)
        @test (1:∞)' .+ 1 ≡ broadcast(+, (1:∞)', 1) ≡ (2:∞)'
        @test transpose(1:∞) .+ 1 ≡ broadcast(+, transpose(1:∞), 1) ≡ transpose(2:∞)
    end

    @testset "near-equal ranges" begin
        @test 0.0:0.1:∞ != 0.0f0:0.1f0:∞
    end

    @testset "comparing InfiniteUnitRanges and OneToInf" begin
        @test 1:2:∞ == 1:2:∞ != 1:3:∞ != 2:3:∞ == 2:3:∞ != 2:∞
        @test 1:1:∞ == 1:∞ == 1:∞ == OneToInf() == OneToInf()
    end

    @testset "Base.OneTo (misleading) overrides" begin
        @test_skip Base.OneTo{BigInt}(∞) isa OneToInf{BigInt}
        @test oneto(∞) isa OneToInf{Int}
    end

    @testset "issue #6973" begin
        r1 = 1.0:0.1:∞
        r2 = 1.0f0:0.2f0:∞
        r3 = 1:2:∞
        @test r1 + r1 == 2*r1
        @test_broken r1 + r2 == 2.0:0.3:∞
        @test (2r1)-3r1 == -1:(2step(r1)-3step(r1)):-∞
        @test_broken (r1 + r2) - r2 == r1
        @test r1 + r3 == 2.0:2.1:∞
        @test r3 + r3 == 2 * r3
    end

    @testset "Preservation of high precision upon addition" begin
        let r = (-0.1:0.1:∞) + broadcast(+, -0.3:0.1:∞, 1e-12)
            @test_broken r[3] == 1e-12
        end
    end

    @testset "issue #8584" begin
        @test (0:1//2:∞)[1:2:3] == 0:1//1:1
    end

    @testset "issue #9962" begin
        @test eltype(0:1//3:∞) <: Rational
        @test (0:1//3:∞)[1] == 0
        @test (0:1//3:∞)[2] == 1//3
    end
    @testset "converting ranges (issue #10965)" begin
        @test promote(0:∞, UInt8(2):∞) === (0:∞, 2:∞)
        @test convert(InfUnitRange{Int}, 0:∞) === 0:∞
        @test convert(InfUnitRange{Int128}, 0:∞) === Int128(0):∞

        @test InfUnitRange{Int16}(1:∞) ≡ AbstractVector{Int16}(1:∞) ≡
                AbstractArray{Int16}(1:∞) ≡ Int16(1):∞

        @test OneToInf{Int16}(OneToInf()) ≡ AbstractVector{Int16}(OneToInf()) ≡
                AbstractArray{Int16}(OneToInf()) ≡ OneToInf{Int16}()

        @test promote(0:1:∞, UInt8(2):UInt8(1):∞) === (0:1:∞, 2:1:∞)
        @test convert(InfStepRange{Int,Int}, 0:1:∞) === 0:1:∞
        @test convert(InfStepRange{Int128,Int128}, 0:1:∞) === Int128(0):Int128(1):∞

        @test promote(0:1:∞, 2:∞) === (0:1:∞, 2:1:∞)
        @test convert(InfStepRange{Int128,Int128}, 0:∞) === Int128(0):Int128(1):∞
        @test convert(InfStepRange, 0:∞) === 0:1:∞
        @test convert(InfStepRange{Int128,Int128}, 0.:∞) === Int128(0):Int128(1):∞

        @test_broken promote(0f0:inv(3f0):∞, 0.:2.:∞) === (0:1/3:∞, 0.:2.:∞)

        @test promote(0:1/3:∞, 0:∞) === (0:1/3:∞, 0.:1.:∞)

        @test AbstractArray{Float64}(1:2:∞) ≡ AbstractVector{Float64}(1:2:∞) ≡
                convert(AbstractVector{Float64}, 1:2:∞) ≡ convert(AbstractArray{Float64}, 1:2:∞)

        @test unitrange(oneto(∞)) ≡ InfUnitRange(oneto(∞)) ≡ InfUnitRange{Int}(oneto(∞)) ≡ InfUnitRange(1)

        @test oneto(RealInfinity()) ≡ oneto(ComplexInfinity()) ≡ OneToInf()
        @test_throws ArgumentError oneto(-ComplexInfinity())
        @test_throws ArgumentError oneto(-∞)
    end

    @testset "inf-range[inf-range]" begin
        @test (1:∞)[1:∞] == 1:∞
        @test (1:∞)[OneToInf()] == 1:∞
        @test (1:∞)[2:∞] == 2:∞
        @test_throws BoundsError (1:∞)[-1:∞]
        @test (1:-1:-∞)[1:∞] == 1:-1:-∞
    end

    @testset "OneToInf" begin
        r = OneToInf()
        @test !isempty(r)
        @test length(r) == ℵ₀
        @test size(r) == (ℵ₀,)
        @test step(r) == 1
        @test first(r) == 1
        @test last(r) == ∞
        @test minimum(r) == 1
        @test maximum(r) == ℵ₀
        @test r[2] == 2
        @test r[2:3] === 2:3
        @test_throws BoundsError r[0]
        @test broadcast(+, r, 1) === 2:∞
        @test 2*r === 2:2:∞
        @test r + r === 2:2:∞

        @test r - r === Zeros{Int}(∞)

        @test intersect(r, Base.OneTo(2)) == Base.OneTo(2)
        @test intersect(r, 0:5) == 1:5
        @test intersect(r, 2) === intersect(2, r) === 2:2

        @test Base.unsafe_indices(Base.Slice(r)) == (r,)

        @testset "iteration with zip + finite iterator" begin
            z = zip(OneToInf(), 1:100)
            @test axes(z) == (Base.OneTo(100),)
            @test size(z) == (100,)
        end
    end

    @testset "show" begin
        # NOTE: Interpolating Int to ensure it's displayed properly across 32- and 64-bit
        @test summary(1:∞) == "ℵ₀-element InfUnitRange{$Int} with indices OneToInf()"
        @test Base.inds2string(axes(1:∞)) == "OneToInf()"
    end

    @testset "end" begin
        @test oneto(∞)[end] ≡ oneto(∞)[∞] ≡ ℵ₀
        @test (1:∞)[end] ≡ (1:∞)[∞] ≡ ℵ₀
        @test (1:2:∞)[end] ≡ (1:2:∞)[∞] ≡ ℵ₀
        @test (1.0:2:∞)[end] ≡ (1.0:2:∞)[∞] ≡ ℵ₀
    end

    @testset "union" begin
        @test @inferred((1:∞) ∪ (3:∞)) ≡ @inferred((3:∞) ∪ (1:∞)) ≡ 1:∞
        @test @inferred((1:∞) ∪ (3:1:∞)) ≡ @inferred((3:1:∞) ∪ (1:∞)) ≡ 1:1:∞
        @test @inferred((2:2:∞) ∪ (4:2:∞)) ≡ 2:2:∞
        @test (2.0:1.5:∞) ∪ (3.5:1.5:∞) ≡ 2.0:1.5:∞
        @test (1:∞) ∪ (2:2:∞) ≡ 1:1:∞
        @test (6:4:∞) ∪ (2:2:∞) ≡ 2:2:∞
        @test_throws ArgumentError (3:∞) ∪ (2:2:∞)
        @test_throws ArgumentError (2:2:∞) ∪ (3:∞)
        @test_throws ArgumentError (2:3:∞) ∪ (2:2:∞)
    end

    @testset "adjoint indexing" begin
        a = (1:∞)'
        @test a[:,:] ≡ a
        @test a[1,:] ≡ 1:∞
        @test a[1,2:2:end] ≡ 2:2:∞
    end

    @testset "big" begin
        @test (big(1):∞)[5] isa BigInt
        @test range(big(1), ∞; step=1) == big(1):1:∞
        @test range(big(1.0), ∞; step=2.5) == range(big(1.0); step=big(2.5), length=ℵ₀) == range(big(1.0), ∞; step=big(2.5)) == big(1.0):2.5:∞
    end

    @testset "maximum/minimum" begin
        @test maximum(2:∞) ≡ ℵ₀
        @test minimum(2:∞) ≡ 2
    end

    @testset "getindex[∞]" begin
        @test_throws BoundsError (2:6)[∞]
        @test (2:∞)[∞] ≡ ℵ₀
        @test oneto(∞)[∞] ≡ ℵ₀
    end

    @testset "show" begin
        @test stringmime("text/plain", 2:∞) ≡ "2:∞"
        @test stringmime("text/plain", OneToInf{BigInt}()) ≡ "OneToInf{BigInt}()"
    end

    @testset "in" begin
        @test ∞ ∉ (1:∞)
    end
end

@testset "fill" begin
    @testset "fill sizes" begin
        for A in (Zeros(∞), Fill(1,∞), Ones(∞),
                    Zeros(5,∞), Ones(5,∞), Fill(1,5,∞),
                    Zeros(∞,5), Ones(∞,5), Fill(1,∞,5),
                    Zeros(∞,∞), Ones(∞,∞), Fill(1,∞,∞))
            @test length(A) ≡ ℵ₀
        end
        @test size(Zeros(∞,5)) ≡ (ℵ₀,5)
        @test size(Zeros(5,∞)) ≡ (5,ℵ₀)
    end

    @testset "Fill indexing" begin
        B = Ones(∞,∞)
        @test IndexStyle(B) == IndexCartesian()
        V = view(B,:,1)
        @test_broken size(V) == (ℵ₀,1)
        V = view(B,1,:)
        @test size(V) == (ℵ₀,)
        V = view(B,1:1,:)
        @test size(V) == (1,ℵ₀)
    end

    @testset "Fill reindex" begin
        F = Fill(2.0,2,∞)
        @test reshape(F,∞) ≡ reshape(F,OneToInf()) ≡ reshape(F,(OneToInf(),)) ≡ reshape(F,Val(1)) ≡ Fill(2.0,∞)
    end

    @testset "adjtrans copy" begin
        @test copy((1:∞)') ≡ (1:∞)'
        @test copy(transpose(1:∞)) ≡ transpose(1:∞)
    end

    @testset "Fill slices" begin
        A = Fill(2,∞,∞)
        Z = Zeros(∞,∞)
        @test A[:,1] ≡ A[1,:] ≡ A[1:∞,1] ≡ Fill(2,∞)
        @test Z[:,1] ≡ Z[1,:] ≡ Z[1:∞,1] ≡ Zeros(∞)
    end

    @testset "maximum/minimum/sum" begin
        c = cache(Fill(2,∞));
        c[1] = 1;
        @test maximum(c) == maximum(Vcat([1], Fill(2,∞))) == 2
        c[1:3] = 1:3;
        @test maximum(c) == maximum(Vcat([1,2,3], Fill(2,∞))) == 3
        @test minimum(c) == minimum(Vcat([1,2,3], Fill(2,∞))) == 1
        @test sum([1; zeros(∞)]) ≡ 1.0
        @test sum([1; ones(∞)]) ≡ 1.0∞
    end
end

@testset "diagonal" begin
    D = Diagonal(1:∞)
    @test D[1:10,1:10] == Diagonal(1:10)
    @test_broken D^2 isa Diagonal
    @test D*D isa Diagonal
    @test MemoryLayout(typeof(D.diag)) == LazyLayout()
    @test MemoryLayout(typeof(D)) == DiagonalLayout{LazyLayout}()
    @test Base.BroadcastStyle(typeof(D)) == LazyArrayStyle{2}()
    @test Base.BroadcastStyle(typeof(permutedims(D.diag))) == LazyArrayStyle{2}()
    bc = broadcasted(*,Ones(∞,∞),permutedims(D.diag))
    @test bc isa Broadcasted{LazyArrayStyle{2}}
    @test instantiate(bc) isa Broadcasted{LazyArrayStyle{2}}
    @test copy(instantiate(bc)) isa BroadcastArray
    @test broadcast(*,Ones(∞,∞),permutedims(D.diag)) isa BroadcastArray
    @test Ones(∞,∞)*D isa BroadcastArray
    @test (Ones(∞,∞)*D)[1:10,1:10] == Ones(10,10)*D[1:10,1:10]
    @test @inferred(broadcast(*,Ones{Int}(∞),D)) ≡ @inferred(broadcast(*,D,Ones{Int}(∞))) ≡ D
    @test @inferred(broadcast(*,Ones(∞),D)) == @inferred(broadcast(*,D,Ones(∞))) == Diagonal(1.0:∞)
    @test @inferred(broadcast(*,Ones{Int}(∞)',D)) == @inferred(broadcast(*,D,Ones{Int}(∞)')) == D
    @test @inferred(broadcast(*,Ones(∞)',D)) == @inferred(broadcast(*,D,Ones(∞)')) == Diagonal(1.0:∞)
    @test @inferred(broadcast(*,Fill(2,∞)',D)) ≡ @inferred(broadcast(*,D,Fill(2,∞)')) ≡ 2D

    @test Eye{Int}(∞, ∞) isa Eye{Int}
    @test Eye{Int}(∞, 5) isa Eye{Int}
    @test Eye{Int}(5, ∞) isa Eye{Int}
    @test Eye(∞, ∞) isa Eye{Float64}
    @test Eye(∞, 5) isa Eye{Float64}
    @test Eye(5, ∞) isa Eye{Float64}

    @test Eye{Int}(∞) * D ≡ Eye{Int}(∞) * D ≡ D
    @test Eye(∞) * D == Eye(∞) * D == D
end

@testset "concat" begin
    @testset "concat indexing" begin
        A = Vcat(1:10, 1:∞)
        @test @inferred(length(A)) == ℵ₀
        @test @inferred(A[5]) == A[15] == 5
        @test A[end] == @inferred(A[∞]) == ∞
        @test_throws BoundsError Vcat(1:10)[∞]

        A = Vcat(Ones(1,∞), Zeros(2,∞))
        @test @inferred(size(A)) == (3,ℵ₀)
        @test @inferred(A[1,5]) == 1
        @test @inferred(A[3,5]) == 0
        @test_throws BoundsError A[4,1]

        A = Vcat(Ones{Int}(1,∞), Diagonal(1:∞))
        @test @inferred(size(A)) ≡ (ℵ₀,ℵ₀)
        @test @inferred(A[1,5]) ≡ 1
        @test @inferred(A[5,5]) ≡ 0
        @test @inferred(A[6,5]) ≡ 5
        @test_throws BoundsError A[-1,1]

        A = Vcat(Ones{Float64}(1,∞), Diagonal(1:∞))
        @test @inferred(size(A)) ≡ (ℵ₀,ℵ₀)
        @test @inferred(A[1,5]) ≡ 1.0
        @test @inferred(A[5,5]) ≡ 0.0
        @test @inferred(A[6,5]) ≡ 5.0
        @test_throws BoundsError A[-1,1]

        A = Vcat(1, Zeros(∞))
        @test @inferred(A[1]) ≡ 1.0
        @test @inferred(A[2]) ≡ 0.0

        A = Hcat(Ones(∞), Zeros(∞,2))
        @test @inferred(size(A)) == (ℵ₀,3)
        @test @inferred(A[5,1]) == 1
        @test @inferred(A[5,3]) == 0
        @test_throws BoundsError A[1,4]

        A = Hcat(Ones{Int}(∞), Diagonal(1:∞))
        @test @inferred(size(A)) ≡ (ℵ₀,ℵ₀)
        @test @inferred(A[5,1]) ≡ 1
        @test @inferred(A[5,5]) ≡ 0
        @test @inferred(A[5,6]) ≡ 5
        @test_throws BoundsError A[-1,1]

        A = Hcat(1, (1:∞)')
        @test A[1,:] isa Vcat{<:Any,1}
        @test A[1,:][1:10] == A[1,1:10]
    end

    # This should be generalized, but it at the moment
    # it is restricted to a single Number. Support smart
    # addition for any number of Number/SVector's would be better
    # allowibng for the tail to be variable lenth
    @testset "Vcat special case" begin
        @test Vcat(1,Zeros{Int}(∞)) + Vcat(3,Zeros{Int}(∞)) ≡
            Vcat(1,Zeros{Int}(∞)) .+ Vcat(3,Zeros{Int}(∞)) ≡
            Vcat(4,Zeros{Int}(∞))

            size(Vcat(1:∞)) ≡ (ℵ₀,)
    end

    @testset "Vcat infrange getindex" begin
        x = Vcat(1, Fill(2,∞))
        @test x[1:end] ≡ x[1:∞] ≡ x
        @test x[3:end] ≡ x[3:∞] ≡ Fill(2,∞)
    end

    @testset "maximum/minimum Vcat" begin
        x = Vcat(1:2, [1,1,1,1,1], 3, Fill(4,∞))
        @test maximum(x) == 4
        @test minimum(x) == 1
    end

    @testset "special vcat" begin
        @test [1; Zeros(∞)][1:10] == [1; zeros(9)]
        @test [[1,2,3]; Zeros(∞)][1:10] == [1;2;3;zeros(7)]
        @test [1; zeros(∞)] isa CachedArray
        @test [[1,2,3]; zeros(∞)] isa CachedArray

        @test [1; 2; zeros(Int,∞)] isa CachedArray
        @test [1; 2; 3; zeros(Int,∞)] isa CachedArray
        @test [[1,2]; 3; zeros(Int,∞)] isa CachedArray
        @test [2; [1,2]; 3; zeros(Int,∞)] isa CachedArray

        @test [[1,2]; [3,4]; Zeros(∞)] isa Vcat{<:Any,1,<:Tuple{Array,Zeros}}
        @test [[1,2]; [3,4]; [5,6]; Zeros(∞)] isa Vcat{<:Any,1,<:Tuple{Array,Zeros}}

        @test [randn(2,2); Zeros(∞,2)] isa Vcat{<:Any,2,<:Tuple{Array,Zeros}}

        a = [[1,2,3]; zeros(Int,∞)]
        @test a[3:∞][1:5] == a[3:7]
        @test cache(1:∞)[2:∞][1:5] == 2:6
    end

    @testset "sparse print" begin
        A = Vcat(1, Zeros(∞))
        @test colsupport(A,1) == 1:1
        @test Base.replace_in_print_matrix(A, 2, 1, "0") == "⋅"
        if VERSION < v"1.6-"
			@test stringmime("text/plain", A; context=(:limit => true)) ==
					"vcat($Int, ℵ₀-element Zeros{Float64,1,Tuple{OneToInf{$Int}}} with indices OneToInf()) with indices OneToInf():\n 1.0\n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n ⋮"
		else
			@test stringmime("text/plain", A; context=(:limit => true)) ==
					"vcat($Int, ℵ₀-element Zeros{Float64, 1, Tuple{OneToInf{$Int}}} with indices OneToInf()) with indices OneToInf():\n 1.0\n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n  ⋅ \n ⋮"
		end
        A = Vcat(Ones{Int}(1,∞), Diagonal(1:∞))
        @test Base.replace_in_print_matrix(A, 2, 2, "0") == "⋅"
    end

    @testset "copymutable" begin
        @test Base.copymutable(Vcat(1., Zeros(∞))) isa CachedArray
        @test Base.copymutable(Vcat([1.], Zeros(∞))) isa CachedArray
        @test Base.copymutable(Vcat([1.,2.], zeros(∞))) isa CachedArray
        @test Base.copymutable(Vcat(1.,2., zeros(∞))) isa CachedArray
    end

    @testset "infinite indexing" begin
        a = Vcat(1, 1:∞)
        @test a[:] isa Vcat
        @test a[3:∞] ≡ 2:∞
        @test a[3:2:∞] isa Vcat

        A = Vcat(Ones(1,∞), Fill(2,1,∞))
        @test A[:,:] == A
        @test A[:,2:∞] isa Vcat

        A = Vcat(Ones(5,5), Fill(2,∞,5))
        @test A[:,:] == A
        @test A[2:∞,:] isa Vcat

        A = Vcat(Ones(1,∞), Fill(2,∞,∞))
        @test A[:,:] == A
        @test A[2:∞,2:∞] isa Vcat
        @test A[2:∞,2:∞][1:10,1:10] == fill(2,10,10)
    end

    @testset "adjoint copy"  begin
        a = Vcat(1,(1:∞))'
        b = transpose(Vcat(1,(1:∞)))
        @test copy(a) ≡ a
        @test copy(b) ≡ b
    end
end

@testset "broadcasting" begin
    @testset "∞ BroadcastArray" begin
        A = 1:∞
        B = BroadcastArray(exp, A)
        @test length(B) == ℵ₀
        @test B[6] == exp(6)
        @test exp.(A) ≡ B
        @test B[2:∞] isa BroadcastArray
        B = Diagonal(1:∞) .+ 1
        @test B isa BroadcastArray{Int}
        @test B[1,5] ≡ 1
        @test B[6,6] == 6+1
        B = Diagonal(1:∞) - Ones{Int}(∞,∞) # lowers to broadcast
        @test B isa BroadcastArray{Int}
        @test B[1,5] ≡ -1
        @test B[6,6] == 6-1
    end

    @testset "Broadcast Fill Lowers" begin
        @test broadcast(+, Zeros{Int}(∞) , Fill(1,∞)) isa Fill
        @test broadcast(+, Zeros{Int}(∞) , Zeros(∞)) isa Zeros
        @test broadcast(*, Ones(∞), Ones(∞)) ≡ Ones(∞)
        @test broadcast(*, Ones{Int}(∞), 1:∞) ≡ broadcast(*, 1:∞, Ones{Int}(∞)) ≡ 1:∞
        @test broadcast(*, Fill(2,∞), 1:∞) ≡ broadcast(*, 1:∞, Fill(2,∞)) ≡ 2:2:∞
        @test broadcast(*, Fill([1,2],∞), 1:∞) isa BroadcastVector
        @test broadcast(*, Fill([1,2],∞), 1:∞)[1:3] == broadcast(*, 1:∞, Fill([1,2],∞))[1:3] == [[1,2],[2,4],[3,6]]

        @test broadcast(*, 1:∞, Ones(∞)') isa BroadcastArray
        @test broadcast(*, 1:∞, Fill(2,∞)') isa BroadcastArray
        @test broadcast(*, Diagonal(1:∞), Ones{Int}(∞)') ≡ broadcast(*, Ones{Int}(∞)', Diagonal(1:∞)) ≡ Diagonal(1:∞)
        @test broadcast(*, Diagonal(1:∞), Fill(2,∞)') ≡ broadcast(*, Fill(2,∞)', Diagonal(1:∞)) ≡ Diagonal(2:2:∞)
    end

    @testset "subview inf broadcast" begin
        b = BroadcastArray(exp, 1:∞)
        v = view(b, 3:∞) .+ 1
        @test v isa BroadcastArray
        @test b[3:10] .+ 1 == v[1:8]
    end

    @testset "views of matrices" begin
        D = Diagonal(1:∞)
        V = Vcat(Ones(2,∞), D)
        @test view(D,:,5) .+ 1 isa BroadcastVector
        @test view(D,5,:) .+ 1 isa BroadcastVector
        @test view(V,:,5) .+ 1 isa BroadcastVector
        @test view(V,5,:) .+ 1 isa BroadcastVector

        @test view(D,2:∞,2:∞) .+ 1 isa BroadcastMatrix
        @test view(V,2:∞,2:∞) .+ 1 isa BroadcastMatrix

        @test view(D,2:∞,[1,2,3]) .+ 1 isa BroadcastMatrix
        @test view(D,[1,2,3],2:∞) .+ 1 isa BroadcastMatrix
        @test view(V,2:∞,[1,2,3]) .+ 1 isa BroadcastMatrix
        @test view(V,[1,2,3],2:∞) .+ 1 isa BroadcastMatrix
    end

    @testset "inf broadcast views" begin
        a = BroadcastArray(cos, 1:∞)
        r = div.(1:∞, 2) .+ 1
        b = SubArray(a, (r,))
        @test b[1:6] == a[r[1:6]]
        @test_broken Base.BroadcastStyle(typeof(b)) isa LazyArrayStyle
        c = SubArray(a, (view(r,2:∞),))
        @test c[1:6] == a[r[2:7]]
        @test Base.BroadcastStyle(typeof(c)) isa LazyArrayStyle
    end
end

@testset "Cumsum and diff" begin
    @test cumsum(Ones(∞)) ≡ 1.0:1.0:∞
    @test cumsum(Fill(2,∞)) ≡ 2:2:∞
    @test cumsum(Ones{Int}(∞)) ≡ oneto(∞)
    @test cumsum(Ones{BigInt}(∞)) ≡ OneToInf{BigInt}()

    @test diff(oneto(∞)) ≡ Ones{Int}(∞)
    @test diff(1:∞) ≡ Fill(1,∞)
    @test diff(1:2:∞) ≡ Fill(2,∞)
    @test diff(1:2.0:∞) ≡ Fill(2.0,∞)

    x = Vcat([3,4], Ones{Int}(5), 3, Fill(2,∞))
    y = @inferred(cumsum(x))
    @test y isa Vcat
    @test y[1:12] == cumsum(x[1:12])

    @test diff(x[1:10]) == diff(x)[1:9]
    @test diff(y)[1:20] == x[2:21]

    @test cumsum(x).args[2] ≡ 8:12
    @test last(y.args) == sum(x[1:9]):2:∞

    for r in (3:4:∞, 2:∞, oneto(∞))
        c = cumsum(r)
        @test c isa InfiniteArrays.RangeCumsum
        @test c[Base.OneTo(20)] == c[1:20] == [c[k] for k=1:20] == cumsum(r[1:20])
        @test c[2:20] == [c[k] for k=2:20] == cumsum(r[1:20])[2:end]
        @test c == c
        @test c[Base.OneTo(20)] isa InfiniteArrays.RangeCumsum
        @test exp.(c)[1:20] == exp.(c[1:20])
    end

    @test cumsum(1:∞)[2:∞][1:5] == cumsum(1:6)[2:end]
end

@testset "Sub-array" begin
    @test Ones(∞)[3:∞] ≡ Ones(∞)
    @test Ones{Int}(∞)[4:6] ≡ Ones{Int}(3)
    @test (1:∞)[3:∞] ≡ 3:∞
end

@testset "conv" begin
    @test conv(1:∞, [2]) ≡ conv([2], 1:∞) ≡ 2:2:∞
    @test conv(1:2:∞, [2]) ≡ conv([2], 1:2:∞) ≡ 2:4:∞
    @test conv(1:∞, Ones(∞))[1:5] == conv(Ones(∞),1:∞)[1:5] == [1,3,6,10,15]
    @test conv(Ones(∞), Ones(∞)) ≡ 1.0:1.0:∞
    @test conv(Ones{Int}(∞), Ones{Int}(∞)) ≡ oneto(∞)
    @test conv(Ones{Bool}(∞), Ones{Bool}(∞)) ≡ oneto(∞)
    @test conv(Fill{Int}(2,∞), Fill{Int}(1,∞)) ≡ conv(Fill{Int}(2,∞), Ones{Int}(∞)) ≡
                conv(Ones{Int}(∞), Fill{Int}(2,∞)) ≡ 2:2:∞
    @test conv(Ones{Int}(∞), [1,5,8])[1:10] == conv([1,5,8], Ones{Int}(∞))[1:10] ==
                    conv(fill(1,100),[1,5,8])[1:10] == conv([1,5,8], fill(1,100))[1:10]
    @test conv(Ones{Int}(∞), 1:4)[1:10] == conv(1:4, Ones{Int}(∞))[1:10] ==
                    conv(fill(1,100),collect(1:4))[1:10] == conv(collect(1:4), fill(1,100))[1:10]
end

@testset "Taylor ODE" begin
    e₁ = Vcat(1, Zeros(∞))
    D = Hcat(Zeros(∞), Diagonal(1:∞))
    I_inf = Eye(∞)
    @test I_inf isa Eye{Float64,OneToInf{Int}}
    @test axes(I_inf) == (OneToInf{Int}(), OneToInf{Int}())
    @test eltype(I_inf) == Float64
    @test Base.BroadcastStyle(typeof(I_inf)) == LazyArrays.LazyArrayStyle{2}()
    L = Vcat(e₁', I_inf + D)
    @test L[1:3,1:3] == [1.0 0.0 0.0;
                         1.0 1.0 0.0;
                         0.0 1.0 2.0]
end

@testset "show" begin
    @test repr(Vcat(1:∞)) == "[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, …]"
    @test repr(Vcat(2,1:∞)) == "[2, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, …]"
end

@testset "*" begin
    A = Fill(1,3,∞)
    B = Diagonal(1:∞)
    C = Vcat([1,2,3], Zeros(∞))
    D = Vcat(Fill(1,3,∞), Zeros(∞,∞))

    AB = A*B
    @test AB isa BroadcastArray
    @test size(AB) == (3,ℵ₀)
    @test (AB)[1:3,1:10] == Fill(1,3,10)*Diagonal(1:10)

    @test A*B*C isa ApplyArray
    @test size(A*B*C) == (3,)
    @test (A*B*C)[1] == 14
    @test A*B*C == fill(14,3)
    @test_throws BoundsError (A*B*C)[1:10,1:10]
    @test A*B*D isa ApplyArray
    @test (A*B*D)[1:3,1:5] == fill(6.0,3,5)
end

@testset "MemoryLayout" begin
    @test MemoryLayout(OneToInf{Int}) == LazyLayout()
    @test MemoryLayout(0:∞) == LazyLayout()
    @test MemoryLayout((0:∞)') == DualLayout{LazyLayout}()
    A = _BandedMatrix((0:∞)', ℵ₀, -1, 1)
    @test MemoryLayout(A) == BandedColumns{LazyLayout}()
end

@testset "Banded" begin
    A = _BandedMatrix((0:∞)', ℵ₀, -1, 1)
    @test (Eye{Int}(∞) * A).data ≡ A.data
    @test 2.0A isa BandedMatrix
    @test (2.0A)[1:10,1:10] == 2.0A[1:10,1:10]
    @test 2.0\A isa BandedMatrix
    @test (2.0\A)[1:10,1:10] == 2.0\A[1:10,1:10]
    @test A/2 isa BandedMatrix
    @test (A/2)[1:10,1:10] == (A/2)[1:10,1:10]

    @test A * Eye(∞) isa BandedMatrix
    @test A / Eye(∞) isa BandedMatrix
    @test Eye(∞) * A isa BandedMatrix
    @test Eye(∞) \ A isa BandedMatrix
end

@testset "reshaped" begin
    @test InfiniteArrays.ReshapedArray(1:6,(2,3)) == [1 3 5; 2 4 6]
    @test InfiniteArrays.ReshapedArray(1:∞,(1,ℵ₀))[1,1:10] == 1:10
    @test reshape(1:∞,1,∞) === InfiniteArrays.ReshapedArray(1:∞,(1,ℵ₀))
    @test permutedims(1:∞) isa InfiniteArrays.ReshapedArray
    @test permutedims(1:∞)[1,1:10] == (1:10)
    a = reshape(Vcat(Fill(1,1,∞),Fill(2,2,∞)),∞)
    @test a[1:7] == [1, 2, 2, 1, 2, 2, 1]
    @test permutedims(permutedims(1:∞)) ≡ 1:∞
end

@testset "norm/dot" begin
    for p in (-Inf, 0, 0.1, 1, 2, 3, Inf)
        @test norm(Zeros(∞), p) == 0.0
        @test norm(Fill(5),p) ≈ norm(Array(Fill(5)),p) # tests tuple bug
        @test norm(Zeros{Float64}(),p) == 0.0 # tests tuple bug
    end
    @test norm([1; zeros(∞)]) ≡ 1.0
    @test dot([1; zeros(∞)], [1; zeros(∞)]) ≡ 1.0
    @test dot([1; zeros(∞)], 1:∞) ≡ 1.0
    @test dot(1:∞, [1; zeros(∞)]) ≡ 1.0
end

@testset "sub-Eye" begin
    @test bandwidths(view(Eye(∞),:,2:∞)) == (1,-1)
end

@testset "findfirst" begin
    @test findfirst(isequal(5), OneToInf()) == 5
    @test isnothing(findfirst(isequal(5.5), OneToInf()))
    @test isnothing(findfirst(isequal(-1), OneToInf()))
    @test findfirst(isequal(5), 2:∞) == 4
    @test isnothing(findfirst(isequal(5.5), 2:∞))
    @test isnothing(findfirst(isequal(-1), 2:∞))
    @test findfirst(isequal(4), 2:2:∞) == 2
    @test isnothing(findfirst(isequal(5), 2:2:∞))
    @test isnothing(findfirst(isequal(5.5), 2:2:∞))
    @test isnothing(findfirst(isequal(-1), 2:2:∞))

    @test searchsorted(Vcat(2,3:∞),10) == 9:9
    @test searchsortedfirst(Vcat(2,3:∞),10) == 9
    @test searchsortedlast(Vcat(2,3:∞),10) == 9
    @test searchsortedlast(Vcat(2,3:∞),0) == 0

    @test searchsorted(factorial.(big(1):∞), 6) == 3:3
    @test searchsortedfirst(factorial.(big(1):∞), 7) == 4
    @test searchsortedlast(factorial.(big(1):∞), 7) == 3
end

@testset "convert infrange" begin
    @test convert(AbstractArray{Float64}, 1:∞) ≡ convert(AbstractArray{Float64}, oneto(∞)) ≡
        convert(AbstractVector{Float64}, 1:∞) ≡ convert(AbstractVector{Float64}, oneto(∞)) ≡
        AbstractVector{Float64}(1:∞) ≡ AbstractVector{Float64}(oneto(∞)) ≡
        AbstractArray{Float64}(1:∞) ≡ AbstractArray{Float64}(oneto(∞)) ≡ InfUnitRange(1.0)

    @test convert(AbstractArray{Float64}, (1:∞)') ≡ convert(AbstractArray{Float64}, oneto(∞)') ≡
        convert(AbstractMatrix{Float64}, (1:∞)') ≡ convert(AbstractMatrix{Float64}, oneto(∞)') ≡
        AbstractMatrix{Float64}((1:∞)') ≡ AbstractMatrix{Float64}(oneto(∞)') ≡
        AbstractArray{Float64}((1:∞)') ≡ AbstractArray{Float64}(oneto(∞)') ≡
        InfUnitRange(1.0)'

    @test convert(AbstractArray{Float64}, transpose(1:∞)) ≡ convert(AbstractArray{Float64}, transpose(oneto(∞))) ≡
        convert(AbstractMatrix{Float64}, transpose(1:∞)) ≡ convert(AbstractMatrix{Float64}, transpose(oneto(∞))) ≡
        AbstractMatrix{Float64}(transpose(1:∞)) ≡ AbstractMatrix{Float64}(transpose(oneto(∞))) ≡
        AbstractArray{Float64}(transpose(1:∞)) ≡ AbstractArray{Float64}(transpose(oneto(∞))) ≡
        transpose(InfUnitRange(1.0))
end

@testset "cached indexing" begin
    @test cache(1:∞)[Fill(2,∞)][1:10] == fill(2,10)
end


text:
"""
* `mpars` : medium parameters for forward or inverse mode
"""
function updateA!(pa::Param, mpars; A=pa.A)
	nz=pa.nz; nx=pa.nx; dz=pa.dz; dx=pa.dx
	nznx=nz*nx
	ddz = inv(pa.dz*pa.dz); # for simplicity below
	ddx = inv(pa.dx*pa.dx);
	updateAcore!(A,mpars,ddx,ddz,nx,nz)
	boundary!(A,nx,nz)

	# update At
	for i in 1:nznx+1
		for j in 1:nznx
			pa.At[j,i]=A[i,j]
		end
	end
	return pa
end

"""
#--------------------  Build explicit operator matrix A=∇⋅(mpars(x,z)∇) for direct Poisson solver --------------------
# Some notes:
# We build the operator matrix explicitly like this to properly deal with the Neumann boundary conditions at all boundaries.
# The standard Kronecker product approach for building sparse operator matrices seemed to generate some errors at the boundaries.
# Medium properties calculated on staggered grid.
# Directions: z-direction row-wise, x-direction column-wise.

#=
# Allocate non-zero operator entries
k = 5*(nz - 2)*(nx - 2)+ # pentadiagonal operator entry allocation: 5, since using 5-point finite-difference stencil (star-shaped)
    2*(2*(nz - 2) + 2*(nx - 2)+  # bidiagonal non-corner 4-boundary allocation # for boundaries (2 points only)
    4); # bidiagonal corner boundary allocation
# The sparse() function is often a handy way to construct sparse matrices.
# It takes as its input a vector I of row indices, a vector J of column indices, and a vector V of nonzero values.

avoid using this subroutine, use `updateA!` instead
"""
function updateAcore!(A::AbstractArray,mpars::AbstractArray,ddx,ddz,nx,nz)
	(length(mpars) ≠ (nz*nx))  && error("error in mpars")
	#-------------- Modify operator A for the interior points (without boundary)---------------------------------------
	for j = 2 : nx - 1 # interior x (column) loop (for amount of internal nodes)
	    k = 2 + nz*(j - 1); # initial interior z for this x. To jump from j (new column), need nz*... to proceed to next column vector k
	    for i = 2 : nz - 1  # interior z (row) loop:
		A[k,k - nz] = .5*(mpars[i+(j-1)*nz] + mpars[i+(j-2)*nz])*ddx; # moving left (and placing actual staggered 2D grid medium value there)
		A[k,k - 1] = .5*(mpars[i+(j-1)*nz] + mpars[i-1+(j-1)*nz])*ddz; # moving up the grid
		A[k,k] = -(mpars[i+(j-1)*nz] + .5*(mpars[i+(j-2)*nz] + mpars[i+(j)*nz]))*ddx -
		(mpars[i+(j-1)*nz] + .5*(mpars[i-1+(j-1)*nz] + mpars[i+1+(j-1)*nz]))*ddz; # diagonal elements: ∇⋅\medpars(x,z)∇ calculated using 5 mpars points
		A[k,k + 1] = .5*(mpars[i+(j-1)*nz] + mpars[i+1+(j-1)*nz])*ddz; # moving down the grid
		A[k,k + nz] = .5*(mpars[i+(j-1)*nz]+ mpars[i+(j)*nz])*ddx; # moving right on grid (need nz values to go to next column)
		k = k + 1; # to next z
	    end # The foregoing A terms are operations by hand ...
	end
	return A

end

"""
transpose of dAdx
"""
function dAdmpars(dx,dz,nx,nz,T)

	ddz = inv(dz*dz); # for simplicity below
	ddx = inv(dx*dx);
	X=[spzeros(T,(nz*nx)+1,nz*nx) for i in 1:nz*nx]  # note the size, it has appended zeros

	nznx1=(nz*nx)+1

	for j = 2 : nx - 1
	    k = 2 + nz*(j - 1);
	    for i = 2 : nz - 1

		imp=i+(j-1)*nz; imp_0_m1=i+(j-2)*nz;
		imp_m1_0=i-1+(j-1)*nz; imp_0_p1=i+(j)*nz;
		imp_p1_0=i+1+(j-1)*nz;

		k1=k+(k-nz-1)*nznx1; k2=k+(k-1-1)*nznx1;
		k3=k+(k-1)*nznx1; k4=k+(k)*nznx1;
		k5=k+(k+nz-1)*nznx1;

		X[imp][k1] += .5*ddx; ; X[imp_0_m1][k1] += .5*ddx; ;
		X[imp][k2] += .5*ddz; ; X[imp_m1_0][k2] += .5*ddz; ;

		X[imp][k3] -= ddx  ; X[imp][k3] -= ddz; ;
		X[imp_0_m1][k3] -= .5*ddx;; X[imp_0_p1][k3] -= .5*ddx;;
		X[imp_m1_0][k3] -= .5*ddz;; X[imp_p1_0][k3] -= .5*ddz;;

		X[imp][k4] += .5*ddz; ; X[imp_p1_0][k4] += .5*ddz; ;
		X[imp][k5] += .5*ddx; ; X[imp_0_p1][k5] += .5*ddx; ;
		k = k + 1;
	    end
	end

	for j = 2 : nx - 1
	    k = 1 + nz*(j - 1);
	    for i = [1 nz]
		o = k - Int(sign(i - .5*nz));
		k1=k+(k-1)*nznx1
		k2=k+(o-1)*nznx1
		k3=o+(k-1)*nznx1
		for chu in 1:nz*nx
			X[chu][k1] = -X[chu][k3];
			X[chu][k2] = -X[chu][k1];
		end
		k = k + nz - 1;
	    end
	end

	for j = [1 nx]
	    k = 2 + nz*(j - 1);
	    for i = 2 : nz - 1
		o = k - Int(sign(j - .5*nx))*nz;
		k1=k+(k-1)*nznx1
		k2=k+(o-1)*nznx1
		k3=o+(k-1)*nznx1
		for chu in 1:nz*nx
			X[chu][k1] = -X[chu][k3];
			X[chu][k2] = -X[chu][k1];
		end
		k = k + 1;
	    end
	 end

	for j = [1 nx]
	    k = 1 + nz*(j - 1);
	    for i = [1 nz]
	       o = k - Int(sign(i - .5*nz))-
		Int(sign(j - .5*nx))*nz;
		k1=k+(k-1)*nznx1
		k2=k+(o-1)*nznx1
		k3=o+(o-1)*nznx1
		for chu in 1:nz*nx
			X[chu][k1] = X[chu][k3];
			X[chu][k2] = -X[chu][k1];
		end
		k = k + nz - 1; #
	    end
	 end

	 return X
end

function boundary!(A, nx, nz)

	#------------- Modify operator A for the top and bottom boundaries including Neumann conditions (without the corners)-------------
	for j = 2 : nx - 1 # interior x (column) loop:
	    k = 1 + nz*(j - 1);# initial boundary z for this x (k = i + m*(j - 1))
	    for i = [1 nz]  # % north & south boundary z (row) loop:
		o = k - Int(sign(i - .5*nz)); #% sign function determines correct direction to: inward from north or south boundary # to define direction(from boundary down or up)
		A[k,k] = -A[o,k]; #% build in symmetry ... --> trick to get symmetry, the next line is leading.
		A[k,o] = -A[k,k]; #% cancelation implies equation # Neumann: value of psi = same. The operator does not know value. Just knows has to subtract two values --> minus sign
		k = k + nz - 1; #% from north to south boundary
	    end
	end

	# Modify A for the east and west boundaries without corner points
	for j = [1 nx] # west & east boundary x (column) loop:
	    k = 2 + nz*(j - 1); #% initial interior z for this x (k = i + m*(j - 1))
	    for i = 2 : nz - 1 # interior z (row) loop:
		o = k - Int(sign(j - .5*nx))*nz; # inward from west or east boundary
		A[k,k] = -A[o,k]; #% build in symmetry ...
		A[k,o] = -A[k,k]; #% cancelation implies equation
		k = k + 1; #% to next z
	    end
	 end

	# Modify A for the corner points
	for j = [1 nx] #% west & east boundary x (column) loop:
	    k = 1 + nz*(j - 1); #% start on north boundary
	    for i = [1 nz] #% north & south boundary y (row) loop:
	       o = k - Int(sign(i - .5*nz))- #% inward from north or south boundary
		Int(sign(j - .5*nx))*nz; # inward from west or east boundary # two terms to move diagonal
		A[k,k] = A[o,o]; #% nonzero, alas asymmetric ...
		A[k,o] = -A[k,k]; #% ... cancelation implies equation
		k = k + nz - 1; #% from north to south boundary
	    end
	 end

	# append A with ones, to remove mean
	# To 'fix' the nullspace, and make the source consistent with the nullspace
	#A[nz*nx+1,:] .= 1.0
	A[nz*nx+1,:] .= 0.0
	# the average of four corners would go to zero
	A[nz*nx+1,1] = 1.0
	A[nz*nx+1,nz] = 1.0
	A[nz*nx+1,nz*(nx-1)+1] = 1.0
	A[nz*nx+1,nz*nx] = 1.0

	return nothing
end

ext: jl


text:
import PredictMD

PredictMD.import_all()
PredictMD.import_all(Main)


text:
function snake11rrr(r::Robot)
        j = 0
        while !isborder(r,Sud)
            if (isborder(r,Ost))
                j += movehole11111rrr(r,West)
            else
                j += movehole11111rrr(r,inverse(West))
            end
            if (isborder(r,Ost) || isborder(r,West))
                move!(r,Sud)
            end
        end
        if (isborder(r,Ost))
            movehole11111rrr(r,West)
        else
            movehole11111rrr(r,inverse(West))
        end
        return(j)
    end
function movehole11111rrr(r::Robot,side::HorizonSide)
    flag = 0
    n = 0
    while (!isborder(r,side))
        move!(r,side)
        if flag == 0
            if isborder(r,Sud)
                n += 1
                flag = 1
            end
        else
            if (!isborder(r,Sud))
                flag = 0
            end
        end
    end
    return(n)
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
"""
`Model` serves as a template for implementing a `compute(...)` function as required
for a `<:MeasurementModel`.
"""
struct Model <: MeasurementModel end
struct XLabel <: Label end
struct YLabel <: Label end
struct ALabel <: Label end
struct BLabel <: Label end

function NeXLUncertainties.compute(
    st::Model,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    # Build input variable labels
    xl, yl = XLabel(), YLabel()
    # Extract input variables
    x, y = inputs[xl], inputs[yl]
    # Compute the values
    a = fa(x, y)
    b = fb(x, y)
    al, bl = ALabel(), BLabel()
    vals = LabeledValues([al, bl], [a, b])
    jac = withJac ? zeros(Float64, length(vals), length(inputs)) : missing
    if withJac
        jac[indexin(vals, al), indexin(inputs, xl)] = δfaδx(x, y)
        jac[indexin(vals, al), indexin(inputs, yl)] = δfaδy(x, y)
        jac[indexin(vals, bl), indexin(inputs, xl)] = δfbδx(x, y)
        jac[indexin(vals, bl), indexin(inputs, yl)] = δfbδy(x, y)
    end
    return (vals, jac)
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()

##
using CSV
using DataFrames
using Statistics
using Distributions
using Pipe

##
using Conda
Conda.pip_interop(true)
Conda.pip("install", "ete3")

ENV["PYTHON"] = ""
Pkg.build("PyCall")
using PyCall

ete3 = pyimport("ete3")

##
data = CSV.read("../data/charMtx.csv", DataFrame)

##
worldGlotF = download("https://osf.io/jyvgt/download", "../data/world_fullGlot.tre")

glot = ete3.Tree(worldGlotF)

glot.prune(data.longname)

##
families = open("../data/glot3.txt") do file
    readlines(file)
end

try
    mkdir("mrbayes/")
catch e
end

##

function mbScript(fm, ngen, append)
    fmTaxa = filter(x -> x.glot_fam == fm, data).longname
    nex = """
#Nexus
\tBegin MrBayes;
\t\tset seed=6789580436154794230;
\t\tset swapseed = 614090213;
\t\texecute ../data/asjpNex/$fm.nex;
\t\tlset rates=gamma coding=all;
"""

fmGlot = glot.copy()
fmGlot.prune(fmTaxa)
constraints = []
if length(fmTaxa) > 5
    for nd in fmGlot.get_descendants()
        if !nd.is_leaf()
            push!(constraints, nd.get_leaf_names())
        end
    end
end

    if length(constraints) > 0
        for (i, cn) in enumerate(constraints)
            nex *= "\t\tconstraint c$i = " * join(cn, " ") * ";\n"
        end

        nex *= "\t\tprset topologypr = constraints("
        nex *= join(["c$i" for i in 1:length(constraints)], ",") * ");\n"
    end

    nex *= """
\t\tprset brlenspr = clock:uniform;
\t\tprset clockvarpr = igr;
\t\tprset treeagepr=Gamma(0.05, 0.005);
\t\tprset shapepr=Exponential(10);
\t\tset beagleprecision=double;
\t\tmcmcp Burninfrac=0.5 stoprule=no stopval=0.01;
\t\tmcmcp filename=../data/asjpNex/output/$fm;
\t\tmcmcp samplefreq=1000 printfreq=5000 append=$append;
\t\tmcmc ngen=$ngen nchains=4 nruns=2;
\t\tsump;
\t\tsumt;
\tend;
"""
    nex
end

##
for fm in families
    mbFile = "mrbayes/$(fm).mb.nex"
    nrun = 1000000
    open(mbFile, "w") do file
        write(file, mbScript(fm, nrun, "no"))
    end
    command = `mpirun -np 8 mb $mbFile`
    run(command)
    function converged()
        pstat = CSV.read(
            "../data/asjpNex/output/$fm.pstat",
            DataFrame,
            header = 2,
            datarow = 3,
            delim = "\t",
            ignorerepeated = true,
        )

        maxPSRF = maximum(pstat.PSRF)

        tstat = CSV.read(
            "../data/asjpNex/output/$fm.tstat",
            DataFrame,
            header = 2,
            datarow = 3,
            delim = "\t",
            ignorerepeated = true,
        )

        meanStdev = mean(tstat[:,4])

        vstat = CSV.read(
            "../data/asjpNex/output/$fm.vstat",
            DataFrame,
            header = 2,
            datarow = 3,
            delim = "\t",
            ignorerepeated = true,
            missingstring="NA",
        ) |> dropmissing

        maxPSRF = maximum([maxPSRF, maximum(vstat.PSRF)])
        maxPSRF <= 1.1 && meanStdev <= 0.01
    end
    while !converged()
        nrun += 1000000
        open(mbFile, "w") do file
            write(file, mbScript(fm, nrun, "yes"))
        end
        run(command)
    end
end

max_stars_repo_name: gerhardJaeger/phylogeneticTypology

max_issues_repo_name: gerhardJaeger/phylogeneticTypology

max_forks_repo_name: gerhardJaeger/phylogeneticTypology


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CraneMove <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CraneMove(; kwargs...)
        obj = new(meta(CraneMove), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CraneMove

const __meta_CraneMove = Ref{ProtoMeta}()

    function meta(::Type{CraneMove})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CraneMove)
            __meta_CraneMove[] = target = ProtoMeta(CraneMove)
            allflds = Pair{Symbol,Union{Type,String}}[:BlockId => Int32, :SourceId => Int32, :TargetId => Int32, :Sequence => Int32, :EmptyMove => Bool]
            meta(target, CraneMove, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CraneMove[]
    end
end

function Base.getproperty(obj::CraneMove, name::Symbol)
    if name === :BlockId
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :SourceId
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :TargetId
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :Sequence
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :EmptyMove
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

mutable struct CraneSchedule <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CraneSchedule(; kwargs...)
        obj = new(meta(CraneSchedule), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CraneSchedule

const __meta_CraneSchedule = Ref{ProtoMeta}()

function meta(::Type{CraneSchedule})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CraneSchedule)
            __meta_CraneSchedule[] = target = ProtoMeta(CraneSchedule)
            allflds = Pair{Symbol,Union{Type,String}}[:Moves => Base.Vector{CraneMove}, :SequenceNr => Int32]
            meta(target, CraneSchedule, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CraneSchedule[]
    end
end

function Base.getproperty(obj::CraneSchedule, name::Symbol)
    if name === :Moves
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CraneMove}
    elseif name === :SequenceNr
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct Performance <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Performance(; kwargs...)
        obj = new(meta(Performance), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Performance

const __meta_Performance = Ref{ProtoMeta}()

function meta(::Type{Performance})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Performance)
            __meta_Performance[] = target = ProtoMeta(Performance)
            allflds = Pair{Symbol,Union{Type,String}}[:CraneManipulations => Int32, :ServiceLevelMean => Float64, :LeadTimeMean => Float64, :DeliveredBlocks => Int32, :TotalBlocksOnTime => Int32, :BlockedArrivalTime => Float64, :TardinessMean => Float64, :BufferUtilizationMean => Float64, :CraneUtilizationMean => Float64, :HandoverUtilizationMean => Float64, :UpstreamUtilizationMean => Float64]
            meta(target, Performance, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Performance[]
        end
end

function Base.getproperty(obj::Performance, name::Symbol)
    if name === :CraneManipulations
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :ServiceLevelMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :LeadTimeMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :DeliveredBlocks
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :TotalBlocksOnTime
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :BlockedArrivalTime
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :TardinessMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :BufferUtilizationMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :CraneUtilizationMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :HandoverUtilizationMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :UpstreamUtilizationMean
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct TimeStamp <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function TimeStamp(; kwargs...)
        obj = new(meta(TimeStamp), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct TimeStamp

const __meta_TimeStamp = Ref{ProtoMeta}()

function meta(::Type{TimeStamp})
    ProtoBuf.metalock() do
        if !isassigned(__meta_TimeStamp)
            __meta_TimeStamp[] = target = ProtoMeta(TimeStamp)
            allflds = Pair{Symbol,Union{Type,String}}[:MilliSeconds => Int64]
            meta(target, TimeStamp, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_TimeStamp[]
    end
end

function Base.getproperty(obj::TimeStamp, name::Symbol)
    if name === :MilliSeconds
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct Block <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Block(; kwargs...)
        obj = new(meta(Block), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Block

const __meta_Block = Ref{ProtoMeta}()

function meta(::Type{Block})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Block)
            __meta_Block[] = target = ProtoMeta(Block)
            allflds = Pair{Symbol,Union{Type,String}}[:Id => Int32, :Release => TimeStamp, :Due => TimeStamp, :Ready => Bool]
            meta(target, Block, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Block[]
    end
end

function Base.getproperty(obj::Block, name::Symbol)
    if name === :Id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :Release
        return (obj.__protobuf_jl_internal_values[name])::TimeStamp
    elseif name === :Due
        return (obj.__protobuf_jl_internal_values[name])::TimeStamp
    elseif name === :Ready
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

    mutable struct Stack <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Stack(; kwargs...)
        obj = new(meta(Stack), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Stack

const __meta_Stack = Ref{ProtoMeta}()

function meta(::Type{Stack})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Stack)
            __meta_Stack[] = target = ProtoMeta(Stack)
            allflds = Pair{Symbol,Union{Type,String}}[:Id => Int32, :MaxHeight => Int32, :BottomToTop => Base.Vector{Block}]
            meta(target, Stack, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Stack[]
    end
end

function Base.getproperty(obj::Stack, name::Symbol)
    if name === :Id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :MaxHeight
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :BottomToTop
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Block}
    else
        getfield(obj, name)
    end
end

mutable struct Crane <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Crane(; kwargs...)
        obj = new(meta(Crane), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
        fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Crane

const __meta_Crane = Ref{ProtoMeta}()

function meta(::Type{Crane})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Crane)
            __meta_Crane[] = target = ProtoMeta(Crane)
            allflds = Pair{Symbol,Union{Type,String}}[:Id => Int32, :LocationId => Int32, :Load => Block, :Schedule => CraneSchedule, :GirderPosition => Float64, :HoistPosition => Float64]
            meta(target, Crane, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Crane[]
    end
end

function Base.getproperty(obj::Crane, name::Symbol)
    if name === :Id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :LocationId
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :Load
        return (obj.__protobuf_jl_internal_values[name])::Block
    elseif name === :Schedule
        return (obj.__protobuf_jl_internal_values[name])::CraneSchedule
    elseif name === :GirderPosition
        return (obj.__protobuf_jl_internal_values[name])::Float64
    elseif name === :HoistPosition
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct Handover <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Handover(; kwargs...)
        obj = new(meta(Handover), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Handover

const __meta_Handover = Ref{ProtoMeta}()
function meta(::Type{Handover})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Handover)
            __meta_Handover[] = target = ProtoMeta(Handover)
            allflds = Pair{Symbol,Union{Type,String}}[:Id => Int32, :Ready => Bool, :Block => Block]
            meta(target, Handover, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Handover[]
    end
end

function Base.getproperty(obj::Handover, name::Symbol)
    if name === :Id
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :Ready
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :Block
        return (obj.__protobuf_jl_internal_values[name])::Block
    else
        getfield(obj, name)
    end
end

mutable struct Uncertainties <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Uncertainties(; kwargs...)
        obj = new(meta(Uncertainties), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct Uncertainties

const __meta_Uncertainties = Ref{ProtoMeta}()

function meta(::Type{Uncertainties})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Uncertainties)
            __meta_Uncertainties[] = target = ProtoMeta(Uncertainties)
            pack = Symbol[:ArrivalIntervals,:CraneMoveTimes,:HandoverReadyIntervals]
            allflds = Pair{Symbol,Union{Type,String}}[:ArrivalIntervals => Base.Vector{Float64}, :CraneMoveTimes => Base.Vector{Float64}, :HandoverReadyIntervals => Base.Vector{Float64}]
            meta(target, Uncertainties, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Uncertainties[]
    end
end

function Base.getproperty(obj::Uncertainties, name::Symbol)
    if name === :ArrivalIntervals
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :CraneMoveTimes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
    elseif name === :HandoverReadyIntervals
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Float64}
        else
        getfield(obj, name)
    end
end

mutable struct World <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function World(; kwargs...)
        obj = new(meta(World), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct World

const __meta_World = Ref{ProtoMeta}()

function meta(::Type{World})
    ProtoBuf.metalock() do
        if !isassigned(__meta_World)
            __meta_World[] = target = ProtoMeta(World)
            allflds = Pair{Symbol,Union{Type,String}}[:Now => TimeStamp, :Production => Stack, :Buffers => Base.Vector{Stack}, :Handover => Handover, :Crane => Crane, :KPIs => Performance, :ObservationData => Uncertainties, :InvalidMoves => Base.Vector{CraneMove}]
            meta(target, World, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_World[]
    end
end

function Base.getproperty(obj::World, name::Symbol)
    if name === :Now
        return (obj.__protobuf_jl_internal_values[name])::TimeStamp
    elseif name === :Production
        return (obj.__protobuf_jl_internal_values[name])::Stack
    elseif name === :Buffers
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Stack}
    elseif name === :Handover
        return (obj.__protobuf_jl_internal_values[name])::Handover
    elseif name === :Crane
        return (obj.__protobuf_jl_internal_values[name])::Crane
    elseif name === :KPIs
        return (obj.__protobuf_jl_internal_values[name])::Performance
    elseif name === :ObservationData
        return (obj.__protobuf_jl_internal_values[name])::Uncertainties
    elseif name === :InvalidMoves
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CraneMove}
    else
        getfield(obj, name)
    end
end

export Block, Crane, CraneMove, CraneSchedule, Handover, Performance, Stack, TimeStamp, Uncertainties, World


text:
"""
    PLLParms

Struct for storing coefficients for the PLL filter.

Fields:

- `T::Float64`: integration time in seconds
- `damping::Float64`: PLL damping coefficient
- `B::Float64`: Pll bandwidth in Hz
- `wn::Float64`:
- `a0::Float64`:
- `a1::Float64`:
- `a2::Float64`:
- `b0::Float64`:
- `b1::Float64`:
- `b2::Float64`:
- `descriminator::String`: string containing descriminator equation expression
"""
struct PLLParms
    T::Float64
    damping::Float64
    B::Float64
    wn::Float64
    a0::Float64
    a1::Float64
    a2::Float64
    b0::Float64
    b1::Float64
    b2::Float64
    descriminator::String
end

"""
    DLLParms

Struct for storing DLL parameters.

Fields:

- `T::Float64`: integration time in seconds
- `B::Float64`: DLL bandwidth in Hz
- `d::Int64`: correlator sample spacing
	* this is the number of samples between descriminators
	* `NOT` the number of chips
- `descriminator::String`: string containing descriminator equation expression
"""
struct DLLParms
    T::Float64
    B::Float64
    d::Int64
    descriminator::String
end

"""
    TrackResults

A struct containing the parameters for tracking and its results.

Fields:

- `prn::Int64`: PRN number to be tracked
- `signaltype::String`: name of signal type
- `dll_parms::DLLParms`: DLL parameters
- `pll_parms::PLLParms`: PLL parameters
- `M::Int64`: the length of the data over the integration time or the number of
              tracking samples
- `T::Float64`: the tracking integration time in seconds
- `N::Int64`: the number of samples corresponding to a single integration time
- `data_file::String`: name of data file processed
- `data_type::String`: the data type (e.g. `sc4` or `sc8`)
- `data_nADC::Int64`: bit depth of data
- `data_start_idx::Int64`: starting index of data
- `data_t_length::Float64`: length of data in seconds loaded into `GNSSData`
                            struct
- `data_total_t_length::Float64`: total length of data file in seconds
- `data_sample_num::Int64`: number of data samples loaded into `GNSSData` struct
- `data_fs::Float64`: data sampling rate in Hz
- `data_fif::Float64`: data IF frequency in Hz
- `data_start_t::T1`: data start UTC time
- `data_site_lla::T2`: latitude, longitude, and height of receiver location,
                       where data was collected
- `data_init_n0::Float64`: index in data where code starts
- `data_init_code_chip::Float64`: fractional chip number at `t=0s`
- `data_init_phi::Float64`: initial phase in rads
- `data_init_fd::Float64`: initial Doppler frequency in Hz
- `t::Array{Float64,1}`: time vector in seconds with length `M`
- `code_err_meas::Array{Float64,1}`: vector containing the meassured code phase
                                     error
- `code_err_filt::Array{Float64,1}`: vector containing filtered code phase error
- `code_phase_meas::Array{Float64,1}`: vector containing the measured code phase
                                       estimate
- `code_phase_filt::Array{Float64,1}`: vector containing the filtered code phase
                                       estimate
- `n0::Array{Float64,1}`: vector containing the estimated chip per time
- `dphi_meas::Array{Float64,1}`: vector containing the measured phase error
                                 in rads
- `phi::Array{Float64,1}`: vector containing the filtered phase estimate in rads
- `delta_fd::Array{Float64,1}`: vector containing the Doppler frequency error in
                                Hz
- `fds::Array{Float64,1}`: vector containing the filtered Doppler frequency in
                           Hz
- `ZP::Array{Complex{Float64},1}`: vector containing the value of the prompt
                                   correlator at each time step
- `SNR::Array{Float64,1}`: vector containing the SNR of the correlation peak at
                           each time step in dB
- `data_bits::Array{Int64,1}`: extracted databits based off value of prompt
                               correlator
	* if `ZP[i] > 0`, `data_bits[i] = 1`
	* if `ZP[i] < 0`, `data_bits[i] = 0`
	* `NOTE:` there may be more samples corresponding to a single data bit
	          depending on the integration time
- `code_length::Int64`: the length of the primary code sequence

Kalman Filter Parameters:

- `Q::Array{Float64,2}`: process noise matrix `Q`
- `A::Array{Float64,2}`: state transition matrix `A`
- `C::Array{Float64,2}`: measurement matrix `C`
- `P::Array{Float64,2}`: state uncertainties per time step
- `x::Array{Float64,2}`: filted state estimates per time step
- `K::Array{Float64,2}`: KF gain per time step
- `R::Array{Float64,1}`: carrier phase measurement variance
- `P_full::Array{Float64,3}`: The full 2x2 or 3x3 state covariance matrix
"""
struct TrackResults{T1,T2}
    prn::Int64
    signaltype::String
    dll_parms::DLLParms
    pll_parms::PLLParms
    M::Int64
    T::Float64
    N::Int64
    data_file::String
    data_type::String
    data_nADC::Int64
    data_start_idx::Int64
    data_t_length::Float64
    data_total_t_length::Float64
    data_sample_num::Int64
    data_fs::Float64
    data_fif::Float64
    data_start_t::T1
    data_site_lla::T2
    data_init_n0::Float64
    data_init_code_chip::Float64
    data_init_phi::Float64
    data_init_fd::Float64
    t::Array{Float64,1}
    code_err_meas::Array{Float64,1}
    code_err_filt::Array{Float64,1}
    code_phase_meas::Array{Float64,1}
    code_phase_filt::Array{Float64,1}
    n0::Array{Float64,1}
    dphi_meas::Array{Float64,1}
    dphi_filt::Array{Float64,1}
    phi::Array{Float64,1}
    delta_fd::Array{Float64,1}
    fds::Array{Float64,1}
    ZP::Array{Complex{Float64},1}
    SNR::Array{Float64,1}
    data_bits::Array{Int64,1}
    code_length::Int64
    # Carrier phase Kalman fiter matrices and results
    Q::Array{Float64,2}
    A::Array{Float64,2}
    C::Array{Float64,2}
    P::Array{Float64,2}
    x::Array{Float64,2}
    dx::Array{Float64,2}
	K::Array{Float64,2}
    R::Array{Float64,1}
    P_full::Array{Float64,3}
end

"""
    definepll(T, B, damping)

Define `PLLParms` struct.

Required Arguments:

- `T`: integration time in seconds
- `B`: PLL bandwidth in Hz
- `ξ`: damping coefficient

Returns:

- `PLLParms` struct
"""
function definepll(T, B, ξ)
    ωₙ = 4*B/(2*ξ + 1/(2*ξ))
    a₀ = ωₙ^2*T^2 + 4*ξ*T*ωₙ + 4
    a₁ = 2*ωₙ^2*T^2 - 8
    a₂ = ωₙ^2*T^2 - 4*ξ*T*ωₙ + 4
    b₀ = ωₙ^2*T^2 + 4*ξ*T*ωₙ
    b₁ = 2*ωₙ^2*T^2
    b₂ = ωₙ^2*T^2 - 4*ξ*T*ωₙ
    descriminator = "ϕ = atan(Imag(ZP)/Real(ZP))"
    return PLLParms(T, ξ, B, ωₙ, a₀, a₁, a₂, b₀, b₁, b₂, descriminator)
end

"""
    definedll(T, B, d)

Define `DLLParms` struct.

Required Arguments:

- `T`: integration time in seconds
- `B`: PLL bandwidth in Hz
- `d`: correlator spacing in number of samples

Returns:

- `DLLParms` struct
"""
function definedll(T, B, d)
    descriminator = "Z4 = 1/d * (abs(ZE) - abs(ZL)) / (abs(ZE) + abs(ZL))"
    return DLLParms(T, B, d, descriminator)
end

"""
    measurephase(ZP)

Calculate the raw phase measurement.

Required Arguments:

- `ZP`: prompt correlator value

Returns:

-
"""
function measurephase(ZP)
    return atan((imag(ZP)/real(ZP)))
end

"""
    filtercarrierphase(pll_parms::PLLParms, ϕ_meas, ϕ_meas_1,
                       ϕ_meas_2, ϕ_filt_1, ϕ_filt_2)

Fiters the raw phase measurement using 2ⁿᵈ order PLL filter.

Required Arguments:

- `pll_parms::PLLParms`: struct containing PLL filter parameters
- `ϕ_meas`: phase measured at `t` step in rads
- `ϕ_meas_1`: phase measured at `t-T` step in rads
- `ϕ_meas_2`: phase measured at `t-2T` step in rads
- `ϕ_filt_1`: filtered phase from `t-T` step in rads
- `ϕ_filt_2`: filtered phase from `t-2T` step in rads

Returns:

- filtered phase from `t` step in rads
"""
function filtercarrierphase(pll_parms::PLLParms, ϕ_meas, ϕ_meas_1,
                            ϕ_meas_2, ϕ_filt_1, ϕ_filt_2)
    a₀ = pll_parms.a0
    a₁ = pll_parms.a1
    a₂ = pll_parms.a2
    b₀ = pll_parms.b0
    b₁ = pll_parms.b1
    b₂ = pll_parms.b2
    return (b₀*ϕ_meas + b₁*ϕ_meas_1 + b₂*ϕ_meas_2 - a₁*ϕ_filt_1 - a₂*ϕ_filt_2)/a₀
end

"""
    Z4(dll_parms::DLLParms, ZE, ZP, ZL)

Calculates the code phase error.

Required Arguments:

- `dll_parms::DLLParms`: struct containing DLL filter parameters
- `ZE`: early correlator value
- `ZP`: prompt correlator value
- `ZL`: late correlator value

Returns:

- unfiltered code phase error
"""
function Z4(dll_parms::DLLParms, ZE, ZP, ZL)
    return (abs(ZE) - abs(ZL)) / (abs(ZE) + abs(ZL))
end

"""
    filtercodephase(dll_parms::DLLParms, current_code_err, last_filt_code_err)

Returns the filtered code phase measusurement.

Required Arguments:

- `dll_parms::DLLParms`: struct containing DLL filter parameters
- `current_code_err`: current code phase error estimated from `Z4` descriminator
- `last_filt_code_err`: filtered code phase error from last time step

Returns:

- filtered code phase error
"""
function filtercodephase(dll_parms::DLLParms, current_code_err, last_filt_code_err)
    return last_filt_code_err+4*dll_parms.T*dll_parms.B*(current_code_err-last_filt_code_err)
end

"""
	shiftandcheck(i, offset, N)

Determines the sample location if shifting the current index `i` by some number,
`offset` in a vector with length `N`. Assumes that the new index can wrap.

Required Arguments:

- `i`: current index
- `offset`: sample offset amount
- `N`: size of vector

Returns:

- `j`: new index location
"""
function shiftandcheck(i, offset, N)
	j = i + offset
	if j > N
		j -= N
	end
	if j < 1
		j  += N
	end
	return j
end

"""
	getcorrelatoroutput!(ZP_array, data, replica, i, N, f_if, f_d,
						 fd_rate, ϕ, d, bin_width=1)

Calculate the early, prompt, and late correlator ouputs. Note that
replica already containts the prompt correlator. Be sure to set
the parameters to `replica` and run `generatereplica!(replica)` before
calling this method.

Required Arguments:

- `datafft`: allocated memory for the FFT result of the wiped off `data` vector
- `data`: vector containin I/Q samples to process
- `replica`: replica generated signal
- `i`: current `iᵗʰ` integration time `T`
- `N`: number of samples in `data` and `replica`
- `f_if`: intermediate frequency (IF) of `data` in Hz
- `f_d`: current Doppler frequency estimate in Hz
- `fd_rate`: current Doppler frequency rate estimate in Hz/s
- `ϕ`: current phase estimate in rads
- `d`: current correlator spacing in sample space
- `pfft`: result from `plan_fft(replica.data` which speeds up FFTs and IFFTs

Optional Arguments:

- `bin_width`: half width of frequency range to sum when estimating the SNR
               `(default = 1)`

Returns:

- `ze`: early correlator value
- `zp`: prompt correlator value
- `zl`: late correlator value
- `SNR`: signal-to-noise ratio of frequency domain correlation peak in dB
"""
function getcorrelatoroutput!(datafft, data, replica, i, N, f_if, f_d,
	                          fd_rate, ϕ, d, pfft, n0_index_start,
                              code_start_idx)
    datasegment = view(data.data, (i-1)*N+1:i*N)
    # datasegment = view(data.data, (i-1)*N+n0_index_start:i*N+n0_index_start-1)
    # Perform carrier wipeoff of intemediate and Doppler frequencies, Doppler
    # frequency rate, and carrier phase
    @inbounds for j in 1:N
        t = calc_t_at_i(j, data.start_t, data.f_s)
        datafft[j] = datasegment[j]*cis(-(2π*(f_if+f_d+0.5*fd_rate*t)*t+ϕ))
    end
    # In-place FFT of `datafft`
    pfft * datafft
    # In-place FFT of `replica.data`
    pfft * replica.data
    # Conjugate multiply `datafft` and `replica.data` and store result in
    # `replica.data`
    conjAmultB1D!(replica.data, datafft, N)
    # In-place IFFT of `replica.data`
    pfft \ replica.data
    # Correlation peak is in 1st index since replica was generated with same
    # code phase offset as the estimate
    ze_idx = N - (d - 1)
    zp_idx = 1
    zl_idx = 1 + d
    # zp_idx = floor(Int, code_start_idx)
    # ze_idx = zp_idx - d
    # if ze_idx < 1
    #     ze_idx = N + ze_idx
    # end
    # zl_idx = zp_idx + d
    # if zl_idx > N
    #     zl_idx = zl_idx%(N+1) + 1
    # end
    # Correlator outputs are the sum. Need to divide by the sample number, `N`.
    ze = replica.data[ze_idx]/N
    zp = replica.data[zp_idx]/N
    zl = replica.data[zl_idx]/N
    # Calculate SNR
    PS = abs2(replica.data[zp_idx])
    PN = (sum(abs2, replica.data) - PS) / (N - 1)
    SNR = calc_snr(PS, PN)
    return (ze, zp, zl, SNR)
end

"""
	calcA(T, state_num=2)

Calculate the state transition matrix, `A`, which is dependent on the
integration time, `T`. Used in the carrier tracking loop. If `state_num` is set
to `2`, then the KF will track the carrier phase and Doppler. If `state_num` is
set to `3`, then the KF will track the carrier phase, Doppler, and Doppler rate.

Required Arguments:

- `T`: integration time in seconds

Optional Arguments:

- `state_num`: number of states to track `(default = 2)`

Returns:

- either `2x2` or `3x3` state transition matrix
"""
function calcA(T, state_num=2)
	if state_num == 3
		return [1   T   T^2/2;
                0   1   T;
                0   0   1]
	elseif state_num == 2
		return [1   T;
                0   1]
	else
		error("Number of states specified must be either 2 or 3.")
	end
end

"""
	calcC(T, state_num=2)

Calculate the measurement matrix, `C`, using the integration time, `T`. Used in
the carrier tracking loop. If `state_num` is set to `2`, then the KF will track
the carrier phase and Doppler. If `state_num` is set to `3`, then the KF will
track the carrier phase, Doppler, and Doppler rate.

Required Arguments:

- `T`: integration time in seconds

Optional Arguments:

- `state_num`: number of states to track `(default = 2)`

Returns:

- either `1x2` or `1x3` measurement matrix
"""
function calcC(T, state_num=2)
	if state_num == 3
		return [1   T/2   T^2/6]
	elseif state_num == 2
		return [1   T/2]
	else
		error("Number of states specified must be either 2 or 3.")
	end
end

"""
	calcQ(T, h₀, h₋₂, qₐ, f_L, state_num=2)

Calculate the process noise covariance matrix, `Q`, which is dependent on the
integration time, `T`, and receiver oscillator h-parameters, `h₀` and `h₋₂`. Can
either specify `state_num` as 2 or 3. `f_L` is the carrier frequency. if
`state_num` is set to `3`, higher order effects from sources such as platform
dynamics, decribed with qₐ, are included. If `state_num` is set to `2`, then
they will not be included.

Required Arguments:

- `T`: integration time in seconds
- `h₀`: white frequency noise
- `h₋₂`: random walk frequency noise
- `qₐ`: platform dynamics in m²/s⁶
- `f_L`: carrier frequency in Hz

Optional Arguments:

- `state_num`: number of states to track `(default = 2)`

Returns:

- either `2x2` or `3x3` process noise matric `Q`
"""
function calcQ(T, h₀, h₋₂, qₐ, f_L, state_num=2)
	qϕ = h₀/2  # oscillator phase PSD
	qω = 2π^2*h₋₂  # oscillator frequency PSD
	if state_num == 3
        # First row
        Q11 = T*qϕ + T^3*qω/3 + T^5*qₐ/(20*c^2)
        Q12 = T^2*qω/2 + T^4*qₐ/(8*c^2)
        Q13 = T^3*qₐ/(6*c^2)
        # Second Row
        Q21 = T^2*qω/2 + T^4*qₐ/(8*c^2)
        Q22 = T*qω + T^3*qₐ/(3*c^2)
        Q23 = T^2*qₐ/(2*c^2)
        # Third row
        Q31 = T^3*qₐ/(6*c^2)
        Q32 = T^2*qₐ/(2*c^2)
        Q33 = T*qₐ/c^2
        return (2π*f_L)^2 .* [Q11 Q12 Q13;
                              Q21 Q22 Q23;
                              Q31 Q32 Q33]
	elseif state_num == 2
        # First row
        Q11 = T*qϕ + T^3*qω/3
        Q12 = T^2*qω/2
        # Second row
        Q21 = T^2*qω/2
        Q22 = T*qω
        return (2π*f_L)^2 .* [Q11 Q12;
                              Q21 Q22]
	else
		error("Number of states specified must be either 2 or 3.")
	end
end

"""
	trackprn(data::GNSSSignal, replica::ReplicaSignal, prn, ϕ_init,
			 fd_init, n0_idx_init, P₀, R; DLL_B=5, PLL_B=15, damping=1.4,
			 fd_rate=0., G=0.2, h₀=1e-21, h₋₂=2e-20, σω=10., qₐ=10.,
			 state_num=2, dynamickf=false, cov_mult=1., q_mult=1.,
			 channel="I")

Perform code and phase tracking on data in `data`.

`replica` decides the signal type. Can pass optional arguments that are minumum
amount to track a given PRN.

Required Arguments:

- `data::GNSSSignal`: either `GNSSData` or `ReplicaSignal` struct
- `replica::ReplicaSignal`: struct to use for replica signal generation
	* `replica.t_length` determines the integration time, `T`, used in tracking
- `prn`: PRN number to track
- `ϕ_init`: initial phase estimate in rads
- `fd_init`: initial Doppler frequency in Hz
- `n0_idx_init`: initial code start index location
- `P₀`: 3x3` diagonal matrix containing the initial uncertainties of the
        carrier phase, Doppler, and Doppler rate
- `R`: initial phase measurement uncertainty

Optional Arguments:

- `DLL_B`: bandwidth of the DLL filter in Hz `(default = 5)`
- `PLL_B`: bandwidth of the PLL filter in Hz `(default = 15)`
- `damping`: PLL filter damping coefficient `(default = 1.4)`
- `fd_rate`: Doppler rate in Hz `(default = 0)`
- `G`: `[DEPRICATED]` carrier phase filter gain
- `h₀`: `(default = 1e-21)`
- `h₋₂`: `(default = 2e-20)`
- `qₐ`: line of site platform dynamics in m²/s⁶ `(default = 1)`
- `state_num`: number of states to track `(default = 3)`
	* if set to `3`, carrier phase, Doppler, and Doppler rate are tracked
	* if set to `2`, only carrier phase and Doppler are tracked
- `dynamickf`: flag to specify if steady state KF gain is used or if KF can
               change from time step to time step `(default = true)`
- `cov_mult`: scalar to inflate the initial covariance matrix `P₀`
              `(default = 1)`
- `q_mult`: scalar to inflate the process noise matrix `Q` `(default = true)`
- `channel`: `[NOT USED]` `(default = "I")`

Returns:

- `TrackResults` struct
"""
function trackprn(data::GNSSSignal, replica::ReplicaSignal, prn, ϕ_init,
	              fd_init, n0_idx_init, P₀, R; DLL_B=0.1, PLL_B=15, damping=1.4,
				  fd_rate=0., G=0.2, h₀=1e-21, h₋₂=2e-20, qₐ=1.,
				  state_num=3, dynamickf=true, cov_mult=1., q_mult=1.,
				  R_mult=1, channel="I")
    # Assign signal specific parameters
    sig_freq = replica.signal_type.sig_freq
	if replica.signal_type.include_I
		chipping_rate = replica.signal_type.I_codes.chipping_rates[1]
    	ranging_code_length = replica.signal_type.I_codes.code_lengths[1]
	elseif replica.signal_type.include_Q
		chipping_rate = replica.signal_type.Q_codes.chipping_rates[1]
    	ranging_code_length = replica.signal_type.Q_codes.code_lengths[1]
	else
		error("Cannot track since there are no codes defined in signal type.")
	end
    # For codes with overlay codes, ensure that the code phase is modded
    # with the number of chips in a single unique sequence.
    #
    # Example: L5Q has 20ms overlay code. If T=20ms, ensure that the code phase
    # of the ranging code wraps around 20×ranging_code_length and not every 1ms,
    # since the Neuman-Hoffman sequence repeats every 20ms.
    code_length_ratio = floor(Int, replica.t_length/(ranging_code_length/chipping_rate))
    code_length = code_length_ratio*ranging_code_length
    # Place holder for fft of data wipeoff calculated at each T interval
    N = replica.sample_num
    datafft = Array{Complex{Float64}}(undef, N)
    # Pre-plan FFTs and IFFTs
    pfft = plan_fft!(replica.data)  # In-place FFT plan
    # pifft = plan_ifft!(replica.data)  # In-place IFFT plan
    # Compute the spacing between the ZE, ZP, and ZL correlators
	d = floor(Int, data.f_s/chipping_rate/2)
    # Initialize common variables and initial conditions
    T = replica.t_length
    t_length = data.t_length
    f_s = data.f_s
    f_if = data.f_if
    f_d = fd_init
    ϕ = ϕ_init
    f_code_d = chipping_rate*(1. + f_d/sig_freq)
    n0_init = calcinitcodephase(code_length,
                                f_code_d, 0.,
                                f_s, n0_idx_init)
    # n0_init = 0.0
    n0 = n0_init
    δn0 = 0.
    M = floor(Int, data.sample_num/replica.sample_num)
    # M = floor(Int, (data.sample_num - (n0_idx_init - 1))/replica.sample_num)
    t = Array(0:T:M*T-T)
    # Define DLL and PLL parameter structs
    dll_parms = definedll(T, DLL_B, d)
    pll_parms = definepll(T, PLL_B, damping)
    # Allocate array space for tracking results
    code_err_meas = Array{Float64}(undef, M)
    code_err_filt = Array{Float64}(undef, M)
    code_phase_meas = Array{Float64}(undef, M)
    code_phase_filt = Array{Float64}(undef, M)
    n0s = Array{Float64}(undef, M)
    dphi_measured = Array{Float64}(undef, M)
    dphi_filtered = Array{Float64}(undef, M)
    phi = Array{Float64}(undef, M)
    delta_fd = Array{Float64}(undef, M)
    fds = Array{Float64}(undef, M)
    ZP = Array{Complex{Float64}}(undef, M)
    SNR = Array{Float64}(undef, M)
    data_bits = Array{Int64}(undef, M)
    A = calcA(T, state_num)
    C = calcC(T, state_num)
    Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num) .* q_mult
    P = Array{Float64}(undef, state_num, M)
    x = Array{Float64}(undef, state_num, M)
    dx = Array{Float64}(undef, state_num, M)
	K = Array{Float64}(undef, state_num, M)
    P_full = Array{Float64}(undef, M, state_num, state_num)
    datafft = Array{Complex{Float64}}(undef, N)
	FFTW.set_num_threads(1)
	if state_num == 3
    	x⁻ᵢ = [ϕ_init; 2π*(f_if+fd_init); 2π*fd_rate]
		P⁻ᵢ = deepcopy(P₀)
        δx⁺ₖ = [0.; 0.; 0.]
	elseif state_num == 2
		x⁻ᵢ = [ϕ_init; 2π*(f_if+fd_init)]
		P⁻ᵢ = deepcopy(P₀)[1:2,1:2]
        δx⁺ₖ = [0.; 0.]
	else
		error("Number of states specified must be either 2 or 3.")
	end
    R = R_mult .* R
	P⁻ᵢ = cov_mult .* P⁻ᵢ
	Kᵢ = zeros(size(x⁻ᵢ))
    if ~dynamickf
	    Kfixed = dkalman(A, C, Q, Diagonal(R))
    end
    # Perform code, carrier phase, and Doppler frequency tracking
    for i in 1:M
		if state_num == 3
			ϕ, ω, ωdot = x⁻ᵢ
		else
			ϕ, ω = x⁻ᵢ
			ωdot = 0.
		end
		f_d = ω/2π - f_if
		fd_rate = ωdot/2π
		f_code_d = chipping_rate*(1. + f_d/sig_freq)
		f_code_dd = chipping_rate*fd_rate/sig_freq
        # Calculate the current code start index
        t₀ = ((code_length-n0)%code_length)/f_code_d
        code_start_idx = t₀*f_s + 1
        # Set signal parameters
        definereplica!(replica;
                       prn=prn, f_d=f_d,
                       fd_rate=fd_rate, phi=0., f_if=0.,
                       code_start_idx=code_start_idx,
                       include_carrier=false)
        # definereplica!(replica;
        #                prn=prn, f_d=f_d,
        #                fd_rate=fd_rate, phi=0., f_if=0.,
        #                code_start_idx=1,
        #                include_carrier=false)
        # Generate prompt correlator
        generatereplica!(replica)
        # Calculate early, prompt, and late correlator outputs
        ze, zp, zl, snr = getcorrelatoroutput!(datafft, data, replica, i, N,
		                                       f_if, f_d, fd_rate, ϕ, d,
                                               pfft, n0_idx_init,
                                               code_start_idx)
        # Estimate code phase error
        n0_err = Z4(dll_parms, ze, zp, zl)  # chips
        # Measure carrier phase and Doopler frequency errors
		# NOTE: `dϕ_meas` is considered the measurement error, `δy`
        dϕ_meas = measurephase(zp)  # rad
		# Estimate Kalman gain
		if dynamickf
            # Kᵢ = (P⁻ᵢ*C')/(C*P⁻ᵢ*C' + R)
            Kᵢ = (P⁻ᵢ*C')*pinv(C*P⁻ᵢ*C' + R)
		else
			Kᵢ = Kfixed
		end
        # Correct state uncertaninty
		# P⁺ᵢ = (I - Kᵢ*C)*P⁻ᵢ
        P⁺ᵢ = (I - Kᵢ*C)*P⁻ᵢ*(I - Kᵢ*C)' + Kᵢ*R*Kᵢ'
        # Correct current state estimate
        correction = Kᵢ.*dϕ_meas
        # correction = Kᵢ.*(dϕ_meas - (C*A*δx⁺ₖ)[1])
        # δx⁺ₖ[:] = correction
		x⁺ᵢ = x⁻ᵢ + correction
        if i > 1
            # Filter raw code phase error measurement
            # f_d_code = chipping_rate*(f_d - f_if)/sig_freq
            # f_dd_code = chipping_rate*fd_rate/sig_freq
            # n0_err_filtered = 4*T*DLL_B*(n0_err)
            n0_err_filtered = 4*T*DLL_B*(n0_err-code_err_filt[i-1])
            # n0_err_filtered = code_err_filt[i-1] +
            #                   4*T*DLL_B*(n0_err-code_err_filt[i-1])
            # n0_err_filtered = filtercodephase(dll_parms, n0_err, code_err_filt[i-1])
        else
            n0_err_filtered = 4*T*DLL_B*n0_err
            # n0_err_filtered = filtercodephase(dll_parms, n0_err, 0)
        end
        # Save to allocated arrays
        code_err_meas[i] = n0_err
        code_err_filt[i] = n0_err_filtered
        code_phase_meas[i] = (n0 + n0_err)
        code_phase_filt[i] = (n0 + n0_err_filtered)
        n0s[i] = n0
        dphi_measured[i] = dϕ_meas
        dphi_filtered[i] = correction[1]
        phi[i] = x⁺ᵢ[1]
        delta_fd[i] = correction[2]/2π
        fds[i] = x⁺ᵢ[2]/2π - f_if
        ZP[i] = zp
		SNR[i] = snr
        K[:,i] = Kᵢ
        x[:,i] = x⁺ᵢ
        dx[:,i] = correction
        P[:,i] = diag(P⁺ᵢ)
        P_full[i,:,:] = P⁺ᵢ
        if real(zp) > 0
            data_bits[i] = 1
        else
            data_bits[i] = 0
        end
        # Calculate main code chipping rate at next `i`
        f_code_d = chipping_rate*(1. + (x⁺ᵢ[2]/2π - f_if)/sig_freq)
        if state_num == 3
            f_code_dd = chipping_rate*(x⁺ᵢ[3]/2π)/sig_freq
        end
        # Update code phase with filtered code phase error and propagate to next `i`
		# This essetially makes the DLL rate aided.
        # δn0 = chipping_rate*(x⁺ᵢ[2]/2π - f_if)*T/sig_freq #+ 0.5*f_code_dd*T^2
        # n0 = (n0 + n0_err_filtered + chipping_rate*T + δn0)%code_length
        n0 += (n0_err_filtered + f_code_d*T + 0.5*f_code_dd*T^2)%code_length
        # n0 += (n0_err_filtered + f_code_d*T)%code_length
		# Propagate x⁺ᵢ to next time step
		x⁻ᵢ = A*x⁺ᵢ
        # Propogate state uncertaninty
		P⁻ᵢ = A*P⁺ᵢ*A' + Q
    end
	FFTW.set_num_threads(nthreads())
    # Return `TrackResults` struct
    # Specify what to store if `data` is a simulated signal struct
    # or `GNSSData` struct
    if typeof(data) == GNSSData
        file_name = data.file_name
        total_data_length = data.total_data_length
        data_type = data.data_type
        start_data_idx = data.start_data_idx
        site_loc_lla = data.site_loc_lla
        data_start_time = data.data_start_time
    else
        file_name = "N/A"
        total_data_length = data.t_length
        data_type = "N/A"
        start_data_idx = 1
        site_loc_lla = "N/A"
        data_start_time = "N/A"
    end
    return TrackResults(prn, replica.name,
                        dll_parms, pll_parms, M, T, N,
                        file_name, data_type, data.nADC,
                        start_data_idx, data.t_length,
                        total_data_length, data.sample_num, f_s,
                        f_if, data_start_time, site_loc_lla,
                        float(n0_idx_init), n0_init, ϕ_init, fd_init, t,
                        code_err_meas, code_err_filt, code_phase_meas,
                        code_phase_filt, n0s, dphi_measured, dphi_filtered,
                        phi, delta_fd, fds, ZP, SNR, data_bits, ranging_code_length,
                        Q, A, C, P, x, dx, K, R, P_full)
end

"""
	plotresults(results::TrackResults; saveto=missing,
				figsize=missing, doppler_curve=missing,
				doppler_t=missing, CN0=missing)

Plots the tracking results from the `trackprn` method.

Required Arguments:

- `results::TrackResults`: struct containing results from signal tracking

Optional Arguments:

- `saveto`: `String` specifying file name to save figure to `(default = missing)`
- `figsize`: `Tuple` of length 2 used to specify figure size in inches
	* format is `(height, width)`
	* `(default = missing)`
- `doppler_curve`: truth Doppler curve that will be plotted against Doppler
                   estimate `(default = missing)`
- `doppler_t`: time vector that must be included with `doppler_curve`
               `(default = missing)`
- `truth_SNR`: expected signal-to-noise `SNR` ratio `(default = missing)`

Plots tracking results figure.
"""
function plotresults(results::TrackResults; saveto=missing,
	                 figsize=missing, doppler_curve=missing,
					 doppler_t=missing, truth_SNR=missing)
	if ismissing(figsize)
		figure()
	else
		figure(figsize=figsize)
	end
    matplotlib.gridspec.GridSpec(3,2)
    # Plot code phase errors
    subplot2grid((3,2), (0,0), colspan=1, rowspan=1)
    plot(results.t, results.code_phase_meas.%results.code_length, "k.", label="Measured code phase")
    plot(results.t, results.code_phase_filt.%results.code_length, "b-", label="Filtered code phase")
    xlabel("Time (s)")
    ylabel("Code Phase (chips)")
    title("DLL Tracking")
    legend()
    # Plot filtered and measured phase errors
    subplot2grid((3,2), (0,1), colspan=1, rowspan=1)
    plot(results.t, results.dphi_meas.*(180/π), "k.", label="Measured Δϕ")
    plot(results.t, results.dphi_filt.*(180/π), "b-", label="Filtered Δϕ")
    xlabel("Time (s)")
    ylabel("ϕ (degrees)")
    title("PLL Tracking")
    legend()
	# Doppler frequency estimate
    subplot2grid((3,2), (1,0), colspan=1, rowspan=1)
	if ~ismissing(doppler_curve) && ~ismissing(doppler_t)
		plot(results.t, results.fds, "k.", label="Estimate")
		plot(doppler_t, doppler_curve, "b-", label="Truth")
		legend()
	else
    	plot(results.t, results.fds, "k.")
	end
    xlabel("Time (s)")
    ylabel("Doppler (Hz)")
    title("Doppler Frequency Estimate")
	# SNR estimate
	subplot2grid((3,2), (1,1), colspan=2, rowspan=1)
	if ~ismissing(truth_SNR)
		plot(results.t, results.SNR, "k.", label="Estimate")
		axhline(y=truth_SNR, color="b",
		        label="Truth")
		legend()
	else
    	plot(results.t, results.SNR, "k.")
	end
    xlabel("Time (s)")
    ylabel("SNR (dB)")
    # twinx()
    # plot(results.t, real(results.ZP), "k:")
    # ylabel("real(ZP)")
    title("SNR Estimate")
    # Plot ZP real and imaginary parts
    subplot2grid((3,2), (2,0), colspan=2, rowspan=1)
    plot(results.t, real(results.ZP), label="real(ZP)")
    plot(results.t, imag(results.ZP), label="imag(ZP)")
    xlabel("Time (s)")
    ylabel("ZP")
    title("Prompt Correlator Output")
    legend()
    suptitle("PRN $(results.prn)\nfd = $(Int(round(results.data_init_fd))) Hz\nn₀ = $(results.data_init_n0) samples")
    # subplots_adjust(hspace=0.4, wspace=0.4)
    subplots_adjust(hspace=0.4, wspace=0.2)
    if ~ismissing(saveto)
        savefig(saveto::String)
    end
end


text:
using EconPDEs, Distributions

Base.@kwdef mutable struct DiTellaModel
  # Utility Function
  γ::Float64 = 5.0
  ψ::Float64 = 1.5
  ρ::Float64 = 0.05
  τ::Float64 = 0.4

  # Technology
  A::Float64 = 200.0
  σ::Float64 = 0.03

  # MoralHazard
  ϕ::Float64 = 0.2

  # Idiosyncratic
  νbar::Float64 = 0.24
  κν::Float64 = 0.22
  σνbar::Float64 = -0.13
end

function (m::DiTellaModel)(state::NamedTuple, y::NamedTuple)
  (; γ, ψ, ρ, τ, A, σ, ϕ, νbar, κν, σνbar) = m
  (; x, ν) = state
  (; pA, pAx_up, pAx_down, pAν_up, pAν_down, pAxx, pAxν, pAνν, pB, pBx_up, pBx_down, pBν_up, pBν_down, pBxx, pBxν, pBνν, p, px_up, px_down, pν_up, pν_down, pxx, pxν, pνν) = y

  # drift and volatility of state variable ν
  g = p / (2 * A)
  i = A * g^2
  μν = κν * (νbar - ν)
  σν = σνbar * sqrt(ν)
  pAν = (μν >= 0) ? pAν_up : pAν_down
  pBν = (μν >= 0) ? pBν_up : pBν_down
  pν = (μν >= 0) ? pν_up : pν_down

  pAx, pBx, px = pAx_up, pBx_up, px_up
  iter = 0
  @label start
  σX = x * (1 - x) * (1 - γ) / (γ * (ψ - 1)) * (pAν / pA - pBν / pB) * σν / (1 - x * (1 - x) * (1 - γ) / (γ * (ψ - 1)) * (pAx / pA - pBx / pB))
  σpA = pAx / pA * σX + pAν / pA * σν
  σpB = pBx / pB * σX + pBν / pB * σν
  σp = px / p * σX + pν / p * σν
  κ = (σp + σ - (1 - γ) / (γ * (ψ - 1)) * (x * σpA + (1 - x) * σpB)) / (1 / γ)
  κν = γ * ϕ * ν / x
  σA = κ / γ + (1 - γ) / (γ * (ψ - 1)) * σpA
  νA = κν / γ
  σB = κ / γ + (1 - γ) / (γ * (ψ - 1)) * σpB
  # Interest rate r
  μX = x * (1 - x) * ((σA * κ + νA * κν - 1 / pA - τ) - (σB * κ -  1 / pB + τ * x / (1 - x)) - (σA - σB) * (σ + σp))
  if (iter == 0) & (μX <= 0)
    iter += 1
    pAx, pBx, px = pAx_down, pBx_down, px_down
    @goto start
  end

  μpA = pAx / pA * μX + pAν / pA * μν + 0.5 * pAxx / pA * σX^2 + 0.5 * pAνν / pA * σν^2 + pAxν / pA * σX * σν
  μpB = pBx / pB * μX + pBν / pB * μν + 0.5 * pBxx / pB * σX^2 + 0.5 * pBνν / pB * σν^2 + pBxν / pB * σX * σν
  μp = px / p * μX + pν / p * μν + 0.5 * pxx / p * σX^2 + 0.5 * pνν / p * σν^2 + pxν / p * σX * σν
  r = (1 - i) / p + g + μp + σ * σp - κ * (σ + σp) - γ / x * (ϕ * ν)^2

  # Market Pricing
  pAt = - pA * (1 / pA  + (ψ - 1) * τ / (1 - γ) * ((pA / pB)^((1 - γ) / (1 - ψ)) - 1) - ψ * ρ + (ψ - 1) * (r + κ * σA + κν * νA) + μpA - (ψ - 1) * γ / 2 * (σA^2 + νA^2) + (2 - ψ - γ) / (2 * (ψ - 1)) * σpA^2 + (1 - γ) * σpA * σA)
  pBt = - pB * (1 / pB - ψ * ρ + (ψ - 1) * (r + κ * σB) + μpB - (ψ - 1) * γ / 2 * σB^2 + (2 - ψ - γ) / (2 * (ψ - 1)) * σpB^2 + (1 - γ) * σpB * σB)
  # algebraic constraint
  pt = - p * ((1 - i) / p - x / pA - (1 - x) / pB)
  return (; pAt, pBt, pt)
end

m = DiTellaModel()
xn = 80
νn = 10
xs = range(0.01, 0.99, length = xn)
distribution = Gamma(2 * m.κν * m.νbar / m.σνbar^2, m.σνbar^2 / (2 * m.κν))
νs = range(quantile(distribution, 0.001), quantile(distribution, 0.999), length = νn)
stategrid = OrderedDict(:x => xs, :ν => νs)
yend = OrderedDict(:pA => ones(xn, νn), :pB => ones(xn, νn), :p => ones(xn, νn))
result = pdesolve(m, stategrid, yend; is_algebraic = OrderedDict(:pA => false, :pB => false, :p => true))
@assert result.residual_norm <= 1e-5