text:
# agent-based economy under covid infection

using Agents,Random,Plots,Serialization, ProgressMeter, Statistics

include("covid_model.jl")

include("covid_snapshot_loader.jl")

path_to_data = "" # for example: data//bailout//duration//12//
include("$(path_to_data)covid_par_ini.jl")

if loadsnapshot
    covidmodel, lochh, locf, unemplist, shorttimelist, empcount, unempcount, shorttimecount, oldcount, unemp, firms, hh, tau, divperhh, weeklyconsumption = restore_snapshot("$snapname")
else
    include("covid_rand_ini.jl")
end

actshoplist = deepcopy(locf) # initialize list of active stores

datapoint = fld(T,datat)
datahhy = zeros(datapoint+1, nhhy, 5) #pos,cor, emp, saving, consbud
datahho = zeros(datapoint+1, nhho, 5) #pos,cor, emp, saving, consbud
dataf = zeros(datapoint+1, nf, 6) #pos,saving,product,stock, demand
castraj = zeros(datapoint+1,k1,k2,2)
unempltraj = zeros(datapoint+1)
shorttimetraj = zeros(datapoint+1)
unempsectraj = zeros(datapoint+1,nsec)
shorttimesectraj = zeros(datapoint+1,nsec)
pubacctraj = zeros(datapoint+1)
tautraj = zeros(datapoint+1)
totfirmtraj = zeros(datapoint+1)
constraj = zeros(datapoint+1,nsec,k1,k2) # total consuption during the data collection interval
conspercaptraj = zeros(datapoint+1)
gdppercaptraj = zeros(datapoint+1)
gdplosstraj = zeros(datapoint+1)
totgdploss = 0
inactivefirmstraj = zeros(datapoint+1)
bailouttraj = zeros(datapoint+1)
totalaccountstraj = zeros(datapoint+1)
R0count = Int[]  # list of number infected by infected agents
R0counttraj = zeros(datapoint+1)
f_av_savings_reg = zeros(datapoint+1,nsec,k1,k2)
f_av_workers_reg = zeros(datapoint+1,nsec+1,k1,k2)
f_av_stock_reg = zeros(datapoint+1,nsec,k1,k2)
f_av_demandexp_reg = zeros(datapoint+1,nsec,k1,k2)
f_av_savings = zeros(datapoint+1,nsec)
f_var_savings = zeros(datapoint+1,nsec)
f_av_workers = zeros(datapoint+1,nsec+1)
f_var_workers = zeros(datapoint+1,nsec+1)
f_av_stock = zeros(datapoint+1,nsec)
f_av_demandexp = zeros(datapoint+1,nsec)
f_var_demandexp = zeros(datapoint+1,nsec)
bankrupttraj = zeros(datapoint+1,nsec)
totinftraj = []
RKIR0traj = zeros(datat*datapoint+1)
curinftraj = []
curinfregtraj = []
contact_count_traj = [] # list of number of contatcs
contact_work_traj = [] # list of number of contatcs in work
contact_social_traj = [] # list of number of social contatcs
contact_shop_traj = [] # list of number of contatcs shopping
bankruptcount = zeros(nsec)

