text:
export PseudoHilbert,Hilbert,SingularIntegral

#############
# Hilbert implements the Hilbert operator as a contour integral:
#
#        H f(z) := 1/π\int_Γ f(t)/(t-z) dt,  z ∈ Γ
#
# SingularIntegral implements the Hilbert operator as a line integral:
#
#       SI f(z) := 1/π\int_Γ f(t)/(t-z) ds(t),  z ∉ Γ
#
# PseudoHilbert corresponds to pseudocauchy, which may not be normalized at
# infinity.
#############

ApproxFun.@calculus_operator(PseudoHilbert)
ApproxFun.@calculus_operator(Hilbert)
ApproxFun.@calculus_operator(SingularIntegral)

for Op in (:PseudoHilbert,:Hilbert,:SingularIntegral)
    ConcOp=parse("Concrete"*string(Op))
    OpWrap=parse(string(Op)*"Wrapper")
    OffOp=parse("Off"*string(Op))
    @eval begin
        ## Convenience routines
        $Op(d::IntervalDomain,n::Int)=$Op(JacobiWeight(-.5,-.5,Chebyshev(d)),n)
        $Op(d::IntervalDomain)=$Op(JacobiWeight(-.5,-.5,Chebyshev(d)))
        $Op(d::PeriodicDomain,n::Int)=$Op(Laurent(d),n)
        $Op(d::PeriodicDomain)=$Op(Laurent(d))

        ## Modifiers for SumSpace, ArraySpace, and PiecewiseSpace

        $Op(AS::ArraySpace,n::Int)=$OpWrap(DiagonalArrayOperator($Op(AS.space,n),size(AS)),n)
        function $Op(S::PiecewiseSpace,n::Int)
            sp=vec(S)
            m=length(sp)
            diag=Any[$Op(sp[k],n) for k=1:m]
            C=Any[k==j?diag[k]:$OffOp(sp[k],rangespace(diag[j]),n) for j=1:m,k=1:m]
            D=Operator{mapreduce(i->eltype(C[i]),promote_type,1:m^2)}[C[j,k] for j=1:m,k=1:m]
            D=promotespaces(D)
            $OpWrap(InterlaceOperator(D,S,PiecewiseSpace(map(rangespace,D[:,1]))),n)
        end

        bandinds{s,DD}(::$ConcOp{Hardy{s,DD}})=0,0
        domainspace{s,DD}(H::$ConcOp{Hardy{s,DD}})=H.space
        rangespace{s,DD}(H::$ConcOp{Hardy{s,DD}})=H.space

        bandinds{DD}(H::$ConcOp{Fourier{DD}})=-H.order,H.order
        domainspace{DD}(H::$ConcOp{Fourier{DD}})=H.space
        rangespace{DD}(H::$ConcOp{Fourier{DD}})=H.space

        function rangespace{DD}(H::$ConcOp{JacobiWeight{Chebyshev{DD},DD}})
            @assert domainspace(H).α==domainspace(H).β==-0.5
            H.order==0?Chebyshev(domain(H)):Ultraspherical(H.order,domain(H))
        end
        function rangespace{DD}(H::$ConcOp{JacobiWeight{Ultraspherical{Int,DD},DD}})
            @assert order(domainspace(H)) == 1
            @assert domainspace(H).α==domainspace(H).β==0.5
            H.order==1?Chebyshev(domain(H)):Ultraspherical(H.order-1,domain(H))
        end
        # bandinds{λ,DD}(H::$ConcOp{JacobiWeight{Ultraspherical{λ,DD},DD}})=-λ,H.order-λ
        bandinds{DD}(H::$ConcOp{JacobiWeight{Chebyshev{DD},DD}}) = 0,H.order
        bandinds{DD}(H::$ConcOp{JacobiWeight{Ultraspherical{Int,DD},DD}}) =
            H.order > 0 ? (-1,H.order-1) : (-2,0)

        choosedomainspace(H::$Op{UnsetSpace},sp::Ultraspherical)=ChebyshevWeight(ChebyshevDirichlet{1,1}(domain(sp)))
        choosedomainspace(H::$Op{UnsetSpace},sp::PiecewiseSpace)=PiecewiseSpace(map(s->choosedomainspace(H,s),sp.spaces))

        # BlockOperator [1 Hilbert()] which allows for bounded solutions
        #TODO: Array values?
        choosedomainspace{T,V}(P::BlockOperator{$ConcOp{UnsetSpace,T,V}},
                               sp::Ultraspherical) =
                TupleSpace(ConstantSpace(),JacobiWeight(0.5,0.5,Ultraspherical(1,domain(sp))))

        choosedomainspace{T,V,W}(P::BlockOperator{ReOperator{$ConcOp{UnsetSpace,T,V},W}},
                               sp::Ultraspherical) =
                TupleSpace(ConstantSpace(),JacobiWeight(0.5,0.5,Ultraspherical(1,domain(sp))))
    end
end

for TYP in (:SumSpace,:PiecewiseSpace,:TupleSpace),(Op,OpWrap) in ((:PseudoHilbert,:PseudoHilbertWrapper),
                          (:Hilbert,:HilbertWrapper),
                          (:SingularIntegral,:SingularIntegralWrapper))
    @eval $Op(S::$TYP,k)=$OpWrap(DiagonalInterlaceOperator(map(s->$Op(s,k),S.spaces),$TYP),k)
end

# Length catch

ConcreteHilbert(sp::Space{ComplexBasis},n) =
    ConcreteHilbert{typeof(sp),typeof(n),Complex{real(eltype(domain(sp)))}}(sp,n)
ConcreteSingularIntegral(sp::Space{ComplexBasis},n) =
    ConcreteSingularIntegral{typeof(sp),typeof(n),Complex{real(eltype(domain(sp)))}}(sp,n)

# Override sumspace

for TYP in (:Hilbert,:SingularIntegral)
    ConcOp=parse("Concrete"*string(TYP))
    @eval function $TYP{DD<:Circle}(F::Fourier{DD},n)
        if !domain(F).orientation
            R=reverseorientation(F)
            Conversion(R,F)*(-$TYP(R,n))*Conversion(F,R)
        else
            $ConcOp(F,n)
        end
    end
end

### Operator Entries

## Circle

Hilbert{s,DD<:Circle}(S::Hardy{s,DD},m::Integer) = m==1?
    HilbertWrapper(SpaceOperator(ConstantOperator((s?1:-1)*im),S,S),m):
    ConcreteHilbert(S,m)

function getindex{DD<:Circle,OT,T}(H::ConcreteHilbert{Hardy{true,DD},OT,T},k::Integer,j::Integer)
    ##TODO: Add scale for different radii.
    m=H.order
    @assert m ≠ 1
    if k==j && m==0
        k==1?-T(2log(2)):one(T)/(k-1)
    elseif k==j && k ≠ 1
        k==im*(one(T)*im*(k-1))^(m-1)
    else
        zero(T)
    end
end

function getindex{DD<:Circle,OT,T}(H::ConcreteHilbert{Hardy{false,DD},OT,T},k::Integer,j::Integer)
    ##TODO: Add scale for different radii.
    m=H.order
    @assert m ≠ 1
    if k==j
        -im*(one(T)*im*k)^(m-1)
    else
        zero(T)
    end
end

function Hilbert{DD<:Circle}(S::Fourier{DD},m::Integer)
    @assert m==0 || m==1
    ConcreteHilbert(S,m)
end

function getindex{DD<:Circle,OT,T}(H::ConcreteHilbert{Fourier{DD},OT,T},k::Integer,j::Integer)
    d = domain(H)
    r = d.radius
    o = d.orientation
    if H.order == 0 && k==j
        T(k==1?2r*log(r):(-r/(k÷2)))
    elseif H.order == 1
        if k==j==1
            T(o?im:-im)
        elseif iseven(k) && j==k+1
            -one(T)
        elseif isodd(k) && j==k-1
            one(T)
        else
            zero(T)
        end
    else
        zero(T)
    end
end

function addentries!{DD<:Circle}(H::ConcreteSingularIntegral{Hardy{true,DD}},A,kr::Range,::Colon)
##TODO: Add scale for different radii.
    m=H.order
    d=domain(H)
    sp=domainspace(H)

    r = domain(H).radius
    if m == 0
        for k=kr
            A[k,k] += k==1?2r*log(r):-r./(k-1)
        end
    elseif m == 1
        for k=kr
            A[k,k] += im
        end
    else
        for k=kr
            A[k,k] += k==1?0.0:1.0im*(1.0im*(k-1))^(m-1)
        end
    end
    A
end

function addentries!{DD<:Circle}(H::ConcreteSingularIntegral{Hardy{false,DD}},A,kr::Range,::Colon)
##TODO: Add scale for different radii.
    m=H.order
    d=domain(H)
    sp=domainspace(H)

    r = domain(H).radius
    if m== 1
        for k=kr
            A[k,k]-= im
        end
    else
        for k=kr
            A[k,k]-=1.0im*(1.0im*k/r)^(m-1)
        end
    end
    A
end

function addentries!{DD<:Circle}(H::ConcreteSingularIntegral{Fourier{DD}},A,kr::Range,::Colon)
    r = domain(H).radius
    if H.order == 0
        for k=kr
            if k==1
                A[1,1]+=2r*log(r)
            else
                j=div(k,2)
                A[k,k]+=-r/j
            end
        end
    elseif H.order == 1
        for k=kr
            if k==1
                A[1,1]+=0
            elseif iseven(k)
                A[k,k+1]-=1
            else   #isodd(k)
                A[k,k-1]+=1
            end
        end
    else
            error("Hilbert order $(H.order) not implemented for Fourier")
    end

    A
end

## JacobiWeight

for (Op,Len) in ((:Hilbert,:complexlength),
                        (:SingularIntegral,:arclength))
    ConcOp=parse("Concrete"*string(Op))
    OpWrap=parse(string(Op)*"Wrapper")

    @eval begin
        function $Op{DD<:Interval}(S::JacobiWeight{Chebyshev{DD},DD},m::Int)
            if S.α==S.β==-0.5
                if m==0
                    $ConcOp(S,m)
                else
                    d=domain(S)
                    C=(4./$Len(d))^(m-1)
                    $OpWrap(SpaceOperator(
                        ToeplitzOperator([0.],[zeros(m);C]),S,Ultraspherical(m,d)),m)
                end
            elseif S.α==S.β==0.5
                d=domain(S)
                if m==1
                    J=JacobiWeight(0.5,0.5,Ultraspherical(1,d))
                    $OpWrap($Op(J,m)*Conversion(S,J),m)
                else
                    J=JacobiWeight(-0.5,-0.5,Chebyshev(d))
                    $OpWrap($Op(J,m)*Conversion(S,J),m)
                end
            else
                error(string($Op)*" not implemented for parameters $(S.α),$(S.β)")
            end
        end

        function getindex{DD<:Interval,OT,T}(H::$ConcOp{JacobiWeight{Chebyshev{DD},DD},OT,T},k::Integer,j::Integer)
            sp=domainspace(H)
            @assert H.order == 0
            @assert sp.α==sp.β==-0.5

            if k==j
                d=domain(H)
                C=$Len(d)/4
                -T(k==1?-2C*log(C):2C/(k-1))
            else
                zero(T)
            end
        end

        # we always have real for n==1
        function $Op{DD<:Interval}(S::JacobiWeight{Ultraspherical{Int,DD},DD},m)
            @assert order(S.space) == 1
            if S.α==S.β==0.5
                if m==1
                    d=domain(S)
                    $OpWrap(SpaceOperator(
                        ToeplitzOperator([-1.],[0.]),S,m==1?Chebyshev(d):Ultraspherical(m-1,d)),m)
                else
                    $ConcOp(sp,m)
                end
            else
                error(string($Op)*" not implemented for parameters $(S.α),$(S.β)")
            end
        end
        function getindex{DD<:Interval,OT,T}(H::$ConcOp{JacobiWeight{Ultraspherical{Int,DD},DD},OT,T},k::Integer,j::Integer)
            # order(domainspace(H))==1
            m=H.order
            d=domain(H)
            sp=domainspace(H)
            @assert sp.α==sp.β==0.5
            @assert m ≠ 1

            C=(4./$Len(d))^(m-1)
            if m == 0
                if k==j==1
                    T(C*log(C))
                elseif k==j
                    -T(C/(k-1))
                elseif j==k-2
                    T(C/(k-1))
                else
                    zero(T)
                end
            elseif j==k+m-2
                -T(.5C*k/(m-1))
            else
                zero(T)
            end
        end
    end
