text:
using IntervalArithmetic
import Base: getindex, ∪
import Base: +, -, *, /, min, max, ^, log, <, >, exp, sin, cos, tan
import IntervalArithmetic: hull

"""

    Multi-Intervals are unions of disjoint intervals.
    This file includes constructors, arithmetic (including intervals and scalars)
    and complement functions

    Empty sets and intersecting intervals are appropriately handled in the constructor:

    julia> a = interval(0,2) ∪ interval(3,4)
    [0, 2] ∪ [3, 4]

    julia> b = interval(1,2) ∪ interval(4,5) ∪ ∅
    [1, 2] ∪ [4, 5]

    julia> c = a * b
    [0, 10] ∪ [12, 20]

    julia> complement(c)
    [-∞, 0] ∪ [10, 12] ∪ [20, ∞]

"""

abstract type MultiInterval{T} <: AbstractInterval{T} end

###
#   Multi-Interval constructor. Consists of a vector of intervals
###
struct IntervalM{T<:Real} <: MultiInterval{T}
    v :: Array{Interval{T}}
end

###
#   Outer constructors
###
function intervalM(x)
    x = IntervalM(x)
    x = removeEmpties(x)
    return condense(x)
end

intervalM(x :: Interval) = IntervalM([x])
∪(x :: Interval) = intervalM(x)

∪(x :: Interval, y :: Interval) = intervalM([x; y])
∪(x :: Array{Interval{T}}) where T <:Real = intervalM(x)

intervalM(x :: Interval, y :: IntervalM) = intervalM([x; y.v])
∪(x :: Interval, y :: IntervalM) = intervalM(x,y)

intervalM(x :: IntervalM, y :: Interval) = intervalM([x.v; y])
∪(x :: IntervalM, y :: Interval) = intervalM(x,y)

intervalM(x :: IntervalM, y :: IntervalM) = intervalM([x.v; y.v])
∪(x :: IntervalM, y :: IntervalM) = intervalM(x,y)

# MultiInterval can act like a vector
getindex(x :: IntervalM, ind :: Integer) = getindex(x.v,ind)
getindex(x :: IntervalM, ind :: Array{ <: Integer}) = getindex(x.v,ind)

# Remove ∅ from Multi-Interval
function removeEmpties(x :: IntervalM)
    v = x.v
    Vnew = v[v .!= ∅]
    return IntervalM(Vnew)
end

# Envolpe intervals which intersect. Recursive condense!
function condense(x :: IntervalM)

    if isCondensed(x); return x; end

    v = sort(x.v)
    v = unique(v)

    Vnew = Interval{Float64}[]
    for i =1:length(v)
        these = findall( intersect.(v[i],v ) .!= ∅)
        push!(Vnew, hull(v[these]))
    end
    return condense( intervalM(Vnew) )
end

function isCondensed(x :: IntervalM)
    v = sort(x.v)
    for i=1:length(v)
        intersects = findall( intersect.(v[i],v[1:end .!= i]) .!= ∅)
        if !isempty(intersects); return false; end
    end
    return true
end

# Envolpe/hull. Keeps it as a multi interval
env(x :: IntervalM) = IntervalM([hull(x.v)])
hull(x :: IntervalM) = IntervalM([hull(x.v)])

# Computes the complement of a Multi-Interval
function complement(x :: IntervalM)

    v = sort(x.v)
    vLo = left.(v)
    vHi = right.(v)

    pushfirst!(vHi, -∞)
    push!(vLo, ∞)

    complements = interval.(vHi,vLo)

    return intervalM(complements)
end

complement(x :: Interval) = complement(intervalM(x))

###
#   Multi-Interval Arithmetic
###

for op in (:+, :-, :*, :/, :min, :max, :^, :log, :<, :>)
    @eval ($op)( x::IntervalM, y::IntervalM) = intervalM([$op(xv, yv) for xv in x.v, yv in y.v][:])

    @eval ($op)( x::IntervalM, y::Interval) = intervalM(broadcast($op, x.v, y))
    @eval ($op)( x::Interval, y::IntervalM) = intervalM(broadcast($op, y, x.v))

    @eval ($op)( x::IntervalM, n::Real) = intervalM(broadcast($op, x.v, n))
    @eval ($op)( n::Real, x::IntervalM) = intervalM(broadcast($op, n, x.v))

end

for op in (:-, :sin, :cos, :tan, :exp, :log)
    @eval ($op)( x::IntervalM) = intervalM(broadcast($op, x.v))
end

###
#   Utilities and show
###

left(x :: Interval) = x.lo
right(x :: Interval) = x.hi

left(x :: IntervalM) = left.(x.v)
right(x :: IntervalM) = right.(x.v)

function Base.show(io::IO, this::IntervalM)
    v = sort(this.v)
    print(io, join(v, " ∪ "));
end


text:
using LinearAlgebra
using Distributions

function kinetic_matrix(num_sites, t)   # one-body matrix, periodic boundary condition
    K_matrix = zeros(num_sites, num_sites)
    for i =  1 : num_sites
        for j = 1 : num_sites
            if abs(i - j) == 1 || abs(i - j) == num_sites - 1
                K_matrix[i, j] = -t
            end
        end
    end
    if flag_background_subtraction == 1
        K_matrix += Diagonal(background_dist) * U
    end
    return K_matrix
end

function auxiliary_field_matrix()
    ϕ = rand(Normal(), num_sites)
    coefficent_vector = ϕ * AF_coefficient
    AF_matrix = complex(zeros(num_sites, num_sites))
    for i = 1 : num_sites
        AF_matrix[i, i] = coefficent_vector[i]
    end
    if flag_background_subtraction == 1
        background_exponent = sum(-coefficent_vector .* background_dist)
        background_weight = exp(background_exponent)
        return background_weight * exp(AF_matrix)
    elseif flag_background_subtraction == 0
        return exp(AF_matrix)
    end
end

ext: jl

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
function schedule_by_day(config)
    # "config" is dict with keys start_day and method
    start_day = config["start_day"]
    method = config["method"]
    employees = config["employees"]
    env = config["env"]
    perfect_optimality = week_sum_coverage(env)

    if (method == "balanced")
            ok, weekly_schedule = schedule_by_day_balanced(employees, env, start_day)
    elseif (method == "greedy")
        ok, weekly_schedule = schedule_by_day_greedy(employees, env, start_day)
    else
        warn("Unknown method $method")
        return {"success" => false}
    end
    if ok
        valid_schedule, message = validate_schedule(employees, env, weekly_schedule)
        if !valid_schedule
            Logging.info( "Unexpected validation issue in $method scheduler - $message")
            return {"success" => false}
        end
        hours_scheduled = week_hours_scheduled(weekly_schedule)
        overage = (hours_scheduled - perfect_optimality) / perfect_optimality
        Logging.info( "$method day scheduler succeeded - overage $overage ")

        return {
            "success" => true,
            "schedule" => weekly_schedule,
        }
    end
    return {"success" => false}
end

function schedule_by_day_balanced(employees, env, start_day=1)
    days_per_week = size(env["coverage"], 1)
    hours_scheduled = Dict()
    weekly_schedule = Dict()
    for e in keys(employees)
        hours_scheduled[e] = 0
        weekly_schedule[e] = Any[]
    end

    end_day = (start_day + days_per_week - 2) % days_per_week + 1

    for day_index in start_day:(days_per_week + start_day - 1)
        # stupid 1-indexing
        # Need to wrap around day so we can support different start days
        day = (day_index - 1) % days_per_week+ 1

        day_env = Dict()
        day_env["coverage"] = env["coverage"][day]
        day_env["cycle_length"] = length(env["coverage"][day])

        day_employees = Dict()
        for e in keys(employees)
            if employees[e]["days_assigned"][day] != 1
                continue
            end

            day_min, day_max = get_day_bounds(
                sub_array(employees[e]["longest_availability"], day, end_day),
                sub_array(employees[e]["days_assigned"], day, end_day),
                hours_scheduled[e],
                employees[e]["hours_min"],
                employees[e]["hours_max"],
                employees[e]["shift_time_min"],
            )

            day_employees[e] = Dict()
            day_employees[e]["min"] = day_min
            day_employees[e]["max"] = day_max
            day_employees[e]["availability"] = employees[e]["availability"][day]
        end

        ok, day_schedule = schedule_day(day_employees, day_env)
        if !ok
            return false, {}
        end

        for e in keys(day_schedule)
            push!(weekly_schedule[e], {
                "day" => day,
                "start" => day_schedule[e]["start"],
                "length" => day_schedule[e]["length"],
            })
            hours_scheduled[e] += day_schedule[e]["length"]
        end

        ## INTERSHIFT
        # Adjust next day's availability
        if day == days_per_week
            # only do it if there is another day
            continue
        end

        for e in keys(day_schedule)
            # Note: this assumes that intershift < day_length
            # TODO - add a check in verify function for this
            next_day = day + 1
            if next_day > days_per_week
                next_day = (next_day - 1) % days_per_week + 1
            end

            if employees[e]["days_assigned"][next_day] != 1
                continue
            end
            # last time slot where the person works
            last_shift = day_schedule[e]["start"] + day_schedule[e]["length"] - 1
            hours_left_in_day = length(env["coverage"][day]) - last_shift

            if "intershift" in keys(employees[e])
                intershift = employees[e]["intershift"]
            else
                intershift = env["intershift"]
            end

            overflow = intershift - hours_left_in_day - env["time_between_coverage"]

            if overflow <= 0
                continue
            end

            # adjust availability
            for t in 1:overflow
                employees[e]["availability"][next_day][t] = 0
                employees[e]["availability"][next_day] = anneal_availability(employees[e]["availability"][next_day], employees[e]["shift_time_min"])
            end
            employees[e]["longest_availability"][next_day] = get_longest_availability(employees[e], next_day)
        end

    end

    return true, weekly_schedule
end

