text:
using PyPlot
using CSV
include("structures.jl")
include("utilitaires.jl")
include("parser.jl")
include("concorde.jl")
include("partitionnement_mono.jl")

function main()
    # Paramètres
    vitesseCamion = 1
    vitesseDrone = 3
    listeInstances = readdir("experimentations")
    nbInstances = length(listeInstances)

    # Première compilation
    instance = parserNousBis(listeInstances[1])
    instance.vitesseCamion = vitesseCamion
    instance.vitesseDrone = vitesseDrone
    concordeFast(instance)
    AEP(instance)
    EP(instance)

    # Expérimentation
    taillesInstance = trunc.(Int64, [1:120;])
    listeTempsTSP = Vector{Float64}(undef, 120)
    listeValeursTSP = Vector{Float64}(undef, 120)
    listeTempsAEP = Vector{Float64}(undef, 120)
    listeValeursAEP = Vector{Float64}(undef, 120)
    listeTempsEP = Vector{Float64}(undef, 120)
    listeValeursEP = Vector{Float64}(undef, 120)

    # Init
    for i in 1:length(taillesInstance)
        listeTempsTSP[i] = 0.0
        listeTempsAEP[i] = 0.0
        listeTempsEP[i] = 0.0
        listeValeursTSP[i] = 0.0
        listeValeursAEP[i] = 0.0
        listeValeursEP[i] = 0.0
    end

    # Début
    for i in 1:nbInstances
        instance = parserNousBis(listeInstances[i])
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone
        println(" -> $i : $(instance.nbPoints)")
        tempsTSP = @elapsed concordeFast(instance)
        valTSP = instance.solution.tempsParcours
        listeTempsTSP[instance.nbPoints] += tempsTSP
        listeValeursTSP[instance.nbPoints] += valTSP
        tempsAEP = @elapsed AEP(instance)
        valAEP = instance.solution.tempsParcours
        listeTempsAEP[instance.nbPoints] += tempsAEP
        listeValeursAEP[instance.nbPoints] += valAEP
        concordeFast(instance)
        tempsEP = @elapsed EP(instance)
        valEP = instance.solution.tempsParcours
        listeTempsEP[instance.nbPoints] += tempsEP
        listeValeursEP[instance.nbPoints] += valEP
    end

    # Lissage
    listeTempsTSP = listeTempsTSP/20
    listeTempsAEP = listeTempsAEP/20
    listeTempsEP = listeTempsEP/20
    listeValeursTSP = listeValeursTSP/20
    listeValeursAEP = listeValeursAEP/20
    listeValeursEP = listeValeursEP/20

    graphique = plt.figure() # Sortie graphique unique
    ax = graphique.gca()
    #ax.title("CPUt for several instance sizes")
    xlabel(L"Instance size ($n$)")
    ylabel(L"Time ($s$)")
    xticks(taillesInstance[10:end])
    ax.plot(taillesInstance[10:end], listeTempsTSP[10:end], label="TSP")
    ax.plot(taillesInstance[10:end], listeTempsAEP[10:end], label="A1")
    ax.plot(taillesInstance[10:end], listeTempsEP[10:end], label="A2")
    ax.grid(true)
    ax.legend(loc="best")
    graphique2 = plt.figure() # Sortie graphique unique
    ax2 = graphique2.gca()
    #ax2.title(L"Value of $f_1$ for several instance sizes")
    xlabel(L"Instance size ($n$)")
    ylabel(L"$z_1$")
    xticks(taillesInstance[10:end])
    ax2.plot(taillesInstance[10:end], listeValeursTSP[10:end], label="TSP")
    ax2.plot(taillesInstance[10:end], listeValeursAEP[10:end], label="A1")
    ax2.plot(taillesInstance[10:end], listeValeursEP[10:end], label="A2")
    ax2.grid(true)
    ax2.legend(loc="best")
end

function expReverse()
    # Paramètres
    vitesseCamion = 1
    vitesseDrone = 3
    listeInstances = readdir("experimentations")
    nbInstances = length(listeInstances)

    taillesInstance = trunc.(Int64, [1:120;])
    A1gapZ1 = Vector{Float64}(undef, 120)
    A1gapZ2 = Vector{Float64}(undef, 120)
    A2gapZ1 = Vector{Float64}(undef, 120)
    A2gapZ2 = Vector{Float64}(undef, 120)

    # Début
    for i in 1:nbInstances
        instance = parserNousBis(listeInstances[i])
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone = vitesseDrone

        println(" -> $i : $(instance.nbPoints)")
        tempsTSP = @elapsed concordeFast(instance)
        tempsParcoursTSP = instance.solution.tempsParcours
        ordreParcoursTSP = instance.solution.ordreVisiteCamion

        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        #println(" Parcours TSP : ", instance.solution.ordreVisiteCamion)
        #println(" Temps parcours TSP : ", instance.solution.tempsParcours)
        tempsAEP = @elapsed AEP(instance)
        #println(" AEP ")
        #println(" Parcours (z1) : ", instance.solution.tempsParcours)
        #println(" Parcours (z2) : ", instance.solution.tempsAttente)
        aepZ1 = instance.solution.tempsParcours
        aepZ2 = instance.solution.tempsAttente

        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        #println(" Parcours TSP : ", instance.solution.ordreVisiteCamion)
        #println(" Temps parcours TSP : ", instance.solution.tempsParcours)
        tempsEP = @elapsed EP(instance)
        #println(" EP ")
        #println(" Parcours (z1) : ", instance.solution.tempsParcours)
        #println(" Parcours (z2) : ", instance.solution.tempsAttente)
        epZ1 = instance.solution.tempsParcours
        epZ2 = instance.solution.tempsAttente

        reverse!(ordreParcoursTSP)
        #println(" REVERSE ")

        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        #println(" Parcours TSP : ", instance.solution.ordreVisiteCamion)
        #println(" Temps parcours TSP : ", instance.solution.tempsParcours)
        tempsAEP = @elapsed AEP(instance)
        #println(" AEP ")
        #println(" Parcours inverse (z1) : ", instance.solution.tempsParcours)
        #println(" Parcours inverse (z2) : ", instance.solution.tempsAttente)
        aepZ1reverse = instance.solution.tempsParcours
        aepZ2reverse = instance.solution.tempsAttente

        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        #println(" Parcours TSP : ", instance.solution.ordreVisiteCamion)
        #println(" Temps parcours TSP : ", instance.solution.tempsParcours)
        tempsEP = @elapsed EP(instance)
        #println(" EP ")
        #println(" Parcours inverse (z1) : ", instance.solution.tempsParcours)
        #println(" Parcours inverse (z2) : ", instance.solution.tempsAttente)
        epZ1reverse = instance.solution.tempsParcours
        epZ2reverse = instance.solution.tempsAttente

        #println(" Gap A1 z1 : ", ( 1 - ( min(aepZ1reverse, aepZ1)/max(aepZ1reverse, aepZ1) ) )*100, "%")
        #println(" Gap A1 z2 : ", (1-(min(aepZ2reverse, aepZ2)/max(aepZ2reverse, aepZ2)))*100, "%")
        #println(" Gap A2 z1 : ", (1-(min(epZ1reverse, epZ1)/max(epZ1reverse, epZ1)))*100, "%")
        #println(" Gap A2 z2 : ", (1-(min(epZ2reverse, epZ2)/max(epZ2reverse, epZ2)))*100, "%")
        A1gapZ1[instance.nbPoints] = (1-(min(aepZ1reverse, aepZ1)/max(aepZ1reverse, aepZ1)))*100
        A1gapZ2[instance.nbPoints] = (1-(min(aepZ2reverse, aepZ2)/max(aepZ2reverse, aepZ2)))*100
        A2gapZ1[instance.nbPoints] = (1-(min(epZ1reverse, epZ1)/max(epZ1reverse, epZ1)))*100
        A2gapZ2[instance.nbPoints] = (1-(min(epZ2reverse, epZ2)/max(epZ2reverse, epZ2)))*100
    end

    # Lissage
    A1gapZ1/20
    A1gapZ2/20
    A2gapZ1/20
    A2gapZ2/20

    meanA1z1 = 0.0
    meanA1z2 = 0.0
    meanA2z1 = 0.0
    meanA2z2 = 0.0
    for i in 10:120
        meanA1z1 += A1gapZ1[i]
        meanA1z2 += A1gapZ2[i]
        meanA2z1 += A2gapZ1[i]
        meanA2z2 += A2gapZ2[i]
    end
    meanA1z1 = meanA1z1/111
    meanA1z2 = meanA1z2/111
    meanA2z1 = meanA2z1/111
    meanA2z2 = meanA2z2/111

    graphique = plt.figure() # Sortie graphique unique
    ax = graphique.gca()
    # ax.title(M"Gap between $R_{TSP}$ and $R'_{TSP}$ (A1)")
    xlabel(L"Instance size ($n$)")
    ylabel(L"Gap ($\%$)")
    xticks(taillesInstance[10:end])
    ax.plot(taillesInstance[10:end], A1gapZ1[10:end], label=L"$z_1$", color="C1") # C1 : orange light
    ax.plot(taillesInstance[10:end], A1gapZ2[10:end], label=L"$z_2$", color="C0") # C0 : blue light
    ax.grid(true)
    axhline(y=meanA1z1, linestyle="--", label=L"mean $z_1$", color="C1")
    axhline(y=meanA1z2, linestyle="--", label=L"mean $z_2$", color="C0")
    ax.legend(loc="best")
    xticks([10,20,30,40,50,60,70,80,90,100,110,120])
    graphique2 = plt.figure() # Sortie graphique unique
    ax2 = graphique2.gca()
    #ax2.title(L"Value of $f_1$ for several instance sizes")
    xlabel(L"Instance size ($n$)")
    ylabel(L"Gap ($\%$)")
    xticks(taillesInstance[10:end])
    ax2.plot(taillesInstance[10:end], A2gapZ1[10:end], label=L"$z_1$", color="C1")
    ax2.plot(taillesInstance[10:end], A2gapZ2[10:end], label=L"$z_2$", color="C0")
    axhline(y=meanA2z1, linestyle="--", label=L"mean $z_1$", color="C1")
    axhline(y=meanA2z2, linestyle="--", label=L"mean $z_2$", color="C0")
    ax2.grid(true)
    ax2.legend(loc="best")
    xticks([10,20,30,40,50,60,70,80,90,100,110,120])
end

function etudeVitesse()
    # Paramètres
    listeInstances = readdir("experimentations")
    nbInstances = length(listeInstances)
    listeAlphas = [1:0.25:50;]
    listeZ1AEP = Vector{Float64}(undef, length(listeAlphas))
    listeZ1EP = Vector{Float64}(undef, length(listeAlphas))

    for j in 1:nbInstances
    println(" --> !!!! $j !!!!")
    vitesseCamion = 1
    vitesseDrone = 1
    instance = parserNousBis(listeInstances[j])
    instance.vitesseCamion = vitesseCamion
    instance.vitesseDrone = vitesseDrone
    concordeFast(instance)
    tempsParcoursTSP = instance.solution.tempsParcours
    ordreParcoursTSP = instance.solution.ordreVisiteCamion
    for i in 1:length(listeAlphas)
        println(" -> alpha : ", instance.vitesseDrone/instance.vitesseCamion)
        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        AEP(instance)
        z1AEP = instance.solution.tempsParcours
        println(" -> z_1 AEP : $z1AEP")
        instance.solution = Solution(ordreParcoursTSP,[],tempsParcoursTSP,0)
        EP(instance)
        z1EP = instance.solution.tempsParcours
        println(" -> z_1 EP : $z1EP")
        listeZ1AEP[i] += z1AEP
        listeZ1EP[i] += z1EP
        instance.vitesseCamion = vitesseCamion
        instance.vitesseDrone += 0.25
    end
    end

    listeZ1AEP = listeZ1AEP/nbInstances
    listeZ1EP = listeZ1EP/nbInstances

    println(listeAlphas)
    println(listeZ1AEP)
    println(listeZ1EP)

    graphique = plt.figure() # Sortie graphique unique
    ax = graphique.gca()
    # ax.title(M"Gap between $R_{TSP}$ and $R'_{TSP}$ (A1)")
    xlabel(L"$\alpha = \frac{speed_{drone}}{speed_{truck}}$")
    ylabel(L"z_1")
    ax.plot(listeAlphas, listeZ1AEP, label=L"A1", color="C1", linestyle="-") # C1 : orange light
    ax.plot(listeAlphas, listeZ1EP, label=L"A2", color="C0", linestyle="--") # C0 : blue light
    ax.grid(true)
    ax.legend(loc="best")