# initial datapoint
for id in hh[1]
    current_agent = getagent(id)
    datahhy[1,id,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
    datahhy[1,id,2] = current_agent.cor
    datahhy[1,id,3] = current_agent.emp
    datahhy[1,id,4] = current_agent.saving
    datahhy[1,id,5] = sum(current_agent.consbud)
end
for id in hh[2]
    current_agent = getagent(id)
    datahho[1,id-nhhy,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
    datahho[1,id-nhhy,2] = current_agent.cor
    datahho[1,id-nhhy,3] = current_agent.emp
    datahho[1,id-nhhy,4] = current_agent.saving
    datahho[1,id-nhhy,5] = sum(current_agent.consbud)
end
for id in firms
    current_agent = getagent(nhhy+nhho+id)
    dataf[1,id,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
    dataf[1,id,2] = current_agent.saving
    dataf[1,id,3] = length(current_agent.worker)
    dataf[1,id,4] = current_agent.stock
    dataf[1,id,5] = current_agent.demandexp
end

#function phasepolicy(t::Int)
#end

#initial unemployment rate
unempltraj[1] = unemp
unempsectraj[1,1] = unemp
unempsectraj[1,2] = unemp
unempsectraj[1,3] = unemp
shorttimetraj[1] = 0
shorttimesectraj[1,1] = 0
shorttimesectraj[1,2] = 0
shorttimesectraj[1,3] = 0
#initial public account
pubacctraj[1] = 0
#initial number of firms
totfirmtraj[1] = size(firms)[1]
# initial tax
tautraj[1] = tau
#initially no infected
push!(totinftraj, 0)
# initial no contact
push!(contact_count_traj, 0)
push!(contact_work_traj, 0)
push!(contact_social_traj, 0)
push!(contact_shop_traj, 0)

# robert koch R0
RKIR0traj[1] =0
ttt=0

currentadaptivepolicy = "NONE"
polswitchcount = 0
badpoltime = 0

#iteration
@showprogress for t = 1:datapoint
    if t*datat+1 > virustime && trigger # activated if t first crosses virustime
        for i=1:k1
            for j=1:k2
#                infected = rand(merge(lochh[1,i,j],lochh[2,i,j]),iniinf[i][j]) # select randomly infected
#                for (id, agent) in infected
#                    agent.cor = 2 # agent is infected
#                    agent.cortime = trec
#                end
                infected = rand(lochh[1,i,j],iniinf[i][j][1]) # select randomly infected
                for (id, agent) in infected
                    agent.cor = 2 # agent is infected
                    agent.cortime = trec
                end
                infected = rand(lochh[2,i,j],iniinf[i][j][2]) # select randomly infected
                for (id, agent) in infected
                    agent.cor = 2 # agent is infected
                    agent.cortime = trec
                end
            end
        end
        global trigger = false
        global virus = true  # activate virus spreading in model_step
        global totgdploss = 0
    end

    # adjust policy variables
    global pinf += poladjfrac*(pinft - pinf)
    global socialmaxyyh += poladjfrac*(socialmaxyyt - socialmaxyyh)
    if socialmaxyyh >= 1
        global socialmaxyy = Int(ceil(socialmaxyyh))
    else
        global socialmaxyy = socialmaxyyh
    end
    global socialmaxoyh += poladjfrac*(socialmaxoyt - socialmaxoyh)
    if socialmaxoyh >= 1
        global socialmaxoy = Int(ceil(socialmaxoyh))
    else
        global socialmaxoy = socialmaxoyh
    end
    global socialmaxooh += poladjfrac*(socialmaxoot - socialmaxooh)
    if socialmaxooh >= 1
        global socialmaxoo = Int(ceil(socialmaxooh))
    else
        global socialmaxoo = socialmaxooh
    end
    global pshop += poladjfrac .* (pshopt .- pshop)
    global phomeoffice += poladjfrac*(phomeofficet - phomeoffice)

    global workmeeth = workmeeth .+ poladjfrac.*(workmeett .- workmeeth)
    global workmeet = Int.(ceil.(workmeeth))

    global shopmeeth = shopmeeth .+ poladjfrac.*(shopmeett .- shopmeeth)
    global shopmeet = Int.(ceil.(shopmeeth))

    if t*datat <= resetfirmsavingstime && (t+1)*datat > resetfirmsavingstime && !loadsnapshot  # reset firm savings after initial burn-in
        for id in firms
            getagent(nhh+id).saving = 500
        end
    end

    if t*datat < enablebankruptciestime && !loadsnapshot
        global bankruptcypossible = false
    else
        global bankruptcypossible = true
    end

    if t*datat <= endfcadj && (t+1)*datat > endfcadj && !loadsnapshot # update firm fixed costs after initial burn-in
        global updatefixedcosts  = false
    end

    for (tpol, policyfile) in policies
        if t*datat <= tpol && (t+1)*datat > tpol
            include("$policyfile")
        end
    end

    global end_bailout_after_weeks
    if t > end_bailout_after_weeks
        global bailoutprogram = false
        global shorttimeprogram = false
    end

    #step!(covidmodel, dummystep, model_step!, datat)
    noinfbow = totinftraj[size(totinftraj, 1)]
    for j = 1:datat
        global ttt +=1
        step!(covidmodel, dummystep, model_step!, 1)
        if ttt > 9 && totinftraj[ttt-8] > 0
            global RKIR0traj[ttt] = (totinftraj[ttt]-totinftraj[ttt-4]) / (totinftraj[ttt-4]-totinftraj[ttt-8])
        end
        push!(curinfregtraj,curinfregcount)
        push!(curinftraj,curinfcount)
    end
    noinfeow = totinftraj[size(totinftraj, 1)]

    newinf = noinfeow - noinfbow

    # Adaptive policy
    global currentadaptivepolicy
    if (t+1)*datat > tadaptivepolicystart && (t+1)*datat < tadaptivepolicyend
        if newinf >= adaptivepolicythreshold && (currentadaptivepolicy == "GOOD" || currentadaptivepolicy == "NONE")
            include("$adaptivepolicybad")
            currentadaptivepolicy = "BAD"
            global polswitchcount += 1
        elseif newinf < adaptivepolicythreshold && currentadaptivepolicy == "BAD"
            include("$adaptivepolicygood")
            currentadaptivepolicy = "GOOD"
        end
    end
    if currentadaptivepolicy == "BAD"
        global badpoltime += datat
    end

    if t == 1
        global gdppercapref = gdp / nhh #size(union(hh[1],hh[2]))[1]
        global firmref = zeros(nsec)
        for tt=1:nsec
             global firmref[tt] = 0
             for i = 1:k1
                 for j=1:k2
                   global firmref[tt] += size(locf[tt,i,j])[1]
                 end
              end
          end
      end

    for id in hh[1]
        current_agent = getagent(id)
        datahhy[t+1,id,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
        datahhy[t+1,id,2] = current_agent.cor
        datahhy[t+1,id,3] = current_agent.emp
        datahhy[t+1,id,4] = current_agent.saving
        datahhy[t+1,id,5] = sum(current_agent.consbud)
    end
    unempltraj[t+1] = unemp
    for tt = 1:nsec
        unempsectraj[t+1,tt] = unempcount[tt] / (size(hh[1])[1] * fracemp[tt])
    end

    for tt = 1:nsec
        shorttimetraj[t+1] += shorttimecount[tt]
        shorttimesectraj[t+1,tt] = shorttimecount[tt] / (size(hh[1])[1] * fracemp[tt])
    end

    pubacctraj[t+1] = pubacc / size(union(hh[1],hh[2]))[1]
    tautraj[t+1] = tau
    totfirmtraj[t+1] = size(firms)[1]
    R0counttraj[t+1] = mean(R0count)

    for id in hh[2]
        current_agent = getagent(id)
        datahho[t+1,id-nhhy,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
        datahho[t+1,id-nhhy,2] = current_agent.cor
        datahho[t+1,id-nhhy,3] = current_agent.emp
        datahho[t+1,id-nhhy,4] = current_agent.saving
        datahho[t+1,id-nhhy,5] = sum(current_agent.consbud)
    end
    for id in firms
        current_agent = getagent(nhh+id)
        dataf[t+1,id,1] = ((current_agent.pos[1]-1)*k2+current_agent.pos[2])
        dataf[t+1,id,2] = current_agent.saving
        dataf[t+1,id,3] = length(current_agent.worker)
        dataf[t+1,id,4] = current_agent.stock
        dataf[t+1,id,5] = current_agent.demandexp
        dataf[t+1,id,6] = 1
    end

    for tt = 1:nsec
             tempsavreg = [[] for i = 1:k1, j = 1:k2]
             tempworkreg = [[] for i = 1:k1, j = 1:k2]
             tempstockreg = [[] for i = 1:k1, j = 1:k2]
             tempdemandexpreg = [[] for i = 1:k1, j = 1:k2]
             tempsav = []
             tempwork = []
             tempstock = []
             tempdemandexp = []
             for id in firms
                 if getagent(nhh+id).type-(nsec+1) == tt
                 push!(tempsavreg[Int(dataf[t+1,id,1])],dataf[t,id,2])
                 push!(tempworkreg[Int(dataf[t+1,id,1])],dataf[t,id,3])
                 push!(tempstockreg[Int(dataf[t+1,id,1])],dataf[t,id,4])
                 push!(tempdemandexpreg[Int(dataf[t+1,id,1])],dataf[t,id,5])
                 push!(tempsav,dataf[t+1,id,2])
                 push!(tempwork,dataf[t+1,id,3])
                 push!(tempstock,dataf[t+1,id,4])
                 push!(tempdemandexp,dataf[t+1,id,5])
                 end
             end
             for i=1:k1
                 for j=1:k2
                     ii = (i-1)*k2+j
                     f_av_savings_reg[t+1,tt,i,j] = mean(tempsavreg[ii])
                     f_av_workers_reg[t+1,tt,i,j] = mean(tempworkreg[ii])
                     f_av_stock_reg[t+1,tt,i,j] = mean(tempstockreg[ii])
                     f_av_demandexp_reg[t+1,tt,i,j] = mean(tempdemandexpreg[ii])
                 end
             end
             f_av_savings[t+1,tt] = mean(tempsav)
             f_var_savings[t+1,tt] = var(tempsav)
             f_av_workers[t+1,tt] = mean(tempwork)
             f_var_workers[t+1,tt] = var(tempwork)
             f_av_stock[t+1,tt] = mean(tempstock)
             f_av_demandexp[t+1,tt] = mean(tempdemandexp)
             f_var_demandexp[t+1,tt] = var(tempdemandexp)
         end
         tempworkreg = [[] for i = 1:k1, j = 1:k2]
         tempwork = []
         for id in firms
             if getagent(nhh+id).type-(nsec+1) == nsec+1
                 push!(tempworkreg[Int(dataf[t,id,1])],dataf[t,id,3])
                 push!(tempwork,dataf[t,id,3])
             end
         end
         for i=1:k1
             for j=1:k2
                 ii = (i-1)*k2+j
                 f_av_workers_reg[t,nsec+1,i,j] = mean(tempworkreg[ii])
             end
         end
         f_av_workers[t,nsec+1] = mean(tempwork)
         f_var_workers[t,nsec+1] = var(tempwork)

    conspercaph = 0
    for i=1:k1
        for j=1:k2
            castraj[t+1,i,j,1] =cas[i,j][1]
            castraj[t+1,i,j,2] =cas[i,j][2]
            for tt =1:nsec
                sumh = 0
                for jj = 1: datat
                    sumh += regcons[(t-1)*datat+jj][tt,i,j]
                end
                constraj[t+1,tt,i,j] = sumh
                conspercaph += sumh
            end
        end
    end
    conspercaptraj[t+1] = conspercaph / size(union(hh[1],hh[2]))[1]
    gdppercaptraj[t+1] = gdp / nhh # size(union(hh[1],hh[2]))[1]
    gdplosstraj[t+1] =  (gdppercapref - gdp / nhh)/ gdppercapref
    global totgdploss += gdplosstraj[t+1]
    for tt = 1:nsec
        bankrupttraj[t+1,tt] = bankruptcount[tt] / firmref[tt]
    end
    inactivefirmstraj[t+1] = inactivefirms
    bailouttraj[t+1] = sumbailouts
    totalaccountstraj[t+1] = totalaccounts
end

curinfregtrajy = [curinfregtraj[i][1,:,:] for i=1:size(curinfregtraj)[1]]
curinfregtrajo = [curinfregtraj[i][2,:,:] for i=1:size(curinfregtraj)[1]]
curinftrajy = [curinftraj[i][1] for i=1:size(curinftraj)[1]]
curinftrajo = [curinftraj[i][2] for i=1:size(curinftraj)[1]]
totcas = 0
for i=1:k1
    for j=1:k2
        global totcas += cas[i,j][1]+cas[i,j][2]
    end
end
totcas = totcas / nhh

if snapshot
    open("$(path_to_data)snapshot","w") do snapshot_file
        serialize(snapshot_file, (covidmodel.agents, lochh, locf, unemplist, shorttimelist, empcount, unempcount, shorttimecount, oldcount, unemp, firms, hh, tau, divperhh, weeklyconsumption))
    end
end

open("$(path_to_data)simdata.dat", "w") do outfile
    serialize(outfile, (datahhy, datahho, dataf, castraj, unempltraj, pubacctraj,constraj,unempsectraj,tautraj,conspercaptraj,R0counttraj,gdppercaptraj,R0count,f_av_savings,f_var_savings,f_av_workers,f_var_workers,f_av_stock,f_av_demandexp,f_var_demandexp,totinftraj,contact_count_traj,contact_work_traj,contact_social_traj, contact_shop_traj,inactivefirmstraj,bailouttraj,shorttimetraj,shorttimesectraj,totalaccountstraj,curinftrajy,curinftrajo,curinfregtrajy,curinfregtrajo,totfirmtraj,gdplosstraj,bankrupttraj))
end

include("covid_plot_reg.jl")

if virus
    println("R0: ", mean(R0count[1:200]))
    totgdploss = totgdploss * 365 / max(1, T - virustime)
    println("Av annual GDPloss of policy: ",totgdploss)
    println("Mortality rate ",totcas)
end


text:
using Test

@testset "ChromeProfileFormat" begin
end

ext: jl


text:
using Random, Distributions, Optim, Plots, Statistics

struct Theta
    z::Float64
    c::Float64
    y1::Float64
    y0::Float64
end

struct SampleTypes
    sample::Array{Theta}
end

function SampleTypes(n)
    z= rand(Normal(0,1), n)
    c = rand(Uniform(1, 10),n)
    θ = z .>0
    y1 = 5 .*(0.75.*θ .+ (1 .- θ))
    y0 = (1 .- θ).*10
    return [Theta(z[i], c[i], y1[i], y0[i]) for i in 1:n]
end

function generatePlots()
        n = 500
        data = SampleTypes(n)
        βstar = computeOptimalBeta(data, true)
        βces = computeOptimalBeta(data, false)
        plotNonStrategic(data, βces)
        savefig("nonstrategic.pdf")
        plotStrategic(data, βstar)
        savefig("strategic_optimal.pdf")
        plotStrategic(data, βces)
        savefig("strategic_ces.pdf")
end

function plotStrategic(data, β)
    x = [reportX(data[i], β) for i in 1:n]
    δ = [allocation(x[i], β) for i in 1:n]
    w = [rand(Bernoulli(δ[i]), 1)[1] for i in 1:n]
    y = [data[i].y1*w[i] + (1-w[i])*data[i].y0 for i in 1:n]
    println(mean(y))
    cate = [(data[i].y1 - data[i].y0) >0  for i in 1:n ]
    plot(x, zeros(length(x)) .+ rand(Normal(-0.2,0.05), length(x)), group=cate, label = ["Negative CATE" "Positive CATE"],
                    xrange = [-3, 3], yrange = [-0.3, 1.4], alpha=0.5, seriestype=:scatter)
    alloc(x) = allocation(x, β)
    plot!(alloc, -3, 3, lw = 2, label = "Assignment Prob.", linecolor=:black, xlabel="X", ylabel="Pr(W_i = 1| X__i = x)")
end

function plotNonStrategic(data, β)
    x = [data[i].z for i in 1:n]
    δ = [allocation(x[i], β) for i in 1:n]
    w = [rand(Bernoulli(δ[i]), 1)[1] for i in 1:n]
    y = [data[i].y1*w[i] + (1-w[i])*data[i].y0 for i in 1:n]
    println(mean(y))
    cate = [(data[i].y1 - data[i].y0) >0  for i in 1:n ]
    plot(x, zeros(length(x)) .+ rand(Normal(-0.2,0.05), length(x)), group=cate, label = ["Negative CATE" "Positive CATE"],
                    xrange = [-3, 3], yrange = [-0.3, 1.4], alpha=0.5, seriestype=:scatter)
    alloc(x) = allocation(x, β)
    plot!(alloc, -3, 3, lw = 2, label = "Assignment Prob.", linecolor=:black, xlabel="X", ylabel="Pr(W_i = 1| X_i = x)")
end

function run_Experiment(n = 2000, T=200)
    types = SampleTypes(10000)
    βopt = computeOptimalBeta(types, true)
    βces = computeOptimalBeta(types, false)

    bhat = 0.0
    Ropt = zeros(T); Rces = zeros(T); Rest = zeros(T)
    step = 0.3
    for t in 1:T
        agents = SampleTypes(n)
        h = 1.0
        ϵ =  h.*rand([-1,1], n)
        βp = bhat .+ ϵ
        data_est = generateObservations(agents, βp)
        Rest[t] = mean(data_est.y)
        Ropt[t] = computeProfit(agents, βopt); Rces[t] = computeProfit(agents, βces)
        γ̂ = linReg(data_est.y, ϵ)[2]
        bhat = bhat + step*γ̂
    end
    println("Estimated: ", bhat, "Optimal: ", βopt, "CES: ", βces)
    println("CES Regret: ", mean(Rces .- Ropt), "Estimated Regret: ", mean(Rest .- Ropt))
    plot(1:T, [Rces Ropt Rest], xlabel="Time Step", ylabel="Profit", labels= ["CES Rule" "Full Information Optimal" "Estimated Optimal"], legend=:bottomright)
end

function computeProfit(agents, β::Float64)
    data = generateObservations(agents, β .+ zeros(length(agents)))
    return mean(data.y)
end

function linReg(y, x)
    x = [ones(length(y)) x ]
    β = inv(x'*x)*x'*y
    return β
end

function generateObservations(data, β)
    n = length(data)
    x = [reportX(data[i], β[i]) for i in 1:n]
    δ = [allocation(x[i], β[i]) for i in 1:n]
    w = [rand(Bernoulli(δ[i]), 1)[1] for i in 1:n]
    y = [data[i].y1*w[i] + (1-w[i])*data[i].y0 for i in 1:n]
    return (y=y, w=w, x=x)
end

function computeOptimalBeta(sampleData, strategic::Bool)
    function objective(β)
        β = β[1]
        y = zeros(n)
        x = zeros(n)

        for i in 1:n
            agent = sampleData[i]
            if strategic
                x[i] = reportX(agent, β)
            else
                x[i] = agent.z
            end
            δ = allocation(x[i], β)
            w = rand(Bernoulli(δ), 1)[1]
            y[i] = agent.y1*w + (1-w)*agent.y0
        end
        return -mean(y)
    end
    if strategic
        res = Optim.optimize(β -> objective(β), -10, 10, GoldenSection(), abs_tol=0.00000001)
    else
        res = Optim.optimize(β -> objective(β), -100, 100, GoldenSection(), abs_tol=0.00000001)
    end
    println(res)
    #δ =allocation(x, Optim.minimizer(res))
    #println(δ[1:5])
    #println((x)[1:5])
    return Optim.minimizer(res)
end

function reportX(type::Theta, β)
    function objective(x)
        return -5*allocation(x, β) + type.c*(x - type.z)^2
    end
    res = Optim.optimize(x -> objective(x), -100, 100, GoldenSection())
    return Optim.minimizer(res)
end

function allocation(x, β)
    return 1 ./ (1 .+ exp.(-1 .* x.*β))
end

struct Sample
    sample::Array{Theta}
end
struct ObservedData
    x::Array{Float64}
    w::Array{Bool}
    y::Array{Float64}
end

max_stars_repo_name: evanmunro/personalized-policy

max_issues_repo_name: evanmunro/personalized-policy

max_forks_repo_name: evanmunro/personalized-policy


text:
≈(a::Number,b::Number) = isapprox(a, b, rtol=1e-3, atol=1e-8)

# with small calculations, such that the finite_sum_small branch is used
@testset "values of riemanntheta (small)" begin

    # using the Wolfram example
    # https://reference.wolfram.com/language/ref/SiegelTheta.html

    Ω = [ im -0.5 ; -0.5 im ]
    zs = [ComplexF64[0.5, 0.]]

    res = riemanntheta(zs, Ω)

    @test real(res[1]) ≈ 1.00748
    @test imag(res[1]) ≈ 0.

    zs = [ ComplexF64[x, 2x] for x in -1:0.1:1 ]
    res = riemanntheta(zs, Ω, eps=1e-3)

    @test maximum(real, res) ≈ 1.165
    @test minimum(real, res) ≈ 0.901
    @test maximum(imag, res) ≈ 0.
    @test minimum(imag, res) ≈ 0.

    # using example in abelfunctions doc
    # https://github.com/abelfunctions/abelfunctions/blob/master/doc/GettingStarted.md

    Ω = [ -1.309017+0.951057im -0.809017+0.587785im ;
          -0.809017+0.587785im -1.000000+1.175571im ]
    z = [0.5, 0.5im]
    res = riemanntheta([z], Ω, eps=1e-3)

    @test real(res[1]) ≈ 1.11415
    @test imag(res[1]) ≈ 0.8824
end

# with larger calculations, such that the finite_sum_large branch is used
@testset "values of riemanntheta (large)" begin

    # using the Wolfram example
    # https://reference.wolfram.com/language/ref/SiegelTheta.html

    Ω = [ im -0.5 ; -0.5 im ]
    zs = fill(ComplexF64[0.5, 0.], 10000)

    res = riemanntheta(zs, Ω)

    @test real(res[1]) ≈ 1.00748
    @test imag(res[1]) ≈ 0.

    zs = repeat( [ ComplexF64[x, 2x] for x in -1:0.01:1 ], outer=100)
    res = riemanntheta(zs, Ω, eps=1e-3)

    @test maximum(real, res) ≈ 1.165
    @test minimum(real, res) ≈ 0.901
    @test maximum(imag, res) ≈ 0.
    @test minimum(imag, res) ≈ 0.

    # using example in abelfunctions doc
    # https://github.com/abelfunctions/abelfunctions/blob/master/doc/GettingStarted.md

    Ω = [ -1.309017+0.951057im -0.809017+0.587785im ;
          -0.809017+0.587785im -1.000000+1.175571im ]
    z = [0.5, 0.5im]
    res = riemanntheta(fill(z,10000), Ω, eps=1e-3)

    @test real(res[1]) ≈ 1.11415
    @test imag(res[1]) ≈ 0.8824
end

################################################################################

g = 3
δ = 1e-8
Random.seed!(0)
tmp = (rand()*10 - 5.) * rand(g,g) ; Ω = Complex.(rand(g, g), tmp*tmp')
z₀ = (rand()*10 - 5.) * rand(ComplexF64, g)

@testset "oscillatory_part derivs are correct (small)" begin

    circvec = [Complex(1.,0.); zeros(ComplexF64, g-1)]

    # calculate function at slightly shifted z₀
    z = [[z₀] ; [z₀ + circshift(δ * circvec, i) for i in 0:g-1]]
    res = oscillatory_part(z, Ω)
    dres = [ (res[i] - res[1]) / δ for i in 2:g+1 ]

    # and compare to calculated derivates
    for i in 0:g-1
        derivs = [ circshift(circvec, i) ]
        res2 = RiemannTheta.oscillatory_part([z₀], Ω, derivs=derivs)
        # println(res2[1], dres[i+1])
        # res2[1]/δ , dres[i+1]/δ
        @test res2[1] ≈ dres[i+1]
    end

end

@testset "oscillatory_part derivs are correct (large)" begin

    circvec = [Complex(1.,0.); zeros(ComplexF64, g-1)]

    # calculate function at slightly shifted z₀
    z = [[z₀] ; [z₀ + circshift(δ * circvec, i) for i in 0:g-1]]
    res = RiemannTheta.oscillatory_part(z, Ω)
    dres = [ (res[i] - res[1]) / δ for i in 2:g+1 ]

    # and compare to calculated derivates
    for i in 0:g-1
        derivs = [ circshift(circvec, i) ]
        res2 = oscillatory_part(repeat([z₀], outer=10000), Ω, derivs=derivs)
        @test res2[1] ≈ dres[i+1]
    end

end

@testset "riemanntheta derivs are correct" begin

    circvec = [Complex(1.,0.); zeros(ComplexF64, g-1)]

    # calculate function at slightly shifted z₀
    z = [[z₀] ; [z₀ + circshift(δ * circvec, i) for i in 0:g-1]]
    res = riemanntheta(z, Ω)
    dres = [ (res[i] - res[1]) / δ for i in 2:g+1 ]

    # and compare to calculated derivates
    for i in 0:g-1
        derivs = [ circshift(circvec, i) ]
        res2 = riemanntheta([z₀], Ω, derivs=derivs)
        # println(res2[1], dres[i+1])
        # res2[1]/δ , dres[i+1]/δ
        @test res2[1] ≈ dres[i+1]
    end
end


text:
	u₀=0
	f(u,p,t) = -F_C*sign(u)+F_ext;
	g(u,p,t) = sqrt(2*D);
	# dt = 1//2^(4)
	tspan = (0.0,totalTime)
	prob = SDEProblem(f,g,u₀,tspan)
	sol = solve(prob,SRIW1(),dt=dt)

	# single sample path
	v = range(minimum(sol.u),maximum(sol.u),length=5000);
	p2 = histogram(sol.u, bins = 50, normalize = :pdf, xlabel = "v", ylabel = "P(v)");
	p1 = plot(sol.t,sol.u,xlabel = "t",ylabel= "v(t)");
	plot!(p2, v,Pst(F_C,F_ext,D,v))

	# ensemble prediction
	using DifferentialEquations.EnsembleAnalysis
	ensembleprob = EnsembleProblem(prob);
	solE = solve(ensembleprob,SRIW1(),EnsembleThreads();
		trajectories=1000);#,dt=totalTime/50) # EM(),EnsembleThreads() is optional
	summ = EnsembleSummary(solE, 0:0.5:totalTime)
	plot(summ)

	ensembleu = map(x-> x.u[end],solE);
	p3 = histogram(ensembleu, bins = 50, normalize = :pdf, xlabel = "v", ylabel = "P(v)");
	v = range(minimum(ensembleu),maximum(ensembleu),length=5000);
	plot!(p3, v,Pst(F_C,F_ext,D,v))
	plot(p1,p2,p3,layout = (1,3), size= (700,200))

	@time solve(ensembleprob,SRIW1(),EnsembleThreads(); trajectories=1000, dt=totalTime/50);
	@time solve(ensembleprob,SRIW1(); trajectories=1000, dt=totalTime/50);


text:
abstract type AbsAgent end

#returns an action
play(agent::AbsAgent, state::AbsState, env::AbsEnvironment; learn=true) = error("play is unimplemented")
observe(agent::AbsAgent, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false) = nothing

struct RandomAgent <: AbsAgent end

function play(agent::RandomAgent, state::AbsState, env::AbsEnvironment; learn=false)
	actionSet = getActions(state, env)
	return shuffle(actionSet)[1]
end

#Q-learning
mutable struct QLearner <: AbsAgent
	Qtable::Dict{Tuple{AbsState, AbsAction}, Float64}
	Ɣ::Float64
	α::Float64
	ε::Float64
	lastState
	lastAction

	function QLearner(env::AbsEnvironment; Ɣ=0.9, α=0.8, ε=0.05)
		#Initialize Qtable
		Qtable = Dict{Tuple{AbsState, AbsAction}, Float64}()
		for s in getAllStates(env)
			terminal = isTerminal(s, env)

			for a in getActions(s, env)
				Qtable[(s, a)] = terminal ? 0 : rand(-5:5)
			end
		end
		new(Qtable, Ɣ, α, ε, nothing, nothing)
	end
end

function maxQ(agent::QLearner, state::AbsState, env::AbsEnvironment)
	m = -1e32
	action = nothing
	for a in getActions(state, env)
		if agent.Qtable[(state, a)] > m
			action = a
			m = agent.Qtable[(state, a)]
		end
	end
	(action, m)
end

function play(agent::QLearner, state::AbsState, env::AbsEnvironment; learn=true)
	r = rand()
	action = nothing
	if r < agent.ε && learn
		actionSet = getActions(state, env)
		action = shuffle(actionSet)[1]
	else
		action, q = maxQ(agent, state, env)
	end

	agent.lastState = state
	agent.lastAction = action
	return action
end

function observe(agent::QLearner, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false)
	if learn
		action, q = maxQ(agent, state, env)
		agent.Qtable[(agent.lastState, agent.lastAction)] = agent.Qtable[(agent.lastState, agent.lastAction)] +
			agent.α * (reward + agent.Ɣ * q - agent.Qtable[(agent.lastState, agent.lastAction)])
		if terminal
			agent.lastState = nothing
			agent.lastAction = nothing
		end
	end
end

mutable struct SarsaLearner <: AbsAgent;
	qlearner::QLearner
	S
	A
	Sprime
	Aprime

	function SarsaLearner(env::AbsEnvironment; Ɣ=0.9, α=0.8, ε=0.05)
		ql = QLearner(env; Ɣ=Ɣ, α=α, ε=ε)
		new(ql, nothing, nothing, nothing, nothing)
	end
end

function play(agent::SarsaLearner, state::AbsState, env::AbsEnvironment; learn=true, observe=false)
	r = rand()
	action = nothing
	if r < agent.qlearner.ε && learn
		actionSet = getActions(state, env)
		action = shuffle(actionSet)[1]
	else
		action, q = maxQ(agent.qlearner, state, env)
	end

	if observe
		agent.Sprime = state
		agent.Aprime = action
	else
		if agent.S == nothing
			agent.S = state
			agent.A = action
		else
			action = learn ? agent.A : action
		end
	end
	return action
end

function observe(agent::SarsaLearner, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false)
	if learn
		action = play(agent, state, env; observe=true)
		agent.qlearner.Qtable[(agent.S, agent.A)] = agent.qlearner.Qtable[(agent.S,agent.A)] +
			agent.qlearner.α * (reward + agent.qlearner.Ɣ * agent.qlearner.Qtable[(agent.Sprime, agent.Aprime)] - agent.qlearner.Qtable[(agent.S,agent.A)])
		if terminal
			agent.S = nothing
			agent.A = nothing
			agent.Sprime = nothing
			agent.Aprime = nothing
		else
			agent.S = agent.Sprime
			agent.A = agent.Aprime
		end
	end
end

mutable struct SarsaLambdaLearner <: AbsAgent;
	qlearner::QLearner
	E::Dict{Tuple{AbsState, AbsAction}, Float64}
	λ
	S
	A
	Sprime
	Aprime

	function SarsaLambdaLearner(env::AbsEnvironment; λ=0.9, Ɣ=0.9, α=0.8, ε=0.05)
		ql = QLearner(env; Ɣ=Ɣ, α=α, ε=ε)
		E = Dict{Tuple{AbsState, AbsAction}, Float64}()
		for s in getAllStates(env)
			for a in getActions(s, env)
				E[(s, a)] = 0.0
			end
		end

		new(ql, E, λ, nothing, nothing, nothing, nothing)
	end
end

function play(agent::SarsaLambdaLearner, state::AbsState, env::AbsEnvironment; learn=true, observe=false)
	r = rand()
	action = nothing
	if r < agent.qlearner.ε && learn
		actionSet = getActions(state, env)
		action = shuffle(actionSet)[1]
	else
		action, q = maxQ(agent.qlearner, state, env)
	end

	if observe
		agent.Sprime = state
		agent.Aprime = action
	else
		if agent.S == nothing
			agent.S = state
			agent.A = action
		else
			action = learn ? agent.A : action
		end
	end
	return action
end

function observe(agent::SarsaLambdaLearner, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false)
	if learn
		action = play(agent, state, env; observe=true)
		delta = reward + agent.qlearner.Ɣ * agent.qlearner.Qtable[(agent.Sprime, agent.Aprime)] - agent.qlearner.Qtable[(agent.S, agent.A)]
		agent.E[(agent.S, agent.A)] = agent.E[(agent.S, agent.A)] + 1

		for s in getAllStates(env)
			for a in getActions(s, env)
				agent.qlearner.Qtable[(s, a)] = agent.qlearner.Qtable[(s, a)] + agent.qlearner.α * delta * agent.E[(s,a)]
				agent.E[(s, a)] = agent.qlearner.Ɣ * agent.λ * agent.E[(s,a)]
			end
		end

		if terminal
			agent.S = nothing
			agent.A = nothing
			agent.Sprime = nothing
			agent.Aprime = nothing

			for s in getAllStates(env)
				for a in getActions(s, env)
					agent.E[(s, a)] = 0.0
				end
			end

		else
			agent.S = agent.Sprime
			agent.A = agent.Aprime
		end
	end
end

mutable struct PolicyAgent <: AbsAgent
	policy::Policy
end

play(agent::PolicyAgent, state::AbsState, env::AbsEnvironment; learn=true) = agent.policy.mapping[state][1][1]

mutable struct MonteCarloAgent <: AbsAgent
	Qtable::Dict{AbsState, Dict{AbsAction, Float64}}
	policy::Policy
	returns::Dict{AbsState, Dict{AbsAction, Array{Float64, 1}}}
	rewards::Array{Float64, 1}
	path::Array{Tuple{AbsState, AbsAction}, 1}

	function MonteCarloAgent(env::AbsEnvironment)
		Qtable = Dict{AbsState, Dict{AbsAction, Float64}}()
		returns = Dict{AbsState, Dict{AbsAction, Array{Float64, 1}}}()
		policy = Policy()
		rewards = Array{Float64, 1}()
		path = Array{Tuple{AbsState, AbsAction}, 1}()

		for s in getAllStates(env)
			Qtable[s] = Dict{AbsAction, Float64}()
			returns[s] = Dict{AbsAction, Array{Float64, 1}}()

			terminal = isTerminal(s, env)
			as = getActions(s, env)
			a = shuffle(collect(as))[1]
			policy.mapping[s] = [(a, 1.0)]

			for a in as
				returns[s][a] = Array{Float64, 1}()
				Qtable[s][a] = rand(-5:5)
			end
		end

		new(Qtable, policy, returns, rewards, path)
	end
end

function play(agent::MonteCarloAgent, state::AbsState, env::AbsEnvironment; learn=true)
	action = nothing
	if length(agent.path) == 0 && learn
		as = getActions(state, env)
		action = shuffle(collect(as))[1]
	else
		action = agent.policy.mapping[state][1][1]
	end
	push!(agent.path, (state, action))
	return action
end

function observe(agent::MonteCarloAgent, state::AbsState, reward::Float64, env::AbsEnvironment; learn=true, terminal=false)
	if learn
		push!(agent.rewards, reward)

		if terminal
			visited = Set{Tuple{AbsState, AbsAction}}()
			g = sum(agent.rewards)
			for t in agent.path
				if ! (t in visited)
					push!(agent.returns[t[1]][t[2]], g)
					agent.Qtable[t[1]][t[2]] = mean(agent.returns[t[1]][t[2]])
				end
			end

			#update the policy
			for t in agent.path
				m = -1e32
				action = nothing
				for a in getActions(t[1], env)
					if agent.Qtable[t[1]][a] > m
						m = agent.Qtable[t[1]][a]
						action = a
					end
				end
				agent.policy.mapping[t[1]] = [(action, 1.0)]
			end

			agent.rewards = Array{Float64, 1}()
			agent.path = Array{Tuple{AbsState, AbsAction}, 1}()
		end
	end
end

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
#=

# 420: Discontinuous Quantities
 ([source code](SOURCE_URL))

Test  jumping species and quantity handling

=#
module Example420_DiscontinuousQuantities
using Printf
using VoronoiFVM
using SparseArrays
using ExtendableGrids
using GridVisualize
using LinearAlgebra

function main(;N=5, Plotter=nothing,unknown_storage=:sparse)
    XX=collect(0:0.1:1)
    xcoord=XX
    for i=1:N-1
	xcoord=glue(xcoord,XX.+i)
    end
    grid2=simplexgrid(xcoord)
    for i=1:N
	cellmask!(grid2,[i-1],[i],i)
    end
    for i=1:N-1
	bfacemask!(grid2,[i],[i],i+2)
    end

    params=zeros(2,num_cellregions(grid2))
    for i=1:num_cellregions(grid2)
        params[1,i]=i
        params[2,i]=10*i
    end

    system=VoronoiFVM.System(grid2,unknown_storage=unknown_storage)

    ## First, we introduce a continuous quantity which we name "cspec". Note that the "species number" can be assigned automatically if not given explicitely.
    cspec=ContinuousQuantity(system,1:N;ispec=1,id=1)

    ## A discontinuous quantity can be introduced as well. by default, each reagion gets a new species number. This can be overwritten by the user.
    dspec=DiscontinuousQuantity(system,1:N; regionspec=[2+i%2 for i=1:N],id=2)

    # check 1D array acces with quantities
    carrierList = [cspec dspec]
    numberCarriers = length(carrierList)

    params2=zeros(1, numberCarriers)

    for icc ∈ carrierList
        params2[icc] = 2
    end

    for i=1:numberCarriers
        @assert params2[i] == 2
    end

    # check 2D array acces with quantities
    for i=1:num_cellregions(grid2)
        @assert params[cspec,i] == i
        @assert params[dspec,i] == 10*i
    end

    for i=1:num_cellregions(grid2)
        params[cspec,i] = -i
        params[dspec,i] = -10*i
    end

    for i=1:num_cellregions(grid2)
        @assert params[1,i] == -i
        @assert params[2,i] == -10*i
    end

    ##For both quantities, we define simple diffusion fluxes:

    function flux(f,u,edge)
	f[dspec]=u[dspec,1]-u[dspec,2]
	f[cspec]=u[cspec,1]-u[cspec,2]
    end

    d1=1
    q1=0.2

    function breaction(f,u,bnode)

        # left outer boundary value for dspec
        if bnode.region == 1
            f[dspec] = u[dspec] + 0.5
        end

        ## Define a thin layer inteface condition for `dspec` and an interface source for `cspec`.
        if bnode.region>2
            react=(u[dspec,1]-u[dspec,2])/d1
            f[dspec,1]= react
            f[dspec,2]= -react
            f[cspec]=-q1*u[cspec]
        end
    end

    physics!(system,VoronoiFVM.Physics(
        flux=flux,
        breaction=breaction
    ))

    ## Set boundary conditions
    boundary_dirichlet!(system,dspec,2,0.1)
    boundary_dirichlet!(system,cspec,1,0.1)
    boundary_dirichlet!(system,cspec,2,1.0)
    subgrids=VoronoiFVM.subgrids(dspec,system)

    U=solve(unknowns(system,inival=0),system)

    dvws=VoronoiFVM.views(U,dspec,subgrids,system)
    cvws=VoronoiFVM.views(U,cspec,subgrids,system)
    vis=GridVisualizer(resolution=(600,300),Plotter=Plotter)
    for i=1:length(dvws)
	scalarplot!(vis,subgrids[i],dvws[i],flimits=(-0.5,1.5),clear=false,color=:red)
	scalarplot!(vis,subgrids[i],cvws[i],flimits=(-0.5,1.5),clear=false,color=:green)
    end
    reveal(vis)
    norm(system,U,2)
end

function test()
    testval=7.812799873197911
    main(unknown_storage=:sparse) ≈ testval &&
        main(unknown_storage=:dense) ≈ testval
end

end


text:
module execute1d
using LPDM
using POMDPs, Parameters, StaticArrays, D3Trees, Distributions, SparseArrays
using Combinatorics
using StatsFuns
using Printf
using Dates

include("LightDarkLPDM.jl")
using LightDarkPOMDPs

include("1D/LPDMBounds1d.jl")

struct LPDMTest
    solver_mode::Symbol
    action_mode::Symbol
    obs_mode::Symbol
    reward_mode::Symbol
end

struct LPDMScenario
    s0::LD1State
end

# reward function options - :quadratic or :fixed
function batch_execute(;n::Int64=1, debug::Int64=1, reward_mode=:quadratic)
    test=Array{LPDMTest}(undef,0)

    push!(test, LPDMTest(:lpdm, :standard, :discrete, reward_mode))
    push!(test, LPDMTest(:lpdm, :standard, :continuous, reward_mode))
    # push!(test, LPDMTest(:lpdm, :adapt, reward_mode)) # simulated annealing
    # push!(test, LPDMTest(:despot, :small, reward_mode))
    # push!(test, LPDMTest(:despot, :large, reward_mode))
    # push!(test, LPDMTest(:lpdm_bv, :bv, reward_mode)) # blind value

    scen=Array{LPDMScenario}(undef,0)
    push!(scen, LPDMScenario(LD1State(-2*π)))
    # push!(scen, LPDMScenario(LD1State(π/2)))
    # push!(scen, LPDMScenario(LD1State(3/2*π)))
    # push!(scen, LPDMScenario(LD1State(2*π)))

    # Dummy execution, just to make sure all the code is compiled and loaded,
    # to improve uniformity of subsequent executions.
    # execute(solv_mode = :lpdm, action_space_type = :adapt, n_sims = 1, s0 = LD1State(π), output = 0)

    # General solver parameters
    steps::Int64                = -1
    time_per_move::Float64      = 1.0
    search_depth::Int64         = 30
    n_particles::Int64          = 50
    max_trials::Int64           = -1

    f = open("results_" * Dates.format(now(),"yyyy-mm-dd_HH_MM") * ".txt", "w")

    Printf.@printf(f,"GENERAL SOLVER PARAMETERS\n")
    Printf.@printf(f,"\treward function:\t\t\t%s\n", reward_func)
    Printf.@printf(f,"\tsteps:\t\t\t%d\n", steps)
    Printf.@printf(f,"\ttime per move:\t\t\t%f\n", time_per_move)
    Printf.@printf(f,"\tsearch depth:\t\t\t%d\n", search_depth)
    Printf.@printf(f,"\tN particles:\t\t\t%d\n", n_particles)
    Printf.@printf(f,"\tmax trials:\t\t\t%d\n\n", max_trials)

    for i in 1:length(scen)
        # write(f,"SCENARIO $i, s0 = $(scen[i].s0)\n\n")
        if debug >= 0
            println("")
            println("SCENARIO $i, s0 = $(scen[i].s0)")
            println("------------------------")
        end

        Printf.@printf(f,"SCENARIO %d, s0 = %f\n", i, scen[i].s0)
        Printf.@printf(f,"=====================================================================\n")
        # Printf.@printf(f,"mode\t\tact. space\t\tsteps(std)\t\treward(std)\n")
        Printf.@printf(f,"SOLVER\t\tACT. MODE\t\tOBS. MODE\t\tSTEPS (STD)\t\t\tREWARD (STD)\n")
        Printf.@printf(f,"=====================================================================\n")
        for t in test
            if debug >= 0
                println("mode: $(t.solver_mode), actions: $(t.action_mode), observations: $(t.obs_mode), rewards: $(t.reward_mode)")
            end
            steps, steps_std, reward, reward_std =
                        execute(solver_mode       = t.solver_mode,
                                action_mode       = t.action_mode,
                                obs_mode          = t.obs_mode,
                                reward_mode       = t.reward_mode,
                                n_sims            = n,
                                s0                = scen[i].s0,
                                output            = debug)

            Printf.@printf(f,"%s\t\t%s\t\t%s\t\t\t%05.2f (%06.2f)\t\t%06.2f (%06.2f)\n",
                            string(t.solver_mode), string(t.action_mode), string(t.obs_mode), steps, steps_std, reward, reward_std)
        end
        Printf.@printf(f,"=====================================================================\n")
        Printf.@printf(f,"%d tests per scenario\n\n", n)
    end
    close(f)
end

function execute(;vis::Vector{Int64}=Int64[],
                solver_mode::Symbol=:despot,
                action_mode::Symbol=:standard,
                obs_mode::Symbol=:discrete,
                reward_mode = :quadratic,
                n_sims::Int64=1,
                steps::Int64=-1,
                s0::LD1State=LD1State(1.9),
                output::Int64=1)#n_sims::Int64 = 100)

    if solver_mode == :despot
        p = LightDark1DDespot(action_mode = action_mode,
                            obs_mode = obs_mode,
                            reward_mode = reward_mode)
    elseif solver_mode == :lpdm
        p = LightDark1DLpdm(action_mode = action_mode,
                            obs_mode = obs_mode,
                            reward_mode = reward_mode)
    end

    sim_rewards = Vector{Float64}(undef,n_sims)
    sim_steps   = Vector{Int64}(undef,n_sims)

    for sim in 1:n_sims
        # world_seed  ::UInt32   = convert(UInt32, 42)

        world_rng = RNGVector(1, UInt32(sim))
        LPDM.set!(world_rng, 1)

        # NOTE: restrict initial state to positive numbers only, for now
        # s::LD1State                  = LD1State(π);    # initial state
        step_rewards::Array{Float64}     = Vector{Float64}(undef,0)
        # println("$(supertype(LightDark1DDespot)))")

        # NOTE: Original version
        # solver = LPDM.LPDMSolver{LD1State, LD1Action, LD1Obs, LDBounds1d{LD1State, LD1Action, LD1Obs}, RNGVector}(
        #                                                                     # rng = sim_rng,
        #                                                                     debug = output,
        #                                                                     time_per_move = 1.0,  #sec
        #                                                                     # time_per_move = 1.0,  #sec
        #                                                                     sim_len = -1,
        #                                                                     search_depth = 50,
        #                                                                     n_particles = 20,
        #                                                                     # seed = UInt32(5),
        #                                                                     seed = UInt32(2*sim+1),
        #                                                                     # max_trials = 10)
        #                                                                     max_trials = -1,
        #                                                                     mode = solv_mode)

        # NOTE: Test version
        solver = LPDM.LPDMSolver{LD1State, LD1Action, LD1Obs, LDBounds1d{LD1State, LD1Action, LD1Obs}, RNGVector}(
                                                                            # rng = sim_rng,
                                                                            debug = output,
                                                                            time_per_move = -1.0,  #sec
                                                                            # time_per_move = 1.0,  #sec
                                                                            sim_len = steps,
                                                                            search_depth = 50,
                                                                            n_particles = 20,
                                                                            # seed = UInt32(2),
                                                                            seed = UInt32(2*sim+1),
                                                                            # max_trials = 10)
                                                                            max_trials = 1,
                                                                            solver_mode = solver_mode,
                                                                            action_mode = action_mode,
                                                                            obs_mode    = obs_mode)

    #---------------------------------------------------------------------------------
        # Belief
        bu = LPDMBeliefUpdater(p,
                               n_particles = solver.config.n_particles,
                               seed = UInt32(3*sim+1));  # initialize belief updater
        initial_states = state_distribution(p, s0, solver.config, world_rng)     # create initial  distribution
        current_belief = LPDM.create_belief(bu)                       # allocate an updated belief object
        LPDM.initialize_belief(bu, initial_states, current_belief)    # initialize belief
        # show(current_belief)
        updated_belief = LPDM.create_belief(bu)
    #---------------------------------------------------------------------------------

        policy::LPDMPolicy = POMDPs.solve(solver, p)
        s = s0
        step::Int64 = 0
        r::Float64 = 0.0

        # output >= 1 && println("=============== SIMULATION # $sim ================")
        if output >= 0
            println("")
            println("*** SIM $sim (solver mode: $solver_mode, action mode: $action_mode, obs mode: $obs_mode, s0: $s0)***")
            println("")
        end

        val, run_time, bytes, gctime, memallocs =
        @timed while !isterminal(p, s) && (solver.config.sim_len == -1 || step < solver.config.sim_len)
            step += 1

            if output >= 1
                println("")
                println("=============== Step $step ================")
                show(current_belief)
            end

            a = POMDPs.action(policy, current_belief)
            if output >= 1
                println("s: $s")
                println("a: $a")
            end

            s, o, r = POMDPs.generate_sor(p, s, a, world_rng)
            push!(step_rewards, r)
            if output >= 1
                println("s': $s")
                println("o': $o")
                println("r: $r")
                println("=======================================")
            end

            # error("$(@which(POMDPs.update(bu, current_belief, a, o, updated_belief)))")
            # update belief
            POMDPs.update(bu, current_belief, a, o, updated_belief)
            current_belief = deepcopy(updated_belief)
            # show(updated_belief) #NOTE: don't show for now

            if LPDM.isterminal(p, current_belief.particles)
                if output >= 1
                    println("Terminal belief. Execution completed.")
                    show(current_belief)
                end
                break
            end

            if step ∈ vis
                t = LPDM.d3tree(solver,
                                detect_repeat=false,
                                title="Step $step",
                                init_expand=10)
                # # show(t)
                inchrome(t)
                # blink(t)
            end
            if output >= 1
                solver_mode == :lpdm && println("root actions: $(solver.root.action_space)")
            end
        end

        sim_steps[sim] = step
        sim_rewards[sim] = sum(step_rewards)
        # # Compute discounted reward
        # discounted_reward::Float64 = 0.0
        # multiplier::Float64 = 1.0
        # for r in rewards
        #     discounted_reward += multiplier * r
        #     multiplier *= p.discount
        # end
        # println("Discounted reward: $discounted_reward")
    end
    # println("SIM_STEPS: $sim_steps, SIM_REWARDS: $sim_rewards")
    return mean(sim_steps), std(sim_steps), mean(sim_rewards), std(sim_rewards)
    # return t
end

function test_move()
    p = LightDarkPOMDPs.LightDark1DDespot()

    coordinates = Vector{Tuple{Float64,Float64}}()
    push!(coordinates,(5.0,0.0))
    push!(coordinates,(-5.0,0.0))
    push!(coordinates,(3.14,4.50))

    for c in coordinates
        r,a = move(p,c[1],c[2])
        println("x1=$(c[1]), x2=$(c[2]), r=$r, a=$a")
    end
end

end #module

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
function Investment_OPF_stage1(optimizer,set_opt_thermalgenerators,set_opt_winds,set_thermalgenerators,set_winds,set_demands,set_nodes,set_nodes_ref,set_nodes_noref,set_scenarios,set_times,P,V,max_demand,R,p_D,D,γ,Τ,wind,wind_opt,Ns_H,links,F_max_dict,B_dict,MapG,MapG_opt,MapD,MapW,MapW_opt,tech_thermal,tech_thermal_opt,tech_wind,tech_wind_opt,capacity_per_unit,varcost,invcost,maxBuilds,ownership,capacity_existingunits,fixedcost,EmissionsRate,HeatRate,fuelprice,life,CRF_thermal_opt,CRF_wind_opt,set_K,Set_Dual_constraint3_st2,Set_Dual_constraint4_st2,Set_Dual_constraint5_st2,Set_Dual_constraint6_st2,Set_Dual_constraint8_st2,Set_Dual_constraint10_st2)

m=Model(optimizer)

#=
@variable(m, f[link in links,t in set_times,s in set_scenarios] )
@variable(m, θ[n in set_nodes,t in set_times,s in set_scenarios])
@variable(m, p_G_e[e in set_thermalgenerators,t in set_times,s in set_scenarios]>= 0)
@variable(m, p_G_e_opt[e in set_opt_thermalgenerators,t in set_times,s in set_scenarios]>= 0)
@variable(m, p_G_w[w in set_winds,t in set_times,s in set_scenarios]>= 0)
@variable(m, p_G_w_opt[w in set_opt_winds,t in set_times,s in set_scenarios]>= 0)
@variable(m, r_d[d in set_demands,t in set_times,s in set_scenarios]>= 0)
@variable(m, r_w[w in set_winds,t in set_times,s in set_scenarios]>= 0)
@variable(m, r_w_opt[w in set_opt_winds,t in set_times,s in set_scenarios]>= 0)
=#

@variable(m, x_e[e in set_opt_thermalgenerators]>= 0)
@variable(m, x_w[w in set_opt_winds]>= 0)
#Auxiliary variable Benders
@variable(m, α >= 0)
#Auxiliary variable to determine the cost per scenario
#@variable(m, problem1_cost_scenario[s in set_scenarios])
@variable(m, investment_cost)

@objective(m, Min, sum(tech_thermal_opt[e,invcost]*CRF_thermal_opt[e]*x_e[e]  for e in set_opt_thermalgenerators)
                  +sum(tech_wind_opt[w,invcost]*CRF_wind_opt[w]*x_w[w]  for w in set_opt_winds)
                  +sum(tech_thermal[e,fixedcost]*(tech_thermal[e,capacity_existingunits])  for e in set_thermalgenerators)
                  +sum(tech_wind[w,fixedcost]*(tech_wind[w,capacity_existingunits])  for w in set_winds)
                  +sum(tech_thermal_opt[e,fixedcost]*(x_e[e])  for e in set_opt_thermalgenerators)
                  +sum(tech_wind_opt[w,fixedcost]*(x_w[w])  for w in set_opt_winds)
                  +α)

@constraint(m, constraint1,  sum(tech_thermal[e,capacity_existingunits] for e in set_thermalgenerators)
                            +sum(tech_wind[w,capacity_existingunits] for w in set_winds)
                            +sum(x_w[w] for w in set_opt_winds)
                            +sum(x_e[e] for e in set_opt_thermalgenerators) >= D*(1+R))

@constraint(m, constraint2[w in set_opt_winds], x_w[w] <= tech_wind_opt[w,maxBuilds])

@constraint(m, constraint3,  α>=α_down)

@constraint(m, constraint4[K in set_K],
                      α>=sum( (sum(sum(tech_thermal[e,capacity_existingunits]*Set_Dual_constraint3_st2[e,t,s,K] for e in set_thermalgenerators) for t in set_times)
                             +sum(sum(x_e[e]*Set_Dual_constraint4_st2[e,t,s,K] for e in set_opt_thermalgenerators) for t in set_times)
                             +sum(sum(tech_wind[w,capacity_existingunits]*wind[w,t]*Set_Dual_constraint5_st2[w,t,s,K] for w in set_winds) for t in set_times)
                             +sum(sum(x_w[w]*wind_opt[w,t]*Set_Dual_constraint6_st2[w,t,s,K] for w in set_opt_winds) for t in set_times)
                             +sum(sum(F_max_dict[Z[2][j]]*Set_Dual_constraint8_st2[Z[1][j],t,s,K] for j in set_number_links) for t in set_times)
                             +sum(sum(sum(p_D[d,t]*max_demand for d in set_demands if n == MapD[d][2])*Set_Dual_constraint10_st2[n,t,s,K] for n in set_nodes) for t in set_times)
                                 )*γ[s] for s in set_scenarios  )  )

#@constraint(m, constraint3[e in set_thermalgenerators,t in set_times,s in set_scenarios], p_G_e[e,t,s] <= tech_thermal[e,capacity_existingunits])

#@constraint(m, constraint4[e in set_opt_thermalgenerators,t in set_times,s in set_scenarios], p_G_e_opt[e,t,s] <= x_e[e])

#@constraint(m, constraint5[w in set_winds,t in set_times,s in set_scenarios], p_G_w[w,t,s] +r_w[w,t,s] == tech_wind[w,capacity_existingunits]*wind[w,t])

#@constraint(m, constraint6[w in set_opt_winds,t in set_times,s in set_scenarios], p_G_w_opt[w,t,s] +r_w_opt[w,t,s] == x_w[w]*wind_opt[w,t])

#@constraint(m, constraint7[j in links,t in set_times,s in set_scenarios], f[j,t,s] == B_dict[j]*(θ[j[1],t,s]-θ[j[2],t,s]))

#@constraint(m, constraint8[j in links,t in set_times,s in set_scenarios], f[j,t,s] <= F_max_dict[j])

#@constraint(m, constraint9[t in set_times,s in set_scenarios], θ[set_nodes_ref,t,s] == 0 )

#=
@constraint(m, constraint10[n in set_nodes, t in set_times,s in set_scenarios],
              +sum(p_G_e[e,t,s] for e in set_thermalgenerators if n == MapG[e][2])
              +sum(p_G_e_opt[e,t,s] for e in set_opt_thermalgenerators if n == MapG_opt[e][2])
              +sum(p_G_w[w,t,s] for w in set_winds if n == MapW[w][2])
              +sum(p_G_w_opt[w,t,s] for w in set_opt_winds if n == MapW_opt[w][2])
              +sum(r_d[d,t,s] for d in set_demands if n == MapD[d][2])
              -sum(p_D[d,t]*max_demand for d in set_demands if n == MapD[d][2])
              -sum(f[j,t,s] for j in links if n == j[1])== 0)
=#

#Auxiliary Constraints= Cost per scenario
#=
@constraint(m, constraint11[s in set_scenarios], problem1_cost_scenario[s] == sum((
                  +sum(tech_thermal[e,varcost]*p_G_e[e,t,s] for e in set_thermalgenerators)
                  +sum(tech_thermal_opt[e,varcost]*p_G_e_opt[e,t,s]  for e in set_opt_thermalgenerators)
                  +sum(tech_wind[w,varcost]*p_G_w[w,t,s]  for w in set_winds)
                  +sum(tech_wind_opt[w,varcost]*p_G_w_opt[w,t,s]  for w in set_opt_winds)
                  +sum(V*r_d[d,t,s]  for d in set_demands)
                  +sum(P*r_w[w,t,s]  for w in set_winds)
                  +sum(P*r_w_opt[w,t,s]  for w in set_opt_winds)
                  +sum(Τ[s]*tech_thermal[e,EmissionsRate]*tech_thermal[e,HeatRate]*p_G_e[e,t,s] for e in set_thermalgenerators)
                  +sum(Τ[s]*tech_thermal_opt[e,EmissionsRate]*tech_thermal_opt[e,HeatRate]*p_G_e_opt[e,t,s] for e in set_opt_thermalgenerators))*Ns_H[t]  for t in set_times)
                  +sum(tech_thermal_opt[e,invcost]*CRF_thermal_opt[e]*x_e[e]  for e in set_opt_thermalgenerators)
                  +sum(tech_wind_opt[w,invcost]*CRF_wind_opt[w]*x_w[w]  for w in set_opt_winds)
                  +sum(tech_thermal[e,fixedcost]*(tech_thermal[e,capacity_existingunits])  for e in set_thermalgenerators)
                  +sum(tech_wind[w,fixedcost]*(tech_wind[w,capacity_existingunits])  for w in set_winds)
                  +sum(tech_thermal_opt[e,fixedcost]*(x_e[e])  for e in set_opt_thermalgenerators)
                  +sum(tech_wind_opt[w,fixedcost]*(x_w[w])  for w in set_opt_winds))
=#
#Auxiliary Constraints= Investment Cost

@constraint(m, constraint12, investment_cost == sum(tech_thermal_opt[e,invcost]*CRF_thermal_opt[e]*x_e[e]  for e in set_opt_thermalgenerators)
+sum(tech_wind_opt[w,invcost]*CRF_wind_opt[w]*x_w[w]  for w in set_opt_winds)
+sum(tech_thermal[e,fixedcost]*(tech_thermal[e,capacity_existingunits])  for e in set_thermalgenerators)
+sum(tech_wind[w,fixedcost]*(tech_wind[w,capacity_existingunits])  for w in set_winds)
+sum(tech_thermal_opt[e,fixedcost]*(x_e[e])  for e in set_opt_thermalgenerators)
+sum(tech_wind_opt[w,fixedcost]*(x_w[w])  for w in set_opt_winds))

@time optimize!(m)

#cd("C:\\Users\\braya\\Desktop\\Doctorado\\Doctorado EEUU\\Semesters\\Spring Semester 2021\\DTU Course\\Project\\Step 4\\Julia")
#Print output
println("Number of variables and Constraints:")
display(m)

status = termination_status(m)
println("The solution status is: $status")

syscost_det=objective_value(m)

println("Total Cost:",syscost_det)

investment_cost_value=JuMP.value.(investment_cost)

x_w_value=zeros(length(set_opt_winds))
for w in set_opt_winds
println("Investment decision for candidate wind unit $w: ",JuMP.value.(x_w[w]))
x_w_value[w]=JuMP.value.(x_w[w])
end

x_e_value=zeros(length(set_opt_thermalgenerators))
for e in set_opt_thermalgenerators
println("Investment decision for candidate thermal unit $e: ",JuMP.value.(x_e[e]))
x_e_value[e]=JuMP.value.(x_e[e])
end

return (syscost_det,x_w_value,x_e_value,investment_cost_value)
end

#=
p_G_e_value=zeros(length(set_thermalgenerators),length(set_times),length(set_scenarios))
for t in set_times, s in set_scenarios, e in set_thermalgenerators
#println("Production level of existing thermal generator $e in time $t under scenario $s:",JuMP.value.(p_G_e[e,t,s]))
  p_G_e_value[e,t,s]=JuMP.value.(p_G_e[e,t,s])
end

p_G_e_opt_value=zeros(length(set_opt_thermalgenerators),length(set_times),length(set_scenarios))
for t in set_times, s in set_scenarios, e in set_opt_thermalgenerators
#println("Production level of candidate thermal generator $e in time $t under scenario $s:",JuMP.value.(p_G_e_opt[e,t,s]))
  p_G_e_opt_value[e,t,s]=JuMP.value.(p_G_e_opt[e,t,s])
end

p_G_w_value=zeros(length(set_winds),length(set_times),length(set_scenarios))
for t in set_times, s in set_scenarios, w in set_winds
#println("Production level of wind unit $w in time $t under scenario $s:",JuMP.value.(p_G_w[w,t,s]))
  p_G_w_value[w,t,s]=JuMP.value.(p_G_w[w,t,s])
end

p_G_w_opt_value=zeros(length(set_opt_winds),length(set_times),length(set_scenarios))
for t in set_times, s in set_scenarios, w in set_opt_winds
#println("Production level of candidare wind unit $w in time $t under scenario $s: ",JuMP.value.(p_G_w_opt[w,t,s]))
  p_G_w_opt_value[w,t,s]=JuMP.value.(p_G_w_opt[w,t,s])
end

p_D_value=zeros(length(set_demands),length(set_scenarios))
for d in set_demands, s in set_scenarios
  p_D_value[d,s]=p_D[d,s]*max_demand
end
#println("Comsuption level of demand d: ",p_D_value)

r_d_value=zeros(length(set_demands),length(set_times),length(set_scenarios))
for t in set_times, s in set_scenarios, d in set_demands
#println("Unserved energy of demand $d in time $t under scenario $s:",JuMP.value.(r_d[d,t,s]))
  r_d_value[d,t,s]=JuMP.value.(r_d[d,t,s])
end

r_w_value=zeros(length(set_winds),length(set_times),length(set_scenarios))

for t in set_times, s in set_scenarios, w in set_winds
#println("Curtailment of wind unit $w in time $t under scenario $s: ",JuMP.value.(r_w[w,t,s]))
  r_w_value[w,t,s]=JuMP.value.(r_w[w,t,s])
end

Dual_constraint10=zeros(length(set_nodes),length(set_times),length(set_scenarios))

for   t in set_times, s in set_scenarios,n in set_nodes
   Dual_constraint10[n,t,s]= JuMP.dual(constraint10[n,t,s])/(Ns_H[t]*γ[s])
#println("Electricity Price in node $n in time $t under scenario $s: ", Dual_constraint10[n,t,s])
end

df1 = DataFrame( Dual_constraint10)
#CSV.write("Electrity_Price.csv", df1)

#for s in set_scenarios, t in set_times, j in links
#println("Power Flow lines $j in time $t under scenario $s: ", JuMP.value.(f[j,t,s]))
#end

f_value=zeros(n_link,length(set_times),length(set_scenarios))
global i=1
for j in links
  for t in set_times
  for s in set_scenarios
global f_value[i,t,s]= JuMP.value.(f[j,t,s])
  end
end
global i=i+1
end

=#

#=
print("")

for n in set_nodes, s in set_scenarios, t in set_times
   println("Electricity Price in node $n in time $t under scenario $s: ", JuMP.dual(constraint10[n,t,s]))
end
=#


text:
# In this file, kernel functions and their gradients are defined

# W(r_ij, h_k)
function W(i, j, k=i)
  h = sph[k].hsml
  invh = 1 / h
  r = dist(i, j) * invh

  if kernel == "CubicSpline"
    norm = 4 / 3 * invh
    if r > 1
      return 0
    elseif r > 0.5
      return norm * 2 * (1 - r)^3
    else
      return norm * (1 - 6 * r^2 + 6 * r^3)
    end
  elseif kernel == "WendlandC2"
    norm = 5 / 4 * invh
    if r > 1
      return 0
    else
      return norm * (1 - r)^3 * (1 + 3 * r)
    end
  else
    @assert 0
  end
end

# nabla_i W(r_ij, h_k)
function gradW(i, j, k=i)
  h = sph[k].hsml
  invh = 1 / h
  r = dist(i, j) * invh

  if gradient == "Standard"
    if kernel == "CubicSpline"
      norm = 4 / 3 * invh
      if r > 1
        gradW = 0
      elseif r > 0.5
        gradW = norm * 6 * (1 - r)^2 * (-invh)
      else
        gradW = norm * (-12 * r * invh + 18 * r^2 * invh)
      end
    elseif kernel == "WendlandC2"
      norm = 5 / 4 * invh
      if r > 1
        gradW = 0
      else
        gradW = norm * (-3 * (1 - r)^2 * (1 + 3 * r) * invh + 3 * (1 - r)^3 * invh)
      end
    else
      @assert 0
    end
    gradW = eij(i, j) * gradW
  elseif gradient == "IntegralApproach"
    gradW = - eij(i, j) * dist(i, j) * W(i, j, k) / sph[k].c11
  end
  return gradW
end

# dW(r_ij, h_i)/dh_i
function dWdh(i, j)
  h = sph[i].hsml
  invh = 1 / h
  r = dist(i, j) * invh

    if kernel == "CubicSpline"
      if r > 1
        return 0
      elseif r > 0.5
        return -8 / 3 * invh * (1 - r)^2 * (1 + r)
      else
        return -4 / 3 * invh * (1 - 18 * r^2 + 24 * r^3)
      end
    elseif kernel == "WendlandC2"
      if r > 1
        return 0
      else
        return 5 / 4 * invh^2 * (1 - r)^2 * (3 * r - 1) * (5 * r + 1)
      end
    else
      @assert 0
    end
end

# averaged for symmetry against exchange of i and j
function W_sym(i, j)
  return 0.5 * (W(i, j) + W(j, i))
end

# nabla_i (W(r_ij, h_i) + W(r_ij, h_j))/2
function gradW_sym(i, j)
  return 0.5 * (gradW(i, j) + gradW(i, j, j))
end

# This correction matrix accounts for local particle distribution.
# In this one-dimentional code, this is not a matrix but a scaler
function correction_matrix(i, j)
  return Z(j)/y(j) * dist(i, j)^2 * W(i, j)
end


text:
################################################################################
###  Introduction to Scientific Programming and Machine Learning with Julia  ###
###                                                                          ###
### Run each script on a new clean Julia session                             ###
### GitHub: https://github.com/sylvaticus/IntroSPMLJuliaCourse               ###
### Licence (apply to all material of the course: scripts, videos, quizes,..)###
### Creative Commons By Attribution (CC BY 4.0), Antonello Lobianco          ###
################################################################################

# ## 0104 Control Flow and Functions

# ## Some stuff to set-up the environment..
cd(@__DIR__)
using Pkg
Pkg.activate(".")
## If using a Julia version different than 1.7 please uncomment and run the following line (reproductibility guarantee will hower be lost)
## Pkg.resolve()
## Pkg.instantiate() # run this if you didn't in Segment 01.01
using Random
Random.seed!(123)
using InteractiveUtils # loaded automatically when working... interactively

# ## Variables scope

# The scope of a variable is the region of code where the variable can be accessed directly (without using prefixes).
# Modules, functions, `for` and other blocks (but notably not "if" blocks) introduce an inner scope that hinerit from the scope where the block or function is defined (but not, for function, from the caller's scope).
# Variables that are defined outside any block or function are _global_ for the module where they are defined (the `Main` module if outside any other module, e.g. on the REPL), the others being _local_.
# Variables defined in a `for` block that already exists as global behave differently depending if we are working interactively or not:

g  = 2
g2 = 20
for i in 1:2
    l1 = g2                                       # l1: local, g2: global (read only)
    l1 += i
    g = i                                         # REPL/INTERACTIVE: global (from Julia 1.5), FILE MODE: local by default (with a warning `g` being already defined)
    g += i
    println("i: $i")
    println("l1: $l1")
    println("g: $g")
    for j in 1:2
        l1 += j                                   # still the local in outer loop, not a new local one
        l2 = j
        g  += j
        println("j: $j")
        println("l1 inside inner loop: $l1")
        println("l2 inside inner loop: $l2")
        println("g inside inner loop: $g")
    end
    ## println("l2 post: $l2")                     # error: l2 not defined in this scope
    println("l1 post: $l1")
    println("g post: $g")
end
## println("l1 global $l1")                        # error; l1 is not defined in the global scope
println("g in global: $g")                        # REPL/INTERACTIVE: "7", FILE MODE: "2"

function foo(i)
    l1 = g2                                       # l1: local, g2: global (read only)
    l1 += i
    g = i                                         # REPL/INTERACTIVE and FILE MODE: local by default (with a warning `g` being already defined)
    g += i
    println("i: $i")
    println("l1: $l1")
    println("g: $g")
    for j in 1:2
        l1 += j                                   # still the local in outer loop, not a new local one
        l2 = j
        g  += j
        println("j: $j")
        println("l1 inside inner loop: $l1")
        println("l2 inside inner loop: $l2")
        println("g inside inner loop: $g")
    end
    ## println("l2 post: $l2")                     # error: l2 not defined in this scope
    println("l1 post: $l1")
    println("g post: $g")
end

println("Calling foo..")
foo(10)
println("g in global: $g")                        # REPL/INTERACTIVE: "7", FILE MODE: "2"

g = 2
include("010401-varScopeExample.jl.txt")            # gives a warning !

# ## Repeated iterations: `for` and `while` loops, List Comprehension, Maps

for i in 1:2, j in 3:4             # j is the inner loop
    println("i: $i, j: $j")
end
a = 1
while true             # or condition, e.g. while a == 10
    global a += 1
    println("a: $a")
    if a == 10
        break
    else
        continue
    end
    println("This is never printed")
end

# ### List Comprehension
[ i+j for i in 1:2, j in 3:4 if j >= 4]

# ### Maps

# Apply a (possible anonymous) function to a list of arguments:

map((name,year) -> println("$name is $year year old"), ["Marc","Anna"], [25,22])

# !!! warninng
#     Don't confuse the single-line arrows used in anonymous functions (`->`) with the double-line arrow used to define a Pair (`=>`)

# We can use maps to substitute values of an array based on a dictionary:
countries = ["US","UK","IT","UK","UK"]
countryNames = Dict("IT" => "Italy", "UK" => "United Kngdom",  "US"=>"United States")
countryLongNames = map(cShortName -> countryNames[cShortName], countries)

# ## Conditional statements: if blocks and ternary operators

a = 10

if a < 4                      # use `!`, `&&` and `||` for "not", "and" and  "or" conditions
    println("a < 4")
elseif a < 8
    println("a < 8")
else
    println("a is big!")
end

a = 10

if a < 5
  b = 100
else
  b = 200
end

b

# Ternary operators
b =  a < 5 ? 100 : 200   # ? condition : if true : if false

# Short-circuit evaluation
b = 100
(a < 5) || (b = 200)    # replace an if: the second part is executed unless the first part is already true
b
(a < 50) || (b = 500)   # here is never executed
b

# !!! warning
#     Don't confuse boolean operators `&&` and `||` with their analogous `&` and `|` bitwise operators

a = 3
b = 2
bitstring(a)
bitstring(b)

##a && b     # error non boolean used in boolean context
a & b
a | b

# ## Functions

function foo(x)             # function definition
    x+2
end
foo(2)                      # function call
inlineFunction(x) = x+2
foo2 = x -> x+2             # anonymous function (aka "lambda function") and assignment to the variable foo2
foo2(2)

# A nested function:
function f1(x)
    function f2(x,y)
        x+y
    end
    f2(x,2)
end

f1(2)

# A recursive function:
function fib(n)                   # This is a naive implementation. Much faster implementations of the Fibonacci numbers exist
    if n == 0  return 0
    elseif n == 1 return 1
    else
     return fib(n-1) + fib(n-2)
    end
end
fib(4)

# ### Function arguments

# #### Positional vs keyword arguments
f(a,b=1;c=1) = a+10b+100c  # `a` and `b` are positional arguments (`b` with a default provided), `c` is a keyword argument
f(2)
f(2,c=3)

foo(a, args...;c=1) = a + length(args) + sum(args) + c
foo(1,2,3,c=4)

# Rules for positional and keyword arguments:
# - keyword arguments follow a semicolon `;` in the parameters list of the function definition
# - a positional argument without a default can not follow a positional argument with a default provided
# - the splat operator to define variable number of arguments must be the last positional argument
# - the function call must use positional arguments by position and keyword arguments by name

## foo(a::String="aaa",b::Int64) = "$a "+string(b) # error! Optional positional argument before a mandatory positional one

# ### Argument types and multiple dispatch

# Simple to understand the usage, complex to understand the deep implications
foo3(a::Int64,b::String) = a + parse(Int64,b)
foo3(2,"3")
foo3(a::String,b::Int64) = parse(Int64,a) + b
foo3("3",2)
methods(foo3)

# Multiple dispatch allows to compile a _specialised_ version JIT at run time, on the first call with the given parameters type
# We will see it again when dealing with type inheritance
# In general, unless we need to write specialised methods, no need to specify the type of the parameters. No influence on performances, this is automatically inferred (and the funciton compiled) based on the run-time type of the argument

# !!! tip Functions performances tip
#     The most important things for performances are (1) that the function is _type stable_, that is, that conditional to a specific combination of the types of the parameters the function returns the same type. This is a condition necessary to have a working chain of type inference across function calls; (2) that no (non constant) global constants are used in the function and indeed all the required information for the functio ndoing its work is embedded in the function parameters

# ### Function templates
foo3(a::T,b::String) where {T<: Number} = a + parse(T,b)             # can use T in the function body
foo3(2,"1")
foo3(1.5,"1.5")
foo4(a::Int64,b::T where T <: Number) = a + b                        # ok not used in functio nbody
foo4(a::Int64,b::Array{T} where T <: Number) = a .+ fill(T,b,2)      # wil lerror, can't use T in the function body
## foo4(2,[1,2])                                                     # run time error, T not defined

# ### Call by reference vs. call by value

# How the variable used as function argument within the function body relates to the variable used in calling the function ?
# - **call by value**: the value of the argument is copied and the function body works on a copy of the value
# - **call by reference**: the function works on the same object being referenced by the caller variable and the function argument
# - **call by sharing** (Julia): the arguments are just new local variables that bind the same object. The effects of "modifications" on the local variable on the caller's one depends on the mutability property of the object as we saw in the _Types and objects_ segment:
#   - immutable objects: we can only have that the argument is rebinded to other objects. No effects on the original caller object
#   - mutable objects: if the argument is rebinded to an other object, no effects on the caller object. If the object is modified, the caller object (being the same object) is also modified

x = 10
foo(y) = y = 1
foo(x)
x
foo(x) = x[1] = 10
x = [1,2]
foo(x)
x
foo3(x) = x = [10,20]
foo3(x)
x

# !!! info
#     Functions that modify at least one of their arguments are named, by convention, with an exclamation mark at the end of their name and the argument(s) that is (are) modified set as the first(s) argument(s)

foo!(x) = x[1] = 10 # to follow the convention

# ## `do` blocks

# Functions that accept an other function as their first parameter can be rewritten with the function itself defined in a `do` block:
using Statistics
pool(f,x,poolSize=3) = [f(x[i:i+poolSize-1]) for i in 1:length(x)-poolSize+1] # a real case, used in neural networks as pooling layer
pool(mean,[1,2,3,4,5,6])
pool(maximum,[1,2,3,4,5,6])
pool([1,2,3,4,5]) do x      # x is a local variable within the do block. We need as many local variables as the number of parameters of the inner function
    sum(x)/length(x)
end
# Using the `do` block we can call the outer function and define the inner function at the same time.
# `do` blocks are frequently used in input/output operations

max_stars_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iotfleethub
using AWS.Compat
using AWS.UUIDs

"""
    create_application(application_name, role_arn)
    create_application(application_name, role_arn, params::Dict{String,<:Any})

Creates a Fleet Hub for AWS IoT Device Management web application.  Fleet Hub for AWS IoT
Device Management is in public preview and is subject to change.

# Arguments
- `application_name`: The name of the web application.
- `role_arn`: The ARN of the role that the web application assumes when it interacts with
  AWS IoT Core.  The name of the role must be in the form AWSIotFleetHub_random_string .

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"applicationDescription"`: An optional description of the web application.
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
- `"tags"`: A set of key/value pairs that you can use to manage the web application
  resource.
"""
function create_application(
    applicationName, roleArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "POST",
        "/applications",
        Dict{String,Any}(
            "applicationName" => applicationName,
            "roleArn" => roleArn,
            "clientToken" => string(uuid4()),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_application(
    applicationName,
    roleArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "POST",
        "/applications",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "applicationName" => applicationName,
                    "roleArn" => roleArn,
                    "clientToken" => string(uuid4()),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_application(application_id)
    delete_application(application_id, params::Dict{String,<:Any})

Deletes a Fleet Hub for AWS IoT Device Management web application.  Fleet Hub for AWS IoT
Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
"""
function delete_application(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "DELETE",
        "/applications/$(applicationId)",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "DELETE",
        "/applications/$(applicationId)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_application(application_id)
    describe_application(application_id, params::Dict{String,<:Any})

Gets information about a Fleet Hub for AWS IoT Device Management web application.  Fleet
Hub for AWS IoT Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

"""
function describe_application(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "GET",
        "/applications/$(applicationId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "GET",
        "/applications/$(applicationId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_applications()
    list_applications(params::Dict{String,<:Any})

Gets a list of Fleet Hub for AWS IoT Device Management web applications for the current
account.  Fleet Hub for AWS IoT Device Management is in public preview and is subject to
change.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"nextToken"`: A token used to get the next set of results.
"""
function list_applications(; aws_config::AbstractAWSConfig=global_aws_config())
    return iotfleethub(
        "GET", "/applications"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_applications(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "GET",
        "/applications",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource_arn)
    list_tags_for_resource(resource_arn, params::Dict{String,<:Any})

Lists the tags for the specified resource.  Fleet Hub for AWS IoT Device Management is in
public preview and is subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.

"""
function list_tags_for_resource(
    resourceArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "GET",
        "/tags/$(resourceArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resourceArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "GET",
        "/tags/$(resourceArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Adds to or modifies the tags of the specified resource. Tags are metadata which can be used
to manage a resource.  Fleet Hub for AWS IoT Device Management is in public preview and is
subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tags`: The new or modified tags for the resource.

"""
function tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return iotfleethub(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tags" => tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    resourceArn,
    tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "POST",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tags" => tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Removes the specified tags (metadata) from the resource.  Fleet Hub for AWS IoT Device
Management is in public preview and is subject to change.

# Arguments
- `resource_arn`: The ARN of the resource.
- `tag_keys`: A list of the keys of the tags to be removed from the resource.

"""
function untag_resource(
    resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resourceArn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "DELETE",
        "/tags/$(resourceArn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_application(application_id)
    update_application(application_id, params::Dict{String,<:Any})

Updates information about a Fleet Hub for a AWS IoT Device Management web application.
Fleet Hub for AWS IoT Device Management is in public preview and is subject to change.

# Arguments
- `application_id`: The unique Id of the web application.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"applicationDescription"`: An optional description of the web application.
- `"applicationName"`: The name of the web application.
- `"clientToken"`: A unique case-sensitive identifier that you can provide to ensure the
  idempotency of the request. Don't reuse this client token if a new idempotent request is
  required.
"""
function update_application(
    applicationId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return iotfleethub(
        "PATCH",
        "/applications/$(applicationId)",
        Dict{String,Any}("clientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_application(
    applicationId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return iotfleethub(
        "PATCH",
        "/applications/$(applicationId)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("clientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
using Base.Test, DataStreams, Nulls

import Base: ==
==(a::DataStreams.Data.Schema, b::DataStreams.Data.Schema) = DataStreams.Data.types(a) == DataStreams.Data.types(b) && DataStreams.Data.header(a) == DataStreams.Data.header(b) && size(a) == size(b)

@static if isdefined(Core, :NamedTuple)

mutable struct Source{T}
    sch::DataStreams.Data.Schema
    nt::T
end

DataStreams.Data.schema(s::Source) = s.sch
DataStreams.Data.isdone(s::Source, row, col, rows, cols) = col > cols || (isnull(rows) ? row > length(s.nt[col]) : row > rows)
DataStreams.Data.streamfrom(s::Source, ::Type{DataStreams.Data.Column}, ::Type{T}, row, col) where {T} = s.nt[col]
DataStreams.Data.streamfrom(s::Source, ::Type{DataStreams.Data.Field}, ::Type{T}, row, col) where {T} = s.nt[col][row]

mutable struct Sink{T}
    nt::T
end

I = (id = Int64[1, 2, 3, 4, 5],
firstname = (Union{String, Null})["Benjamin", "Wayne", "Sean", "Charles", null],
lastname = String["Chavez", "Burke", "Richards", "Long", "Rose"],
salary = (Union{Float64, Null})[null, 46134.1, 45046.2, 30555.6, 88894.1],
rate = Float64[39.44, 33.8, 15.64, 17.67, 34.6],
hired = (Union{Date, Null})[Date("2011-07-07"), Date("2016-02-19"), null, Date("2002-01-05"), Date("2008-05-15")],
fired = DateTime[DateTime("2016-04-07T14:07:00"), DateTime("2015-03-19T15:01:00"), DateTime("2006-11-18T05:07:00"), DateTime("2002-07-18T06:24:00"), DateTime("2007-09-29T12:09:00")]
)
J = (; :_0=>["0"], (Symbol("_$i")=>[i] for i = 1:501)...);
K = (; (Symbol("_$i")=>[i] for i = 1:501)...);

nms(::NamedTuple{names}) where {names} = names
I_L = Source(DataStreams.Data.Schema(collect(map(eltype, I)), nms(I), 5), I);
I_M = Source(DataStreams.Data.Schema(collect(map(eltype, I)), nms(I), null), I);

J_L = Source(DataStreams.Data.Schema(collect(map(eltype, J)), nms(J), 1), J);
J_M = Source(DataStreams.Data.Schema(collect(map(eltype, J)), nms(J), null), J);

K_L = Source(DataStreams.Data.Schema(collect(map(eltype, K)), nms(K), 1), K);
K_M = Source(DataStreams.Data.Schema(collect(map(eltype, K)), nms(K), null), K);

Sink(sch::DataStreams.Data.Schema, S, append, args...; reference::Vector{UInt8}=UInt8[]) = Sink(NamedTuple(sch, S, append, args...; reference=reference))
Sink(sink, sch::DataStreams.Data.Schema, S, append; reference::Vector{UInt8}=UInt8[]) = Sink(NamedTuple(sink.nt, sch, S, append; reference=reference))
DataStreams.Data.streamtypes(::Type{Sink}) = [DataStreams.Data.Column, DataStreams.Data.Field]
DataStreams.Data.weakrefstrings(::Type{Sink}) = true
DataStreams.Data.streamto!(sink::Sink, ::Type{DataStreams.Data.Field}, val, row, col) =
(C = getfield(sink.nt, col); row > length(C) ? push!(C, val) : setindex!(C, val, row); return)
DataStreams.Data.streamto!(sink::Sink, ::Type{DataStreams.Data.Column}, column, col) =
append!(getfield(sink.nt, col), column)

end # isdefined

@testset "DataStreams" begin

@testset "Data.Schema" begin

sch = DataStreams.Data.Schema()
@test size(sch) == (0, 0)
@test DataStreams.Data.header(sch) == String[]
@test DataStreams.Data.types(sch) == ()

sch = DataStreams.Data.Schema((Int,))
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["Column1"]
@test DataStreams.Data.types(sch) == (Int,)

sch = DataStreams.Data.Schema([Int])
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["Column1"]
@test DataStreams.Data.types(sch) == (Int,)

sch = DataStreams.Data.Schema((Int,), ["col1"])
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)

@test_throws ArgumentError DataStreams.Data.Schema((String,), ["col1", "col2"])

getR(::DataStreams.Data.Schema{R}) where {R} = R
getT(::DataStreams.Data.Schema{R, T}) where {R, T} = T

sch = DataStreams.Data.Schema((Int,), ["col1"], 1)
@test size(sch) == (1, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)
@test getR(sch) == true
@test getT(sch) == Tuple{Int}

sch = DataStreams.Data.Schema((String,), ["col1"], null)
@test size(sch) == (null, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (String,)
@test getR(sch) == false
@test getT(sch) == Tuple{String}

@test_throws ArgumentError DataStreams.Data.Schema((String,), ["col1"], -1)

sch = DataStreams.Data.Schema((Int,), ["col1"], 0, Dict("hey"=>"ho"))
@test size(sch) == (0, 1)
@test DataStreams.Data.header(sch) == String["col1"]
@test DataStreams.Data.types(sch) == (Int,)
@test DataStreams.Data.metadata(sch) == Dict("hey"=>"ho")

sch = DataStreams.Data.Schema((Int,), ["col1"])
@test sch["col1"] == 1

end # @testset "Data.Schema"

@testset "Data.transform" begin

# knownrows vs. not
sch = DataStreams.Data.Schema((Int,), ["col1"], 1)
sch2, trans = DataStreams.Data.transform(sch, Dict{Int, Function}(), true)
@test sch == sch2
@test trans == (identity,)

sch = DataStreams.Data.Schema((Int, Float64, String,), ["col1", "col2", "col3"], null)
sch2, trans = DataStreams.Data.transform(sch, Dict{Int, Function}(), true)
@test sch == sch2
@test trans == (identity, identity, identity)

sch = DataStreams.Data.Schema((Int,), ["col1"], 1)
f = x->string(x)
sch2, trans = DataStreams.Data.transform(sch, Dict(1=>f), true)
@test DataStreams.Data.header(sch) == DataStreams.Data.header(sch2)
@test DataStreams.Data.types(sch2) == (String,)
@test trans == (f,)

sch = DataStreams.Data.Schema((Int, Float64, String), ["col1", "col2", "col3"], null)
f = x->parse(Int, x)
sch2, trans = DataStreams.Data.transform(sch, Dict("col3"=>f), true)
@test DataStreams.Data.header(sch) == DataStreams.Data.header(sch2)
@test DataStreams.Data.types(sch2) == (Int, Float64, Int)
@test trans == (identity, identity, f)

sch = DataStreams.Data.Schema((Int, Float64, String), ["col1", "col2", "col3"], null)
sch2, trans = DataStreams.Data.transform(sch, Dict{Int, Function}(), false)
@test DataStreams.Data.header(sch) == DataStreams.Data.header(sch2)
@test DataStreams.Data.types(sch2) == (Int, Float64, String)
@test trans == (identity, identity, identity)

sch = DataStreams.Data.Schema((Union{String, Null},), ["col1"])
sch2, trans = DataStreams.Data.transform(sch, Dict{Int, Function}(), false)
@test DataStreams.Data.header(sch) == DataStreams.Data.header(sch2)
@test DataStreams.Data.types(sch2) == (Union{String, Null},)
@test trans == (identity,)

sch = DataStreams.Data.Schema((Int,), ["col1"], 1)
sch2, trans = DataStreams.Data.transform(sch, Dict("col1"=>sin), false)
@test DataStreams.Data.header(sch) == DataStreams.Data.header(sch2)
@test DataStreams.Data.types(sch2) == (Float64,)
@test trans == (sin,)

end # @testset "Data.transform"

@static if isdefined(Core, :NamedTuple)

@testset "Data.stream!" begin

### ALL PAIRS TESTING

# constructed Sink vs. constructed on-the-fly: A, B
# append vs. not; C, D
# w/ transform functions vs. not: E, F
# Data.Field vs. Data.Column: G, H
# <500, >500, homogenous types code path: I, J, K
# knownrows vs. not: L, M

# incompatible source => constructed Sink
@test_throws ArgumentError DataStreams.Data.stream!(I_L, Sink(I))

# incompatible source => otf Sink
@test_throws ArgumentError DataStreams.Data.stream!(I_L, Sink, I)

## Data.Field
DataStreams.Data.streamtype(::Type{Source}, ::Type{DataStreams.Data.Field}) = true
# A, C, E, G, I, L: append to constructed Sink w/ transforms via Data.Field w/ Source=I_J
source = I_L
sink = Sink(deepcopy(I))
transforms = Dict(1=>x->x+1)
DataStreams.Data.stream!(source, sink; append=true, transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (10, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Null}, String, Union{Float64, Null}, Float64, Union{Date, Null}, DateTime)
@test sink.nt.id == [1,2,3,4,5,2,3,4,5,6]

# A, C, F, G, J, L: append to constructed Sink w/o transforms via Data.Field w/ Source=J_L
source = J_L;
sink = Sink(deepcopy(J));
Data.stream!(source, sink; append=true);

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (2, 502)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:501)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:501)...)
@test sink.nt._0 == ["0", "0"]
@test sink.nt._1 == [1, 1]

# A, D, F, G, J, M: replace constructed Sink w/o transforms via Data.Field w/ Source=J_M
source = J_M;
sink = Sink(deepcopy(J));
Data.stream!(source, sink)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 502)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:501)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:501)...)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [1]

# B, D, E, G, J, L: replace otf Sink w/ transforms via Data.Field w/ Source=J_L
source = J_L
sink = Sink(deepcopy(J))
transforms = Dict(2=>x->x+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 502)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:501)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:501)...)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [2]

# B, D, E, G, K, L: replace otf Sink w/ transforms via Data.Field w/ Source=K_L
source = K_L
sink = Sink(deepcopy(K))
transforms = Dict(2=>x->x+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 501)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:501)
@test DataStreams.Data.types(sch) == ((Int for i = 1:501)...)
@test sink.nt._1 == [1]
@test sink.nt._2 == [3]

## Data.Column
# A, D, E, H, K, M: replace constructed Sink w/ transforms via Data.Column w/ Source=K_M
DataStreams.Data.streamtype(::Type{<:Source}, ::Type{DataStreams.Data.Column}) = true
source = K_M
sink = Sink(deepcopy(K))
transforms = Dict(2=>x->x+1)
Data.stream!(source, sink; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 501)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:501)
@test DataStreams.Data.types(sch) == ((Int for i = 1:501)...)
@test sink.nt._1 == [1]
@test sink.nt._2 == [3]

# B, C, F, H, K, M: append to otf Sink w/o transforms via Data.Column w/ Source=K_M
source = K_M
sink = Sink(deepcopy(K))
Data.stream!(source, Sink, sink.nt; append=true)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (2, 501)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 1:501)
@test DataStreams.Data.types(sch) == ((Int for i = 1:501)...)
@test sink.nt._1 == [1, 1]
@test sink.nt._2 == [2, 2]

# B, D, F, H, I, M: replace otf Sink w/o transforms via Data.Column w/ Source=I_M
source = I_M
sink = Sink(deepcopy(I))
transforms = Dict()
Data.stream!(source, Sink, sink.nt)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (5, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Null}, String, Union{Float64, Null}, Float64, Union{Date, Null}, DateTime)
@test sink.nt.id == [1,2,3,4,5]

# B, D, F, H, I, L: replace otf Sink w/o transforms via Data.Column w/ Source=I_L
source = I_L
sink = Sink(deepcopy(I))
Data.stream!(source, Sink, sink.nt)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (5, 7)
@test DataStreams.Data.header(sch) == ["id","firstname","lastname","salary","rate","hired","fired"]
@test DataStreams.Data.types(sch) == (Int64, Union{String, Null}, String, Union{Float64, Null}, Float64, Union{Date, Null}, DateTime)
@test sink.nt.id == [1,2,3,4,5]

# B, D, E, H, J, L: replace otf Sink w/ transforms via Data.Field w/ Source=J_L
source = J_L
sink = Sink(deepcopy(J))
transforms = Dict(2=>x->x+1)
Data.stream!(source, Sink, sink.nt; transforms=transforms)

sch = DataStreams.Data.schema(sink.nt)
@test size(sch) == (1, 502)
@test DataStreams.Data.header(sch) == collect("_$i" for i = 0:501)
@test DataStreams.Data.types(sch) == (String, (Int for i = 1:501)...)
@test sink.nt._0 == ["0"]
@test sink.nt._1 == [2]

end # @testset "Data.stream!"
end # isdefined

# @testset "DataStreams NamedTuple" begin
#
# end
end # @testset "DataStreams"


text:
module MakeOnlineStatsDocs
reload("OnlineStats")
using OnlineStats

rootdir = Pkg.dir("OnlineStats")

#-------------------------------------------------------------------# Generate api.md
api = rootdir * "/doc/api.md"
rm(api)
touch(api)
file = open(api, "r+")
write(file, "<!--- Generated at " * string(now()) * ".  Don't edit --->\n")
write(file, "# API\n\n")
info("The following items are included in the output file:\n")
nms = setdiff(names(OnlineStats), [:OnlineStats])  # Vector{Symbol} of names
for nm in nms
    @eval obj = OnlineStats.$nm
    d = Docs.doc(obj)
    if d != nothing
        println(nm)
        write(file, "## " * string(nm) * "\n" * Markdown.plain(d) * "\n")
    end
end
close(file)

#----------------------------------------------------------------# Figure for Weights
# info("Figure for Weights")
# using Plots; pyplot()
# o1 = EqualWeight()
# o2 = ExponentialWeight(.2)
# o3 = BoundedEqualWeight(.2)
# o4 = LearningRate(.5)
# ovec = [o1, o2, o3, o4]
# map(OnlineStats.updatecounter!, ovec)
#
# plt_wt = plot([0], map(OnlineStats.weight, ovec)', w = 3,
# label = ["EqualWeight" "ExponentialWeight(.2)" "BoundedEqualWeight(.2)" "LearningRate(.5)"],
# xlabel = "nobs", ylabel = "weight value", layout = 4, ylims = (0, 1))
# for i in 1:50
#     for o in ovec
#         OnlineStats.updatecounter!(o)
#     end
#     push!(plt_wt, nobs(o1), map(OnlineStats.weight, ovec))
# end
# png(plt_wt, rootdir * "/doc/images/weights.png")

#-------------------------------------------------------------# push site to gh-pages
cd(rootdir)
run(`mkdocs gh-deploy --clean`)
end  # MakeOnlineStatsDocs


text:
using QMDP
using POMDPs
using POMDPModels
using POMDPModelTools
using BeliefUpdaters
using POMDPTesting
using Test
using Random

pomdp = TigerPOMDP()
solver = QMDPSolver()

@requirements_info solver pomdp

solver.verbose = true
policy = solve(solver, pomdp)

solver.verbose = false
rng = MersenneTwister(11)

bu = updater(policy)
sd = initialstate_distribution(pomdp)
b = initialize_belief(bu, sd)

a = action(policy, b)
v = value(policy, b)

s = initialstate(pomdp, rng)
sp, o = generate_so(pomdp, s, a, rng)
bp = update(bu, b, a, o)
@test isa(bp, DiscreteBelief)

r = test_solver(solver, pomdp)
@test isapprox(r, 17.711, atol=1e-2)

println("There should be a warning here: ")
solve(solver, pomdp, verbose=true)