#=
This scheduler does not account for future availability of employees.
Instead, when it reaches an infeasibility it tries to *drop* employee
shifts. This may behoove our algorithm because it tends to assign people
to the upper bound of their shift limits.
=#
function schedule_by_day_greedy(employees, env, start_day=1)
    days_per_week = size(env["coverage"], 1)
    hours_scheduled = Dict()
    weekly_schedule = Dict()
    shifts_scheduled = Dict()
    for e in keys(employees)
        hours_scheduled[e] = 0
        weekly_schedule[e] = Any[]
        shifts_scheduled[e] = sum(employees[e]["days_assigned"])
    end

    end_day = (start_day + days_per_week - 2) % days_per_week + 1

    for day_index in start_day:(days_per_week + start_day - 1)
        # stupid 1-indexing
        # Need to wrap around day so we can support different start days
        day = (day_index - 1) % days_per_week+ 1

        day_env = Dict()
        day_env["coverage"] = env["coverage"][day]
        day_env["cycle_length"] = length(env["coverage"][day])

        day_employees = Dict()
        for e in keys(employees)
            if employees[e]["days_assigned"][day] != 1
                continue
            end

            shift_time_min = employees[e]["shift_time_min"]
            shift_time_max = employees[e]["shift_time_max"]
            hours_min = employees[e]["hours_min"]
            hours_max = employees[e]["hours_max"]

            if (
                # Worked enough hours already
                hours_scheduled[e] >= hours_min &&
                hours_scheduled[e] < hours_max
                )
                # Either adjust the shift
                remaining_hours = hours_max - hours_scheduled[e]
                if (remaining_hours >= shift_time_min)
                    # only adjust if necessary
                    if (remaining_hours < shift_time_max)
                        shift_time_max = remaining_hours
                    end
                    # otherwise don't touch
                else # Try dropping  the shift
                    if shifts_scheduled[e] > employees[e]["shift_time_min"]
                        shifts_scheduled[e] -= 1
                        continue # will this work?
                    else
                        # We're fucked
                        return false, {}
                    end
                end
            end

            day_employees[e] = Dict()
            day_employees[e]["min"] = shift_time_min
            day_employees[e]["max"] = shift_time_max
            day_employees[e]["availability"] = employees[e]["availability"][day]
        end

        ok, day_schedule = schedule_day(day_employees, day_env)
        if !ok
            return false, {}
        end

        for e in keys(day_schedule)
            push!(weekly_schedule[e], {
                "day" => day,
                "start" => day_schedule[e]["start"],
                "length" => day_schedule[e]["length"],
            })
            hours_scheduled[e] += day_schedule[e]["length"]
        end

        ## INTERSHIFT
        # Adjust next day's availability
        if day == days_per_week
            # only do it if there is another day
            continue
        end

        for e in keys(day_schedule)
            # Note: this assumes that intershift < day_length
            # TODO - add a check in verify function for this
            next_day = day + 1
            if next_day > days_per_week
                next_day = (next_day - 1) % days_per_week + 1
            end

            if employees[e]["days_assigned"][next_day] != 1
                continue
            end
            # last time slot where the person works
            last_shift = day_schedule[e]["start"] + day_schedule[e]["length"] - 1
            hours_left_in_day = length(env["coverage"][day]) - last_shift

            if "intershift" in keys(employees[e])
                intershift = employees[e]["intershift"]
            else
                intershift = env["intershift"]
            end

            overflow = intershift - hours_left_in_day - env["time_between_coverage"]

            if overflow <= 0
                continue
            end

            # adjust availability
            for t in 1:overflow
                employees[e]["availability"][next_day][t] = 0
                employees[e]["availability"][next_day] = anneal_availability(employees[e]["availability"][next_day], employees[e]["shift_time_min"])
            end
            employees[e]["longest_availability"][next_day] = get_longest_availability(employees[e], next_day)
        end

    end

    return true, weekly_schedule
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
module Passes

include("SyntaxHighlighter.jl")
include("BracketHighlighter.jl")
include("RainbowBrackets.jl")
end # module


text:
#=
This is a template for a development/debugging script.
=#

# Wrap things up in a module so we can reload stuff as needed.
module M

using Pseudospectra, Test

# To run like a test:
# displayplots = true # If you want to see results immediately.
include("../test/testplotter.jl")

# Alternatively, to run like a user:
#=
using PyPlot
setpsplotter(:PyPlot)
gs = setgs()
=#

A = # something new and challenging

opts = Dict{Symbol,Any}(:verbosity => 2,
                        # whatever
                        )
ps_data = new_matrix(A,opts)
driver!(ps_data,opts,gs)

end # module


text:
# Solves a few basic eigenvalue problems to test various aspects of NLEIGS

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using NonlinearEigenproblems.RKHelper
using Test
using LinearAlgebra

function nleigs_basic()
    n = 2
    B = Vector{Matrix{Float64}}([[1 3; 5 6], [3 4; 6 6], [1 0; 0 1]])
    pep = PEP(B)

    Σ = [-10.0-2im, 10-2im, 10+2im, -10+2im]

    @bench @testset "Polynomial only" begin
        lambda, X = nleigs(pep, Σ, maxit=10, v=ones(n).+0im, blksize=5)
        verify_lambdas(4, pep, lambda, X)
    end

    @bench @testset "Non-convergent linearization" begin
        @test_logs (:warn, r".*Linearization not converged.*") match_mode = :any begin
            lambda, X = nleigs(pep, Σ, maxit=10, v=ones(n).+0im, maxdgr=5, blksize=5)
            verify_lambdas(4, pep, lambda, X)
        end
    end

    @bench @testset "Non-convergent linearization (static)" begin
        @test_logs (:warn, r".*Linearization not converged.*") match_mode = :any begin
            lambda, X = nleigs(pep, Σ, maxit=10, v=ones(n).+0im, maxdgr=5, blksize=5, static=true)
            verify_lambdas(4, pep, lambda, X)
        end
    end

    @bench @testset "Non-convergent linearization (return_details)" begin
        @test_logs (:warn, r".*Linearization not converged.*") match_mode = :any begin
            lambda, X, _ = nleigs(pep, Σ, maxit=5, v=ones(n).+0im, blksize=5, return_details=true)
            verify_lambdas(0, pep, lambda, X)
        end
    end

    @bench @testset "Complex-valued matrices" begin
        complex_B = map(X -> X + im*I, B)
        complex_pep = PEP(complex_B)
        lambda, X, _ = nleigs(complex_pep, Σ, maxit=10, v=ones(n).+0im, blksize=5, return_details=true)
        verify_lambdas(3, complex_pep, lambda, X)
    end

    @bench @testset "Complex-valued start vector" begin
        lambda, X, _ = nleigs(pep, Σ, maxit=10, v=ones(n) * (1+0.1im), blksize=5, return_details=true)
        verify_lambdas(4, pep, lambda, X)
    end

    @bench @testset "return_details" begin
        lambda, X, res, details = nleigs(pep, Σ, maxit=10, v=ones(n).+0im, blksize=5, return_details=true)
        verify_lambdas(4, pep, lambda, X)

        lam = details.Lam[:,end]
        res = details.Res[:,end]
        conv = map(i -> res[i] < 1e-12 && inpolygon(real(lam[i]), imag(lam[i]), real(Σ), imag(Σ)), 1:size(lam, 1))
        lamconv = lam[conv, end]

        # test that eigenvalues in the info are the same as those returned by nleigs
        @test length(lamconv) == 4
        @test length(union(lambda, lamconv)) == 4
    end
end

@testset "NLEIGS: Basic functionality" begin
    nleigs_basic()
end


text:
module Unwinder

using DWARF
using DWARF.CallFrameInfo
import DWARF.CallFrameInfo: realize_cie, RegStates, CIE
using ObjFileBase
using ObjFileBase: handle
using Gallium
using ..Registers
using ..Registers: ip
using ObjFileBase
using ObjFileBase: Sections, mangle_sname
using ELF
using MachO
using COFF
using Gallium: find_module, Module, load, make_fdetab, make_inverse_symtab
import Gallium: symbolicate

typealias CFICacheEntry Tuple{CallFrameInfo.RegStates,CallFrameInfo.CIE,UInt}
type CFICache
    sz :: Int
    values :: Dict{RemotePtr{Void}, CFICacheEntry}
end

CFICache(sz::Int) = CFICache(sz, Dict{RemotePtr{Void}, Tuple{CallFrameInfo.RegStates,CallFrameInfo.CIE,UInt}}())

function get_word(s::Gallium.LocalSession, ptr::RemotePtr)
    Gallium.Hooking.mem_validate(UInt(ptr), sizeof(Ptr{Void})) || error("Invalid load")
    UInt64(load(s, RemotePtr{Gallium.intptr(Gallium.getarch(s))}(ptr)))
end
get_word(s, ptr::RemotePtr) = UInt64(load(s, RemotePtr{Gallium.intptr(Gallium.getarch(s))}(ptr)))

function find_fde(mod, modrel)
    slide = 0
    eh_frame = Gallium.find_ehframes(mod)[1]
    if isa(mod, Module) && isnull(mod.ehfr)
        isnull(mod.FDETab) &&
            (mod.FDETab = Nullable(make_fdetab(mod.base, mod)))
        return CallFrameInfo.search_fde_offset(eh_frame, get(mod.FDETab),
            modrel, slide, is_eh_not_debug = mod.is_eh_not_debug)
    else
        tab = Gallium.find_ehfr(mod)
        modrel = Int(modrel)-Int(sectionoffset(tab.hdr_sec))
        slide = Int(sectionoffset(tab.hdr_sec)) - Int(sectionoffset(eh_frame))
        loc, fde = CallFrameInfo.search_fde_offset(eh_frame, tab, modrel, slide)
        loc = loc + Int(sectionoffset(tab.hdr_sec))
        return (loc, fde)
    end
end

function probably_in_entrypoint(h, ip)
    start = Gallium.X86_64.BasicRegs(deref(first(filter(x->isa(deref(x),MachO.thread_command),LoadCmds(h))))).rip[]
    start -= Gallium.first_actual_segment(h).vmaddr
    # Default OS X _start is 63 bytes
    start <= ip <= start+63
