text:
# Unit 01 - Linear Classifiers and Generalizations

function perceptron(set, T=1000)
    n = length(set); d = length(set[1][1])
    θ = zeros(d); θ₀ = 0;
    for t in 1:T
        ϵ = 0
        for i in 1:n
            x = set[i][1]; y = set[i][2]
            if (y*(θ' * x + θ₀) <= 0)
                θ  = θ + y * x
                θ₀ = θ₀ + y
                ϵ += 1
            end
        end
        if (ϵ == 0)
            return (θ,θ₀,ϵ,t)
        end
    end
    return (θ,θ₀,ϵ,T)
end

T = 5000
set = [([7,8],-1),([4,2],-1),([2,7],1),([3,-2],1),([-3,-2],1)]
(θ,θ₀,ϵ,t) = perceptron(set,T)

using Plots, LinearAlgebra

θₙ = θ/norm(θ)
θ₀ₙ = θ₀/norm(θ) *  θₙ
x⁺ = vcat([([tuple[1][1] tuple[1][2]])  for tuple in set if tuple[2] == 1]...)
x⁻ = vcat([([tuple[1][1] tuple[1][2]])  for tuple in set if tuple[2] == -1]...)

scatter(x⁺[:,1],x⁺[:,2], label="+1")
scatter!(x⁻[:,1],x⁻[:,2], label="-1")
#plot!([0,θ[1]],[0,θ[2]])
plot!([-θ₀ₙ[1],θₙ[1]-θ₀ₙ[1]],[-θ₀ₙ[2],θₙ[2]-θ₀ₙ[2]], arrow=1., label = "theta" )

set = [([3,3],1),([4,2],1),([2,7],-1),([3,-2],-1)]

function perceptronOrigin(set, T=30)
    println("*** New Lookup ***")
    n  = length(set); d = length(set[1][1])
    θ  = zeros(d)
    ϵt = 0
    for t in 1:T
        ϵ = 0
        for i in 1:n
            x = set[i][1]; y = set[i][2]
            error = false
            θpre = θ
            if (y*(θ' * x) <= 0)
                θ  = θ + y * x
                ϵ += 1; error = true; ϵt += 1
            end
            println("θpre: $θpre x: $x y: $y error: $error θpost=$θ")
        end
        if (ϵ == 0)
            return (θ,ϵ,t, ϵt)
        end
    end
    return (θ,ϵ,T,ϵt)
end

set = [([-1,-1],1),([1,0],-1),([-1,1.5],1)]
perceptronOrigin(set)

set = [([1,0],-1),([-1,1.5],1),([-1,-1],1)]
perceptronOrigin(set)

set = [([-1,-1],1),([1,0],-1),([-1,10],1)]
perceptronOrigin(set)

set = [([1,0],-1),([-1,10],1),([-1,-1],1)]
perceptronOrigin(set)

function perceptronFull(set; T=30, θ = zeros(length(set[1])), θ₀ = 0)
    println("*** New Classifier Lookup ***")
    println("Data: $set")
    println("Parameters: T = $T, θ=$θ, θ₀=$θ₀")
    n = length(set); d = length(set[1][1])
    for t in 1:T
        ϵ = 0
        for i in 1:n
            x = set[i][1]; y = set[i][2]
            error = false
            θpre = θ
            θ₀pre = θ₀
            if (y*(θ' * x + θ₀) <= 0)
                θ  = θ + y * x
                θ₀ = θ₀ + y
                ϵ += 1
                error = true
            end
            println("θpre: $θpre θ₀pre: $θ₀pre x: $x y: $y error: $error θpost=$θ θ₀post=$θ₀")
        end
        if (ϵ == 0)
            return (θ,θ₀,ϵ,T)
        end
    end
    return (θ,θ₀,ϵ,T)
end

set = [([-4,2],1),([-2,1],1),([-1,-1],-1),([2,2],-1),([1,-2],-1)]
(θ,θ₀,ϵ,T) = perceptronFull(set)

## Homework 1

### Tab 2

scatter([-1,1],[1,-1],label="positive")
scatter!([1,2],[1,2],label="negative")

## Homework 2

### Tab 1

λ = 0.5
x = [1,0]
y = 1

θ =
loss(θ,y=1,x=[1,0],λ=0.5) = max(0,1-y*(θ' * x)) + λ/2 * norm(θ)^2

loss([0,0])
loss([1,0])
loss([0.5,0])
loss([-1,0])
loss([1,-1])

max(0,1-1*([1,0]' * [1,0]))
max(0,1-1*([0.5,0]' * [1,0]))

norm()

### Tab 3

cos(0)
cos(π)
cos(2*π)
cos(3*π)
cos(4*π)

function setFactory(d,labels=ones(d))
    out = Tuple{Array{Int64,1},eltype(labels)}[] # or just out=[]
    #out = []
    for i in 1:d
      x = zeros(d)
      x[i] = cos(i*π)
      push!(out,(x,labels[i]))
    end
    return out
end

set = setFactory(20)
using Random
shuffle!(set)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)

randLabel = rand(20)
set = setFactory(20,randLabel)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)
shuffle!(set)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)

randLabel = rand(3)
set = setFactory(3,randLabel)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)

set = setFactory(3)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)

set=setFactory(2)
(θ,ϵ,t,ϵt) = perceptronOrigin(set,400)

max_stars_repo_path: Unit 01 - Linear Classifiers and Generalizations/Unit 01 - Linear Classifiers and Generalizations.jl


text:
using Clang.Generators
using MacroTools
using GEOS_jll

function rewriter(dag::ExprDAG)
    for n in dag.nodes
        map!(rewriter, n.exprs, n.exprs)
    end
end

"Functions that return a Cstring are wrapped in unsafe_string to return a String"
function rewriter(x::Expr)
    if @capture(x,
        function f_(fargs__)
            ccall(fname_, rettype_, argtypes_, argvalues__)
        end
    )
        # it is a function wrapper around a ccall

        # bind the ccall such that we can easily wrap it
        cc = :(ccall($fname, $rettype, $argtypes, $(argvalues...)))

        cc2 = if rettype == :Cstring
            :(unsafe_string($cc))
        else
            cc
        end

        # stitch the modified function expression back together
        x2 = :(function $f($(fargs...))
            $cc2
        end) |> prettify
        x2
    else
        # do not modify expressions that are no ccall function wrappers
        x
    end
end

cd(@__DIR__)

include_dir = joinpath(GEOS_jll.artifact_dir, "include") |> normpath

# wrapper generator options
options = load_options(joinpath(@__DIR__, "generate.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

header_dir = include_dir
headers = [joinpath(header_dir, "geos_c.h")]

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx, BUILDSTAGE_NO_PRINTING)
rewriter(ctx.dag)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)

ext: jl


text:
function probabilityM1(s₁::BigInt, t₁::BigFloat, n::BigInt)
    pM1 = factorial(s₁) ./ (exp(t₁) * (n + 1) ^(s₁ + 1))
    pM2 = (factorial(n) * factorial(s₁)) ./ factorial(n + s₁ + 1)
    p = pM1 ./(pM1 + pM2)
    p = Float64(p)
    return p
end

s1calc(x) = sum(x)
t1calc(x::Array{BigInt, 1}) = Float64(sum(log.(map(a -> factorial(a), x))))
t1calc(x::Array{Int64, 1}) = sum(log.(map(a -> factorial(a), x)))

function poissonsimulator(params, cst, targetdata)

    λ = params[1]
    simdata = rand(Poisson(λ), 100)
    s₁ = s1calc(simdata)
    if sum(simdata.>19) > 0
        t₁ = t1calc(map(x -> BigInt(x), simdata))
    else
        t₁ = t1calc(simdata)
    end
    euclidean(targetdata, [s₁, t₁]), 1
    #euclidean(targetdata[2], [s₁, t₁][2]), 1
end

function geometricsimulator(params, cst, targetdata)

    μ = params[1]
    simdata = rand(Geometric(μ), 100)
    s₁ = s1calc(simdata)
    if sum(simdata.>19) > 0
        t₁ = t1calc(map(x -> BigInt(x), simdata))
    else
        t₁ = t1calc(simdata)
    end
    euclidean(targetdata, [s₁, t₁]), 1
    #euclidean(targetdata[2], [s₁, t₁][2]), 1
end

function generatedata()
    λ = 0.5
    n = 100
    data = rand(Poisson(λ), n)
    s₁ = s1calc(data)
    t₁ = t1calc(data)
    pM1 = probabilityM1(BigInt(s₁), BigFloat(t₁), BigInt(n))

  return [s₁, t₁], pM1
end

################################################################
println("Test Bayes factors are calculated correctly (within 5% of true value)")
srand(1234)

cst = [[i] for i in 1:2]

ABCsetupsmc = ABCSMCModel([poissonsimulator, geometricsimulator], [1, 1], 1.0,
[Prior([Exponential(1)]), Prior([Uniform(0.0, 1.0)])];
nparticles = 500, maxiterations = 10^7)

ABCsetuprej = ABCRejectionModel([poissonsimulator, geometricsimulator], [1, 1], 1.0,
[Prior([Exponential(1)]), Prior([Uniform(0.0, 1.0)])];
nparticles = 500, maxiterations = 10^7)

td, pM1 = generatedata()

println("\t Checking ABC SMC model selection")
@time abcressmc = runabc(ABCsetupsmc, td);
@test isapprox(pM1, abcressmc.modelprob[1], rtol = 0.05)

println("\t Checking ABC rejection model selection")
@time abcresrej = runabc(ABCsetuprej, td);
@test isapprox(pM1, abcresrej.modelfreq[1], rtol = 0.05)

println("\t Check no errors arising from plotting")

plotmodelposterior(abcressmc, save = true)
@test isfile("ABCSMCmodelposteriors.pdf")
rm("ABCSMCmodelposteriors.pdf")

plotmodelposterior(abcresrej, save = true)
@test isfile("ABCRejectionmodelposteriors.pdf")
rm("ABCRejectionmodelposteriors.pdf")

plotparameterposterior(abcressmc, 1, save = true)
@test isfile("ABCSMCparameterposteriors-model1.pdf")
rm(("ABCSMCparameterposteriors-model1.pdf"))

plotparameterposterior(abcresrej, 1, save = true)
@test isfile("ABCRejectionparameterposteriors-model1.pdf")
rm("ABCRejectionparameterposteriors-model1.pdf")

writeoutput(abcressmc)
@test isfile("SMCModel-outputmodel1.txt")
@test isfile("SMCModel-outputmodel2.txt")
@test isfile("SMCModel-outputmodelprobabilities.txt")
rm("SMCModel-outputmodel1.txt")
rm("SMCModel-outputmodel2.txt")
rm("SMCModel-outputmodelprobabilities.txt")

writeoutput(abcresrej)
@test isfile("RejectionModel-outputmodel1.txt")
@test isfile("RejectionModel-outputmodel2.txt")
@test isfile("RejectionModel-outputmodelprobabilities.txt")
rm("RejectionModel-outputmodel1.txt")
rm("RejectionModel-outputmodel2.txt")
rm("RejectionModel-outputmodelprobabilities.txt")


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#* $Id$
#*--------------------------------------------------------------------
type Indutor
	Indutor()=begin
		new(
			indutance(),
			wire(Dict{Symbol,Any}(
				:Brief=>"Inlet",
				:PosX=>0.4638,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			wire(Dict{Symbol,Any}(
				:Brief=>"Outlet",
				:PosX=>0.4638,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			[
				:(inlet.V - outlet.V = L * diff(inlet.i)),
				:(outlet.i = inlet.i),
			],
			[
				"","",
			],
			[:L,],
			[:inlet,:outlet,]
		)
	end
	L::indutance
	inlet::wire
	outlet::wire
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Indutor
function setEquationFlow(in::Indutor)
	addEquation(1)
	addEquation(2)
end
function atributes(in::Indutor,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Indutor"
	fields[:Brief]="Electrical Indutor."
	drive!(fields,_)
	return fields
end
Indutor(_::Dict{Symbol,Any})=begin
	newModel=Indutor()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"""
`jutrie` is a pure Julia implementation of the
`trie` (prefix tree) data structure.
"""

import Base

mutable struct Node
    value
    # If the children is `Vector`
    # we can map it by Char's ASCII
    # Else if we want to support Unicode
    # we can Hierarchical it by `0` and `1`, deepen the Trie.
    # To simply,
    # we use the map structure
    children::Dict{Char, Node}

    Node(value=0) = new(value, Dict{Char, Node}())
end

struct __NULL end
NULL = __NULL()
isNULL(x) = (typeof(x) == __NULL)
# When we use the node, we need to check if it is NULL

struct Trie
    root::Node
    keyDict::Dict{String, Node}

    Trie() = new(Node(), Dict{String, Node}())
end

"""
    update!(trie::Trie, dataDict::Dict{String, Number})

Update the key and value by dataDict.
"""
function update!(trie::Trie, dataDict::Dict{String, T}) where {T <: Number}
    for (key, value) in dataDict
        set_value!(trie, key, value)
    end
    return trie
end

"""
    increase!(trie::Trie, key::String, value::Number)

It will increase VALUE on key,
if the key didn't exist,
we will initalize it by value.
"""
function increase!(trie::Trie, key::String, value::Number)
    set_value!(trie, key, exist_key(trie, key) ? get_value(trie, key)+value : value )
end

function set_value!(trie::Trie, key::String, value::Number)
    if ( !exist_key(trie, key) )
        node = trie.root
        for part in key
            if !(part in keys(node.children))
                node.children[part] = Node()
            end
            node = node.children[part]
        end
    else
        node = find_prefix_node(trie, key)
    end
    node.value = value
end

function get_value(trie::Trie, key::String)
    if ( !exist_key(trie, key) )
        error("KeyError: key $(key) not found")
    else
        return find_prefix_node(trie, key).value
    end
end

"""
    keys(trie::Trie) -> Array{String, 1}
Return the elements which not zero
"""
function Base.keys(trie::Trie)
    trieKeys = []

    function traversal(node, prefix::String)
        for (key, child) in node.children
            if child.value != 0
                append!(trieKeys, [string(prefix, key)])
            end
            traversal(child, string(prefix, key))
        end
    end

    traversal(trie.root, "")
    return trieKeys
end

function items(trie::Trie)

end

"""
    test_key(trie::Trie, key::String) -> Bool

The external interface to check if the key is valid.
"""
function test_key(trie::Trie, key::String)
    node = find_prefix_node(trie, key)
    return !isNULL(node) && (node.value > 0)
end

"""
    exist_key(trie::Trie, key::String) -> Bool

The internal interface to check if the key node is NULL.
"""
function exist_key(trie::Trie, key::String)
    return !isNULL( find_prefix_node(trie, key) )
end

function prefix_items(trie::Trie, prefix::String)
    node = find_prefix_node(trie, prefix)
    if isNULL(node) return [] end

    # Return all node's children element
    items = []
    function itemAppend(node, prefix)
        for (key, child) in node.children
            append!(items, (prefix+key, node.value))
            itemAppend(child)
        end
    end
    itemAppend(node)

    # Find the Real Keys
    return filter(p -> p[2] > 0, items)
end

function find_prefix_node(trie::Trie, prefix::String)
    # Find the prefix node
    node = trie.root
    for part in prefix
        if (part in keys(node.children))
            node = node.children[part]
        else return NULL end
    end
    return node
end

function longest_prefix_item(trie::Trie, prefix::String)
    longest_prefix_key, longest_prefix_value = NULL, NULL
    for (key, value) in prefix_items(trie, prefix)
        if (isNULL(longest_prefix_key)) ||
            (length(key) > length(longest_prefix_key))
            longest_prefix_value = value
            longest_prefix_key = key
        end
    end
    return ( isNULL(longest_prefix_key) ? () :
            (longest_prefix_key, longest_prefix_value) )
end

function test_before()
    trie = Trie()
    update!(trie, Dict("abc" => 3, "all" => 1, "bad" => 4))
end

function test_keys()
    trie = test_before()
    show(keys(trie))
end

function test_get_value()
    trie = test_before()
    get_value(trie, "all") |> println
    try
        get_value(trie, "agc")
    catch
        println("agc not found")
    end
    increase!(trie, "all", 10) |> println
    get_value(trie, "all") |> println
end

function main()
    test_keys()
    test_get_value()
end

main()

max_stars_repo_name: gajanlee/toys

max_issues_repo_name: gajanlee/toys

max_forks_repo_name: gajanlee/toys


text:
function fibonacci(limit)
    println(0)
    x = BigInt(0)
    y = BigInt(1)
    for i = 1:limit - 1
        result = BigInt(x + y)
        y = BigInt(x)
        x = BigInt(result)
        println(result)
    end
end
fibonacci(100)

max_stars_repo_name: ellygaytor/Projects

max_issues_repo_name: ellygaytor/Projects

max_forks_repo_name: ellygaytor/Projects


text:
@testset "Low Pass Tests" begin
    @testset "Agreement with explicit smoothing for random data" begin
        ts = cumsum(rand(1000))
        ys = randn(1000)
        dys = 0.1 .+ rand(1000)

        i = 100
        tau = pi

        wts = exp.(-abs.(ts.-ts[i])./tau)
        yout = sum(ys.*wts)/sum(wts)
        dyout = sqrt(sum(wts.*wts.*dys.*dys)/(sum(wts)*sum(wts)))

        ys_filt, dys_filt = LowPass.exp_lowpass(ts, ys, dys, tau)

        @test abs(ys_filt[i] - yout) < 1e-10
        @test abs(dys_filt[i] - dyout) < 1e-10
    end
end


text:
struct FenwickTree{T}
    bi_tree::Vector{T} #bi_tree is shorthand for Binary Indexed Tree, an alternative name for Fenwick Tree
    n::Int
end

"""
    FenwickTree{T}(n::Integer)

Constructs a [`FenwickTree`](https://en.wikipedia.org/wiki/Fenwick_tree) of length `n`.

"""
FenwickTree{T}(n::Integer) where T = FenwickTree{T}(zeros(T, n), n)

"""
    FenwickTree(counts::AbstractArray)

Constructs a [`FenwickTree`](https://en.wikipedia.org/wiki/Fenwick_tree) from an array of `counts`

"""
function FenwickTree(a::AbstractVector{U}) where U
    n = length(a)
    tree = FenwickTree{U}(n)
    @inbounds for i = 1:n
        inc!(tree, i, a[i])
    end
    tree
end

length(ft::FenwickTree) = ft.n
Base.eltype(::Type{FenwickTree{T}}) where T = T

"""
    inc!(ft::FenwickTree{T}, ind::Integer, val)

Increases the value of the [`FenwickTree`] by `val` from the index `ind` upto the length of the Fenwick Tree.

"""
function inc!(ft::FenwickTree{T}, ind::Integer, val = 1) where T
    val0 = convert(T, val)
    i = ind
    n = ft.n
    @boundscheck 1 <= i <= n || throw(ArgumentError("$i should be in between 1 and $n"))
    @inbounds while i <= n
        ft.bi_tree[i] += val0
        i += i&(-i)
    end
end

"""
    dec!(ft::FenwickTree, ind::Integer, val)

Decreases the value of the [`FenwickTree`] by `val` from the index `ind` upto the length of the Fenwick Tree.

"""
dec!(ft::FenwickTree, ind::Integer, val = 1 ) = inc!(ft, ind, -val)

"""
    incdec!(ft::FenwickTree{T}, left::Integer, right::Integer, val)

Increases the value of the [`FenwickTree`] by `val` from the indices from `left` and decreases it from the `right`.

"""
function incdec!(ft::FenwickTree{T}, left::Integer, right::Integer, val = one(T)) where T
    val0 = convert(T, val)
    inc!(ft, left, val0)
    dec!(ft, right, val0)
end

"""
    prefixsum(ft::FenwickTree{T}, ind::Integer)

Return the cumulative sum from index 1 upto `ind` of the [`FenwickTree`](@ref)

# Examples
```
julia> f = FenwickTree{Int}(6)
julia> inc!(f, 2, 5)
julia> prefixsum(f, 1)
 0
julia> prefixsum(f, 3)
 5
```
"""
function prefixsum(ft::FenwickTree{T}, ind::Integer) where T
    sum = zero(T)
    i = ind
    n = ft.n
    @boundscheck 1 <= i <= n || throw(ArgumentError("$i should be in between 1 and $n"))
    @inbounds while i > 0
        sum += ft.bi_tree[i]
        i -= i&(-i)
    end
    sum
end

getindex(ft::FenwickTree{T}, ind::Integer) where T = prefixsum(ft, ind)


text:
import GR

GR.selntran(0)
GR.setcharheight(0.024)

function main()
    GR.settextfontprec(2, 0)
    y = 0
    for i in 0:1
        GR.text(0.05, 0.85-y, " !\"#\$%&'()*+,-./")
        GR.text(0.05, 0.80-y, "0123456789:;<=>?")
        GR.text(0.05, 0.75-y, "@ABCDEFGHIJKLMNO")
        GR.text(0.05, 0.70-y, "PQRSTUVWXYZ[\\]^_")
        GR.text(0.05, 0.65-y, "`abcdefghijklmno")
        GR.text(0.05, 0.60-y, "pqrstuvwxyz{|}~ ")

        GR.text(0.5, 0.85-y, " ¡¢£¤¥¦§¨©ª«¬®¯")
        GR.text(0.5, 0.80-y, "°±²³´µ¶·¸¹º»¼½¾¿")
        GR.text(0.5, 0.75-y, "ÀÁÂÃÄÅÆÇÈÉÊËÌÍÎÏ")
        GR.text(0.5, 0.70-y, "ÐÑÒÓÔÕÖ×ØÙÚÛÜÝÞß")
        GR.text(0.5, 0.65-y, "àáâãäåæçèéêëìíîï")
        GR.text(0.5, 0.60-y, "ðñòóôõö÷øùúûüýþÿ")

        GR.settextfontprec(233, 3)
        y = 0.4
    end
    GR.updatews()
end

main()


text:
# ============================================================================ #
# Line
# ============================================================================ #
abstract type AbstractLine{T} end

struct Line{PointType <: AbstractPoint} # <: AbstractLine
    p1::PointType
    p2::PointType
end

Line(a::A, b::A) where {A <: AbstractArray} = Line( Point(a[1], a[2]), Point(b[1], b[2]) )

# ---------------------------------------------------------------------------- #
# Methods
# ---------------------------------------------------------------------------- #

Base.:(==)(line1::Line, line2::Line) = line1.p1 == line2.p1 && line1.p2 == line2.p2

function shift(line::T, v) where {T <: AbstractLine}
    return T( line.p1 + v, line.p2 + v )
end
function scale(line::T, factor) where {T <: AbstractLine}
    return T( line.p1 * factor, line.p2 * factor )
end

function Base.getindex(line::AbstractLine, index::Int)
    if index == 1
        return line.p1
    elseif index == 2
        return line.p2
    else
        return NaN
    end
end

function vector(line::Line)
    line.p2 - line.p1
end

function line_intersection(line1::Line, line2::Line, abs_tol::Real=0, rel_tol::Real=abs_tol>0 ? 0 : sqrt(eps))
    rhs = line1.p1 - line1.p2
    n1 = orthogonal_vector( vector(line1) )
    n2 = orthogonal_vector( vector(line2) )
    lamdas = n1 ./ n2
    if isapprox(lamdas[1], lamdas[2], atol = abs_tol, rtol = rel_tol)   # lamdas[1] == lamdas[2]
        return Point(v)
    end
    A = [n1, n2]
    t12 = A \ rhs
    v = t12[1] > t12[2] ? n1 .* t12[1] : n2 .* t12[2]
    return Point(v)
end

function segment_intersection(line1_a::P, line1_b::P, line2_a::P, line2_b::P) where {P <: AbstractPoint}

    d1 = point_orientation(line1_a, line1_b, line2_a)
    d2 = point_orientation(line1_a, line1_b, line2_b)
    d3 = point_orientation(line2_a, line2_b, line1_a)
    d4 = point_orientation(line2_a, line2_b, line1_b)

    b12 = (d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)
    b34 = (d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0)

    if b12 && b34
        return true

    elseif d1 == 0 && on_segment_paralell(line2_a, line2_b, line1_a)
        return true

    elseif d2 == 0 && on_segment_paralell(line2_a, line2_b, line1_b)
        return true

    elseif d3 == 0 && on_segment_paralell(line1_a, line1_b, line2_a)
        return true

    elseif d4 == 0 && on_segment_paralell(line1_a, line1_b, line2_b)
        return true

    else
        return false
    end

end

function segment_intersection(line1::Line, line2::Line)
    return segment_intersection(line1.p1, line1.p2, line2.p1, line2.p2)
end

function on_segment_paralell(line::Line, point::AbstractPoint)
# a point which lies on the line but not known whether on the segment
    return  (min(line.p1.x, line.p2.x) <= point.x <= max(line.p1.x, line.p2.x)) &&
            (min(line.p1.y, line.p2.y) <= point.y <= max(line.p1.y, line.p2.y))
end

function on_segment_paralell(line_a, line_b, point::AbstractPoint)
# a point which lies on the line but not known whether on the segment
    return  (min(line_a.x, line_b.x) <= point.x <= max(line_a.x, line_b.x)) &&
            (min(line_a.y, line_b.y) <= point.y <= max(line_a.y, line_b.y))
end

function point_orientation(p_base::P, p2::P, p3::P) where {P <: AbstractPoint}
    return cross2(p3 - p_base, p2 - p_base)
end

function discretize_for_plot(line::Line)
    x = [line.p1.x, line.p2.x]
    y = [line.p1.y, line.p2.y]
    return (x, y)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Test
using Distributions
using Caesar_VariableFactorExamples
using JSON

@testset "ExamplePrior" begin
    include("ExamplePrior.jl")
end


text:
#=
Type and methods facilitating simple but effective syntax interface for defining trading rules
=#

#TODO: figure out how to make this a function that interfaces with the portfolio & account objects
struct Rule{S,F,T}
    trigger::S
    action::F
    args::Tuple{Vararg{T}}
    function Rule(trigger::S, action::F, args::Tuple{Vararg{T}}) where {S<:Signal, F<:Function, T}
        return new{S,F,T}(trigger, action, args)
    end
end

macro rule(logic::Expr, args...)
    trigger = :($(logic.args[2]))
    #action = :($(logic.args[3])$((args...)))
    action = :($(logic.args[3]))
    args = :($(args...))
    return esc(:(Rule($trigger, $action, $args)))
end

→(a,b) = a ? b() : nothing


text:
"""
    Variant(b::Bgen, offset::Integer)
    Variant(io, offset, compression, layout, expected_n)
Parse information of a single variant beginning from `offset`.
"""
function Variant(io::IOStream, offset::Integer,
        compression::Integer, layout::Integer,
        expected_n::Integer)
    seek(io, offset)
    if eof(io)
        @error "reached end of file"
    end
    if layout == 1
        n_samples = read(io, UInt32)
    else
        n_samples = expected_n
    end
    if n_samples != expected_n
        @error "number of samples does not match"
    end
    varid_len = read(io, UInt16)
    varid = String(read(io, varid_len))
    rsid_len = read(io, UInt16)
    rsid = String(read(io, rsid_len))
    chrom_len = read(io, UInt16)
    chrom = String(read(io, chrom_len))
    pos = read(io, UInt32)

    if layout == 1
        n_alleles = 2
    else
        n_alleles = read(io, UInt16)
    end

    alleles = Array{String}(undef, n_alleles)
    for i in 1:n_alleles
        allele_len = read(io, UInt32)
        alleles_bytes = read(io, allele_len)
        alleles[i] = String(alleles_bytes)
    end

    if compression == 0 && layout == 1
        geno_block_size = 6 * n_samples
    else
        geno_block_size = read(io, UInt32)
    end
    geno_offset = position(io)
    next_var_offset = geno_offset + geno_block_size

    Variant(offset, geno_offset, next_var_offset, geno_block_size, n_samples,
        varid, rsid, chrom, pos, n_alleles, alleles, nothing)
end

function Variant(b::Bgen, offset::Integer)
    h = b.header
    if offset >= b.fsize
        @error "reached end of file"
    end
    Variant(b.io, offset, h.compression, h.layout, h.n_samples)
end

@inline n_samples(v::Variant)::Int = v.n_samples
@inline varid(v::Variant) = v.varid
@inline rsid(v::Variant) = v.rsid
@inline chrom(v::Variant) = v.chrom
@inline pos(v::Variant)::Int = v.pos
@inline n_alleles(v::Variant)::Int = v.n_alleles
@inline alleles(v::Variant) = v.alleles

# The following functions are valid only after calling `probabilities!()`
# or `minor_allele_dosage!()`
@inline phased(v::Variant) = v.genotypes.preamble.phased
@inline min_ploidy(v::Variant) = v.genotypes.preamble.min_ploidy
@inline max_ploidy(v::Variant) = v.genotypes.preamble.max_ploidy
@inline ploidy(v::Variant) = v.genotypes.preamble.ploidy
@inline bit_depth(v::Variant) = v.genotypes.preamble.bit_depth
@inline missings(v::Variant) = v.genotypes.preamble.missings

# The below are valid after calling `minor_allele_dosage!()`
@inline function minor_allele(v::Variant)
    midx = v.genotypes.minor_idx
    if midx == 0
        @error "`minor_allele_dosage!()` must be called before `minor_allele()`"
    else
        v.alleles[v.genotypes.minor_idx]
    end
end

@inline function major_allele(v::Variant)
    midx = v.genotypes.minor_idx
    if midx == 0
        @error "`minor_allele_dosage!()` must be called before `minor_allele()`"
    else
        v.alleles[3 - midx]
    end
end

"""
    destroy_genotypes!(v::Variant)
Destroy any parsed genotype information.
"""
function clear!(v::Variant)
    v.genotypes = nothing
    return
end


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer

import SimpleDirectMediaLayer.LoadBMP

SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLEBUFFERS, 16)
SDL2.GL_SetAttribute(SDL2.GL_MULTISAMPLESAMPLES, 16)

SDL2.init()

win = SDL2.CreateWindow("Hello World!", Int32(100), Int32(100), Int32(800), Int32(600),
    UInt32(SDL2.WINDOW_SHOWN))
SDL2.SetWindowResizable(win,true)

renderer = SDL2.CreateRenderer(win, Int32(-1),
    UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))

import Base.unsafe_convert
unsafe_convert(::Type{Ptr{SDL2.RWops}}, s::String) = SDL2.RWFromFile(s, "rb")

LoadBMP(src::String) = SDL2.LoadBMP_RW(src,Int32(1))

bkg = SDL2.Color(200, 200, 200, 255)

# Create text
font = TTF_OpenFont(joinpath(@__DIR__,"../../assets/fonts/FiraCode/ttf/FiraCode-Regular.ttf"), 14)
#txt = "@BinDeps.install Dict([ (:glib, :libglib) ])"
txt = "hi"
text = TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
tex = SDL2.CreateTextureFromSurface(renderer,text)

fx,fy = Int[1], Int[1]
TTF_SizeText(font, txt, pointer(fx), pointer(fy))
fx,fy = fx[1],fy[1]

#img = SDL2.LoadBMP("LB2951.jpg")
#tex = SDL2.CreateTextureFromSurface(ren, img)
#SDL2.FreeSurface(img)

for i = 1:200
    x,y = Int[1], Int[1]
    SDL2.PumpEvents()
    SDL2.GetMouseState(pointer(x), pointer(y))

    # Set render color to red ( background will be rendered in this color )
    SDL2.SetRenderDrawColor(renderer, 200, 200, 200, 255)
    SDL2.RenderClear(renderer)

    SDL2.SetRenderDrawColor(renderer, 20, 50, 105, 255)
    SDL2.RenderDrawLine(renderer,0,0,800,600)

    rect = SDL2.Rect(1,1,200,200)
    SDL2.RenderFillRect(renderer, pointer_from_objref(rect) )
    SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(SDL2.Rect(x[1],y[1],fx,fy)))

    SDL2.RenderPresent(renderer)
    sleep(0.001)
end
SDL2.Quit()


text:
"""
```
compute_system(m; tvis::Bool = false, verbose = :high)
compute_system_helper(m; tvis::Bool = false, verbose = :high)
```

Given the current model parameters, compute the state-space system
corresponding to model `m`. Returns a `System` or `RegimeSwitchingSystem` object.
The keyword `apply_altpolicy` indicates whether the state-space system
should reflect an alternative policy, and the keyword `tvis` indicates
whether the state-space system involves time-varying information sets.
To use `tvis = true`, at least the setting `:tvis_information_set` must
exist. See `?DSGE.compute_tvis_system` for more information about
computing state-space systems with time-varying information sets.

compute_system just runs compute_system_helper (which actually computes
the system) and then in the case of imperfect but positive credibility,
adjusts the anticipated observables and pseudo-observables' measurement
equations.
"""
function compute_system(m::AbstractDSGEModel{T}; tvis::Bool = false,
                        set_regime_eqcond_info!::Union{Function, Nothing} = nothing,
                        verbose::Symbol = :high) where {T <: Real}

    # do NOT delete the second part of the Boolean for apply_altpolicy; it is needed for
    # the scenarios code, which continues to use :alternative_policy as a Setting
    # AND for implementing imperfect awareness via uncertain_altpolicy
    apply_altpolicy = haskey(get_settings(m), :regime_eqcond_info) || haskey(get_settings(m), :alternative_policy)

    # When regime-switching, we need to do some extra work
    if haskey(get_settings(m), :regime_switching) && get_setting(m, :regime_switching)
        # Grab these settings
        has_uncertain_altpolicy = haskey(get_settings(m), :uncertain_altpolicy)
        has_uncertain_temporary_altpolicy = haskey(get_settings(m), :uncertain_temporary_altpolicy)
        has_regime_eqcond_info = haskey(get_settings(m), :regime_eqcond_info)
        uncertain_altpolicy = has_uncertain_altpolicy && get_setting(m, :uncertain_altpolicy)
        uncertain_temporary_altpolicy = has_uncertain_temporary_altpolicy && get_setting(m, :uncertain_temporary_altpolicy)
        is_gensys2 = haskey(get_settings(m), :gensys2) && get_setting(m, :gensys2)

        # Grab regime info dictionary, if one exists
        regime_eqcond_info = has_regime_eqcond_info ? get_setting(m, :regime_eqcond_info) : Dict{Int, EqcondEntry}()

        # Ensure :alternative_policies exists as a Setting
        @assert (!uncertain_altpolicy && !uncertain_temporary_altpolicy) ||
            haskey(get_settings(m),
                   :alternative_policies) "To use uncertain_altpolicy/uncertain_temporary_altpolicy, the setting :alternative_policies must be set"

        # # If uncertain_temporary_altpolicy is false, want to make sure temp altpol period is treated as certain
        # # in the case that uncertain_altpolicy is true (hence imperfect awareness is one
        if has_uncertain_temporary_altpolicy && !uncertain_temporary_altpolicy && has_regime_eqcond_info
            for reg in keys(regime_eqcond_info)
                if regime_eqcond_info[reg].alternative_policy.key in get_setting(m, :temporary_altpolicy_names)
                    if ismissing(regime_eqcond_info[reg].weights) && haskey(get_settings(m), :alternative_policies)
                        # We can infer the correct weight vector from alternative_policies.
                        # If :alternative_policies is not a setting and weights is missing,
                        # then nothing will happen b/c either the lack of an :alternative_policies
                        # will trigger an error, or the weights won't be needed.
                        altpol_vec = zeros(length(get_setting(m, :alternative_policies)) + 1)
                        altpol_vec[1] = 1.0
                        regime_eqcond_info[reg] = EqcondEntry(regime_eqcond_info[reg].alternative_policy, altpol_vec)
                    elseif !ismissing(regime_eqcond_info[reg].weights)
                        # if the weights vector is not missing, then we can infer the correct length from it
                        altpol_vec = zeros(length(regime_eqcond_info[reg].weights))
                        altpol_vec[1] = 1.0
                        regime_eqcond_info[reg].weights = altpol_vec
                    end
                end
            end
        end

        # Same for uncertain_altpolicy (note: this step is unnecessary for compute_system_helper,
        ## but it is helpful for combining historical and alternative policies with
        ## the right weights, e.g. in the time_varying_credibility.jl test)
        if has_uncertain_altpolicy && !uncertain_altpolicy && has_regime_eqcond_info
            for reg in keys(regime_eqcond_info)
                if regime_eqcond_info[reg].alternative_policy.key == alternative_policy(m).key # TODO: update this to allow for a Vector{Symbol}
                    # Same logic as previous block for these if conditions                     # as with temporary policies
                    if ismissing(regime_eqcond_info[reg].weights) && haskey(get_settings(m), :alternative_policies)
                        altpol_vec = zeros(length(get_setting(m, :alternative_policies)))
                        altpol_vec[1] = 1.0
                        regime_eqcond_info[reg] = EqcondEntry(regime_eqcond_info[reg].alternative_policy, altpol_vec)
                    elseif !ismissing(regime_eqcond_info[reg].weights)
                        altpol_vec = zeros(length(regime_eqcond_info[reg].weights))
                        altpol_vec[1] = 1.0
                        regime_eqcond_info[reg].weights = altpol_vec
                    end
                end
            end
        end

        system_main = if uncertain_altpolicy
            # Handle this special case on its own
            compute_uncertain_altpolicy_system_helper(m; tvis = tvis, verbose = verbose)
        else
            compute_system_helper(m; tvis = tvis, verbose = verbose)
        end
    else
        system_main = compute_system_helper(m; tvis = tvis, verbose = verbose)
    end

    # If correcting measurement eqs for anticipated (pseudo) observables is unnecessary
    # (eg. running Taylor or no regime switching or no uncertainty in temp altpol or altpolicy or
    # either perfect or zero credibility - invariant perfect or zero credibility in the case of time varying),
    # then return system now.
    # The !apply_altpolicy check may be problematic after refactoring altpolicy.
    if !apply_altpolicy || !haskey(get_settings(m), :regime_switching) || !get_setting(m, :regime_switching) ||
        !has_regime_eqcond_info || # if regime_eqcond_info is not defined, then no alt policies occur
        (has_uncertain_temporary_altpolicy && !uncertain_temporary_altpolicy && has_uncertain_altpolicy && !uncertain_altpolicy) ||
        (!has_uncertain_temporary_altpolicy && !has_uncertain_altpolicy)

        if haskey(get_settings(m), :regime_switching) && get_setting(m, :regime_switching) && !apply_altpolicy
            if has_regime_eqcond_info
                m <= Setting(:regime_eqcond_info, regime_info_copy)
            end
            m <= Setting(:gensys2, is_gensys2)
        end
    end

    return system_main
end

function compute_system_helper(m::AbstractDSGEModel{T}; tvis::Bool = false, verbose::Symbol = :high) where {T <: Real}

    if tvis
        @assert haskey(get_settings(m), :tvis_information_set) "The setting :tvis_information_set is not defined"
        n_tvis = haskey(get_settings(m), :tvis_regime_eqcond_info) ? length(get_setting(m, :tvis_regime_eqcond_info)) : 1
        if n_tvis == 1 && haskey(get_settings(m), :tvis_regime_eqcond_info)
            if haskey(get_settings(m), :regime_eqcond_info)
                if get_setting(m, :tvis_regime_eqcond_info)[1] != get_setting(m, :regime_eqcond_info)
                    warn_str = "The dictionary of functions in the Setting :regime_eqcond_info does not match the one specified " *
                    "by the length-one Setting :tvis_regime_eqcond_info. Replacing :regime_eqcond_info with the dictionary " *
                    "of functions contained in :tvis_regime_eqcond_info."
                    @warn warn_str
                    m <= Setting(:regime_eqcond_info, get_setting(m, :tvis_regime_eqcond_info)[1])
                end
            else
                m <= Setting(:regime_eqcond_info, get_setting(m, :tvis_regime_eqcond_info)[1])
            end
        end

        if n_tvis > 1 # case of n_tvis = 1 handled below to avoid constructing redundant TimeVaryingInformationSetSystem
            @assert haskey(get_settings(m), :tvis_select_system) "The setting :tvis_select_system is not defined"
            tvis_sys = compute_tvis_system(m; verbose = verbose)
            transition_eqns = Transition{T}[tvis_sys[select, reg, :transition] for (reg, select) in enumerate(tvis_sys[:select])]

            return RegimeSwitchingSystem(transition_eqns, tvis_sys[:measurements], tvis_sys[:pseudo_measurements])
        end
    end

    solution_method  = get_setting(m, :solution_method)
    regime_switching = haskey(get_settings(m), :regime_switching) && get_setting(m, :regime_switching)

    # Solve model
    if regime_switching
        n_regimes      = (regime_switching && haskey(get_settings(m), :n_regimes)) ? get_setting(m, :n_regimes) : 1
        n_hist_regimes = (regime_switching && haskey(get_settings(m), :n_hist_regimes)) ? get_setting(m, :n_hist_regimes) : 1

        if solution_method == :gensys
            # Determine which regimes should use gensys2/gensys
            gensys_regimes, gensys2_regimes = compute_gensys_gensys2_regimes(m)

            # Solve!
            TTTs, RRRs, CCCs = solve(m; regime_switching = regime_switching,
                                     regimes = collect(1:n_regimes),
                                     gensys_regimes = gensys_regimes,
                                     gensys2_regimes = gensys2_regimes,
                                     verbose = verbose)

            if length(gensys2_regimes) > 1 && haskey(get_settings(m), :use_forward_expectations_memo) &&
               get_setting(m, :use_forward_expectations_memo)
                warn_str = "There are multiple gensys2 regimes. Double check that this should be the case." *
                    " If so, note that the current implementation of use_forward_expectations_memo" *
                    " may not be robust to having 2+ gensys2_regimes. In particular, it is assumed" *
                    " that we only need to compute the memo once, but if the TVIS infoset" *
                    " implies multiple gensys2_regimes, it is likely that you need to recompute the memo" *
                    " to ensure the correct sequence of products of the TTTs is calculated"
                @warn warn_str
            end

            memo_gensys2_regimes = isempty(gensys2_regimes) ? (1:1) : gensys2_regimes[1]

            return RegimeSwitchingSystem(m, TTTs, RRRs, CCCs, n_regimes, tvis;
                                         gensys2_regimes = memo_gensys2_regimes) # handles formation of transition and measurement equations
        else
            error("Regime switching with solution algorithms other than gensys has not been implemented.")
        end
    else
        if solution_method == :gensys

            TTT, RRR, CCC = solve(m; verbose = verbose)
            transition_equation = Transition(TTT, RRR, CCC)

            # Solve measurement equation
            measurement_equation = measurement(m, TTT, RRR, CCC)

        elseif solution_method == :klein
            # Unpacking the method from solve to hang on to TTT_jump
            if m.spec == "het_dsge"
                TTT_jump, TTT_state, eu = klein(m)
            else
                TTT_jump, TTT_state, eu = klein(m)
            end
            if eu==-1
                throw(KleinError())
            end

            TTT, RRR = klein_transition_matrices(m, TTT_state, TTT_jump)
            CCC = zeros(n_model_states(m))

            if m.spec == "real_bond_mkup"
                GDPeqn = construct_GDPeqn(m, TTT_jump)
                TTT, RRR, CCC = augment_states(m, TTT, TTT_jump, RRR, CCC, GDPeqn)
                # Measurement (needs the additional TTT_jump argument)
                measurement_equation = measurement(m, TTT, TTT_jump, RRR, CCC, GDPeqn)
            elseif m.spec == "het_dsge" || m.spec == "rep_dsge"
                TTT, RRR, CCC = augment_states(m, TTT, RRR, CCC)
                measurement_equation = measurement(m, TTT, RRR, CCC)
            else
                TTT, RRR, CCC        = augment_states(m, TTT, RRR, CCC)
                measurement_equation = measurement(m, TTT, RRR, CCC)
            end

            transition_equation = Transition(TTT, RRR, CCC)

        else
            throw("solution_method provided does not exist.")
        end

        type_tuple = (typeof(m), Matrix{T}, Matrix{T}, Vector{T})
        if hasmethod(pseudo_measurement, type_tuple)
            # Solve pseudo-measurement equation
            pseudo_measurement_equation = pseudo_measurement(m, TTT, RRR, CCC)
            return System(transition_equation, measurement_equation, pseudo_measurement_equation)
        else
            return System(transition_equation, measurement_equation)
        end
    end
end

function RegimeSwitchingSystem(m::AbstractDSGEModel{T}, TTTs::Vector{<: AbstractMatrix{T}}, RRRs::Vector{<: AbstractMatrix{T}},
                               CCCs::Vector{<: AbstractVector{T}}, n_regimes::Int, tvis::Bool;
                               gensys2_regimes::UnitRange{Int64} = 1:1) where {T <: Real}

    transition_equations = Vector{Transition{T}}(undef, n_regimes)
    for i = 1:n_regimes
        transition_equations[i] = Transition(TTTs[i], RRRs[i], CCCs[i])
    end

    # Compute memo for forward expectations if requested
    use_fwd_exp_sum = haskey(get_settings(m), :use_forward_expected_sum_memo) && get_setting(m, :use_forward_expected_sum_memo)
    use_fwd_exp     = haskey(get_settings(m), :use_forward_expectations_memo) && get_setting(m, :use_forward_expectations_memo)
    perm_pol_reg    = haskey(get_settings(m), :memo_permanent_policy_regime) ?
        get_setting(m, :memo_permanent_policy_regime) : gensys2_regimes[end] # make a guess that the permanent policy is gensys2_regimes[end]
    memo = if use_fwd_exp_sum                                                # but this guess may not work
        # TODO: update the hard-coded 40 to allow for other cases, but currently, the only use case
        # requires power up to 40.
        ForwardMultipleExpectationsMemo(TTTs, gensys2_regimes[1] + 1, gensys2_regimes[1] + 2, perm_pol_reg, perm_pol_reg, 1, 40)
    elseif use_fwd_exp
        ForwardMultipleExpectationsMemo(TTTs, gensys2_regimes[1] + 1, gensys2_regimes[1] + 2, perm_pol_reg,
                                        perm_pol_reg, 1, n_mon_anticipated_shocks(m))
    else
        nothing
    end

    # Infer which measurement and pseudo-measurement equations to use
    model_type = typeof(m)
    type_tuple = (model_type, Vector{Matrix{T}}, Vector{Matrix{T}}, Vector{Vector{T}})
    has_pseudo = hasmethod(pseudo_measurement, type_tuple) ||
        hasmethod(pseudo_measurement, (model_type, Matrix{T}, Matrix{T}, Vector{T}))

    if tvis
        if hasmethod(measurement, type_tuple)
            measurement_equations =
                measurement(m, TTTs, RRRs, CCCs; information_set = get_setting(m, :tvis_information_set)[reg], memo = memo)
        else
            empty_meas_eqn = haskey(get_settings(m), :empty_measurement_equation) ?
                get_setting(m, :empty_measurement_equation) : falses(n_regimes)

            measurement_equations = Vector{Measurement{T}}(undef, n_regimes)
            for reg in 1:n_regimes
                measurement_equations[reg] = if empty_meas_eqn[reg]
                    Measurement(Matrix{T}(undef, 0, 0), Vector{T}(undef, 0),
                                Matrix{T}(undef, 0, 0), Matrix{T}(undef, 0, 0))
                else
                    measurement(m, TTTs[reg], RRRs[reg], CCCs[reg], reg = reg,
                                TTTs = TTTs, CCCs = CCCs,
                                information_set = get_setting(m, :tvis_information_set)[reg], memo = memo)
                end
            end
        end

        if hasmethod(pseudo_measurement, type_tuple)
            pseudo_measurement_equations =
                pseudo_measurement(m, TTTs, RRRs, CCCs;
                                   information_set = get_setting(m, :tvis_information_set), memo = memo)
        elseif hasmethod(pseudo_measurement, (typeof(m), Matrix{T}, Matrix{T}, Vector{T}))
            empty_pseudo_meas_eqn = haskey(get_settings(m), :empty_pseudo_measurement_equation) ?
                get_setting(m, :empty_pseudo_measurement_equation) : falses(n_regimes)

            pseudo_measurement_equations = Vector{PseudoMeasurement{T}}(undef, n_regimes)
            for reg in 1:n_regimes
                pseudo_measurement_equations[reg] = if empty_pseudo_meas_eqn[reg]
                    PseudoMeasurement(Matrix{T}(undef, 0, 0), Vector{T}(undef, 0))
                else
                    pseudo_measurement(m, TTTs[reg], RRRs[reg], CCCs[reg];
                                       reg = reg, TTTs = TTTs, CCCs = CCCs,
                                       information_set = get_setting(m, :tvis_information_set)[reg], memo = memo)
                end
            end
        end
    else
        if hasmethod(measurement, type_tuple)
            measurement_equations = measurement(m, TTTs, RRRs, CCCs; memo = memo)
        else
            empty_meas_eqn = haskey(get_settings(m), :empty_measurement_equation) ?
                get_setting(m, :empty_measurement_equation) : falses(n_regimes)

            measurement_equations = Vector{Measurement{T}}(undef, n_regimes)
            for reg in 1:n_regimes
                measurement_equations[reg] = if empty_meas_eqn[reg]
                    Measurement(Matrix{T}(undef, 0, 0), Vector{T}(undef, 0),
                                Matrix{T}(undef, 0, 0), Matrix{T}(undef, 0, 0))
                else
                    measurement(m, TTTs[reg], RRRs[reg], CCCs[reg]; reg = reg, memo = memo)
                end
            end
        end

        if hasmethod(pseudo_measurement, type_tuple)
            pseudo_measurement_equations = pseudo_measurement(m, TTTs, RRRs, CCCs; memo = memo)
        elseif hasmethod(pseudo_measurement, (typeof(m), Matrix{T}, Matrix{T}, Vector{T}))
            empty_pseudo_meas_eqn = haskey(get_settings(m), :empty_pseudo_measurement_equation) ?
                get_setting(m, :empty_pseudo_measurement_equation) : falses(n_regimes)

            pseudo_measurement_equations = Vector{PseudoMeasurement{T}}(undef, n_regimes)
            for reg in 1:n_regimes
                pseudo_measurement_equations[reg] = if empty_pseudo_meas_eqn[reg]
                    PseudoMeasurement(Matrix{T}(undef, 0, 0), Vector{T}(undef, 0))
                else
                    pseudo_measurement(m, TTTs[reg], RRRs[reg], CCCs[reg]; reg = reg, memo = memo)
                end
            end
        end
    end

    if has_pseudo
        return RegimeSwitchingSystem(transition_equations, measurement_equations, pseudo_measurement_equations)
    else
        return RegimeSwitchingSystem(transition_equations, measurement_equations)
    end
end

"""
```
compute_system(m::PoolModel{T})
```

Given the current model parameters, compute the state-space system
corresponding to the PoolModel model `m`.

Outputs

```
Φ: state transition function
Ψ: likelihood function, given weights on underlying models (the states) and predictive densities
F_ϵ: structural shock distribution
F_u: likelihood function measurement error distribution
F_λ: initial distribution of λ for state transition function
"""
function compute_system(m::PoolModel{T}; tvis::Bool = false, verbose::Symbol = :high) where {T <: Real}
    Φ, F_ϵ, F_λ = transition(m)
    Ψ, F_u = measurement(m)
    return Φ, Ψ, F_ϵ, F_u, F_λ
end

"""
```
compute_uncertain_altpolicy_system_helper(m::AbstractDSGEModel{T}; tvis::Bool = false, verbose::Symbol = :high) where {T <: Real}
```

calculates a state-space system with imperfect awareness and possibly multi-period alternative policies.
This calculation is performed separately to ensure computational efficiency.

Note that the endogenous states across all policies (implemented and alternative) must be the same, and any augmentations
to the state space after calling `gensys`/`gensys2` should be done within the default `augment_states` function. The user
can typically add if-else conditions within `augment_states` to handle state space augmentations for different alternative policy.
"""
function compute_uncertain_altpolicy_system_helper(m::AbstractDSGEModel{T}; tvis::Bool = false,
                                                   verbose::Symbol = :high) where {T <: Real}

    # Set up
    solution_method = get_setting(m, :solution_method)
    @assert solution_method == :gensys "No solution algorithm except gensys is allowed with multi-period alternative policies"

    assert_str = "To use imperfect awareness, " *
        "the setting :replace_eqcond must be true, and " *
        "either regime_eqcond_info or tvis_regime_eqcond_info is a setting"
    @assert haskey(get_settings(m), :replace_eqcond) &&
        get_setting(m, :replace_eqcond) && (haskey(get_settings(m), :regime_eqcond_info) ||
                                            haskey(get_settings(m), :tvis_regime_eqcond_info)) assert_str

    n_regimes = get_setting(m, :n_regimes)

    if tvis
        @assert haskey(get_settings(m), :tvis_information_set) "The setting :tvis_information_set is not defined"
        n_tvis = haskey(get_settings(m), :tvis_regime_eqcond_info) ? length(get_setting(m, :tvis_regime_eqcond_info)) : 1
        if n_tvis == 1 && haskey(get_settings(m), :tvis_regime_eqcond_info)
            if haskey(get_settings(m), :regime_eqcond_info)
                if get_setting(m, :tvis_regime_eqcond_info)[1] != get_setting(m, :regime_eqcond_info)
                    warn_str = "The dictionary of functions in the Setting :regime_eqcond_info does not match the one specified " *
                    "by the length-one Setting :tvis_regime_eqcond_info. Replacing :regime_eqcond_info with the dictionary " *
                    "of functions contained in :tvis_regime_eqcond_info."
                    @warn warn_str
                    m <= Setting(:regime_eqcond_info, get_setting(m, :tvis_regime_eqcond_info)[1])
                end
            else
                m <= Setting(:regime_eqcond_info, get_setting(m, :tvis_regime_eqcond_info)[1])
            end
        end

        if n_tvis > 1 # case of n_tvis = 1 handled below to avoid constructing redundant TimeVaryingInformationSetSystem
            assert_str =  "Multiple equilibrium conditions with time-varying information sets and multi-period " *
                "alternative policies has not been implemented"
            @assert false assert_str
        end
    end

    # First calculate the perfectly credible transition matrices for the policy actually being implemented
    is_altpol = [isa(x, AltPolicy) for x in get_setting(m, :alternative_policies)]
    if all(is_altpol)
        system_perfect_cred_totpolicies = perfect_cred_altpolicy_systems(m; verbose = verbose)
    else
        system_perfect_cred_totpolicies = perfect_cred_multiperiod_altpolicy_systems(m, is_altpol; verbose = verbose)
    end
    n_tot_policies = length(system_perfect_cred_totpolicies)

    # Checks if pseudo measurement is required
    type_tuple = (typeof(m), Vector{Matrix{T}}, Vector{Matrix{T}}, Vector{Vector{T}})
    has_pseudo = hasmethod(pseudo_measurement, type_tuple) ||
    hasmethod(pseudo_measurement, (typeof(m), Matrix{T}, Matrix{T}, Vector{T}))

    has_fwd_looking_obs = haskey(get_settings(m), :forward_looking_observables)
    if has_pseudo
        has_fwd_looking_pseudo = haskey(get_settings(m), :forward_looking_pseudo_observables)
    end

    # Form measurement and pseudo-measurement equations for imperfect awareness state space system
    gensys_regimes, gensys2_regimes = compute_gensys_gensys2_regimes(m)
    TTTs, RRRs, CCCs = solve_uncertain_altpolicy(m, system_perfect_cred_totpolicies, is_altpol;
                                                 gensys_regimes = gensys_regimes, gensys2_regimes = gensys2_regimes,
                                                 regimes = collect(1:n_regimes), verbose = verbose)

    if length(gensys2_regimes) > 1 && haskey(get_settings(m), :use_forward_expectations_memo) &&
        get_setting(m, :use_forward_expectations_memo)
        warn_str = "There are multiple gensys2 regimes. Double check that this should be the case." *
            " If so, note that the current implementation of use_forward_expectations_memo" *
            " may not be robust to having 2+ gensys2_regimes. In particular, it is assumed" *
            " that we only need to compute the memo once, but if the TVIS infoset" *
            " implies multiple gensys2_regimes, it is likely that you need to recompute the memo" *
            " to ensure the correct sequence of products of the TTTs is calculated"
        @warn warn_str
    end
    memo_gensys2_regimes = isempty(gensys2_regimes) ? (1:1) : gensys2_regimes[1]
    system_main = RegimeSwitchingSystem(m, TTTs, RRRs, CCCs, n_regimes, tvis;
                                        gensys2_regimes = memo_gensys2_regimes)

    # Correct the measurement equations for anticipated observables via weighted average
    regime_eqcond_info = get_setting(m, :regime_eqcond_info)
    which_is_system = vcat(false, is_altpol) # need to add a false to the start of is_altpol for implemented policy
    empty_meas_eqn = haskey(get_settings(m), :empty_measurement_equation) ?
        get_setting(m, :empty_measurement_equation) : falses(n_regimes)
    empty_pseudo_meas_eqn = haskey(get_settings(m), :empty_pseudo_measurement_equation) ?
        get_setting(m, :empty_pseudo_measurement_equation) : falses(n_regimes)
    for reg in sort!(collect(keys(regime_eqcond_info)))

        new_wt = regime_eqcond_info[reg].weights

        if !empty_meas_eqn[reg]
            if has_fwd_looking_obs
                for k in get_setting(m, :forward_looking_observables)
                    system_main.measurements[reg][:ZZ][m.observables[k], :] =
                        sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:ZZ][m.observables[k], :] :
                                           system_perfect_cred_totpolicies[i].measurements[min(length(system_perfect_cred_totpolicies[i].measurements), reg)][:ZZ][m.observables[k], :]) for i in 1:length(new_wt)])
                    system_main.measurements[reg][:DD][m.observables[k]] =
                        sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:DD][m.observables[k]] :
                                           system_perfect_cred_totpolicies[i].measurements[min(length(system_perfect_cred_totpolicies[i].measurements), reg)][:DD][m.observables[k]]) for i in 1:length(new_wt)])
                end
            else
                system_main.measurements[reg][:ZZ] .=
                    sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:ZZ] :
                                       system_perfect_cred_totpolicies[i].measurements[min(length(system_perfect_cred_totpolicies[i].measurements), reg)][:ZZ]) for i in 1:length(new_wt)])
                system_main.measurements[reg][:DD] .=
                    sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:DD] :
                                       system_perfect_cred_totpolicies[i].measurements[min(length(system_perfect_cred_totpolicies[i].measurements), reg)][:DD]) for i in 1:length(new_wt)])
            end
        end

        if has_pseudo && !empty_pseudo_meas_eqn[reg]
            if has_fwd_looking_pseudo
                for k in get_setting(m, :forward_looking_pseudo_observables)
                    system_main.pseudo_measurements[reg][:ZZ_pseudo][m.pseudo_observables[k], :] =
                        sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:ZZ_pseudo][m.pseudo_observables[k], :] :
                                           system_perfect_cred_totpolicies[i].pseudo_measurements[min(length(system_perfect_cred_totpolicies[i].pseudo_measurements), reg)][:ZZ_pseudo][m.pseudo_observables[k], :])
                             for i in 1:length(new_wt)])
                    system_main.pseudo_measurements[reg][:DD_pseudo][m.pseudo_observables[k]] =
                        sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:DD_pseudo][m.pseudo_observables[k]] :
                                           system_perfect_cred_totpolicies[i].pseudo_measurements[min(length(system_perfect_cred_totpolicies[i].pseudo_measurements), reg)][:DD_pseudo][m.pseudo_observables[k]])
                             for i in 1:length(new_wt)])
                end
            else
                system_main.pseudo_measurements[reg][:ZZ_pseudo] .=
                    sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:ZZ_pseudo] :
                                       system_perfect_cred_totpolicies[i].pseudo_measurements[min(length(system_perfect_cred_totpolicies[i].pseudo_measurements), reg)][:ZZ_pseudo])
                         for i in 1:length(new_wt)])
                system_main.pseudo_measurements[reg][:DD_pseudo] .=
                        sum([new_wt[i] .* (which_is_system[i] ? system_perfect_cred_totpolicies[i][:DD_pseudo] :
                                           system_perfect_cred_totpolicies[i].pseudo_measurements[min(length(system_perfect_cred_totpolicies[i].pseudo_measurements), reg)][:DD_pseudo])
                             for i in 1:length(new_wt)])
            end
        end
    end

    return system_main