end

#### Expe Bi

function expBi()
    listeInstances = readdir("experimentations")
    nbInstances = length(listeInstances)
    cpt = 0
    cpt2 = 0
    for i in 1:1000
        print(" --> $i : ")
        instance = parserNousBis(listeInstances[i])
        instance.vitesseCamion = 1
        instance.vitesseDrone = 0.5
        concordeFast(instance)
        cptTemp, cpt2Temp = tempBi(instance)
        cpt += cptTemp
        cpt2 += cpt2Temp
    end
    println(" Same : $cpt")
    println(" Not same : $cpt2")
end

function tempBi(instance)
    nb1 = tempV1(instance)
    nb2 = tempV2(instance)
    if nb1 == nb2
        println(" Same")
        cpt = 1
        cpt2 = 0
    else
        println(" Not same")
        cpt = 0
        cpt2 = 1
    end
    return cpt, cpt2
end

function tempV1(instance)
    tempsMinIJ, operationMinIJ, attenteIJ = calculMeilleuresOperationsV1(instance)
    #display(tempsMinIJ)
    #display(attenteIJ)
    #display(operationMinIJ)

    V = []
    A = []
    P = []
    for i in 1:length(instance.solution.ordreVisiteCamion)
        push!(V, [])
        push!(A, [])
        push!(P, [])
    end

    push!(V[1], 0.0)
    push!(A[1], 0.0)
    push!(P[1], -1) # A changer
    for i in 2:length(V)
        for k in 1:(i-1)
            for m in 1:length(V[k])
                for n in 1:length(tempsMinIJ[k,i])
                    temps = V[k][m] + (tempsMinIJ[k,i])[n]
                    attente = A[k][m] + (attenteIJ[k,i])[n]
                    pred = k
                    push!(V[i], temps)
                    push!(A[i], attente)
                    push!(P[i], pred)
                end
            end
        end
    #plot(V[i], A[i], marker="x", linestyle="", color="C0")
    #grid(true)
    #title("2-objectifs")
    #xlabel("Temps de parcours total")
    #ylabel("Temps d'attente total")
    #clf()
    end

    #println(" V = $V")
    #println(" A = $A")
    #println(" P = $P")

    #plot(V[end], A[end], marker="x", linestyle="", color="C0")
    #grid(true)
    #title("2-objectifs")
    #xlabel("Temps de parcours total")
    #ylabel("Temps d'attente total")
    #println("\n Nombre de points (génération exhaustive) : ", length(V[end]))
    #clf()

    ND = []
    compteur = 0
    for i in 1:length(V[end]) # On veut tester si le point i est dominé par un point j
        estDomine = false
        for j in 1:length(V[end])
            if j != i
                if round(V[end][j], digits=5) <= round(V[end][i], digits=5) && round(A[end][j], digits=5) <= round(A[end][i], digits=5) # le point j est meilleur sur les 2 objectifs -> il domine le point i
                    estDomine = true
                    break
                end
            end
        end
        if estDomine == false
            push!(ND, i)
            compteur += 1
            #println(" $(compteur)) z_1 = $(V[end][i]); z_2 = $(A[end][i])")
            #plot(V[end][i], A[end][i], marker="x", linestyle="", color="C1")
        end
    end
    return length(ND)
end

function tempV2(instance)
    tempsMinIJ, operationMinIJ, attenteIJ = calculMeilleuresOperationsV2(instance)
    #display(tempsMinIJ)
    #display(attenteIJ)
    #display(operationMinIJ)

    V = []
    A = []
    P = []
    for i in 1:length(instance.solution.ordreVisiteCamion)
        push!(V, [])
        push!(A, [])
        push!(P, [])
    end

    push!(V[1], 0.0)
    push!(A[1], 0.0)
    push!(P[1], -1) # A changer
    for i in 2:length(V)
        for k in 1:(i-1)
            for m in 1:length(V[k])
                for n in 1:length(tempsMinIJ[k,i])
                    temps = V[k][m] + (tempsMinIJ[k,i])[n]
                    attente = A[k][m] + (attenteIJ[k,i])[n]
                    pred = k
                    push!(V[i], temps)
                    push!(A[i], attente)
                    push!(P[i], pred)
                end
            end
        end
    #plot(V[i], A[i], marker="x", linestyle="", color="C0")

    end

    ND = []
    compteur = 0
    for i in 1:length(V[end]) # On veut tester si le point i est dominé par un point j
        estDomine = false
        for j in 1:length(V[end])
            if j != i
                if round(V[end][j], digits=5) <= round(V[end][i], digits=5) && round(A[end][j], digits=5) <= round(A[end][i], digits=5) # le point j est meilleur sur les 2 objectifs -> il domine le point i
                    estDomine = true
                    break
                end
            end
        end
        if estDomine == false
            push!(ND, i)
            compteur += 1
            #println(" $(compteur)) z_1 = $(V[end][i]); z_2 = $(A[end][i])")
            #plot(V[end][i], A[end][i], marker="x", linestyle="", color="C1")
        end
    end
    return length(ND)
end

expReverse()

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
## TODO: pointmass rules
## TODO: save sigma vectors in meta to limit allocations

@rule Flow(:out, Marginalisation) (m_in::MvNormalMeanCovariance, meta::FlowMeta{M,Linearization}) where { M } = begin

    # extract parameters
    μ_in, Σ_in = mean_cov(m_in)

    # extract model
    model = getmodel(meta)

    # calculate new parameters
    (μ_out, J)   = forward_jacobian(model, μ_in)
    Σ_out   = J * Σ_in * J'

    # return distribution
    return MvNormalMeanCovariance(μ_out, Σ_out)

end

@rule Flow(:out, Marginalisation) (m_in::MvNormalMeanPrecision, meta::FlowMeta{M,Linearization}) where { M } = begin

    # extract parameters
    μ_in, Λ_in = mean_precision(m_in)

    # extract model
    model = getmodel(meta)

    # calculate new parameters
    μ_out   = forward(model, μ_in)
    Ji      = inv_jacobian(model, μ_in)
    Λ_out   = Ji' * Λ_in * Ji

    # return distribution
    return MvNormalMeanPrecision(μ_out, Λ_out)

end

@rule Flow(:out, Marginalisation) (m_in::MvNormalWeightedMeanPrecision, meta::FlowMeta{M,Linearization}) where { M } = begin

    # extract parameters
    μ_in, Λ_in = mean_precision(m_in)

    # extract model
    model = getmodel(meta)

    # calculate new parameters
    μ_out   = forward(model, μ_in)
    Ji      = inv_jacobian(model, μ_in)
    Λ_out   = Ji' * Λ_in * Ji

    # return distribution
    return MvNormalMeanPrecision(μ_out, Λ_out)

end

@rule Flow(:out, Marginalisation) (m_in::MultivariateNormalDistributionsFamily, meta::FlowMeta{M,Unscented}) where { M } = begin

    # extract parameters
    μ_in, Σ_in = mean_cov(m_in)

    # extract model
    model = getmodel(meta)
    T = eltype(model)

    # extract parameters of linearization
    approximation = getapproximation(meta)
    λ  = getλ(approximation)
    L  = getL(approximation)
    Wm = getWm(approximation)
    Wc = getWc(approximation)

    # calculate sigma points/vectors
    sqrtΣ = sqrt((L + λ)*Σ_in)
    χ = Vector{Vector{T}}(undef, 2*L + 1)
    for k = 1:length(χ)
        χ[k] = copy(μ_in)
    end
    for l = 2:L+1
        χ[l]     .+= sqrtΣ[l-1,:]
        χ[L + l] .-= sqrtΣ[l-1,:]
    end

    # transform sigma points
    Y = forward.(model, χ)

    # calculate new parameters
    μ_out = zeros(T, L)
    Σ_out = zeros(T, L, L)
    for k = 1:2*L+1
        μ_out .+= Wm[k] .* Y[k]
    end
    for k = 1:2*L+1
        Σ_out .+= Wc[k] .* ( Y[k] - μ_out ) *  ( Y[k] - μ_out )'
    end

    # return distribution
    return MvNormalMeanCovariance(μ_out, collect(Hermitian(Σ_out)))

end

ext: jl


text:
using Plots
pgfplots()

import NaNMath
# import BenchmarkProfiles

# Customized performance_profile method, derived from BenchmarkProfiles.jl
# Copyright (c) 2016: Dominique Orban.
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

function performance_profile(T0::Array{Float64,2}, labels::Vector{String}; logscale::Bool=true, title::String="", ymax = 1.1, xmax = NaN, linewidth=1, kwargs...)

    T0[isinf.(T0)] = NaN
    T0[T0 .< 0] = NaN

    # Remove instances for which all solvers have NaN
    T = T0[.!(isnan.(T0[:,1]) .& isnan.(T0[:,2])), :]

    # Number of instances and number of solvers
    (np, ns) = size(T)

    # Minimal (i.e., best) performance per solver
    minperf = mapslices(NaNMath.minimum, T, 2)

    # Compute ratios and divide by smallest element in each row
    r = zeros(np, ns)
    for p = 1 : np
        r[p, :] = T[p, :] / minperf[p]
    end

    logscale && (r = log2.(r))
    max_ratio = NaNMath.maximum(r)

    # Replace failures with twice the max_ratio and sort each column of r
    failures = isnan.(r)
    r[failures] = 2 * max_ratio
    r = sort(r, 1)

    (np, ns) = size(r)

    ratios = [r; 2.0 * max_ratio * ones(1, ns)]
    xs = collect(1:np+1)/np
    if length(labels) == 0
        labels = [@sprintf("column %d", col) for col = 1 : ns]
    end

    profile = Plots.plot(; kwargs...)
    for s = 1:ns
        Plots.plot!(ratios[:, s], xs, t=:steppost, label=labels[s], linewidth=linewidth)
    end
    if isfinite(xmax)
        Plots.xlims!(logscale ? 0.0 : 1.0, xmax)
    else
        Plots.xlims!(logscale ? 0.0 : 1.0, 1.1 * max_ratio)
    end
    Plots.ylims!(0, ymax)
    Plots.xlabel!("\$F\$" * (logscale ? " (log scale)" : ""))
    Plots.ylabel!("\$P\$")
    Plots.title!(title)

    return profile
end

max_stars_repo_name: chriscoey/PajaritoSupplement

max_issues_repo_name: chriscoey/PajaritoSupplement

max_forks_repo_name: chriscoey/PajaritoSupplement


text:
include("Main.jl")
include("Taskgraph.jl")
include("Architecture.jl")
include("Dump.jl")


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Computations
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, RData, CSV, XLSX, LinearAlgebra, Statistics

dir = "X:/VIVES/1-Personal/Florian/git/OECD_ICIO/src/"
dir_raw = "C:/Users/u0148308/data/raw/" # location of raw data

# other scripts
include(dir * "functions.jl") # Script with functions to import and transform raw data

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# OECD rev. 2021
source = "OECD"
revision = "2021"
year = 2018 # specified year
N = 67 # number of countries (originally MEX/CHN three times so 71)
S = 45 # number of industries

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Notes:
#   - Z: NS×NS, intermediate demand matrix
#   - Y: NS×N, final demand matrix
#   - W: NS×1, value added vector (notice: in OECD its a row vector)
#   - X: NS×1, output vector (notice: in OECD its a row vector)

Z, Y, W, X = import_data(dir_raw, source, revision, year)

Z, Y, W, X = remove_reorder(Z, Y, W, X, N, S)

any(Z .< 0.0)
any(Y .< 0.0) # possible to have negative values due to inventory!
any(W .< 0.0)
any(X .< 0.0)

any(iszero.(Z)) # will introduce NaN when used as divisor!
any(iszero.(Y)) # what to do in this case?
any(iszero.(W))
any(iszero.(X))

