text:
export TDLearner

using LinearAlgebra: dot
using Distributions: pdf

Base.@kwdef struct TDLearner{A} <: Any
    approximator::A
    γ::Float64 = 1.0
    method::Symbol
    n::Int = 0
end

(L::TDLearner)(env::AbstractEnv) = L.approximator(state(env))
(L::TDLearner)(s) = L.approximator(s)
(L::TDLearner)(s, a) = L.approximator(s, a)

## update policies

function RLBase.update!(
    p::QBasedPolicy{<:TDLearner},
    t::Any,
    e::AbstractEnv,
    s::AbstractStage,
)
    if p.learner.method === :ExpectedSARSA && s === PRE_ACT_STAGE
        # A special case
        update!(p.learner, (t, pdf(prob(p, e))), e, s)
    else
        update!(p.learner, t, e, s)
    end
end

function RLBase.update!(L::TDLearner, t::Any, ::AbstractEnv, s::PreActStage)
    _update!(L, L.approximator, Val(L.method), t, s)
end

function RLBase.update!(L::TDLearner, t::Any, ::AbstractEnv, s::PostEpisodeStage)
    _update!(L, L.approximator, Val(L.method), t, s)
end

# for ExpectedSARSA
function RLBase.update!(
    L::TDLearner,
    t::Tuple,
    ::AbstractEnv,
    s::Union{PreActStage,PostEpisodeStage},
)
    _update!(L, L.approximator, Val(L.method), t, s)
end

## update trajectories