end

function perfect_cred_altpolicy_systems(m::AbstractDSGEModel{T}; verbose::Symbol = :high) where {T <: Real}

    # Set up
    m <= Setting(:uncertain_altpolicy, false)
    if haskey(get_settings(m), :uncertain_temporary_altpolicy) && get_setting(m, :uncertain_temporary_altpolicy)
        m <= Setting(:uncertain_temporary_altpolicy, false)
        orig_temp_altpol = true
    else
        orig_temp_altpol = false
    end
    n_totpol = length(get_setting(m, :alternative_policies)) + 1 # + 1 for the implemented policy
    sys_totpolicies = Vector{Union{RegimeSwitchingSystem, System}}(undef, n_totpol)

    # Save the following elements.
    # Note deepcopy is necessary for regime_eqcond_info b/c weights will be updated to [1., 0., ...]
    # when we compute the perfectly credible implemented policy
    orig_n_regimes          = get_setting(m, :n_regimes)
    orig_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
    orig_altpol             = haskey(get_settings(m), :alternative_policy) ? get_setting(m, :alternative_policy) : nothing
    is_gensys2              = haskey(get_settings(m), :gensys2) && get_setting(m, :gensys2)

    ## Now calculate the perfectly credible transition matrices for all alternative policies
    sys_totpolicies[1] = compute_system(m; tvis = haskey(get_settings(m), :tvis_information_set), verbose = verbose)

    # m <= Setting(:regime_eqcond_info, Dict{Int64, EqcondEntry}()) # TODO: maybe also delete this setting
    delete!(get_settings(m), :regime_eqcond_info)
    delete!(get_settings(m), :alternative_policy)   # does nothing if alternative_policy is not a key in get_settings(m)

    for i in 2:n_totpol # want to populate TTTs, RRRs, and CCCs, and there are n_totpol total policies
        new_altpol = get_setting(m, :alternative_policies)[i - 1] # decrement by 1 b/c :alternative_policies only includes
        # We assume that the user only wants the permanent altpolicy system
        # AND there is no parameter regime-switching that affects the TTT matrix or CCC vector.
        # If there is parameter regime-switching or other kinds of regime-switching
        # that affect the TTT matrix/CCC vector in the alternative policies which
        # people are using to form expectations, then the user needs to pass a
        # MultiPeriodAltPolicy type (and set the gensys2 field to false if
        # calling gensys2 is unnecessary)
        m <= Setting(:regime_switching, false) # turn off regime-switching
        delete!(get_settings(m), :regime_eqcond_info)
        if new_altpol.key != :default_policy # in this case, we can save some extra time
            m <= Setting(:alternative_policy, new_altpol)
        end
        sys_totpolicies[i] = compute_system(m; verbose = verbose)
        m <= Setting(:regime_switching, true) # turn off regime-switching
        if new_altpol.key != :default_policy
            delete!(get_settings(m), :alternative_policy)
        end
    end

    # Now add original settings back
    m <= Setting(:n_regimes, orig_n_regimes)
    m <= Setting(:regime_eqcond_info, orig_regime_eqcond_info)
    m <= Setting(:gensys2, is_gensys2)
    m <= Setting(:uncertain_altpolicy, true)
    if !isnothing(orig_altpol)
        m <= Setting(:alternative_policy, orig_altpol)
    end
    if haskey(get_settings(m), :uncertain_temporary_altpolicy) && orig_temp_altpol
        m <= Setting(:uncertain_temporary_altpolicy, orig_temp_altpol)
    end

    return sys_totpolicies