end

function entry_cfa(mod, r)
    rs = DWARF.CallFrameInfo.RegStates()
    regs = Gallium.X86_64.inverse_dwarf
    rs.cfa = DWARF.CallFrameInfo.Offset(
        DWARF.CallFrameInfo.RegNum(regs[:rbp]), 0)
    rs[regs[:rip]] = DWARF.CallFrameInfo.Load(DWARF.CallFrameInfo.RegCFA, 0)
    rs, DWARF.CallFrameInfo.CIE(0,0,0,regs[:rip],UInt8[]), UInt64(0)
end

function modulerel(mod, base, ip)
    ret = (ip - base)
end

get_ciecache(_) = nothing
function get_ciecache(mod::Module)
    if isnull(mod.ciecache) && !isnull(mod.eh_frame)
        CallFrameInfo.precompute(get(mod.eh_frame), mod.is_eh_not_debug)
    end
    isnull(mod.ciecache) && return nothing
    get(mod.ciecache)
end

realize_cie(mod::Module, fde) = realize_cie(get_ciecache(mod), fde)

function compute_register_states(s, base, mod, r, stacktop, ::Void)
    modrel = UInt(modulerel(mod, base, UInt(ip(r))))
    loc, fde = try
        find_fde(mod, modrel)
    catch e
        # As a special case, if we're in a MachO executable's entry point,
        # we probably don't have unwind info. TODO: Remove this once we support
        # compact unwind info which the entry point does have.
        if isa(handle(mod), MachO.MachOHandle) && readheader(handle(mod)).filetype == MachO.MH_EXECUTE
            probably_in_entrypoint(handle(mod), modrel) && return entry_cfa(mod, r)
        end
        rethrow(e)
    end
    ciecache = get_ciecache(mod)
    cie::CIE, ccoff = realize_cieoff(fde, ciecache)
    # Compute CFA
    target_delta::UInt64 = modrel - loc - (stacktop?0:1)
    @assert target_delta < UInt(CallFrameInfo.fde_range(fde, cie))
    #out = STDOUT #IOContext(STDOUT, :reg_map => Gallium.X86_64.dwarf_numbering)
    #drs = CallFrameInfo.RegStates()
    #CallFrameInfo.dump_program(out, cie, target = UInt(target_delta), rs = drs); println(out)
    #CallFrameInfo.dump_program(out, fde, cie = cie, target = UInt(target_delta), rs = drs)
    CallFrameInfo.evaluate_program(fde, UInt(target_delta), cie = cie, ciecache = ciecache, ccoff=ccoff)::RegStates, cie, target_delta
end

function compute_register_states(s, base, mod, r, stacktop, cfi_cache::CFICache)
    pc = RemotePtr{Void}(UInt(ip(r)))
    if haskey(cfi_cache.values, pc)
        return cfi_cache.values[pc]::CFICacheEntry
    else
        result = compute_register_states(s, base, mod, r, stacktop, nothing)::CFICacheEntry
        if length(cfi_cache.values) < cfi_cache.sz
            cfi_cache.values[pc] = result
        end
        result::CFICacheEntry
    end
end

immutable Frame
    cfa_addr::RemotePtr{Void}
    rs::RegStates
    cie::CIE
    target_delta::UInt64
end

function evaluate_cfi_expr(opcodes, s, r, rs)
    sm = DWARF.Expressions.StateMachine{UInt64}() # typeof(unsigned(ip(r)))
    getreg(reg) = get_dwarf(r, reg)
    getword(addr) = get_word(s, RemotePtr{UInt64}(addr))[]
    addr_func(addr) = addr
    loc = DWARF.Expressions.evaluate_simple_location(sm, opcodes, getreg, getword, addr_func, :NativeEndian)
    if isa(loc, DWARF.Expressions.RegisterLocation)
        addr = RemotePtr{Void}(get_dwarf(r, loc.i))
    else
        addr = RemotePtr{Void}(loc.i)
    end
    addr
end
evaluate_cfa_expr(s, r, rs) = evaluate_cfi_expr(rs.cfa_expr.opcodes, s, r, rs)

function compute_cfa_addr(s, r, rs)
    local cfa_addr
    if !CallFrameInfo.isdwarfexpr(rs.cfa)
        if rs.cfa.flag != CallFrameInfo.Flag.Val
            error("invalid CFA value $(rs.cfa)")
        end
        cfa_addr = RemotePtr{Void}(convert(Int, get_dwarf(r, Int(rs.cfa.base)) + rs.cfa.offset))
    else
        cfa_addr = evaluate_cfa_expr(s, r, rs)
    end
end

function frame(s, base, mod, r, stacktop :: Bool, cfi_cache)
    rs, cie, target_delta = compute_register_states(s, base, mod, r, stacktop, cfi_cache)::CFICacheEntry
    Frame(compute_cfa_addr(s, r, rs), rs, cie, UInt64(target_delta))
end

immutable TransformedArray{A,T,F} <: AbstractArray{T,1}
    arr::A
    func::F
end
TransformedArray{T,F}(a::AbstractArray{T,1}, func::F) =
    TransformedArray{typeof(a), T, F}(a, func)
Base.length(a::TransformedArray) = length(a.arr)
Base.size(a::TransformedArray) = size(a.arr)
Base.getindex(a::TransformedArray, idx) = a.func(a.arr[idx])

symbolicate(modules, ip) = symbolicate(Gallium.LocalSession(), modules, ip)
function symbolicate(session, modules, ip)
    base, mod = try
        find_module(session, modules, ip)
    catch err
        if !isa(err, ErrorException) || !contains(err.msg, "not found")
            rethrow(err)
        end
        return (false, "<Unknown Module>")
    end
    symbolicate(session, base, mod, ip)