end

## PseudoHilbert
# The default is Hilbert

# getindex(H::ConcretePseudoHilbert,k::Integer,j::Integer)=Hilbert(H.space,H.order)[k,j]
# bandinds(H::ConcretePseudoHilbert)=bandinds(Hilbert(H.space,H.order))


text:
using Test
using Random: seed!

using NearestNeighborsPBC

function neighbor_sets_equal(S, T)
    return all(Set(s) == Set(t) for (s, t) in zip(S, T))
end

@testset "test NearestNeighborsPBC against naive N^2 implementation. with r_max passed" begin
    seed!(2354523)
    L = 10.0
    res = []
    for _i in 1:20
        for r_max in [0.1, 0.3, 1.0, 1.5, 9.]
            for ndim in [2, 3]
                for n  in [4, 10]
                    X = L .* rand(Float64, (ndim, n,))

                    ind1 = inrangePBC_naive(X, L, r_max)

                    kdt = KDTreePBC(X, L, r_max=r_max)
                    ind2 = inrange(kdt, X, r_max)

                    push!(res, neighbor_sets_equal(ind1, ind2))
                end
            end
        end
    end
    @test all(res)
end

@testset "test NearestNeighborsPBC against naive N^2 implementation. without passing r_max and using Float32's." begin
    seed!(2354524)
    L = 6.4f0
    res = []
    for _i in 1:20
        for r_max in [0.1f0, 0.3f0, 1.0f0, 1.5f0, 4.0f0]
            for ndim in [2, 3]
                for n  in [4, 10]
                    X = L .* rand(Float32, (ndim, n,))

                    ind1 = inrangePBC_naive(X, L, r_max)

                    kdt = KDTreePBC(X, L)
                    ind2 = inrange(kdt, X, r_max)

                    push!(res, neighbor_sets_equal(ind1, ind2))
                end
            end
        end
    end
    @test all(res)
end

ext: jl

max_stars_repo_name: gregrf/Geo

max_issues_repo_name: gregrf/Geo

max_forks_repo_name: gregrf/Geo


text:
using NLsolve, Distributions

@doc raw"""
	CapitalAdjustModel2(; kwargs...)

Construct capital investment model with two types of capital and one
productivity shock and convex adjustment costs. The recursive formulation is:

```math
\begin{equation*}
V(K,N,z) = \max_{a_K,a_N} K^\alpha e^z - i_K K - \frac{\gamma_K}{2} i_K^2 K - i_N N - \frac{\gamma_N}{2} i_N^2 N + \beta \mathbb{E} V(K',N',z') \\
i_K \equiv \frac{a_K}{K} - (1-\delta_K) \\
i_N \equiv \frac{a_N}{N} - (1-\delta_N) \\
\text{where } K' = a_K \\
N' = a_N \\
z' = \rho z + \sigma \varepsilon, \quad \varepsilon \sim \mathcal{N}(0,1)
\end{equation*}
```

The firm solves the following problem to initialize at t=0:

```math
\begin{equation*}
V_0(z_0) = \max_{K_0,N_0} V(K_0,N_0,z_0) - (1 + C0) (K_0 + N_0) \\
\text{where } z_0 = \sqrt{\frac{\sigma^2}{1-\rho^2}} ε_0 \\
ε_0 \sim \mathcal{N}(0,1)
\end{equation*}
```

Example:

```julia
using StructuralDDP
prob = CapitalAdjustModel2(nK=30, nN=25, nz=5, ρ=0.5, σ=0.3, δ_K=0.1)
```

""" CapitalAdjustModel2
function CapitalAdjustModel2(;
	β   = 0.90,
	α  = 0.67,
	η  = 0.6,
	δ_K  = 0.1,
	δ_N  = 0.2,
	γ_K  = 1.,
	γ_N  = 2.,
	ρ  = 0.6,
	σ  = 0.3,
	C0  = 0.4,
	F = 0.,

	nK = 50,
	nN = 50,
	nz = 3,

	# options for testing
	intdim = :Separable_ExogStates,
	initialize_exact = true)

	# initial standard deviation of transitiory shock
    stdz = sqrt(σ^2/(1-ρ^2))
    minz = -3*stdz
    maxz =  3*stdz
    vz = collect(LinRange(minz, maxz, nz))

	function steadystate()
        # check https://www3.nd.edu/~esims1/investment_notes.pdf for details

        function f!(F,x)
            (K, N, I_K, I_N, V_K, V_N) = x
            K = abs(K)
            N = abs(N)

            F[1] = I_K - δ_K * K
            F[2] = I_N - δ_N * N
            F[3] = - 1 - γ_K * I_K/(K+N) + β*V_K
			F[4] = - (1 + γ_N * I_N/(K+N)) + β*V_N
            adjcosts = γ_K/2 * (I_K/(K+N))^2 + γ_N/2 * (I_N/(K+N))^2
            F[5] = - V_K + β*( (exp(stdz^2/2) * α*η*K^(α*η-1)*N^((1-α)*η) +
                                adjcosts) + (1-δ_K)*V_K)
            F[6] = - V_N + β*( (exp(stdz^2/2) * (1-α)*η*K^(α*η)*N^((1-α)*η-1) +
                                adjcosts) + (1-δ_N)*V_N)
        end

        out = nlsolve(f!, [ 10.; 10.; 0.1; 0.1; 1.1; 1.1])
        return (abs(out.zero[1]), abs(out.zero[2]))
    end

    (K_ss, N_ss) = steadystate()

	minK = 0.1 * K_ss
    maxK = 3.5 * K_ss
    vK   = exp.(collect(LinRange(log(minK), log(maxK), nK)))

	minN = 0.1 * N_ss
    maxN = 3.5 * N_ss
 	vN   = exp.(collect(LinRange(log(minN), log(maxN), nN)))

    tStateVectors = (vK, vN, vz)

	function rewardfunc_partial(vStates)
		(K, N, z) = vStates
		return (K^α * N^(1-α))^η * exp(z)
	end

	function adjustcosts(Kprime, K, Nprime, N)

	    capx = Kprime - (1-δ_K)*K
	    intx = Nprime - (1-δ_N)*N

	    Ktot = K + N

	    adj_K = γ_K/2*(capx/Ktot)^2 * Ktot
	    adj_N = γ_N/2*(intx/Ktot)^2 * Ktot

		action = Kprime != K

	    return adj_K + adj_N + F*K*action
	end

	function reward(vStates, vChoices)
		(K, N, z) = vStates
		(Kprime, Nprime) = vChoices

		capx = Kprime - (1-δ_K)*K
		intx = Nprime - (1-δ_N)*N

		return rewardfunc_partial(vStates) - capx - intx -
									adjustcosts(Kprime, K, Nprime, N)
	end

	function reward(Output, vStates, vChoices)
		K, N = vStates
		Kprime, Nprime = vChoices
	    capx = Kprime - (1-δ_K)*K
	    intx = Nprime - (1-δ_N)*N
	    adj = adjustcosts(Kprime, K, Nprime, N)
	    return Output - adj - intx - capx
	end

	if intdim == :All
	    transfunc = function mytransfunc(vStates, vChoices, Shock)
	        z = vStates[3]
	        zprime  = ρ*z + σ * Shock
	        return vChoices[1], vChoices[2], zprime
	    end
		tChoiceVectors = (vK, vN)

	elseif intdim == :Separable
	    transfunc = function mytransfunc1(vStates, vChoices, Shock)
	        z = vStates[3]
	        zprime  = ρ*z + σ * Shock
	        return zprime
	    end
		tChoiceVectors = (1, 2)

	elseif intdim == :Separable_States
	    transfunc = function mytransfunc2(vStates, Shock)
	        z = vStates[3]
	        zprime  = ρ*z + σ * Shock
	        return zprime
	    end
		tChoiceVectors = (1, 2)

	elseif intdim == :Separable_ExogStates
	    transfunc = function mytransfunc3(z, Shock)
	        zprime  = ρ*z + σ * Shock
	        return zprime
	    end
		tChoiceVectors = (1, 2)

	else
		error("$intdim wrong intdim")
	end

	if initialize_exact
		tChoiceVectorsZero = (1,2)

		initializationproblem(value, vChoices) = value - (1+C0)*(vChoices[1] + vChoices[2])

		function initialize(shock)
		    z0 = shock * sqrt(σ^2 / (1-ρ^2))
			# z0 = inbounds(z0, tStateVectors[3][1], tStateVectors[3][end])
			return z0
		end
		shockdist_initial = Normal()
	else
		tChoiceVectorsZero = nothing
		initializationproblem = nothing
		initialize = nothing
		shockdist_initial = nothing
	end

	DDP(tStateVectors,
		tChoiceVectors,
		reward,
        transfunc,
		Normal(), # distribution of shocks: standard normal
		β;
		intdim = intdim,
        rewardfunc_partial = rewardfunc_partial,
        initializationproblem = initializationproblem,
        initializefunc = initialize,
		shockdist_initial = shockdist_initial,
		tChoiceVectorsZero = tChoiceVectorsZero,
        )
end


text:
"""
$(TYPEDEF)

Struct to hold sparse matrix in the linked list format.

Modeled after the linked list sparse matrix format described in
the  [whitepaper](https://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps)
and the  [SPARSEKIT2 source code](https://www-users.cs.umn.edu/~saad/software/SPARSKIT/SPARSKIT2.tar.gz)
by Y. Saad. He writes "This is one of the oldest data structures used for sparse matrix computations."

The relevant source [formats.f](https://salsa.debian.org/science-team/sparskit/blob/master/FORMATS/formats.f)
is also available in the debian/science gitlab.

The advantage of the linked list structure is the fact that upon insertion
of a new entry, the arrays describing the structure can grow at their respective ends and
can be conveniently updated via `push!`.  No copying of existing data is necessary.

$(TYPEDFIELDS)
"""
mutable struct SparseMatrixLNK{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}

    """
    Number of rows
    """
    m::Ti

    """
    Number of columns
    """
    n::Ti

    """
    Number of nonzeros
    """
    nnz::Ti

    """
    Length of arrays
    """
    nentries::Ti

    """
    Linked list of column entries. Initial length is n,
    it grows with each new entry.

    colptr[index] contains the next
    index in the list or zero, in the later case terminating the list which
    starts at index 1<=j<=n for each column j.
    """
    colptr::Vector{Ti}

    """
    Row numbers. For each index it contains the zero (initial state)
    or the row numbers corresponding to the column entry list in colptr.

    Initial length is n,
    it grows with each new entry.
    """
    rowval::Vector{Ti}

    """
    Nonzero entry values correspondin to each pair
    (colptr[index],rowval[index])

    Initial length is n,  it grows with each new entry.
    """
    nzval::Vector{Tv}
end

"""
$(SIGNATURES)

Constructor of empty matrix.
"""
SparseMatrixLNK{Tv,Ti}(m, n)  where {Tv,Ti<:Integer} =    SparseMatrixLNK{Tv,Ti}(m,n,0,n,zeros(Ti,n),zeros(Ti,n),zeros(Tv,n))

"""
$(SIGNATURES)

Constructor of empty matrix.
"""
SparseMatrixLNK(valuetype::Type{Tv},indextype::Type{Ti},m, n)  where {Tv,Ti<:Integer} =SparseMatrixLNK{Tv,Ti}(m,n)

"""
$(SIGNATURES)

Constructor of empty matrix.
"""
SparseMatrixLNK(valuetype::Type{Tv},m, n)  where {Tv} =SparseMatrixLNK(Tv,Int,m,n)

"""
$(SIGNATURES)

Constructor of empty matrix.
"""
SparseMatrixLNK(m, n)  where {Tv} =SparseMatrixLNK(Float64,m,n)

"""
$(SIGNATURES)

Constructor from SparseMatrixCSC.

"""
function SparseMatrixLNK(csc::SparseArrays.SparseMatrixCSC{Tv,Ti})  where {Tv,Ti<:Integer}
    lnk=SparseMatrixLNK{Tv,Ti}(csc.m,csc.n)
    for j=1:csc.n
        for k=csc.colptr[j]:csc.colptr[j+1]-1
            lnk[csc.rowval[k],j]=csc.nzval[k]
        end
    end
    lnk
end