end

function perfect_cred_multiperiod_altpolicy_systems(m::AbstractDSGEModel{T}, is_altpol::Vector{Bool};
                                                    verbose::Symbol = :high) where {T <: Real}

    # Set up
    m        <= Setting(:uncertain_altpolicy,           false)
    m        <= Setting(:uncertain_temporary_altpolicy, false)
    n_totpol = length(get_setting(m, :alternative_policies)) + 1 # + 1 for the implemented policy
    sys_totpolicies = Vector{Union{System, RegimeSwitchingSystem}}(undef, n_totpol)

    # Save the following elements.
    # Note deepcopy is necessary for regime_eqcond_info b/c weights will be updated to [1., 0., ...]
    # when we compute the perfectly credible implemented policy
    orig_n_regimes          = get_setting(m, :n_regimes)
    orig_regime_eqcond_info = deepcopy(get_setting(m, :regime_eqcond_info))
    orig_altpol             = haskey(get_settings(m), :alternative_policy) ? get_setting(m, :alternative_policy) : nothing
    orig_tvis_infoset       = haskey(get_settings(m), :tvis_information_set) ? get_setting(m, :tvis_information_set) : nothing
    orig_temp_altpol_names  = haskey(get_settings(m), :temporary_altpolicy_names) ? get_setting(m, :temporary_altpolicy_names) : nothing
    orig_temp_altpol_len    = haskey(get_settings(m), :temporary_altpolicy_length) ? get_setting(m, :temporary_altpolicy_length) : nothing
    orig_perf_cred_same_reg = haskey(get_settings(m), :perfect_credibility_identical_transitions) ? get_setting(m, :perfect_credibility_identical_transitions) : nothing
    orig_iden_eqcond_reg    = haskey(get_settings(m), :identical_eqcond_regimes) ? get_setting(m, :identical_eqcond_regimes) : nothing
    is_gensys2              = haskey(get_settings(m), :gensys2) && get_setting(m, :gensys2)

    ## Now calculate the perfectly credible transition matrices for all alternative policies
    sys_totpolicies[1] = compute_system(m; tvis = haskey(get_settings(m), :tvis_information_set), verbose = verbose)

    # m <= Setting(:regime_eqcond_info, Dict{Int64, EqcondEntry}()) # TODO: maybe also delete this setting
    delete!(get_settings(m), :regime_eqcond_info)
    delete!(get_settings(m), :alternative_policy)   # does nothing if alternative_policy is not a key in get_settings(m)
    delete!(get_settings(m), :tvis_information_set) # does nothing if tvis_information_set is not a key in get_settings(m)
    delete!(get_settings(m), :temporary_altpolicy_names) # does nothing if temporary_altpolicy_names is not a key in get_settings(m)
    delete!(get_settings(m), :perfect_credibility_identical_transitions) # does nothing if temporary_altpolicy_names is not a key in get_settings(m)
    delete!(get_settings(m), :identical_eqcond_regimes) # does nothing if temporary_altpolicy_names is not a key in get_settings(m)

    for i in 2:n_totpol # want to populate TTTs, RRRs, and CCCs, and there are n_totpol total policies
        new_altpol = get_setting(m, :alternative_policies)[i - 1] # decrement by 1 b/c :alternative_policies only includes
        if is_altpol[i - 1]                                       # policies which aren't being implemented, as is the case for is_altpol
            # If AltPolicy, we assume that the user only wants the permanent altpolicy system
            # AND there is no parameter regime-switching that affects the TTT matrix or CCC vector.
            # If there is parameter regime-switching or other kinds of regime-switching
            # that affect the TTT matrix/CCC vector in the alternative policies which
            # people are using to form expectations, then the user needs to pass a
            # MultiPeriodAltPolicy type (and set the gensys2 field to false if
            # calling gensys2 is unnecessary)
            m <= Setting(:regime_switching, false) # turn off regime-switching
            delete!(get_settings(m), :regime_eqcond_info)
            if new_altpol.key != :default_policy # in this case, we can save some extra time
                m <= Setting(:alternative_policy, new_altpol)
            end
            sys_totpolicies[i] = compute_system(m; verbose = verbose)
            m <= Setting(:regime_switching, true) # turn off regime-switching
            if new_altpol.key != :default_policy
                delete!(get_settings(m), :alternative_policy)
            end
        else # Then it's a MultiPeriodAltPolicy
            m <= Setting(:n_regimes, new_altpol.n_regimes)
            m <= Setting(:regime_eqcond_info, new_altpol.regime_eqcond_info)
            m <= Setting(:gensys2, new_altpol.gensys2)
            if !isnothing(new_altpol.temporary_altpolicy_names)
                m <= Setting(:temporary_altpolicy_names, new_altpol.temporary_altpolicy_names)
            end
            if !isnothing(new_altpol.temporary_altpolicy_length)
                m <= Setting(:temporary_altpolicy_length, new_altpol.temporary_altpolicy_length)
            end
            if !isnothing(new_altpol.identical_eqcond_regimes)
                m <= Setting(:identical_eqcond_regimes, new_altpol.identical_eqcond_regimes)
            end
            if !isnothing(new_altpol.perfect_credibility_identical_transitions)
                m <= Setting(:perfect_credibility_identical_transitions, new_altpol.perfect_credibility_identical_transitions)
            end
            if !isnothing(new_altpol.infoset)
                m <= Setting(:tvis_information_set, new_altpol.infoset)
                sys_totpolicies[i] = compute_system(m; tvis = true, verbose = verbose)
            else
                sys_totpolicies[i] = compute_system(m; tvis = false, verbose = verbose)
            end
        end
    end

    # Now add original settings back
    m <= Setting(:n_regimes, orig_n_regimes)
    m <= Setting(:regime_eqcond_info, orig_regime_eqcond_info)
    m <= Setting(:gensys2, is_gensys2)
    m <= Setting(:uncertain_altpolicy,           true)
    m <= Setting(:uncertain_temporary_altpolicy, true)
    if !isnothing(orig_altpol)
        m <= Setting(:alternative_policy, orig_altpol)
    end
    if !isnothing(orig_tvis_infoset)
        m <= Setting(:tvis_information_set, orig_tvis_infoset)
    end
    if !isnothing(orig_temp_altpol_names)
        m <= Setting(:temporary_altpolicy_names, orig_temp_altpol_names)
    end
    if !isnothing(orig_temp_altpol_len)
        m <= Setting(:temporary_altpolicy_length, orig_temp_altpol_len)
    end
    if !isnothing(orig_iden_eqcond_reg)
        m <= Setting(:identical_eqcond_regimes, orig_iden_eqcond_reg)
    end
    if !isnothing(orig_perf_cred_same_reg)
        m <= Setting(:perfect_credibility_identical_transitions, orig_perf_cred_same_reg)
    end
    m <= Setting(:uncertain_altpolicy,           true)
    m <= Setting(:uncertain_temporary_altpolicy, true)

    return sys_totpolicies
