text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(sigmoid(-0.7172866311917714 + 0.8530094465195428*sigmoid(0.49975691401319544 + 0.817850285997455*sigmoid(-0.9044530867886094 + -0.3504940249833268*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*sigmoid(-0.49866226232151334 + 0.31623027699891626*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + -0.2302663292873568*sigmoid(0.6056562038214484 + 0.017530899400252764*sigmoid(-0.9044530867886094 + -0.3504940249833268*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*sigmoid(-0.49866226232151334 + 0.31623027699891626*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + sigmoid(0.11802092343800519 + -0.23198020992542823*sigmoid(0.49975691401319544 + 0.817850285997455*sigmoid(-0.9044530867886094 + -0.3504940249833268*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*sigmoid(-0.49866226232151334 + 0.31623027699891626*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + 0.3694431087223533*sigmoid(0.6056562038214484 + 0.017530899400252764*sigmoid(-0.9044530867886094 + -0.3504940249833268*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*sigmoid(-0.49866226232151334 + 0.31623027699891626*sigmoid(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*sigmoid(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
dims(H::Hermitian{<:Number,<:Tensor}) = dims(parent(H))

blockdims(H::Hermitian{<:Number,<:Tensor}, b) = blockdims(parent(H), b)

dim(H::Hermitian{<:Number,<:Tensor}, i::Int) = dim(parent(H), i)

matrix(H::Hermitian{<:Number,<:Tensor}) = Hermitian(matrix(parent(H)))

inds(H::Hermitian{<:Number,<:Tensor}) = inds(parent(H))

ind(H::Hermitian{<:Number,<:Tensor}, i::Int) = ind(parent(H), i)

nnzblocks(H::Hermitian{<:Number,<:Tensor}) = nnzblocks(parent(H))

nzblocks(H::Hermitian{<:Number,<:Tensor}) = nzblocks(parent(H))

eachnzblock(H::Hermitian{<:Number,<:Tensor}) = eachnzblock(parent(H))

eachblock(H::Hermitian{<:Number,<:Tensor}) = eachblock(parent(H))

eachdiagblock(H::Hermitian{<:Number,<:Tensor}) = eachdiagblock(parent(H))

nblocks(H::Hermitian{<:Number,<:Tensor}) = nblocks(parent(H))

function blockview(H::Hermitian{<:Number,<:Tensor}, block)
  return _blockview(H, blockview(parent(H), block))
end
_blockview(::Hermitian{<:Number,<:Tensor}, blockviewH) = Hermitian(blockviewH)
_blockview(::Hermitian{<:Number,<:Tensor}, ::Nothing) = nothing

ext: jl


text:
using Plots, ProgressMeter, LinearAlgebra
gr()

m, d  = 1000, 2
coef3 = 2 * rand(d,d,d) .- 1
coef2 = 2 * rand(d,d) .- 1
coef1 = 2 * rand(d,1)
coef0 = rand(1)

struct Stump

    j :: Int64
    theta :: Float64
    b :: Float64

    function Stump(j, theta, b)
        @assert j > 0
        new(j, theta, b)
    end

end

function (self :: Stump)( x )
    self.b * sign.(self.theta .- x[:,self.j])
end

function erreur(self :: Stump, x, y, D)
    sum( D .* (y .!= self(x)))
end

mutable struct Hypotheses

    weights :: Vector{Float64}
    stumps  :: Array{Stump,1}
    D       :: Vector{Vector{Float64}}

    function Hypotheses( )
        weights = Float64[]
        stumps  = Stump[]
        D       = Vector{Float64}[]
        new( weights, stumps, D)
    end

end

function (self::Hypotheses)(x)
    T = length(self.weights)
    m = size(x)[1]
    s = zeros(Float64, m )
    for t in 1:T
        s .+= self.weights[t] .* self.stumps[t](x)
    end
    sign.(s)
end

function erreur(self,x,y)
    m = length(y)
    sum(self(x) .!= y)/m
end

function plot_results( x, y, ytest)

    p    = plot(layout=(1,2); )
    blue = y .< 0
    red  = y .> 0
    scatter!(p[1,1],x[red,1],x[red,2],label=:Xb, title=:test)
    scatter!(p[1,1],x[blue,1],x[blue,2],label=:or)
    blue = ytest .< 0
    red  = ytest .> 0
    scatter!(p[1,2],x[red,1],x[red,2],label=:Xb, title=:test)
    scatter!(p[1,2],x[blue,1],x[blue,2],label=:or)
    p

end

function ERM_Stump( x :: Array{Float64,2},
                    y :: Vector{Float64},
                    D :: Vector{Float64})
    m, d       = size(x)
    Fstar      = 1e15
    F          = Fstar
    bstar      = 1
    thetastar  = 0.
    jstar      = 1
    indice_tri = hcat([sortperm(x[:,j]) for j in 1:d]...) #np.argsort(x,axis=0)
    xsort      = zeros(Float64,m+1)
    for j in 1:d
        xsort[1:m] = x[indice_tri,j][:,j]
        xsort[end] = xsort[m]+1
        ysort = y[indice_tri[:,j]]
        Dsort = D[indice_tri[:,j]]
        F = sum(Dsort .* (ysort .+ 1)/2)
        if F < Fstar
            bstar = 1
            Fstar = F
            thetastar = xsort[1] - 1
            jstar = j
        end
        for i in 1:m
            F = F - ysort[i] * Dsort[i]
            if F < Fstar && xsort[i] < xsort[i+1]
                bstar=1
                Fstar=F
                thetastar=(xsort[i]+xsort[i+1])/2
                jstar=j
            end
        end
        F = sum(Dsort .* (1 .- ysort)/2)
        if F < Fstar
            bstar = -1
            Fstar = F
            thetastar = xsort[1]-1
            jstar = j
        end
        for i in 1:m
            F = F + ysort[i] * Dsort[i]
            if F < Fstar && xsort[i]<xsort[i+1]
                bstar = -1
                Fstar = F
                thetastar = (xsort[i]+xsort[i+1])/2
                jstar = j
            end
        end
    end
    Stump(jstar, thetastar, bstar)
end

function generate_data(m, d)
    x = 2 * rand(m,d) .- 1
    y = zeros(Float64,m)
    for i in 1:d, j in 1:d, k in 1:d
        y .+= coef3[i,j,k] .* x[:,i] .* x[:,j] .* x[:,k]
    end
    y = sign.(y .+ diag((x * coef2) * transpose(x)) .+ x * coef1 .+ coef0);
    x, vec(y)
end

x, y = generate_data(m, d)

#d = 2
#m = 1000
#f = open("test1.txt")
#data = [ transpose(map(x->parse(Float64,x),split(line," "))) for line in eachline(f)]
#data = vcat(data...)
#
#x = data[:,1:2]
#y = data[:,3]

D = 1 / m .* ones(Float64, m)
h = ERM_Stump(x,y,D)
plot_results( x, y, h(x))

function adaboost( x, y , T)

    m = size(x)[1]
    D = 1 / m * ones(Float64,m)
    H = Hypotheses()
    epsilon = Float64[]

    @showprogress 1 for t in 1:T
        ht = ERM_Stump(x, y, D)
        hx = ht(x)
        e  = erreur(ht, x, y, D)
        w0 = 0.5 * log(1 / e-1)
        push!(epsilon, e)
        push!(H.weights, w0)
        push!(H.stumps, ht)
        #push!(H.D, D)
        D = D .* exp.(-w0 * y .* hx)
        D = D / sum(D)
    end
    H, epsilon
end

H, epsilon = adaboost( x, y, 100)
yresult = H(x)
plot_results( x, y, yresult)

a=100*(1-erreur(H,x,y))
println("Performance sur les données de test :",a[1],"%")

xtest, ytest = generate_data(m, d)
plot_results( xtest, ytest, H(xtest))

a=100*(1-erreur(H,xtest,ytest))
println("Performance sur les données de test :",a[1],"%")

err       = Float64[]
err0      = zeros(Float64,m)
err_test  = Float64[]
err_test0 = zeros(Float64,m)
for (weight, stump)  in zip(H.weights, H.stumps)
    err0 .+= weight .* stump(x)
    push!(err, sum( sign.(err0) .!= y) / m)
    err_test0 .+= weight .* stump(xtest)
    push!(err_test, sum(sign.(err_test0) .!= ytest) / m)
end

plot(err; marker=:o)
plot!(err_test; marker=:s)
xlabel!("Itérations")
ylabel!("Erreur en %")

scatter(epsilon)
xlabel!("steps")
ylabel!("Epsilon")
ylims!(minimum(epsilon), 0.5)

f = open("data_banknote_authentication.txt")
databrute = [ transpose(map(x->parse(Float64,x),split(line,","))) for line in eachline(f)]
data = vcat(databrute...);

@show n = size(data)[1]
print("Pourcentage de vrais billets :",(n-sum(data[:,5]))/n*100,"%")

using Random
taille       = n
taille_train = trunc(Int64,taille*0.8)
taille_test  = taille-taille_train
indice       = randperm(taille)
train        = data[indice[1:taille_train],:]
test         = data[indice[taille_train+1:end],:]
xtrain       = train[:,1:4]
ytrain       = 2 * train[:,5] .- 1
xtest        = test[:,1:4]
ytest        = 2 * test[:,5] .- 1
d            = 4
m            = taille_train

H, epsilon = adaboost(xtrain, ytrain, 50)

yresult     =  H(xtrain)
yresulttest =  H(xtest);

a=100*(1 .- erreur(H,xtrain,ytrain))
print("Performance sur les données d'entraînement :",a[1],"%")

a=100*(1-erreur(H,xtest,ytest))
print("Performance sur les données de test :",a[1],"%")

nt        = length(H.stumps)
err       = Float64[]
err0      = zeros(Float64,length(ytrain))
err_test  = Float64[]
err_test0 = zeros(Float64,length(ytest))

for i in 1:nt
    err0 .+= H.weights[i] .* H.stumps[i](xtrain)
    push!(err, sum(sign.(err0) .!= ytrain) / length(ytrain))
    err_test0 .+= H.weights[i] .* H.stumps[i](xtest)
    push!(err_test, sum(sign.(err_test0) .!= ytest)/length(ytest))
end

plot(err)
plot!(err_test)
xlabel!("Itérations ")
ylabel!("Erreur en %")

plot(epsilon)
xlabel!("Itérations")
ylabel!("Epsilon")
ylims!(minimum(epsilon), 0.5)

print("borne d'erreur théorique ",exp(-2*sum((1/2 .- epsilon).^2))*100,"%")

plot_results(xtrain, ytrain, yresult)

plot_results(xtest, ytest, yresulttest)

f = open("EEG_Eye_State.txt")

databrute = [ transpose(map(x->parse(Float64,x),split(line,","))) for line in eachline(f)]
data = vcat(databrute...)

n = size(databrute)[1]
println("Nombre de données : ", n)

count = sum( data[:,15] .== 0)
println("Pourcentage de classe 0 :",count/n*100,"%")

n_train  = trunc(Int64,n/3)
n_valid  = trunc(Int64,n/3)
n_test   = n-n_train-n_valid
indice   = randperm(n)
train    = data[indice[1:n_train],:]
valid    = data[indice[n_train+1:n_train+n_valid],:]
test     = data[indice[n_train+n_valid+1:end],:]
x        = train[:,1:14]
y        = 2 * train[:,15] .- 1
xvalid   = valid[:,1:14]
yvalid   = 2 * valid[:,15] .- 1
xtest    = test[:,1:14]
ytest    = 2 * test[:,15] .- 1
d        = 14
m        = n_train

(H,epsilon)=adaboost(x,y,2000)

yresult     = H(x)
yresulttest = H(xtest)

a = sum(yresult .== y) / m * 100
println("Performance sur les données d'entraînement :",a[1],"%")

a = sum(yresulttest .== ytest) / n_test*100
println("Performance sur les données de test :",a[1],"%")

err        = Float64[]
err0       = zeros(Float64,length(y))
err_valid  = Float64[]
err_valid0 = zeros(Float64,length(yvalid))
err_test   = []
err_test0  = zeros(Float64,length(ytest))
nt = length(H.stumps)
for i in 1:nt
    err0 .+= H.weights[i] .* H.stumps[i](x)
    push!(err, sum(sign.(err0) .!= y) / length(y))
    err_test0 .+= H.weights[i] .* H.stumps[i](xtest)
    push!(err_test, sum(sign.(err_test0) .!= ytest) / length(ytest))
    err_valid0 .+= H.weights[i]*H.stumps[i](xvalid)
    push!(err_valid,sum(sign.(err_valid0) .!= yvalid) / length(yvalid))
end

plot(err)
plot!(err_test)
plot!(err_valid)
xlabel!("Itérations")
ylabel!("Erreur en %")

plot(epsilon)
xlabel!("Itérations")
ylabel!("Epsilon")
ylims!(minimum(epsilon), 0.5)

println("borne d'erreur théorique : ",exp(-2*sum((1/2 .- epsilon).^2))*100,"%")

i = argmin(err_valid)
println(i)
println(err_test[i])
println(minimum(err_valid))


text:
using InverseIsing
using StatsBase:sample
using JSON

function estim_test1()
    samples = [1 -1 -1;-1 1 1]
    model = GBM(size(samples, 2))
    fit(model, samples)
    W_fit = infer(model)
    @test W_fit == [0 -1 -1;-1 0 1;-1 1 0]
    @test decode(W_fit) == Dict((1, 2) => -1, (1, 3) => -1, (2, 3) => 1)
end

function estim_test2()
    n_unit = 5 # Number of units
    rate = 0.9 # zero rate (in all interaction-pairs)
    h = Dict()
    J_idx = [(i,j) for i in 1:n_unit for j in i+1:n_unit]
    len = length(J_idx)
    J_val = Int.(ones(len))
    J_val[sample(1:len, Int(ceil(len * rate)), replace=false)] .= 0
    J = Dict(zip(J_idx, J_val))

    # annealing and create samples
    resp = anneal(h, J, n_read=1000)
    samples = hcat(resp.states...)'

    # inverse ising inference
    model = GBM(n_unit)
    fit(model, samples)
    pred = decode(infer(model))
    @test J == pred
end

@testset "Test estim" begin
    estim_test1()
    estim_test2()
end

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
@testset "Genetic Programming" begin
    Random.seed!(9874984737486);
    pop = 10
    terms = Terminal[:x, :y, rand]
    funcs = Function[+,-,*,/]

    t = TreeGP(pop, terms, funcs, maxdepth=2)
    @test Evolutionary.population_size(t) == pop
    @test sort(Evolutionary.terminals(t)) == [:x, :y]
    @testset for (term, dim) in t.terminals
        @test dim == 1
    end
    @testset for (func, arity) in t.functions
        @test arity == 2
    end
    @test_skip summary(t) == "TreeGP[P=10,Parameter[x,y],Function[*, +, /, -]]"

    # population initialization
    popexp = Evolutionary.initial_population(t);
    @test length(popexp) == pop
    popexp = Evolutionary.initial_population(t, :(x + 1));
    @test popexp[1] == :(x + 1)

    # recursive helper functions
    Random.seed!(9874984737482)
    gt = rand(TreeGP(pop, terms, funcs, maxdepth=2, initialization=:grow), 3)
    @test Evolutionary.nodes(gt) < 15
    @test Evolutionary.height(gt) <= 3
    @test length(gt) < 15
    ft = rand(TreeGP(pop, terms, funcs, maxdepth=2, initialization=:full), 3)
    @test Evolutionary.nodes(ft) == 15
    @test Evolutionary.height(ft) == 3
    @test length(ft) == 15
    # @test Evolutionary.depth(ft, :x) == 3
    # ft[3] = :z
    # @test Evolutionary.depth(ft, :z) == 3
    @test Evolutionary.depth(ft, ft) == 0
    @test Evolutionary.depth(ft, ft[3]) > 0
    @test Evolutionary.depth(ft, :w) == -1
    @test Evolutionary.evaluate([1.0, 2.0], :y, [:y, :z]) == 1.0
    copyto!(ft, gt)
    @test ft == gt
    # @test Evolutionary.symbols(ft) |> sort == [:x, :y]

    # simplification
    @test Expr(:call, -, :x, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, :x, :x) |> Evolutionary.simplify! == 1
    @test Expr(:call, *, 0, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, *, :x, 0) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, 0, :x) |> Evolutionary.simplify! == 0
    @test Expr(:call, /, 1, 0) |> Evolutionary.simplify! == 1
    @test Expr(:call, +, 0, :x) |> Evolutionary.simplify! == :x
    @test Expr(:call, +, :x, 0) |> Evolutionary.simplify! == :x
    @test Expr(:call, -, :x, 0) |> Evolutionary.simplify! == :x
    @test Expr(:call, +, :x, :x) |> Evolutionary.simplify! == Expr(:call, *, 2, :x)
    @test Expr(:call, -, 2, 1) |> Evolutionary.simplify! == 1

    # evaluation
    ex = Expr(:call, +, 1, :x) |> Evolutionary.Expression
    xs = rand(10)
    @test ex(xs[1]) == xs[1]+1
    @test ex.(xs) == xs.+1
    io = IOBuffer()
    show(io, "text/latex", ex)
    @test String(take!(io)) == "\\left(1.0+x\\right)"

    depth = 5
    fitfun(x) = x*x + x + 1.0
    function fitobj(expr)
        rg = -5.0:0.5:5.0
        ex = Evolutionary.Expression(expr)
        sum(v->isnan(v) ? 1.0 : v, abs2.(fitfun.(rg) - ex.(rg)) )/length(rg) |> sqrt
    end

    Random.seed!(9874984737426);
    res = Evolutionary.optimize(fitobj,
        TreeGP(50, Terminal[:x, randn], Function[+,-,*,Evolutionary.pdiv],
            mindepth=1,
            maxdepth=depth,
            optimizer = GA(
                selection = tournament(5),
                ɛ = 0.1,
                mutationRate = 0.8,
                crossoverRate = 0.2,
            ),
        )
    )
    @test minimum(res) < 1.1

end


text:
module Tabular

using ..FastAI
using ..FastAI:
    # blocks
    Block, WrapperBlock, AbstractBlock, OneHotTensor, OneHotTensorMulti, Label,
    LabelMulti, wrapped, Continuous, getencodings, getblocks, encodetarget, encodeinput,
    # encodings
    Encoding, StatefulEncoding, OneHot,
    # visualization
    ShowText,
    # other
    Context, Training, Validation
import ..FastAI: Datasets
using ..FastAI.Datasets

# for tests
using ..FastAI: testencoding

# extending
import ..FastAI:
    Datasets,
    blockmodel, blockbackbone, blocklossfn, encode, decode, checkblock,
    encodedblock, decodedblock, showblock!, mockblock, setup
using ..FastAI.Datasets

import DataAugmentation
import DataFrames: DataFrame
import MLUtils: MLUtils, eachobs, getobs, numobs
import Flux
import Flux: Embedding, Chain, Dropout, Dense, Parallel, BatchNorm
import PrettyTables
import Requires: @require
import ShowCases: ShowCase
import Tables
import Statistics

using InlineTest

# Blocks
include("blocks/tablerow.jl")

# Encodings
include("encodings/tabularpreprocessing.jl")

include("models.jl")

const _tasks = Dict{String, Any}()
include("tasks/classification.jl")
include("tasks/regression.jl")
include("recipes.jl")

function __init__()
    _registerrecipes()
    foreach(values(_tasks)) do t
        if !haskey(learningtasks(), t.id)
            push!(learningtasks(), t)
        end
    end
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        import .Makie
        import .Makie: @recipe, @lift
        import .FastAI: ShowMakie
        include("makie.jl")
    end
end

export TableRow, TabularPreprocessing, TabularClassificationSingle, TabularRegression

end


text:
using HSARuntime
using Test

if HSARuntime.configured
    agents = get_agents()
    if length(agents) > 0
        agent = get_default_agent()
        @test agent !== nothing

        HSARuntime.set_default_agent!(:gpu)
        agent = get_default_agent()
        @test HSARuntime.device_type(agent) == :gpu

        agent_name = HSARuntime.get_name(agent)
        @test length(agent_name) > 0

        agent_isa = get_first_isa(agent)
        @test length(agent_isa) > 0

        @info "Testing with agent $agent_name with ISA $agent_isa"

        arr = HSAArray(agent, UInt32, (8, ))
        fill!(arr, UInt32(1))
        _arr = Array(arr)
        @test length(arr) == length(_arr)
        for idx in 1:length(arr)
            _arr[idx] = arr[idx]
        end
        @test arr == _arr
        fill!(_arr, UInt32(2))
        for idx in 1:length(arr)
            arr[idx] = _arr[idx]
        end
        @test arr == _arr

        arr1 = HSAArray(agent, UInt32, (8,))
        arr2 = HSAArray(agent, UInt32, (8,))
        fill!(arr1, 1)
        fill!(arr2, 2)
        arr3 = similar(arr1)
        for i in 1:length(arr1)
            arr3[i] = arr1[i] + arr2[i]
        end
        @test all(arr3 .== 3)
        @test all(isa.((arr1, arr2, arr3), Ref(HSAArray)))
        @test size(arr3) == size(arr1)

        @testset "API wrappers" begin
            include("memory.jl")
        end
    else
        @warn("No devices detected; skipping on-device tests")
    end
else
    @warn("HSARuntime.jl has not been configured; skipping on-device tests.")
end

#=
queue = HSAQueue(agent)
kernel = HSAKernel(agent, executable, "&__vector_copy_kernel")
inp = zeros(HSAArray{UInt32}, 1024, 1024)
out = ones(HSAArray{UInt32}, 1024, 1024)
args = HSAKernelArgs(inp, out)
signal = HSASignal()
launch(queue, kernel, args, signal;
       workgroup_size=(256,1,1),
       grid_size=(1024*1024,1,1))
wait(signal)
=#


text:
module DefineSingletons

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end DefineSingletons

export @def_singleton

"""
    @def_singleton singleton_name
    @def_singleton singleton_name::SingletonType
    @def_singleton singleton_name::SingletonType <: SuperType
    @def_singleton singleton_name isa SuperType
    @def_singleton singleton_name = SingletonType()

Define a singleton named `singleton_name` and its two-argument
`show(::IO, ::SingletonType)` method.

With the form `singleton_name = SingletonType()`, the type
`SingletonType` can be a parametric type as long as
`Base.issingletontype(SingletonType)` is `true`.

# Examples
```jldoctest
julia> using DefineSingletons

julia> @def_singleton mysingleton::MySingletonType;

julia> mysingleton
mysingleton

julia> mysingleton isa MySingletonType
true

julia> mysingleton === MySingletonType()
true
```

With supertype:

```jldoctest; setup = :(using DefineSingletons)
julia> abstract type MySuperType end;

julia> @def_singleton mysingleton2::MySingletonType2 <: MySuperType;

julia> mysingleton2
mysingleton2

julia> mysingleton2 isa MySingletonType2
true

julia> MySingletonType2 <: MySuperType
true

julia> @def_singleton mysingleton3 isa MySuperType;

julia> mysingleton3 isa MySuperType
true
```

With pre-existing parametric type:

```jldoctest; setup = :(using DefineSingletons)
julia> struct MyParametricType{T} end;

julia> @def_singleton P1 = MyParametricType{1}();

julia> P1
P1

julia> P1 isa MyParametricType{1}
true
```
"""
macro def_singleton(ex)
    ans = handle_gentype(ex)
    ans === nothing || return ans
    ans = handle_predefined(ex)
    ans === nothing || return ans
    throw(ArgumentError("invalid input: $ex"))
end

function handle_gentype(ex)
    super_type = Any
    type_name = nothing
    if ex isa Symbol
        singleton_name = ex
    elseif Meta.isexpr(ex, :(::), 2)
        singleton_name, type_name = ex.args
        if Meta.isexpr(type_name, :(<:), 2)
            type_name, super_type = type_name.args
        end
    elseif Meta.isexpr(ex, :(<:), 2)
        x, super_type = ex.args
        if Meta.isexpr(x, :(::), 2)
            singleton_name, type_name = x.args
        else
            return nothing
        end
    elseif Meta.isexpr(ex, :call, 3) && ex.args[1] === :isa
        _, singleton_name, super_type = ex.args
    else
        return nothing
    end
    if type_name === nothing
        type_name = gensym(string("typeof_", singleton_name))
    end
    show_def = define_show(singleton_name, type_name)
    singleton_name = esc(singleton_name)
    type_name = esc(type_name)
    super_type = super_type isa Type ? super_type : esc(super_type)
    quote
        struct $type_name <: $super_type end
        const $singleton_name = $type_name()
        $show_def
        $singleton_name
    end
end

function handle_predefined(ex)
    if Meta.isexpr(ex, :(=), 2)
        singleton_name, type_call = ex.args
        if Meta.isexpr(type_call, :call, 1)
            type_name, = type_call.args
        else
            return nothing
        end
    else
        return nothing
    end
    show_def = define_show(singleton_name, type_name)
    singleton_name = esc(singleton_name)
    type_name = esc(type_name)
    quote
        const $singleton_name = $type_name()
        $show_def
        $singleton_name
    end
end

function define_show(singleton_name::Symbol, type_name)
    singleton_name = string(singleton_name)
    type_name = esc(type_name)
    quote
        function Base.show(io::IO, x::$type_name)
            if !get(io, :limit, false)
                # Don't show full name in REPL etc.:
                print(io, parentmodule($type_name), '.')
            end
            print(io, $singleton_name)
        end
    end
end

end # module


text:
mutable struct SkyRoom2
    led_arduino::LEDArduino
    camera::PiCamera
    frame::Observable
    function SkyRoom2()
        led_arduino = LEDArduino(led_port["skyroom2"])
        camera = PiCamera(30, 67, 67, 4)
        frame = Observable(snap(camera))
        new(led_arduino, camera, frame)
    end
end

connect!(a::SkyRoom2) = @async while isopen(a.camera)
    a.frame[] = snap(a.camera)
    sleep(0.01)
end

function restart(a::SkyRoom2)
    empty!(a.frame.listeners)
    restart(a.led_arduino)
    restart(a.camera)
    connect!(a)
end

mutable struct SkyRoom1
    wind_arduinos::Vector{FanArduino}
    led_arduino::LEDArduino
    camera::PiCamera
    data::Observable
    function SkyRoom1()
        wind_arduinos = [FanArduino(id, port) for (id, port) in enumerate(fan_ports) if isconnected(port)]
        led_arduino = LEDArduino(led_port["skyroom"])
        camera = PiCamera(30, 67, 67, 4)
        data = Observable((; frame = snap(camera), trpms = get_rpms(wind_arduinos)))
        new(wind_arduinos, led_arduino, camera, data)
    end
end

connect!(a::SkyRoom1) = @async while all(isopen, a.wind_arduinos) && isopen(a.camera)
    a.data[] = (; frame = snap(a.camera), trpms = get_rpms(a.wind_arduinos))
    sleep(0.01)
end

function restart(a::SkyRoom1)
    empty!(a.data.listeners)
    restart.(a.wind_arduinos)
    restart(a.led_arduino)
    restart(a.camera)
    connect!(a)
end

function dropdown(options, option)
    dropdown_onchange = js"update_obs($option, this.options[this.selectedIndex].text);"
    DOM.select(DOM.option.(options); class="bandpass-dropdown", onclick=dropdown_onchange)
end

function recordfans(trpms, folder, ids)
    fan_io = open(folder / "fans.csv", "w")
    println(fan_io, "time,", join([join(["fan$(id)_speed$j" for j in 1:3], ",") for id in ids], ","))
    on(trpms) do (t, rpms)
        println(fan_io, t, ",",join(Iterators.flatten(rpms), ","))
    end
end

# label_setup(x) = string("fans=", Int[i.pwm for i in x.fans], "; stars=", [string(i.cardinality, " ", i.elevation, " ", i.intensity, " ", i.radius) for i in x.stars])

function update_arena!(wind_arduinos, led_arduino, setup)
    for a in wind_arduinos
        a.pwm[] = setup.fans[a.id].pwm
    end
    led_arduino.pwm[] = parse2arduino(setup.stars)
end

function backup(left2upload, bucket)
    left2upload[] = 1.0
    n = length(readpath(datadir))
    for (i, folder) in enumerate(readpath(datadir))
        tmp = Tar.create(string(folder))
        rm(folder, recursive = true)
        name = basename(folder)
        run(`aws s3 mv $tmp s3://$bucket/$name.tar --quiet`)
        left2upload[] = (n - i)/n
    end
end

typedict(x::T) where {T} = Dict(fn=>getfield(x, fn) for fn ∈ fieldnames(T))

function todict(setup)
    x = Dict(pairs(setup))
    if haskey(x, :fans)
        x[:fans] = typedict.(x[:fans])
    end
    x[:stars] = typedict.(x[:stars])
    return x
end

function dom_handler(sr::SkyRoom1, left2upload, session, request)

    restart(sr)

    frame = map(sr.data) do x
        x.frame
    end
    trpms = map(sr.data) do x
        x.trpms
    end

    rpms = map(trpms) do (_, x)
        [ismissing(x) ? 0.0 : x for x in Iterators.flatten(x)]
    end

    on(trpms) do (t, rpms)
        println(fan_io[], t, ",",join(Iterators.flatten(rpms), ","))
    end

    md = Dict()

    setup_file = download(setupsurl["skyroom"])
    df = CSV.File(setup_file, header = 1:2, types = Dict(1 => String))  |> TableOperations.transform(setup_label = strip) |> DataFrame

    setuplog = similar(df[1:1,:])
    setuplog[:, :time] .= now()
    empty!(setuplog)

    setups = select(df, :setup_label => identity => :label, r"fan" => ByRow(parse2wind ∘ tuple) => :fans, r"star" => ByRow(parse2stars ∘ tuple) => :stars)
    options = collect(df.setup_label)
    option = Node(first(options))
    setup = map(option) do o
        i = findfirst(==(o), options)
        push!(setuplog, merge((; time = now()), NamedTuple(df[i,:])))
        setups[i, :]
    end

    on(setup) do x
        for _ in 1:3
            update_arena!(sr.wind_arduinos, sr.led_arduino, x)
            sleep(0.05)
        end
    end

    # GC.gc(true)

    fanrecording = Observable(recordfans(trpms, tmpdir(), [a.id for a in sr.wind_arduinos]))
    off(trpms, fanrecording[])

    recordingtime = Node(now())
    timestamp = map(string, recordingtime)
    timestamp[] = ""
    recording = JSServe.Checkbox(false)

    on(recording) do tf
        if tf
            delete!(setuplog, 1:nrow(setuplog) - 1)
            recordingtime[] = now()
            folder = datadir / timestamp[]
            mkdir(folder)
            sr.camera.cam.start_recording(string(folder / "video.h264"))
            fanrecording[] = recordfans(trpms, folder, [a.id for a in sr.wind_arduinos])
        else
            off(trpms, fanrecording[])
            sr.camera.cam.stop_recording()
        end
    end
    comment = JSServe.TextField("")
    beetleid = JSServe.TextField("")
    save = JSServe.Button("Save")
    on(save) do _
        if recording[]
            recording[] = false
        end

        md["setup"] = todict(setup[])
        md["recording_time"] = recordingtime[]
        md["beetle_id"] = beetleid[]
        md["comment"] = comment[]
        folder = datadir / timestamp[]
        open(folder / "metadata.toml", "w") do io
            TOML.print(io, md)
        end
        CSV.write(folder / "setups_log.csv", setuplog)

        timestamp[] = ""
        comment[] = ""
        beetleid[] = ""
    end
    upload = JSServe.Button("Backup")
    on(_ -> backup(left2upload, "nicolas-cage-skyroom"), upload)
    left2upload_label = map(left2upload) do i
        if i > 0
            string(round(Int, 100i), "% left to upload...")
        else
            "all backed up"
        end
    end

    # GC.gc(true)

    colors = repeat(1:5, inner = [3])
    rpmx = vcat(((i - 1)*4 + 1 : 4i - 1 for i in 1:5)...)
    rpmplot = Scene(show_axis = false, resolution = (540, round(Int, 3*5 + 3*540/(3*5+4))))
    barplot!(rpmplot, rpmx, top_rpm*ones(3*5), color = :white, strokecolor = :black, strokewidth = 1)
    barplot!(rpmplot, rpmx, rpms, color = colors, strokecolor = :transparent, strokewidth = 0)
    disconnect!(AbstractPlotting.camera(rpmplot))

    frameplot = image(frame, scale_plot = false, show_axis = false)
    disconnect!(AbstractPlotting.camera(frameplot))

    # GC.gc(true)

    # print_sizes()

    empty!(WGLMakie.SAVE_POINTER_IDENTITY_FOR_TEXTURES)

    return DOM.div(
        DOM.div(rpmplot),
        DOM.div(frameplot),
        DOM.div("Setup: ", dropdown(options, option)),
        DOM.div("Recording: ", recording, timestamp),
        DOM.div("ID: ", beetleid),
        DOM.div("Comment: ", comment),
        DOM.div(save),
        DOM.div(upload, left2upload_label)
    )

end

function dom_handler(sr::SkyRoom2, left2upload, session, request)

    restart(sr)

    md = Dict()
    setup_file = download(setupsurl["skyroom2"])
    df = CSV.File(setup_file, header = 1:2, types = Dict(1 => String))  |> TableOperations.transform(setup_label = strip) |> DataFrame
    setups = select(df, :setup_label => identity => :label, r"star" => ByRow(parse2stars ∘ tuple) => :stars)
    options = collect(df.setup_label)
    option = Node(first(options))
    setup = map(option) do o
        i = findfirst(==(o), options)
        setups[i, :]
    end

    on(setup) do x
        for _ in 1:3
            sr.led_arduino.pwm[] = parse2arduino(x.stars)
            sleep(0.05)
        end
    end

    # GC.gc(true)

    recordingtime = Node(now())
    timestamp = map(string, recordingtime)
    timestamp[] = ""
    recording = JSServe.Checkbox(false)
    on(recording) do tf
        if tf
            recordingtime[] = now()
            folder = datadir / timestamp[]
            mkdir(folder)
            sr.camera.cam.start_recording(string(folder / "video.h264"))
        else
            sr.camera.cam.stop_recording()
        end
    end
    comment = JSServe.TextField("")
    beetleid = JSServe.TextField("")
    save = JSServe.Button("Save")
    on(save) do _
        if recording[]
            recording[] = false
        end

        md["setup"] = todict(setup[])
        md["recording_time"] = recordingtime[]
        md["beetle_id"] = beetleid[]
        md["comment"] = comment[]
        folder = datadir / timestamp[]
        open(folder / "metadata.toml", "w") do io
            TOML.print(io, md)
        end

        timestamp[] = ""
        comment[] = ""
        beetleid[] = ""
    end
    upload = JSServe.Button("Backup")
    on(_ -> backup(left2upload, "top-floor-skyroom2"), upload)
    left2upload_label = map(left2upload) do i
        if i > 0
            string(round(Int, 100i), "% left to upload...")
        else
            "all backed up"
        end
    end

    # GC.gc(true)

    frameplot = image(sr.frame, scale_plot = false, show_axis = false)
    disconnect!(AbstractPlotting.camera(frameplot))

    # GC.gc(true)

    # print_sizes()

    empty!(WGLMakie.SAVE_POINTER_IDENTITY_FOR_TEXTURES)

    return DOM.div(
        DOM.div(frameplot),
        DOM.div("Setup: ", dropdown(options, option)),
        DOM.div("Recording: ", recording, timestamp),
        DOM.div("ID: ", beetleid),
        DOM.div("Comment: ", comment),
        DOM.div(save),
        DOM.div(upload, left2upload_label)
    )
end

# function print_sizes()
#     sizes = sort(filter(x-> x[2] > 5*(10^6), map(collect(values(Base.loaded_modules))) do m
#                             m=>Base.summarysize(m)
#                         end), rev=true, by=last)
#     foreach(((m, s),)-> println(m, ": ", Base.format_bytes(s)), sizes)
#     println("Free: ", Base.format_bytes(Sys.free_memory()))
# end

# function print_sizes()
#     mem = Pair{String, Int}[]
#     for m in values(Base.loaded_modules), vs in names(m, all = true)
#         if isdefined(m, vs)
#             v = getfield(m, vs)
#             x = Base.summarysize(v)
#             if x > 10^6
#                 push!(mem, string(m, "; ", vs, ": ") => x)
#             end
#         end
#     end
#     sort!(mem, by = last, rev = true)
#     n = length(mem)
#     if n > 15
#         deleteat!(mem, 16:n)
#     end
#     println("Memory usage:")
#     for (txt, i) in mem
#         println(txt, Base.format_bytes(i))
#     end
#     println("")
# end


text:
# Default implementations of ONNX operators

import NNlib
import Flux

flipweights(w::AbstractArray{T, N}) where {T,N} = w[(size(w,i):-1:1 for i in 1:(N-2))..., :, :]

conv(x, w; kw...) = NNlib.conv(x, flipweights(w); kw...)

function conv(x, w, b; kw...)
    d = ndims(x) - 2
    bias_size = (ntuple(_ -> 1, d)..., :, 1)
    b = reshape(b, bias_size)
    return conv(x, w; kw...) .+ b
end

function onnx_gemm(A, B, C; tA=0, tB=false, α=1, β=1)
    A = Bool(tA) ? A' : A
    B = Bool(tB) ? B' : B
    # note: order of arguments reversed due to row-major layout
    return α * B * A .+ β * C
end

function onnx_gemm(A, B; tA=0, tB=0, α=1)
    A = Bool(tA) ? A' : A
    B = Bool(tB) ? B' : B
    # note: order of arguments reversed due to row-major layout
    return α * B * A
end

# Julia-friendly flatten
function flatten(x; dim=ndims(x) - 1)
    sz = size(x)
    keep = dim < ndims(x) ? sz[dim + 1:end] : 1
    return reshape(x, :, keep...)
end

# ONNX-specific flatten
function onnx_flatten(x; axis=1)
    dim = axis >= 0 ? ndims(x) - axis + 1 : axis + 1
    return flatten(x; dim=dim)

end

add(xs...) = +(xs...)
mul(xs...) = .*(xs...)
relu(x) = NNlib.relu.(x)
maxpool(x; kernel, pad=0, stride=1) = NNlib.maxpool(x, kernel; pad=pad, stride=stride)

# mutable struct BatchNorm{F,V,N,W}
#     λ::F  # activation function
#     β::V  # bias
#     γ::V  # scale
#     μ::W     # moving mean
#     σ²::W    # moving var
#     ϵ::N
#     momentum::N
#     affine::Bool
#     track_stats::Bool
#     active::Union{Bool, Nothing}
#     chs::Int # number of channels
#   end

function batch_norm(x, γ, β, μ, σ², ϵ, momentum, training_mode)
    bn = Flux.BatchNorm(
        identity, β, γ, μ, σ², ϵ, momentum, true, true, nothing, length(γ))
    y = bn(x)
    if training_mode
        return y, bn.μ, bn.σ²  # TODO: are bn.μ and bn.σ² actually updated?
    else
        return y
    end
end

function global_average_pool(x)
    return Flux.GlobalMeanPool()(x)
end


text:
# A delay eigenvalue problem
function dep0(n::Int=5)
    msws_rng = MSWS_RNG()
    A0 = gen_rng_mat(msws_rng,n,n)
    A1 = gen_rng_mat(msws_rng,n,n)
    tau = 1.0
    nep = DEP([A0,A1],[0,tau])
    return nep
end

# A delay eigenvalue problem with sparse matrices
function dep0_sparse(n::Integer=100, p::Real=0.25)
    msws_rng = MSWS_RNG()
    A0 = sparse(1:n,1:n,vec(gen_rng_mat(msws_rng,n,1)))+gen_rng_spmat(msws_rng,n,n,p)
    A1 = sparse(1:n,1:n,vec(gen_rng_mat(msws_rng,n,1)))+gen_rng_spmat(msws_rng,n,n,p)
    tau = 1.0
    nep = DEP([A0,A1],[0,tau])
    return nep
end

# A delay eigenvalue problem with sparse tridiagonal matrices
function dep0_tridiag(n::Integer=100)
    msws_rng = MSWS_RNG()
    K = [1:n;2:n;1:n-1]; J=[1:n;1:n-1;2:n] # sparsity pattern of tridiag matrix
    A0 = sparse(K, J, vec(gen_rng_mat(msws_rng,3*n-2,1)))
    A1 = sparse(K, J, vec(gen_rng_mat(msws_rng,3*n-2,1)))
    tau = 1.0
    nep = DEP([A0,A1],[0,tau])
    return nep
end

# A polynomial eigenvalue problem
function pep0(n::Integer=200)
    msws_rng = MSWS_RNG()
    A0=gen_rng_mat(msws_rng,n,n)
    A1=gen_rng_mat(msws_rng,n,n)
    A2=gen_rng_mat(msws_rng,n,n)
    A=[A0,A1,A2]
    nep=PEP(A)
    return nep
end

# A polynomial eigenvalue problem with symmetric matrices
function pep0_sym(n::Integer=200)
    msws_rng = MSWS_RNG()
    A0 = Symmetric(gen_rng_mat(msws_rng,n,n))
    A1 = Symmetric(gen_rng_mat(msws_rng,n,n))
    A2 = Symmetric(gen_rng_mat(msws_rng,n,n))
    A = [A0, A1, A2]
    nep = PEP(A)
    return nep
end

# A sparse polynomial eigenvalue problem
function pep0_sparse(n::Integer=200, p::Real=0.03)
    msws_rng = MSWS_RNG()
    A0=gen_rng_spmat(msws_rng,n,n,p)
    A1=gen_rng_spmat(msws_rng,n,n,p)
    A2=gen_rng_spmat(msws_rng,n,n,p)
    A=[A0,A1,A2]
    nep=PEP(A)
    return nep
end

# Helper for random matrices, stable over releases
# B. Widynski, Middle Square Weyl Sequence RNG, arXiv 1704.00358
mutable struct MSWS_RNG
    x::UInt128
    w::UInt128
    s::UInt128
    function MSWS_RNG(seed::UInt128 = zero(UInt128))
        base = 0x9ef09a97ac0f9ecaef01c4f2db0958c9
        s = (seed << 1) + base
        x = 0x1de568e1a1ca1b593cbf13f7407cf43e
        w = 0xd4ac5c288559e14a5fafc1b7df9f9e0e
        return new(x, w, s)
    end
end

function gen_rng_int(rng::MSWS_RNG)
    rng.x *= rng.x
    rng.x += (rng.w += rng.s)
    rng.x = (rng.x>>64) | (rng.x<<64)
    return UInt64(rng.x & typemax(UInt64))
end

function gen_rng_float(rng::MSWS_RNG)
    return Float64(gen_rng_int(rng)/typemax(UInt64))
end

function gen_rng_mat(rng::MSWS_RNG, n::Integer, m::Integer)
    A = zeros(Float64,n,m)
    for c = 1:m
        for r = 1:n
            A[r,c] = 1-2*gen_rng_float(rng)
        end
    end
    return A
end

function gen_rng_spmat(rng::MSWS_RNG, n::Integer, m::Integer, p::Real)
    n=Int64(n)
    m=Int64(m)
    nonzeros = round(p*m*n)
    dict = Dict{Tuple{Int64,Int64},Float64}()
    for i = 1:nonzeros
        r = mod(gen_rng_int(rng),n) + 1
        c = mod(gen_rng_int(rng),m) + 1
        dict[r,c] = 1-2*gen_rng_float(rng)
    end
    idxes = collect(keys(dict))
    nnzs = length(idxes)
    r = zeros(Int64,nnzs)
    c = zeros(Int64,nnzs)
    vals = zeros(Float64,nnzs)
    for i = 1:nnzs
        r[i] = idxes[i][1]
        c[i] = idxes[i][2]
        vals[i] = dict[idxes[i]]
    end
    return sparse(r,c,vals,n,m)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ e89c8374-8a79-11eb-3029-7d8170bde727
begin
	using PlutoUI
	using NativeSVG
	using GDBM
end

# ╔═╡ fa2d4204-8a79-11eb-0a0a-69a684f6c422
md"""# Files

This chapter introduces the idea of persistent programs that keep data in permanent storage, and shows how to use different kinds of permanent storage, like files and databases.
"""

# ╔═╡ 087bcda8-8a7a-11eb-0bdf-758ab26f7fae
md"""## Persistence

Most of the programs we have seen so far are transient, in the sense that they run for a short time and produce some output, but when they end, their data disappears. If you run the program again, it starts with a clean slate.

Other programs are *persistent*: they run for a long time (or all the time); they keep at least some of their data in permanent storage (a hard drive, for example); and if they shut down and restart, they pick up where they left off.

Examples of persistent programs are operating systems, which run pretty much whenever a computer is on, and web servers, which run all the time, waiting for requests to come in on the network.

One of the simplest ways for programs to maintain their data is by reading and writing text files. We have already seen programs that read text files; in this chapter we will see programs that write them.

An alternative is to store the state of the program in a database. In this chapter I will also present how to use a simple database.
"""

# ╔═╡ 23dd5184-8a7a-11eb-1dd0-35b05d248988
md"""## Reading and Writing

A *text file* is a sequence of characters stored on a permanent medium like a hard drive or flash memory. We saw how to open and read a file in “Reading Word Lists”.

To write a file, you have to open it with mode `"w"` as a second parameter:

```julia
julia> fout = open("output.txt", "w")
IOStream(<file output.txt>)
```

If the file already exists, opening it in write mode clears out the old data and starts fresh, so be careful! If the file doesn’t exist, a new one is created. open returns a file object and the write function puts data into the file:

```julia
julia> line1 = "This here's the wattle,\n";

julia> write(fout, line1)
24
```

The return value is the number of characters that were written. The file object keeps track of where it is, so if you call write again, it adds the new data to the end of the file:

```julia
julia> line2 = "the emblem of our land.\n";

julia> write(fout, line2)
24
```

When you are done writing, you should close the file:

```julia
julia> close(fout)
```
If you don’t close the file, it gets closed for you when the program ends.
"""

# ╔═╡ 6b77c68c-8a7a-11eb-2ba0-87a6effbfd81
md"""## Formatting

The argument of `write` has to be a string, so if we want to put other values in a file, we have to convert them to strings. The easiest way to do that is with `string` or string interpolation:

```julia
julia> fout = open("output.txt", "w")
IOStream(<file output.txt>)
julia> write(fout, string(150))
3
```

An alternative is to use the print or print(ln) functions:

```julia
julia> camels = 42
42
julia> println(fout, "I have spotted $camels camels.")
```

!!! tip
    A more powerful alternative is the `@printf` macro of the `Printf` package, which prints using a C-style format specification string.
"""

# ╔═╡ e880ce50-8a7a-11eb-3d3c-31a12cd844d3
md"""## Filenames and Paths

Files are organized into *directories* (also called “folders”). Every running program has a “current directory,” which is the default directory for most operations. For example, when you open a file for reading, Julia looks for it in the current directory.

The function `pwd` returns the name of the current directory:

```julia
julia> cwd = pwd()
"/home/ben"
```

`cwd` stands for “current working directory.” The result in this example is */home/ben*, which is the home directory of a user named *ben*.

A string like `"/home/ben"` that identifies a file or directory is called a *path*.

A simple filename like *memo.txt* is also considered a path, but it is a *relative path* because it relates to the current directory. If the current directory is */home/ben*, the filename memo.txt would refer to */home/ben/memo.txt*.

A path that begins with */* does not depend on the current directory; it is called an *absolute path*. To find the absolute path to a file, you can use `abspath`:

```julia
julia> abspath("memo.txt")
"/home/ben/memo.txt"
```

Julia provides other functions for working with filenames and paths. For example, `ispath` checks whether a file or directory exists:

```julia
julia> ispath("memo.txt")
true
```

If it exists, `isdir` checks whether it’s a directory:

```julia
julia> isdir("memo.txt")
false
julia> isdir("/home/ben")
true
```

Similarly, `isfile` checks whether it’s a file.

`readdir` returns an array of the files (and other directories) in the given directory:

```julia
julia> readdir(cwd)
3-element Array{String,1}:
 "memo.txt"
 "music"
 "photos"
```

To demonstrate these functions, the following example “walks” through a directory, prints the names of all the files, and calls itself recursively on all the directories:

```julia
function walk(dirname)
	for name in readdir(dirname)
		path = joinpath(dirname, name)
		if isfile(path)
			println(path)
		else
			walk(path)
		end
	end
end
```

`joinpath` takes a directory and a filename and joins them into a complete path.

!!! tip
    Julia provides a function called `walkdir` that is similar to this one but more versatile. As an exercise, read the documentation and use it to print the names of the files in a given directory and its subdirectories.
"""

# ╔═╡ 9ebf14b8-8a7b-11eb-3347-dbc6f1972fb2
md"""## Catching Exceptions

A lot of things can go wrong when you try to read and write files. If you try to open a file that doesn’t exist, you get a `SystemError`:

```julia
julia> fin = open("bad_file")
ERROR: SystemError: opening file "bad_file": No such file or directory
```

The same thing happens if you don’t have permission to access a file:

```julia
julia> fout = open("/etc/passwd", "w")
ERROR: SystemError: opening file "/etc/passwd": Permission denied
```

To avoid these errors you could use functions like `ispath` and `isfile`, but it would take a lot of time and code to check all the possibilities.

It’s easier to go ahead and try, and deal with problems if they happen—which is exactly what the `try` statement does. The syntax is similar to an `if` statement:

```julia
try
	fin = open("bad_file.txt")
catch exc
	println("Something went wrong: $exc")
end
```

Julia starts by executing the `try` clause. If all goes well, it skips the `catch` clause and proceeds. If an exception occurs, it jumps out of the `try` clause and runs the `catch` clause.

Handling an exception with a `try` statement is called *catching* an exception. In this example, the except clause prints an error message that is not very helpful. In general, catching an exception gives you a chance to fix the problem, or try again, or at least end the program gracefully.

In code that performs state changes or uses resources like files, there is typically cleanup work (such as closing files) that needs to be done when the code is finished. Exceptions potentially complicate this task, since they can cause a block of code to exit before reaching its normal end. The `finally` keyword provides a way to run some code when a given block of code exits, regardless of how it exits:

```julia
f = open("output.txt")
try
	line = readline(f)
	println(line)
finally
	close(f)
end
```

The function `close` will always be executed.
"""

# ╔═╡ 22787ce0-8a7c-11eb-2c30-17f9710f7862
md"""## Databases

A *database* is a file that is organized for storing data. Many databases are organized like dictionaries, in the sense that they map from keys to values. The biggest difference between a database and a dictionary is that the database is on disk (or other permanent storage), so it persists after the program ends.

The `GDBM` package (`pkg"add https://github.com/BenLauwens/GDBM.jl.git"`) provides an interface to the *GDBM* (GNU dbm) library of functions for creating and updating database files. As an example, I’ll create a database that contains captions for image files.

Opening a database is similar to opening other files:

```julia
julia> using GDBM

julia> db = DBM("captions", "c")
DBM(<captions>)
```

The mode `"c"` means that the database should be created if it doesn’t already exist. The result is a database object that can be used (for most operations) like a dictionary.

When you create a new item, `GDBM` updates the database file:

```julia
julia> db["cleese.png"] = "Photo of John Cleese."
"Photo of John Cleese."
```

When you access one of the items, `GDBM` reads the file:

```julia
julia> db["cleese.png"]
"Photo of John Cleese."
```

If you make another assignment to an existing key, `GDBM` replaces the old value:

```julia
julia> db["cleese.png"] = "Photo of John Cleese doing a silly walk."
"Photo of John Cleese doing a silly walk."
julia> db["cleese.png"]
"Photo of John Cleese doing a silly walk."
```

Some functions that take a dictionary as an argument, like `keys` and `values`, don’t work with database objects. But iteration with a `for` loop works:

```julia
for (key, value) in db
	println(key, ": ", value)
end
```

As with other files, you should close the database when you are done:

```julia
julia> close(db)
```
"""

# ╔═╡ a350e348-8a7c-11eb-22a5-2b82d4b9ee63
md"""## Serialization

A limitation of `GDBM` is that the keys and the values have to be strings or byte arrays. If you try to use any other type, you get an error.

The functions `serialize` and `deserialize` can help. The `serialize` function can translate almost any type of object into a byte array (an `IOBuffer`) suitable for storage in a database:

```julia
julia> using Serialization

julia> io = IOBuffer();

julia> t = [1, 2, 3];

julia> serialize(io, t)
24
julia> print(take!(io))
UInt8[0x37, 0x4a, 0x4c, 0x07, 0x04, 0x00, 0x00, 0x00, 0x15, 0x00, 0x08, 0xe2, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]
```

The format isn’t obvious to human readers; it is meant to be easy for Julia to interpret. `deserialize` reconstitutes the object:

```julia
julia> io = IOBuffer();

julia> t1 = [1, 2, 3];

julia> serialize(io, t1)
24
julia> s = take!(io);

julia> t2 = deserialize(IOBuffer(s));

julia> print(t2)
[1, 2, 3]
```

`serialize` and `deserialize` write to and read from an `IOBuffer` object that represents an in-memory I/O stream. The function `take!` fetches the contents of the `IOBuffer` as a byte array and resets the `IOBuffer` to its initial state.

Although the new object has the same value as the old one, it is not (in general) the same object:

```julia
julia> t1 == t2
true
julia> t1 ≡ t2
false
```

In other words, `serialization` and then `deserialization` has the same effect as copying the object.

You can use this to store nonstrings in a database.
!!! tip
    In fact, the storage of nonstrings in a database is so common that
this functionality has been encapsulated in a package called `JLD`.
"""

# ╔═╡ 72f18e34-8a7d-11eb-2c21-311d92932b2a
md"""## Command Objects

Most operating systems provide a command-line interface, also known as a shell. Shells usually provide commands to navigate the filesystem and launch applications. For example, in Unix you can change directories with cd, display the contents of a directory with ls, and launch a web browser by typing (for example) **`firefox`**.

Any program that you can launch from the shell can also be launched from Julia using a *command object*:

```julia
julia> cmd = `echo hello`
`echo hello`
```

Backticks are used to delimit the command. The function `run` executes the command:

```julia
julia> run(cmd);
hello
```

`hello` is the output of the `echo` command, sent to `STDOUT`. The run function itself returns a process object, and throws an `ErrorException` if the external command fails to run successfully.

If you want to read the output of the external command, `read` can be used instead:

```julia
julia> a = read(cmd, String)
"hello\n"
```

For example, most Unix systems provide a command called `md5sum` or `md5` that reads the contents of a file and computes an [MD5 checksum](https://en.wikipedia.org/wiki/MD5). This command provides an efficient way to check whether two files have the same contents. The probability that different contents will yield the same checksum is very small (i.e., unlikely to happen before the universe collapses).

You can use a command object to run `md5` from Julia and get the result:

```julia
julia> filename = "output.txt"
"output.txt"
julia> cmd = `md5 $filename`
`md5 output.txt`
julia> res = read(cmd, String)
"MD5 (output.txt) = d41d8cd98f00b204e9800998ecf8427e\n"
```
"""

# ╔═╡ f9d3aed4-8a7d-11eb-03cb-c11448a93747
md"""## Modules
Suppose you have a file named *wc.jl* with the following code:

```julia
function linecount(filename)
	count = 0
	for line in eachline(filename)
		count += 1
	end
	count
end

print(linecount("wc.jl"))
```

If you run this program, it reads itself and prints the number of lines in the file, which is 9. You can also include it in the REPL like this:

```julia
julia> include("wc.jl")
9
```

What if you don’t want the `linecount` function to be directly available in `Main`, but you’d still like to use the function in different parts of your code? Julia introduces modules to create separate variable workspaces.

A module starts with the keyword `module` and ends with `end`. Using modules allows you to avoid naming conflicts between your own top-level definitions and those found in somebody else’s code. `import` allows you to control which names from other modules are visible, and `export` specifies which of your names are public (i.e., can be used outside the module without being prefixed with the name of the module):

```julia
module LineCount
	export linecount
	function linecount(filename)
		count = 0
		for line in eachline(filename)
			count += 1
		end
		count
	end
end
```

The `using` statement allows you to make use of a module’s public names from elsewhere, so you can use the `linecount` function that `LineCount` provides outside that module as follows:

```julia
julia> using LineCount

julia> linecount("wc.jl")
11
```

#### Exercise 14-1

Type this example into a file named *wc.jl*, include it into the REPL, and enter **`using LineCount`**.

!!! danger
    If you import a module that has already been imported, Julia does nothing. It does not reread the file, even if it has changed.

    If you want to reload a module, you have to restart the REPL. If you want to avoid this, you can use the `Revise` package to keep your sessions running longer.

"""

# ╔═╡ 8c1051c6-8a7e-11eb-313e-c321a1f68d13
md"""## Debugging

When you are reading and writing files, you might run into problems with whitespace. These errors can be hard to debug because spaces, tabs, and newlines are normally invisible:

```julia
julia> s = "1 2\t 3\n 4";

julia> println(s)
1 2 	3
 4
```

The built-in functions `repr` and `dump` can help. They take any object as an argument and return a string representation of the object:

```julia
julia> repr(s)
"\"1 2\\t 3\\n 4\""
julia> dump(s)
String "1 2\t 3\n 4"
```

This can be helpful for debugging.

One other problem you might run into is that different systems use different characters to indicate the end of a line. Some systems use a newline, represented by `\n`. Others use a Return character, represented by `\r`. Some use both. If you move files between different systems, these inconsistencies can cause problems.

For most systems, there are applications to convert from one format to another. You can read more about newline characters and conversion applications or, of course, you could write one yourself.
"""

# ╔═╡ e93ab6a2-8a7e-11eb-1631-53da88e645bc
md"""## Glossary

*persistent*:
Pertaining to a program that runs indefinitely and keeps at least some of its data in permanent storage.

*text file*:
A sequence of characters stored in permanent storage, like a hard drive.

*directory*:
A named collection of files, also called a folder.

*path*:
A string that identifies a file.

*relative path*:
A path that starts from the current directory.

*absolute path*:
A path that starts from the topmost directory in the filesystem.

*catch*:
To prevent an exception from terminating a program using the `try ... catch ... finally` statements.

*database*:
A file whose contents are organized like a dictionary with keys that correspond to values.

*shell*:
A program that allows users to type commands and then executes them by start‐ ing other programs.

*command object*:
An object that represents a shell command, allowing a Julia program to run com‐ mands and read the results.

*module*:
A separate global variable workspace used to avoid naming conflicts.
"""

# ╔═╡ 21f358fa-8a7f-11eb-0630-4f3d9dbc6869
md"""## Exercises

#### Exercise 14-2

Write a function called `sed` that takes as arguments a pattern string, a replacement string, and two filenames; it should read the first file and write the contents into the second file (creating it if necessary). If the pattern string appears anywhere in the file, it should be replaced with the replacement string.

If an error occurs while opening, reading, writing, or closing files, your program should catch the exception, print an error message, and exit.
"""

# ╔═╡ 35895b6c-8a7f-11eb-19e3-f96884b91d48
md"""#### Exercise 14-3

If you have done “Exercise 12-3”, you’ll see that a dictionary is created that maps from a sorted string of letters to the array of words that can be spelled with those letters. For example, `"opst"` maps to the array `["opts", "post", "pots", "spot", "stop", "tops"]`.

Write a module that imports `anagramsets` and provides two new functions: `storeanagrams` should store the anagram dictionary using `JLD`, and `readanagrams` should look up a word and return an array of its anagrams.
"""

# ╔═╡ 5d0288be-8a7f-11eb-0519-dd123ae3db3c
md"""#### Exercise 14-4

In a large collection of MP3 files, there may be more than one copy of the same song, stored in different directories or with different filenames. The goal of this exercise is to search for duplicates.

1. Write a program that searches a directory and all of its subdirectories, recur‐ sively, and returns an array of complete paths for all files with a given suffix (like *.mp3*).
2. To recognize duplicates, you can use `md5sum` or `md5` to compute a “checksum” for each file. If two files have the same checksum, they probably have the same contents.
3. To double-check, you can use the Unix command `diff`.
"""

# ╔═╡ Cell order:
# ╟─e89c8374-8a79-11eb-3029-7d8170bde727
# ╟─fa2d4204-8a79-11eb-0a0a-69a684f6c422
# ╟─087bcda8-8a7a-11eb-0bdf-758ab26f7fae
# ╟─23dd5184-8a7a-11eb-1dd0-35b05d248988
# ╟─6b77c68c-8a7a-11eb-2ba0-87a6effbfd81
# ╟─e880ce50-8a7a-11eb-3d3c-31a12cd844d3
# ╟─9ebf14b8-8a7b-11eb-3347-dbc6f1972fb2
# ╟─22787ce0-8a7c-11eb-2c30-17f9710f7862
# ╟─a350e348-8a7c-11eb-22a5-2b82d4b9ee63
# ╟─72f18e34-8a7d-11eb-2c21-311d92932b2a
# ╟─f9d3aed4-8a7d-11eb-03cb-c11448a93747
# ╟─8c1051c6-8a7e-11eb-313e-c321a1f68d13
# ╟─e93ab6a2-8a7e-11eb-1631-53da88e645bc
# ╟─21f358fa-8a7f-11eb-0630-4f3d9dbc6869
# ╟─35895b6c-8a7f-11eb-19e3-f96884b91d48
# ╟─5d0288be-8a7f-11eb-0519-dd123ae3db3c


text:
using Distributed
import Nanosoldier, GitHub, Sockets

nodes = Dict(Any => addprocs(1))
@everywhere import Nanosoldier

auth = GitHub.authenticate("GITHUB_AUTH00000000000000000000000000000")
secret = GITHUB_SECRET
port = Int(0xffff)

config = Nanosoldier.Config("nanosoldier-worker", nodes, auth, secret;
                            workdir = joinpath(@__DIR__, "workdir"),
                            trackrepo = "JuliaLang/julia",
                            reportrepo = "JuliaCI/NanosoldierReports",
                            testmode = false)

server = Nanosoldier.Server(config)
run(server, Sockets.IPv4(0), port)


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

const config_file_name_glob = "*-config.json"
const config_file_name_expr = r"^(.*)-config.json$"

const raw_file_suffix = Dict(
    :phy => "",
    :cal => "-calib",
    :pca => "-pcali",
)

const calib_categories = [ :cal, :pca ]

export GERDAData
mutable struct GERDAData
    config::DataConfig
end

raw_data_location(data::GERDAData, setup::Symbol) = data.config.setups[setup].data.raw
export raw_data_location

gen_data_location(data::GERDAData, setup::Symbol) = data.config.setups[setup].data.gen
export gen_data_location

meta_data_location(data::GERDAData, setup::Symbol) = data.config.setups[setup].data.meta
export meta_data_location

function data_dirname(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    runstr = filekey_run_str(key)
    if tier == :tier0
        runstr
    else
        joinpath(string(tier), string(system), string(key.category), runstr)
    end
end
export data_dirname

const tier0_timestamp_format = dateformat"yyyymmdd-HHMMSS"

function data_basename(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    if tier == :tier0
        timestamp = Dates.format(Dates.unix2datetime(key.time), tier0_timestamp_format)
        suffix = raw_file_suffix[key.category]
        "$(timestamp).events$(suffix)"
    else
        "$key-$system-$tier"
    end
end
export data_basename

function data_base_filename(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    data_location = if tier == :tier0
        raw_data_location(data, key.setup)
    else
        gen_data_location(data, key.setup)
    end

    joinpath(
        data_location,
        data_dirname(data, key, system, tier),
        data_basename(data, key, system, tier)
    )
end
export data_base_filename

function data_filename(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    if tier == :tier0
        data_base_filename(data, key, system, tier) * ".bz2"
    else
        data_base_filename(data, key, system, tier) * ".root"
    end
end
export data_filename

function log_filename(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    if tier == :tier0
        throw(ArgumentError("Can't provide log file names for tier0"))
    else
        data_base_filename(data, key, system, tier) * ".log"
    end
end
export log_filename

function calib_filename(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    if !(key.category in calib_categories)
        throw(ArgumentError("Can't generate calib file path for file key of category \"$(key.category)\""))
    end

    joinpath(
        meta_data_location(data, key.setup),
        "calib",
        filekey_run_str(key),
        "$key-$system-$tier-calib.json"
    )
end
export calib_filename

function calib_catalog_filename(data::GERDAData, setup::Symbol)
    joinpath(
        meta_data_location(data, setup),
        "calib",
        "$setup-calibrations.jsonl"
    )
end
export calib_catalog_filename

function calib_catalog(data::GERDAData, setup::Symbol)
    read(CalibCatalog, calib_catalog_filename(data, setup))
end
export calib_catalog

function calib_available_for(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    calib_available(calib_catalog(data, key.setup), key, system)
end
export calib_available_for

function calib_filename_for(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    calib_key = calibfor(calib_catalog(data, key.setup), key, system)
    calib_filename(data, calib_key, system, tier)
end
export calib_filename_for

function calib_props_for(data::GERDAData, key::FileKey, system::Symbol, tier::Symbol)
    filename = calib_filename_for(data, key, system, tier)
    read(PropDict, filename, subst_pathvar = true, subst_env = true, trim_null = true)
end
export calib_props_for

function config_filenames(data::GERDAData, key::FileKey)
    cfg_files = glob(config_file_name_glob, joinpath(meta_data_location(data, key.setup), "config"))

    function matches_key(filename)
        bname = basename(filename)
        m = match(config_file_name_expr, bname)
        if m != nothing
            occursin(key, m.captures[1])
        else
            false
        end
    end

    sort([f for f in cfg_files if matches_key(f)])
end
export config_filenames

function config_props_for(data::GERDAData, key::FileKey)
    read(PropDict, config_filenames(data, key), subst_pathvar = true, subst_env = true, trim_null = true)
end
export config_props_for


text:
using TOML
using ConfParser
import Mongoc

mutable struct MetaPackage
    name
    git_url
    dependencies
    function MetaPackage(path::String)
        package_path = "$(homedir())/.julia/registries/General/$path/"
        deps_file = package_path * "Deps.toml"
        package_file = package_path * "Package.toml"
        fh = open(package_file, "r")
        package_info = TOML.parse(fh)
        close(fh)

        git_url = package_info["repo"]
        name = package_info["name"]

        dependencies = []
        if isfile(deps_file)
            fh = open(deps_file, "r")
            deps = TOML.parse(fh)
            close(fh)
            for key in keys(deps)
                dict = deps[key]
                for key in keys(dict)
                    push!(dependencies, key)
                end
            end
        end
        new(name, git_url, dependencies)
    end
end

function db_package_save(collection, meta)
    document = Mongoc.BSON()
    document["name"] = meta.name
    document["git_url"] = meta.git_url
    document["dependencies"] = meta.dependencies
    push!(collection, document)
end

conf = ConfParse("$(homedir())/.config/julia-package-registry-tool.jl.conf")
parse_conf!(conf)
user     = retrieve(conf, "database", "user")
password = retrieve(conf, "database", "password")
url      = retrieve(conf, "database", "url")

connect_string = "mongodb+srv://$user:$password@$url"
client = Mongoc.Client(connect_string)
db = client["julia_pakages_regitry"]

# load julia package registry file
registry_file_path = "$(@__DIR__)/resources/julia_registry.toml"
fh = open(registry_file_path, "r")
registry = TOML.parse(fh)
close(fh)

# building MetaPackage Dict
packages = registry["packages"]
metas = Dict{String,MetaPackage}()
for package in values(packages)
    path = package["path"]
    meta = MetaPackage(path)
    metas[meta.name] = meta
end

# save registry to db
collection = db["packages"]
for meta in values(metas)
    db_package_save(collection, meta)
end

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_name: leroy-walton/playgrounds


text:
module OpenGL

import OpenGLStd.glEnable, OpenGLStd.glDisable

import Images.imread

export glimread

#Enabling lists of stuff.
type _GlEnable
    things::Vector
end

function glEnable(things::Vector)
    for thing in things
        glEnable(thing)
    end
    return _GlEnable(things)
end
glEnable(things...) = glEnable(things)

#Disabling lists of stuff
function glDisable(things::Vector)
    for thing in things
        glDisable(thing)
    end
end
glDisable(things...) = glDisable(things)

# Going with a basic OpenGL RGB 1D-interleaved image-texture
function glimread(filename::String)
    img = imread(filename)
    w = size(img,2)
    h = size(img,3)
    img1D = Array(Uint8,w*h*3)
    for j=1:w
        for k=1:h
            for i=1:3
                img1D[3*w*(j-1)+3*(k-1)+i] = img[i,w-j+1,k]
            end
        end
    end
    return img1D, w, h
end

end


text:
# Kth Order Markov Model

"""
    KOM{SymbolType}(context_length::Int)

Creates a K-th Order Markov Model predictor with a context depth of `context_length` for
`SymbolType`. `SymbolType` can be any valid type including `Char`, `Int` etc.,

## Examples
```julia-repl
julia> p = KOM{Char}(3)
DiscretePredictors.KOM{Char}([*] (0)
, Char[], 3)

julia> p = KOM{Int}(5)
DiscretePredictors.KOM{Int}([*] (0)
, Int[], 5)

julia> p = KOM{String}(4)
DiscretePredictors.KOM{String}([*] (0)
, String[], 4)
```

Reference:

"""
mutable struct KOM{T} <: BasePredictor{T}
    model::Trie{T,Int}
    context::Vector{T}
    cxt_length::Int

    KOM{T}( _c::Int ) where {T} = new( Trie{T,Int}(), Vector{T}(), _c )
    # (::Type{KOM{T}}){T}(_c::Int) = new{T}( Trie{T,Int}(), Vector{T}(), _c )
end

function add!( p::KOM{T}, sym::T ) where {T}
    p.model.value += 1

    # Create node path
    push!( p.context, sym )
    buffer  = p.context[1:end]
    while !isempty(buffer)
        p.model[buffer] = (haskey(p.model,buffer)) ? p.model[buffer]+1 : 1
        popfirst!( buffer )
    end

    if length(p.context) > p.cxt_length
        popfirst!( p.context )
    end

    nothing
end

function predict( p::KOM{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    buffer  = Vector{T}()
    for i = length(p.context):-1:1
        pushfirst!( buffer, p.context[i] )
        # println( "    Buffer: ", buffer, " -> ", p.model[buffer] );
        # Apply escape probability
        esc_prob = 1/p.model[buffer]
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        list_of_children    = children( p.model, buffer )
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            symbols[k] += list_of_children[k].value/p.model[buffer]
        end
    end

    return symbols;
end

function info_string( p::KOM{T} ) where {T}
    return @sprintf( "KOM(%d)", p.cxt_length )
end

function unique_string( p::KOM{T} ) where {T}
    return @sprintf( "KOM_%02d", p.cxt_length )
end


text:
# # Replication of Mackowiak, Matejka and Wiederholt (2018)

# This example replicates [Mackowiak, Matejka and Wiederholt (2018)](https://www.sciencedirect.com/science/article/abs/pii/S002205311830139X) using the [DRIPs](https://github.com/afrouzi/DRIPs.jl) package.

# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/afrouzi/DRIPs.jl/binder?filepath=examples) to run and modify the following code (no software is needed on the local machine).

# See [Afrouzi and Yang (2020)](http://www.afrouzi.com/dynamic_inattention.pdf) for background on the theory.

# ## Contents
# * [Ex. 1: AR(2) Process](@ref ex1)
#     * [Initialization](@ref ex1_initialize)
#     * [Replication of Figures (2) and (3)](@ref ex1_solution)
#     * [Measure Performance](@ref ex1_performance)
# * [Ex. 2: AR(3) and ARMA(2,1) Processes](@ref ex2)
#     * [Solve and Measure Performance](@ref ex2_solve)
#     * [Replication of Figure (5)](@ref ex2_fig5)
# * [Ex. 3: Price-setting with Rational Inattention](@ref ex3)
#     * [The Case with No Strategic Complementarity](@id ex3_no_strcomp)
#     * [The Case with Strategic Complementarity](@id ex3_strcomp)
#     * [Replication of Figure (6)](@ref ex3_fig6)
# * [Ex. 4: Business Cycle Model with News Shocks](@id ex4)
#     * [Setup](@ref ex4_setup)
#     * [Initialization](@ref ex4_initialize)
#     * [Solve and Measure Performance](@ref ex4_solve)
#     * [Replication of Figure (7)](@ref ex4_fig7)

# ## [Ex. 1: AR(2) Process](@id ex1)
# In this example, authors assume that the optimal action follows an AR(2) process,
# ```math
# \begin{aligned}
#     x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \theta_0 \varepsilon_t
# \end{aligned}
# ```
# Then, we can write a state-space form as:
# ```math
# \begin{aligned}
#         \left[\begin{array}{c} x_t \\ x_{t-1} \end{array}\right]
#         & =
#        \underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}
#             \phi_1 & \phi_2 \\
#             1 & 0\end{array}\right]}}\,
#             \left[\begin{array}{c} x_{t-1} \\ x_{t-2} \end{array}\right]
#         + \underset{\mathbf{Q}}{\underbrace{\left[\begin{array}{c}
#             \theta_0 \\
#             0\end{array}\right]}}\, \varepsilon_t
# \end{aligned}
# ```
# We now characterize the optimal signal, $S_{t} = h_1 X_t + h_2 X_{t-1} + \psi_t$, as a function of $\phi_2$ and the capacity for information processing ($\kappa$).

# ### [Initialization](@id ex1_initialize)
# Include the package:
using DRIPs, LinearAlgebra;

# Assign value to deep parameters and define the structure of the problem

β   = 1.0   ; # Time preference
θ0  = 1.0   ;

ϕ2_seq  = 0:0.02:0.9        ;      # sequence of values for the AR(2) parameter
n_ϕ     = length(ϕ2_seq)    ;

κ_seq   = 0.008:0.01:3.5    ;      # sequence of values for the information processing capacity
n_κ     = length(κ_seq)     ;

H   = [1; 0]    ;
Q   = [θ0 ; 0]    ;

# ### [Replication of Figures (2) and (3)](@id ex1_solution)
# Solve for different values of $\phi_2$:

κ   = 0.2   ; # fix κ for this exercise so we can vary ϕ2

h1  = zeros(n_ϕ,1)    ;
h2  = zeros(n_ϕ,1)    ;
h2norm = zeros(n_ϕ,1) ;

for i in 1:n_ϕ
    ϕ2  = ϕ2_seq[i] ;
    ϕ1  = 0.9 - ϕ2 ;
    A   = [ϕ1 ϕ2 ; 1.0 0.0] ;

    ex1a = Drip(κ,β,A,Q,H,fcap=true);

    h1[i] = ex1a.ss.Y[1]
    h2[i] = ex1a.ss.Y[2]
    h2norm[i] = ex1a.ss.Y[2]/ex1a.ss.Y[1] # normalize the first signal weight h1 to one
end

# Here is the replication of Figure (2):
using Printf, LaTeXStrings, Plots; pyplot();
plot(ϕ2_seq,h2norm[:,1],
    title       = L"Optimal Signal, AR(2) processes with $\phi_1 + \phi_2 = 0.9$, $\kappa$ constant.",
    xlabel      = L"$\phi_2$",
    ylabel      = L"$h_2$(with $h_1$ normalized to 1)",
    label       = L"Signal Weight on $X_{t-1}$, $h_2$ (with $h_1$ normalized to 1)",
    legend      = :topleft,
    lw          = 2,
    color       = :black,
    xlim        = (0,0.9),
    xticks      = (0:0.1:0.9),
    titlefont   = font(10), legendfont = font(7), tickfont = font(7),
    size        = (650,300),
    grid        = :off, framestyle = :box)

# Now, we solve for the signal weight and standard deviation of noise, for different values of $\kappa$

ϕ2   = 0.4           ;  # fix ϕ2 in this exercise so we can vary κ
ϕ1   = 0.99 - ϕ2     ;

A    = [ϕ1 ϕ2 ; 1.0 0.0] ;

ω_sol   = zeros(n_κ,1) ;
σz_sol  = zeros(n_κ,1) ;

for i in 1:n_κ
    κ = κ_seq[i] ;

    ex1b = Drip(κ,β,A,Q,H,fcap=true);

    h1temp = ex1b.ss.Y[1]
    h2temp = ex1b.ss.Y[2]

    ω_sol[i]  = 1.0 - (h2temp/h1temp)/(ϕ2 + (1.0-ϕ1)*(h2temp/h1temp)) ;
    ρ = (ω_sol[i] + (1.0-ω_sol[i])*ϕ1)/h1temp
    σz_sol[i] = ρ*ex1b.ss.Σ_z[1,1] ;
end

# Below is the replication of Figure (3):
a = fill(NaN, n_κ, 1)

plot(κ_seq,[ω_sol[:,1] a],
    xlabel      = L"$\kappa$",
    label       = [L"Signal weight on $X_{t}$, $\omega$ (left axis)" L"St. dev. of noise  $\sigma_{\psi}$ (right axis)"],
    title       = L"Optimal signal as function of $\kappa$, AR(2) example.",
    linestyle   = [:solid :dash], color = [:black :gray40],
    xticks      = (0:1:3),
    xlim        = (-0.02,3.02),
    yticks      = (0:0.5:1),
    ylim        = (-0.01,1.01),
    lw          = 2, grid = :off, legend = :right,
    titlefont   = font(10), legendfont = font(8), tickfont = font(8),
    framestyle  = :box)
plot!(twinx(),κ_seq,σz_sol[:,1],
    linestyle   = :dash, color = :gray40,
    label       = "",
    xlim        = (-0.02,3.02),
    xticks      = false,
    yticks      = (0:6:12),
    ylim        = (-0.12,12.12),
    lw          = 2, tickfont = font(7),
    size        = (650,300),
    grid        = :off, framestyle  = :box)

# ### [Measure Performance](@id ex1_performance)
# Benchmark the solution for random values of $\phi_2$:
using BenchmarkTools;
@benchmark Drip(κ,β,[0.9-ϕ2 ϕ2; 1.0 0.0],Q,H,fcap = true) setup = (ϕ2 = 0.9*rand())

# Benchmark the solution for random values of $κ$:
@benchmark Drip(κ,β,[0.5 0.4; 1.0 0.0],Q,H,fcap = true) setup = (κ = rand())

# ## [Ex. 2: AR(3) and ARMA(2,1) Processes](@id ex2)

# Here, we replicate the AR(3) and ARMA(2,1) examples in Mackowiak, Matejka and Wiederholt (2018).
# Using a similar method to the AR(2) case, we can write the law of motion of optimal actions as a state-space form.

# ### [Solve and Meausre Performance](@id ex2_solve)

# #### The AR(3) Example

# Initialize:

β       = 1.0           ;
θ0      = 1.0           ;
κ       = 10.8          ;

ϕ1      = 1.5           ;
ϕ2      = -0.9          ;
ϕ3      = 0.1           ;

A   = [ϕ1 ϕ2 ϕ3 ; 1 0 0 ; 0 1 0] ;
Q   = [θ0 ; 0; 0]   ;
H   = [1  ; 0; 0]   ;

# Solve:

ex2a  = Drip(κ,β,A,Q,H;tol=1e-8);

# Measure performance for different values of $\kappa$:

@benchmark Drip(κ,β,A,Q,H;tol=1e-8) setup = (κ = 10*rand())

# To be consistent with MMW(2018), scale the signal vector so that the weight on the first element is one ($h_1 = 1$):

h1 = 1;
h2 = ex2a.ss.Y[2]/ex2a.ss.Y[1]
h3 = ex2a.ss.Y[3]/ex2a.ss.Y[1]

# Print the weights:

s = @sprintf("  h1 = %5.3f, h2 = %5.3f, h3 = %5.3f", h1, h2, h3)  ;
println(s) ;

# Since we scaled the signal vector, we also need to adjust the noise in the signal accordingly:
AdjNoise   = 3.879
AdjPara_ex2a= AdjNoise*ex2a.ss.Y[1]

# Calculate IRFs:

irf_ex2a = irfs(ex2a; T=30) ;
xirf_ex2a   = irf_ex2a.x[1,1,:]
xhatirf_ex2a= irf_ex2a.x_hat[1,1,:]

x     = zeros(3,30);
xhat_noise_ex2a = zeros(3,30);

for ii in 1:30
    if ii==1
        xhat_noise_ex2a[:,ii] = ex2a.ss.K;
    else
        xhat_noise_ex2a[:,ii] = A*xhat_noise_ex2a[:,ii-1]+(ex2a.ss.K*ex2a.ss.Y')*(x[:,ii]-A*xhat_noise_ex2a[:,ii-1]);
    end
end
xhat_noise_ex2a = xhat_noise_ex2a[1,:]*AdjPara_ex2a ;

# #### The ARMA(2,1) Example:

# Initialize:

β       = 1.0           ;
θ0      = 0.5           ;
θ1      = -0.1          ;
κ       = 1.79          ;

ϕ1      = 1.3           ;
ϕ2      = -0.4          ;

A   = [ϕ1 ϕ2 θ1 ; 1.0 0.0 0.0 ; 0.0 0.0 0.0] ;
Q   = [θ0 ; 0.0; 1.0]   ;
H   = [1.0; 0.0; 0.0]   ;

# Solve:

ex2b = Drip(κ,β,A,Q,H,tol=1e-8);

# Measure performance for different values of $\kappa$:

@benchmark Drip(κ,β,A,Q,H;tol=1e-8) setup = (κ = 2*rand())

# To be consistent with MMW(2018), scale the signal vector so that the weight on the first element is one ($h_1 = 1$):

h1 = 1;
h2 = ex2b.ss.Y[2]/ex2b.ss.Y[1]
h3 = ex2b.ss.Y[3]/ex2b.ss.Y[1]

# Print the weights:

s = @sprintf("  h1 = %5.3f, h2 = %5.3f, h3 = %5.3f", h1, h2, h3)  ;
println(s) ;

# Since we scaled the signal vector, we also need to adjust the noise in the signal accordingly:

AdjNoise    = 1.349
AdjPara_ex2b= AdjNoise*ex2b.ss.Y[1]

# Calculate the IRFs

irf_ex2b     = irfs(ex2b; T=30) ;
xirf_ex2b    = irf_ex2b.x[1,1,:]
xhatirf_ex2b = irf_ex2b.x_hat[1,1,:]

x     = zeros(3,30);
xhat_noise_ex2b = zeros(3,30);

for ii in 1:30
    if ii==1
        xhat_noise_ex2b[:,ii] = ex2b.ss.K;
    else
        xhat_noise_ex2b[:,ii] = A*xhat_noise_ex2b[:,ii-1]+(ex2b.ss.K*ex2b.ss.Y')*(x[:,ii]-A*xhat_noise_ex2b[:,ii-1]);
    end
end
xhat_noise_ex2b = xhat_noise_ex2b[1,:]*AdjPara_ex2b ;

# ### [Replication of Figure (5)](@id ex2_fig5)

p1 = plot(1:30,[xirf_ex2a xhatirf_ex2a xhat_noise_ex2a],
    title       = "AR(3) example, optimal signal with i.i.d. noise",
    color             = [:black :gray20 :gray40], markerstrokecolor = [:black :gray20 :gray40],
    yticks      = (-0.8:0.4:1.6),
    ylim        = (-0.82,1.62))
p2 = plot(1:30,[xirf_ex2b xhatirf_ex2b xhat_noise_ex2b],
    title       = "ARMA(2,1) example, optimal signal with i.i.d. noise",
    color       = [:black :gray20 :gray40], markerstrokecolor = [:black :gray20 :gray40],
    yticks      = (0:0.1:0.6),
    ylim        = (-0.05,0.65))
Plots.plot(p1,p2,
    layout      = (1,2),
    label       = [L"$X_{t}$ to $\varepsilon_t$" L"$Y_t=E[X_t|I_t]$ to $\varepsilon_t$" L"$Y_t=E[X_t|I_t]$ to $\psi_t$"],
    marker      = [:circle :circle :star8], markercolor = [:black :false :gray40], markersize  = [3 5 5],
    legend      = :topright,
    xticks      = (0:2:30),
    xlim        = (0,30),
    lw          = 1.5,
    legendfont  = font(8), guidefont=font(9), titlefont=font(9), tickfont=font(8),
    size        = (700,300),
    grid        = :off, framestyle = :box)

# ## [Ex. 3: Price-setting with Rational Inattention](@id ex3)

# We now replicate the price-setting exercise in MMW (2018) and its comparison with the Woodford (2002) example. This corresponds to Figure (6) in their paper.
# The model structure is identifcal to our [Example 1](https://github.com/afrouzi/DRIPs.jl/blob/master/examples/notebooks/ex1_pricing_pe_nofeedback.ipynb) (without strategic complementarity) and [Example 2](https://github.com/afrouzi/DRIPs.jl/blob/master/examples/notebooks/ex2_pricing_pe_with_feedback.ipynb) (with strategic complementarity).

# ### [The Case with No Strategic Complementarity](@id ex3_no_strcomp)

# Initialize:
ρ   = 0.9;        #persistence of money growth
σ_u = 0.1;          #std. deviation of shocks to money growth

# Primitives of the DRIP:
κ   = 0.62;
β   = 1.0;
A   = [1 ρ; 0 ρ];
Q   = σ_u*[1; 1];
H   = [1; 0];
ω   = 0.1

# Solve:

ex_opt   = Drip(κ,β,A,Q,H,fcap=true);
capa_opt = DRIPs.capacity(ex_opt) #returns capacity utilized in bits

# Measure performance for different values of $\kappa$:

@benchmark Drip(κ,β,A,Q,H,fcap=true) setup = (κ = 2*rand())

# Calculate IRFs:

irfs_ex_opt = irfs(ex_opt, T = 12);
output_opt  = (irfs_ex_opt.x[1,1,:] - irfs_ex_opt.a[1,1,:]) ;
output_opt  = [0;output_opt] ;

# Now, to compare with Woodford (2002), assume that firms observe a noisy signal, $S_t = q_t + \zeta_t$ where $\zeta_t$ is an idiosyncratic noise. We first define a function to solve the corresponding Kalman filtering problem.

function K_filtering(A,Q,Ysignal,Σz,Σ0 ; maxit=10000,tol=1e-10,w=1)
    err = 1
    iter = 0
    while (err > tol) & (iter < maxit)
        global Knew = Σ0*Ysignal*inv(Ysignal'*Σ0*Ysignal .+ Σz)
        global Σp_temp = Σ0 - Knew*Ysignal'*Σ0
        global Σ1 = A*Σp_temp*A' + Q*Q'

        err     = norm(Σ1 - Σ0,2)/norm(Σ0,2)
        Σ0 = w*Σ1 + (1-w)*Σ0

        #println("Iteration $iter. Difference: $err")
        iter += 1
    end

    return(Knew,Σ0,Σp_temp)
end;

# Now find the capacity utilized under Woodford’s formulation such that it yields the same information flow as the optimal signal under ratinoal inattention.

Ywoodford = [1;0]
Σ1_init = ex_opt.ss.Σ_1

Σz_new_b  = 0.01
Σz_new_u  = 0.1
Σz_new    = (Σz_new_b+Σz_new_u)/2
for i in 1:10000
    (Knew,Σ1_new,Σp_temp) = K_filtering(A,Q,Ywoodford,Σz_new,Σ1_init;w=0.5)
    capa_woodford = 0.5*log(det(Σ1_new)/det(Σp_temp))/log(2)

    if capa_woodford > capa_opt
        global Σz_new_b  = Σz_new
    else
        global Σz_new_u  = Σz_new
    end
    global Σz_new    = (Σz_new_b+Σz_new_u)/2
    err = abs(capa_woodford - capa_opt)
    if err < 1e-5
        break
    end
end

# Calculate impulse responses under Woodford's formulation
e_k = 1;
x   = zeros(2,12);
xhat= zeros(2,12);
a   = zeros(2,12);

for ii in 1:12
    if ii==1
        x[:,ii]     = Q*e_k;
        xhat[:,ii]  = (Knew*Ywoodford')*(x[:,ii]);
    else
        x[:,ii]     = A*x[:,ii-1];
        xhat[:,ii]  = A*xhat[:,ii-1]+(Knew*Ywoodford')*(x[:,ii]-A*xhat[:,ii-1]);
    end
    a[:,ii]  .= H'*xhat[:,ii];
end

output_woodford  = (x[1,:] - a[1,:]) ;
output_woodford  = [0;output_woodford] ;

# Before plotting the IRFs, we also solve the case with strategic complementarity.

# ### [The Case with Strategic Complementarity](@id ex3_strcomp)

# We now turn to the example with strategic complementarity.
# As in our [Example 2](https://github.com/afrouzi/DRIPs.jl/blob/master/examples/notebooks/ex2_pricing_pe_with_feedback.ipynb), we first define a function to solve the fixed point with endogenous feedback.

function ge_drip(ω,β,A,Q,          #primitives of drip except for H because H is endogenous
                 α,                #strategic complementarity
                 Hq,               #state space rep. of Δq
                 L;                #length of truncation
                 H0       = Hq,    #optional: initial guess for H (Hq is the true solution when α=0)
                 maxit    = 200,   #optional: max number of iterations for GE code
                 tol      = 1e-4)  #optional: tolerance for iterations
    err   = 1;
    iter  = 0;
    M     = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
    while (err > tol) & (iter < maxit)
        if iter == 0
            global ge  = Drip(ω,β,A,Q,H0; w=0.9);
        else
            global ge  = Drip(ω,β,A,Q,H0; Ω0=ge.ss.Ω, Σ0=ge.ss.Σ_1, maxit=100);
        end

        XFUN(jj) = ((I-ge.ss.K*ge.ss.Y')*ge.A)^jj * (ge.ss.K*ge.ss.Y') * (M')^jj
        X = DRIPs.infinitesum(XFUN; maxit=L, start = 0);  #E[x⃗]=X×x⃗

        XpFUN(jj) = α^jj * X^(jj)
        Xp = DRIPs.infinitesum(XpFUN; maxit=L, start = 0);

        H1 = (1-α)*Xp'*Hq;
        err= 0.5*norm(H1-H0,2)/norm(H0)+0.5*err;
        H0 = H1;

        iter += 1;
        if iter == maxit
            print("GE loop hit maxit\n")
        elseif mod(iter,10) == 0
            println("Iteration $iter. Difference: $err")
        end

    end
    print(" Iteration Done.\n")
    return(ge)
end;

# Now, we solve for the optimal signal structure under rational inattention.

# Initialize:

ρ   = 0.9;        #persistence of money growth
σ_u = 0.1;        #std. deviation of shocks to money growth
α   = 0.85;        #degree of strategic complementarity
L   = 40;         #length of truncation
Hq  = ρ.^(0:L-1); #state-space rep. of Δq

ω   = 0.08;
β   = 1 ;
A   = [1 zeros(1,L-2) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
M   = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)]; # shift matrix
Q   = [σ_u; zeros(L-1,1)];

# Solve:

ex_ge   = ge_drip(ω,β,A,Q,α,Hq,L) ;

# Print capacity utilized in with strategic complementarity:

capa_ge = DRIPs.capacity(ge)

# Measure performance for random values of ω

using Suppressor
@suppress @benchmark ge_drip(ω,β,A,Q,α,Hq,L) setup = (ω = 0.1*rand())

# Calculate IRFs
geirfs  = irfs(ex_ge,T = L) ;

dq = diagm(Hq)*geirfs.x[1,1,:];
q  = inv(I-M)*dq ;
output_ge_opt  = q - geirfs.a[1,1,:] ;
output_ge_opt  = [0;output_ge_opt] ;

# Finally, to compare with the IRFs under Woodford (2002)’s specification, we find signal noise such that it yields the same information flow as the optimal signal structure.

Ywoodford_ge = Hq
Σ1_init  = ex_ge.ss.Σ_1

Σz_new_b = 0.05
Σz_new_u = 0.12
Σz_new   = (Σz_new_b+Σz_new_u)/2

for i in 1:10000
    (Knew,Σ1_new,Σp_temp) = K_filtering(A,Q,Ywoodford_ge,Σz_new,Σ1_init;w=0.5)
    capa_woodford = 0.5*log(det(Σ1_new)/det(Σp_temp))/log(2)

    if capa_woodford > capa_ge
        global Σz_new_b  = Σz_new
    else
        global Σz_new_u  = Σz_new
    end
    global Σz_new    = (Σz_new_b+Σz_new_u)/2
    err = abs(capa_woodford - capa_ge)
    if err < 1e-5
        break
    end
end

XFUN(jj) = ((I-Knew*Ywoodford_ge')*A)^jj * (Knew*Ywoodford_ge') * (M')^jj
X = DRIPs.infinitesum(XFUN; maxit=L, start = 0);  #E[x⃗]=X×x⃗

XpFUN(jj) = α^jj * X^(jj)
Xp = DRIPs.infinitesum(XpFUN; maxit=L, start = 0);

H1 = (1-α)*Xp'*Hq;

# Calculate impurse responses under Woodford's signal
e_k = 1;
x   = zeros(40,40);
xhat= zeros(40,40);
a   = zeros(1,40);

for ii in 1:40
    if ii==1
        x[:,ii]     = Q*e_k;
        xhat[:,ii]  = (Knew*Ywoodford_ge')*(x[:,ii]);
    else
        x[:,ii]     = A*x[:,ii-1];
        xhat[:,ii]  = A*xhat[:,ii-1]+(Knew*Ywoodford_ge')*(x[:,ii]-A*xhat[:,ii-1]);
    end
    a[:,ii]  .= H1'*xhat[:,ii];
end
dq = diagm(Hq)*x[1,:];
q  = inv(I-M)*dq ;
output_ge_woodford  = q - a[1,:] ;
output_ge_woodford = [0;output_ge_woodford] ;

# We now have all the IRFs to replicate Figure (6)

# ### [Replication of Figure (6)](@id ex3_fig6)

p1 = plot(0:12,[output_woodford output_opt],
    title       = L"The case of $\xi=1$",
    color       = [:black :gray40], markerstrokecolor = [:black :gray40],
    ylim        = (-0.003,0.053),
    ytick       = (0:0.01:0.05))
p2 = plot(0:12,[output_ge_woodford[1:13] output_ge_opt[1:13]],
    title       = L"The case of $\xi=0.15$",
    color       = [:black :gray40], markerstrokecolor = [:black :gray40],
    ylim        = (-0.005,0.105),
    ytick       = (0:0.02:0.1))
Plots.plot(p1,p2,
    layout      = (2,1),
    label       = ["Woodford Model" "Model with optimal signals"],
    legend      = :topright,
    marker      = [:circle :star8], markersize = [5 5], markercolor = [:false :gray40],
    xlim        = (0,12),
    xtick       = (0:1:12),
    xlabel      = "Time",
    lw          = 1.5,
    legendfont  = font(8), titlefont = font(10), guidefont = font(9), tickfont = font(8),
    size        = (700,400),
    grid        = :off, framestyle = :box)

# ## [Ex. 4: Business Cycle Model with News Shocks](@id ex4)
# In this section, we replicate the business cycle model with news shocks in Section 7 in Mackowiak, Matejka and Wiederholt (2018).
#
# ### [Setup](@id ex4_setup)
#
# #### Full-Information
# The techonology shock follows AR(1) process:
# ```math
# \begin{aligned}
# z_{t} = \rho z_{t-1} + \sigma \varepsilon_{t-k}
# \end{aligned}
# ```
# and the total labor input is:
# ```math
# \begin{aligned}
# n_{t} = \int_0^1 n_{i,t} di .
# \end{aligned}
# ```
#
# Under perfect information, the households chooses the utility-maximizing labor supply, all firms choose the profit-maximizing labor input, and the labor market clearing condition is:
# ```math
# \begin{aligned}
# \frac{1-\gamma}{\psi + \gamma}w_{t} = \frac{1}{\alpha}(z_t - w_t).
# \end{aligned}
# ```
# Then, the market clearing wages and the equilibrium labor input are:
# ```math
# \begin{aligned}
#     w_{t} & = \frac{\frac{1}{\alpha}}{\frac{1-\gamma}{\psi+\gamma} + \frac{1}{\alpha}} z_t \equiv \xi z_t \\
#     n_t & = \frac{1}{\alpha}(1-\xi) z_t.
# \end{aligned}
# ```
#
#
# #### Rational Inattention
#
# Firms wants to keep track of their ideal price,
# ```math
# \begin{aligned}
#     n_{t}^* = \frac{1}{\alpha} z_t - \frac{1}{\alpha} \frac{\psi + \gamma}{1 - \gamma} n_t
# \end{aligned}
# ```
# where $n_{t} = \int_0^1 n_{i,t} di$. Then, firm $i$'s choice depends on its information set at time $t$:
# ```math
# \begin{aligned}
#     n_{i,t} = E_{i,t} [n_{t}^*].
# \end{aligned}
# ```
#
# Note that now the state space representation for $n_{t}^*$ is determined in the equilibrium. However, we know that this is a Guassian process and by Wold's theorem we can decompose it to its $MA(\infty)$ representation:
# ```math
# \begin{aligned}
#     n_{t}^*=\Phi(L)\varepsilon_t
# \end{aligned}
# ```
# where $\Phi(.)$ is a lag polynomial and $\varepsilon_t$ is the shock to technology. Here, we have basically guessed that the process for $p_{i,t}^*$ is determined uniquely by the history of monetary shocks which requires that rational inattention errors of firms are orthogonal (See [Afrouzi (2020)](http://www.afrouzi.com/strategic_inattetion.pdf)). Our objective is to find $\Phi(.)$.
#
# Now, as in our [Example 2](https://github.com/afrouzi/DRIPs.jl/blob/master/examples/notebooks/ex2_pricing_pe_with_feedback.ipynb), we can represent the problem in a matrix notation.
#
# ### [Initialization](@id ex4_initialize)

β   = 1     ;   # Time preference
γ   = 1/3   ;   # Inverse of intertemporal elasticity of substitution
ψ   = 0     ;   # Inverse of Frisch elasticity
α   = 3/4   ;   # Labor share in production function
θ   = -1/α*(ψ+γ)/(1-γ)  ;
ξ   = θ/(θ-1)   ;

ρ   = 0.9   ;   #persistence of technology shocks
σ   = 1     ;   #std. deviation of technology shocks
ω   = 6.5   ;   # Information cost

L   = 40    ;   #length of truncation
k   = 8     ;   #news horizon

M   = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)]; # shift matrix
Hz  = ρ.^(0:L-1)
Hz = (M^k)*Hz

A   = M     ;
Q   = [σ; zeros(L-1,1)] ;

# Also, define a function that solves the GE problem and returns the solution in a `Drip` structure:

function ge_drip(ω,β,A,Q,          #primitives of drip except for H because H is endogenous
                 α,                #strategic complementarity
                 θ,
                 Hz,               #state space rep. of z
                 L;                #length of truncation
                 w_out    = 0.5,    #optional: initial guess for H (Hq is the true solution when α=0)
                 H0       = Hz,    #optional: initial guess for H (Hq is the true solution when α=0)
                 maxit    = 200,   #optional: max number of iterations for GE code
                 tol      = 1e-6)  #optional: tolerance for iterations
    err   = 1;
    iter  = 0;
    M     = [zeros(1,L-1) 0; Matrix(I,L-1,L-1) zeros(L-1,1)];
    eye   = Matrix(I,L,L)
    while (err > tol) & (iter < maxit)
        if iter == 0
            global ge  = Drip(ω,β,A,Q,H0;w=0.5, tol=1e-8);
        else
            global ge  = Drip(ω,β,A,Q,H0;w=0.9, tol=1e-8, Ω0=ge.ss.Ω, Σ0=ge.ss.Σ_1, maxit=1000);
        end

        XFUN(jj) = ((eye-ge.ss.K*ge.ss.Y')*ge.A)^jj * (ge.ss.K*ge.ss.Y') * (M')^jj
        X = DRIPs.infinitesum(XFUN; maxit=L, start = 0);  #E[x⃗]=X×x⃗

        H1 = (1/α)*Hz + θ*X'*H0 ;

        err= 0.5*norm(H1-H0,2)/norm(H0)+0.5*err;

        H0 = w_out*H1 + (1.0-w_out)*H0 ;

        iter += 1;
        if iter == maxit
            print("GE loop hit maxit\n")
        elseif mod(iter,10) == 0
            println("Iteration $iter. Difference: $err")
        end

    end
    print(" Iteration Done.\n")
    return(ge)
end;

# ### [Solve and Measure Performance](@id ex4_solve)

# Solve:

ge      = ge_drip(ω,β,A,Q,α,θ,Hz,L)  ;

# Measure performance by solving the model for different values of ω:

@suppress @benchmark ge_drip(ω,β,A,Q,α,θ,Hz,L) setup = (ω = 6.5*rand())

# Calculate IRFs and profit loss:

geirfs  = irfs(ge,T = L)         ;
profit_loss = sum((geirfs.x[1,1,:]/100 - geirfs.x_hat[1,1,:]/100).^2) ;

s = @sprintf("==: Profit loss from rational inattention = %6.5f", profit_loss)  ;
println(s) ;

# ### [Replication of Figure (7)](@id ex4_fig7)

n_opt   = geirfs.a[1,1,:]   ; # Optimal labor input under rational inattention
n_fullinfo = σ*1/α*(1-ξ)*Hz ; # Optimal labor input under full information

plot(1:30,[n_fullinfo[1:30] n_opt[1:30]],
    title       = "The impulse response of labor input to a productivity shock.",
    ylabel      = "Percent",
    label       = ["Equilibrium under perfect information" "Equilibirum when firms are subject to rational inattention"],
    legend      = :topright,
    color       = [:black :gray40], markerstrokecolor = [:black :gray40],
    marker      = [:circle :star8], markercolor = [:false :gray40], markersize = [5 5],
    ylim        = (-0.05,0.85),
    ytick       = (0:0.1:0.8),
    xlim        = (0,30),
    xtick       = (0:2:30),
    lw          = 1.5,
    legendfont  = font(8), titlefont = font(10), tickfont = font(8), guidefont = font(9),
    size        = (650,400),
    grid        = :off, framestyle = :box)


text:
import FLOWFarm; const ff = FLOWFarm

# based on https://backend.orbit.dtu.dk/ws/portalfiles/portal/146793996/Pena18_wes_3_191_2018.pdf

diam = 120.0

locdata = readdlm("inputfiles/AnholtOffshoreWindFarmLocations.txt", ' ', skipstart=1)
turbine_x = locdata[:, 1]*diam
turbine_y = locdata[:, 2]*diam
nturbines = length(turbine_x)

turbine_z = zeros(nturbines)
rotor_diameter = zeros(nturbines).+diam
hub_height = zeros(nturbines).+81.6
turbine_yaw = zeros(nturbines).+0.0
turbine_ai = zeros(nturbines).+(1.0/3.0)

generator_efficiency = zeros(nturbines).+0.944
cut_in_speed = zeros(nturbines).+3.5  # m/s
cut_out_speed = zeros(nturbines).+25.  # m/s
rated_speed = zeros(nturbines).+12.0  # m/s
rated_power = zeros(nturbines).+3.6E6  # W

# load cp data
cpdata = readdlm("inputfiles/power_curve_siemens_swp3.6-120.txt",  ',', skipstart=1)
vel_points = cpdata[:, 1]
power_points = cpdata[:, 2]

power_model = ff.PowerModelPowerPoints(vel_points, power_points)
power_models = Vector{typeof(power_model)}(undef, nturbines)
for i = 1:nturbines
    power_models[i] = power_model
end

# rotor sample points
rotor_points_y = [0.0]
rotor_points_z = [0.0]

winddirections = deg2rad.(range(0,stop=350,length=36))
windspeeds = zeros(36).+(9.23)
windprobabilities = [1.3,1.1,1.3,1.1,1.0,1.3,1.7,2.2,2.1,3.3,3.5,4.1,4.5,5.0,
                            3.6,3.3,2.8,2.6,2.8,3.7,4.7,4.3,3.6,3.5,3.1,3.6,4.0,3.6,
                            4.2,3.8,1.8,1.6,1.6,1.4,1.2,1.3]
windprobabilities = windprobabilities./(sum(windprobabilities))

air_density = 1.1716  # kg/m^3
ambient_ti = 0.1
measurementheight = copy(hub_height)
ambient_tis = ones(length(winddirections)).*ambient_ti
shearexponent = 0.15
wind_shear_model = ff.PowerLawWindShear(shearexponent, 0.0, "nothing")
windresource = ff.DiscretizedWindResource(winddirections, windspeeds, windprobabilities, measurementheight, air_density, ambient_tis, wind_shear_model)

alpha = 0.04
wakedeficitmodel = ff.JensenTopHat(alpha)
wakedeflectionmodel = ff.JiminezYawDeflection(alpha)
wakecombinationmodel = ff.SumOfSquaresFreestreamSuperposition()
localtimodel = ff.LocalTIModelNoLocalTI()
model_set = ff.WindFarmModelSet(wakedeficitmodel, wakedeflectionmodel, wakecombinationmodel, localtimodel)

ct = 0.689
ct_model = ff.ThrustModelConstantCt(ct)
ct_models = Vector{typeof(ct_model)}(undef, nturbines)
for i = 1:nturbines
    ct_models[i] = ct_model
end


text:
"""
    PODsvd(X)

Uses the SVD based decomposition technique to calculate the POD basis of X.
"""
function PODsvd(X; subtractmean::Bool = false)

    Xcop = deepcopy(X)
    PODsvd!(Xcop, subtractmean = subtractmean)

end

"""
    PODsvd(X,W)

Same as `PODsvd(X)` but uses weights for each data point. The weights are equal to the
cell volume for a volume mesh.
"""
function PODsvd(X,W::AbstractVector; subtractmean::Bool = false)

    Xcop = deepcopy(X)
    PODsvd!(Xcop, W, subtractmean = subtractmean)

end

"""
    PODsvd!(X)

Same as `PODsvd(X)` but overwrites memory.
"""
function PODsvd!(X; subtractmean::Bool = false)

    if subtractmean
        X .-= mean(X,dims=2)
    end

    # Economy sized SVD
    F = svd!(X)

    # Mode coefficients
    a = Diagonal(F.S)*F.Vt

    POD = PODBasis(a, F.U)

    return POD, F.S

end

"""
    PODsvd!(X,W)

Same as `PODsvd!(X)` but uses weights for each data point. The weights are equal to the
cell volume for a volume mesh.
"""
function PODsvd!(X,W::AbstractVector; subtractmean::Bool = false)

    if subtractmean
        X .-= mean(X,dims=2)
    end

    # Take into account the weights
    X = sqrt.(W).*X

    # Economy sized SVD
    F = svd!(X)

    # Mode coefficients
    a = Diagonal(F.S)*F.Vt
    Φ = 1 ./sqrt.(W).*F.U

    POD = PODBasis(a, Φ)

    return POD, F.S

end


text:
function double_integrator_dynamics!(ẋ,x,u)
    ẋ[1] = x[2]
    ẋ[2] = u[1]
end

n = 2
m = 1

model = Model(double_integrator_dynamics!,n,m)

# initial and goal states
x0 = [1.;0.]
xf = [0.;0.]

# costs
Q = (1e-2)*Diagonal(I,n)
Qf = 100.0*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

# simulation
tf = 5.0
dt = 0.1

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)

double_integrator = [model, obj_uncon]


text:
include("kernels.jl)
a = CuArray(rand(Float32, (3, 4)))
b = CuArray(zeros(Float32, (3, 4)))

display(a)
softmax!(a, b)
display(b)

a = CuArray(rand(Float32, (16384, 64*32)));
b = CuArray(zeros(Float32, (16384, 64*32)));
@btime begin softmax!(a, b); synchronize() end
# 795.205 μs (24 allocations: 704 bytes)