function findindex(lnk::SparseMatrixLNK,i,j)
    if !((1 <= i <= lnk.m) & (1 <= j <= lnk.n))
        throw(BoundsError(lnk, (i,j)))
    end

    k=j
    k0=j
    while k>0
        if lnk.rowval[k]==i
            return k,0
        end
        k0=k
        k=lnk.colptr[k]
    end
    return 0,k0
end

"""
$(SIGNATURES)

Return value stored for entry or zero if not found
"""
function Base.getindex(lnk::SparseMatrixLNK{Tv,Ti},i,j) where {Tv,Ti}
    k, k0=findindex(lnk,i,j)
    if k==0
        return zero(Tv)
    else
        return lnk.nzval[k]
    end
end

function addentry!(lnk::SparseMatrixLNK,i,j,k,k0)
    # increase number of entries
    lnk.nentries+=1
    if length(lnk.nzval)<lnk.nentries
        newsize=Int(ceil(5.0*lnk.nentries/4.0))
        resize!(lnk.nzval,newsize)
        resize!(lnk.rowval,newsize)
        resize!(lnk.colptr,newsize)
    end

    # Append entry if not found
    lnk.rowval[lnk.nentries]=i

    # Shift the end of the list
    lnk.colptr[lnk.nentries]=0
    lnk.colptr[k0]=lnk.nentries

    # Update number of nonzero entries
    lnk.nnz+=1
    return lnk.nentries
end

"""
$(SIGNATURES)

Update value of existing entry, otherwise extend matrix if v is nonzero.
"""
function Base.setindex!(lnk::SparseMatrixLNK,v,i,j)
    if !((1 <= i <= lnk.m) & (1 <= j <= lnk.n))
        throw(BoundsError(lnk, (i,j)))
    end

    # Set the first  column entry if it was not yet set.
    if lnk.rowval[j]==0 && !iszero(v)
        lnk.rowval[j]=i
        lnk.nzval[j]=v
        lnk.nnz+=1
        return lnk
    end

    k,k0=findindex(lnk,i,j)
    if k>0
        lnk.nzval[k]=v
        return lnk
    end
    if !iszero(v)
        k=addentry!(lnk,i,j,k,k0)
        lnk.nzval[k]=v
    end
    return lnk
end

"""
$(SIGNATURES)

Update element of the matrix  with operation `op`.
It assumes that `op(0,0)==0`. If `v` is zero, no new
entry is created.
"""
function updateindex!(lnk::SparseMatrixLNK{Tv,Ti},op, v, i,j) where {Tv,Ti}
    # Set the first  column entry if it was not yet set.
    if lnk.rowval[j]==0 && !iszero(v)
        lnk.rowval[j]=i
        lnk.nzval[j]=op(lnk.nzval[j],v)
        lnk.nnz+=1
        return lnk
    end
    k,k0=findindex(lnk,i,j)
    if k>0
        lnk.nzval[k]=op(lnk.nzval[k],v)
        return lnk
    end
    if !iszero(v)
        k=addentry!(lnk,i,j,k,k0)
        lnk.nzval[k]=op(zero(Tv),v)
    end
    lnk
end

"""
$(SIGNATURES)

Update element of the matrix  with operation `op`.
It assumes that `op(0,0)==0`. If `v` is zero a new entry
is created nevertheless.
"""
function rawupdateindex!(lnk::SparseMatrixLNK{Tv,Ti},op, v, i,j) where {Tv,Ti}
    # Set the first  column entry if it was not yet set.
    if lnk.rowval[j]==0
        lnk.rowval[j]=i
        lnk.nzval[j]=op(lnk.nzval[j],v)
        lnk.nnz+=1
        return lnk
    end
    k,k0=findindex(lnk,i,j)
    if k>0
        lnk.nzval[k]=op(lnk.nzval[k],v)
    else
        k=addentry!(lnk,i,j,k,k0)
        lnk.nzval[k]=op(zero(Tv),v)
    end
    lnk
end

"""
$(SIGNATURES)

Return tuple containing size of the matrix.
"""
Base.size(lnk::SparseMatrixLNK) = (lnk.m, lnk.n)

"""
$(SIGNATURES)

Return number of nonzero entries.
"""
SparseArrays.nnz(lnk::SparseMatrixLNK)=lnk.nnz

"""
$(SIGNATURES)

Dummy flush! method for SparseMatrixLNK. Just
used in test methods
"""
function flush!(lnk::SparseMatrixLNK{Tv, Ti}) where{Tv, Ti}
    return lnk
end

# Struct holding pair of value and row
# number, for sorting
struct ColEntry{Tv,Ti<:Integer}
    rowval::Ti
    nzval::Tv
end

# Comparison method for sorting
Base.isless(x::ColEntry, y::ColEntry)= (x.rowval<y.rowval)

"""
$(SIGNATURES)

Add SparseMatrixCSC matrix and [`SparseMatrixLNK`](@ref)  lnk, returning a SparseMatrixCSC
"""
function Base.:+(lnk::SparseMatrixLNK{Tv,Ti},csc::SparseMatrixCSC)::SparseMatrixCSC where {Tv,Ti<:Integer}
    @assert(csc.m==lnk.m)
    @assert(csc.n==lnk.n)

    # overallocate arrays in order to avoid
    # presumably slower push!
    xnnz=nnz(csc)+nnz(lnk)
    colptr=Vector{Ti}(undef,csc.n+1)
    rowval=Vector{Ti}(undef,xnnz)
    nzval=Vector{Tv}(undef,xnnz)

    # Detect the maximum column length of lnk
    lnk_maxcol=0
    for j=1:csc.n
        lcol=zero(Ti)
        k=j
        while k>0
            lcol+=1
            k=lnk.colptr[k]
        end
        lnk_maxcol=max(lcol,lnk_maxcol)
    end

    # pre-allocate column  data
    col=[ColEntry{Tv,Ti}(0,0) for i=1:lnk_maxcol]

    inz=1 # counts the nonzero entries in the new matrix

    in_csc_col(jcsc,j)=(nnz(csc)>zero(Ti)) && (jcsc<csc.colptr[j+1])

    in_lnk_col(jlnk,l_lnk_col)=(jlnk<=l_lnk_col)

    # loop over all columns
    for j=1:csc.n
        # Copy extension entries into col and sort them
        k=j
        l_lnk_col=0
        while k>0
            if lnk.rowval[k]>0
                l_lnk_col+=1
                col[l_lnk_col]=ColEntry(lnk.rowval[k],lnk.nzval[k])
            end
            k=lnk.colptr[k]
        end
        sort!(col,1,l_lnk_col, Base.QuickSort, Base.Forward)

        # jointly sort lnk and csc entries  into new matrix data
        # this could be replaced in a more transparent manner by joint sorting:
        # make a joint array for csc and lnk col, sort them.
        # Will this be faster?

        colptr[j]=inz
        jlnk=one(Ti) # counts the entries in col
        jcsc=csc.colptr[j]  # counts entries in csc

        while true
            if in_csc_col(jcsc,j) &&  (in_lnk_col(jlnk,l_lnk_col)  && csc.rowval[jcsc]<col[jlnk].rowval  ||  !in_lnk_col(jlnk,l_lnk_col))
                # Insert entries from csc into new structure
                rowval[inz]=csc.rowval[jcsc]
                nzval[inz]=csc.nzval[jcsc]
                jcsc+=1
                inz+=1
            elseif in_csc_col(jcsc,j) &&  (in_lnk_col(jlnk,l_lnk_col)  && csc.rowval[jcsc]==col[jlnk].rowval)
                # Add up entries from csc and lnk
                rowval[inz]=csc.rowval[jcsc]
                nzval[inz]=csc.nzval[jcsc]+col[jlnk].nzval
                jcsc+=1
                inz+=1
                jlnk+=1
            elseif in_lnk_col(jlnk,l_lnk_col)
                # Insert entries from lnk res. col into new structure
                rowval[inz]=col[jlnk].rowval
                nzval[inz]=col[jlnk].nzval
                jlnk+=1
                inz+=1
            else
                break
            end
        end
    end
    colptr[csc.n+1]=inz
    # Julia 1.7 wants this correct
    resize!(rowval,inz-1)
    resize!(nzval,inz-1)
    SparseMatrixCSC{Tv,Ti}(csc.m,csc.n,colptr,rowval,nzval)
end

Base.:+(csc::SparseMatrixCSC,lnk::SparseMatrixLNK)=lnk+csc

"""
$(SIGNATURES)

Constructor from SparseMatrixLNK.

"""
function SparseArrays.SparseMatrixCSC(lnk::SparseMatrixLNK)::SparseMatrixCSC
    csc=spzeros(lnk.m,lnk.n)
    lnk+csc
end

rowvals(S::SparseMatrixLNK) = getfield(S, :rowval)
getcolptr(S::SparseMatrixLNK) = getfield(S, :colptr)
nonzeros(S::SparseMatrixLNK) = getfield(S, :nzval)

copy(S::SparseMatrixLNK) = SparseMatrixLNK(size(S, 1), size(S, 2), S.nnz, S.nentries, copy(getcolptr(S)), copy(rowvals(S)), copy(nonzeros(S)))


text:
using AbstractFFTs
using Test
using Statistics
using StaticRanges
using NamedDims
using IntervalSets
using FFTW
using LinearAlgebra
using MappedArrays
using Dates
using Documenter
using StaticArrays
using ArrayInterface
using ArrayInterface: indices, known_length, StaticInt

#=
using Dates,MappedArrays,Statistics,LinearAlgebra,Base,Core
pkgs = (Dates,MappedArrays,Statistics,LinearAlgebra,Base,Core);
ambs = detect_ambiguities(pkgs...);
using AxisIndices
ambs = setdiff(detect_ambiguities(AxisIndices, pkgs...), ambs);

unique([i[1].name for i in ambs])
for amb_i in ambs
    if amb_i[1].name == :CartesianIndices
        print(amb_i)
    end
end

inds = [i[1].name == :copyto! for i in ambs]
ambs[[i[1].name == :mapreduce for i in ambs]]

itr[1].name
unique([itr[1].name for itr in ambs])

if !isempty(ambs)
    println("Ambiguities:")
    for a in ambs
        println(a)
    end
=#

using DelimitedFiles
using AxisIndices
using AxisIndices: is_key, cat_axis, hcat_axes, vcat_axes, matmul_axes
using AxisIndices: CenteredAxis, IdentityAxis, OffsetAxis
using StaticRanges: can_set_first, can_set_last, can_set_length, parent_type
using StaticRanges: grow_last, grow_last!, grow_first, grow_first!
using StaticRanges: shrink_last, shrink_last!, shrink_first, shrink_first!

using ArrayInterface: to_axes, to_index
using Base: step_hp, OneTo
using Base.Broadcast: broadcasted
bstyle = Base.Broadcast.DefaultArrayStyle{1}()

@test Base.to_shape(SimpleAxis(1)) == 1

include("simple_axis.jl")
include("axis.jl")
include("abstract_axis.jl")
include("centered_axis.jl")
include("identity_axis.jl")
include("offset_axis.jl")
include("struct_axis.jl")
include("padded_axis.jl")
include("arrays.jl")
include("indexing.jl")
include("append_tests.jl")
include("similar_tests.jl")
include("broadcasting_tests.jl")
include("linear_algebra.jl")
include("copyto_tests.jl")
include("permutedims.jl")
include("closest.jl")

@testset "CartesianAxes" begin
    cartaxes = CartesianAxes((2.0:5.0, 1:4))
    cartinds = CartesianIndices((1:4, 1:4))
    for (axs, inds) in zip(collect(cartaxes), collect(cartinds))
        @test axs == inds
    end

    for (axs, inds) in zip(cartaxes, cartinds)
        @test axs == inds
    end

    @test collect(cartaxes) == cartaxes[1:4, 1:4]
end

@testset "LinearAxes" begin
    linaxes = LinearAxes((2.0:5.0, 1:4))
    lininds = LinearIndices((1:4, 1:4))
    @test @inferred(linaxes[10]) == 10
    for (axs, inds) in zip(collect(linaxes), collect(lininds))
        @test axs == inds
    end

    for (axs, inds) in zip(linaxes, lininds)
        @test axs == inds
    end
    @test collect(linaxes) == linaxes[1:4,1:4]
end

@test first(axes(@inferred(filter(isodd, AxisArray(1:7, (2:8,)))), 1)) == 2
@test axes(filter(isodd, AxisArray([11 12; 21 22], (2:3, 3:4))), 1) isa AbstractAxis
include("mapped_arrays.jl")
#include("offset_array_tests.jl")
include("resize_tests.jl")
include("fft.jl")