end

"""
```
compute_tvis_system(m; apply_altpolicy = false, verbose = :high)
```

Given the current model parameters, compute the regime-switching
state-space system corresponding to model `m` with time-varying
information sets.

Aside from settings required for regime-switching,
additional required settings that must exist in `m` are:
- `:tvis_information_set`: Vector of `UnitRange{Int}` specifying
which regimes should be included in the information set
corresponding each regime.
- `:tvis_regime_eqcond_info`: Vector of `regime_eqcond_info`
to specify different sets of equilibrium conditions, which
generate different state space systems (usually).
- `:tvis_select_system`: Vector of `Int` to specify which
state space system to use when calculating the measurement
and pseudo measurement equations. These state space systems
correspond to different sets of equilibrium conditions (usually).
"""
function compute_tvis_system(m::AbstractDSGEModel{T}; verbose::Symbol = :high) where {T <: Real}

    # TODO: update this compute_tvis_system to compute the average over forward-looking variables
    #       when using uncertain altpol and uncertain zlb. This means that, for each
    #       vector of transition equations, we need to compute the associated average of measurement equations,
    #       which requires us to calculate the perfect credibility solutions for each vector of transition equations
    #       that are meant to have uncertainty in them.

    @assert get_setting(m, :solution_method) ==
        :gensys "Currently, the solution method must be :gensys to calculate a state-space system with time-varying information sets"

    @assert haskey(get_settings(m), :replace_eqcond) "The setting :replace_eqcond must be true to calculate a state-space system with time-varying information sets"

    # :regime_dates should have the same number of possible regimes. Any differences in eqcond
    # should be specified by tvis_regime_eqcond_info
    tvis_infoset            = get_setting(m, :tvis_information_set)
    tvis_regime_eqcond_info = get_setting(m, :tvis_regime_eqcond_info)
    tvis_select             = get_setting(m, :tvis_select_system)
    regime_switching        = get_setting(m, :regime_switching)

    n_tvis         = length(tvis_regime_eqcond_info)
    n_regimes      = regime_switching && haskey(get_settings(m), :n_regimes) ? get_setting(m, :n_regimes) : 1
    n_hist_regimes = regime_switching && haskey(get_settings(m), :n_hist_regimes) ? get_setting(m, :n_hist_regimes) : 1

    apply_altpolicy = any(.!isempty.(tvis_regime_eqcond_info)) ||
        haskey(get_settings(m), :regime_eqcond_info) || (haskey(get_settings(m), :alternative_policy) &&
        get_setting(m, :alternative_policy).key != :historical)

    # Solve model
    transitions = Vector{Vector{Transition{T}}}(undef, n_tvis)
    TTTs_vec    = Vector{Vector{Matrix{T}}}(undef, n_tvis)
    RRRs_vec    = Vector{Vector{Matrix{T}}}(undef, n_tvis)
    CCCs_vec    = Vector{Vector{Vector{T}}}(undef, n_tvis)

    for (i, regime_eqcond_info) in enumerate(tvis_regime_eqcond_info) # For each set of equilibrium conditions,
        # Update regime_eqcond_info and compute implied gensys/gensys2 regimes
        if apply_altpolicy
            m <= Setting(:regime_eqcond_info, regime_eqcond_info) # calculate the implied regime-switching system
            first_gensys2_regime = minimum(collect(keys(get_setting(m, :regime_eqcond_info))))
        else
            first_gensys2_regime = n_hist_regimes + 1
        end
        last_gensys2_regime = haskey(get_settings(m), :temporary_altpolicy_length) ?
            first_gensys2_regime + get_setting(m, :temporary_altpolicy_length) : n_regimes
        if get_setting(m, :gensys2)
            gensys_regimes = [1:first_gensys2_regime-1]
            if last_gensys2_regime != n_regimes
                append!(gensys_regimes, [last_gensys2_regime+1:n_regimes])
            end
        else
            gensys_regimes = [1:n_regimes]
        end
        gensys2_regimes = [first_gensys2_regime-1:last_gensys2_regime]

        TTTs_vec[i], RRRs_vec[i], CCCs_vec[i] = solve(m; regime_switching = regime_switching,
                                                      regimes = collect(1:n_regimes),
                                                      gensys_regimes = gensys_regimes,
                                                      gensys2_regimes = gensys2_regimes, verbose = verbose)
        transitions[i] = Vector{Transition{T}}(undef, n_regimes)
        for j in 1:n_regimes # Compute vector of Transition for constructing the TimeVaryingInformationSetSystem
            transitions[i][j] = Transition(TTTs_vec[i][j], RRRs_vec[i][j], CCCs_vec[i][j])
        end
    end

    # Infer which measurement and pseudo-measurement equations to use
    measurement_eqns = Vector{Measurement{T}}(undef,       n_regimes)
    has_pseudo       = hasmethod(pseudo_measurement, (typeof(m), Matrix{T}, Matrix{T}, Vector{T}))
    if has_pseudo # Only calculate PseudoMeasurement equation if the method exists
        pseudo_measurement_eqns = Vector{PseudoMeasurement{T}}(undef, n_regimes)
    end

    for (reg, i) in enumerate(tvis_select)
        measurement_eqns[reg] = measurement(m, TTTs_vec[i][reg], RRRs_vec[i][reg], CCCs_vec[i][reg],
                                            reg = reg, TTTs = TTTs_vec[i], CCCs = CCCs_vec[i],
                                            information_set = tvis_infoset[reg])

        if has_pseudo
            pseudo_measurement_eqns[reg] = pseudo_measurement(m, TTTs_vec[i][reg], RRRs_vec[i][reg], CCCs_vec[i][reg], reg = reg,
                                                              TTTs = TTTs_vec[i], CCCs = CCCs_vec[i],
                                                              information_set = tvis_infoset[reg])
        end
    end

    if has_pseudo
        return TimeVaryingInformationSetSystem(transitions, measurement_eqns, pseudo_measurement_eqns,
                                               tvis_infoset, tvis_select)
    else
        return TimeVaryingInformationSetSystem(transitions, measurement_eqns, tvis_infoset, tvis_select)
    end