end
function symbolicate(session, base, mod, ip)
    modrel = UInt(modulerel(mod, base, ip))
    loc = modrel
    approximate = true
    try
        if !isnull(mod.eh_frame)
            loc, fde = find_fde(mod, modrel)
            approximate = false
        elseif !isnull(mod.xpdata)
            loc = find_seh_entry(mod, modrel).start
            approximate = false
        end
    end
    sections = Sections(Gallium.dhandle(mod))
    syms = Gallium.get_syms(mod)
    function correct_symbol(x)
        isundef(x) && return (false, UInt64(0))
        !isa(handle(mod), COFF.COFFHandle) || COFF.isfunction(x) || return (false, UInt64(0))
        isa(handle(mod), ELF.ELFHandle) &&
            (ELF.st_type(x) != ELF.STT_FUNC) && return (false, UInt64(0))
        value = symbolvalue(x, sections)
        #@show value
        (true, value)
    end
    if isa(mod, Gallium.Module)
        isnull(mod.inverse_symtab) &&
            (mod.inverse_symtab = Nullable(make_inverse_symtab(Gallium.dhandle(mod))))
        idx = searchsortedfirst(TransformedArray(get(mod.inverse_symtab),
            idx->symbolvalue(syms[idx], sections)), loc)
        if !approximate
            while idx <= length(syms)
                ok, value = correct_symbol(syms[get(mod.inverse_symtab)[idx]])
                (mod.is_jit_dobj) && (value -= base)
                (ok && value == loc) && break
                ok && value > loc && (#=idx = 0;=# break)
                idx += 1
            end
        end
        (idx == length(get(mod.inverse_symtab))+1) && (idx = 0)
        idx != 0 && (idx = get(mod.inverse_symtab)[idx])
    else
        @assert !approximate
        idx = findfirst(syms) do sym
            ok, value = correct_symbol(sym)
            ok && value == loc
        end
    end
    idx == 0 && return "???"
    name = symname(syms[idx]; strtab = StrTab(syms))
    (!approximate, name)
end

function fetch_cfi_val_value(s, r, resolution, cfa_addr)
    if resolution.base == CallFrameInfo.RegCFA
        return (convert(UInt64,cfa_addr)%Int64 + resolution.offset) % UInt64
    else
        return convert(UInt64, get_dwarf(r, Int(resolution.base))%Int64 + resolution.offset)
    end
end

function fetch_cfi_value(s, r, rs, reg, cfa_addr)
    resolution = rs.values[reg]
    if CallFrameInfo.isdwarfexpr(resolution)
        ve = rs.values_expr[reg]
        val = evaluate_cfi_expr(ve.opcodes, s, r, rs)
        ve.is_val || (val = get_word(s,RemotePtr{UInt64}(val)))
        return val
    elseif CallFrameInfo.issame(resolution)
        return get_dwarf(r, reg)
    elseif resolution.flag == CallFrameInfo.Flag.Val#isa(resolution, CallFrameInfo.Offset)
        return fetch_cfi_val_value(s, r, resolution, cfa_addr)
    elseif resolution.flag == CallFrameInfo.Flag.Deref
        new_resolution = CallFrameInfo.RegState(resolution.base, resolution.offset, CallFrameInfo.Flag.Val)
        addr = RemotePtr{Void}(fetch_cfi_val_value(s, r, new_resolution, cfa_addr))
        return get_word(s, addr)
    else
        error("Unknown resolution $resolution")
    end

end

function unwind_step_frame_pointer!(new_registers, s)
    # We don't really know anything about how to unwind here, but let's
    # try frame pointer based unwinding and hope FPO isn't enabled
    # TODO: Maybe try assembly profiling
    if isa(getarch(s),X86_64.X86_64Arch)
        old_rbp = get_dwarf(new_registers, :rbp)
        set_dwarf!(new_registers, :rbp, get_word(s, RemotePtr{UInt64}(old_rbp)))
        set_ip!(new_registers, get_word(s, RemotePtr{UInt64}(old_rbp+8)))
    else
        old_ebp = get_dwarf(new_registers, :ebp)
        set_dwarf!(new_registers, :ebp, get_word(s, RemotePtr{UInt32}(old_ebp)))
        set_ip!(new_registers, get_word(s, RemotePtr{UInt32}(old_ebp+4)))
    end
    new_registers
end

using Gallium: X86_64
function unwind_step(s, modules, r, cfi_cache = nothing; stacktop = false, ip_only = false, allow_frame_based = true)
    new_registers = copy(r)
    # A priori the registers in the new frame will not be valid, we copy them
    # over from above still and propagate as usual in case somebody wants to
    # look at them.
    invalidate_regs!(new_registers)

    # First, find the module we're currently in
    base, mod = find_module(s, modules, UInt(ip(r)))
    modrel = UInt64(ip(r)) - base

    # Determine if we have windows or DWARF unwind info
    if !isnull(mod.eh_frame)
        cf = try
            frame(s, base, mod, r, stacktop, cfi_cache)
        catch e
            allow_frame_based &&
                return (true, unwind_step_frame_pointer!(new_registers, s))
            rethrow(e)
            return (false, r)
        end

        # By definition, the next frame's stack pointer is our CFA
        set_sp!(new_registers, UInt(cf.cfa_addr))
        CallFrameInfo.isundef(cf.rs[cf.cie.return_reg]) && return (false, r)
        # Find current frame's return address, (i.e. the new frame's ip)
        set_ip!(new_registers, fetch_cfi_value(s, r, cf.rs, cf.cie.return_reg, cf.cfa_addr))
        # Now set other registers recorded in the CFI
        if !ip_only
            for reg in keys(cf.rs.values)
                resolution = cf.rs.values[reg]
                reg == cf.cie.return_reg && continue
                set_dwarf!(new_registers, reg, fetch_cfi_value(s, r, cf.rs, reg, cf.cfa_addr))
            end
        end
    elseif !isnull(mod.xpdata)
        entry = find_seh_entry(mod, modrel)
        offs = modrel - entry.start
        frameaddr = 0
        if entry.info.FrameRegOff != 0
            reg = entry.info.FrameRegOff & 0xf
            regoff = (entry.info.FrameRegOff & 0xf0) >> 4
            frameaddr = get_dwarf(r, X86_64.seh_numbering[reg]) -
                16*regoff
        end
        set_dwarf!(new_registers, :rsp, get_dwarf(r, :rsp))
        i = 1
        codes = reinterpret(UNWIND_CODE, entry.opcodes)
        while i <= length(entry.opcodes)
            op = codes[i]
            run = op.PrologueOffset <= offs
            old_rsp = get_dwarf(new_registers, :rsp)
            opcode, info = Gallium.opinfo(op)
            if opcode == UWOP_PUSH_NONVOL
                if run
                    # Pop the register
                    set_dwarf!(new_registers, X86_64.seh_numbering[info],
                        get_word(s, RemotePtr{UInt64}(old_rsp)))
                    set_dwarf!(new_registers, :rsp, old_rsp+8)
                end
                i += 1
            elseif opcode == UWOP_ALLOC_LARGE || opcode == UWOP_ALLOC_SMALL
                # Undo the allocation
                size, nskip = Gallium.compute_alloc_op_size(entry.opcodes, i)
                run && set_dwarf!(new_registers, :rsp, old_rsp + size)
                i += nskip
            elseif opcode == UWOP_SAVE_XMM128
                # TODO
                i += 2
            elseif opcode == UWOP_SET_FPREG
                set_dwarf!(new_registers, :rsp, frameaddr)
                i += 1
            else
                error("Unknown operation $(op.Operation&0xf)")
            end
        end
        # The return address should now be on the stack. Pop it as well
        old_rsp = get_dwarf(new_registers, :rsp)
        set_ip!(new_registers, get_word(s, RemotePtr{UInt64}(old_rsp)))
        set_dwarf!(new_registers, :rsp, old_rsp+8)
    else
        allow_frame_based &&
            return (true, unwind_step_frame_pointer!(new_registers, s))
        error("Ununwindable module")
    end
    UInt(ip(new_registers)) == 0 &&  return (false, r)
    (true, new_registers)
end

# Win64 unwinding
using Gallium: PData, RUNTIME_FUNCTION, UNWIND_CODE,
    UWOP_PUSH_NONVOL, UWOP_ALLOC_LARGE, UWOP_ALLOC_SMALL,
    UWOP_SET_FPREG, UWOP_SAVE_NONVOL, UWOP_SAVE_NONVOL_FAR,
    UWOP_SAVE_XMM128, UWOP_SAVE_XMM128_FAR, UWOP_PUSH_MACHFRAME

function find_seh_entry(mod, modrel)
    slide = 0
    xpref = get(mod.xpdata)
    idx = searchsortedlast(PData(xpref.pdata), RUNTIME_FUNCTION(modrel,0,0), by = x->Int64(x.startoff))
    return xpref[idx]
end

end


text:
import Pkg
Pkg.add("Documenter")
using Documenter
using CellListMap
DocMeta.setdocmeta!(CellListMap, :DocTestSetup, :(using CellListMap); recursive=true)
push!(LOAD_PATH,"../src/")
makedocs(
    modules=[CellListMap],
    sitename="CellListMap.jl",
    pages = [
        "Overview" => "index.md",
        "Examples" => "examples.md",
        "Periodic conditions" => "pbc.md",
        "Parallelization" => "parallelization.md",
        "Units, autodiff, etc." => "units_etc.md",
        "Performance" => "performance.md",
        "Options" => "options.md",
        "From Python" => "python.md",
        "Help entries" => "help.md",
        "Reference" => "reference.md",
    ]
)
deploydocs(
    repo = "github.com/m3g/CellListMap.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "main",
    versions = ["stable" => "v^", "v#.#" ],
)


text:
@testset "devices" begin

name(dev)
totalmem(dev)
attribute(dev, CUDAdrv.DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK)
@test warpsize(dev) == Int32(32)
capability(dev)
@grab_output show(stdout, "text/plain", dev)

@test eval(Meta.parse(repr(dev))) == dev

@test eltype(devices()) == CuDevice
@grab_output show(stdout, "text/plain", CUDAdrv.DEVICE_CPU)
@grab_output show(stdout, "text/plain", CUDAdrv.DEVICE_INVALID)
end


text:
using DSGE, ModelConstructors, Dates, CSV, DataFrames, Plots, Test, HDF5
include("tvcred_parameterize.jl")

regenerate_reference_forecasts = false

# Forecast settings
fcast_date = DSGE.quartertodate("2020-Q4")
date_fcast_end  = iterate_quarters(fcast_date, 60)

# Flexible AIT Rule
imperfect_cred_flexait = .33
imperfect_cred_taylor  = 1. - imperfect_cred_flexait
φ_π                    = 10.
φ_y                    = 6.
Thalf                  = 10.
ρ_smooth               = 0.7
pgap_init              = 0.125
ygap_init              = 12.
pgap_ygap_init_date    = Date(2020, 6, 30)
start_zlb_date         = Date(2021, 3, 31)
end_zlb_date           = Date(2022, 12, 31)

# additional settings to implement Flexible AIT rule
flexait_custom = Dict{Symbol, Setting}(:add_initialize_pgap_ygap_pseudoobs => Setting(:add_initialize_pgap_ygap_pseudoobs, true),
                                       :add_pgap => Setting(:add_pgap, true), :add_ygap => Setting(:add_ygap, true),
                                       :add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true),
                                       :add_altpolicy_ygap => Setting(:add_altpolicy_ygap, true),
                                       :forecast_horizons => Setting(:forecast_horizons,
                                                                     subtract_quarters(date_fcast_end, fcast_date)),
                                       :forecast_smoother => Setting(:forecast_smoother, :durbin_koopman),
                                       :contemporaneous_and_proportional_antshocks =>
                                       Setting(:contemporaneous_and_proportional_antshocks, Symbol[:biidc]),
                                       :antshocks =>
                                       Setting(:antshocks, Dict{Symbol, Int}(:biidc => 1,
                                                                             :φ => 1, :ziid => 1)),
                                       :ant_eq_mapping =>
                                       Setting(:ant_eq_mapping, Dict{Symbol, Symbol}(:biidc => :biidc,
                                                                                :φ => :φ,
                                                                                :ziid => :ziid)),
                                       :ant_eq_E_mapping =>
                                       Setting(:ant_eq_E_mapping, Dict{Symbol, Symbol}(:φ => :Eφ)),
                                       :proportional_antshocks =>
                                       Setting(:proportional_antshocks, [:biidc, :φ, :ziid]),
                                       :n_anticipated_obs_gdp => Setting(:n_anticipated_obs_gdp, 1),
                                       :add_anticipated_obs_gdp => Setting(:add_anticipated_obs_gdp, true),
                                        :meas_err_anticipated_obs_gdp =>
                                        Setting(:meas_err_anticipated_obs_gdp, 1.),
                                       :flexible_ait_policy_change =>
                                       Setting(:flexible_ait_policy_change, false))

# Initialize model object
m = Model1002("ss59"; custom_settings = flexait_custom)
usual_model_settings!(m, "200001", cdvt = "200001", fcast_date = fcast_date)
m <= Setting(:time_varying_trends, true)
get_setting(m, :regime_dates)[5] = Date(2020, 12, 31)
setup_regime_switching_inds!(m, cond_type = :full)

m10 = Model1002("ss10") # for help initializing parameters of m
θ10 = h5read(joinpath(dirname(@__FILE__), "..", "reference", "tvcred_reference_forecast.h5"), "para")
DSGE.update!(m10, θ10)
update_vals = zeros(Float64, length(m.parameters))
for k in map(x-> x.key, m10.parameters)
    update_vals[m.keys[k]] = m10[k].value
end
DSGE.update!(m, update_vals) # make sure parameters match Model 1002 ss10, w/zeros for other parameters

# Parameterize
tvcred_parameterize!(m)