X_1 = [sum(Z[i,:]) + sum(Y[i,:]) for i in 1:N*S] # exports
X_2 = [sum(Z[:,i]) + W[i] for i in 1:N*S] # imports + VA
X_comp = [X_1 X_2 X X_1./X_2] # for comparison: seems to coincide better with exports computations

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# other basic matrices
A = Z ./ repeat(X', N*S) # NS×NS, input coefficients
any(isnan.(A))
A .= ifelse.(isnan.(A), 0.0, A) # in case X = 0 (otherwise we cannot compute Leontief inverse since NaN in A)

B = inv(I - A) # NS×NS, Leontief inverse
any(isnan.(B))

count(broadcast(in, B*[sum(Y[i,:]) for i in 1:N*S] ./ X, 0.95..1.05) .== 0) # gives number of values lying outside of interval

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Crosscheck computations with TiVA database: https://stats.oecd.org/Index.aspx?datasetcode=TIVA_2021_C1#
# reporter:
#   - country: AUT (n = 2)
#   - industry: D16 (wood and wood products) (s = 8)
#   - row number: (n-1)*S+s = 1*S+8 = 53
# partner:
#   - country: DEU (n = 13)
#   - industry: D41T43 (construction) (s = 25)
#   - column number: (n-1)*S+s = 12*S+25 = 565
#   - column number: (n-1)*S+s = 12*S+8 = 548 (for imports, i.e. D16 as origin industry)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.1 PROD: Production (gross output), USD million
X

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.2 VALU: Value added, USD million
W

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.3 PROD_VASH: Value added as a share of Gross Output, by industry, percentage
V = W ./ X # NS×1, (row vector in OECD)
any(isnan.(V))
V .= ifelse.(isnan.(V), 1.0, V) # in case X = 0

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.4 EXGR | EXGR_INT | EXGR_FNL: Gross exports, by industry and by partner country, USD million (f.o.b.)
#   - to obtain total => [sum(EXGR[i, :]) for i in 1:N]
#   - ctry: AUT exports to ctry: DEU => EXGR[2,13]
#   - ctry: AUT, industry: D16 exports to ctry: DEU => EXGR[53,13]
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # N×N
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.5. IMGR | IMGR_INT | IMGR_FNL: Gross imports, by industry and by partner country, USD million (f.o.b.)
#   - to obtain total => [sum(IMGR[:, i]) for i in 1:N]
#   - ctry: DEU imports from ctry: AUT => IMGR[13,2]
#   - ctry: DEU imports from ctry: AUT, industry: D16 => IMGR[13,53]
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "country", "imports") # N×N
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.6. BALGR | BALGR_INT | BALGR_FNL: Gross trade balance, by partner country, USD million (f.o.b.)
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # N×N
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "country", "imports") # N×N

#   - to obtain total => [sum(BALGR[i, :]) for i in 1:N], [sum(BALGR[:, i]) for i in 1:N]
#   - look at one row: negative/positive => import more than export (trade deficit) / export more than import (trade surplus)
#   - look at one column: positive/negative => import more than export (trade deficit) / export more than import (trade surplus)
BALGR = EXGR .- IMGR # N×N
BALGR_INT = EXGR_INT .- IMGR_INT # N×N
BALGR_FNL = EXGR_FNL .- IMGR_FNL # N×N

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.6.1 BALGR | BALGR_INT | BALGR_FNL: Gross trade balance, by partner country and by industry, USD million (f.o.b.)
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

#   - ctry: AUT, industry: D16 exports to ctry: DEU => EXGR[53,13] = EXGR[(n_r-1)*S+s, n_c]
#   - ctry: AUT imports from ctry: DEU, industry: D16 => IMGR[2,548] = IMGR[n_r, (n_c-1)*S+s]
#       + trade balance form perspective of AUT with DEU in industry D16 => BALGR[2, 13, 8] = BALGR[n_r, n_c, s]

BALGR = [EXGR[(n_r-1)*S+s, n_c] - IMGR[n_r, (n_c-1)*S+s] for n_r in 1:N, n_c in 1:N, s in 1:S]
BALGR_INT = [EXGR_INT[(n_r-1)*S+s, n_c] - IMGR_INT[n_r, (n_c-1)*S+s] for n_r in 1:N, n_c in 1:N, s in 1:S]
BALGR_FNL = [EXGR_FNL[(n_r-1)*S+s, n_c] - IMGR_FNL[n_r, (n_c-1)*S+s] for n_r in 1:N, n_c in 1:N, s in 1:S]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.7. EXGRpSH: Gross exports, partner shares, by industry, percentage
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - total exports by industry
EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S]

#   - ctry: AUT, industry: D16 export to ctry: DEU as percentage of total industry exports => EXGRpSH[53,13]
EXGRpSH = EXGR ./ repeat(EXGR_TOTAL, 1, N) .* 100
EXGRpSH_INT = EXGR_INT ./ repeat(EXGR_TOTAL, 1, N) .* 100 # what to do with NaN?
EXGRpSH_FNL = EXGR_FNL ./ repeat(EXGR_TOTAL, 1, N) .* 100

all([sum(EXGRpSH[i, :]) for i in 1:N] .≈ 100) # NaN still in there!
all([sum(EXGRpSH_INT[i, :]) .+ sum(EXGRpSH_FNL[i, :]) for i in 1:N] .≈ [sum(EXGRpSH[i, :]) for i in 1:N])

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.7.1 EXGRpSH: Gross exports, partner shares, percentage
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # NS×N

#   - total exports by industry
EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N]

#   - ctry: AUT, industry: D16 exports to ctry: DEU as percentage of total exports => EXGRpSH[2,13]
EXGRpSH = EXGR ./ repeat(EXGR_TOTAL, 1, N) .* 100
EXGRpSH_INT = EXGR_INT ./ repeat(EXGR_TOTAL, 1, N) .* 100 # what to do with NaN?
EXGRpSH_FNL = EXGR_FNL ./ repeat(EXGR_TOTAL, 1, N) .* 100

all([sum(EXGRpSH[i, :]) for i in 1:N] .≈ 100) # NaN still in there!
all([sum(EXGRpSH_INT[i, :]) .+ sum(EXGRpSH_FNL[i, :]) for i in 1:N] .≈ [sum(EXGRpSH[i, :]) for i in 1:N])

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.8. IMGRpSH: Gross imports, partner shares, by industry, percentage
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

#   - total imports
IMGR_TOTAL = [sum(IMGR[i, :]) for i in 1:N]

#   - ctry: AUT imports from ctry: DEU, industry: D16 as percentage of total imports => IMGRpSH[2,548]
#       + notice difference here to exports => percentage of total industry exports vs percentage of total imports
IMGRpSH = IMGR ./ repeat(IMGR_TOTAL, 1, N*S) .* 100
IMGRpSH_INT = IMGR_INT ./ repeat(IMGR_TOTAL, 1, N*S) .* 100
IMGRpSH_FNL = IMGR_FNL ./ repeat(IMGR_TOTAL, 1, N*S) .* 100

all([sum(IMGRpSH[i, :]) for i in 1:N] .≈ 100)
all([sum(IMGRpSH_INT[i, :]) .+ sum(IMGRpSH_FNL[i, :]) for i in 1:N] .≈ [sum(IMGRpSH[i, :]) for i in 1:N])

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 4.8.1 IMGRpSH: Gross imports, partner shares, percentage
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "country", "imports") # N×NS

#   - total imports
IMGR_TOTAL = [sum(IMGR[i, :]) for i in 1:N]

#   - ctry: AUT imports from ctry: DEU as percentage of total imports  => IMGRpSH[2,13]
IMGRpSH = IMGR ./ repeat(IMGR_TOTAL, 1, N) .* 100
IMGRpSH_INT = IMGR_INT ./ repeat(IMGR_TOTAL, 1, N) .* 100
IMGRpSH_FNL = IMGR_FNL ./ repeat(IMGR_TOTAL, 1, N) .* 100

all([sum(IMGRpSH[i, :]) for i in 1:N] .≈ 100)
all([sum(IMGRpSH_INT[i, :]) .+ sum(IMGRpSH_FNL[i, :]) for i in 1:N] .≈ [sum(IMGRpSH[i, :]) for i in 1:N])

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.1. EXGR_DVA: Domestic value added content of gross exports, USD million
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - ctry: AUT, industry: D16 domestic value added content of gross exports to ctry: DEU => EXGR_DVA[53, 13]
EXGR_DVA = fill(0.0, N*S, N) # NS×N, initialize

for c in 1:N
    v = V[(c-1)*S+1:c*S] # V_c
    b = B[(c-1)*S+1:c*S, (c-1)*S+1:c*S] # B_c,c
    exgr = EXGR[(c-1)*S+1:c*S, :] # all exports of country "c"
    for p in 1:N
        for i in 1:S
            e = fill(0.0, S) # initialize zero vector
            e[i] = exgr[i, p] # export value of country "c", industry "i" and partner country "p"
            EXGR_DVA[(c-1)*S+i, p] = v'*b*e
        end
    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.2. EXGR_DVASH: Domestic value added share of gross exports, percentage
#   - for EXGR_DVA computation look refer to 5.1.

EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - total exports by industry
EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1