function RLBase.update!(
    t::Any,
    ::Union{
        QBasedPolicy{<:TDLearner},
        NamedPolicy{<:QBasedPolicy{<:TDLearner}},
        VBasedPolicy{<:TDLearner},
    },
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    empty!(t)
end

## implementations

function _update!(
    L::TDLearner,
    ::Union{TabularQApproximator,LinearQApproximator},
    ::Union{Val{:SARSA},Val{:ExpectedSARSA},Val{:SARS}},
    t::Trajectory,
    ::PostEpisodeStage,
)
    S, A, R, T = [t[x] for x in SART]
    n, γ, Q = L.n, L.γ, L.approximator
    G = 0.0
    for i in 1:min(n + 1, length(R))
        G = R[end-i+1] + γ * G
        s, a = S[end-i], A[end-i]
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function _update!(
    L::TDLearner,
    ::Union{TabularQApproximator,LinearQApproximator},
    ::Val{:SARSA},
    t::Trajectory,
    ::PreActStage,
)
    S, A, R, T = [t[x] for x in SART]
    n, γ, Q = L.n, L.γ, L.approximator

    if length(R) >= n + 1
        s, a, s′, a′ = S[end-n-1], A[end-n-1], S[end], A[end]
        G = discount_rewards_reduced(@view(R[end-n:end]), γ) + γ^(n + 1) * Q(s′, a′)
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function _update!(
    L::TDLearner,
    ::TabularQApproximator,
    ::Val{:ExpectedSARSA},
    experience,
    ::PreActStage,
)
    t, p = experience

    S = t[:state]
    A = t[:action]
    R = t[:reward]

    n, γ, Q = L.n, L.γ, L.approximator

    if length(R) >= n + 1
        s, a, s′ = S[end-n-1], A[end-n-1], S[end]
        G = discount_rewards_reduced(@view(R[end-n:end]), γ) + γ^(n + 1) * dot(Q(s′), p)
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function _update!(
    L::TDLearner,
    ::Union{TabularQApproximator,LinearQApproximator},
    ::Val{:SARS},
    t::Any,
    ::PreActStage,
)
    S = t[:state]
    A = t[:action]
    R = t[:reward]

    n, γ, Q = L.n, L.γ, L.approximator

    if length(R) >= n + 1
        s, a, s′ = S[end-n-1], A[end-n-1], S[end]
        G = discount_rewards_reduced(@view(R[end-n:end]), γ) + γ^(n + 1) * maximum(Q(s′))
        update!(Q, (s, a) => Q(s, a) - G)
    end
end

function _update!(
    L::TDLearner,
    ::Union{TabularVApproximator,LinearVApproximator},
    ::Val{:SRS},
    t::Trajectory,
    ::PostEpisodeStage,
)
    S, R = t[:state], t[:reward]
    n, γ, V = L.n, L.γ, L.approximator
    G = 0.0
    w = 1.0
    for i in 1:min(n + 1, length(R))
        G = R[end-i+1] + γ * G
        s = S[end-i]
        if haskey(t, :weight)
            w *= t[:weight][end-i]
        end
        update!(V, s => w * (V(s) - G))
    end
end

function _update!(
    L::TDLearner,
    ::Union{TabularVApproximator,LinearVApproximator},
    ::Val{:SRS},
    t::Any,
    ::PreActStage,
)
    S = t[:state]
    R = t[:reward]

    n, γ, V = L.n, L.γ, L.approximator
    if length(R) >= n + 1
        s, s′ = S[end-n-1], S[end]
        G = discount_rewards_reduced(@view(R[end-n:end]), γ) + γ^(n + 1) * V(s′)
        if haskey(t, :weight)
            W = t[:weight]
            @views w = reduce(*, W[end-n-1:end-1])
        else
            w = 1.0
        end
        update!(V, s => w * (V(s) - G))
    end
end

#####
# DynaAgent
#####

function RLBase.update!(
    p::QBasedPolicy{<:TDLearner},
    m::Union{ExperienceBasedSamplingModel,TimeBasedSamplingModel},
    ::Any,
    env::AbstractEnv,
    ::Union{PreActStage,PostEpisodeStage},
)
    if p.learner.method == :SARS
        transition = sample(m)
        if !isnothing(transition)
            s, a, r, t, s′ = transition
            traj = VectorSARTTrajectory()
            push!(traj; state = s, action = a, reward = r, terminal = t)
            push!(traj; state = s′, action = a)  # here a is a dummy one
            update!(p.learner, traj, env, t ? POST_EPISODE_STAGE : PRE_ACT_STAGE)
        end
    else
        @error "unsupported method $(p.learner.method)"
    end
end

function RLBase.update!(
    p::QBasedPolicy{<:TDLearner},
    m::PrioritizedSweepingSamplingModel,
    ::Any,
    env::AbstractEnv,
    ::Union{PreActStage,PostEpisodeStage},
)
    if p.learner.method == :SARS
        transition = sample(m)
        if !isnothing(transition)
            s, a, r, t, s′ = transition
            traj = VectorSARTTrajectory()
            push!(traj; state = s, action = a, reward = r, terminal = t)
            push!(traj; state = s′, action = a)  # here a is a dummy one
            update!(p.learner, traj, env, t ? POST_EPISODE_STAGE : PRE_ACT_STAGE)

            # update priority
            for (s̄, ā, r̄, d̄) in m.predecessors[s]
                P = RLBase.priority(p.learner, (s̄, ā, r̄, d̄, s))
                if P ≥ m.θ
                    m.PQueue[(s̄, ā)] = P
                end
            end
        end
    else
        @error "unsupported method $(p.learner.method)"
    end
end

function RLBase.priority(L::TDLearner, transition::Tuple)
    if L.method == :SARS
        s, a, r, d, s′ = transition
        γ, Q = L.γ, L.approximator
        Δ = d ? (r - Q(s, a)) : (r + γ^(L.n + 1) * maximum(Q(s′)) - Q(s, a))
        Δ = [Δ]  # must be broadcastable in Flux.Optimise
        Flux.Optimise.apply!(Q.optimizer, (s, a), Δ)
        abs(Δ[])
    else
        @error "unsupported method"
    end
end

#####
# TDλReturnLearner
#####

export TDλReturnLearner

Base.@kwdef struct TDλReturnLearner{Tapp} <: Any
    approximator::Tapp
    γ::Float64 = 1.0
    λ::Float64
end

(L::TDλReturnLearner)(env::AbstractEnv) = L(state(env))
(L::TDλReturnLearner)(s) = L.approximator(s)
(L::TDλReturnLearner)(s, a) = L.approximator(s, a)

function RLBase.update!(L::TDλReturnLearner, t::Any, ::AbstractEnv, ::PreActStage) end

function RLBase.update!(L::TDλReturnLearner, t::Any, ::AbstractEnv, ::PostEpisodeStage)
    λ, γ, V = L.λ, L.γ, L.approximator
    R = t[:reward]
    S = @view t[:state][1:end-1]
    S′ = @view t[:state][2:end]
    T = length(R)
    for t in 1:T
        G = 0.0
        for n in 1:(T-t)
            G +=
                λ^(n - 1) *
                (discount_rewards_reduced(@view(R[t:t+n-1]), γ) + γ^n * V(S′[t+n-1]))
        end
        G *= 1 - λ
        G +=
            λ^(T - t) *
            (discount_rewards_reduced(@view(R[t:T]), γ) + γ^(T - t + 1) * V(S′[T]))
        sₜ = S[t]
        update!(V, sₜ => V(sₜ) - G)
    end
end

function RLBase.update!(
    t::Any,
    ::VBasedPolicy{<:TDλReturnLearner},
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    empty!(t)
end


text:
## Tests for compatiability with the old form
# test for not loaded
@test_throws UndefVarError Poly([1,2,3])
@test_throws UndefVarError poly([1,2,3])
@test_throws UndefVarError polyval(Polynomial([1,2,3]), 1)
@test_throws UndefVarError polyder(Polynomial([1,2,3]))
@test_throws UndefVarError polyint(Polynomial([1,2,3]))
a = Polynomial(1 .//convert(Vector{BigInt},map(gamma,BigFloat(1):BigFloat(17))),"x")
@test_throws UndefVarError Pade(a,8,8)

# must load the compatability module to have access
using Polynomials.PolyCompat

## the methods  `polyval`, `polyint` and `polyder` are only for the
## legacy `Poly`  type.
# check method are defined
p = Poly([1,2,3])
q = poly([1,2,3])
polyval(p, 1)
polyint(p)
polyder(p)
a = Poly(1 .//convert(Vector{BigInt},map(gamma,BigFloat(1):BigFloat(17))),"x")
pp = Pade(a,8,8)
padeval(pp, 1)

p = Polynomial([1,2,3])
@test_throws ErrorException polyval(p,1)
@test_throws ErrorException polyint(p)
@test_throws ErrorException polyder(p)

## ---------
## Tests for  Poly, and Pade types once PolyCompat is loaded

pNULL = Poly(Float32[])
p0 = Poly([0])
p1 = Poly([1,0,0,0,0,0,0,0,0,0,0,0,0,0])
p2 = Poly([1,1,0,0])
p3 = Poly([1,2,1,0,0,0,0])
p4 = Poly([1,3,3,1,0,0])
p5 = Poly([1,4,6,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0])
pN = Poly([276,3,87,15,24,0])
pR = Poly([3//4, -2//1, 1//1])
X = Poly([0.0, 1.0])
T = Int64
# Poly{T}([zero(T), one(T)])
# Poly{T}([zero(T), one(T)], :y)

p1000 = Poly(randn(1000))

@test length(pNULL) == 1
@test length(p1000-p1000) == 1
@test length(p1000^0) == 1
@test length(0*p1000) == 1
@test length(p1000) == 1000
sprint(show, p1000)
sprint(show, pNULL)

@test Poly([1, -2, 1]) == poly(diagm(0=>[1., 1.]))

@test p3 == Poly([1,2,1])
@test pN*10 == Poly([2760, 30, 870, 150, 240])
@test pN/10 == Poly([27.6, 0.3, 8.7, 1.5, 2.4])
@test 10*pNULL + pN == pN
@test 10*p0 + pN == pN
@test p5 + 2*p1 == Poly([3,4,6,4,1])
@test 10*pNULL - pN == -pN
@test p0 - pN == -pN
@test p5 - 2*p1 == Poly([-1,4,6,4,1])
@test p2*p2*p2 == p4
@test p2^4 == p5
@test pNULL^3 == pNULL
@test pNULL*pNULL == pNULL

@test map(degree, [pNULL,p0,p1,p2,p3,p4,p5,pN,pR,p1000]) == [-1,-1,0,1,2,3,4,4,2,999]

@test polyval(poly(Int[]), 2.) == 1.
@test polyval(pN, -.125) == 276.9609375
@test polyval(pNULL, 10) == 0
@test polyval(p0, -10) == 0
@test isa(polyval(p0, X), Poly)
@test isa(polyval(p1, X), Poly)
@test polyval(poly([1//2, 3//2]), 1//2) == 0//1
@test polyder(polyint(pN)) == pN
@test polyder(pR) == Poly([-2//1,2//1])
@test polyder(p3) == Poly([2,2])
@test polyder(p1) == polyder(p0) == polyder(pNULL) == pNULL
@test_throws ErrorException polyder(pR, -1)
@test polyint(pNULL,1) == p1
@test polyint(Poly(Rational[1,2,3])) == Poly(Rational[0, 1, 1, 1])
@test polyint(p2, 0, 2) == 4.0

@test pN(-.125) == 276.9609375
@test pN.([0.1, 0.2, 0.3]) == polyval.(pN, [0.1, 0.2, 0.3])

@test poly([-1,-1]) == p3
@test roots(p0)==roots(p1)==roots(pNULL)==[]
@test roots(Poly([0,1,0])) == [0.0]
r = roots(Poly([0,1,0]))
@test roots(p2) == [-1]
a_roots = copy(coeffs(pN))
@test all(map(abs,sort(roots(poly(a_roots))) - sort(a_roots)) .< 1e6)
@test length(roots(p5)) == 4
@test roots(pNULL) == []
@test sort(roots(pR)) == [1//2, 3//2]
# default type in variable() is Float64
x = variable()
plarge = 8.362779449448982e41 - 2.510840694154672e57x + 4.2817430781178795e44x^2 - 1.6225927682921337e31x^3 + 1.0x^4  # #120
@test length(roots(plarge)) == 4

@test pNULL + 2 == p0 + 2 == 2 + p0 == Poly([2])
@test p2 - 2 == -2 + p2 == Poly([-1,1])
@test 2 - p2 == Poly([1,-1])

p0 = Poly([0])
p1 = Poly([1])
p2 = Poly([5, 6, -3, 2 ,4])
p3 = Poly([7, -3, 2, 6])
p4 = p2 * p3
@test divrem(p4, p2) == (p3, zero(p3))
@test p3%p2 == p3
@test all((map(abs,(p2 ÷ p3 - Poly([1/9,2/3]))|> coeffs)) .< eps())
@test divrem(p0,p1) == (p0,p0)
@test divrem(p1,p1) == (p1,p0)
@test divrem(p2,p2) == (p1,p0)
@test divrem(pR, pR) == (one(pR), zero(pR))
@test_throws DivideError p1 ÷ p0
@test_throws DivideError divrem(p0,p0)

# test chop
pchop = Poly([1, 2, 3, 0, 0, 0])
pchopped = chop(pchop)
@test roots(pchop) == roots(pchopped)

#Tests for multivariable support
pX = Poly([1, 2, 3, 4, 5])
pS1 = Poly([1, 2, 3, 4, 5], "s")
pS2 = Poly([1, 2, 3, 4, 5], 's')
pS3 = Poly([1, 2, 3, 4, 5], :s)
@test pX != pS1
@test pS1 == pS2
@test pS1 == pS3
@test_throws ErrorException pS1 + pX
@test_throws ErrorException pS1 - pX
@test_throws ErrorException pS1 * pX
@test_throws ErrorException pS1 ÷ pX
@test_throws ErrorException pS1 % pX

#Testing copying.
pcpy1 = Poly([1,2,3,4,5], :y)
pcpy2 = copy(pcpy1)
@test pcpy1 == pcpy2

#Tests for Pade approximants

a = Poly(1 .//convert(Vector{BigInt},map(gamma,BigFloat(1):BigFloat(17))),"x")
PQexp = Pade(a,8,8)
@test isapprox(convert(Float64, padeval(PQexp,1.0)), exp(1.0))
@test isapprox(convert(Float64, padeval(PQexp,-1.0)), exp(-1.0))

b = Poly(convert(Vector{BigInt},map(sinpi,(0:16)/2)).//convert(Vector{BigInt},map(gamma,BigFloat(1):BigFloat(17))),"x")
PQsin = Pade(b,8,7)
@test isapprox(convert(Float64, padeval(PQsin,1.0)), sin(1.0))
@test isapprox(convert(Float64, padeval(PQsin,-1.0)),sin(-1.0))

c = Poly(convert(Vector{BigInt},map(sinpi,(1:17)/2)).//convert(Vector{BigInt},map(gamma,BigFloat(1):BigFloat(17))),"x")
PQcos = Pade(c,8,8)
@test isapprox(convert(Float64, padeval(PQcos,1.0)), cos(1.0))
@test isapprox(convert(Float64, padeval(PQcos,-1.0)), cos(-1.0))

const _γ = 0.5772156649015

d = Poly(convert(Vector{BigInt},(-1).^(0:60).* map(gamma,BigFloat(1):BigFloat(61.0))).//1,"x")
PQexpint = Pade(d,30,30)

@test isapprox(convert(Float64, padeval(PQexpint,1.0)),
               exp(1)*(-_γ-sum([(-1).^k/k./gamma(k+1) for k=1:20])))

## polyfit(xs,ys) -> fit(Poly, xs, ys)
xs = range(0, stop=pi, length=10)

ys = map(sin,xs)
p = fit(Poly, xs, ys)
p = fit(Poly, xs, ys, var=:t)
p = fit(Poly, xs, ys, 2)
@test maximum(map(abs,map(x->polyval(p, x), xs) - ys)) <= 0.03
#https://stackoverflow.com/questions/50832823/error-with-polyfit-function-julia
# relax type assumptions on x, y
fit(Poly, Number[1,2,3], Number[2,3,1])

## truncation
p1 = Poly([1,1]/10)
p2 = Poly([1,2]/10)
p3 = Poly([1,3]/10)
psum = p1 + p2 - p3
@test degree(psum) == 1         # will have wrong degree
@test degree(truncate(psum)) == 0 # the degree should be correct after truncation

@test truncate(Poly([2,1]),rtol=1/2,atol=0) == Poly([2])
@test truncate(Poly([2,1]),rtol=1,atol=0)   == Poly([0])
@test truncate(Poly([2,1]),rtol=0,atol=1)   == Poly([2])

@test norm(Poly([1., 2.])) == norm([1., 2.])
@test norm(Poly([1., 2.]), 1) == norm([1., 2.], 1)

## setindex!
p1    = Poly([1,2,1])
p1[5] = 1
@test p1[5] == 1
@test p1 == Poly([1,2,1,0,0,1])

p1[:] = 0
@test p1 ≈ zero(p1)
@test zero(Poly{Int}) == Poly(Int[])
@test one(Poly{Int}) == Poly([1])
p1[:] = [1,2,1,0,0,1]
@test p1 == Poly([1,2,1,0,0,1])

## elementwise operations #52
p1  = Poly([1, 2])
p2  = Poly([3, 1.])
p   = [p1, p2]
q   = [3, p1]
@test isa(q,Vector{Poly{Int64}})
psum  = p .+ 3
pprod = p .* 3
pmin  = p .- 3
@test isa(psum, Vector{Poly{Float64}})
@test isa(pprod,Vector{Poly{Float64}})
@test isa(pmin, Vector{Poly{Float64}})

## getindex with ranges #43
p1 = Poly([4,5,6])
@test all(p1[0:1] .== [4,5])
@test all(p1[0:end] .== [4,5,6])
p1[0:1] = [7,8]
@test all(p1[0:end] .== [7,8,6])

@test p1[:] == coeffs(p1)

## conjugate of poly (issue #59)
as = [im, 1, 2]
bs = [1, 1, 2]
@test conj(Poly(as)) == Poly(conj.(as))
#@test conj(Poly(bs)) == Poly(conj(bs)) # conj gives warning as no defn on T
## and transpose
@test transpose(Poly(as)) == Poly(as)

## unnecessary copy in convert #65
p1 = Poly([1,2])
p2 = convert(Poly{Int64}, p1)
p2[3] = 3
@test p1[3] == 3

## Polynomials with non-Real type
import Base: +, *, -
struct Mod2 <: Number
  v::Bool
end
+(x::Mod2,y::Mod2) = Mod2(xor(x.v, y.v))
*(x::Mod2,y::Mod2) = Mod2(x.v&y.v)
-(x::Mod2,y::Mod2) = x+y
-(x::Mod2) = x
Base.one(::Type{Mod2}) = Mod2(true)
Base.zero(::Type{Mod2}) = Mod2(false)
Base.convert(::Type{Mod2},x::Integer) = Mod2(convert(Bool,x))
Base.convert(::Type{Bool},x::Mod2) = x.v

# Test that none of this throws
p = Poly([Mod2(true),Mod2(false), Mod2(true)])
repr(p)
@test p(Mod2(false)) == p[0]

## changes to show
p = Poly([1,2,3,1])  # leading coefficient of 1
@test repr(p) == "Poly(1 + 2*x + 3*x^2 + x^3)"
p = Poly([1.0, 2.0, 3.0, 1.0])
@test repr(p) == "Poly(1.0 + 2.0*x + 3.0*x^2 + 1.0*x^3)"
p = Poly([1, im])
@test repr(p) == "Poly(1 + im*x)"
p = Poly([1+im, 1-im, -1+im, -1 - im])# minus signs
@test repr(p) == "Poly((1 + 1im) + (1 - 1im)*x - (1 - 1im)*x^2 - (1 + 1im)*x^3)"
p = Poly([1.0, 0 + NaN*im, NaN, Inf, 0 - Inf*im]) # handle NaN or Inf appropriately
@test repr(p) == "Poly(1.0 + NaN*im*x + NaN*x^2 + Inf*x^3 - Inf*im*x^4)"

p = Poly([1,2,3])

@test repr("text/latex", p) == "\$1 + 2\\cdot x + 3\\cdot x^{2}\$"
p = Poly([1//2, 2//3, 1])
@test repr("text/latex", p) == "\$\\frac{1}{2} + \\frac{2}{3}\\cdot x + x^{2}\$"

# customized printing with printpoly
function printpoly_to_string(args...; kwargs...)
    buf = IOBuffer()
    printpoly(buf, args...; kwargs...)
    return String(take!(buf))
end
@test printpoly_to_string(Poly([1,2,3], "y")) == "1 + 2*y + 3*y^2"
@test printpoly_to_string(Poly([1,2,3], "y"), descending_powers=true) == "3*y^2 + 2*y + 1"
@test printpoly_to_string(Poly([2, 3, 1], :z), descending_powers=true, offset=-2) == "1 + 3*z^-1 + 2*z^-2"
@test printpoly_to_string(Poly([-1, 0, 1], :z), offset=-1, descending_powers=true) == "z - z^-1"

## want to be able to copy and paste
## check hashing
p = poly([1,2,3])
q = poly([1,2,3])
@test hash(p) == hash(q)

## Check for Inf/NaN operations
p1 = Poly([Inf, Inf])
p2 = Poly([0, Inf])
@test p1(Inf) == Inf
@test isnan(p1(-Inf))
@test isnan(p1(0))
@test p2(-Inf) == -Inf

## Check for isequal
p1 = Poly([-0., 5., Inf])
p2 = Poly([0., 5., Inf])
p3 = Poly([0, NaN])

@test p1 == p2 && !isequal(p1, p2)
@test p3 === p3 && p3 ≠ p3 && isequal(p3, p3)

## Handling of `NaN`s
p     = Poly([NaN, 1, 5])
pder  = polyder(p)
pint  = polyint(p)

@test isnan(p(1))                 # p(1) evaluates to NaN
@test isequal(pder, Poly([NaN]))  # derivative will give Poly([NaN])
@test isequal(pint, Poly([NaN]))  # integral will give Poly([NaN])

pint  = polyint(p, Complex(0.))
@test isequal(pint, Poly([NaN]))  # integral will give Poly([NaN])

## proper conversions in arithmetic with different element-types #94
p = Poly([0,one(Float64)])
@test Poly{Complex{Float64}} == typeof(p+1im)
@test Poly{Complex{Float64}} == typeof(1im-p)
@test Poly{Complex{Float64}} == typeof(p*1im)

## comparison between `Number`s and `Poly`s
p1s = Poly([1,2], :s)
p1x = Poly([1,2], :x)
p2s = Poly([1], :s)

@test p1s == p1s
@test p1s ≠ p1x
@test p1s ≠ p2s

@test_throws ErrorException p1s ≈ p1x
@test p1s ≉ p2s
@test p1s ≈ Poly([1,2.], :s)

@test p2s ≈ 1. ≈ p2s
@test p2s == 1. == p2s
@test p2s ≠ 2. ≠ p2s
@test p1s ≠ 2. ≠ p1s

@test nnz(map(Poly, sparse(1.0I, 5,5))) == 5

@test Poly([0.5]) + 2 == Poly([2.5])
@test 2 - Poly([0.5]) == Poly([1.5])

# test size
@test size(Poly([0.5, 0.2])) == (2,)
@test size(Poly([0.5, 0.2]), 1) == 2
# @test size(Poly([0.5, 0.2]), 1, 2) == (2,1) # deprecated in v0.7

# test iteration
p1 = Poly([1,2,0,3])
for term in p1
  @test isa(term, Poly)
end

@test eltype(typeof(p1)) == typeof(p1)

# after implementing the iteration interface, i.e., `start` and its friends,
# previously throwing `collect` function (and maybe, other similar functions)
# started giving different errors due to `eltype{T}(::Poly{T}) = T`. Changing
# this definition has resulted in JuliaLang/METADATA.jl#8528. One small fix
# was to direct `collect{T}(p::Poly{T})` to `collect(Poly{T}, p)`.

@test eltype(p1) == Int
@test eltype(collect(p1)) == Poly{Int}
@test eltype(collect(Poly{Float64}, p1)) == Poly{Float64}
@test_throws InexactError collect(Poly{Int}, Poly([1.2]))

@test length(collect(p1)) == degree(p1)+1

@test [p1[idx] for idx in eachindex(p1)] == [1,2,0,3]

p1 = Poly([2.,5.,1.])
p2 = Poly([1.,2.,3.])

@test degree(gcd(p1, p2))                       == 0          # no common roots
@test degree(gcd(p1, Poly(5)))                  == 0          # ditto
@test degree(gcd(p1, Poly(eps(0.))))            == 0          # ditto
@test degree(gcd(p1, Poly(0)))                  == degree(p1) # Poly(0) has the roots of p1
@test degree(gcd(p1+p2*170.10734737144486, p2)) == 0          # see, c.f., #122

p1 = poly([1.,2.,3.])
p2 = poly([1.,2.,6.])

@test (res = roots(gcd(p1, p2)); 1. ∈ res && 2. ∈ res)

## Getting error on passing Real arrays to polyfit #146
xx = Real[20.0, 30.0, 40.0]
yy = Real[15.7696, 21.4851, 28.2463]
fit(Poly, xx,yy,2)

## Issue with overflow and polyder Issue #159
@test !iszero(polyder(Poly(BigInt[0, 1])^100, 100))

@testset "Pade" begin
    # exponential
    coeffs = 1 .// BigInt.(gamma.(1:17))
    a = Polynomial(coeffs)
    PQexp = Pade(a, 8, 8)
    @test PQexp(1.0) ≈ exp(1.0)
    @test PQexp(-1.0) ≈ exp(-1.0)

    # sine
    coeffs = BigInt.(sinpi.((0:16) ./ 2)) .// BigInt.(gamma.(1:17))
    p = Polynomial(coeffs)
    PQsin = Pade(p, 8, 7)
    @test PQsin(1.0) ≈ sin(1.0)
    @test PQsin(-1.0) ≈ sin(-1.0)

    # cosine
    coeffs = BigInt.(sinpi.((1:17) ./ 2)) .// BigInt.(gamma.(1:17))
    p = Polynomial(coeffs)
    PQcos = Pade(p, 8, 8)
    @test PQcos(1.0) ≈ cos(1.0)
    @test PQcos(-1.0) ≈ cos(-1.0)

    # summation of a factorially divergent series
    γ = 0.5772156649015
    s = BigInt.(gamma.(BigInt(1):BigInt(61)))
    coeffs = (BigInt(-1)).^(0:60) .* s .// 1
    d = Polynomial(coeffs)
    PQexpint = Pade(d, 30, 30)
    @test Float64(PQexpint(1.0)) ≈ exp(1) * (-γ - sum([(-1)^k / k / gamma(k + 1) for k = 1:20]))
end

ext: jl


text:
module Pitchjx

export pitchjx

using DataFrames, Dates

include("extractor.jl")

"""
Scrape pitching data.

# How to use

`pitchjx("2018-10-20")`
"""
function pitchjx(start, fin=start)
    @info "Initialize: Start"
    params = Dict(
        "startDate" => start,
        "endDate" => fin,
        "sportId" => 1
    )
    @info "Initialize: Finish!"
    @info "Extract dataset: Start"
    result = extract(params)
    @info "Extract dataset: Finish!"
    return result
end

end


text:
using SolidStateDetectors
using IntervalSets
using SpecialFunctions
T=Float32

struct CustomChargeDensity{T, N} <: SolidStateDetectors.AbstractChargeDensity{T}
	G1::SolidStateDetectors.AbstractGeometry{T,N}
	G2::SolidStateDetectors.AbstractGeometry{T,N}
	L1::SolidStateDetectors.AbstractChargeDensity{T}
	L2::SolidStateDetectors.AbstractChargeDensity{T}
	G3::SolidStateDetectors.AbstractGeometry{T,N}
	L3::SolidStateDetectors.AbstractChargeDensity{T}
end

function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::SolidStateDetectors.AbstractCoordinatePoint{T})::T where {T}
	(pt in cdm.G1) * SolidStateDetectors.get_charge_density(cdm.L1, pt) +
	(pt in cdm.G2) * SolidStateDetectors.get_charge_density(cdm.L2, pt) +
	(pt in cdm.G3) * SolidStateDetectors.get_charge_density(cdm.L3, pt)
end

G1 = SolidStateDetectors.Tube{T}(0..0.00725,0..2π,0..0.002,CartesianPoint{T}(0.0, 0.0, 0))
G2 = SolidStateDetectors.Tube{T}(0..0.00725,0..2π,0..1e-05,CartesianPoint{T}(0.0, 0.0, 0))
CD1 = SolidStateDetectors.CylindricalChargeDensity{T}((0.0f0,0.0f0,50000000000.0),(3e+24,0.0f0,0.0f0))
CD2 = SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,-1.9e+16),(0.0f0,0.0f0,0.0f0))
G3 = SolidStateDetectors.HexagonalPrism{T}(0, 0.00515, 2e-05,CartesianPoint{T}(0, 0, 0.00199),0)
CD3= SolidStateDetectors.LinearChargeDensity{T}((0.0f0,0.0f0,1.9e+16),(0.0f0,0.0f0,0.0f0))

ccdm = CustomChargeDensity(G1,G2,CD1,CD2,G3,CD3)

max_stars_repo_name: leenderthayen/NDP

max_issues_repo_name: leenderthayen/NDP

max_forks_repo_name: leenderthayen/NDP


text:
# This is part of Calendars.jl. See the copyright note there.
# ===================== European calendar ===================

# The Gregorian calendar is used for dates on and after
# CE 1582-10-15, otherwise the Julian calendar is used.

# Is the date a valid European date?
function isValidDateEuropean(cd::CDate, warn=true)
    cal, year, month, day = cd
    if cal != EC
        warn && @warn(Warning(cd))
        return false
    end

    dn = DayNumberGregorian(year, month, day)
    if dn < GregorysBreak
        return isValidDateJulian((JD, year, month, day), warn)
    end
    return isValidDateGregorian((CE, year, month, day), warn)
end

# Computes the day number from a date.
# All dates prior to CE 1582-10-15 are interpreted as Julian dates.
function DayNumberEuropean(year::DPart, month::DPart, day::DPart)
    dn = DayNumberGregorian(year, month, day)
    if dn < GregorysBreak
        return DayNumberJulian(year, month, day)
    end
    return dn
end

# Computes the day number from a valid date.
function DayNumberEuropean(cd::CDate)
    cal, year, month, day = cd
    return DayNumberEuropean(year, month, day)
end

# Computes the historic date from a day number.
# All dates prior to CE 1582-10-15 are returned as Julian dates.
function DateEuropean(dn::DPart)
    d = dn < GregorysBreak ? DateJulian(dn) : DateGregorian(dn)
    cal, year, month, day = d
    return (EC, year, month, day)
end

# Return the days this year so far.
function DayOfYearEuropean(year::DPart, month::DPart, day::DPart)
    dn = DayNumberGregorian(year, month, day)
    if dn < GregorysBreak
        return DayOfYearJulian(year, month, day)
    end
    return DayOfYearGregorian(year, month, day)
end

# Return the day number of the first day in the given European year.
function YearStartEuropean(year::DPart)
    if year <= 1582
        return YearStartJulian(year)
    end
    return YearStartGregorian(year)
end

# Return the day number of the last day in the given European year.
function YearEndEuropean(year::DPart)
    if year < 1582
        return YearEndJulian(year)
    end
    return YearEndGregorian(year)
end

function MonthsInYearEuropean(year::DPart)
    return 12
end

function DaysInYearEuropean(year::DPart)
    if year < 1582
        return DaysInYearJulian(year)
    end
    return DaysInYearGregorian(year)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using ValueShapes, Distributions

@testset "variate_transform" begin
    @testset "identity_vt" begin
        @test @inferred(BAT.IdentityVT(BAT.InfiniteSpace(), ScalarShape{Real}())) isa BAT.VariateTransform{Univariate,BAT.InfiniteSpace,BAT.InfiniteSpace}
        @test @inferred(BAT.IdentityVT(BAT.UnitSpace(), ArrayShape{Int}(2))) isa BAT.VariateTransform{Multivariate,BAT.UnitSpace,BAT.UnitSpace}
        @test @inferred(BAT.IdentityVT(BAT.MixedSpace(), ArrayShape{Float32}(2, 3))) isa BAT.VariateTransform{Matrixvariate,BAT.MixedSpace,BAT.MixedSpace}

        ntshape = NamedTupleShape(a = ScalarShape{Real}(), b = ArrayShape{Int}(2))
        @test @inferred(BAT.IdentityVT(BAT.MixedSpace(), ntshape)) isa BAT.VariateTransform{<:ValueShapes.NamedTupleVariate{(:a, :b)},BAT.MixedSpace,BAT.MixedSpace}

        nttrafo = BAT.IdentityVT(BAT.MixedSpace(), ntshape)
        ntvalue = (a = 4.2, b = [5, 7])

        @test @inferred(varshape(nttrafo)) == ntshape
        @test @inferred(inv(nttrafo)) === nttrafo

        @test @inferred(BAT.apply_vartrafo(nttrafo, ntvalue, 0)).v === ntvalue
        @test BAT.apply_vartrafo(nttrafo, ntvalue, 0.79).ladj == 0.79
        @test isnan(BAT.apply_vartrafo(nttrafo, ntvalue, NaN).ladj)
        @test @inferred(BAT.apply_vartrafo(inv(nttrafo), ntvalue, 0.79)).v === ntvalue
        @test BAT.apply_vartrafo(inv(nttrafo), ntvalue, 0.79).ladj == 0.79
        @test isnan(BAT.apply_vartrafo(inv(nttrafo), ntvalue, NaN).ladj)

        @test @inferred(nttrafo ∘ nttrafo) === nttrafo
    end
end


text:
using SimplePolynomials, BigCombinatorics

function sine_coeffs(n::Int)
    T = Rational{BigInt}
    result = zeros(T,2n)
    for j=1:n
        result[2j] = (-1)^(j-1)//Factorial(2j-1)
    end

    return result
end

"""
`sine_poly(n)` create exact sine power series of degree `2n`
(actually `2n-1` because `sin` is odd).
"""
function sine_poly(n::Int)
    return SimplePolynomial(sine_coeffs(n))
end


text:
module Searches

using App, SearchLight

export Search

mutable struct Search <: AbstractModel
  ### internals
  _table_name::String
  _id::String

  ### fields
  id::Nullable{SearchLight.DbId}
  search::String

  ### constructor
  Search(;
    id = Nullable{SearchLight.DbId}(),
    search = ""
  ) = new("searches", "id", id, search)
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
export dTtransposedMultiplication, QtransposedMultiplication, computeParametricHessian, hessianFunction

function dTtransposedMultiplication(N::Int,dX::Array{Float64,1},dY::Array{Float64,1},r::Array{Float64,1})
    result = zeros(2*N)
    @simd for i=1:N
        @inbounds result[i] = dX[i] * r[i];
        @inbounds result[N+i] = dY[i] * r[i];
    end
    return result
end

# Q' * hdTr with hdTr = h * dT' * residual
function QtransposedMultiplication(N::Int,hdTr::Array{Float64,1},centeredGrid::Array{Float64,1})
    result = zeros(6)
    @simd for i=1:N
        @inbounds result[1] += centeredGrid[i]*hdTr[i]
        @inbounds result[2] += centeredGrid[N+i]*hdTr[i]
        @inbounds result[3] += hdTr[i]
        @inbounds result[4] += centeredGrid[i]*hdTr[N+i]
        @inbounds result[5] += centeredGrid[N+i]*hdTr[N+i]
        @inbounds result[6] += hdTr[N+i]
    end
    return result
end

# compute parametric hessian
# result = h*(dT*Q)'*(dT*Q) = h* Q' * dT' dT * Q
function computeParametricHessian(N::Int,h::Float64,dTX::Array{Float64,1},dTY::Array{Float64,1},centeredGrid::Array{Float64,1})
    result = zeros(6,6)
    help = zeros(2*N); help[1:N] = centeredGrid[1:N]
    result[:,1] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    help[1:N] = centeredGrid[N+1:end]
    result[:,2] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    help[1:N] = ones(N)
    result[:,3] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    help = zeros(2*N); help[N+1:end] = centeredGrid[1:N]
    result[:,4] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    help[N+1:end] = centeredGrid[N+1:end]
    result[:,5] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    help[N+1:end] = ones(N)
    result[:,6] = QtransposedMultiplication(N,hessianFunction(N,h,dTX,dTY,help),centeredGrid)
    return result
end

# result = h*H*x, H = h*dT'*dT
function hessianFunction(N::Int,h::Float64,dTX::Array{Float64,1},dTY::Array{Float64,1},x::Array{Float64,1})
    result = zeros(2*N)
    @simd for i=1:N
        @inbounds result[i] = h*dTX[i]*dTX[i]*x[i] + h*dTX[i]*dTY[i]*x[N+i]
        @inbounds result[N+i] = h*dTX[i]*dTY[i]*x[i] + h*dTY[i]*dTY[i]*x[N+i]
    end
    return result
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module PeaksAndValleys

using Compat

export PEAK, VALLEY
export identify_initial_pivot, peak_valley_pivots, max_drawdown
export compute_segment_returns, pivots_to_modes

@compat const PEAK, VALLEY = Int8(1), Int8(-1)

"""Quickly identify the xs[1] as a peak or valley."""
function identify_initial_pivot(xs, up_thresh, down_thresh)
    x_1 = xs[1]
    max_x = x_1
    max_t = 1
    min_x = x_1
    min_t = 1
    up_thresh += 1
    down_thresh += 1

    for (t, x_t) in enumerate(xs[2:end])
        t += 1  # I wish Julia had an enumerate(*,from=X)

        if x_t / min_x >= up_thresh
            return min_t == 1 ? VALLEY : PEAK
        end

        if x_t / max_x <= down_thresh
            return max_t == 1 ? PEAK : VALLEY
        end

        if x_t > max_x
            max_x = x_t
            max_t = t
        end

        if x_t < min_x
            min_x = x_t
            min_t = t
        end
    end
    x_1 < xs[length(xs)] ? VALLEY : PEAK
end

"""
Finds the peaks and valleys of a series.

Parameters
---
- **xs** : This is your series.
- **up_thresh** : The minimum relative change necessary to define a peak.
- **down_thesh** : The minimum relative change necessary to define a valley.

Returns
---
an array with 0 indicating no pivot and -1 and 1 indicating valley and peak
respectively.

The First and Last Elements
---
The first and last elements are guaranteed to be annotated as peak or
valley even if the segments formed do not have the necessary relative
changes. This is a tradeoff between technical correctness and the
propensity to make mistakes in data analysis. The possible mistake is
ignoring data outside the fully realized segments, which may bias analysis.
"""
function peak_valley_pivots(xs, up_thresh, down_thresh)
    assert(down_thresh < 0) # "It must be negative.

    initial_pivot = identify_initial_pivot(xs, up_thresh, down_thresh)

    t_n = length(xs)
    pivots = zeros(Int8, t_n)
    pivots[1] = initial_pivot

    # Adding one to the relative change thresholds saves operations. Instead
    # of computing relative change at each point as x_j / x_i - 1, it is
    # computed as x_j / x_1. Then, this value is compared to the threshold + 1.
    # This saves (t_n - 1) subtractions.
    up_thresh += 1
    down_thresh += 1

    trend = -initial_pivot
    last_pivot_t = 1
    last_pivot_x = xs[1]
    for (t, x) in enumerate(xs[2:end])
        t += 1
        r = x / last_pivot_x

        if trend == -1
            if r >= up_thresh
                pivots[last_pivot_t] = trend
                trend = 1
                last_pivot_x = x
                last_pivot_t = t
            elseif x < last_pivot_x
                last_pivot_x = x
                last_pivot_t = t
            end
        else
            if r <= down_thresh
                pivots[last_pivot_t] = trend
                trend = -1
                last_pivot_x = x
                last_pivot_t = t
            elseif x > last_pivot_x
                last_pivot_x = x
                last_pivot_t = t
            end
        end
    end

    pivots[t_n] == 0 && (pivots[t_n] = trend)

    return pivots
end

"""
Return the absolute value of the maximum drawdown of sequence X.

Note
---
If the sequence is increasing, 0 is returned."""
function max_drawdown(xs)
    mdd = 0
    peak = xs[1]
    for x in xs
        if x > peak
            peak = x
        end
        dd = (peak - x) / peak
        if dd > mdd
            mdd = dd
        end
    end
    return mdd
end

"""Return a numpy array of the pivot-to-pivot returns for each segment."""
function compute_segment_returns(xs, pivots)
    pivot_points = xs[pivots .!= 0]
    return pivot_points[2:end] ./ pivot_points[1:end-1] - 1.0
end

"""
Translate pivots into trend modes.

Parameters
---
pivots : the result of calling `peak_valley_pivots`

Returns
---
A vector of trend modes. That is, between (VALLEY, PEAK] it is 1 and
between (PEAK, VALLEY] it is -1.
"""
function pivots_to_modes(pivots)
    modes = zeros(Int8, length(pivots))
    modes[1] = pivots[1]
    mode = -modes[1]
    for (t, x) in enumerate(pivots[2:end])
        if x != 0
            modes[t+1] = mode
            mode = -x
        else
            modes[t+1] = mode
        end
    end
    modes
end

end


text:
using Test
using LogicCircuits
import .LogicCircuits.LoadSave:
   parse_sdd_file, CircuitFormatLines, zoo_sdd_file

@testset "SDD file parser tests" begin

  @test parse_sdd_file(zoo_sdd_file("random.sdd")) isa CircuitFormatLines

  circuit = load_logic_circuit(zoo_sdd_file("random.sdd"))

  @test circuit isa PlainLogicCircuit
  @test num_nodes(circuit) == 1676
  @test isdecomposable(circuit)
  @test !issmooth(circuit)
  @test any(isfalse, linearize(circuit))
  @test any(istrue, linearize(circuit))
  @test num_variables(circuit) == 30
  @test issomething(canonical_literals(circuit))
  @test issomething(canonical_constants(circuit))

end


text:
using Calendar
perf = ref(String)
timing = ref(Float64)
function test()
	tic()
	for i = 1:1000000
		t = ymd_hms(1982,6,30,23,59,59)
	end
	return toq()
end
push!(perf,"Create Default CalendarTime")
push!(timing,test())
function test()
	tic()
	for i = 1:1000000
		t = ymd_hms(1982,6,30,23,59,59,"UTC")
	end
	return toq()
end
push!(perf,"Create UTC CalendarTime")
push!(timing,test())
function test()
	tic()
	for i = 1:1000000
		t = ymd_hms(1982,6,30,23,59,59,"VET")
	end
	return toq()
end
push!(perf,"Create VET CalendarTime")
push!(timing,test())
function test()
	tic()
	for i = 1:1000000
		t = ymd_hms(1982,6,30,23,59,59,"PST")
	end
	return toq()
end
push!(perf,"Create PST CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8)
	for i = 1:1000000
		r = dayofyear(t)
	end
	return toq()
end
push!(perf,"dayofyear CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8)
	for i = 1:1000000
		r = dayofweek(t)
	end
	return toq()
end
push!(perf,"dayofweek CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8)
	for i = 1:1000000
		r = week(t)
	end
	return toq()
end
push!(perf,"week CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = years(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ years(1) CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = months(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ months(1) CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = days(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ days(1) CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = hours(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ hours(1) CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = minutes(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ minutes(1) CalendarTime")
push!(timing,test())
function test()
	tic()
	t = ymd(2013,7,8,"UTC")
	tt = seconds(1)
	for i = 1:1000000
		r = t + tt
	end
	return toq()
end
push!(perf,"+ seconds(1) CalendarTime")
push!(timing,test())
#results
results = [perf timing]
println(results)


text:
# Indexing and dimensions (B.4)

export
    threadIdx, blockDim, blockIdx, gridDim,
    warpsize

@generated function _index(::Val{name}, ::Val{range}) where {name, range}
    JuliaContext() do ctx
        T_int32 = LLVM.Int32Type(ctx)

        # create function
        llvm_f, _ = create_function(T_int32)
        mod = LLVM.parent(llvm_f)

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry", ctx)
            position!(builder, entry)

            # call the indexing intrinsic
            intr_typ = LLVM.FunctionType(T_int32)
            intr = LLVM.Function(mod, "llvm.nvvm.read.ptx.sreg.$name", intr_typ)
            idx = call!(builder, intr)

            # attach range metadata
            range_metadata = MDNode([ConstantInt(Int32(range.start), ctx),
                                    ConstantInt(Int32(range.stop), ctx)],
                                    ctx)
            metadata(idx)[LLVM.MD_range] = range_metadata

            ret!(builder, idx)
        end

        call_function(llvm_f, UInt32)
    end
end

# XXX: these depend on the compute capability
#      https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities
const max_block_size = (x=1024, y=1024, z=64)
const max_grid_size  = (x=2^31-1, y=65535, z=65535)

for dim in (:x, :y, :z)
    # Thread index
    fn = Symbol("threadIdx_$dim")
    intr = Symbol("tid.$dim")
    @eval @inline $fn() = Int(_index($(Val(intr)), $(Val(0:max_block_size[dim]-1)))) + 1

    # Block size (#threads per block)
    fn = Symbol("blockDim_$dim")
    intr = Symbol("ntid.$dim")
    @eval @inline $fn() = Int(_index($(Val(intr)), $(Val(1:max_block_size[dim]))))

    # Block index
    fn = Symbol("blockIdx_$dim")
    intr = Symbol("ctaid.$dim")
    @eval @inline $fn() = Int(_index($(Val(intr)), $(Val(0:max_grid_size[dim]-1)))) + 1

    # Grid size (#blocks per grid)
    fn = Symbol("gridDim_$dim")
    intr = Symbol("nctaid.$dim")
    @eval @inline $fn() = Int(_index($(Val(intr)), $(Val(1:max_grid_size[dim]))))
end

"""
    gridDim()::CuDim3

Returns the dimensions of the grid.
"""
@inline gridDim() =   (x=gridDim_x(),   y=gridDim_y(),   z=gridDim_z())

"""
    blockIdx()::CuDim3

Returns the block index within the grid.
"""
@inline blockIdx() =  (x=blockIdx_x(),  y=blockIdx_y(),  z=blockIdx_z())

"""
    blockDim()::CuDim3

Returns the dimensions of the block.
"""
@inline blockDim() =  (x=blockDim_x(),  y=blockDim_y(),  z=blockDim_z())

"""
    threadIdx()::CuDim3

Returns the thread index within the block.
"""
@inline threadIdx() = (x=threadIdx_x(), y=threadIdx_y(), z=threadIdx_z())

"""
    warpsize()::UInt32

Returns the warp size (in threads).
"""
@inline warpsize() = Int(ccall("llvm.nvvm.read.ptx.sreg.warpsize", llvmcall, UInt32, ()))


text:
using ShrinkageEstimators
using Base.Test

d = 1
n = 10000
b = sqrt(n) # getfrequencies function uses sqrt(n) bins

arr = rand(n, d)
arrX = arr
arrY = rand(10000, 1)
arrZ = rand(10000, 1)

# Since the estimator is only an estimator, the tests are only approximate.
# In general a tolerance of 0.5 is used, which is relatively large. This
# should pass even if the estimate is slightly off, but fail if something
# is going drastically wrong.

# Check entropy is roughly log_base(b) for uniform distribution with b bins
@test_approx_eq_eps get_shrinkage_entropy(arr) log2(b) 0.5
println("Entropy passed")

# Check change of base works
@test_approx_eq_eps get_shrinkage_entropy(arr, base=e) log(b) 0.5
println("Entropy with change of base passed")

# Check specifying lambda works
@test_approx_eq_eps get_shrinkage_entropy(arr, lambda=0) log2(b) 0.5
println("Entropy with with lambda 0 passed")

# Check uniform width discretization works
@test_approx_eq_eps get_shrinkage_entropy(arr, mode="uniformcount") log2(b) 0.5
println("Entropy with with uniform count passed")

# TODO: Work out how to test "bayesianblocks" mode

# Check joint entropy between multiple identical distributions is roughly
# the same as entropy for one of them
@test_approx_eq_eps get_shrinkage_entropy(arr) get_shrinkage_entropy(arr, arr) 0.5
println("Joint entropy between 2 variables passed")
@test_approx_eq_eps get_shrinkage_entropy(arr) get_shrinkage_entropy(arr, arr, arr) 0.5
println("Joint entropy between 3 variables passed")

# Check I(X; Y) is H(X) + H(Y) - H(X, Y)
a = get_shrinkage_entropy(arrX) + get_shrinkage_entropy(arrY) - get_shrinkage_entropy(arrX, arrY)
@test a == get_shrinkage_mi(arrX, arrY)
println("Mutual information passed")

# Check I(X; Y | Z) is H(X, Z) + H(Y, Z) - H(X, Y, Z) - H(Z)
a = get_shrinkage_entropy(arrX, arrZ) + get_shrinkage_entropy(arrY, arrZ) - get_shrinkage_entropy(arrX, arrY, arrZ) - get_shrinkage_entropy(arrZ)
@test a == get_shrinkage_cmi(arrX, arrY, arrZ)
println("Conditional mutual information passed")


text:
function generatempo(model::UnitCellQModel;
                     symmetry::Symbol=:AUTO,
                     verbose::Bool=false)
    if model.symmetry == Trivial
        return _generatempo_nosym(model, verbose=verbose)

    elseif model.symmetry == U1
        return _generatempo_sym(model, verbose=verbose)

    elseif eltype(model.inters) <: FermionQModelInteraction
        error("JW not yet implemented! Do it manually and input the spin model for the mpo!")
    else
        error("Dont' recognize $(eltype(model.inters)) for MPO generation!")
    end
end

function _generatempo_nosym(model::UnitCellQModel;
                            verbose::Bool=false)
    typeof(model.qtype) == SpinType ||
        error("Can generate MPO only for Bosonic models.")
    D = dimension(model)
    model.lattice.bcs[D] == :OBC ||
        @warn "MPO generated periodic in X direction! discouraged!"

    d = model.qtype.d
    n_sites = prod(model.lattice.sizes)
    T = eltype(model.inters[1])

    sz,sp,sm = spinoperators(1/2)
    dict = Dict{Symbol, Matrix{T}}()
    dict[:Z] = sz
    dict[:P] = sp
    dict[:M] = sm
    dict[:I] = I(2)

    allterms = Vector{QTerm{Matrix{T}}}()
    for is in Iterators.product([1:l for l in model.lattice.sizes]...)
        for interaction in model.inters
            ns = interaction.ucidxs
            offs = interaction.offsets
            indexes = zeros(Int, support(interaction))
            isinside = true
            for i in 1:support(interaction)
                index, crossings = sitelinearindex(model.lattice, ns[i], offs[i] .+ is)
                if isnothing(index)
                    isinside = false
                    break
                end
                indexes[i] = index
            end
            if isinside
                perm = sortperm(indexes)
                verbose && println("Adding $(interaction.amp) between $(indexes[perm])")
                for qaterm in interaction.terms
                    ops = [dict[qaterm.ops[i]] for i in eachindex(qaterm.ops)]
                    permutedops = ops[perm]
                    termops = (permutedops[1:end-1]..., qaterm.amp * interaction.amp * permutedops[end])
                    push!(allterms, QTerm(Tuple(indexes[perm]), termops))
                end
            end
        end
    end
    sort!(allterms, by=x->x.sites[1])

    dims = ones(Int, n_sites+1)
    pointer=1
    mpo = MPO{T}()
    ldim = 1
    n_pcols = 0
    pterms = Vector{QTerm}()
    nextterms = Vector{QTerm}()
    for n in 1:n_sites
        lterms = Vector{QTerm}()
        sterms = Vector{QTerm}()
        n_ncols = 0
        while pointer <= length(allterms) && allterms[pointer].sites[1] == n
            if support(allterms[pointer]) > 1
                push!(sterms, allterms[pointer])
                n_ncols += length(allterms[pointer].ops)
            else
                push!(lterms, allterms[pointer])
            end
            pointer += 1
        end

        n_pcols = 0
        for pterm in pterms
            if pterm.sites[1] > n || support(pterm) > 1
                n_pcols += length(pterm.ops)
            end
        end

        if n < n_sites
            rdim = 2 + n_pcols + n_ncols
        else
            n_ncols == 0 || error()
            n_pcols == 0 || error()
            rdim = 1
        end

        dims[n+1] = rdim
        # MPO matrix generation
        W = zeros(T, ldim, d, rdim, d)

        # better logic is needed for this!
        if n > 1
            W[1,:,1,:] = I(d)
        end
        if n < n_sites
            W[ldim, :, rdim, :] = I(d)
        end

        row = 2
        col = 2
        for pterm in pterms
            if pterm.sites[1] == n
                if support(pterm) == 1
                    W[row, :, 1, :] = pterm.ops[1]
                        row+=1
                else
                        W[row, :, col, :] = pterm.ops[1]
                        row+=1
                        col+=1
                    push!(nextterms, removehead(pterm))
                end
            else
                    W[row, :, col, :] = I(d)
                    row+=1
                    col+=1
                push!(nextterms, pterm)
            end
        end
        n_pcols = col - 2

        for lterm in lterms
                W[ldim, :, 1, :] += lterm.ops[1]
        end

        for sterm in sterms
                W[ldim, :, col, :] = sterm.ops[1]
                col += 1
            push!(nextterms, removehead(sterm))
        end

        push!(mpo, W)
        pterms = nextterms
        nextterms = Vector{QTerm}()
        ldim = rdim
    end
    mpo
end

function _generatempo_sym(model::UnitCellQModel;
                          verbose::Bool=false)
    typeof(model.qtype) == SpinType ||
        error("Can generate MPO only for Bosonic models.")
    D = dimension(model)
    model.lattice.bcs[D] == :OBC ||
        @warn "MPO generated periodic in X direction! discouraged!"

    # typeof(model.inters[1]) <: SymQModelInteraction ||
    #     error("Only symmetric (U1) is allowed!")

    d = model.qtype.d
    d == 2 || error("only works for spin 1/2 for now!")
    n_sites = prod(model.lattice.sizes)
    T = eltype(model.inters[1])

    sz,sp,sm = spinoperators(1/2, symmetry=U1)
    dict = Dict{Symbol, SymMatrix{U1, T}}()
    dict[:Z] = sz
    dict[:P] = sp
    dict[:M] = sm
    dict[:I] = eye(T, U1Space(0=>1, 1=>1))

    allterms = Vector{QTerm{SymMatrix{U1,T}}}()
    for is in Iterators.product([1:l for l in model.lattice.sizes]...)
        for interaction in model.inters
            ns = interaction.ucidxs
            offs = interaction.offsets
            indexes = zeros(Int, support(interaction))
            isinside = true
            for i in 1:support(interaction)
                index, crossings = sitelinearindex(model.lattice, ns[i], offs[i] .+ is)
                if isnothing(index)
                    isinside=false
                    break
                end
                indexes[i] = index
            end
            if isinside
                perm = sortperm(indexes)
                for qaterm in interaction.terms
                    ops = [dict[qaterm.ops[i]] for i in eachindex(qaterm.ops)]
                    permutedops = ops[perm]
                    termops = (permutedops[1:end-1]..., qaterm.amp * interaction.amp * permutedops[end])
                    push!(allterms, QTerm(Tuple(indexes[perm]), termops))
                end
            end
        end
    end
    sort!(allterms, by=x->x.sites[1])

    pointer=1
    mpo = MPOperator{SymTensor{U1, T, 4}}()
    chrs = [zero(U1)]
    Vw = U1Space(0=>1)
    Vo = U1Space(0=>1, 1=>1)
    symI = eye(T, Vo)
    pchrs = Dict{U1Charge, Vector{QTerm}}()
    for n in 1:n_sites
        # find the new charge for each previous-term that still
        # continues to the next site
        nextchrs = U1[zero(U1)]
        for (pchr, pterms) in pchrs
            for pterm in pterms
                if pterm.sites[1] > n
                    push!(nextchrs, pchr)
                elseif support(pterm) > 1
                    push!(nextchrs, pchr + pterm.ops[1].charge)
                end
            end
        end

        # making local-terms and starting-terms and for starting terms
        # that continue find the charge the go to
        lterms = Vector{QTerm}()
        sterms = Vector{QTerm}()
        while pointer <= length(allterms) && allterms[pointer].sites[1] == n
            if support(allterms[pointer]) > 1
                push!(sterms, allterms[pointer])
                push!(nextchrs, allterms[pointer].ops[1].charge)
            else
                push!(lterms, allterms[pointer])
            end
            pointer += 1
        end

        sort!(nextchrs)
        chrs = U1[]
        chrdims = Int[]
        push!(chrs, nextchrs[1])
        push!(chrdims, 1)
        for c in nextchrs[2:end]
            if c == chrs[end]
                chrdims[end] += 1
            else
                push!(chrs, c)
                push!(chrdims, 1)
            end
        end

        nextterms = Dict(chr => Vector{QTerm}() for chr in chrs)
        # make the new leg
        if n < n_sites
            idx = searchsortedfirst(chrs, U1(0))
            chrdims[idx] += 1
            Vr = dual(U1Space([chrs[i]=>chrdims[i] for i in eachindex(chrs)]))
        else
            chrs == [U1(0)] || error("$chrs")
            chrdims == [1] || error("$chrdims")
            Vr = dual(U1Space(zero(U1)=>1))
        end

        W = fill(zero(T), U1(0), (Vw, Vo, Vr, dual(Vo)))

        l00, r00 = size(W.blocks[Sector{U1}(0,0,0,0)])[[1,3]]
        # better logic is needed for this!
        if n > 1
            W.blocks[Sector{U1}(0,0,0,0)][1,1,1,1] = one(T)
            W.blocks[Sector{U1}(0,1,0,1)][1,1,1,1] = one(T)
        end
        if n < n_sites
            W.blocks[Sector{U1}(0,0,0,0)][l00,1,r00,1] = one(T)
            W.blocks[Sector{U1}(0,1,0,1)][l00,1,r00,1] = one(T)
        end

        for lterm in lterms
            op = lterm.ops[1]
            op.charge == 0 || error()
            W.blocks[Sector{U1}(0,0,0,0)][l00, 1, 1, 1] += op[Sector{U1}(0,0)]
            W.blocks[Sector{U1}(0,1,0,1)][l00, 1, 1, 1] += op[Sector{U1}(1,1)]
        end

        rows = Dict(c=>1 for c in charges(Vw))
        rows[0] = 2
        cols = Dict(U1(c)=>1 for c in chrs)
        cols[0] = 2

        for (pchr, pterms) in pchrs
            for pterm in pterms
                if pterm.sites[1] == n
                    op = pterm.ops[1]
                    if support(pterm) == 1
                        # pchr goes to zero
                        row ,col = rows[pchr], cols[U1(0)]
                        for (s,blk) in op.blocks
                            c1, c2 = s.charges
                            sect = Sector{U1}(pchr, c1, 0, c2)
                            size(blk) == (1, 1) || error()
                            W.blocks[sect][row, 1, 1, 1] = blk[1,1]
                        end
                        rows[pchr] += 1
                    else
                        #pchr goes to pchr+op.charge
                        row ,col = rows[pchr], cols[pchr+op.charge]
                        for (s,blk) in op.blocks
                            c1, c2 = s.charges
                            sect = Sector{U1}(pchr, c1, pchr+op.charge, c2)
                            size(blk) == (1, 1) || error()
                            W.blocks[sect][row, 1, col, 1] = blk[1,1]
                        end
                        rows[pchr] += 1
                        cols[pchr+op.charge] += 1
                        push!(nextterms[pchr+op.charge], removehead(pterm))
                    end
                else
                    # pchr goes to pchr
                    row, col = rows[pchr], cols[pchr]
                    for (s,blk) in symI.blocks
                        c1, c2 = s.charges
                        sect = Sector{U1}(pchr, c1, pchr, c2)
                        size(blk) == (1, 1) || error()
                        W.blocks[sect][row, 1, col, 1] = blk[1,1]
                    end
                    rows[pchr] += 1
                    cols[pchr] += 1
                    push!(nextterms[pchr], pterm)
                end
            end
        end

        for sterm in sterms
            op = sterm.ops[1]
            row = rows[U1(0)]
            col = cols[op.charge]
            for (s,blk) in op.blocks
                c1, c2 = s.charges
                sect = Sector{U1}(0, c1, op.charge, c2)
                size(blk) == (1, 1) || error()
                W.blocks[sect][dim(Vw, U1(0)), 1, col, 1] = blk[1,1]
            end
            cols[op.charge] += 1
            push!(nextterms[op.charge], removehead(sterm))
        end

push!(mpo, W)
pchrs = nextterms
Vw = dual(Vr)
end
mpo
end

function _generatempo_infinite(model::UnitCellQModel;
                               verbose::Bool=false)
    typeof(model.qtype) == SpinType ||
        error("Can generate MPO only for Bosonic models!")
    D = dimension(model)
    lx = model.lattice.sizes[D]
    T = eltype(model.inters[1])
    model.lattice.bcs[D] in [:INF] ||
        error("boundary conitions inf not :INF!")

    typeof(model.inters[1]) <: SymQModelInteraction ||
        error("Only symmetric (U1) is allowed!")

    n_sites = prod(model.lattice.sizes[1:D-1])
    xrange = largestxrange(model)
    println(xrange)
    mpo = generatesymmpo(changeboundary(changesize(model, D, 2*(xrange+1)),
                                        D, :OBC))

    b = n_sites * xrange
    for i = 1:n_sites
        mpo.tensors[b+i] == mpo.tensors[b+i+n_sites] ||
            error("mpo not converged!")
    end

    tensors = Vector{SymTensor{T, 4}}()

    dims = Int[]
    for l = 1:lx
        for i=1:n_sites
            push!(dims, size(mpo.tensors[b+i], 1))
            push!(tensors, mpo.tensors[b+i])
        end
    end
    push!(dims, size(tensors[end], 3))
    return MPOperator(n_sites * lx, mpo.d, dims, tensors)
end
# function generateinfinitempo(model::UnitCellQModel)
#     typeof(model.qtype) == SpinType ||
#         error("Can generate MPO only for SpinType models.")
#     model.lattice.bc in [:PBCYAPBCX, :PBCX, :PBCYPBCX] ||
#         error("unrecognized boundary condition $boundary!")

#     typeof(model.inters[1]) <: SymQModelInteraction ||
#         error("Only symmetric (U1) is allowed!")

#
#     d = model.qtype.d
#     d == 2 || error("only works for spin 1/2 for now!")
#
#     T = eltype(model.inters[1])

#     ## notes: each term defined in the unit cell should be ordered by
#     ## indexes on the -inf to +inf index system. Then is should be
#     ## translated along the x-axis (last axis) so that the current
#     ## x-positions can be at any position from the beginnning to the
#     ## end, and all of those terms should be included in the infinte MPO
#     lattice = changeboundary(model.lattice, D, :INF)
#     allterms = Vector{QTerm{SymMatrix{T}}}()
#     for is in Iterators.product([1:l for l in model.lattice.sizes[1:D-1]]...)
#         for interaction in model.inters
#             ns = interaction.ucidxs
#             offs = interaction.offsets
#             indexes = zeros(Int, support(interaction))
#             isinside = true
#             for i in 1:support(interaction)
#                 index, crossings = sitelinearindex(lattice, ns[i], offs[i] .+ is)
#                 if isnothing(index)
#                     isinside = false
#                     break
#                 end
#                 indexes[i] = index
#             end
#             if isinside
#             perm = sortperm(indexes)
#             permindexes = indexes[perm]
#             if support(interaction) == 1
#                 for ops in interaction.terms
#                     n = mod1(indexes[1], n_sites)
#                     push!(lterms[n], interaction.amp * ops[1])
#                 end
#                 break
#             end
#             for ops in interaction.terms
#                 permutedops = ops[perm]
#                 termops = (permutedops[1:end-1]..., interaction.amp * permutedops[end])
#                 pointer = 1
#                 n = mod1(permindexes[pointer], n_sites)
#                 push!(sterms[n], termops[pointer])
#                 index = permindexes[pointer] + 1
#                 pointer += 1
#                 while pointer < support(interation)
#                     if index == permindexes[pointer]
#                         pointer += 1
#                     else
#                         push!(cterms[n], I)
#                     end
#                     index += 1
#                 end
#                 n = mod1(permindexes[pointer], n_sites)
#                 push!(eterms[n], termops[pointer])

#                         if n == permindexes[pointer]
#                             push!(cterms[n], termops[pointer])
#                             pointer +=1
#                         else
#                             # Add I
#                         end
#                     end
#                 end
#             end
#         end
#         end
#     end
#     sort!(allterms, by=x->x.sites[1])

#     dims = ones(Int, n_sites+1)
#     pointer=1
#     tensors = Vector{SymTensor{Float64, 4}}(undef, n_sites)
#     chrs = [0]
#     lleg = STLeg(+1, chrs, [1])
#     o1leg = STLeg(+1, [0,1], [1,1])
#     o2leg = STLeg(-1, [0,1], [1,1])
#     symI = eye(T, collect(0:d-1), ones(Int, d))
#     pchrs = Dict{Int, Vector{QTerm}}()
#     for n in 1:n_sites
#         # find the new charge for each previous-term that still
#         # continues to the next site
#         nextchrs = Int[0]
#         for (pchr, pterms) in pchrs
#             for pterm in pterms
#                 if pterm.sites[1] > n
#                     push!(nextchrs, pchr)
#                 elseif support(pterm) > 1
#                     push!(nextchrs, pchr + pterm.ops[1].charge)
#                 end
#             end
#         end

#         # making local-terms and starting-terms and for starting terms
#         # that continue find the charge the go to
#         lterms = Vector{QTerm}()
#         sterms = Vector{QTerm}()
#         while pointer <= length(allterms) && allterms[pointer].sites[1] == n
#             if support(allterms[pointer]) > 1
#                 push!(sterms, allterms[pointer])
#                 push!(nextchrs, allterms[pointer].ops[1].charge)
#             else
#                 push!(lterms, allterms[pointer])
#             end
#             pointer += 1
#         end

#         sort!(nextchrs)
#         chrs = Int[]
#         chrdims = Int[]
#         push!(chrs, nextchrs[1])
#         push!(chrdims, 1)
#         for c in nextchrs[2:end]
#             if c == chrs[end]
#                 chrdims[end] += 1
#             else
#                 push!(chrs, c)
#                 push!(chrdims, 1)
#             end
#         end

#         nextterms = Dict(chr => Vector{QTerm}() for chr in chrs)
#         # make the new leg
#         if n < n_sites
#             idx = searchsortedfirst(chrs, 0)
#             chrdims[idx] += 1
#             rleg = STLeg(-1, chrs, chrdims)
#         else
#             chrs == [0] || error("$chrs")
#             chrdims == [1] || error("$chrdims")
#             rleg = STLeg(-1, [0], [1])
#         end

#         dims[n+1] = fulldims(rleg)

#         W = fill(zero(T), 0, (lleg, o1leg, rleg, o2leg))

#         index0000 = index_sector(W, (0,0,0,0))
#         index0101 = index_sector(W, (0,1,0,1))
#         l00, r00 = size(W.blocks[index0000])[[1,3]]
#         # better logic is needed for this!
#         if n > 1
#             W.blocks[index0000][1,1,1,1] = one(T)
#             W.blocks[index0101][1,1,1,1] = one(T)
#         end
#         if n < n_sites
#             W.blocks[index0000][l00,1,r00,1] = one(T)
#             W.blocks[index0101][l00,1,r00,1] = one(T)
#         end

#         for lterm in lterms
#             op = lterm.ops[1]
#             op.charge == 0 || error()
#             W.blocks[index0000][l00, 1, 1, 1] += get_sector(op, (0,0))
#             W.blocks[index0101][l00, 1, 1, 1] += get_sector(op, (1,1))
#         end

#         rows = Dict(c=>1 for c in lleg.chrs)
#         rows[0] = 2
#         cols = Dict(c=>1 for c in chrs)
#         cols[0] = 2

#         for (pchr, pterms) in pchrs
#             for pterm in pterms
#                 if pterm.sites[1] == n
#                     op = pterm.ops[1]
#                     if support(pterm) == 1
#                         # pchr goes to zero
#                         row ,col = rows[pchr], cols[0]
#                         for i in 1:length(op.sects)
#                             c1, c2 = op.sects[i]
#                             idx = index_sector(W, (pchr, c1, 0, c2))
#                             size(op.blocks[i]) == (1, 1) || error()
#                             W.blocks[idx][row, 1, 1, 1] = op.blocks[i][1,1]
#                         end
#                         rows[pchr] += 1
#                     else
#                         #pchr goes to pchr+op.charge
#                         row ,col = rows[pchr], cols[pchr+op.charge]
#                         for i in 1:length(op.sects)
#                             c1, c2 = op.sects[i]
#                             idx = index_sector(W, (pchr, c1, pchr+op.charge, c2))
#                             size(op.blocks[i]) == (1, 1) || error()
#                             W.blocks[idx][row, 1, col, 1] = op.blocks[i][1,1]
#                         end
#                         rows[pchr] += 1
#                         cols[pchr+op.charge] += 1
#                         push!(nextterms[pchr+op.charge], removehead(pterm))
#                     end
#                 else
#                     # pchr goes to pchr
#                     row, col = rows[pchr], cols[pchr]
#                     for i in 1:length(symI.sects)
#                         c1, c2 = symI.sects[i]
#                         idx = index_sector(W, (pchr, c1, pchr, c2))
#                         size(symI.blocks[i]) == (1, 1) || error()
#                         W.blocks[idx][row, 1, col, 1] = symI.blocks[i][1,1]
#                     end
#                     rows[pchr] += 1
#                     cols[pchr] += 1
#                     push!(nextterms[pchr], pterm)
#                 end
#             end
#         end

#         for sterm in sterms
#             op = sterm.ops[1]
#             row, col = rows[0], cols[op.charge]
#             for i in 1:length(op.sects)
#                 c1, c2 = op.sects[i]
#                 idx = index_sector(W, (0, c1, op.charge, c2))
#                 size(op.blocks[i]) == (1, 1) || error()
#                 W.blocks[idx][getdim(lleg, 0), 1, col, 1] = op.blocks[i][1,1]
#             end
#             cols[op.charge] += 1
#             push!(nextterms[op.charge], removehead(sterm))
#         end

# tensors[n] =  W
# pchrs = nextterms
# lleg = STLeg(+1, rleg.chrs, rleg.dims)
# end
# MPOperator(n_sites, d, dims, tensors)
# end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
"""
```
AnSchorfheide{T} <: AbstractRepModel{T}
```

The `AnSchorfheide` type defines the structure of the simple New Keynesian DSGE
model described in 'Bayesian Estimation of DSGE Models' by Sungbae An and Frank
Schorfheide.

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::OrderedDict{Symbol,Int}`: Maps lagged states to
  their columns in the measurement and equilibrium condition equations. These
  are added after `gensys` solves the model.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"an_schorfheide\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct AnSchorfheide{T} <: AbstractRepModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                       # model steady-state values
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
                                                           # parameters and steady-states

    endogenous_states::OrderedDict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::OrderedDict{Symbol,Int}              # measurement and equilibrium condition equations.
    expected_shocks::OrderedDict{Symbol,Int}               #
    equilibrium_conditions::OrderedDict{Symbol,Int}        #
    endogenous_states_augmented::OrderedDict{Symbol,Int}   #
    observables::OrderedDict{Symbol,Int}                   #
    pseudo_observables::OrderedDict{Symbol,Int}            #

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::AnSchorfheide) = "Julia implementation of model defined in 'Bayesian Estimation of DSGE Models' by Sungbae An and Frank Schorfheide: AnSchorfheide, $(m.subspec)"

"""
`init_model_indices!(m::AnSchorfheide)`

Arguments:
`m:: AnSchorfheide`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::AnSchorfheide)
    # Endogenous states
    endogenous_states = collect([
        :y_t, :π_t, :R_t, :y_t1, :g_t, :z_t, :Ey_t1, :Eπ_t1])

    # Exogenous shocks
    exogenous_shocks = collect([
        :z_sh, :g_sh, :rm_sh])

    # Expectations shocks
    expected_shocks = collect([
        :Ey_sh, :Eπ_sh])

    # Equilibrium conditions
    equilibrium_conditions = collect([
        :eq_euler, :eq_phillips, :eq_mp, :eq_y_t1, :eq_g, :eq_z, :eq_Ey, :eq_Eπ])

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = []

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(exogenous_shocks);            m.exogenous_shocks[k]            = i end
    for (i,k) in enumerate(expected_shocks);             m.expected_shocks[k]             = i end
    for (i,k) in enumerate(equilibrium_conditions);      m.equilibrium_conditions[k]      = i end
    for (i,k) in enumerate(endogenous_states);           m.endogenous_states[k]           = i end
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(endogenous_states) end
    for (i,k) in enumerate(observables);                 m.observables[k]                 = i end
    for (i,k) in enumerate(pseudo_observables);          m.pseudo_observables[k]          = i end
end

function AnSchorfheide(subspec::String="ss0";
                       custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                       testing = false)

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = AnSchorfheide{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), OrderedDict{Symbol,Int}(),

            # model indices
            OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(), OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Set settings
    model_settings!(m)
    default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable and pseudo-observable transformations
    init_observable_mappings!(m)
    init_pseudo_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)

    return m
end

"""
```
init_parameters!(m::AnSchorfheide)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::AnSchorfheide)
    # Initialize parameters
    m <= parameter(:τ, 1.9937, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), GammaAlt(2., 0.5), fixed=false,
                   description="τ: The inverse of the intemporal elasticity of substitution.",
                   tex_label="\\tau")

    m <= parameter(:κ, 0.7306, (1e-20, 1-1e-7), (1e-20, 1-1e-7), ModelConstructors.SquareRoot(), Uniform(0,1), fixed=false,
                   description="κ: Composite parameter in New Keynesian Phillips Curve.",
                   tex_label="\\kappa")

    m <= parameter(:ψ_1, 1.1434, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), GammaAlt(1.5, 0.25), fixed=false,
                   description="ψ_1: The weight on inflation in the monetary policy rule.",
                   tex_label="\\psi_1")
    m <= parameter(:ψ_2, 0.4536, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), GammaAlt(0.5, 0.25), fixed=false,
                   description="ψ_2: The weight on the output gap in the monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:rA, 0.0313, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), GammaAlt(0.5, 0.5), fixed=false,
                   description="rA: β (discount factor) = 1/(1+ rA/400).",
                   tex_label="rA")

    m <= parameter(:π_star, 8.1508, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), GammaAlt(7., 2.), fixed=false,
                   description="π_star: Target inflation rate.",
                   tex_label="\\pi*")

    m <= parameter(:γ_Q, 1.5, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), Normal(0.40, 0.20), fixed=false,

                   description="γ_Q: Steady state growth rate of technology.",
                   tex_label="\\gamma_Q")

    m <= parameter(:ρ_R, 0.3847, (1e-20, 1-1e-7), (1e-20, 1-1e-7), ModelConstructors.SquareRoot(), Uniform(0,1), fixed=false,
                   description="ρ_R: AR(1) coefficient on interest rate.",
                   tex_label="\\rho_R")

    m <= parameter(:ρ_g, 0.3777, (1e-20, 1-1e-7), (1e-20, 1-1e-7), ModelConstructors.SquareRoot(), Uniform(0,1), fixed=false,
                   description="ρ_g: AR(1) coefficient on g_t = 1/(1 - ζ_t), where ζ_t is government spending as a fraction of output.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_z, 0.9579, (1e-20, 1-1e-7), (1e-20, 1-1e-7), ModelConstructors.SquareRoot(), Uniform(0,1), fixed=false,
                   description="ρ_z: AR(1) coefficient on shocks to the technology growth rate.",
                   tex_label="\\rho_z")

    m <= parameter(:σ_R, 0.4900, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), RootInverseGamma(4, .4), fixed=false,
                   description="σ_R: Standard deviation of shocks to the nominal interest rate.",
                   tex_label="\\sigma_R")

    m <= parameter(:σ_g, 1.4594, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), RootInverseGamma(4, 1.), fixed=false,
                   description="σ_g: Standard deviation of shocks to the government spending process.",
                   tex_label="\\sigma_g")

    m <= parameter(:σ_z, 0.9247, (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(), RootInverseGamma(4, 0.5), fixed=false,
                   description="σ_z: Standard deviation of shocks to the technology growth rate process.",
                   tex_label="\\sigma_z")

    m <= parameter(:e_y, 0.20*0.579923, fixed=true,
                   description="e_y: Measurement error on GDP growth.",
                   tex_label="e_y")

    m <= parameter(:e_π, 0.20*1.470832, fixed=true,
                   description="e_π: Measurement error on inflation.",
                   tex_label="e_\\pi")

    m <= parameter(:e_R, 0.20*2.237937, fixed=true,
                   description="e_R: Measurement error on the interest rate.",
                   tex_label="e_R")
end

"""
```
steadystate!(m::AnSchorfheide)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
the parameters of `m` are updated.
"""
function steadystate!(m::AnSchorfheide)
    return m
end

function model_settings!(m::AnSchorfheide)
    default_settings!(m)

    # Data
    m <= Setting(:data_id, 0, "Dataset identifier")
    m <= Setting(:cond_full_names, [:obs_gdp, :obs_nominalrate],
        "Observables used in conditional forecasts")
    m <= Setting(:cond_semi_names, [:obs_nominalrate],
        "Observables used in semiconditional forecasts")

    # Metropolis-Hastings
    m <= Setting(:mh_cc, 0.27,
                 "Jump size for Metropolis-Hastings (after initialization)")

    # Forecast
    m <= Setting(:use_population_forecast, true,
                 "Whether to use population forecasts as data")
    m <= Setting(:forecast_zlb_value, 0.13,
        "Value of the zero lower bound in forecast periods, if we choose to enforce it")
end

function shock_groupings(m::AnSchorfheide)
    gov = ShockGroup("g", [:g_sh], RGB(0.70, 0.13, 0.13)) # firebrick
    tfp = ShockGroup("z", [:z_sh], RGB(1.0, 0.55, 0.0)) # darkorange
    pol = ShockGroup("pol", vcat([:rm_sh], [Symbol("rm_shl$i") for i = 1:n_anticipated_shocks(m)]),
                     RGB(1.0, 0.84, 0.0)) # gold
    det = ShockGroup("dt", [:dettrend], :gray40)

    return [gov, tfp, pol, det]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RunCommandResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{InstanceViewStatus} } # spec name: value

    function RunCommandResult(;value=nothing)
        o = new()
        validate_property(RunCommandResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type RunCommandResult

const _property_map_RunCommandResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_RunCommandResult = Dict{Symbol,String}(Symbol("value")=>"Vector{InstanceViewStatus}")
Base.propertynames(::Type{ RunCommandResult }) = collect(keys(_property_map_RunCommandResult))
Swagger.property_type(::Type{ RunCommandResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RunCommandResult[name]))}
Swagger.field_name(::Type{ RunCommandResult }, property_name::Symbol) =  _property_map_RunCommandResult[property_name]

function check_required(o::RunCommandResult)
    true
end

function validate_property(::Type{ RunCommandResult }, name::Symbol, val)
end


text:
# You cannot start or import this file directly, use runtest.jl instead
using InteractiveUtils
include("../src/SimpleCache.jl")

@noinline function multi(x::Int)
    println("printing stuff 1")
    @debug "debug macro 1"
    return VALUE_A
end

@noinline function multi(x)
    println("printing stuff 2")
    @debug "debug macro 2"
    return VALUE_B
end

@noinline function caller(x)
    multi(x)
    return 5
end

global glob = VALUE_C
@noinline function use_global()
    return glob
end

const constglob = VALUE_D
@noinline function use_const()
    return constglob
end

@noinline function ya_func(x)
    return VALUE_E
end

@noinline function pointer(x)
    if rand() > 0.5
        ptr = multi
    else
        ptr = ya_func
    end
    return ptr(x)
end

println(SimpleCache.hash_function_recursive([], @code_typed(caller(4)), true))
println(SimpleCache.hash_function_recursive([], @code_typed(caller(4.0)), true))
println(SimpleCache.hash_function_recursive([], @code_typed(use_global()), true))
println(SimpleCache.hash_function_recursive([], @code_typed(use_const()), true))
println(SimpleCache.hash_function_recursive([], @code_typed(pointer(4)), true))


text:
function perform_stft(x,w,q,n)
    """
        perform_stft - compute a local Fourier transform

        Forward transform:
          MF = perform_stft(M,w,q, options);
        Backward transform:
          M  = perform_stft(MF,w,q, options);

          w is the width of the window used to perform local computation.
          q is the spacing betwen each window.

          MF(:,i) contains the spectrum around point (i-1)*q

          A typical use, for an redundancy of 2 could be w=2*q+1

          options.bound can be either 'per' or 'sym'

          No multichannel, no Gabor stft

          options.normalization can be set to
              'tightframe': tight frame transform, with energy conservation.
              'unit': unit norm basis vectors, usefull to do thresholding

          Copyright (c) 2008 Gabriel Peyre
    """
    if ndims(x) == 1
        dir = 1
    else
        dir = -1
    end

    # perform sampling
    X = Array(1:q:n+2)

    p = length(X)
    eta = 1

    if w%2 == 1
        w = ceil((w-1)/2)*2+1
        w1 = Base.div(w-1,2)
        dX = Array(-w1:w1)
    else
        dX = Array((-Base.div(w,2)+1):Base.div(w,2))
    end

    X1 = repeat(X,outer=(1,w)) + repeat(dX, outer=(1,p))'
    #periodic boundary conditions
    X1 = ((X1-1).%n)+1;

    I = X1'
    I[I.<=0]+=n

    # build a sin weight function
    W = .5 *(1 - cos(2*pi*Array(0:(w-1))./(w-1)))

    #renormalize the windows
    weight = zeros(n)

    for i in 1:p
        weight[I[:,i]] = weight[I[:,i]] + W.^2
    end

    weight = sqrt(weight)
    Weight = repeat(W, outer=(1,p))

    for i in 1:p
        Weight[:,i] = Weight[:,i]./weight[I[:,i]]
    end

    #compute the transform
    if dir == 1
        y = zeros(eta*w,p)
        if w%2 == 1
            m = Base.div(eta*w+1,2)
            w1 = Base.div((w-1),2)
            sel = Array(m-w1:m+w1) - 1
        else
            m = Base.div((eta*w),2)+1
            w1 = Base.div(w,2)
            sel = Array(m-w1:m+w1-1) - 1
        sel[sel.<=0] += eta*w
        y[sel,:] = x[I].*Weight

        #perform the transform
        y = my_transform(y,+1)
        end
    else
        x = my_transform(x,-1)
        x = real.(x.*Weight)
        y = zeros(n)
        for i in 1:p
            y[I[:,i]] = y[I[:,i]] + x[:,i]
        end
    end
    return y
end
########################################################################
function my_transform(x,dir)

    # my_transform - perform either FFT with energy conservation.
    # Works on array of size (w,w,a,b) on the 2 first dimensions.
    w = size(x)[1]
    if dir == 1
        y = fft(x,1)/sqrt(w)
    else
        y = ifft(x*sqrt(w),1)
    end
    return y
end


text:
struct TimeInvariantFullTimeGF{T, scalar} <: AbstractTimeGF{T, scalar}
  grid::FullTimeGrid
  gtr::AntiHermitianToeplitzStorage{T,scalar}
  les::AntiHermitianToeplitzStorage{T,scalar}
  rm::GenericStorage{T,scalar}
  mat::PeriodicStorage{T,scalar}
  ξ::GFSignEnum
end

norbitals(G::TimeInvariantFullTimeGF) = G.gtr.norb

function TimeInvariantFullTimeGF(::Type{T}, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  nt = grid.nt
  ntau = grid.ntau

  gtr = AntiHermitianToeplitzStorage(T, nt, norb, scalar)
  les = AntiHermitianToeplitzStorage(T, nt, norb, scalar)
  rm = GenericStorage(T, ntau, nt, norb, scalar)
  mat = PeriodicStorage(T, ntau, norb, scalar)

  TimeInvariantFullTimeGF(grid, gtr, les, rm, mat, ξ)
end
TimeInvariantFullTimeGF(grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = TimeInvariantFullTimeGF(ComplexF64, grid, norb, ξ, scalar)

@inline function Base.getindex(G::TimeInvariantFullTimeGF, k, l, t1::TimeGridPoint, t2::TimeGridPoint, greater=true)
  greater = t1 == t2 ? greater : heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx
  ξ = Int(G.ξ)

  real_branches = (forward_branch, backward_branch)
  if (t1.bpoint.domain ∈ real_branches && t2.bpoint.domain ∈ real_branches)
    return greater ? G.gtr[k,l,i,j] : G.les[k,l,i,j]
  elseif (t1.bpoint.domain == imaginary_branch && (t2.bpoint.domain == forward_branch || t2.bpoint.domain == backward_branch))
    return G.rm[k,l,i,j]
  elseif ((t1.bpoint.domain == forward_branch || t1.bpoint.domain == backward_branch) && t2.bpoint.domain == imaginary_branch)
    ntau = G.grid.ntau
    return copy(-ξ * adjoint(G.rm[l,k,ntau+1-j,i])) # Aoki 19c
  else
    greater ? G.mat[k,l,i,j,greater] : ξ * G.mat[k,l,i,j,greater]
  end
end

function Base.setindex!(G::TimeInvariantFullTimeGF, v, k, l, t1::TimeGridPoint, t2::TimeGridPoint)
  greater = heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx
  ξ = Int(G.ξ)

  real_branches = (forward_branch, backward_branch)
  if (t1.bpoint.domain ∈ real_branches && t2.bpoint.domain ∈ real_branches)
    return greater ? G.gtr[k,l,i,j] = v : G.les[k,l,i,j] = v
  elseif (t1.bpoint.domain == imaginary_branch && (t2.bpoint.domain == forward_branch || t2.bpoint.domain == backward_branch))
    return G.rm[k,l,i,j] = v
  elseif ((t1.bpoint.domain == forward_branch || t1.bpoint.domain == backward_branch) && t2.bpoint.domain == imaginary_branch)
    ntau = G.grid.ntau
    return G.rm[l,k,ntau+1-j,i] = -ξ * adjoint(v) # Aoki 19c
  else
    if greater
      G.mat[k,l,i,j] = v
    else
      G.mat[k,l,i,j] = ξ * v
    end
  end
end

function Base.similar(G::T, ξ::GFSignEnum) where T <: TimeInvariantFullTimeGF
  T(G.grid, similar(G.gtr), similar(G.les), similar(G.rm), similar(G.mat), ξ)
end

function Base.zero(G::T, ξ::GFSignEnum) where T <: TimeInvariantFullTimeGF
  T(G.grid, zero(G.gtr), zero(G.les), zero(G.rm), zero(G.mat), ξ)
end

function Base.similar(G::T) where T <: TimeInvariantFullTimeGF
  similar(G, G.ξ)
end

function Base.zero(G::T) where T <: TimeInvariantFullTimeGF
  zero(G, G.ξ)
end

function TimeDomain(G::TimeInvariantFullTimeGF)
  grid = G.grid
  tplus = grid[forward_branch]
  tminus = reverse(grid[backward_branch])
  tau = grid[imaginary_branch]

  nt = length(tplus)
  ntau = length(tau)

  points = Tuple{TimeGridPoint, TimeGridPoint}[]

  for i in 1:nt
    push!(points, (tminus[1], tminus[i])) # greater
    push!(points, (tplus[1], tminus[i])) # lesser
  end

  for i in 1:nt
    for j in 1:ntau
      push!(points, (tau[j], tminus[i])) # right mixing
    end
  end

  for i in 1:ntau
    push!(points, (tau[i], tau[1])) # matsubara
  end

  return TimeDomain(points)
end

function TimeInvariantFullTimeGF(f::Function, ::Type{T}, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  G = TimeInvariantFullTimeGF(T, grid, norb, ξ, scalar)
  D = TimeDomain(G)

  for (t1, t2) in D.points
    G[t1, t2] = f(t1, t2)
  end

  return G
end

TimeInvariantFullTimeGF(f::Function, grid::FullTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = TimeInvariantFullTimeGF(f, ComplexF64, grid, norb, ξ, scalar)

function TimeInvariantFullTimeGF(dos::AbstractDOS, grid::FullTimeGrid)
  β = grid.contour.β
  TimeInvariantFullTimeGF(grid, 1, fermionic, true) do t1, t2
    Keldysh.dos2gf(dos, β, t1.bpoint, t2.bpoint)
  end
end