# Load data
df = CSV.read(joinpath(dirname(@__FILE__), "..", "reference", "timevaryingcred_df.csv"), DataFrame)
df[!, :obs_pgap] .= NaN # NaN out just to be safe
df[!, :obs_ygap] .= NaN
ind_init = findfirst(df[!, :date] .== pgap_ygap_init_date)
df[end, :obs_gdp] = df[end - 1, :obs_gdp1] # To match reference forecast

# Set up Flexible AIT with a temporary ZLB. Note that df will be altered
θ = (φ_π = φ_π, φ_y = φ_y, Thalf = Thalf, pgap = pgap_init,
     ygap = ygap_init, ρ_smooth = ρ_smooth,
     cred = imperfect_cred_flexait,
     historical_policy = DSGE.taylor_rule()) # NamedTuple with parameters of the Flexible AIT policy
df[ind_init, :obs_pgap] = -θ[:pgap] # initialize pgap and ygap in the date
df[ind_init, :obs_ygap] = -θ[:ygap] # pgap_ygap_init_date, default to 2020:Q2

m <= Setting(:flexible_ait_φ_π, θ[:φ_π])
m <= Setting(:flexible_ait_φ_y, θ[:φ_y])
m <= Setting(:ait_Thalf, θ[:Thalf])
m <= Setting(:gdp_Thalf, θ[:Thalf])
m <= Setting(:pgap_value, θ[:pgap])
m <= Setting(:pgap_type, :flexible_ait)
m <= Setting(:ygap_value, θ[:ygap])
m <= Setting(:ygap_type, :flexible_ait)
m <= Setting(:flexible_ait_ρ_smooth, θ[:ρ_smooth])
m <= Setting(:alternative_policies, AltPolicy[θ[:historical_policy]])
m <= Setting(:skip_altpolicy_state_init, true)

## Set up temporary ZLB
m <= Setting(:gensys2, true)
m <= Setting(:uncertain_temporary_altpolicy, true)
m <= Setting(:uncertain_altpolicy, true)
m <= Setting(:temporary_altpolicy_names, [:zero_rate])

# Set regime_dates and regime_eqcond_info for temp ZLB to flexible AIT
gensys2_first_regime = get_setting(m, :n_regimes) + 1
get_setting(m, :regime_dates)[gensys2_first_regime] = start_zlb_date
n_zlb_reg = DSGE.subtract_quarters(end_zlb_date, start_zlb_date) + 1

# Set up regime_eqcond_info
m <= Setting(:replace_eqcond, true)
reg_dates = deepcopy(get_setting(m, :regime_dates))
regime_eqcond_info = Dict{Int, DSGE.EqcondEntry}()
for (regind, date) in zip(gensys2_first_regime:(n_zlb_reg - 1 + gensys2_first_regime), # See comments starting at line 57
                          DSGE.quarter_range(reg_dates[gensys2_first_regime],
                                             DSGE.iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg - 1)))
    reg_dates[regind] = date
    regime_eqcond_info[regind] = DSGE.EqcondEntry(DSGE.zero_rate(), [0.5, 0.5])
end
reg_dates[n_zlb_reg + gensys2_first_regime] = DSGE.iterate_quarters(reg_dates[gensys2_first_regime], n_zlb_reg)
regime_eqcond_info[n_zlb_reg + gensys2_first_regime] = DSGE.EqcondEntry(DSGE.flexible_ait(), [0.5, 0.5])
nreg0 = length(reg_dates)
m <= Setting(:regime_dates,             reg_dates)
m <= Setting(:regime_eqcond_info, regime_eqcond_info)
setup_regime_switching_inds!(m; cond_type = :full)
set_regime_vals_fnct(m, get_setting(m, :n_regimes))
m <= Setting(:temporary_altpolicy_length, n_zlb_reg)
m <= Setting(:zero_rate_zlb_value, 0.)
# Now add additional regimes of flexible AIT to allow time-varying credibility
for i in nreg0:(nreg0 + 15)
    reg_dates[i] = DSGE.iterate_quarters(reg_dates[nreg0], i - nreg0)
    regime_eqcond_info[i] = DSGE.EqcondEntry(DSGE.flexible_ait(), [0.5, 0.5])
end
m <= Setting(:regime_dates,             reg_dates)
m <= Setting(:regime_eqcond_info, regime_eqcond_info)
setup_regime_switching_inds!(m; cond_type = :full)
set_regime_vals_fnct(m, get_setting(m, :n_regimes))
# Set up TVIS information set
m <= Setting(:tvis_information_set, vcat([i:i for i in 1:(gensys2_first_regime - 1)],
                                         [i:get_setting(m, :n_regimes) for i in
                                          gensys2_first_regime:get_setting(m, :n_regimes)]))

# Run forecast with 2 policies as reference output
output_vars = [:forecastobs, :forecastpseudo]
modal_params = map(x -> x.value, m.parameters)
# sys_2pol = compute_system(m; tvis = true)
out1 = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                              regime_switching = true, n_regimes = get_setting(m, :n_regimes))
# Test multiple alternative policies with multi-period altpolicies and temporary policies

## First start with fake temporary policy
temp_taylor_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
temp_default_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
for k in keys(temp_taylor_regime_eqcond_info)
    get_setting(m, :regime_eqcond_info)[k].weights =
        [get_setting(m, :regime_eqcond_info)[k].weights[1], get_setting(m, :regime_eqcond_info)[k].weights[2] / 2.,
         get_setting(m, :regime_eqcond_info)[k].weights[2] / 2.] # update the weights vector length
    temp_taylor_regime_eqcond_info[k] = DSGE.EqcondEntry(taylor_rule())
    temp_default_regime_eqcond_info[k] = DSGE.EqcondEntry(default_policy())
end
temp_taylor = MultiPeriodAltPolicy(:temporary_taylor, get_setting(m, :n_regimes), temp_taylor_regime_eqcond_info, gensys2 = true,
                                   temporary_altpolicy_names = [:taylor_rule],
                                   temporary_altpolicy_length = get_setting(m, :temporary_altpolicy_length),
                                   infoset = copy(get_setting(m, :tvis_information_set))) # also test w/ tvis and w/out
temp_default = MultiPeriodAltPolicy(:temporary_default, get_setting(m, :n_regimes), temp_default_regime_eqcond_info,
                                    gensys2 = false, temporary_altpolicy_names = [:default]) # also check without gensys2 on

delete!(DSGE.get_settings(m), :alternative_policies) # delete to update typing
m <= Setting(:alternative_policies, [DSGE.taylor_rule(), temp_taylor])
# sys_taylor = compute_system(m; tvis = true)
out_taylor_temp_taylor = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                                regime_switching = true, n_regimes = get_setting(m, :n_regimes))
m <= Setting(:alternative_policies, [DSGE.default_policy(), temp_default])
# sys_default = compute_system(m; tvis = true)
out_default_temp_default = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                                  regime_switching = true, n_regimes = get_setting(m, :n_regimes))

## Use implemented policy as nontrivial third policy
temp_flexait_zlb_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
for k in keys(temp_taylor_regime_eqcond_info)
    get_setting(m, :regime_eqcond_info)[k].weights =
        [get_setting(m, :regime_eqcond_info)[k].weights[1] / 2., 2 * get_setting(m, :regime_eqcond_info)[k].weights[2],
         get_setting(m, :regime_eqcond_info)[k].weights[1] / 2.] # need to multiply by 2 b/c of how weights are set previously
    if get_setting(m, :regime_eqcond_info)[k].alternative_policy.key == :zero_rate
        temp_flexait_zlb_regime_eqcond_info[k] = EqcondEntry(zero_rate())
    else
        temp_flexait_zlb_regime_eqcond_info[k] = EqcondEntry(flexible_ait())
    end
end
temp_flexait_zlb = MultiPeriodAltPolicy(:temporary_flexait_zlb, get_setting(m, :n_regimes), temp_flexait_zlb_regime_eqcond_info, gensys2 = true,
                                        temporary_altpolicy_names = [:zero_rate],
                                        temporary_altpolicy_length = get_setting(m, :temporary_altpolicy_length),
                                        infoset = copy(get_setting(m, :tvis_information_set)))
m <= Setting(:alternative_policies, [DSGE.default_policy(), temp_flexait_zlb])
# sys_flexait = compute_system(m; tvis = true)
out_flexait_zlb_temp_flexait_zlb = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                                          regime_switching = true, n_regimes = get_setting(m, :n_regimes))

## Now add nontrivial temporary policy that is different from implemented policy
temp_flexible_ait_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
for k in keys(temp_taylor_regime_eqcond_info)
    if get_setting(m, :regime_eqcond_info)[k].alternative_policy.key == :zero_rate
        temp_flexible_ait_regime_eqcond_info[k] = DSGE.EqcondEntry(flexible_ait()) # temporary flexible ait during ZLB periods
    else
        temp_flexible_ait_regime_eqcond_info[k] = DSGE.EqcondEntry(taylor_rule()) # followed by Taylor rule in the end
    end
end
temp_flexible_ait = MultiPeriodAltPolicy(:temporary_flexible_ait, get_setting(m, :n_regimes), temp_flexible_ait_regime_eqcond_info, gensys2 = true,
                                         temporary_altpolicy_names = [:flexible_ait],
                                         temporary_altpolicy_length = get_setting(m, :temporary_altpolicy_length),
                                         infoset = copy(get_setting(m, :tvis_information_set))) # also test w/ tvis and w/out

m <= Setting(:alternative_policies, [DSGE.default_policy(), temp_flexible_ait])
out_temp_flexible_ait = DSGE.forecast_one_draw(m, :mode, :full, output_vars, modal_params, df,
                                               regime_switching = true, n_regimes = get_setting(m, :n_regimes))

if regenerate_reference_forecasts
    h5open(joinpath(dirname(@__FILE__), "../reference/multiple_altpol_imperfect_awareness_output.h5"), "w") do file
        write(file, "forecastobs", out_temp_flexible_ait[:forecastobs])
        write(file, "forecastpseudo", out_temp_flexible_ait[:forecastpseudo])
    end
end