if VERSION >= v"1.6.0-DEV.421"
    @testset "docs" begin
        doctest(AxisIndices)
    end
end

#include("NamedMetaAxisArray_tests.jl")


text:
using ClusterManagers, DSGE, ModelConstructors, FileIO, Plots, Dates, OrderedCollections, Distributed
using Plots.PlotMeasures
GR.inline("pdf")
fn = dirname(@__FILE__)

# This script shows how to use regime-switching
# to enact a temporary ZLB with imperfect awareness, followed by an
# alternative policy with imperfect awareness (in this case, AIT).
# We sometimes also call policies with imperfect awareness as "uncertain" policies.
# This script demonstrates a specific case of the general regime switching code
# demonstrated in regime_switching.jl.

# What do you want to do?
zlb_altpolicy     = true   # run a forecast
make_plots        = false  # Make plots
save_plots        = false  # Save plots to figurespath(m, "forecast")
endo_zlb          = false  # Run forecast with endogenous zlb
add_workers       = false  # Run in parallel
n_workers         = 10

pol_str = "Uncertain ZLB, Uncertain AIT"

# Initialize model objects and desired settings
custom_settings = Dict{Symbol, Setting}(:n_mon_anticipated_shocks =>
                                        Setting(:n_mon_anticipated_shocks, 6, "Number of anticipated policy shocks"),
                                        :add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true),
                                        :add_altpolicy_yagp => Setting(:add_altpolicy_ygap, true),
                                        :flextible_ait_policy_change =>
                                        Setting(:flexible_ait_policy_change, false))
m = Model1002("ss10"; custom_settings = custom_settings)   # We will directly construct the matrix of parameter draws
m <= Setting(:flexible_ait_2020Q3_policy_change, false)

data_vint = "201117"
cond_vint = "201119"
fcast_date = DSGE.quartertodate("2020-Q4")
cond_date = DSGE.quartertodate("2020-Q4")
date_fcast_end = iterate_quarters(fcast_date, 40)

m10 = Model1002("ss10") # b/c loading from a saved estimation file has not been fully implemented

for models in [m, m10]
    # Not strictly necessary to have the same settings, but it makes it easier for us
    usual_model_settings!(m, "201117"; cdvt = "201119", fcast_date = quartertodate("2020-Q4"))
    m <= Setting(:use_population_forecast, false)
    m <= Setting(:saveroot, joinpath("$(fn)", "..", "save"))
    m <= Setting(:dataroot, joinpath("$(fn)", "..", "save", "input_data"))
    m <= Setting(:date_forecast_end, date_fcast_end)
    m <= Setting(:forecast_block_size, 40)
    if get_setting(m, :sampling_method) == :SMC
        m <= Setting(:forecast_jstep, 1)
    end
end

# Set up parameters for regime-switching.
overrides = forecast_input_file_overrides(m10)
overrides[:full] = "$(fn)/../test/reference/mhsave_vint=181115.h5"
modal_params = map(x -> x.value, m.parameters) # The default parameters will be our "modal" parameters for this exercise

forecast_string = ""

# Set up regime switching and forecast settings
m <= Setting(:n_regimes, 9, true, "reg", "")   # How many regimes?
m <= Setting(:regime_switching, true)          # Need to set that the model does have regime-switching
m <= Setting(:forecast_block_size, 5000)
m <= Setting(:use_parallel_workers, true)
baseline_regime_dates = Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2020, 9, 30))
m <= Setting(:regime_dates, baseline_regime_dates)
setup_regime_switching_inds!(m)

# Some settings for alternative policy
if zlb_altpolicy
    # parametrize the alternative policy (AIT, in this case)
    m <= Setting(:pgap_value, 0.)
    m <= Setting(:pgap_type, :smooth_ait_gdp_alt)
    m <= Setting(:ygap_value, 12.)
    m <= Setting(:ygap_type, :smooth_ait_gdp_alt)
    m <= Setting(:smooth_ait_gdp_alt_ρ_smooth, 0.)
    m <= Setting(:ait_Thalf, 10.)
    m <= Setting(:gdp_Thalf, 10.)
    m <= Setting(:smooth_ait_gdp_alt_φ_y, 6.)
    m <= Setting(:smooth_ait_gdp_alt_φ_π, 6.)
end

# Some settings for the endogenous zlb (if we're using one)
if endo_zlb
    # Set the maximum endogenous zlb length (in the forecast)
    m <= Setting(:max_temporary_altpolicy_length, 20)
    # Use only two regimes of ZLB if nominal rates do not drop below the zlb
    # (in the forecast)
    m <= Setting(:min_temporary_altpolicy_length, 2)
    # If we have zlb in the history, we indicate the length of pre-forecast zlb using
    # the following setting (not applicable in this example)
    # m <= Setting(:min_temporary_altpolicy_length, length_prefcast_zlb)
    # If we have time varying credibility, we specify the length of the credibility
    # variation as follows:
    m <= Setting(:cred_vary_until, 9)
end

for i in 1:2
    adj = (i == 2) ? .9 : 1. # adjustment to the size of the standard deviation

    # For each parameter, we need to instruct it that a regime exists.
    # Note that the regimes MUST be entered in order of the regimes
    # because regime values are stored in an OrderedDict, which is sorted
    # by insertion order.
    ModelConstructors.set_regime_val!(m[:σ_g], i, adj * m10[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_b], i, adj * m10[:σ_b].value)
    ModelConstructors.set_regime_val!(m[:σ_μ], i, adj * m10[:σ_μ].value)
    ModelConstructors.set_regime_val!(m[:σ_ztil], i, adj * m10[:σ_ztil].value)
    ModelConstructors.set_regime_val!(m[:σ_λ_f], i, adj * m10[:σ_λ_f].value)
    ModelConstructors.set_regime_val!(m[:σ_λ_w], i, adj * m10[:σ_λ_w].value)
    ModelConstructors.set_regime_val!(m[:σ_r_m], i, adj * m10[:σ_r_m].value)
    ModelConstructors.set_regime_val!(m[:σ_σ_ω], i, adj * m10[:σ_σ_ω].value)
    ModelConstructors.set_regime_val!(m[:σ_μ_e], i, adj * m10[:σ_μ_e].value)
    ModelConstructors.set_regime_val!(m[:σ_γ], i, adj * m10[:σ_γ].value)
    ModelConstructors.set_regime_val!(m[:σ_π_star], i, adj * m10[:σ_π_star].value)
    ModelConstructors.set_regime_val!(m[:σ_lr], i, adj * m10[:σ_lr].value)
    ModelConstructors.set_regime_val!(m[:σ_z_p], i, adj * m10[:σ_z_p].value)
    ModelConstructors.set_regime_val!(m[:σ_tfp], i, adj * m10[:σ_tfp].value)
    ModelConstructors.set_regime_val!(m[:σ_gdpdef], i, adj * m10[:σ_gdpdef].value)
    ModelConstructors.set_regime_val!(m[:σ_corepce], i, adj * m10[:σ_corepce].value)
    ModelConstructors.set_regime_val!(m[:σ_gdp], i, adj * m10[:σ_gdp].value)
    ModelConstructors.set_regime_val!(m[:σ_gdi], i, adj * m10[:σ_gdi].value)

    for j = 1:DSGE.n_mon_anticipated_shocks(m)
        ModelConstructors.set_regime_val!(m[Symbol("σ_r_m$(j)")], i, adj * m10[Symbol("σ_r_m$(j)")])
    end
end

df = load_data(m; check_empty_columns = false)

output_vars = [:histpseudo, :histobs, :histstdshocks,
               :hist4qpseudo, :hist4qobs, :histutpseudo,
               :forecastpseudo, :forecastobs, :forecastutpseudo,
               :forecast4qpseudo, :forecast4qobs, :forecaststdshocks]
if endo_zlb
    append!(output_vars, [:bddforecastobs, :bddforecastpseudo])
end

if zlb_altpolicy
    Hbar = 6 # number of zlb regimes
    n_tempZLB_regimes = 3 + Hbar - 1 # 2 historical regimes, 6 quarters of zlb, and a switch to uncertain flexible ait in the last regime
    # Add parameter values for additional regimes
    for i in 3:(n_tempZLB_regimes+1)
        set_regime_val!(m[:σ_g], i, m10[:σ_g].value)
        set_regime_val!(m[:σ_b], i, m10[:σ_b].value)
        set_regime_val!(m[:σ_μ], i, m10[:σ_μ].value)
        set_regime_val!(m[:σ_ztil], i, m10[:σ_ztil].value)
        set_regime_val!(m[:σ_λ_f], i, m10[:σ_λ_f].value)
        set_regime_val!(m[:σ_λ_w], i, m10[:σ_λ_w].value)
        set_regime_val!(m[:σ_r_m], i, m10[:σ_r_m].value)
        set_regime_val!(m[:σ_σ_ω], i, m10[:σ_σ_ω].value)
        set_regime_val!(m[:σ_μ_e], i, m10[:σ_μ_e].value)
        set_regime_val!(m[:σ_γ], i, m10[:σ_γ].value)
        set_regime_val!(m[:σ_π_star], i, m10[:σ_π_star].value)
        set_regime_val!(m[:σ_lr], i, m10[:σ_lr].value)
        set_regime_val!(m[:σ_z_p], i, m10[:σ_z_p].value)
        set_regime_val!(m[:σ_tfp], i, m10[:σ_tfp].value)
        set_regime_val!(m[:σ_gdpdef], i, m10[:σ_gdpdef].value)
        set_regime_val!(m[:σ_corepce], i, m10[:σ_corepce].value)
        set_regime_val!(m[:σ_gdp], i, m10[:σ_gdp].value)
        set_regime_val!(m[:σ_gdi], i, m10[:σ_gdi].value)

        for j = 1:DSGE.n_mon_anticipated_shocks(m)
            set_regime_val!(m[Symbol("σ_r_m$(j)")], i, m10[Symbol("σ_r_m$(j)")].value)
        end
    end

    # then, set up new set of regime dates to implement a temporary zlb and uncertain policy change
    # Set up regime dates
    temp_regime_dates = deepcopy(get_setting(m, :regime_dates))
    temp_regime_dates[1] = date_presample_start(m)
    end_date = Date(2020, 12, 31) + Dates.Year(1) + Dates.Month(9)
    for (i, date) in zip(3:(n_tempZLB_regimes+1),
                         Date(2020, 12, 31):Dates.Month(3):end_date) # add forecast dates
        temp_regime_dates[i] = date
    end
    m <= Setting(:regime_dates, temp_regime_dates)
    m <= Setting(:regime_switching, true)
    setup_regime_switching_inds!(m)

    # Following two settings turn imperfect awareness
    m <= Setting(:uncertain_altpolicy, true) # turns imperfect awareness on
    m <= Setting(:uncertain_temporary_altpolicy, true) # needed to allow a temporary alternative policy with imperfect awareness

    # set the "alternative" rule for imperfect awareness
    # taylor_rule corresponds to the default monetary policy rule we use,
    # which has coefficients on the output gap, inflation gap, and output growth
    m <= Setting(:alternative_policies, AltPolicy[taylor_rule()])

    # Now set up settings for temp alt policy
    m <= Setting(:gensys2, true) # Temporary alternative policies use a special gensys algorithm
    m <= Setting(:replace_eqcond, true) # The gensys2 algo replaces eqcond matrices, so this step is required
    if endo_zlb
        m <= Setting(:temporary_altpolicy_length, get_setting(m, :min_temporary_altpolicy_length)) # specifies number of ZLB regimes; required if there is
    else
        m <= Setting(:temporary_altpolicy_length, n_tempZLB_regimes - 2) # specifies number of ZLB regimes; required if there is
    end
    # further regime-switching after the ZLB ends (aside from
    # the extra regime for the "lift-off" from ZLB)
    credvec = range(0., stop = 1., length = n_tempZLB_regimes - 2) # Credibility of ZLB increases as time goes on
    replace_eqcond = Dict{Int, EqcondEntry}() # Which eqcond to use in which periods
    for (i, reg) in enumerate(3:n_tempZLB_regimes)
        #if using endozlb, only use zlb_rule in the eqcond for the minimum ZLB
        # length (in this example, two quarters)
        if endo_zlb && reg > 3+ get_setting(m, :min_temporary_altpolicy_length)
            reg_policy = DSGE.smooth_ait_gdp_alt()
        else
            reg_policy = DSGE.zlb_rule()
        end
        replace_eqcond[reg] = EqcondEntry(reg_policy, [credvec[i], 1. - credvec[i]]) # Temp ZLB rule in these regimes
    end
    replace_eqcond[n_tempZLB_regimes + 1] = EqcondEntry(DSGE.smooth_ait_gdp_alt(), [1., 0.]) # switch to AIT in the final regime
    m <= Setting(:regime_eqcond_info, replace_eqcond) # Add mapping of regimes to new eqcond matrices
    m <= Setting(:temporary_altpolicy_names, [:zlb_rule])

    # set up information sets
    m <= Setting(:tvis_information_set, vcat([1:1, 2:2],
                                             [i:get_setting(m, :n_regimes) for i in
                                              3:get_setting(m, :n_regimes)]))

    # pgap and ygap initialization
    df[!, :obs_pgap] .= NaN
    df[!, :obs_ygap] .= NaN
    ind_init = findfirst(df[!, :date] .== Date(2020, 6, 30))
    df[ind_init, :obs_pgap] = 0.125
    df[ind_init, :obs_ygap] = 12.

    # adjust nominal rates data during temporary ZLB
    tempzlb_quarters = DSGE.quarter_range(Date(2020, 12, 31),end_date)
    start_ind = findfirst(df[!, :date] .== Date(2020, 12, 31))
    if !isnothing(start_ind)
        inds_tempzlb = start_ind:findfirst(date_forecast_start(m))
        df[inds_tempzlb, :obs_nominalrate] .= NaN
        df[inds_tempzlb, [Symbol("obs_nominalrate$i") for i in 1:n_mon_anticipated_shocks(m)]] .= NaN
    end

    # Run forecast
    if endo_zlb
        # helper function used to set up parameter values when adjusting the length of the zlb
        # in the forecast
        function baseline_set_regime_vals(m::AbstractDSGEModel, n::Int; start_regime::Int = 3)
            if n >= start_regime
                for p in m.parameters
                    if !isempty(p.regimes)
                        if haskey(p.regimes, :value)
                            for i in start_regime:n
                                if !haskey(p.regimes[:value], i)
                                    ModelConstructors.set_regime_val!(p, i, ModelConstructors.regime_val(p, 1))
                                end
                            end
                        end
                        if haskey(p.regimes, :prior)
                            for i in start_regime:n
                                if !haskey(p.regimes[:prior], i)
                                    ModelConstructors.set_regime_prior!(p, i, ModelConstructors.regime_prior(p, 1))
                                end
                            end
                        end
                    end
                end
            end

            m
        end

        # Run the endogenous ZLB forecast
        fcast_tempzlb = DSGE.forecast_one_draw(m, :mode, :none, output_vars, modal_params,
                                               df; regime_switching = true, n_regimes = get_setting(m, :n_regimes), zlb_method = :temporary_altpolicy, rerun_smoother = true,
                                               set_regime_vals_altpolicy = (x, n) ->
                                               baseline_set_regime_vals(x, n;
                                                                              start_regime = 3))
    else
        fcast_tempzlb = DSGE.forecast_one_draw(m, :mode, :none, output_vars, modal_params,
                                               df; regime_switching = true, n_regimes = get_setting(m, :n_regimes))
    end