#   - total domestic value added contained in gross exports by industry
EXGR_DVA_TOTAL = [sum(EXGR_DVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 share of domestic value added as percentage of total industry exports => EXGR_DVASH[53]
EXGR_DVASH = EXGR_DVA_TOTAL ./ EXGR_TOTAL .* 100 # NS×1, what to do with NaN?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.3. EXGR_TDVAIND: Industry domestic value added contribution to gross exports, as a percentage of total gross exports
#   - for EXGR_DVA computation look refer to 5.1.
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - total exports
EXGR_TOTAL = [sum(EXGR[(i-1)*S+1:i*S, :]) for i in 1:N] # N×1

#   - total domestic value added contained in gross exports by industry
EXGR_DVA_TOTAL = [sum(EXGR_DVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 share of domestic value added as percentage of total exports => EXGR_TDVAIND[53]
EXGR_TDVAIND = EXGR_DVA_TOTAL ./ repeat(EXGR_TOTAL, inner=S) .* 100 # NS×1, what to do with NaN?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.4. EXGR_DVApSH: Domestic value added in gross exports, partner shares, percentage
#   - for EXGR_DVA computation look refer to 5.1.
#   - same as 5.2. but now on the differentiated along partners
#   - [sum(EXGR_DVApSH[i, :]) for i in 1:N*S] .== EXGR_DVASH

#   - total domestic value added contained in gross exports by industry
EXGR_DVA_TOTAL = [sum(EXGR_DVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 domestic value added in exports to ctry: DEU as percentage of total industry domestic value added exports => EXGR_DVApSH[53,13]
#       + 20% of domestic value exported from AUT in D16 go to DEU
EXGR_DVApSH = EXGR_DVA ./ repeat(EXGR_DVA_TOTAL, 1, N) .* 100 # NS×N, what to do with NaN?

EXGR_DVApSH .= ifelse.(isnan.(EXGR_DVApSH), 0.0, EXGR_DVApSH)

[sum(EXGR_DVApSH[i,:]) for i in 1:N*S] .≈ 100.0

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.5. EXGR_DDC: Direct domestic industry value added content of gross exports, USD million
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

EXGR_DDC = Vector{Float64}[] # initialize
for c in 1:N
    v_hat = Diagonal(V[(c-1)*S+1:c*S])
    b_diag = Diagonal(B[(c-1)*S+1:c*S, (c-1)*S+1:c*S])
    e = [sum(EXGR[i, :]) for i in (c-1)*S+1:c*S]
    push!(EXGR_DDC, v_hat * b_diag * e)
end

#   - ctry: AUT, industry: D16 direct domestic value added content to total gross exports => EXGR_DDC[53]
EXGR_DDC = reduce(vcat, EXGR_DDC) # NS×1, collapse vector

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.6. EXGR_IDC: Indirect domestic content of gross exports (originating from domestic intermediates), USD million
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

EXGR_IDC_i = fill(0.0, N*S, S) # initialize
for c in 1:N
    v_hat = Diagonal(V[(c-1)*S+1:c*S])
    b_offdiag = B[(c-1)*S+1:c*S, (c-1)*S+1:c*S]
    for i in 1:S b_offdiag[i, i] = 0.0 end
    e_tot = [sum(EXGR[i, :]) for i in (c-1)*S+1:c*S]
    for s in 1:S
        e = fill(0.0, S)
        e[s] = e_tot[s]
        EXGR_IDC_i[(c-1)*S+s, :] = v_hat * b_offdiag * e
    end
end

#   - ctry: AUT, industry: D16 indirect domestic value added content to total gross exports from ind: D41T43 => EXGR_IDC_i[53, 25]
EXGR_IDC_i # NS×S

#   - ctry: AUT, industry: D16 indirect domestic value added content to total gross exports => EXGR_IDC[53]
EXGR_IDC = [sum(EXGR_IDC_i[i, :]) for i in 1:N*S] # NS×1

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.7. EXGR_RIM: Re-imported domestic value added content of gross exports, USD million

#   - total domestic value added in gross exports per industry
EXGR_DVA_TOTAL = [sum(EXGR_DVA[i, :]) for i in 1:N*S]

#   - ctry: AUT, industry: D16 re-imported domestic value added content to total gross exports => EXGR_RIM[53]
EXGR_RIM = EXGR_DVA_TOTAL .- EXGR_DDC .- EXGR_IDC

# wrong => no discrepancy?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.8. EXGR_FVA: Foreign value added content of gross exports, by industry, USD million
#   - EXGR_RIM is included in EXGR_FVA!
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

v_hat = Diagonal(V) # NS×NS
EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total exports per industry

EXGR_FVA_all = fill(0.0, N*S, N*S) # initialize

for c in 1:N
    for s in 1:S
        b_column = B[:, (c-1)*S+s] # all input requirements of country: c, industry: s from each other country-industry pair
        b_column[(c-1)*S+1:c*S] .= 0 # all domestic input requirements are set to zero => since captured by EXGR_DVA

        fva = v_hat * b_column * EXGR_TOTAL'
        EXGR_FVA_all[(c-1)*S+s, :] .= fva[:, (c-1)*S+s]
    end
end

#   - ctry: AUT, industry: D16 foreign value added content in gross exports => EXGR_FVA[53]
EXGR_FVA = [sum(EXGR_FVA_all[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 foreign value added content in gross exports per partner ctry: DEU => EXGR_FVA_p[53, 13]
EXGR_FVA_p = [sum(EXGR_FVA_all[i, (j-1)*S+1:j*S]) for i in 1:N*S, j in 1:N] # N×NS

any([sum(EXGR_FVA_p[i, :]) for i in 1:N*S] ≈ EXGR_FVA)

#   - ctry: AUT, ind: D16 foreign value added content in gross exports per partner ctry: DEU, ind: D41T43 => EXGR_FVA_all[53, 565]
EXGR_FVA_all # NS×NS

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.9. EXGR_FVASH: Foreign value added share of gross exports, percentage
#   - EXGR_RIM is included in EXGR_FVA!
#   - see 5.8. for EXGR_FVA_all computation
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total exports per industry

#   - ctry: AUT percentage share of foreign value added content in gross exports in ind: D16 => EXGR_FVASH[53]
EXGR_FVASH = EXGR_FVA ./ EXGR_TOTAL .* 100 # NS×1, what to do with NaN?

#   - ctry: AUT, ind: D16 percentage points of foreign value added content in gross exports by partner ctry: DEU => EXGR_FVApSH[53, 13]
#       + EXGR_FVApSH[53, 13] percentage points of foreign value added in ctry: AUT, ind: D16 originates in ctry: DEU
#       + EXGR_FVASH[53] = 30% of VA comes is foreign, EXGR_FVApSH[53, 13] = 8 p.p. originate in DEU
#           i.e. EXGR_FVApSH[53, 13]/EXGR_FVASH[53] = 8/30 = 27% of FVA in ctry: AUT, ind: D16 comes from DEU
EXGR_FVApSH = EXGR_FVA_p ./ repeat(EXGR_TOTAL, 1, N) .* 100

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.10. EXGR_TFVAIND Industry foreign value added contribution to gross exports, a as a percentage of total gross exports
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N] # N×1, total exports

#   - ctry: AUT, ind: D16 percentage share of foreign value added content in total gross exports => EXGR_TFVAIND[53]
#       + EXGR_TFVAIND[53] = 0.67 p.p. of foreign value added contained in AUT exports are contained in industry D16
#           i.e. to compare across industries
EXGR_TFVAIND = EXGR_FVA ./ repeat(EXGR_TOTAL, inner=S) .* 100

#   - ctry: AUT percentage share of foreign value added content in total gross exports => EXGR_TFVAIND_c[2]
EXGR_TFVAIND_c = [sum(EXGR_TFVAIND[(i-1)*S+1:i*S]) for i in 1:N]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.1. EXGR_BSCI: Origin of value added in gross exports, USD million
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total exports per industry
v_hat = Diagonal(V) # NS×NS, diagonalized value added coefficients

EXGR_BSCI = fill(0.0, N*S, N*S) # initialize
for c in 1:N
    for s in 1:S
        e = fill(0.0, N*S)
        e[(c-1)*S+s] = EXGR_TOTAL[(c-1)*S+s]

        EXGR_BSCI[(c-1)*S+s, :] = v_hat * B * e
    end
end

#   - ctry: AUT, ind: D16 value added in gross exports originating from ctry: DEU, ind: D41T43 => EXGR_BSCI[53, 565]
EXGR_BSCI

#   - EXGR_BSCI gathers all value added information (why not compute other statistics from there?)
#       + domestic value added in ctry: AUT, ind: D16 => sum(EXGR_BSCI[53, 45:90]) ≈ sum(EXGR_DVA[53,:])
#       + direct domestic value added in ctry: AUT, ind: D16 => EXGR_BSCI[53, 53] == EXGR_DDC[53]
#       + indirect domestic value added in ctry: AUT, ind: D16 => sum(EXGR_BSCI[53, [46:52; 54:90]]) == EXGR_IDC[53]
#       + foreign value added in ctry: AUT, ind: D16 => sum(EXGR_BSCI[53, [1:45; 91:N*S]]) == EXGR_FVA[53]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.11. DEXFVApSH: Backward participation in GVCs, percentage
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N] # N×1, total exports

#   - ctry: AUT value added in total gross exports from ctry: DEU => EXGR_BSCI_c[2, 13]
EXGR_BSCI_c = [sum(EXGR_BSCI[(c-1)*S+1:c*S, (p-1)*S+1:p*S]) for c in 1:N, p in 1:N] # N×N

#   - ctry: AUT percentage share of value added in total gross exports originating from ctry: DEU => DEXFVApSH[2, 13]
#       + DEXFVApSH[2, 13] = 8.3% of value added exports from AUT originates from DEU
DEXFVApSH = EXGR_BSCI_c ./ repeat(EXGR_TOTAL, 1, N) .* 100

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.12. EXGR_DVAFXSH: Domestic value added embodied in foreign exports as share of gross exports, percentage
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N] # N×1, total exports

#   - ctry: AUT, ind: D16 value added in total gross exports from ctry: DEU => EXGR_BSCI_c_i[53, 13]
EXGR_BSCI_c_i = [sum(EXGR_BSCI[i, (p-1)*S+1:p*S]) for i in 1:N*S, p in 1:N] # NS×N

EXGR_DVAFXSH = Float64[] # initialize
for c in 1:N
    va_c = EXGR_BSCI_c_i[:, c] # column gives domestic VA per ctry c in contained in row country-industry
    va_c[(c-1)*S+1:c*S] .= 0.0 # set domestic VA to c's own exports to zero
    for s in 1:S
        va_c_i = va_c[s:S:N*S] # subset to industry i
        share = sum(va_c_i) / EXGR_TOTAL[c] * 100 # sum for total (i.e. WORLD as partner) and divide by c's total exports
        push!(EXGR_DVAFXSH, share)
    end
end

#   - VA originating from ctry: AUT, ind: D16 as percentage share of total exports => EXGR_DVAFXSH[53]
#       + AUT domestic value added content embodied in the gross exports of industry D16 in foreign countries
#           as a percentage of total gross exports of country AUT

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.13. FEXDVApSH: Forward participation in GVCs, percentage
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "country", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N] # N×1, total exports

#   - ctry: AUT value added in total gross exports from ctry: DEU => EXGR_BSCI_c[2, 13]
EXGR_BSCI_c = [sum(EXGR_BSCI[(c-1)*S+1:c*S, (p-1)*S+1:p*S]) for c in 1:N, p in 1:N] # N×N

#   - ctry: AUT percentage share of domestic value added in total gross exports of ctry: DEU => FEXDVApSH[2, 13]
#       + FEXDVApSH[2, 13] = 1.2% of DEU exports is domestic value added originating in AUT
FEXDVApSH = EXGR_BSCI_c ./ repeat(EXGR_TOTAL', N) .* 100

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.14. EXGR_INTDVASH: Domestic value added in exports of intermediate products, as a share of total gross exports, percentage
#   - same calculation as for 5.1. EXGR_DVA
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - ctry: AUT, industry: D16 domestic value added content of gross exports in intermediate input to ctry: DEU => EXGR_INTDVA[53, 13]
EXGR_INTDVA = fill(0.0, N*S, N) # NS×N, initialize

for c in 1:N
    v = V[(c-1)*S+1:c*S] # V_c
    b = B[(c-1)*S+1:c*S, (c-1)*S+1:c*S] # B_c,c
    exgr = EXGR_INT[(c-1)*S+1:c*S, :] # all exports of country "c"
    for p in 1:N
        for i in 1:S
            e = fill(0.0, S) # initialize zero vector
            e[i] = exgr[i, p] # export value of country "c", industry "i" and partner country "p"
            EXGR_INTDVA[(c-1)*S+i, p] = v'*b*e
        end
    end
end

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total industry exports

#   - total domestic value added contained in gross exports intermediate inputs by industry
EXGR_INTDVA_TOTAL = [sum(EXGR_INTDVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 share of domestic value added in gross exports of intermediate inputs as percentage of total industry exports => EXGR_INTDVASH[53]
EXGR_INTDVASH = EXGR_INTDVA_TOTAL ./ EXGR_TOTAL .* 100 # NS×1, what to do with NaN?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.15. EXGR_FNLDVASH: Domestic value added in exports of final products, as a share of total gross exports, percentage
#   - same calculation as for 5.1. EXGR_DVA and 5.14. EXGR_INTDVA
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

#   - ctry: AUT, industry: D16 domestic value added content of gross final demand exports to ctry: DEU => EXGR_FNLDVA[53, 13]
EXGR_FNLDVA = fill(0.0, N*S, N) # NS×N, initialize

for c in 1:N
    v = V[(c-1)*S+1:c*S] # V_c
    b = B[(c-1)*S+1:c*S, (c-1)*S+1:c*S] # B_c,c
    exgr = EXGR_FNL[(c-1)*S+1:c*S, :] # all exports of country "c"
    for p in 1:N
        for i in 1:S
            e = fill(0.0, S) # initialize zero vector
            e[i] = exgr[i, p] # export value of country "c", industry "i" and partner country "p"
            EXGR_FNLDVA[(c-1)*S+i, p] = v'*b*e
        end
    end
end

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total industry exports

#   - total domestic value added contained in gross final demand exports by industry
EXGR_FNLDVA_TOTAL = [sum(EXGR_FNLDVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 share of domestic value added in gross exports of final demand as percentage of total industry exports => EXGR_FNLDVASH[53]
EXGR_FNLDVASH = EXGR_FNLDVA_TOTAL ./ EXGR_TOTAL .* 100 # NS×1, what to do with NaN?

any(EXGR_INTDVASH .+ EXGR_FNLDVASH .== EXGR_DVASH)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.16. EXGR_INTDVApSH: Domestic value added in exports of intermediate products, partner shares, percentage

#   - total domestic value added contained in gross exports by industry
EXGR_INTDVA_TOTAL = [sum(EXGR_INTDVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 domestic value added in gross exports of intermediate inputs to ctry: DEU
#           as percentage of total industry domestic value added exports of intermediate inputs => EXGR_INTDVApSH[53,13]
#       + EXGR_INTDVApSH[53,13] = 19% of domestic value added of intermediate input exports from AUT in D16 go to DEU
EXGR_INTDVApSH = EXGR_INTDVA ./ repeat(EXGR_INTDVA_TOTAL, 1, N) .* 100 # NS×N, what to do with NaN?

EXGR_INTDVApSH .= ifelse.(isnan.(EXGR_INTDVApSH), 0.0, EXGR_INTDVApSH)
count([sum(EXGR_INTDVApSH[i,:]) for i in 1:N*S] .≈ 100.0)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.16.1 EXGR_FNLDVApSH: Domestic value added in exports of final products, partner shares, percentage

#   - total domestic value added contained in gross exports by industry
EXGR_FNLDVA_TOTAL = [sum(EXGR_FNLDVA[i, :]) for i in 1:N*S] # NS×1

#   - ctry: AUT, industry: D16 domestic value added in gross exports of final products to ctry: DEU
#           as percentage of total industry domestic value added exports of final products => EXGR_FNLDVApSH[53,13]
#       + EXGR_FNLDVApSH[53,13] = 38% of domestic value added of final product exports from AUT in D16 go to DEU
EXGR_FNLDVApSH = EXGR_FNLDVA ./ repeat(EXGR_FNLDVA_TOTAL, 1, N) .* 100 # NS×N, what to do with NaN?

EXGR_FNLDVApSH .= ifelse.(isnan.(EXGR_FNLDVApSH), 0.0, EXGR_FNLDVApSH)
count([sum(EXGR_FNLDVApSH[i,:]) for i in 1:N*S] .≈ 100.0)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.19. IMGR_DVA: Domestic value added embodied in gross imports, USD million
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

#   - ctry: AUT domestic value added content of gross imports from ctry: DEU, ind: D16 => IMGR_DVA[2, 548]
IMGR_DVA = fill(0.0, N, N*S) # NS×N, initialize

for c in 1:N
    v = V[(c-1)*S+1:c*S] # V_c
    for p in 1:N
        b = B[(c-1)*S+1:c*S, (p-1)*S+1:p*S] # B_c,p
        imgr = Diagonal(IMGR[c, (p-1)*S+1:p*S])
        IMGR_DVA[c, (p-1)*S+1:p*S] = v'*b*imgr
    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.20. IMGR_DVASH: Domestic value added share of gross imports, percentage
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

#   - ctry: AUT share of domestic value added in gross imports from ctry: DEU, ind: 16 => IMGR_DVASH[2, 548]
IMGR_DVASH = IMGR_DVA ./ IMGR .* 100 # N×NS, what to do with NaN? (i.e. domestic va imported from own ctry)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.21. REII: Re-exported intermediate imports, USD million
EXGR, EXGR_INT, EXGR_FNL = trade(Z, Y, N, S, "industry", "exports") # NS×N

EXGR_TOTAL = [sum(EXGR[i, :]) for i in 1:N*S] # NS×1, total industry exports

REII_p = fill(0.0, N*S, N)

for c in 1:N
    for p in 1:N
        a = A[(p-1)*S+1:p*S, (c-1)*S+1:c*S]
        #for i in 1:S a[i, i] = 0 end # written in TiVA guide but not true
        e = EXGR_TOTAL[(c-1)*S+1:c*S]
        b = B[(c-1)*S+1:c*S, (c-1)*S+1:c*S]
        REII_p[(c-1)*S+1:c*S, p] = a*b*e
    end
end

#   - Total intermediate products imported in ctry: AUT, ind: D16 originated from ctry: DEU => REII_p[53, 13]
REII_p # NS×N

REII = copy(REII_p)
for c in 1:N*S
    REII[c, ceil(Int, c/S)] = 0.0 # set all domestic contributions to zero
end

#   - Total intermediate products absorbed in ctry: AUT, ind: D16 originated from all foreign countries => REII_p[53]
#       + sum(REII_p[53, :]) - REII_p[53, 2] ≈ REII[53]
REII = [sum(REII[i, :]) for i in 1:N*S] # NS×1

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 5.22. IMGRINT_REII: Re-exported intermediate imports as a % of total intermediate imports, percentage
#   - ctry: DEU imports from ctry: AUT, industry: D16 => IMGR[13,53]
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS

IMGR_INT_TOTAL = [sum(IMGR_INT[i, s:S:N*S]) for i in 1:N, s in 1:S] # N×S, total intermediate industry imports

#   - ctry: AUT, ind: D16 re-exported intermediate imports as a share of total intermediate imports from all foreign countries => IMGRINT_REII[53]
IMGRINT_REII = REII ./ reshape(IMGR_INT_TOTAL', N*S) .* 100 # NS×1

#   - ctry: AUT, ind: D16 re-exported intermediate imports as a share of total intermediate imports from ctry: DEU => IMGRINT_REII_p[53, 13]
IMGRINT_REII_p = REII_p ./ IMGR_INT' .* 100 # NS×N, what to do with Inf?

IMGRINT_REII_p .= ifelse.(isinf.(IMGRINT_REII_p), 0.0, IMGRINT_REII_p)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.1. FFD_DVA: Domestic value added embodied in foreign final demand, USD million
v_hat = Diagonal(V)
FD = copy(Y)

#   - Value added originating in ctry:AUT, ind: D16 embodied in final demand of ctry: DEU => FFD_DVA[53, 13]
FFD_DVA = v_hat * B * FD # NS×N

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.2. FFD_DVApSH: Domestic value added embodied in foreign final demand, partner shares, percentage
FFD_DVA_c = copy(FFD_DVA)
for c in 1:N*S
    FFD_DVA_c[c, ceil(Int, c/S)] = 0.0 # do not count domestic value added to domestic FD in total
end
FFD_DVA_TOTAL = [sum(FFD_DVA_c[i, :]) for i in 1:N*S] # NS×1

#   - Value added originating in ctry:AUT, ind: D16 embodied in final demand of ctry: DEU
#       as percentage of total domestic value added originating in ctry: AUT, ind: D16 demanded abroad => FFD_DVApSH[53, 13]
FFD_DVApSH = FFD_DVA_c ./ repeat(FFD_DVA_TOTAL, 1, N) .* 100 # NS×N, what to do with NaN?

count([sum(FFD_DVApSH[i, :]) for i in 1:N*S] .≈ 100)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.3. VALU_FFDDVA: Domestic value added embodied in foreign final demand as a % of total value added, percentage
VALU = copy(W) # NS×1, industry value added

FFD_DVA_c = copy(FFD_DVA)
for c in 1:N*S
    FFD_DVA_c[c, ceil(Int, c/S)] = 0.0 # do not count domestic value added to domestic FD in total
end

#   - Domestic value added from ctry: AUT, ind: D16 embodied in foreign final demand in ctry: DEU
#       as percentage of total value added of ctry: AUT, ind: D16 => VALU_FFDDVA_p[53, 13]
VALU_FFDDVA_p = FFD_DVA_c ./ repeat(VALU, 1, N) .* 100

#   - Domestic value added from ctry: AUT, ind: D16 embodied in total foreign final demand
#       as percentage of total value added of ctry: AUT, ind: D16 => VALU_FFDDVA[53]
VALU_FFDDVA = [sum(VALU_FFDDVA_p[i, :]) for i in 1:N*S]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.4. DFD_FVA: Foreign value added embodied in domestic final demand, USD million

#   - Value added originating from ctry: AUT, ind: D16 embodied in final demand of ctry: DEU => DFD_FVA[13, 53]
DFD_FVA = FFD_DVA'

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.5. DFD_FVApSH: Foreign value added embodied in domestic final demand, partner shares, percentage

DFD_FVA_TOTAL = copy(DFD_FVA)
for p in 1:N*S
    DFD_FVA_TOTAL[ceil(Int, p/S), p] = 0.0 # do not count domestic value added
end
DFD_FVA_TOTAL = [sum(DFD_FVA_TOTAL[i, s:S:N*S]) for i in 1:N, s in 1:S]

#   - Foreign value added originating from ctry: AUT, ind: D16 embodied in final demand of ctry: DEU
#       as percentage of total foreign value added of ind: D16 => DFD_FVApSH[13, 53]
DFD_FVApSH = DFD_FVA ./ repeat(DFD_FVA_TOTAL, 1, N) .* 100
for p in 1:N*S
    DFD_FVApSH[ceil(Int, p/S), p] = 0.0 # set domestic VA to zero
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.6. BALVAFD: Value added embodied in final demand, balance, USD million
BALVAFD_all = [FFD_DVA[i, ceil(Int, j/S)] - DFD_FVA[ceil(Int, i/S), j] for i in 1:N*S, j in 1:N*S] # NS×NS, balance of all ctry-ind pairs

#   - Domestic value added embodied in foreign final demand minus foreign value added in domestic final demand by the value added industry
#   - Balance from the point of ctry: AUT with ctry: DEU in ind: D16 => BALVAFD[53, 13]
#   - Balance from the point of ctry: DEU with ctry: AUT in ind: D16 => BALVAFD[548, 2]
BALVAFD = fill(0.0, N*S, N)
for c in 1:N
    for s in 1:S
        BALVAFD[(c-1)*S+s, :] = BALVAFD_all[(c-1)*S+s, s:S:N*S] # NS×N, filter only same industries (i.e. balance between D16 and D41T43 doesnt make sense)
    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

#6.7. FD_VA | CONS_VA | GFCF_VA: Value added embodied in final demand, consumption and GFCF, USD million
v_hat = Diagonal(V)
FD = copy(Y)

#   - Value added originating in ctry: AUT, ind: D16 embodied in final demand of ctry: DEU => FD_VA[53, 13]
#       + same as FFD_DVA but we do not take out domestic VA later
FD_VA = v_hat * B * FD # NS×N

#   - Value added originating in ctry: AUT, ind: D16 embodied in consumption of ctry: DEU => CONS_VA[53, 13]
#       + Consumption (CONS) is made out of household consumption, government expenditure and NPISHs (subset of FD)
# CONS_VA = v_hat * B * CONS

#   - Value added originating in ctry: AUT, ind: D16 embodied in gross fixed capital formation of ctry: DEU => GFCF_VA[53, 13]
#       + Gross fixed capital formation (GFCF) covers demand for private and government investment (subset of FD)
# GFCF_VA = v_hat * B * GFCF

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 6.8. FD_VASH | CONS_VASH | GFCF_VASH: Value added shares in final demand, consumption and GFCF, percentage
FD_VA_TOTAL = [sum(FD_VA[s:S:N*S, j]) for s in 1:S, j in 1:N] # S×N, total VA embodied in final demand per industry and destination ctry

#   - Value added originating in ctry: AUT, ind: D16 embodied in final demand of ctry: DEU
#       as percentage share of total value added used in ctry: DEU, ind: D16 => FD_VASH[53, 13]
FD_VASH = FD_VA ./ repeat(FD_VA_TOTAL, N) .* 100 # NS×N, what to do with NaN?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.2. IMGR_BSCI: Origin of value added in gross imports, USD million
IMGR, IMGR_INT, IMGR_FNL = trade(Z, Y, N, S, "industry", "imports") # N×NS
v_hat = Diagonal(V)

IMGR_BSCI = fill(0.0, N*S, N*S, N) # NS×NS×N
for c in 1:N
    IMGR_BSCI[:, :, c] = v_hat * B * Diagonal(IMGR[c, :])
end

#   - Value added in gross imports in ctry: AUT from ctry: DEU, ind: D16 where VA originated in source ctry: USA, ind: D29  => IMGR_BSCI[1685, 548, 2]
#       + rows: source country and industry, columns: exporing country and industry, 3rd dimension: importing ctry
IMGR_BSCI

#   - IMGR_BSCI gathers most (all?) value added information (why not compute other statistics from there?)
#       + total gross imports by ctry: AUT => sum(IMGR_BSCI[:, :, 2]) ≈ sum(IMGR[2,:])
#       + total gross imports by ctry: AUT per industry => [sum(IMGR_BSCI[:, j, 2]) for j in 1:N*S] ≈ IMGR[2, :]
#       + domestic value added by ctry: AUT importied via ctry: DEU, ind: D16 => sum(IMGR_BSCI[46:90, 548, 2]) ≈ IMGR_DVA[2, 548]

# Note:
#   - Note that the same value added originating from source country s can be present in the gross imports of more
#       than one importing country c (as embodied value added, from upstream production, may cross national borders many times).
#       In general, therefore, these estimates should be viewed from the perspective of an importing country c.
#   => reason why the sums are usually slightly higher than when using the other method?

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.3. FDVA_BSCI: Origin of value added in final demand, USD million
FD = copy(Y)
FD_TOTAL = [sum(FD[i, :]) for i in 1:N*S]

#   - Value added entailed in final demand of ctry: AUT, ind: D16 generated by ctry: DEU, ind: D41T43 => FDVA_BSCI_2[53, 553]
FDVA_BSCI = zeros(N*S, N*S)
for ph in 1:N*S
    fd = zeros(N*S)
    fd[ph] = FD_TOTAL[ph]
    FDVA_BSCI[ph, :] = v_hat * B * fd # NS×NS
end

c = 2
s = 8
fd = zeros(N*S)
fd[(c-1)*S+s] = FD_TOTAL[(c-1)*S+s]
e = v_hat * B * fd
sum(e)

sum(e[541:585])

sum(FDVA_BSCI_2[46:90, 548])

# WRONG!

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.4. FD_EXGRINT_VA: Gross exports of intermediate products by origin of value added and final destination, USD million
FD = copy(Y)
FD_TOTAL = [sum(FD[i, :]) for i in 1:N*S]

γ = zeros(N*S, N*S)
for c in 1:N
    γ[(c-1)*S+1:c*S, (c-1)*S+1:c*S] .= 1.0
    for s in 1:S
        γ[(c-1)*S+s, (c-1)*S+s] = 0.0
    end
end

F = Diagonal( ( γ .*  A * Diagonal(B * FD_TOTAL) ) * ones(N*S) )

# NO IDEA

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.5. FD_EXGRFNL_VA: Gross exports of final products by origin of value added and final destination, USD million

# NO IDEA

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# 7.6. FD_EXGR_VA: Gross exports by origin of value added and final destination, USD million

# NO IDEA

# -------------------------------------------------------------------------------------------------------------------------------------------------------------


text:
using Statistics

include("update.jl")

function randomwalk(num::Integer, nsteps::Integer, prob::Real=0.5)
    means, vars = zeros(nsteps+1), zeros(nsteps+1)

    xs = zeros(num)
    means[1] = mean(xs)
    vars[1] = var(xs)
    for t in 1:nsteps
        map!(x -> update(x, prob), xs, xs)
        means[t+1] = mean(xs)
        vars[t+1] = var(xs)
    end
    return means, vars
end

max_stars_repo_name: isuke/JuliaBook-Samples

max_issues_repo_name: isuke/JuliaBook-Samples

max_forks_repo_name: isuke/JuliaBook-Samples


text:
using Base64
using Serialization

OVERRIDE_MSG = "This should be overridden in implemention of AbstractEnvState"
num_rows(::Type{<:AbstractEnvState})              = @assert false OVERRIDE_MSG
num_cols(::Type{<:AbstractEnvState})              = @assert false OVERRIDE_MSG
player1_markers(::Type{<:AbstractEnvState})       = @assert false OVERRIDE_MSG
player2_markers(::Type{<:AbstractEnvState})       = @assert false OVERRIDE_MSG
get_pieces(::Type{<:AbstractEnvState})            = @assert false OVERRIDE_MSG
image_scale(::Type{<:AbstractEnvState})           = @assert false OVERRIDE_MSG
image_title_offset(::Type{<:AbstractEnvState})    = @assert false OVERRIDE_MSG
get_winner(::AbstractEnvState)                    = @assert false OVERRIDE_MSG
actions(::Type{<:AbstractEnvState})               = @assert false OVERRIDE_MSG
num_players(::Type{<:AbstractEnvState})           = @assert false OVERRIDE_MSG

key(s::AbstractEnvState) = encode(s)
svg(s::AbstractEnvState) = BoardStateSVG(s)
image_scale(x::AbstractEnvState) = image_scale(typeof(x))
image_title_offset(x::AbstractEnvState) = image_title_offset(typeof(x))
player1_markers(x::AbstractEnvState) = player1_markers(typeof(x))
player2_markers(x::AbstractEnvState) = player2_markers(typeof(x))
get_pieces(x::AbstractEnvState) = get_pieces(typeof(x))
num_rows(x::AbstractEnvState) = num_rows(typeof(x))
num_cols(x::AbstractEnvState) = num_cols(typeof(x))
num_cells(T::Type{<:AbstractEnvState}) = GridGames.getnumrows(T)*GridGames.getnumcols(T)
num_cells(x::AbstractEnvState) = num_cells(typeof(x))
Base.getindex(x::AbstractEnvState, idx::Integer) = x.board_state[idx]
Base.iterate(x::AbstractEnvState, state=1) = iterate(x.board_state, state)
action_space_size(T::Type{<:AbstractEnvState}) = getactions(T) |> length
action_space_size(x::AbstractEnvState) = action_space_size(typeof(x))
piece(s::AbstractEnvState, name::Symbol) = piece(typeof(s), name)
piece(::AbstractEnvState, ::Nothing) = nothing
piece(::Type{<:AbstractEnvState}, ::Nothing) = nothing
game_is_won(s::AbstractEnvState) = !isnothing(get_winner(s))
actions(s::AbstractEnvState) = actions(typeof(s))
to_play(s::AbstractEnvState) = s.to_play
dynamics(s::AbstractEnvState, a::Integer) = get_reward(s, a), next_state(s, a)
get_action_image(s::AbstractEnvState, a::Integer) = get_action_image(typeof(s), a)

function get_action_image(T::Type{<:AbstractEnvState}, a::Integer)
    channel = reshape(Flux.onehot(a, 1:num_cells(T)), GridGames.getnumrows(T), GridGames.getnumcols(T))
    Image(convert(Array{Float32}, Flux.unsqueeze(channel, 3)))
end

# function next_to_play(s::AbstractEnvState)
#     if num_players(typeof(s)) == 1
#         return PLAYER_1
#     end
#     to_play(s) == PLAYER_1 ? PLAYER_2 : PLAYER_1
# end

function get_image(s::AbstractEnvState)
    #pieces = piece_matrix(s)
    pieces = getpiecematrix(s)
    piece_channels = map(piece_idx->begin
        pieces .== piece_idx
    end, 1:length(getpieces(s)))
    to_play_multiplier = to_play(s).idx == PLAYER_1 ? Float32(1) : Float32(-1)
    to_play_channel = to_play_multiplier*ones(Float32, GridGames.getnumrows(s), GridGames.getnumcols(s))
    channels = reduce((a,b)->cat(a, b; dims=3), piece_channels; init=to_play_channel)
    Image(channels)
end
function get_reward(s::AbstractEnvState, action::Integer)
    winner = next_state(s, action) |> get_winner
    if isnothing(winner)
        return Float32(0)
    end
    winner == to_play(s) ? Float32(1) : Float32(-1)
end
function next_state(s::AbstractEnvState, action::Integer)
    typeof(s)(
        next_board_state(s, action),
        next_to_play(s)
    )
end
# function next_board_state(s::AbstractEnvState, action::Integer)
#     marker = get_marker(s, action)
#     coord = marker |> get_coordinate
#     byte_idx = coord_to_byte_idx(s, coord)
#     new_cell_byte = UInt8(get_piece_idx(s, action))
#     BoardState(s.board_state, byte_idx, new_cell_byte)
# end
get_marker(s::AbstractEnvState, a::Integer) = get_marker(s, to_play(s), a)
get_piece(s::AbstractEnvState, a::AbstractAction) = get_marker(s, a.i) |> get_piece

function get_piece_idx(s::AbstractEnvState, action::Integer)
    marker = get_marker(s, action)
    if isnothing(marker)
        return nothing
    end
    get_piece_idx(typeof(s), marker.piece.name)
end

function piece(T::Type{<:AbstractEnvState}, name::Symbol)
    get_pieces(T)[get_piece_idx(T, name)]
end

function get_piece(s::AbstractEnvState, name::Symbol)
    get_pieces(s)[get_piece_idx(s, name)]
end

function get_piece(s::AbstractEnvState, cell_idx::Integer)
    byte = s.board_state[cell_idx]
    if iszero(byte)
        return nothing
    end
    get_pieces(s)[byte]
end

function piece_count(s::AbstractEnvState, symbol)
    filter(byte->begin
        if iszero(byte)
            return false
        end
        get_pieces(s)[byte] == get_piece(s, symbol)
    end, s.board_state.bytes) |> length
end

get_piece_idx(s::AbstractEnvState, name::Symbol) = get_piece_idx(typeof(s), name)
function get_piece_idx(T::Type{<:AbstractEnvState}, name::Symbol)
    findfirst(x->x.name == name, get_pieces(T))
end

# @todo Remove these functions
hex_str(x::AbstractEnvState) = encode(x)[1:7]
#

# @todo Remove this. Needs to be at the game level to support
# more flexible game play outcomes.
function last_reward(s::AbstractEnvState)
    isnothing(get_winner(s)) ? Float32(0) : Float32(1)
end

function Base.show(io::IO, s::AbstractEnvState)
    encoded = encode(s)
    if length(encoded) >= 40
        encoded = join([first(encoded, 20), "...", last(encoded, 20)])
    end
    print(io, "$(string(typeof(s))){\"$encoded\"}")
end

function Base.:-(s1::AbstractEnvState, s2::AbstractEnvState)
    bs1 = s1.board_state
    bs2 = s2.board_state
    byte_indices = findall(x->x > 0, bs1 - bs2)
    byte_idx_to_coord.(typeof(s1), byte_indices)
end

function byte_idx_to_coord(T::Type{<:AbstractEnvState}, byte_idx::Integer)
    row_idx = mod1(byte_idx, num_rows(T))
    col_idx = ((byte_idx - 1) ÷ num_rows(T)) + 1
    row_idx, col_idx
end

function coord_to_byte_idx(s::AbstractEnvState, coord::Tuple{<:Integer, <:Integer})
    coord_to_byte_idx(typeof(s), coord)
end

function coord_to_byte_idx(T::Type{<:AbstractEnvState}, coord::Tuple{<:Integer, <:Integer})
    row_idx, col_idx = coord
    (col_idx - 1)*num_rows(T) + row_idx
end

# function Base.getindex(s::AbstractEnvState, row_idx::Integer, col_idx::Integer)
#     byte_idx = coord_to_byte_idx(s, (row_idx, col_idx))
#     s[byte_idx]
# end

# function Base.getindex(s::AbstractEnvState, row_idx::Integer, ::Colon)
#     s[row_idx, 1:num_cols(s)]
# end

# function Base.getindex(s::AbstractEnvState, row_idx::Integer, col_range::UnitRange{<:Integer})
#     start_byte_idx = coord_to_byte_idx(s, (row_idx, col_range[1]))
#     end_byte_idx = coord_to_byte_idx(s, (row_idx, col_range[end]))
#     by = num_rows(s)
#     [s[idx] for idx in start_byte_idx:by:end_byte_idx]
# end

# function Base.getindex(s::AbstractEnvState, ::Colon, col_idx::Integer)
#     s[1:num_rows(s), col_idx]
# end

# @todo row_range should support step size (the intermediate value) as well
# function Base.getindex(s::AbstractEnvState, row_range::UnitRange{<:Integer}, col_idx::Integer)
#     start_byte_idx = coord_to_byte_idx(s, (row_range[1], col_idx))
#     end_byte_idx = coord_to_byte_idx(s, (row_range[end], col_idx))
#     [s[idx] for idx in start_byte_idx:end_byte_idx]
# end

# @todo
# function Base.getindex(board_state::AbstractEnvState, ::Colon, col_idx::Integer)

#Base.copy(s::AbstractEnvState) = typeof(s)(copy(g.board_state))

function Base.show(io::IO, m::MIME"text/plain", s::AbstractEnvState)
    function cell_fn(_,_,row,col)
        piece_idx = s[row,col]
        piece_idx == 0 ? nothing : plaintext(get_pieces(s)[piece_idx])
    end
    g = Grid(num_cols(s), num_rows(s))
    show(io, "text/plain", g(cell_fn))
end

# function piece_matrix(s::AbstractEnvState)
#     reshape(
#         collect(s.board_state.bytes),
#         num_rows(s),
#         num_cols(s))
# end

function diff_coord(s1::AbstractEnvState, s2::AbstractEnvState)
    diff = (s1 - s2)
    isempty(diff) ? nothing : diff[1]
end

function targetsvg(input_state::AbstractEnvState, target_state::AbstractEnvState, target::Target)
    r = @sprintf "%.2f" reward(target)
    v = @sprintf "%.2f" value(target)
    highlight_coord = diff_coord(target_state, input_state)
    gridsvg(
        target_state;
        title="target",
        title_offset=(-210,-195),
        header="r: $r",
        header_color="#333",
        br_header="v: $v",
        br_header_color="#333",
        highlight=highlight_coord
    )
end

function gridsvg(
    s::AbstractEnvState;
    highlight=nothing,
    highlight_color="#aacef0",
    main_title="state",
    aspect_ratio=(num_cols(s), num_rows(s)),
    scale=image_scale(s),
    main_title_offset=image_title_offset(s),
    kwargs...)
    function cell_fn(_,_,row,col)
        piece_idx = s[row,col]
        if piece_idx == 0
            return nothing
        end
        p = get_pieces(s)[piece_idx]
        bg_color = highlight == (row,col) ? highlight_color : nothing
        (value=p.icon_symbol, color=p.icon_color, size=p.icon_size, bg_color=bg_color)
    end
    GridSVG(
        Grid(num_cols(s), num_rows(s), cell_fn);
        main_title=main_title,
        main_title_offset=main_title_offset,
        aspect_ratio=aspect_ratio,
        scale=scale,
        kwargs...
    )
end

function Base.show(io::IO, m::MIME"image/svg+xml", s::AbstractBoardStateSVG)
    show(io, m, gridsvg(s.x))
end


text:
export SCIPSolver

# Linear Quadratic Model

type SCIPMathProgModel <: AbstractLinearQuadraticModel
    ptr_model::Ptr{Void}
    options
    lazy_userdata
    heur_userdata

    function SCIPMathProgModel(options...)
        _arr = Array(Ptr{Void}, 1)
        # TODO: check return code (everywhere!)
        ccall((:CSIPcreateModel, libcsip), Cint, (Ptr{Ptr{Void}}, ), _arr)
        m = new(_arr[1], options)
        # QUESTION: Why is _arr not garbage-collected?
        setparams!(m)
        return m
    end
end

# Nonlinear Model

type SCIPNonlinearModel <: AbstractNonlinearModel
    m::SCIPMathProgModel
end

# Solver

immutable SCIPSolver <: AbstractMathProgSolver
    options
end
SCIPSolver(kwargs...) = SCIPSolver(kwargs)
LinearQuadraticModel(s::SCIPSolver) = SCIPMathProgModel(s.options)
NonlinearModel(s::SCIPSolver) = SCIPMathProgModel(s.options)


text:
# The zlib Compression
# ====================

function compress!(dst::Vector{UInt8}, src::Vector{UInt8})::UInt64
    sz = Ref(Culong(length(dst)))
    code = ccall(
        (:compress, Libz.zlib),
        Cint,
        (Ptr{Cvoid}, Ref{Culong}, Ptr{Cvoid}, Culong),
        dst, sz, src, sizeof(src))
    if code != Libz.Z_OK
        Libz.zerror(code)
    end
    return sz[]
end

function uncompress!(dst::Vector{UInt8}, src::Vector{UInt8})::UInt64
    sz = Ref(Culong(length(dst)))
    code = ccall(
        (:uncompress, Libz.zlib),
        Cint,
        (Ptr{Cvoid}, Ref{Culong}, Ptr{Cvoid}, Culong),
        dst, sz, src, sizeof(src))
    if code != Libz.Z_OK
        Libz.zerror(code)
    end
    return sz[]
end


text:
module operators

export Operator, length, basis, dagger, ishermitian, tensor, embed,
        trace, ptrace, normalize, normalize!, expect, variance,
        expm, permutesystems, identityoperator

import Base: ==, +, -, *, /, length, trace, one, ishermitian, expm
import ..bases: basis, tensor, ptrace, permutesystems,
            samebases, check_samebases, multiplicable
import ..states: dagger, normalize, normalize!

using Compat
using ..sortedindices, ..bases, ..states

"""
Abstract base class for all operators.

All deriving operator classes have to define the fields
`basis_l` and `basis_r` defining the left and right side bases.

For fast time evolution also at least the function
`gemv!(alpha, op::Operator, x::Ket, beta, result::Ket)` should be
implemented. Many other generic multiplication functions can be defined in
terms of this function and are provided automatically.
"""
abstract type Operator end

# Common error messages
arithmetic_unary_error(funcname, x::Operator) = throw(ArgumentError("$funcname is not defined for this type of operator: $(typeof(x)).\nTry to convert to another operator type first with e.g. full() or sparse()."))
arithmetic_binary_error(funcname, a::Operator, b::Operator) = throw(ArgumentError("$funcname is not defined for this combination of types of operators: $(typeof(a)), $(typeof(b)).\nTry to convert to a common operator type first with e.g. full() or sparse()."))
addnumbererror() = throw(ArgumentError("Can't add or subtract a number and an operator. You probably want 'op + identityoperator(op)*x'."))

length(a::Operator) = length(a.basis_l)::Int*length(a.basis_r)::Int
basis(a::Operator) = (check_samebases(a); a.basis_l)

# Arithmetic operations
+(a::Operator, b::Operator) = arithmetic_binary_error("Addition", a, b)
+(a::Number, b::Operator) = addnumbererror()
+(a::Operator, b::Number) = addnumbererror()

-(a::Operator) = arithmetic_unary_error("Negation", a)
-(a::Operator, b::Operator) = arithmetic_binary_error("Subtraction", a, b)
-(a::Number, b::Operator) = addnumbererror()
-(a::Operator, b::Number) = addnumbererror()

*(a::Operator, b::Operator) = arithmetic_binary_error("Multiplication", a, b)

dagger(a::Operator) = arithmetic_unary_error("Hermitian conjugate", a)

"""
    ishermitian(op::Operator)

Check if an operator is Hermitian.
"""
ishermitian(op::Operator) = arithmetic_unary_error(ishermitian, op)

"""
    tensor(x::Operator, y::Operator, z::Operator...)

Tensor product ``\\hat{x}⊗\\hat{y}⊗\\hat{z}⊗…`` of the given operators.
"""
tensor(a::Operator, b::Operator) = arithmetic_binary_error("Tensor product", a, b)
tensor(op::Operator) = op
tensor(operators::Operator...) = reduce(tensor, operators)

"""
    embed(basis1[, basis2], indices::Vector, operators::Vector)

Tensor product of operators where missing indices are filled up with identity operators.
"""
function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
               indices::Vector{Int}, operators::Vector{T}) where T<:Operator
    N = length(basis_l.bases)
    @assert length(basis_r.bases) == N
    @assert length(indices) == length(operators)
    sortedindices.check_indices(N, indices)
    tensor([i ∈ indices ? operators[findfirst(indices, i)] : identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i=1:N]...)
end
embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, op::Operator) = embed(basis_l, basis_r, Int[index], [op])
embed(basis::CompositeBasis, index::Int, op::Operator) = embed(basis, basis, Int[index], [op])
embed(basis::CompositeBasis, indices::Vector{Int}, operators::Vector{T}) where {T<:Operator} = embed(basis, basis, indices, operators)

"""
    embed(basis1[, basis2], operators::Dict)

`operators` is a dictionary `Dict{Vector{Int}, Operator}`. The integer vector
specifies in which subsystems the corresponding operator is defined.
"""
function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
               operators::Dict{Vector{Int}, T}) where T<:Operator
    @assert length(basis_l.bases) == length(basis_r.bases)
    N = length(basis_l.bases)
    if length(operators) == 0
        return identityoperator(T, basis_l, basis_r)
    end
    indices, operator_list = zip(operators...)
    operator_list = [operator_list...;]
    indices_flat = [indices...;]
    start_indices_flat = [i[1] for i in indices]
    complement_indices_flat = Int[i for i=1:N if i ∉ indices_flat]
    operators_flat = T[]
    if all([minimum(I):maximum(I);]==I for I in indices)
        for i in 1:N
            if i in complement_indices_flat
                push!(operators_flat, identityoperator(T, basis_l.bases[i], basis_r.bases[i]))
            elseif i in start_indices_flat
                push!(operators_flat, operator_list[findfirst(start_indices_flat, i)])
            end
        end
        return tensor(operators_flat...)
    else
        complement_operators = [identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in complement_indices_flat]
        op = tensor([operator_list; complement_operators]...)
        perm = sortperm([indices_flat; complement_indices_flat])
        return permutesystems(op, perm)
    end
end
embed(basis_l::CompositeBasis, basis_r::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:Operator} = embed(basis_l, basis_r, Dict([i]=>op_i for (i, op_i) in operators); kwargs...)
embed(basis::CompositeBasis, operators::Dict{Int, T}; kwargs...) where {T<:Operator} = embed(basis, basis, operators; kwargs...)
embed(basis::CompositeBasis, operators::Dict{Vector{Int}, T}; kwargs...) where {T<:Operator} = embed(basis, basis, operators; kwargs...)

"""
    trace(x::Operator)

Trace of the given operator.
"""
trace(x::Operator) = arithmetic_unary_error("Trace", x)

ptrace(a::Operator, index::Vector{Int}) = arithmetic_unary_error("Partial trace", a)

"""
    normalize(op)

Return the normalized operator so that its `trace(op)` is one.
"""
normalize(op::Operator) = op/trace(op)

"""
    normalize!(op)

In-place normalization of the given operator so that its `trace(x)` is one.
"""
normalize!(op::Operator) = throw(ArgumentError("normalize! is not defined for this type of operator: $(typeof(op)).\n You may have to fall back to the non-inplace version 'normalize()'."))

"""
    expect(op, state)

Expectation value of the given operator `op` for the specified `state`.

`state` can either be a (density) operator or a ket.
"""
expect(op::Operator, state::Ket) = dagger(state)*(op*state)
expect(op::Operator, state::Operator) = trace(op*state)

"""
    expect(index, op, state)

If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number.
"""
function expect(indices::Vector{Int}, op::Operator, state::Operator)
    N = length(state.basis_l.shape)
    indices_ = sortedindices.complement(N, indices)
    expect(op, ptrace(state, indices_))
end
function expect(indices::Vector{Int}, op::Operator, state::Ket)
    N = length(state.basis.shape)
    indices_ = sortedindices.complement(N, indices)
    expect(op, ptrace(state, indices_))
end
expect(index::Int, op::Operator, state) = expect([index], op, state)
expect(op::Operator, states::Vector) = [expect(op, state) for state=states]
expect(indices::Vector{Int}, op::Operator, states::Vector) = [expect(indices, op, state) for state=states]

"""
    variance(op, state)

Variance of the given operator `op` for the specified `state`.

`state` can either be a (density) operator or a ket.
"""
function variance(op::Operator, state::Ket)
    x = op*state
    stateT = dagger(state)
    stateT*op*x - (stateT*x)^2
end
function variance(op::Operator, state::Operator)
    expect(op*op, state) - expect(op, state)^2
end

"""
    variance(index, op, state)

If an `index` is given, it assumes that `op` is defined in the subsystem specified by this number
"""
function variance(indices::Vector{Int}, op::Operator, state::Operator)
    N = length(state.basis_l.shape)
    indices_ = sortedindices.complement(N, indices)
    variance(op, ptrace(state, indices_))
end
function variance(indices::Vector{Int}, op::Operator, state::Ket)
    N = length(state.basis.shape)
    indices_ = sortedindices.complement(N, indices)
    variance(op, ptrace(state, indices_))
end
variance(index::Int, op::Operator, state) = variance([index], op, state)
variance(op::Operator, states::Vector) = [variance(op, state) for state=states]
variance(indices::Vector{Int}, op::Operator, states::Vector) = [variance(indices, op, state) for state=states]

"""
    expm(op::Operator)

Operator exponential.
"""
expm(op::Operator) = throw(ArgumentError("expm() is not defined for this type of operator: $(typeof(op)).\nTry to convert to dense operator first with full()."))

permutesystems(a::Operator, perm::Vector{Int}) = arithmetic_unary_error("Permutations of subsystems", a)

"""
    identityoperator(a::Basis[, b::Basis])

Return an identityoperator in the given bases.
"""
identityoperator(::Type{T}, b1::Basis, b2::Basis) where {T<:Operator} = throw(ArgumentError("Identity operator not defined for operator type $T."))
identityoperator(::Type{T}, b::Basis) where {T<:Operator} = identityoperator(T, b, b)
identityoperator(op::T) where {T<:Operator} = identityoperator(T, op.basis_l, op.basis_r)

one(b::Basis) = identityoperator(b)
one(op::Operator) = identityoperator(op)

# Fast in-place multiplication
"""
    gemv!(alpha, a, b, beta, result)

Fast in-place multiplication of operators with state vectors. It
implements the relation `result = beta*result + alpha*a*b`.
Here, `alpha` and `beta` are complex numbers, while `result` and either `a`
or `b` are state vectors while the other one can be of any operator type.
"""
gemv!() = error("Not Implemented.")

"""
    gemm!(alpha, a, b, beta, result)

Fast in-place multiplication of of operators with DenseOperators. It
implements the relation `result = beta*result + alpha*a*b`.
Here, `alpha` and `beta` are complex numbers, while `result` and either `a`
or `b` are dense operators while the other one can be of any operator type.
"""
gemm!() = error("Not Implemented.")

# Helper functions to check validity of arguments
function check_ptrace_arguments(a::Operator, indices::Vector{Int})
    if !isa(a.basis_l, CompositeBasis) || !isa(a.basis_r, CompositeBasis)
        throw(ArgumentError("Partial trace can only be applied onto operators with composite bases."))
    end
    rank = length(a.basis_l.shape)
    if rank != length(a.basis_r.shape)
        throw(ArgumentError("Partial trace can only be applied onto operators wich have the same number of subsystems in the left basis and right basis."))
    end
    if rank == length(indices)
        throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use trace() instead."))
    end
    sortedindices.check_indices(length(a.basis_l.shape), indices)
    for i=indices
        if a.basis_l.shape[i] != a.basis_r.shape[i]
            throw(ArgumentError("Partial trace can only be applied onto subsystems that have the same left and right dimension."))
        end
    end
end
function check_ptrace_arguments(a::StateVector, indices::Vector{Int})
    if length(basis(a).shape) == length(indices)
        throw(ArgumentError("Partial trace can't be used to trace out all subsystems - use trace() instead."))
    end
    sortedindices.check_indices(length(basis(a).shape), indices)
end

samebases(a::Operator) = samebases(a.basis_l, a.basis_r)::Bool
samebases(a::Operator, b::Operator) = samebases(a.basis_l, b.basis_l)::Bool && samebases(a.basis_r, b.basis_r)::Bool
check_samebases(a::Operator) = check_samebases(a.basis_l, a.basis_r)

multiplicable(a::Operator, b::Ket) = multiplicable(a.basis_r, b.basis)
multiplicable(a::Bra, b::Operator) = multiplicable(a.basis, b.basis_l)
multiplicable(a::Operator, b::Operator) = multiplicable(a.basis_r, b.basis_l)

end # module


text:
using HDF5
using Test

macro gcvalid(args...)
    Expr(:block, quote
        GC.enable(true)
        GC.gc()
        GC.enable(false)
    end,
    [:(@test HDF5.isvalid($(esc(x)))) for x in args]...)
end

macro closederror(x)
    quote
        try
            $(esc(x))
        catch e
            isa(e, ErrorException) || rethrow(e)
            e.msg == "File or object has been closed" || error("Attempt to access closed object did not throw")
        end
    end
end

@testset "gc" begin

GC.enable(false)
fn = tempname()
for i = 1:10
    file = h5open(fn, "w")
    memtype_id = HDF5.API.h5t_create(HDF5.API.H5T_COMPOUND, 2*sizeof(Float64))
    HDF5.API.h5t_insert(memtype_id, "real", 0, HDF5.hdf5_type_id(Float64))
    HDF5.API.h5t_insert(memtype_id, "imag", sizeof(Float64), HDF5.hdf5_type_id(Float64))
    dt = HDF5.Datatype(memtype_id)
    commit_datatype(file, "dt", dt)
    ds = dataspace((2,))
    d = create_dataset(file, "d", dt, ds)
    g = create_group(file, "g")
    a = create_attribute(file, "a", dt, ds)
    @gcvalid dt ds d g a
    close(file)

    @closederror read(d)
    for obj in (d, g)
        @closederror read_attribute(obj, "a")
        @closederror write_attribute(obj, "a", 1)
    end
    for obj in (g, file)
        @closederror open_dataset(obj, "d")
        @closederror read_dataset(obj, "d")
        @closederror write_dataset(obj, "d", 1)
        @closederror read(obj, "x")
        @closederror write(obj, "x", "y")
    end
end
for i = 1:10
    file = h5open(fn, "r")
    dt = file["dt"]
    d = file["d"]
    ds = dataspace(d)
    g = file["g"]
    a = attributes(file)["a"]
    @gcvalid dt ds d g a
    close(file)
end
GC.enable(true)

let plist = create_property(HDF5.API.H5P_FILE_ACCESS)  # related to issue #620
    HDF5.API.h5p_close(plist)
    @test_nowarn finalize(plist)
end

rm(fn)

end # testset gc


text:
include("ttv_nplanet.jl")
function ttv_wrapper(tt0,nplanet::Int64,ntrans::Vector{Int64},params::Vector{T},jmax::Integer,EM::Bool) where T<:Real
  # These lines need modification for different choices of parameters:
  if nplanet == 2
    n1,n2 = ntrans
  elseif nplanet == 3
    n1,n2,n3 = ntrans
  elseif nplanet == 4
    n1,n2,n3,n4 = ntrans
  end

  # Call ttv_nplanet:
  ttv = ttv_nplanet(nplanet,jmax,ntrans,params[1:5*nplanet])
  # We measure transit times,not TTVs,so add back in the linear ephemeris:
  t01 = params[3]
  per1 = params[2]
  ttv1 = collect(range(t01,stop = t01+per1*(n1-1),length = n1)) #this doesnt account for skipped transits
  for i=1:n1
    ttv1[i]+= ttv[1,i]
  end
  t02 = params[8]
  per2 = params[7]
  ttv2 = collect(range(t02,stop = t02+per2*(n2-1),length = n2)) #check memory allocation >>>>>>>>>>>>
  for i=1:n2
    if EM
      ttv2[i] += ttv[2,i]
    else
      ts = params[end-2] #tmax sinphi0
      tc = params[end-1] #tmax cosphi0
      deltaphi = params[end]
      ttv2[i] += ttv[2,i] + ts*cos((i-1)*deltaphi) + tc*sin((i-1)*deltaphi)
    end
  end
  # If transit times of additional planets were observable these would need to be added in.
  #println("param2: ",param)
  return [ttv1;ttv2]  #check memory allocation >>>>>>>>>>>>
end

function chisquare(tt0,nplanet,ntrans,params,tt,sigtt,jmax,EM)#,fixp3::Bool = false,p3_cur::Float64 = 0.0)
  chisq = 0.0  #check memory allocation >>>>>>>>>>>>
  # println(params,tt[1],sigtt[1])
  tt_model = ttv_wrapper(tt0,nplanet,ntrans,params,jmax,EM) #,fixp3,p3_cur)
  for j=1:length(tt)
    chisq += (tt[j]-tt_model[j])^2/sigtt[j]^2
  end
  # println(nplanet)
  return chisq
end

max_stars_repo_name: bmlindor/TTVFaster

max_issues_repo_name: bmlindor/TTVFaster

max_forks_repo_name: bmlindor/TTVFaster


text:
export defΩ, SimpleΩ, LazyΩ

using ..IDS, ..Util, ..Tagging
using Distributions: Distribution
using Random: AbstractRNG
import ..Basis: replacetags

"Simplest, immutable Omega"
struct SimpleΩ{TAGS <: Tags, T} <: AbstractΩ
  data::T
  tags::TAGS
end

idtype(::SimpleΩ{TAGS, Dict{T, V}}) where {TAGS, T, V} = T

SimpleΩ(data) = SimpleΩ(data, Tags())
replacetags(ω::SimpleΩ, tags) = SimpleΩ(ω.data, tags)
(T::Type{<:Distribution})(π::SimpleΩ, args...) = ω.data[scope(ω)]

recurse(d::D, ω::SimpleΩ) where {D<:Distribution} = getindex(ω.data, scope(ω))::eltype(D)

traithastag(t::Type{SimpleΩ{TAGS, T}}, tag) where {TAGS, T} = traithastag(TAGS, tag)
traits(::Type{SimpleΩ{TAGS, T}}) where {TAGS, T} = traits(TAGS)

function resolve(dist, id, ω::SimpleΩ)
  # @assert false
  @show ω
  @show id
  @show ω.data[id]
  d, val = ω.data[id]
  val::eltype(dist)
end


text:
module AxelFraud

    using Agents
    using StatsBase
    using DataFrames
    using Feather
    using Query
    using Pipe
    using Random
    using ProgressMeter

    include("functions.jl")

    export initialize_model
    export run_config
    export run_random
    export run_batch
    export create_config_table

end


text:
#==============================================================================#
# EC2/test/runtests.jl
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

using AWSEC2
using Base.Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config()

#-------------------------------------------------------------------------------
# EC2 tests
#-------------------------------------------------------------------------------

@test_throws AWSCore.AWSException ec2_id(aws, "Not a real server name!!")

r = ec2(aws, Dict("Action"           => "DescribeImages",
                  "Filter.1.Name"    => "image-id",
                  "Filter.1.Value.1" => "ami-1ecae776"))

@test r["imagesSet"]["item"]["description"] ==
      "Amazon Linux AMI 2015.03.0 x86_64 HVM GP2"

#==============================================================================#
# End of file.
#==============================================================================#


text:
# Julia wrapper for header: /usr/include/scip/scip_cut.h
# Automatically generated using Clang.jl wrap_c

function SCIPgetCutEfficacy(scip, sol, cut)
    ccall((:SCIPgetCutEfficacy, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_SOL}, Ptr{SCIP_ROW}), scip, sol, cut)
end

function SCIPisCutEfficacious(scip, sol, cut)
    ccall((:SCIPisCutEfficacious, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_SOL}, Ptr{SCIP_ROW}), scip, sol, cut)
end

function SCIPisEfficacious(scip, efficacy)
    ccall((:SCIPisEfficacious, libscip), UInt32, (Ptr{SCIP_}, Cdouble), scip, efficacy)
end

function SCIPgetVectorEfficacyNorm(scip, vals, nvals)
    ccall((:SCIPgetVectorEfficacyNorm, libscip), Cdouble, (Ptr{SCIP_}, Ptr{Cdouble}, Cint), scip, vals, nvals)
end

function SCIPisCutApplicable(scip, cut)
    ccall((:SCIPisCutApplicable, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, cut)
end

function SCIPaddCut(scip, sol, cut, forcecut, infeasible)
    ccall((:SCIPaddCut, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SOL}, Ptr{SCIP_ROW}, UInt32, Ptr{UInt32}), scip, sol, cut, forcecut, infeasible)
end

function SCIPaddRow(scip, row, forcecut, infeasible)
    ccall((:SCIPaddRow, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}, UInt32, Ptr{UInt32}), scip, row, forcecut, infeasible)
end

function SCIPisCutNew(scip, row)
    ccall((:SCIPisCutNew, libscip), UInt32, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPaddPoolCut(scip, row)
    ccall((:SCIPaddPoolCut, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPdelPoolCut(scip, row)
    ccall((:SCIPdelPoolCut, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetPoolCuts(scip)
    ccall((:SCIPgetPoolCuts, libscip), Ptr{Ptr{SCIP_CUT}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNPoolCuts(scip)
    ccall((:SCIPgetNPoolCuts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetGlobalCutpool(scip)
    ccall((:SCIPgetGlobalCutpool, libscip), Ptr{SCIP_CUTPOOL}, (Ptr{SCIP_},), scip)
end

function SCIPcreateCutpool(scip, cutpool, agelimit)
    ccall((:SCIPcreateCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CUTPOOL}}, Cint), scip, cutpool, agelimit)
end

function SCIPfreeCutpool(scip, cutpool)
    ccall((:SCIPfreeCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CUTPOOL}}), scip, cutpool)
end

function SCIPaddRowCutpool(scip, cutpool, row)
    ccall((:SCIPaddRowCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CUTPOOL}, Ptr{SCIP_ROW}), scip, cutpool, row)
end

function SCIPaddNewRowCutpool(scip, cutpool, row)
    ccall((:SCIPaddNewRowCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CUTPOOL}, Ptr{SCIP_ROW}), scip, cutpool, row)
end

function SCIPdelRowCutpool(scip, cutpool, row)
    ccall((:SCIPdelRowCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CUTPOOL}, Ptr{SCIP_ROW}), scip, cutpool, row)
end

function SCIPseparateCutpool(scip, cutpool, result)
    ccall((:SCIPseparateCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CUTPOOL}, Ptr{SCIP_RESULT}), scip, cutpool, result)
end

function SCIPseparateSolCutpool(scip, cutpool, sol, result)
    ccall((:SCIPseparateSolCutpool, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CUTPOOL}, Ptr{SCIP_SOL}, Ptr{SCIP_RESULT}), scip, cutpool, sol, result)
end

function SCIPaddDelayedPoolCut(scip, row)
    ccall((:SCIPaddDelayedPoolCut, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPdelDelayedPoolCut(scip, row)
    ccall((:SCIPdelDelayedPoolCut, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_ROW}), scip, row)
end

function SCIPgetDelayedPoolCuts(scip)
    ccall((:SCIPgetDelayedPoolCuts, libscip), Ptr{Ptr{SCIP_CUT}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNDelayedPoolCuts(scip)
    ccall((:SCIPgetNDelayedPoolCuts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPgetDelayedGlobalCutpool(scip)
    ccall((:SCIPgetDelayedGlobalCutpool, libscip), Ptr{SCIP_CUTPOOL}, (Ptr{SCIP_},), scip)
end

function SCIPseparateSol(scip, sol, pretendroot, allowlocal, onlydelayed, delayed, cutoff)
    ccall((:SCIPseparateSol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_SOL}, UInt32, UInt32, UInt32, Ptr{UInt32}, Ptr{UInt32}), scip, sol, pretendroot, allowlocal, onlydelayed, delayed, cutoff)
end

function SCIPgetCuts(scip)
    ccall((:SCIPgetCuts, libscip), Ptr{Ptr{SCIP_ROW}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNCuts(scip)
    ccall((:SCIPgetNCuts, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPclearCuts(scip)
    ccall((:SCIPclearCuts, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPremoveInefficaciousCuts(scip)
    ccall((:SCIPremoveInefficaciousCuts, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end