end


text:
# =============================================================
#		MODULE: table
# =============================================================
module table

	using DelimitedFiles

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	#		FUNCTION : Output(path.Output, ∑Psd, Id, N_∑Psd, Clay, Sand, Clay, SoilClass) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	function Output(Path, ∑Psd, Id, N_∑Psd, Sand, Silt, Clay, SoilClass)

		Header =  ["Id" "Sand" "Silt" "Clay" "SoilClass"]

		if isfile(Path) # If the file is there than delete it before saving figure
			  rm(Path)
		end

		DelimitedFiles.writedlm(Path, [Header; Id Sand Silt Clay SoilClass], ",")

	end

end  # module table
# ............................................................

max_stars_repo_name: joseph-pollacco/TextureTriangleNZ

max_issues_repo_name: joseph-pollacco/TextureTriangleNZ

max_forks_repo_name: joseph-pollacco/TextureTriangleNZ


text:
#!/usr/bin/julia

using Random

function passwd(name :: String, site :: String, seed :: Int64, n :: Int64)
    Random.seed!(seed)
    nums = map((x)->Int64(x), collect(name))
    seed = atan(Int64(sum(nums)) * 131);
    for i=1:n
        res = 213*rand()*sin(seed * (Int64(name[(i % length(name)) + 1]) / 13 + Int64(site[(i % length(site)) + 1]) / 31))
        seed *= 3*sin(res/5)+5
        nextch = Int64(round(abs(res))) % 26 + 'a'
	print(rand() < 0.5 ? nextch : nextch - 32)
    end