end

if make_plots
    plots_dict = Dict()

    # Define some plotting functions
    function make_comp_plot(m::AbstractDSGEModel, plot_mat::OrderedDict{String, Dict{Symbol, Array{Float64}}},
                            obs::Symbol, product::Symbol,
                            dates_plot, data_plot,
                            fcast_date, nfcast;
                            styles::Vector{Symbol} = repeat([:solid], length(keys(plot_mat))),
                            colors::Vector{Symbol} = repeat([:red], length(keys(plot_mat))))
        series_ind = if product ==:forecastobs
            m.observables[obs]
        else
            m.pseudo_observables[obs]
        end

        if product == :forecastobs
            dates_plot = data_plot[!, :date]
            histobs_plot = obs == :obs_hours ? data_plot[!, obs] : 4 .* data_plot[!, obs]
            fcast_date_tmp = vcat(dates_plot[end - 1], fcast_date)
            dates_plot = map(x -> DSGE.prev_quarter(x), dates_plot)
            p = Plots.plot(dates_plot, histobs_plot, color = :black,
                           label = "Data", linewidth = 3, legend = :bottomright, left_margin = 20px)
        else
            fcast_date_tmp = fcast_date
            p = Plots.plot()
        end

        for (i, (key, val)) in enumerate(plot_mat)

            if obs == :obs_hours || product == :forecastpseudo
                forecast_plot = val[product][series_ind, :][1:nfcast]
            else
                forecast_plot = 4 .* val[product][series_ind, :][1:nfcast]
            end
            if product == :forecastobs
                forecast_plot = vcat(histobs_plot[end], forecast_plot)
            end

            plot!(fcast_date_tmp, forecast_plot, color = colors[i],
                  label = key, linewidth = 1, linestyle = styles[i])
            if obs == :obs_nominalrate
                ylims!((-2, 4))
            end
        end
        return p
    end

    data_plot  = df[df[:, :date] .>= Date(2016, 3, 31), :]
    fcast_date = map(x -> DSGE.prev_quarter(x), DSGE.get_quarter_ends(date_forecast_start(m), DSGE.quartertodate("2031-Q4")))
    nfcast     = length(fcast_date)
    dates_plot = map(x -> DSGE.prev_quarter(x), data_plot[!, :date])

    for obs in [:obs_gdp, :obs_corepce, :obs_gdpdeflator, :obs_nominalrate, :obs_investment, :obs_consumption,
                :obs_spread, :obs_wages, :obs_hours]
        plots_dict[obs] = make_comp_plot(m, OrderedDict("$(pol_str)" => fcast_tempzlb),
                                         obs, :forecastobs,
                                         dates_plot, data_plot, fcast_date, nfcast,
                                         colors = [:red, :blue, :purple])
    end
end

if save_plots
    if !isdir(figurespath(m, "forecast"))
        mkdir(figurespath(m, "forecast"))
    end

    for (k, p) in plots_dict
        Plots.savefig(p, joinpath(figurespath(m, "forecast"), "$(pol_str)_$(string(k)).pdf"))
    end
end

nothing


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD

function InitialNetwork(dim, P)
    Network = sample([-1,0], Weights([1-P, P]),(dim,dim))
    return Network
end

function FindNeighbors(i, j, Network, dim)
    Numbereds = 0
    colorList = []
    for neighbor in ([0,-1], [-1,0])
        if j + neighbor[2] >= 1 && i + neighbor[1] >= 1 && Network[([i, j] + neighbor)...] > 0
            push!(colorList,Network[([i, j] + neighbor)...])
            Numbereds += 1
        end
    end
    return colorList, Numbereds #colorList contain the node and its neighbors
end

function RecursionMapping(x, L)
    while L[x] != L[L[x]]
        x = L[x]
    end
    return L[x]
end

function CheckPercolation(dim, L, Network)
    fside = []
    lside = []
    for i in 1:dim
        if Network[i,dim] != -1
            push!(lside, RecursionMapping(Network[i,dim],L))
        end
        if Network[i,1] != -1
            push!(fside, RecursionMapping(Network[i,1],L))
        end
    end
    RandomPoint = rand(Network)
    if RandomPoint != -1 && RecursionMapping(RandomPoint,L) in intersect(fside,lside)
        return 1
    else
        return 0
    end
end

function HKNetworkDynamic(dim, P)
    Network = InitialNetwork(dim, P)
    S = []
    L = []
    k = 1

    for j in 1:dim
        for i in 1:dim
            if Network[i,j] != 0
                continue
            end
            NColor, NNum = FindNeighbors(i,j,Network,dim)
            if NNum == 0
                Network[i,j] = k
                push!(L, k)
                push!(S, 1)
                k += 1
            elseif NNum == 1
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            elseif NNum == 2 && RecursionMapping(NColor[1],L) == RecursionMapping(NColor[2],L)
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            else
                S[RecursionMapping(NColor[1],L)] += 1 + S[RecursionMapping(NColor[2],L)]
                S[RecursionMapping(NColor[2],L)] = 0
                Network[i,j] = RecursionMapping(NColor[1],L)
                L[NColor[2]] = RecursionMapping(NColor[1],L)
            end
        end
    end
    # return Network
    return CheckPercolation(dim, L, Network)
end

HKNetworkDynamic(200, 0.6)

Plist = hcat(0.5:0.01:1)
Avglist = []
STDlist = []
TotData = []

for P in Plist
    print("\r$P")
    AllRuns = []
    for i in 1:10000
        push!(AllRuns, HKNetworkDynamic(200, P))
    end
    push!(TotData, AllRuns)
    push!(STDlist, std(AllRuns))
    push!(Avglist, mean(AllRuns))
end

save("200x200-HK.jld", "data", TotData)

scatter(Plist,
    Avglist,
    yerr = STDlist,
    legend = nothing,
    xlabel = L"P",
    ylabel = L"Avg\ Q_{\infty}",
    title = L"Avg \ Q_{\infty} \ for \ 200\times200 \ Network\ (10000 \ runs)")
savefig("../../Figs/Q6/200x200-10000.pdf")


text:
module day01

include("inputs.jl")
ex = inputs.ex_d1()
input = inputs.input_d1()

# Part 1
function sum_two_to_2020(input)
    for a in 1:length(input)-1, b in 1:a
        input[a] + input[b] == 2020 && return input[a] * input[b]
    end
end

@assert sum_two_to_2020(ex) == 514579
@info "Day 1, Part 1 answer: $(sum_two_to_2020(input))"

# Part 2
function sum_three_to_2020(input)
    local i = 0
    for a in 1:length(input)-1, b in 1:(a-1), c in 1:b
        input[a] + input[b] + input[c] == 2020 && return input[a]*input[b]*input[c]
    end
end

@assert sum_three_to_2020(ex) == 241861950
@info "Day 1, Part 2 answer: $(sum_three_to_2020(input))"

end


text:
using Revise, ProfileView
using nMDS, Distances
using Test
x = rand(10,10)
dist = pairwise(BrayCurtis(),x)
@profview nmds(dist,3)
@profview nmds(dist,3)

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


text:
function evolve(parent::String, target::String, mutrate::Float64, nchild::Int)
    fitness(a, b::String=target) = count(l == t for (l, t) in zip(a, b))
    function mutate(str::String, rate::Float64=mutrate)
        L = split(" ABCDEFGHIJKLMNOPQRSTUVWXYZ", "")
        r = collect(rand() < rate ? rand(L) : c for c in str)
        return string(r...)
    end
    println("Initial parent is $parent, its fitness is $(fitness(parent))")
    gens = 0
    while parent != target
        children = collect(mutate(parent, mutrate) for i in 1:nchild)
        bestfit, best = findmax(fitness.(children))
        parent = children[best]
        gens += 1
        if gens % 10 == 0
            println("After $gens generations, the new parent is $parent and its fitness is $(fitness(parent))")
        end
    end
    println("After $gens generations, the parent evolved into the target $target")
end

evolve("IU RFSGJABGOLYWF XSMFXNIABKT", "METHINKS IT IS LIKE A WEASEL", 0.08998, 100)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
abstract type AbstractCallback end

# source: https://youtu.be/EpNeNCGmyZE?t=633
# default fucntion blabla(cb;kwargs...) and fucntion blabla!(cb;kwargs...)
#for sym in [:before_fit, :before_epoch, :before_all_batches, :before_batch, :after_batch, :before_validate, :after_all_batches, :after_epoch, :after_fit]
#    @eval function $sym(cb::AbstractCallback;kwargs...)
#    end
#end

for sym in [:before_fit, :before_epoch, :before_all_batches, :before_batch, :after_batch, :before_validate, :after_all_batches, :after_epoch, :after_fit]
    @eval function $(Symbol(string(sym,"!")))(cb::AbstractCallback;kwargs...)
    end
end

@with_kw mutable struct TrainEvalCallback <: AbstractCallback
    order::Int = -1
    n_epochs::Int = 0
    current_epoch::Float64 = 0.0
    n_batches::Int = 0
    current_batch::Int = 0
    in_train::Bool = true
    use_cuda::Bool = false
    num_recorder::Int = 0
end

init_c_epoch!(cb::TrainEvalCallback) = cb.current_epoch = 0
init_c_batch!(cb::TrainEvalCallback) = cb.current_batch = 0
incr_c_batch!(cb::TrainEvalCallback) = cb.current_batch += 1
incr_c_epoch!(cb::TrainEvalCallback) = cb.current_epoch += 1. /cb.n_batches