@testset "Multiple alternative policies, incl. MultiPeriodAltPolicy, with imperfect awareness" begin
    for k in keys(out_taylor_temp_taylor) # check implementations that should yield same results as a forecast
        @test out_taylor_temp_taylor[k] ≈ out1[k] # with 2 alternative policies and imperfect awareness
        @test out_default_temp_default[k] ≈ out1[k]
        @test out_flexait_zlb_temp_flexait_zlb[k] ≈ out1[k]
        @test maximum(abs.(out_temp_flexible_ait[k] -
                          h5read(joinpath(dirname(@__FILE__), "../reference/multiple_altpol_imperfect_awareness_output.h5"), string(k)))) < 1e-3
    end
end
nothing


text:
########################################################################
#
#    Reverse diff on graph
#
#########################################################################

###### creates reverse mode diff graph ######
function reversegraph(g::ExGraph, exitnode::ExNode, diffsym::Array{Symbol})
    g2     = ExGraph()  # graph that will contain the derivative evaluation

    # create starting nodes
    dnodes = Dict()     # map linking nodes of g to their derivative node in g2
    for n in filter(n-> !isa(n,NFor), g.nodes)
        if n == exitnode
            dnodes[n] = addnode!(g2, NConst(1.0))
        else
            dnodes[n] = addgraph!( zeronode(n), g2, Dict( :tv => n) )
            # dnodes[n] = createzeronode!(g2, n)
        end
    end

    # builds the graph for derivatives calculations
    reversepass!(g2, g, dnodes)

    # store in setmap the nodes containing the derivatives of diffsym
    for (k,v) in filter((k,v) -> isa(k, NExt) & in(k.main, diffsym), dnodes)
        g2.seti[v] = dprefix(k.main)
    end

    g2
end

#  climbs the reversed evaluation tree
function reversepass!(g2::ExGraph, g::ExGraph, dnodes::Dict)
    # TODO : have a generic treatment of mutating nodes (NSRef / NSDot / NFor ...)

    rev(n::ExNode) = nothing  # do nothing by default

	function rev(n::NCall)
		op = n.parents[1].main
		for (index, arg) in enumerate(n.parents)
			if !isa(arg, Union{NConst, NComp})
				# haskey(drules, (op, index-1)) || error("no derivation rule for $(op) at arg #$(index-1)")
				# ddict = drules[(op, index-1)]
                ddict = getrule(op, index-1)

                targs = Tuple{ Type[ typeof(x.val) for x in n.parents[2:end]]... }

                sk = tmatch( targs, collect(keys(ddict)) )
                (sk == nothing) && error("no derivation rule for $(op) at arg #$(index-1) for signature $targs")

                # dg, dd = drules[(op, index-1)][sk]
                dg, dd = ddict[sk]
            	smap = Dict( zip(dd, [n.parents[2:end]; dnodes[n]]) )

                exitnode = addgraph!(dg, g2, smap)

                vp = addnode!(g2, NConst(+))
                dnodes[arg] = addnode!(g2, NCall(:call, [vp, dnodes[arg], exitnode]) )
            end
        end
    end

    function rev(n::NRef)
        v2 = addnode!(g2, NRef(:getidx,  [ dnodes[n.parents[1]]; n.parents[2:end] ]) )
        vp = addnode!(g2, NConst(+))
        v3 = addnode!(g2, NCall(:call, [vp, v2, dnodes[n]]) )

		v4 = addnode!(g2, NSRef(:setidx, [ dnodes[n.parents[1]]; v3; n.parents[2:end] ]) )
		# TODO : update precedence of v4 here ? can 'dnodes[n.parents[1]' be already a parent elsewhere ?
		dnodes[n.parents[1]] = v4
	end

	function rev(n::NSRef)
		if length(n.parents) >= 3   # regular setindex
			v2 = addnode!(g2, NRef(:getidx, [ dnodes[n] ; n.parents[3:end] ]) )

			# treat case where a single value is allocated to several array elements
			if length(n.parents[2].val) == 1
				sz = mapreduce(x -> length(x.val), *, n.parents[3:end])
				if sz > 1
			       	vp = addnode!(g2, NConst(sum))
					v2 = addnode!(g2, NCall(:call, [ vp, v2]))
				end
			end

	       	vp = addnode!(g2, NConst(+))
			v3 = addnode!(g2, NCall(:call, [ vp, dnodes[n.parents[2]], v2 ]) )
			dnodes[n.parents[2]] = v3

			# shut down the influence of these indices
			zn = addnode!(g2, NConst(0.))
			v4 = addnode!(g2, NSRef(:setidx, [ dnodes[n]; zn; n.parents[3:end] ]) )
			v4.precedence = filter(n2 -> dnodes[n] in n2.parents && n2 != v4, g2.nodes)
			dnodes[n.parents[1]] = v4
		else   # simple assignment
	       	vp = addnode!(g2, NConst(+))
			v3 = addnode!(g2, NCall(:call, [ vp, dnodes[n.parents[2]], dnodes[n] ]) )
			dnodes[n.parents[2]] = v3

			# shut down the influence of the variable
			zn = addnode!(g2, NConst(0.))
			v4 = addnode!(g2, NSRef(:setidx, [ dnodes[n], zn ]) )
			v4.precedence = filter(n2 -> dnodes[n] in n2.parents && n2 != v4, g2.nodes)
			dnodes[n.parents[1]] = v4
		end
	end

	function rev(n::NDot)
		fsym = isa(n.main, Expr) ? n.main.args[1] : n.main.value  # can be Expr or QuoteNode
		idx = findfirst( fieldnames(typeof(n.parents[1].val)) .== fsym )
		(idx == 0) && error("[reversegraph] field $(n.main) not found in $(typeof(n.val))")

		v1 = addnode!(g2, NConst(idx) )
        v2 = addnode!(g2, NRef(:getidx,  [ dnodes[n.parents[1]], v1 ]) )
        vp = addnode!(g2, NConst(+))
        v3 = addnode!(g2, NCall(:call, [vp, v2, dnodes[n]]) )

        v4 = addnode!(g2, NSRef(:setidx, [ dnodes[n.parents[1]], v3, v1 ]) )
        # TODO : update precedence of v4 here ? can 'dnodes[n.parents[1]' be already a parent elsewhere ?
        dnodes[n.parents[1]] = v4
    end

    function rev(n::NSDot)
        fsym = isa(n.main, Expr) ? n.main.args[1] : n.main.value  # can be Expr or QuoteNode
        idx = findfirst( fieldnames(typeof(n.parents[1].val)) .== fsym )
        (idx == 0) && error("[reversegraph] field $(n.main) not found in $(typeof(n.val))")

        v1 = addnode!(g2, NConst(idx) )
        v2 = addnode!(g2, NRef(:getidx, [ dnodes[n] , v1 ]) )

        vp = addnode!(g2, NConst(+))
        v3 = addnode!(g2, NCall(:call, [ vp, dnodes[n.parents[2]], v2 ]) )
        dnodes[n.parents[2]] = v3

        # shut down the influence of these indices
        zn = addnode!(g2, NConst(0.))
        v4 = addnode!(g2, NSRef(:setidx, [ dnodes[n], zn, v1 ]) )
        v4.precedence = filter(n2 -> dnodes[n] in n2.parents && n2 != v4, g2.nodes)
        dnodes[n.parents[1]] = v4
    end

    function rev(n::NIn)
        isa(n.parents[1], NFor) && return nothing  # do nothing in the case of for loops

        vp = addnode!(g2, NConst(+))
        v2 = addnode!(g2, NCall(:call, [vp, dnodes[n], dnodes[n.parents[1]]]) )
        dnodes[n.parents[1]] = v2
    end

    function rev(n::NFor)
        fg  = copy(n.main[2])      # subgraph of for loop, copied to make new loop
        fg2 = ExGraph()            # will contain dnodes
        is  = n.main[1]            # symbol of loop index

        fdnodes = Dict()
        ndmap = Dict()

        nexti = NSMap()
        nexto = NSMap()
        # nexti = fg.exti
        # nexto = fg.exto
        # outgoing nodes generate ingoing dnodes
        for (n2, sym) in fg.seti
            hassym(fg.seto, sym)  || continue
            on = getnode(fg.seto, sym)
            dsym = newvar(:_dtmp)  # dprefix(sym)
            #  derivative of var
            nn = addnode!(fg2, NExt(dsym))
            nn.val = "seti"
            nexti[nn] = dsym
            nexto[dnodes[on]] = dsym
            # println("seti : $sym / $dsym : $nn")
            fdnodes[n2] = nn
        end

        # ingoing nodes become potential outgoing dnodes
        for (n2, sym) in fg.exti
            hassym(fg.exto, sym)  || continue
            on = getnode(fg.exto, sym)

            ## derivative accumulator
            # if haskey(nexti.vk, dsym)  # already mapped ?
            if hassym(fg.seto, sym)  # hassym(fg.seti, sym)  # already mapped ?
                on2 = getnode(fg.seti, sym)
                nn   = fdnodes[on2]  # nexti.vk[dsym]
                dsym = nexti[nn]
                # println("exti (refused) : $sym / $dsym : $nn")
            else
                dsym = newvar(:_dtmp)  # dprefix(sym)
                nn = addnode!(fg2, NExt(dsym))
                nn.val = "exti"
                nexti[nn] = dsym
                nexto[dnodes[on]] = dsym

                # println("exti : $sym / $dsym : $nn")
            end
            ndmap[n2] = (dsym, on)
            fdnodes[n2] = nn
        end

        fg.exti = NSMap(merge(fg.exti.kv, nexti.kv))
        fg.exto = NSMap(merge(fg.exto.kv, nexto.kv))

        # create regular zeronodes for the remaining fg nodes
        for n2 in filter(n-> !isa(n,NFor) & !haskey(fdnodes, n), fg.nodes)
            nn = addgraph!( zeronode(n2), fg2, Dict( :tv => n2) )
            # nn = createzeronode!(fg2, n2)
            fdnodes[n2] = nn
        end

        # builds the graph for derivatives calculations
        reversepass!(fg2, fg, fdnodes)
        append!(fg.nodes, fg2.nodes)

        # println("after reverse\n$fg")
        # variables of interest are derivatives only
        fg.seti = NSMap()
        for (ni, (sym, on)) in ndmap
            fg.seti[ fdnodes[ni] ] = sym
        end
        # println("after seti update\n$fg")

        prune!(fg) # reduce to derivatives evaluation only
        # println("after pruning\n$fg")

        # create for loop
        nr = addgraph!(:( reverse( x ) ), g2,
                       Dict( :x => n.parents[1] ) ) # range in reverse order
        v2 = addnode!(g2, NFor(Any[ n.main[1], fg ]) )
        v2.parents = [nr; collect( nodes( fg.exto)) ]

        # seto = dnodes of fg's ingoing variables
        fg.seto = NSMap()
        for (ns2, (sym, on)) in ndmap
            rn = addnode!(g2, NIn(sym, [v2]))  # external node, receiving loop result
            fg.seto[rn] = sym

            append!(v2.precedence, filter(n -> dnodes[on] in n.parents && n != v2, g2.nodes))
            dnodes[on] = rn
        end
    end

    evalsort!(g)
    for n2 in reverse(g.nodes)
        rev(n2)
    end