end

passwd(ARGS[1], ARGS[2], parse(Int64, ARGS[3]), parse(Int64, ARGS[4]))
println()


text:
###
# This script extracts roughly the last 2,000 songs from your "Liked" songs
## It includes all the default track attributes and dates added_at
## YMMV, I'm not sure if I am messing up the indexing, I was never good at counting past 20 anyways, you might see some duplicates which get
## filtered via the `unique` command on `id`
#
# Hope this helps get you started on using `Spotify.jl`!
###

## Libraries
using Spotify
using Spotify.Library
using Spotify.Tracks
using DataFrames

## Get tracks loop

function define_df()
    @info Defining the dataframe
    temp = library_get_saved_tracks(1)[1]["items"][1]
    global tracks_df = DataFrame(;Dict(temp["track"])...)
    tracks_df["added_at"] = temp["added_at"]
    delete!(tracks_df,1)
end

define_df()
##

@warn "Attempting to retrieve the last 2,000 songs from Spotify \n This may take some time"
for batch = 0:20:2000
    println("Getting batch: ", batch/20)
    temp = library_get_saved_tracks(20,(batch+1))[1]["items"]

    for i in 1:20
        temp2 = DataFrame(;Dict(temp[i]["track"])...)
        temp2["added_at"] = temp[i]["added_at"]
        println("Adding Track: ", temp[i]["track"]["name"])
        append!(tracks_df, temp2, cols=:union)

    end

    #sleep(rand(1:11))