function calculate_n_nbatches!(cb::TrainEvalCallback, d::DataLoader)
    cb.n_batches =  ceil(d.nobs / d.batchsize)
end

function before_fit!(cb::TrainEvalCallback; n_epochs::Int, model, loss_func, kwargs...)
    init_c_epoch!(cb)
    init_c_batch!(cb)
    cb.n_epochs = n_epochs
    if CUDA.has_cuda_gpu()
        cb.use_cuda = true
    end
end

function after_batch!(cb::TrainEvalCallback; kwargs...)
    incr_c_batch!(cb)
    incr_c_epoch!(cb)
end

function before_epoch!(cb::TrainEvalCallback; epoch, kwargs...)
    cb.current_epoch = epoch-1# start at zero
    cb.in_train = true
end

function before_all_batches!(cb::TrainEvalCallback; d::DataLoader, kwargs...)
    calculate_n_nbatches!(cb,d)
end

function before_validate!(cb::TrainEvalCallback; kwargs...)
    cb.in_train = false
end

@with_kw mutable struct AvgStatsCallback{F<:Function} <: AbstractCallback
    order::Int = 0
    current_sum_stat::Float64 = 1000.0
    num_samples::Int = 0
    average_stat::Float64 = 1000.0
    f::F = logitcrossentropy
    name::String = "loss"
    recorder::Bool = true
    rec_ind::Int = 0
    rec_num::Int = 0
    recorder_initialzed::Bool = false
end

function AvgStatsCallback(order::Int,f::Function; kwargs...)
    AvgStatsCallback{typeof(f)}(order=order, f=f, name=String(Symbol(f)))
end

function before_fit!(cb::AvgStatsCallback; te_cb::TrainEvalCallback, kwargs...)
    if cb.recorder && !cb.recorder_initialzed
        te_cb.num_recorder += 1
        cb.rec_num = te_cb.num_recorder
        cb.rec_ind = 0
        cb.recorder_initialzed = true
    elseif cb.recorder_initialzed && cb.recorder_initialzed
        cb.rec_ind = 0
    end
end

function after_batch!(cb::AvgStatsCallback; rec, batch_loss, batch_size, batch_pred, batch_label, kwargs...)
    if cb.name == "loss"
         incr = batch_loss
    else
        incr = cb.f(batch_pred,batch_label)
    end
    cb.current_sum_stat += incr*batch_size
    cb.num_samples += batch_size
    if cb.recorder
        cb.rec_ind += 1
        rec[cb.rec_num,cb.rec_ind] = incr
    end
end

function after_all_batches!(cb::AvgStatsCallback; te_cb::TrainEvalCallback, epoch::Int, kwargs...)
    cb.average_stat = cb.current_sum_stat/cb.num_samples
    cb.current_sum_stat = 0
    cb.num_samples = 0
    if te_cb.in_train
        @assert round(Int, te_cb.current_epoch) == epoch "pb not passing over all the train set"
        println("Train - epoch: ", epoch, " ", cb.name, ": ", cb.average_stat)
    else
        @assert round(Int, te_cb.current_epoch) == epoch+1 "pb not passing over all the validation set"
        println("Val - epoch: ", epoch, " ", cb.name, ": ", cb.average_stat)
    end
end


text:
module TetrahedronIntegration

include("tetrahedron.jl")
include("gaussian.jl")
include("sampling.jl")

end

max_stars_repo_name: jaemolihm/TetrahedronIntegration

max_issues_repo_name: jaemolihm/TetrahedronIntegration

max_forks_repo_name: jaemolihm/TetrahedronIntegration


text:
using SearchLight, SearchLight.QueryBuilder

try
  SearchLight.Configuration.load() |> SearchLight.Database.connect!
  SearchLight.load_resources()
catch ex
  @error "Failed loading SearchLight database configuration. Please make sure you have a valid connection.yml file."
end


text:
#=
Hidden layer activation functions and gradients.
Options are:
:logistic, :tanh, :rectifier
=#

################################################################################
#
# Activation Functions and Derivatives
#
################################################################################

########################################
# Logistic
########################################

# logistic function of element z
# f(z) = 1 / (1 + exp(-z))
function logistic{F<:AbstractFloat}(z::F)
    one(z) / (one(z) + exp(-z))
end

# logistic function of array Z
# replacing array A in-place
# A[i] = f(Z[i]) = 1 / (1 + exp(-Z[i])) ∀ i
function logistic!{F<:AbstractFloat}(A::Array{F}, Z::Array{F})
    size(A) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    @threads for i in eachindex(Z)
        @inbounds A[i] = logistic(Z[i])
    end
end

# logistic function derivative of element z
# df/dz = f(z) (1 - f(z))
function logistic_deriv{F<:AbstractFloat}(z::F)
    s = logistic(z)
    s * (F(1) - s)
end

# logistic function derivative of array Z
# ∂f/∂Z[i] = f(Z[i]) (1 - f(Z[i]))
# replacing dest in-place with Hademaard product of dest and derivative i.e.
# dest[i] = dest[i] * ∂f/∂Z[i]
function logistic_deriv!{F<:AbstractFloat}(dest::Array{F}, Z::Array{F})
    size(dest) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    @threads for i in eachindex(Z)
        @inbounds dest[i] *= logistic_deriv(Z[i])
    end
end

########################################
# Tanh
########################################

# Elementwise tanh function of array Z
# replacing array A in-place
function tanh!{F<:AbstractFloat}(A::Array{F}, Z::Array{F})
    size(A) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    @threads for i in eachindex(Z)
        @inbounds A[i] = tanh(Z[i])
    end
end

# Elementwise tanh function derivative of array Z
# ∂f/∂Z[i] = 1 - f(Z[i])^2
# replacing dest in-place with Hademaard product of dest and derivative i.e.
# dest[i] = dest[i] * ∂f/∂Z[i]
function tanh_deriv!{F<:AbstractFloat}(dest::Array{F}, Z::Array{F})
    size(dest) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    @threads for i in eachindex(Z)
        @inbounds dest[i] *= (F(1) - tanh(Z[i])^2)
    end
end

########################################
# Rectified linear units
########################################

# Elementwise rectifier function of array Z
# replacing array A in-place
function relu!{F<:AbstractFloat}(A::Array{F}, Z::Array{F})
    size(A) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    @threads for i in eachindex(Z)
        @inbounds A[i] = max(F(0), Z[i])
    end
end

# Elementwise rectifier function derivative of array Z
# ∂f/∂Z[i] = 0 (if Z[i] <= 0) or 1 (if Z[i] > 0)
# replacing dest in-place with Hademaard product of dest and derivative i.e.
# dest[i] = dest[i] * ∂f/∂Z[i]
function relu_deriv!{F<:AbstractFloat}(dest::Array{F}, Z::Array{F})
    size(dest) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    for i in eachindex(Z)
        @inbounds (Z[i] <= F(0)) && (dest[i] = F(0))
        #@inbounds dest[i] *= ifelse(Z[i] > F(0), F(1), F(0))
    end
end

########################################
# Softmax
########################################

# Softmax function of matrix Z replacing matrix A in-place
# Applies function to each column of Z i.e. each column will sum to 1
function softmax!{F<:AbstractFloat}(A::Matrix{F}, Z::Matrix{F}, u::Matrix{F})
    #size(A) == size(Z) || throw(DimensionMismatch("Inconsistent array sizes"))
    maximum!(u, Z)
    for j in indices(Z, 2)
        s = F(0)
        for i in indices(Z, 1)
            @inbounds s += (A[i, j] = exp(Z[i, j] - u[j]))
        end
        for i in indices(Z, 1)
            @inbounds A[i, j] /= s
        end
    end
end

################################################################################
#
# Activation Dictionaries
#
################################################################################

activation_dict = Dict{Symbol, Function}(:logistic => logistic!,
                                         :tanh => tanh!,
                                         :relu => relu!)

activation_deriv_dict = Dict{Symbol, Function}(:logistic => logistic_deriv!,
                                               :tanh => tanh_deriv!,
                                               :relu => relu_deriv!)


text:
using DistributedReinforcementLearning
using Test
using Distributed
using ReinforcementLearningBase
using ReinforcementLearningCore
using ReinforcementLearningZoo
using ReinforcementLearningEnvironments
using Flux

@testset "DistributedReinforcementLearning.jl" begin

include("actor.jl")
include("core.jl")

end


text:
if length(ARGS) >= 1
  nw = parse(ARGS[1])
  if nw >1
     addprocs(nw)
  end
else
  addprocs()
end

import Distributions; @everywhere using Distributions
import ABC; @everywhere using ABC

@everywhere include(joinpath(Pkg.dir("ABC"),"src/composite.jl"))
@everywhere using CompositeDistributions
@everywhere import Compat.view

# Set Prior for Population Parameters
@everywhere d1 = Rayleigh(0.3)
@everywhere d2 = MvNormal([0.0, 1.0], ones(2))
@everywhere param_prior = CompositeDist( ContinuousDistribution[d1,d2] )

# Code to generate simulated data given array of model parameters
@everywhere num_data_default = 100
@everywhere num_outputs = 2
@everywhere function gen_data(theta::Array, n::Integer = num_data_default)
  data = Array(Float64,(3,n))
  data[1,:] = rand(Rayleigh(theta[1]),n)
  data[2:3,:] = rand(MvNormal(theta[2:3],ones(2)),n)
  return data
end

# Function to adjust originally proposed model parameters, so that they will be valid
@everywhere function normalize_theta13_pos!(theta::Array)
 theta[1] = abs(theta[1])
 theta[3] = abs(theta[3])
 return theta
end

# Function to test if the proposed model parameter are valid
@everywhere function is_valid_theta13_pos(theta::Array)
  theta[1]>0.0  && theta[3]>0.0
end

# True Population Parameters
@everywhere theta_true = [0.3, 0.0, 1.0]

# Tell ABC what it needs to know for a simulation
@everywhere in_parallel = length(workers()) > 1 ? true : false
@everywhere abc_plan = abc_pmc_plan_type(gen_data,ABC.calc_summary_stats_mean_var,ABC.calc_dist_max, param_prior, is_valid=is_valid_theta13_pos,num_max_attempt=10000, in_parallel=in_parallel);

# Generate "true/observed "data"
data_true = gen_data(theta_true)   # Draw "real" data from same model as for analysis
ss_true = abc_plan.calc_summary_stats(data_true)
#println("theta= ",theta_true," ss= ",ss_true, " d= ", 0.)

# Run ABC simulation
@time pop_out = run_abc(abc_plan,ss_true,in_parallel=in_parallel);

#=
# Optional plotting of results
using PyPlot

hist(pop_out.weights*length(pop_out.weights));
hist(pop_out.dist);

num_grid_x = 100
num_grid_y = 100
limit = 1.0
x = collect(linspace(theta_true[1]-limit,theta_true[1]+limit,num_grid_x));
y = collect(linspace(theta_true[2]-limit,theta_true[2]+limit,num_grid_y));
z = zeros(Float64,(num_param,length(x),length(y)))
for i in 1:length(x), j in 1:length(y)
    z[1,i,j] = x[i]
    z[2,i,j] = y[j]