end


text:
# addions
"""
    function genr_ions(data::Data, ion_type::String; num_pairs=false, conc=false)

This will return a variable in "Data_Unit" type containing information of ions. Two options are offered, determining by which parameters, `num_pairs` or `conc`, has been specified.
- `num_pairs` will create ion composite with the number of `num_ions`.
- `conc` will create ion composite with the number which make the solution with a concentration of `conc`.

# Supported List
- Type          Mass(g/mol)
- ``Na⁺``      22.99
- ``K⁺``         39.09
- ``Al³⁺``     26.98
- ``Cl⁻``        35.453
Notcice: `ion_type` is case sensetive and space is needed for recogonition

# Example
```julia-repl
cell = genr_cell([30,20,3])
str = Tip3p()
data = genr_atom(cell, str)
data = addions(data, "Na Cl", conc=0.5)
```
"""
function genr_ions(data::Data, ion_type::String; num_pairs=false, conc=false)
    # Supported List
    list_ion = split("K Na Fe Al La Cl")
    list_residue = split("POT SOD FE AL LAN CLA")
    list_charge = [1, 1, 2, 3, 3, -1]
    list_mass = [39.09, 22.99, 55.845, 26.98, 138.04, 35.453]

    # Reading Input
    if isa(num_pairs, Bool) & isa(conc, Bool)
        error("One of `num_pairs` or `conc` should be specified")
    end

    if !(isa(num_pairs, Bool) | isa(conc, Bool))
        error("Only one of `num_pairs` or `conc` can be specified")
    end

    ion = split(ion_type)
    len = length(ion)
    ion_id =  [findall(x->x==ion[n], list_ion)[1] for n = 1:len]
    ion_residue_name = [list_residue[ion_id[n]] for n = 1:len]
    ion_charge = [list_charge[ion_id[n]] for n = 1:len]
    ion_mass = [list_mass[ion_id[n]] for n = 1:len]

    # Genrate ion_mode
    num_unit_ions = lcm(abs.(ion_charge))
    num_ions = convert.(Int64, num_unit_ions ./ abs.(ion_charge))
    num_unit_ions = sum(num_ions)

    unit_type = vcat([[j for i=1:num_ions[j]] for j=1:len]...) # Genrate type vec like [1 2 2 2]
    unit_charge = vcat([[ion_charge[j] for i=1:num_ions[j]] for j=1:len]...)

    center = (max(data.vec_atom, "coord")-min(data.vec_atom, "coord")) ./ 2
    ion_mode = [Atom(i, i, unit_type[i], unit_charge[i], center) for i = 1:num_unit_ions]

    # Calculate # of solute molecules
    if conc!= false
        ratio_sol2wat = conc / (Const_density_wat * Const_kg2g / Const_m2dm^3 / sum(data.data_str.para_mass))
        num_wat = data.data_basic.num_atoms
        num_sol = convert(Int64, round(num_wat * ratio_sol2wat))
    else
        num_sol = num_pairs
    end
    num_atoms = num_sol * num_unit_ions

    # Generate new Atoms
    vec_atom = data.vec_atom
    vec_atom_new = Vector{Atom}(undef, num_atoms)
    box_vec = data.data_basic.box_size[:,2] - data.data_basic.box_size[:,1]
    #max_mol = max(vec_atom, "mol")
    #max_type = max(vec_atom, "typ")
    for sol = 1:num_sol
        for ion = 1:num_unit_ions
            id_now = (sol-1) * num_unit_ions + ion
            atom_now = Atom(ion_mode[ion])
            atom_now.atom = id_now
            atom_now.mol = id_now
            atom_now.coord .+= (rand(3).-0.5) .* box_vec
            vec_atom_new[id_now] = atom_now
        end
    end

    # Generate Data instance
    vec_atom_cation = vec_atom_new[findall(x->x>0, get_data(vec_atom_new, "charge"))]
    for (id, atom) in enumerate(vec_atom_cation)
        atom.atom = id
        atom.mol = id
    end
    data_str = Structure_Ion(ion[1], ion_charge[1], ion_mass[1], 1, [1], ion_residue_name[1])
    data_basic = Data_Basic(data.data_basic)
    data_basic.num_atoms = length(vec_atom_cation)
    data_basic.num_atom_types = 1
    data_basic.num_bonds = 0
    data_basic.num_bond_types = 0
    data_basic.num_angles = 0
    data_basic.num_angle_types = 0
    data_basic.num_dihedrals = 0
    data_basic.num_dihedral_types = 0
    data_basic.num_impropers = 0
    data_basic.num_improper_types = 0
    data_cation = Data_Unit(data_basic, data_str, vec_atom_cation, 0, 0, 0, 0)

    vec_atom_anion = vec_atom_new[findall(x->x<0, get_data(vec_atom_new, "charge"))]
    for (id, atom) in enumerate(vec_atom_anion)
        atom.atom = id
        atom.mol = id
        atom.typ = 1
    end
    data_str = Structure_Ion(ion[2], ion_charge[2], ion_mass[2], 1, [1], ion_residue_name[2])
    data_basic = Data_Basic(data.data_basic)
    data_basic.num_atoms = length(vec_atom_anion)
    data_basic.num_atom_types = 1
    data_basic.num_bonds = 0
    data_basic.num_bond_types = 0
    data_basic.num_angles = 0
    data_basic.num_angle_types = 0
    data_basic.num_dihedrals = 0
    data_basic.num_dihedral_types = 0
    data_basic.num_impropers = 0
    data_basic.num_improper_types = 0
    data_anion = Data_Unit(data_basic, data_str, vec_atom_anion, 0, 0, 0, 0)

    # Output
    return data_cation, data_anion
end


text:
@testset verbose = true "Integration Tests" begin

    @testset "Clear board" begin
        # Write your tests here.
        env = NatureEnv(
            num_starting_players=2,
            world_size=(8, 8),
            window=0,
            player_starting_food=100f0,
            food_generators=[
                FoodGen([5,5],[1,1], max_pos=(8,8)),
                FoodGen([2,2],[1,1], max_pos=(8,8)),
               ])

        RLBase.reset!(env)
        Nature.move_and_collect(env, 1, 1)
        Nature.move_and_collect(env, 2, 2)
        @test env.food_frames[1:2] |> sum |> unique == [0]
    end

    @testset "Pick/Place" begin
        env = NatureEnv(
            num_starting_players=1,
            world_size=(1, 1),
            window=0,
            food_generators=[
                FoodGen([1,1],[1,1], max_pos=(1,1)),
               ])

        RLBase.reset!(env)
        @test env.players[1].pos == (1,1)
        @test env.players[1].food_counts == (env.player_starting_food,)
        f = env.food_frames[1][env.players[1].pos...]
        env(PICK_FOOD_1,1)
        @test env.players[1].food_counts[1] == env.player_starting_food + f - 0.1f0
        env(PLACE_FOOD_1,1)
        @test 0.1f0 > abs(env.players[1].food_counts[1] - (env.player_starting_food + f - 0.2f0 - 0.5f0))
    end

    @testset "Exchange" begin
        env = NatureEnv(
            num_starting_players=2,
            world_size=(1, 1),
            window=0,
            food_generators=[
                FoodGen([1,1],[1,1], max_pos=(1,1)),
               ])

        RLBase.reset!(env)
        ff = env.food_frames[1][1,1]
        # Test exchange works between agents
        env(PLACE_FOOD_1, 1)
        @test env.food_frames[1][1,1] == ff + 0.5f0
        env(PICK_FOOD_1, 2)
        @test env.exchanges[1] == 0.5f0
        # Test exchange doesn't work between same agent
        env(PLACE_FOOD_1, 1)
        env(PICK_FOOD_1, 1)
        @test env.exchanges[1] == 0.5f0
        RLBase.reset!(env)
        @test sum(env.exchanges) == 0f0
    end

    @testset "Comm" begin
        env = NatureEnv(
            num_starting_players=2,
            world_size=(1, 1),
            window=0,
            vocab_size=2,
            num_frames=1,
            food_generators=[
                FoodGen([1,1],[1,1], max_pos=(1,1)),
               ])

        RLBase.reset!(env)
        COMM_1 = 4+(2*env.food_types) + 1
        COMM_1_CHANNEL = 7
        COMM_2 = 4+(2*env.food_types) + 2
        COMM_2_CHANNEL = 8
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 0f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 0f0
        env.step += 1
        env(COMM_1, 1)
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 1f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 0f0
        env(COMM_2, 2)
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 1f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 1f0
        env.step += 1
        env(COMM_1, 1)
        env(COMM_2, 2)
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 1f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 1f0
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 1f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 1f0
        env.step += 1
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 0f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 0f0
        @test state(env, 1)[:,:,COMM_1_CHANNEL][1] == state(env, 2)[:, :, COMM_1_CHANNEL][1] == 0f0
        @test state(env, 1)[:,:,COMM_2_CHANNEL][1] == state(env, 2)[:, :, COMM_2_CHANNEL][1] == 0f0
        RLBase.reset!(env)
        @test length(env.comms) == env.episode_len
        @test all([isempty(comm_step) for comm_step in env.comms])
    end