end

##
size(tracks_df)

unique!(tracks_df, "id")

##
using CSV
something.(tracks_df, missing) |> CSV.write("track_data.csv")


text:
@testset "Verify evaluation of Probabilistic hermite functions" begin

    # Unscaled

    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = false), 0, x)

    for i=0:5
        @test norm(dV[:,i+1] - FamilyProPolyHermite[i+1].(x) .* exp.(-x.^2/4))<1e-8
    end

    # Scaled
    dV = vander(ProHermite(m; scaled = true), 0, x)

    for i=0:5
        @test norm(dV[:,i+1] - FamilyScaledProPolyHermite[i+1].(x).*exp.(-x.^2/4))<1e-8
    end
end

@testset "Verify 1st derivative of Probabilistic hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = false), 1, x)

    for i=0:5
        F = ProHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(F, xi), x) )<1e-8

        @test norm(dV[:,i+1] - (FamilyProPolyHermite[i+1].(x) .* (-0.5*x) .* exp.(-x.^2/4) +
            derivative(FamilyProPolyHermite[i+1],1).(x) .* exp.(-x.^2/4)))<1e-8
    end

    # Scaled
    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = true), 1, x)

    for i=0:5
        F = ProHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(F, xi), x) )<1e-8

        @test norm(dV[:,i+1] - (1/Cpro(i))*(FamilyProPolyHermite[i+1].(x) .* (-0.5*x) .* exp.(-x.^2/4) +
            derivative(FamilyProPolyHermite[i+1],1).(x) .* exp.(-x.^2/4)))<1e-8
    end
end

@testset "Verify second derivative of Probabilistic hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = false), 2, x)

    for i=0:5
        Hen = FamilyProPolyHermite[i+1].(x)
        Henp = derivative(FamilyProPolyHermite[i+1],1).(x)
        Henpp = derivative(FamilyProPolyHermite[i+1],2).(x)
        E = exp.(-x.^2/4)

        F = ProHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), x))<1e-8

        @test norm(dV[:,i+1] - (Henpp .*E - x .* Henp.*E + 0.5*Hen .* (0.5*x.^2 .- 1.0) .* E))<1e-8
    end

    # Scaled
    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = true), 2, x)

    for i=0:5
        Hen = FamilyProPolyHermite[i+1].(x)
        Henp = derivative(FamilyProPolyHermite[i+1],1).(x)
        Henpp = derivative(FamilyProPolyHermite[i+1],2).(x)
        E = exp.(-x.^2/4)

        F = ProHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), x))<1e-8

        @test norm(dV[:,i+1] - (1/Cpro(i))*(Henpp .*E - x .* Henp.*E + 0.5*Hen .* (0.5*x.^2 .- 1.0) .* E))<1e-8
    end
end