end
z = reshape(z,(num_param,length(x)*length(y)))
zz = [ ABC.pdf(ABC.GaussianMixtureModelCommonCovar(pop_out.theta,pop_out.weights,ABC.cov_weighted(pop_out.theta',pop_out.weights)),vec(z[:,i])) for i in 1:size(z,2) ]
zz = reshape(zz ,(length(x),length(y)));
levels = [exp(-0.5*i^2)/sqrt(2pi^num_param) for i in 0:5];
PyPlot.contour(x,y,zz',levels);
plot(pop_out.theta[1,:],pop_out.theta[2,:],".");
=#


text:
plusnumber(t::Int) = t>0 ? "+$t" : "$t"

#treat thsi table as a one-off since it has no common elements w/ other tables
function maketable7(bb::DataFrame)
  T7_PANEL_A_COLKEYS = [:StdDev, :Mean, :Range, :Return, :AbsReturn]
  T7_PANEL_A_NAME::String = "t7panela-$(REPLICATION_TYPE[])"
  T7_PANEL_A_KEYROW_RANGES::Dict = Dict(-3=>-6:-3, 0=>-2:2, 3=>3:6)
  local days::Int = SEO_DAYS

  table7panela(bb,
      days=days,
      tablename= T7_PANEL_A_NAME,
      colkeys = T7_PANEL_A_COLKEYS,
      keyrowranges = T7_PANEL_A_KEYROW_RANGES,
      )
end

function table7panela(bb::DataFrame;
    days = SEO_DAYS,
    tablename::String = error("tablename is required"),
    tablepath::String = TABLE_PATH,
    decimals::Int = DEFAULT_DECIMALS,
    colkeys::Vector{Symbol} = error("colkeys is required"),
    colnames::Vector{Vector{String}} = [["Daily Cross-", "", ""],
      ["Sectional Returns", "", "Average"],
      (string).(colkeys)],
    widthcolnames::Vector{Vector{Int}} = [[2,1,2], [2,1,2], ones(Int, length(colnames[2]))],
    descrownames::Vector{String} = (plusnumber).(-days:days),
    keyrowranges::Dict = error("keyrow ranges is required"))

  local parallel::Bool = PARALLEL[] && (!DISABLE_FM_PARALLEL)

  #################################
  #basic lookups
  i2t::Vector{Int} = collect(-days:days)
  dayfields::Vector{Symbol} = getdayfields(days)
  Fₜ::Dict{Int, Symbol} = Dict(t=>dayfields[i] for (i,t) ∈ enumerate(i2t))

  #preallocation and housekeeping

  Ncols::Int = sum(widthcolnames[1])
  Nrows::Int = 2*days+1
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, Ncols)).(1:Nrows)

  #create an index to make this part of the code more readable
  colidx::Dict = Dict(colkey=>i for (i, colkey) ∈ enumerate(colkeys))

  #now populate the the table
  for (i,rowkey) ∈ enumerate(descrownames)
    f::Symbol = dayfields[i]
    t::Int = i2t[i] #get the period offset

    #compute the following stats for each row
    daystd::Float64 = std(skipmissing(bb[!, f]))*100
    daymean::Float64 = mean(skipmissing(bb[!, f]))*100

    desccontent[i][colidx[:StdDev]] = num2str(daystd, decimals)
    desccontent[i][colidx[:Mean]] = num2str(daymean, decimals)

    if haskey(keyrowranges, t)

      tₗ::Int = minimum(keyrowranges[t])
      tₕ::Int = maximum(keyrowranges[t])
      rngfields::Vector{Symbol} = (i->Fₜ[i]).(keyrowranges[t])
      desccontent[i][colidx[:Range]] = "\\text{$(plusnumber(tₗ)) to $(plusnumber(tₕ))}"

      println(rngfields)
      rngmat::Matrix{Float64} = Matrix(bb[:, rngfields])
      rets::Vector{Float64} = vec(prod(rngmat .+ 1., dims=2)) .- 1
      retmean::Float64 = mean(rets)*100
      retabsmean::Float64 = mean((abs).(rets))*100

      desccontent[i][colidx[:Return]] = num2str(retmean, decimals)
      desccontent[i][colidx[:AbsReturn]] = num2str(retabsmean, decimals)
    else
      desccontent[i][colidx[:Range]] = ""
      desccontent[i][colidx[:Return]] = ""
      desccontent[i][colidx[:AbsReturn]] = ""
    end
  end

  alignmentstring::String = "rrr c rr"

  #println(desccontent)

  #display(desccontent)
  #construct the tex table
  tt::String = textable(colnames=colnames,
    descrownames=descrownames,
    desccontent=desccontent,
    alignmentstring=alignmentstring,
    widthcolnames=widthcolnames)

  writetextable(tt, path=tablepath, outname="$tablename.tex")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Symbolics, Test, SparseArrays
@variables a b c

# Auxiliary Functions
get_sparsity_pattern(h::Union{SparseVector{Num}, SparseMatrixCSC{Num,Int}}) = get_sparsity_pattern(Array(h))
get_sparsity_pattern(h::Array{Num}) = sparse(Int.(.!isequal.(h, 0)))

input = [1, 2, 3]

# ===== Dense tests =====
# Arrays of Matrices
h_dense_arraymat = [[a 1; b 0], [0 0; 0 0], [a c; 1 0]] # empty array support required
function h_dense_arraymat_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1] 1; b[1] 0]
    out[2] .= [0 0; 0 0]
    out[3] .= [a[1] c[1]; 1 0]
end

h_dense_arraymat_ip! = eval(Symbolics.build_function(h_dense_arraymat, [a, b, c])[2])
h_dense_arraymat_ip_skip! = eval(Symbolics.build_function(h_dense_arraymat, [a, b, c], skipzeros=true, fillzeros=false)[2])
out_1_arraymat = [fill(42, 2, 2) for i in 1:3]
out_2_arraymat = deepcopy(out_1_arraymat)
h_dense_arraymat_julia!(out_1_arraymat, input)
h_dense_arraymat_ip_skip!(out_2_arraymat, input)
@test all(isequal(42), out_2_arraymat[2])
foreach(mat->fill!(mat, 0), out_2_arraymat)
h_dense_arraymat_ip_skip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat
h_dense_arraymat_ip!(out_2_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Arrays of Matrices, Heterogeneous element types
test_exp = [exp(a) * exp(b), a]
h_dense_arraymat_het = [Symbolics.hessian(t, [a, b]) for t in test_exp]
function h_dense_arraymat_het_julia!(out, x)
    a, b, c = x
    out[1] .= [exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1]); exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1])]
    out[2] .= [0 0; 0 0]
end

h_dense_arraymat_het_str = Symbolics.build_function(h_dense_arraymat_het, [a, b, c])
h_dense_arraymat_het_ip! = eval(h_dense_arraymat_het_str[2])
out_1_arraymat_het = [Array{Float64}(undef, 2, 2) for i in 1:2]
out_2_arraymat_het = [similar(x) for x in out_1_arraymat_het]
h_dense_arraymat_het_julia!(out_1_arraymat_het, input)
h_dense_arraymat_het_ip!(out_2_arraymat_het, input)
@test out_1_arraymat_het == out_2_arraymat_het

# Arrays of 1D Vectors
h_dense_arrayvec = [[a, 0, c], [0, 0, 0], [1, a, b]] # same for empty vectors, etc.
function h_dense_arrayvec_julia!(out, x)
    a, b, c = x
    out[1] .= [a[1], 0, c[1]]
    out[2] .= [0, 0, 0]
    out[3] .= [1, a[1], b[1]]
end

h_dense_arrayvec_str = Symbolics.build_function(h_dense_arrayvec, [a, b, c])
h_dense_arrayvec_ip! = eval(h_dense_arrayvec_str[2])
out_1_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
out_2_arrayvec = [Vector{Int64}(undef, 3) for i in 1:3]
h_dense_arrayvec_julia!(out_1_arrayvec, input)
h_dense_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_dense_arrayNestedMat = [[[a 1; b 0], [0 0; 0 0]], [[b 1; a 0], [b c; 0 1]]]
function h_dense_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1] .= [a[1] 1; b[1] 0]
    out[1][2] .= [0 0; 0 0]
    out[2][1] .= [b[1] 1; a[1] 0]
    out[2][2] .= [b[1] c[1]; 0 1]
end

h_dense_arrayNestedMat_str = Symbolics.build_function(h_dense_arrayNestedMat, [a, b, c])
h_dense_arrayNestedMat_ip! = eval(h_dense_arrayNestedMat_str[2])
out_1_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]] # avoid undef broadcasting issue
out_2_arrayNestedMat = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]]
h_dense_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_dense_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# Arrays of Arrays of Matrices, Heterogeneous element types
test_exp = [exp(a) * exp(b), a]
h_dense_arrayNestedMat_het = [[Symbolics.hessian(t, [a, b]) for t in test_exp], [Num[0 0; 0 0], Num[0 0; 0 0]]]
function h_dense_arrayNestedMat_het_julia!(out, x)
    a, b, c = x
    out[1][1] .= [exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1]); exp(a[1]) * exp(b[1]) exp(a[1]) * exp(b[1])]
    out[1][2] .= [0 0; 0 0]
    out[2][1] .= [0 0; 0 0]
    out[2][2] .= [0 0; 0 0]
end
h_dense_arrayNestedMat_het_str = Symbolics.build_function(h_dense_arrayNestedMat_het, [a, b, c])
h_dense_arrayNestedMat_het_ip! = eval(h_dense_arrayNestedMat_het_str[2])
out_1_arrayNestedMat_het = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]] # avoid undef broadcasting issue
out_2_arrayNestedMat_het = [[rand(Int64, 2, 2), rand(Int64, 2, 2)], [rand(Int64, 2, 2), rand(Int64, 2, 2)]]
h_dense_arrayNestedMat_julia!(out_1_arrayNestedMat_het, input)
h_dense_arrayNestedMat_ip!(out_2_arrayNestedMat_het, input)
@test out_1_arrayNestedMat_het == out_2_arrayNestedMat_het

# ===== Sparse tests =====
# Array of Matrices
h_sparse_arraymat = sparse.([[a 1; b 0], [0 0; 0 0], [a c; 1 0]])
function h_sparse_arraymat_julia!(out, x)
    a, b, c = x
    out[1][1, 1] = a[1]
    out[1][1, 2] = 1
    out[1][2, 1] = b[1]
    out[2] = sparse([0 0; 0 0]) # no undef constructor for SparseMatrixCSC
    out[3][1, 1] = a[1]
    out[3][1, 2] = c[1]
    out[3][2, 1] = 1
end

h_sparse_arraymat_str = Symbolics.build_function(h_sparse_arraymat, [a, b, c])
h_sparse_arraymat_ip! = eval(h_sparse_arraymat_str[2])
h_sparse_arraymat_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arraymat)
out_1_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns]
out_2_arraymat = [similar(h) for h in h_sparse_arraymat_sparsity_patterns] # can't do similar() because it will just be #undef, with the wrong sparsity pattern
h_sparse_arraymat_ip!(out_2_arraymat, input)
h_sparse_arraymat_julia!(out_1_arraymat, input)
@test out_1_arraymat == out_2_arraymat

# Array of 1D Vectors
h_sparse_arrayvec = sparse.([[a, 0, c], [0, 0, 0], [1, a, b]])
function h_sparse_arrayvec_julia!(out, x)
    a, b, c = x
    out[1][1] = a[1]
    out[1][3] = c[1]
    out[2] = sparse([0, 0, 0]) # necessary because sparsity pattern is 3 elements with 0 stored, not 0 elements
    out[3][1] = 1
    out[3][2] = a[1]
    out[3][3] = b[1]
end

h_sparse_arrayvec_str = Symbolics.build_function(h_sparse_arrayvec, [a, b, c])
h_sparse_arrayvec_ip! = eval(h_sparse_arrayvec_str[2])
h_sparse_arrayvec_sparsity_patterns = map(get_sparsity_pattern, h_sparse_arrayvec)
out_1_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
out_2_arrayvec = [similar(h) for h in h_sparse_arrayvec_sparsity_patterns]
h_sparse_arrayvec_julia!(out_1_arrayvec, input)
h_sparse_arrayvec_ip!(out_2_arrayvec, input)
@test out_1_arrayvec == out_2_arrayvec

# Arrays of Arrays of Matrices
h_sparse_arrayNestedMat = [sparse.([[a 1; b 0], [0 0; 0 0]]), sparse.([[b 1; a 0], [b c; 0 1]])]
function h_sparse_arrayNestedMat_julia!(out, x)
    a, b, c = x
    out[1][1][1, 1] = a[1]
    out[1][1][1, 2] = 1
    out[1][1][2, 1] = b[1]
    out[1][2] = sparse([0 0; 0 0])
    out[2][1][1, 1] = b[1]
    out[2][1][1, 2] = 1
    out[2][1][2, 1] = a[1]
    out[2][2][1, 1] = b[1]
    out[2][2][1, 2] = c[1]
    out[2][2][2, 2] = 1
end

h_sparse_arrayNestedMat_str = Symbolics.build_function(h_sparse_arrayNestedMat, [a, b, c])
h_sparse_arrayNestedMat_ip! = eval(h_sparse_arrayNestedMat_str[2])
h_sparse_arrayNestedMat_sparsity_patterns = [map(get_sparsity_pattern, h) for h in h_sparse_arrayNestedMat]
out_1_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
out_2_arrayNestedMat = [[similar(h_sub) for h_sub in h] for h in h_sparse_arrayNestedMat_sparsity_patterns]
h_sparse_arrayNestedMat_julia!(out_1_arrayNestedMat, input)
h_sparse_arrayNestedMat_ip!(out_2_arrayNestedMat, input)
@test out_1_arrayNestedMat == out_2_arrayNestedMat