end


text:
using ShellModel
using Test

@testset "ShellModel.jl" begin
    main("../snts/x_mass.snt","Be8",1,[];mdimmode=true)
end


text:
"""
PTCUpdate(FPS, FS, x, ItRules, step, residm, dt)

Updates the PTC iteration. This is a much simpler algorithm that Newton.
We update the Jacobian every iteration and there is no line search
to manage.

Do not mess with this function!
In particular do not touch the line that calls PrepareJac!.
        FPF = PrepareJac!(FPS, FS, x, ItRules,dt)
PrePareJac! builds F'(u) + (1/dt)I, factors it, and sends ptcsol that
factorization.

FPF is not the same as FPS (the storage you allocate for the Jacobian)
for a reason. FPF and FPS do not have the same type, even though they
share storage. So, FPS=PrepareJac!(FPS, FS, ...) will break things.

"""
function PTCUpdate(FPS, FS, x, ItRules, step, residm, dt)
    T = eltype(FPS)
    F! = ItRules.f
    pdata = ItRules.pdata
    #
    FPF = PrepareJac!(FPS, FS, x, ItRules, dt)
    #
    #    step .= -(FPF \ FS)
    T == Float64 ? (step .= -(FPF \ FS)) : (step .= -(FPF \ T.(FS)))
    #
    # update solution/function value
    #
    x .= x + step
    EvalF!(F!, FS, x, pdata)
    resnorm = norm(FS)
    #
    # Update dt
    #
    dt *= (residm / resnorm)
    return (x, dt, FS, resnorm)
end

"""
PTCUpdate(df::Real, fval, x, ItRules, step, residm, dt)

PTC for scalar equations.
"""

function PTCUpdate(df::Real, fval, x, ItRules, step, resnorm, dt)
    f = ItRules.f
    dx = ItRules.dx
    fp = ItRules.fp
    pdata = ItRules.pdata
    df = fpeval_newton(x, f, fval, fp, dx, pdata)
    idt = 1.0 / dt
    step = -fval / (idt + df)
    x = x + step
    #    fval = f(x)
    fval = EvalF!(f, fval, x, pdata)
    # SER
    residm = resnorm
    resnorm = abs(fval)
    dt *= (residm / resnorm)
    return (x, dt, fval, resnorm)
end


text:
function my_sum(a,b)
    return a + b;
end


text:
# ---
# title: 217. Contains Duplicate
# id: problem217
# author: zhwang
# date: 2022-02-10
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/contains-duplicate/description/>
# hidden: true
# ---
#
# Given an array of integers, find if the array contains any duplicates.
#
# Your function should return true if any value appears at least twice in the
# array, and it should return false if every element is distinct.
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,1]
#     Output: true
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,4]
#     Output: false
#
# **Example 3:**
#
#
#
#     Input: [1,1,1,3,3,4,3,2,4,2]
#     Output: true
#
#
## @lc code=start
using LeetCode

function contains_duplicate(nums::Vector{Int})::Bool
    hash_table = Int[]
    for i in nums
        !(i in hash_table) ? push!(hash_table, i) : return true
    end
    return false
end
## @lc code=end


text:
"""
    TupleView{T, N, Skip, A}

TupleView, groups elements of an array as tuples.
N is the dimension of the tuple, M tells how many elements to skip to the next tuple.
By default TupleView{N} defaults to skip N items.
# a few examples:

```julia

x = [1, 2, 3, 4, 5, 6]
TupleView{2, 1}(x):
> [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]

TupleView{2}(x):
> [(1, 2), (3, 4), (5, 6)]

TupleView{2, 3}(x):
> [(1, 2), (4, 5)]

TupleView{3, 1}(x):
> [(1, 2, 3), (2, 3, 4), (3, 4, 5), (4, 5, 6)]
```

TupleView can be used together with reinterpret:

```julia
x = [1, 2, 3, 4]
y = reinterpret(Point{2, Int}, TupleView{2, 1}(x))
> [Point(1, 2), Point(2, 3), Point(3, 4)]
```

"""
struct TupleView{T,N,Skip,A} <: AbstractVector{T}
    data::A
    connect::Bool
end

coordinates(tw::TupleView) = coordinates(tw.data)

function Base.size(x::TupleView{T,N,M}) where {T,N,M}
    nitems = length(x.data) ÷ (N - (N - M))
    nitems = nitems - max(N - M, 0)
    return (nitems + x.connect,) # plus one item if we connect
end

function Base.getindex(x::TupleView{T,N,M}, index::Integer) where {T,N,M}
    return ntuple(i -> x.data[mod1(((index - 1) * M) + i, length(x.data))], N)
end

function TupleView{N}(x::AbstractVector; connect=false) where {N}
    return TupleView{N,N}(x, connect=connect)
end

function TupleView{N,M}(x::AbstractVector{T}; connect=false) where {T,N,M}
    return TupleView{NTuple{N,T},N,M,typeof(x)}(x, connect)
end

@inline function connected_line(points::AbstractVector{<:AbstractPoint{N}},
                                skip=N) where {N}
    return connect(points, Line, skip)
end

"""
    connect(points::AbstractVector{<: AbstractPoint}, P::Type{<: Polytope{N}}, skip::Int = N)

Creates a view that connects a number of points to a Polytope `P`.
Between each polytope, `skip` elements are skipped untill the next starts.
Example:
```julia
x = connect(Point[(1, 2), (3, 4), (5, 6), (7, 8)], Line, 2)
x == [Line(Point(1, 2), Point(3, 4)), Line(Point(5, 6), Point(7, 8))]
"""
@inline function connect(points::AbstractVector{Point},
                         P::Type{<:Polytope{N,T} where {N,T}},
                         skip::Int=length(P)) where {Point <: AbstractPoint}
    return reinterpret(Polytope(P, Point), TupleView{length(P),skip}(points))
end

@inline function connect(points::AbstractVector{T}, P::Type{<:Point{N}},
                         skip::Int=N) where {T <: Real,N}
    return reinterpret(Point{N,T}, TupleView{N,skip}(points))
end

@inline function connect(points::AbstractVector{T}, P::Type{<:AbstractFace{N}},
                         skip::Int=N) where {T <: Real,N}
    return reinterpret(Face(P, T), TupleView{N,skip}(points))
end

@inline function connect(points::AbstractMatrix{T},
                         P::Type{<:AbstractPoint{N}}) where {T <: Real,N}
    return if size(points, 1) === N
        return reinterpret(Point{N,T}, points)
    elseif size(points, 2) === N
        seglen = size(points, 1)
        columns = ntuple(N) do i
            return view(points, ((i - 1) * seglen + 1):(i * seglen))
        end
        return StructArray{Point{N,T}}(columns)
    else
        error("Dim 1 or 2 must be equal to the point dimension!")
    end
end

"""
    FaceView{Element, Point, Face, P, F}

FaceView enables to link one array of points via a face array, to generate one
abstract array of elements.
E.g., this becomes possible:
```
x = FaceView(rand(Point3f, 10), TriangleFace[(1, 2, 3), (2, 4, 5), ...])
x[1] isa Triangle == true
x isa AbstractVector{<: Triangle} == true
# This means we can use it as a mesh:
Mesh(x) # should just work!
# Can also be used for Points:

linestring = FaceView(points, LineFace[...])
Polygon(linestring)
```
"""
struct FaceView{Element,Point <: AbstractPoint,Face <: AbstractFace,P <: AbstractVector{Point},F <: AbstractVector{Face}} <: AbstractVector{Element}
    elements::P
    faces::F
end

const SimpleFaceView{Dim,T,NFace,IndexType,PointType <: AbstractPoint{Dim,T},FaceType <: AbstractFace{NFace,IndexType}} = FaceView{Ngon{Dim,T,NFace,PointType},PointType,FaceType,Vector{PointType},Vector{FaceType}}

function Base.getproperty(faceview::FaceView, name::Symbol)
    return getproperty(getfield(faceview, :elements), name)
end

function Base.propertynames(faceview::FaceView)
    return propertynames(getfield(faceview, :elements))
end

Tables.schema(faceview::FaceView) = Tables.schema(getfield(faceview, :elements))

Base.size(faceview::FaceView) = size(getfield(faceview, :faces))

@propagate_inbounds function Base.getindex(x::FaceView{Element}, i) where {Element}
    return Element(map(idx -> coordinates(x)[idx], faces(x)[i]))
end

@propagate_inbounds function Base.setindex!(x::FaceView{Element}, element::Element,
                                            i) where {Element}
    face = faces(x)[i]
    for (i, f) in enumerate(face) # TODO unroll!?
        coordinates(x)[face[i]] = element[i]
    end
    return element
end

function connect(points::AbstractVector{P},
                 faces::AbstractVector{F}) where {P <: AbstractPoint,F <: AbstractFace}
    return FaceView{Polytope(P, F),P,F,typeof(points),typeof(faces)}(points, faces)
end

coordinates(mesh::FaceView) = getfield(mesh, :elements)
faces(mesh::FaceView) = getfield(mesh, :faces)


text:
# Launch with `julia --project --track-allocation=user`
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
import Profile

include("common.jl")
case_name = ENV["ALLOCATION_CASE_NAME"]
@info "Recording allocations for $case_name"
sim = init_sim(case_name)
tendencies = copy(sim.state.prog)
update_n(sim, tendencies, 1) # compile first
Profile.clear_malloc_data()
update_n(sim, tendencies, 1)

# Quit julia (which generates .mem files), then call
#=
import Coverage
allocs = Coverage.analyze_malloc("src")
=#