@testset "Verify third derivative of Probabilistic hermite functions" begin
    # Unscaled
    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = false), 3, x)

    for i=0:5
        F = ProHermite(i; scaled = false)
        @test norm(dV[:,i+1] - map(xi->Diff(F, 3)(xi),x))<1e-8
    end

    # Scaled

    x = randn(100)
    m = 5
    dV = vander(ProHermite(m; scaled = true), 3, x)

    for i=0:5
        F = ProHermite(i; scaled = true)
        @test norm(dV[:,i+1] - map(xi->Diff(F, 3)(xi),x))<1e-8
    end
end

# @testset "Verify integration of Probabilistic hermite functions" begin
#     Nx = 10
#     x = randn(10)
#
#     # Unscaled
#     dV = vander(ProHermite(8; scaled = false), -1, x)
#
#     for i=0:8
#         for j=1:Nx
#             @test abs(dV[j,i+1] - quadgk(y-> ProHermite(i)(y), 0.0, x[j], rtol=1e-12)[1])<1e-8
#         end
#     end
#
#     # Scaled
#     dV = vander(ProHermite(8; scaled = true), -1, x)
#
#     for i=0:8
#         for j=1:Nx
#             @test abs(dV[j,i+1] - quadgk(y-> ProHermite(i; scaled=true)(y), 0.0, x[j], rtol=1e-12)[1])<1e-8
#         end
#     end
#
# end


text:
include("parallel.jl")

# Load restart file and store solution in solver
load_restart_file!(dg, restart_filename) = load_restart_file!(dg, restart_filename, mpi_parallel())
function load_restart_file!(dg::AbstractDg, restart_filename, mpi_parallel::Val{false})
  # Create variables to be returned later
  time = NaN
  step = -1

  # Open file
  h5open(restart_filename, "r") do file
    # Read attributes to perform some sanity checks
    if read(attrs(file)["ndims"]) != ndims(dg)
      error("restart mismatch: ndims in solver differs from value in restart file")
    end
    if read(attrs(file)["equations"]) != get_name(equations(dg))
      error("restart mismatch: equations in solver differs from value in restart file")
    end
    if read(attrs(file)["polydeg"]) != polydeg(dg)
      error("restart mismatch: polynomial degree in solver differs from value in restart file")
    end
    if read(attrs(file)["n_elements"]) != dg.n_elements
      error("restart mismatch: polynomial degree in solver differs from value in restart file")
    end

    # Read time and time step
    time = read(attrs(file)["time"])
    step = read(attrs(file)["timestep"])

    # Read data
    varnames = varnames_cons(equations(dg))
    for v in 1:nvariables(dg)
      # Check if variable name matches
      var = file["variables_$v"]
      if (name = read(attrs(var)["name"])) != varnames[v]
        error("mismatch: variables_$v should be '$(varnames[v])', but found '$name'")
      end

      # Read variable
      dg.elements.u[v, .., :] = read(file["variables_$v"])
    end
  end

  return time, step
end

# Save current DG solution with some context information as a HDF5 file for
# restarting.
save_restart_file(dg, mesh, time, dt, timestep) = save_restart_file(dg, mesh, time, dt, timestep,
                                                                    mpi_parallel())
function save_restart_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep,
                           mpi_parallel::Val{false})
  # Create output directory (if it does not exist)
  output_directory = parameter("output_directory", "out")
  mkpath(output_directory)

  # Filename without extension based on current time step
  filename = joinpath(output_directory, @sprintf("restart_%06d", timestep))

  # Convert time and time step size to floats
  time = convert(Float64, time)
  dt = convert(Float64, dt)

  # Open file (clobber existing content)
  h5open(filename * ".h5", "w") do file
    # Add context information as attributes
    attrs(file)["ndims"] = ndims(dg)
    attrs(file)["equations"] = get_name(equations(dg))
    attrs(file)["polydeg"] = polydeg(dg)
    attrs(file)["n_vars"] = nvariables(dg)
    attrs(file)["n_elements"] = dg.n_elements_global
    attrs(file)["mesh_file"] = splitdir(mesh.current_filename)[2]
    attrs(file)["time"] = time
    attrs(file)["dt"] = dt
    attrs(file)["timestep"] = timestep

    # Restart files always store conservative variables
    data = dg.elements.u
    varnames = varnames_cons(equations(dg))

    # Store each variable of the solution
    for v in 1:nvariables(dg)
      # Convert to 1D array
      file["variables_$v"] = vec(data[v, .., :])

      # Add variable name as attribute
      var = file["variables_$v"]
      attrs(var)["name"] = varnames[v]
    end
  end
end

# Save current DG solution with some context information as a HDF5 file for
# postprocessing.
function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system="")
  return save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system,
                            mpi_parallel())
end
function save_solution_file(dg::AbstractDg, mesh::TreeMesh, time, dt, timestep, system,
                            mpi_parallel::Val{false})
  # Create output directory (if it does not exist)
  output_directory = parameter("output_directory", "out")
  mkpath(output_directory)

  # Filename without extension based on current time step
  if isempty(system)
    filename = joinpath(output_directory, @sprintf("solution_%06d", timestep))
  else
    filename = joinpath(output_directory, @sprintf("solution_%s_%06d", system, timestep))
  end

  # Convert time and time step size to floats
  time = convert(Float64, time)
  dt = convert(Float64, dt)

  # Open file (clobber existing content)
  h5open(filename * ".h5", "w") do file
    # Add context information as attributes
    attrs(file)["ndims"] = ndims(dg)
    attrs(file)["equations"] = get_name(equations(dg))
    attrs(file)["polydeg"] = polydeg(dg)
    attrs(file)["n_vars"] = nvariables(dg)
    attrs(file)["n_elements"] = dg.n_elements
    attrs(file)["mesh_file"] = splitdir(mesh.current_filename)[2]
    attrs(file)["time"] = time
    attrs(file)["dt"] = dt
    attrs(file)["timestep"] = timestep

    # Convert to primitive variables if requested
    solution_variables = parameter("solution_variables", "primitive",
                                   valid=["conservative", "primitive"])
    if solution_variables == "conservative"
      data = dg.elements.u
      varnames = varnames_cons(equations(dg))
    else
      # Reinterpret the solution array as an array of conservative variables,
      # compute the primitive variables via broadcasting, and reinterpret the
      # result as a plain array of floating point numbers
      data = Array(reinterpret(eltype(dg.elements.u),
             cons2prim.(reinterpret(SVector{nvariables(dg),eltype(dg.elements.u)}, dg.elements.u),
                        Ref(equations(dg)))))
      varnames = varnames_prim(equations(dg))
    end

    # Store each variable of the solution
    for v in 1:nvariables(dg)
      # Convert to 1D array
      file["variables_$v"] = vec(data[v, .., :])

      # Add variable name as attribute
      var = file["variables_$v"]
      attrs(var)["name"] = varnames[v]
    end

    # Store element variables
    for (v, (key, element_variables)) in enumerate(dg.element_variables)
      # Add to file
      file["element_variables_$v"] = element_variables

      # Add variable name as attribute
      var = file["element_variables_$v"]
      attrs(var)["name"] = string(key)
    end
  end
end

# Save current mesh with some context information as an HDF5 file.
save_mesh_file(mesh, timestep=-1) = save_mesh_file(mesh, timestep, mpi_parallel())
function save_mesh_file(mesh::TreeMesh, timestep, mpi_parallel::Val{false})
  # Create output directory (if it does not exist)
  output_directory = parameter("output_directory", "out")
  mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep >= 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d", timestep))
  else
    filename = joinpath(output_directory, "mesh")
  end

  # Create output directory (if it does not exist)
  # Open file (clobber existing content)
  h5open(filename * ".h5", "w") do file
    # Add context information as attributes
    n_cells = length(mesh.tree)
    attrs(file)["ndims"] = ndims(mesh)
    attrs(file)["n_cells"] = n_cells
    attrs(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
    attrs(file)["minimum_level"] = minimum_level(mesh.tree)
    attrs(file)["maximum_level"] = maximum_level(mesh.tree)
    attrs(file)["center_level_0"] = mesh.tree.center_level_0
    attrs(file)["length_level_0"] = mesh.tree.length_level_0
    attrs(file)["periodicity"] = collect(mesh.tree.periodicity)

    # Add tree data
    file["parent_ids"] = @view mesh.tree.parent_ids[1:n_cells]
    file["child_ids"] = @view mesh.tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view mesh.tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view mesh.tree.levels[1:n_cells]
    file["coordinates"] = @view mesh.tree.coordinates[:, 1:n_cells]
  end

  return filename * ".h5"
end


text:
# mock implementation for Julia versions where the Debugger doesn't work
module Debugger
isdebugging() = false
end


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

immutable JuMPArray{T,N,NT<:NTuple} <: JuMPContainer{T}
    innerArray::Array{T,N}
    indexsets::NT
    lookup::NTuple{N,Dict}
    meta::Dict{Symbol,Any}
end

@generated function JuMPArray{T,N}(innerArray::Array{T,N}, indexsets::NTuple{N})
    dicttuple = Expr(:tuple)
    for i in 1:N
        inner = quote
            idxset = indexsets[$i]
            ret = Dict{eltype(idxset), Int}()
        end
        tupelem = indexsets.parameters[i]
        if !(tupelem == UnitRange{Int} || tupelem == StepRange{Int})
            inner = quote
                $inner
                cnt = 1
                for x in idxset
                    ret[x] = cnt
                    cnt += 1
                end
                ret
            end
        end
        push!(dicttuple.args, inner)
    end
    :(JuMPArray(innerArray, indexsets, $dicttuple, Dict{Symbol,Any}()))
end

Base.getindex(d::JuMPArray, ::Colon) = d.innerArray[:]
@generated function Base.getindex{T,N,NT<:NTuple}(d::JuMPArray{T,N,NT}, idx...)
    if N != length(idx)
        error("Indexed into a JuMPArray with $(length(idx)) indices (expected $N indices)")
    end
    Expr(:call, :getindex, :(d.innerArray), _to_cartesian(d,NT,idx)...)
end

@generated function Base.setindex!{T,N,NT<:NTuple}(d::JuMPArray{T,N,NT}, v::T, idx...)
    if N != length(idx)
        error("Indexed into a JuMPArray with $(length(idx)) indices (expected $N indices)")
    end
    Expr(:call, :setindex!, :(d.innerArray), :v, _to_cartesian(d,NT,idx)...)
end

function _to_cartesian(d,NT,idx...)
    indexing = Any[]
    for (i,S) in enumerate(NT.parameters)
        idxtype = idx[1][i]
        if S == UnitRange{Int}
            if idxtype == Colon
                # special stuff
                push!(indexing, Colon())
            elseif idxtype <: Range
                push!(indexing, quote
                    rng = d.indexsets[$i]
                    I = idx[$i]
                    I - (start(rng) - 1)
                end)
            else
                push!(indexing, quote
                    rng = d.indexsets[$i]
                    I = idx[$i]
                    first(rng) <= I <= last(rng) || error("Failed attempt to index JuMPArray along dimension $($i): $I ∉ $(d.indexsets[$i])")
                    I - (start(rng) - 1)
                end)
            end
        elseif S == StepRange{Int}
            if idx[1][i] == Colon
                push!(indexing, Colon())
            else
                push!(indexing, quote
                    rng = $(d.indexsets[i])
                    I = idx[$i]
                    first(rng) <= I <= last(rng) || error("Failed attempt to index JuMPArray along dimension $($i): $I ∉ $(d.indexsets[$i])")
                    dv, rv = divrem(I - start(rng), step(rng))
                    rv == 0 || error("Failed attempt to index JuMPArray along dimension $($i): $I ∉ $(d.indexsets[$i])")
                    dv + 1
                end)
            end
        else
            push!(indexing, quote
                if !haskey(d.lookup[$i],idx[$i])
                    error("Failed attempt to index JuMPArray along dimension $($i): $(idx[$i]) ∉ $(d.indexsets[$i])")
                end
                d.lookup[$i][idx[$i]]::Int
            end)
        end
    end
    indexing
end