# Additional Tests
# Returning 0-element structures (corresponding to empty Jacobians)
# Arrays of Matrices
h_empty = [[a b; c 0], Array{Num,2}(undef, 0,0)]
h_empty_str = Symbolics.build_function(h_empty, [a, b, c])
h_empty_ip! = eval(h_empty_str[2])
out = [Matrix{Int64}(undef, 2, 2), Matrix{Int64}(undef, 0, 0)]
h_empty_ip!(out, input) # should just not fail

# Array of Vectors
h_empty_vec = [[a, b, c, 0], Vector{Num}(undef,0)]
h_empty_vec_str = Symbolics.build_function(h_empty_vec, [a, b, c])
h_empty_vec_ip! = eval(h_empty_vec_str[2])
out = [Vector{Int64}(undef, 4), Vector{Int64}(undef, 0)]
h_empty_vec_ip!(out, input) # should just not fail

# Arrays of Arrays of Matrices
h_emptyNested = [[[a b; c 0]], Array{Array{Num, 2}}(undef, 0)] # emptyNested array of arrays
h_emptyNested_str = Symbolics.build_function(h_emptyNested, [a, b, c])
h_emptyNested_ip! = eval(h_emptyNested_str[2])
out = [[[1 2;3 4]], Array{Array{Int64,2},1}(undef, 0)]
h_emptyNested_ip!(out, input) # should just not fail


text:
"""
get detailed GDAL information

## return
- `file`     :
- `range`    : [lon_min, lon_max, lat_min, lat_max]
- `cellsize` : [dx, dy]
- `lon`      : longitudes with the length of nlon
- `lat`      : latitudes with the length of nlat
- `dim`      : [width, height, ntime]
"""
function gdalinfo(ga::GeoArray)
    # mid = [1, 1];
    # if length(mid) == 1; mid = [mid, mid]; end
    dx = ga.f.linear[1]
    dy = abs(ga.f.linear[4])
    dy2 = ga.f.linear[4]

    mid = [1, 1] # default mid is true
    delta = [dx, dy]/2 .* mid

    b = st_bbox(ga)
    lon = b.xmin + delta[1]:dx:b.xmax
    lat = b.ymin + delta[2]:dy:b.ymax
    if dy2 < 0; lat = reverse(lat); end

    Dict(
        "bbox"     => b,
        "cellsize" => [dx, dy],
        "coords"    => [lon, lat],
        "dim"      => size(ga.A))
end

function gdalinfo(file::AbstractString)
    ds = ArchGDAL.read(file)
    gt = ArchGDAL.getgeotransform(ds)
    # band = ArchGDAL.getband(ds, 1)
    w, h = ArchGDAL.width(ds), ArchGDAL.height(ds)
    dx, dy = gt[2], -gt[end]
    x0 = gt[1] #+ dx/2
    x1 = x0 + w* dx
    y1 = gt[4] #- dy/2
    y0 = y1 - h*dy
    b = bbox(x0, y0, x1, y1)

    lon = x0 + dx/2 : dx: x1
    lat = reverse(y0 + dy/2 : dy: y1)
    nband = ArchGDAL.nraster(ds)

    Dict(
        "file"     => basename(file),
        "bbox"     => b,
        "cellsize" => [dx, dy],
        "coords"    => [lon, lat],
        "dim"      => [w, h, nband])
end

gdal_open(file::AbstractString) = ArchGDAL.read(file)

function nband(file::AbstractString)
    # ArchGDAL.unsafe_read(file) do ds
    ArchGDAL.read(file) do ds
        ArchGDAL.nraster(ds)
    end
end
nratser = nband

export gdalinfo, nband, gdal_open


text:
using Raylib

const screenWidth = 800
const screenHeight = 450

function main()
    Raylib.InitWindow(
        screenWidth, screenHeight,
        "raylib [core] example - input mouse wheel"
    )

    boxPositionY = floor(Int, screenHeight/2) - 40
    scrollSpeed = 4 # Scrolling speed in pixels

    Raylib.SetTargetFPS(60)
    while !Raylib.WindowShouldClose()    # Detect window close button or ESC key
        # Update
        boxPositionY -= Int(Raylib.GetMouseWheelMove() * scrollSpeed)

        # Draw
        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE);
        Raylib.DrawRectangle(floor(Int, screenWidth/2) - 40, boxPositionY, 80, 80, Raylib.MAROON)
        Raylib.DrawText("Use mouse wheel to move the cube up and down!", 10, 10, 20, Raylib.GRAY)
        Raylib.DrawText("Box position Y: $boxPositionY", 10, 40, 20, Raylib.LIGHTGRAY)
        Raylib.EndDrawing()
    end

    Raylib.CloseWindow()
end


text:
module ParametricContractor

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOParametricInterval

# sets options for contractor
opt1 = Any[100    #Int64: Number of iterations
       1.0E-30 #Float64: Tolerance for equality of
       1.0E-6 #Float64: Add Interval(1E-8,1E8) to add to M[i,i] when
              #         processing extended interval division.
      ]

# Test Problem #1 Out of place
P1 = [Interval(5.0,7.0),Interval(5.0,7.0)]
Z1 = [Interval(-1.5, 0.0),Interval(0.0, 0.5)]
h1(z,p) = [z[1]^2+z[2]^2+p[1]*z[1]+4;
           z[1]+p[2]*z[2]]
hj1(z,p) = [(2*z[1]+p[1]) (2*z[2]);
              1              p[2]]
Eflag = false
Iflag = false
eDflag = false
newtonGS1 = PI_NewtonGS(Z1,P1,hj1,h1,opt1,Eflag,Iflag,eDflag)
krawczykCW1 = PI_KrawczykCW(Z1,P1,hj1,h1,opt1,Eflag,Iflag)

@test -1E-4 <= newtonGS1[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= newtonGS1[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= newtonGS1[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= newtonGS1[1][2].hi - 0.208486 <= 1E-4
@test -1E-4 <= krawczykCW1[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= krawczykCW1[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= krawczykCW1[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= krawczykCW1[1][2].hi - 0.208486 <= 1E-4

# Test Problem #1 In place
function h1!(h,z,p)
   h[1] = z[1]^2+z[2]^2+p[1]*z[1]+4
   h[2] = z[1]+p[2]*z[2]
end

function hj1!(h,z,p)
   h[1,1] = (2*z[1]+p[1])
   h[1,2] = 2*z[2]
   h[2,1] = 1
   h[2,2] = p[2]
end

Eflag = false
Iflag = false
eDflag = false
newtonGS2 = PIn_NewtonGS(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag,eDflag)
krawczykCW2 = PIn_KrawczykCW(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag)

@test -1E-4 <= newtonGS2[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= newtonGS2[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= newtonGS2[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= newtonGS2[1][2].hi - 0.208486 <= 1E-4
@test -1E-4 <= krawczykCW2[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= krawczykCW2[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= krawczykCW2[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= krawczykCW2[1][2].hi - 0.208486 <= 1E-4

# Starts testing problem # 2
P2 = [Interval(6.0,9.0)]
Z2 = [Interval(-1.78, -0.1)]

# Test Problem #2 Out of Place
h2(x,p) = [x[1]^2 + p[1]*x[1] + 4.0]
hj2(x,p) = [2*x[1]+p[1]]
Eflag = false
Iflag = false
eDflag = false
newtonGS3 = PI_NewtonGS(Z2,P2,hj2,h2,opt1,Eflag,Iflag,eDflag)
krawczykCW3 = PI_KrawczykCW(Z2,P2,hj2,h2,opt1,Eflag,Iflag)
@test -0.772413 - 1E-4 <= krawczykCW3[1][1].lo <= -0.772413 + 1E-4
@test -0.383299 - 1E-4 <= krawczykCW3[1][1].hi <= -0.383299 + 1E-4
@test krawczykCW3[2] == false
@test krawczykCW3[3] == true
@test -0.769225 - 1E-4 <= newtonGS3[1][1].lo <= -0.769225 + 1E-4
@test -0.461725 - 1E-4 <= newtonGS3[1][1].hi <= -0.461725 + 1E-4
@test newtonGS3[3] == false
@test newtonGS3[4] == true

# Test Problem #2 Inplace
function h2!(hout,x,p)
   hout[:] = [x[1]^2 + p[1]*x[1] + 4.0]
end
function hj2!(hout,x,p)
   hout[:] = [2*x[1]+p[1]]
end
Eflag = false
Iflag = false
eDflag = false
newtonGS4 = PIn_NewtonGS(Z2,P2,hj2!,h2!,opt1,Eflag,Iflag,eDflag)
krawczykCW4 = PIn_KrawczykCW(Z2,P2,hj2!,h2!,opt1,Eflag,Iflag)
@test -0.772413 - 1E-4 <= krawczykCW4[1][1].lo <= -0.772413 + 1E-4
@test -0.383299 - 1E-4 <= krawczykCW4[1][1].hi <= -0.383299 + 1E-4
@test krawczykCW4[2] == false
@test krawczykCW4[3] == true
@test -0.769225 - 1E-4 <= newtonGS4[1][1].lo <= -0.769225 + 1E-4
@test -0.461725 - 1E-4 <= newtonGS4[1][1].hi <= -0.461725 + 1E-4
@test newtonGS4[3] == false
@test newtonGS4[4] == true

# Test Problem #3 (Exclusion)
Z2ex = [Interval(-30, -25)]
Eflag = false
Iflag = false
eDflag = false
newtonGS5 = PI_NewtonGS(Z2ex,P2,hj2,h2,opt1,Eflag,Iflag,eDflag)
krawczykCW5 = PI_KrawczykCW(Z2ex,P2,hj2,h2,opt1,Eflag,Iflag)
@test newtonGS5[3] == true
@test krawczykCW5[2] == true

Eflag = false
Iflag = false
eDflag = false
newtonGS6 = PIn_NewtonGS(Z2ex,P2,hj2!,h2!,opt1,Eflag,Iflag,eDflag)
krawczykCW6 = PIn_KrawczykCW(Z2ex,P2,hj2!,h2!,opt1,Eflag,Iflag)
@test newtonGS6[3] == true
@test krawczykCW6[2] == true

# Test Problem #4 (Extended Division)
Z3 = [Interval(-10, 10)]
newtonGS7 = PI_NewtonGS(Z3,P2,hj2,h2,opt1,Eflag,Iflag,eDflag)
newtonGS8 = PIn_NewtonGS(Z3,P2,hj2!,h2!,opt1,Eflag,Iflag,eDflag)
@test 0.285713 - 1E-4 <= newtonGS7[1][1].lo <= 0.285713 + 1E-4
@test 10 - 1E-4 <= newtonGS7[1][1].hi <= 10 + 1E-4
@test -10 - 1E-4 <= newtonGS7[2][1].lo <= -10 + 1E-4
@test -0.13793 - 1E-4 <= newtonGS7[2][1].hi <= -0.13793 + 1E-4
@test newtonGS7[3] == false
@test newtonGS7[4] == false
@test newtonGS7[5] == true
@test 0.285713 - 1E-4 <= newtonGS8[1][1].lo <= 0.285713 + 1E-4
@test 10 - 1E-4 <= newtonGS8[1][1].hi <= 10 + 1E-4
@test -10 - 1E-4 <= newtonGS8[2][1].lo <= -10 + 1E-4
@test -0.13793 - 1E-4 <= newtonGS8[2][1].hi <= -0.13793 + 1E-4
@test newtonGS8[3] == false
@test newtonGS8[4] == false
@test newtonGS8[5] == true

#=
Eflag = false
Iflag = false
eDflag = false
newtonGS9 = PId_NewtonGS(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag,eDflag)
Eflag = false
Iflag = false
krawczykCW9 = PId_KrawczykCW(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag)

@test -1E-4 <= newtonGS9[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= newtonGS9[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= newtonGS9[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= newtonGS9[1][2].hi - 0.208486 <= 1E-4
@test -1E-4 <= krawczykCW9[1][1].lo + 1.04243 <= 1E-4
@test -1E-4 <= krawczykCW9[1][1].hi + 0.492759 <= 1E-4
@test -1E-4 <= krawczykCW9[1][2].lo - 0.0473789 <= 1E-4
@test -1E-4 <= krawczykCW9[1][2].hi - 0.208486 <= 1E-4
=#

end