text:
using LoopVectorization, OffsetArrays, Test

@testset "copy" begin

    function copyavx1!(x, y)
        @turbo for i ∈ eachindex(x)
            x[i] = y[i]
        end
    end
    function copy_avx1!(x, y)
        @_avx for i ∈ eachindex(x)
            x[i] = y[i]
        end
    end
    function copyavx2!(x, y)
        @turbo for i ∈ eachindex(x)
            yᵢ = y[i]
            x[i] = yᵢ
        end
    end
    function copy_avx2!(x, y)
        @_avx for i ∈ eachindex(x)
            yᵢ = y[i]
            x[i] = yᵢ
        end
    end
    function offset_copy!(A, B)
        @inbounds for i=1:size(A,1), j=1:size(B,2)
            A[i,j+2] = B[i,j]
        end
    end
    function offset_copyavx1!(A, B)
        @turbo for i=1:size(A,1), j=1:size(B,2)
           A[i,j+2] = B[i,j]
        end
    end
    function offset_copy_avx1!(A, B)
        @_avx for i=1:size(A,1), j=1:size(B,2)
              @inbounds A[i,j+2] = B[i,j]
        end
    end
    function offset_copyavx2!(A, B)
        @turbo for i=1:size(A,1), j=1:size(B,2)
            Bᵢⱼ = B[i,j]
            A[i,j+2] = Bᵢⱼ
        end
    end
    function offset_copy_avx2!(A, B)
        @_avx for i=1:size(A,1), j=1:size(B,2)
            Bᵢⱼ = B[i,j]
            A[i,j+2] = Bᵢⱼ
        end
    end
    function make2point3avx!(x)
        a = 1.742416161578685
        b = 1.5
        @turbo for i ∈ eachindex(x)
            x[i] = a ^ b
        end
    end
    function make2point3_avx!(x)
        @_avx for i ∈ eachindex(x)
            x[i] = 2.3
        end
    end
    function make23avx!(x)
        @turbo for i ∈ eachindex(x)
            @inbounds x[i] = 23
        end
    end
    function make23_avx!(x)
        @_avx for i ∈ eachindex(x)
            x[i] = 23
        end
    end
    function myfillavx!(x, a)
        @turbo for i ∈ eachindex(x)
            x[i] = a
        end
    end
    function myfill_avx!(x, a)
        @_avx for i ∈ eachindex(x)
            x[i] = a
        end
    end
    function issue_256!(A)
        @turbo for i = 1:size(A,1)
            A[i,(x=1,).x] = 0
        end
        A
    end

    function reversecopy1!(B, A)
        for i in eachindex(A)
            B[i] = A[11-i]
        end
        B
    end
    function reversecopy1avx!(B, A)
        @turbo for i in eachindex(A)
            B[i] = A[11-i]
        end
        B
    end
    function reversecopy2!(B, A)
        for i in eachindex(B)
            B[i] = A[-i]
        end
        B
    end
    function reversecopy2avx!(B, A)
        @turbo for i in eachindex(B)
            B[i] = A[-i]
        end
        B
    end
    function reversecopy3avx!(B, A)
        n = length(A)
        @turbo for i in 1:n
            B[i] = A[n+1-i]
        end
        B
    end
    function copy3!(B, A)
        @assert (length(B) ≥ 3) && (length(A) ≥ 3)
        @turbo for i in 1:3
            B[i] = A[i]
        end
        B
    end
    function copyselfdot!(s, x)
        m = zero(eltype(x))
        @turbo for i ∈ 1:2
            sᵢ = x[i]
            s[i] = sᵢ
            m += sᵢ * sᵢ
        end
        m
    end
  function scattercopyavx!(H,a,j)
    @turbo for i ∈ eachindex(j), k ∈ eachindex(a)
      H[j[i],k] = a[k]
    end
    H
  end
  function scattercopy!(H,a,j)
    @inbounds for i ∈ eachindex(j), k ∈ eachindex(a)
      H[j[i],k] = a[k]
    end
    H
  end

  function issue279!(output, input, idx)
    @turbo for j in axes(output, 2), i in axes(output, 1)
      output[i, j, idx] = input[1, 1, i, j, idx]
    end
    output
  end

  function issue279_plain!(output, input, idx)
    for j in axes(output, 2), i in axes(output, 1)
      output[i, j, idx] = input[1, 1, i, j, idx]
    end
    output
  end
  function issue280!(dest, src)
    @turbo for i in indices((dest, src), (2, 2))
      dest[1, i] = src[2, i]
      dest[2, i] = src[1, i]
    end
  end

  function fill_with_3rd!(dest, src)
    p = 1
    @turbo for i ∈ eachindex(dest)
      dest[i] = src[3*p]
    end
  end
  function fill_with_r1c3!(dest, src::AbstractMatrix)
    p = 1
    @turbo for i ∈ eachindex(dest)
      dest[i] = src[p, 3*p]
    end
  end

    for T ∈ (Float32, Float64, Int32, Int64)
        @show T, @__LINE__
        R = T <: Integer ? (-T(100):T(100)) : T
        x = rand(R, 237);
        q1 = similar(x); q2 = similar(x);

        fill!(q2, -999999); copyavx1!(q2, x);
        @test x == q2
        fill!(q2, -999999); copy_avx1!(q2, x);
        @test x == q2
        fill!(q2, -999999); copyavx2!(q2, x);
        @test x == q2
        fill!(q2, -999999); copy_avx2!(q2, x);
        @test x == q2
        fill!(q2, -999999); @turbo q2 .= x;
        @test x == q2
        @test all(iszero, issue_256!(reshape(x,(length(x),1))))

        B = rand(R, 79, 83);
        A1 = zeros(T, 79, 85);
        A2 = zeros(T, 79, 85);
        offset_copy!(A1, B);
        fill!(A2, 0); offset_copyavx1!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copyavx2!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copy_avx1!(A2, B);
        @test A1 == A2
        fill!(A2, 0); offset_copy_avx2!(A2, B);
        @test A1 == A2
        fill_with_3rd!(A1, B)
        @test all(==(B[3]), A1)
        fill_with_3rd!(A1, eachindex(B))
        @test all(==(3), A1)

        fill_with_r1c3!(A1, B)
        @test all(==(B[1,3]), A1)
        fill_with_r1c3!(A1, LinearIndices(B))# check args fails
        @test all(==(LinearIndices(B)[1,3]), A1)

        a = rand(R)
        myfillavx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfill_avx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfill_avx!(x, a);
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        myfillavx!(x, a);
        fill!(q2, a);
        @test x == q2
        q2 .= 23;
        fill!(q1, -99999); make23_avx!(q1);
        @test q2 == q1
        fill!(q1, -99999); make23avx!(q1);
        @test q2 == q1
        if T <: Union{Float32,Float64}
            make2point3avx!(x)
            fill!(q2, 2.3)
            @test x == q2
            fill!(x, -999999); make2point3_avx!(x)
            @test x == q2
        end
        a = rand(R)
        @turbo x .= a;
        fill!(q2, a);
        @test x == q2
        a = rand(R)
        @turbo x .= a;
        fill!(q2, a);
        @test x == q2

        @test reversecopy1!(zeros(T, 10), collect(1:10)) == reversecopy1avx!(zeros(T, 10), collect(1:10))
        @test reversecopy2!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1)) == reversecopy2avx!(zeros(T, 10), OffsetArray(collect(1:10), -10:-1))
        @test 177:-1:1 == reversecopy3avx!(zeros(T, 177), collect(1:177))

        x = rand(R, 3); y = similar(x);
        @test copy3!(y, x) == x
        fill!(y,0);
        @test copyselfdot!(y, x) ≈ x[1]^2 + x[2]^2
      @test view(x, 1:2) == view(y, 1:2)

      H0 = zeros(T,10,10); H1 = zeros(T,10,10);
      j = [1,2,5,8]; a = rand(R,10);
      @test scattercopyavx!(H0, a, j) == scattercopy!(H1, a, j)

      input = rand(R, 2, 2, 5, 5, 1); output = Array{T}(undef, size(input)[3:end]...); output_plain = similar(output);

      @test issue279!(output, input, 1) ≈ issue279_plain!(output_plain, input, 1)

      src = rand(R, 2, 17); dest = similar(src);
      issue280!(dest, src)
      @test dest ≈ vcat(view(src,2,:)',view(src,1,:)')
      if VERSION ≥ v"1.6"
        src2 = reinterpret(reshape,T,Vector{Tuple{T,T}}(undef, 17)); src2 .= src;
        dest2 = reinterpret(reshape,T,Vector{Tuple{T,T}}(undef, 17));
        issue280!(dest2, src2)
        @test dest2 ≈ vcat(view(src,2,:)',view(src,1,:)')
      end

      rng = 10:20
      arr1 = zeros(30);
      arr2 = zeros(30);
      ifirst = 5
      @turbo for i in rng
        arr1[i-ifirst] += 1
      end
      @. arr2[rng-ifirst] += 1;
      @test arr1 == arr2
    end
end


text:
module IOlib
using CSV, DataFrames, LinearAlgebra, DelimitedFiles, YAML, Printf

export supplyusedata, inputoutputcalc, prices_init, parse_input_file, write_matrix_to_csv, write_vector_to_csv,
       IOdata, PriceData, ExogParams, ϵ

"Global variable to store the result of parsing the configuration file"
global_params = nothing

"Small number for avoiding divide-by-zero problems"
const ϵ = 1.0e-11

"Values calculated from supply-use and associated tables"
mutable struct IOdata
    D::Array{Float64,2} # np x ns
    S::Array{Float64,2} # ns x np
	Vnorm::Array{Float64,2} # ns x np
    F::Array{Float64,1} # np
    I::Array{Float64,1} # np
    X::Array{Float64,1} # np
	M::Array{Float64,1} # np
    W::Array{Float64,1} # ns
    g::Array{Float64,1} # ns
    marg_pos_ratio::Array{Float64,1} # np
    marg_neg_share::Array{Float64,1} # np
    m_frac::Array{Float64,1} # np
    μ::Array{Float64,1} # ns
    τd::Array{Float64,1} # np
end

"Price indices"
mutable struct PriceData
    pb::Array{Float64,1}
    pd::Array{Float64,1}
    pp::Array{Float64,1}
    pw::Array{Float64,1}
    Pg::Float64
    Pc::Float64
end

"User-specified parameters"
mutable struct ExogParams
    πw::Array{Float64,1} # World inflation rate x year
    world_grs::Array{Float64,1} # Real world GDP growth rate x year
	working_age_grs::Array{Float64,1} # Growth rate of the working age population x year
	αKV::Array{Float64,1} # Kaldor-Verdoorn coefficient x year
	βKV::Array{Float64,1} # Kaldor-Verdoorn intercept x year
	xr::Array{Float64,1} # Nominal exchange rate x year
    δ::Array{Float64,1} # ns
    export_elast_demand::Array{Any,1} # np x year
    wage_elast_demand::Array{Any,1} # np x year
end

"""
    write_matrix_to_csv(filename, array, colnames, rownames)

Write a matrix to a CSV file.
"""
function write_matrix_to_csv(filename, array, rownames, colnames)
	open(filename, "w") do io
        # Write header
		write(io, string("", ',', join(colnames, ','), "\r\n"))
        # Write each row
        for i in 1:length(rownames)
            write(io, string(rownames[i], ',', join(array[i,:], ','), "\r\n"))
        end
	end
end

"""
    write_vector_to_csv(filename, vector, varname, rownames)

Write a vector to a CSV file.
"""
function write_vector_to_csv(filename, vector, varname, rownames)
	open(filename, "w") do io
        # Write header
		write(io, string("", ',', varname, "\r\n"))
        for i in 1:length(rownames)
            write(io, string(rownames[i], ',', vector[i], "\r\n"))
        end
	end
end

"""
    string_to_float(str)

Convert strings in Dataframe to floats and fill in missing values with 0.
"""
function string_to_float(str)
    try return(parse(Float64, str)) catch
        return(0.0) # Missing values set to zero
    end
end

"""
    char_index_to_int(str)

Convert an index in the form of characters as used in Excel (e.g., "AC") into an integer index
"""
function char_index_to_int(str)
	str_up = uppercase(str)
	mult = 26^length(str_up) # 26 is number of characters in the Roman alphabet
	ndx = 0
	for c in str_up
		mult /= 26
		ndx += mult * (Int(c) - Int('A') + 1)
	end
	return(Int(ndx))
end

"""
    excel_ref_to_rowcol(str)

Convert a cell reference in Excel form (e.g., "AB3") to [row, col] (eg., [3,28])
"""
function excel_ref_to_rowcol(str)
	m = match(r"([A-Za-z]+)([0-9]+)", strip(str)) # Be tolerant of leading and trailing spaces
	if m === nothing
		error("Cell reference '$str' is not in the correct format (e.g., 'AB3')")
		return nothing
	else
		c = m.captures
		return[parse(Int64,c[2]),char_index_to_int(c[1])]
	end
end

"""
    excel_cell_to_float(df, str)

Find the value of a cell in dataframe df referenced in Excel form (e.g., "AB3") and convert to float
"""
function excel_cell_to_float(df, str)
	rowcol = excel_ref_to_rowcol(str)
	return string_to_float(df[rowcol[1],rowcol[2]])
end

"""
    excel_range_to_rowcol_pair(str)

Convert a cell range in Excel form (e.g., "AB3:BG10") to [[row, col],[row,col]]
"""
function excel_range_to_rowcol_pair(str)
	range_ends = split(str,":")
	if length(range_ends) != 2
		error("Cell range '$str' is not in the correct format (e.g., 'AB3:BG10')")
		return nothing
	else
		return excel_ref_to_rowcol.(range_ends)
	end
end

"""
    df_to_mat(df)

First convert strings in a Dataframe to floats and then convert the Dataframe to a Matrix.
"""
function df_to_mat(df)
    for x = 1:size(df,2)
        df[!,x] = map(string_to_float, df[!,x]) # converts string to float
    end
    mat = Matrix(df) # convert dataframe to matrix
    return mat
end

"""
    excel_range_to_mat(df, str)

Extract a range in Excel format from a dataframe and return a numeric matrix
"""
function excel_range_to_mat(df, str)
	rng = excel_range_to_rowcol_pair(str)
	return df_to_mat(df[rng[1][1]:rng[2][1],rng[1][2]:rng[2][2]])
end

"""
    parse_input_file(YAML_file::String; force::Bool = false, include_energy_sectors::Bool = false)

Read in YAML input file and convert ranges if necessary.
Force: reset global_params
"""
function parse_input_file(YAML_file::String; force::Bool = false, include_energy_sectors::Bool = false)
	global global_params

	if !force & !isnothing(global_params)
		return global_params
	end

    global_params = YAML.load_file(YAML_file)

    if include_energy_sectors
        output_folder_name = string(global_params["output_folder"], "_full")
    else
        output_folder_name = global_params["output_folder"]
    end
    global_params["results_path"] = joinpath("outputs/", output_folder_name, "results")
    global_params["calibration_path"] = joinpath("outputs/", output_folder_name, "calibration")
    global_params["diagnostics_path"] = joinpath("outputs/", output_folder_name, "diagnostics")

	# First, check if any sectors or products should be excluded because production is zero (or effectively zero, for products)
	SUT_df = CSV.read(joinpath("inputs",global_params["files"]["SUT"]), header=false, DataFrame)
	M = vec(sum(excel_range_to_mat(SUT_df, global_params["SUT_ranges"]["imports"]), dims=2))
	supply_table = excel_range_to_mat(SUT_df, global_params["SUT_ranges"]["supply_table"])
	qs = vec(sum(supply_table, dims=2))
	g = vec(sum(supply_table, dims=1))
    M_equiv = transpose(supply_table) * Diagonal(1 ./ (qs .+ ϵ)) * M
    # Remove all sectors for which there is negligible domestic production relative to imports or where g = 0
	θ = global_params["domestic_production_share_threshold"]/100
	zero_domprod_ndxs = findall(x -> x <= 0, min(g, (1 - θ) * g - θ * M_equiv))
    # Remove all products for which there is negligible total supply (domestic + imported) relative to total domestic output
	zero_prod_ndxs = findall(x -> abs(x) < ϵ, (qs + M)/sum(g))

	all_sectors = CSV.read(joinpath("inputs",global_params["files"]["sector_info"]), header=1, types=Dict(:code => String, :name => String), select=[:code,:name], NamedTuple)
	all_products = CSV.read(joinpath("inputs",global_params["files"]["product_info"]), header=1, types=Dict(:code => String, :name => String), select=[:code,:name], NamedTuple)
	sector_codes = all_sectors[:code]
    product_codes = all_products[:code]

	excluded_sectors = vcat(global_params["excluded_sectors"]["territorial_adjustment"],
							global_params["excluded_sectors"]["others"])
    excluded_products = vcat(global_params["excluded_products"]["territorial_adjustment"],
                            global_params["excluded_products"]["others"])
   if !include_energy_sectors
        excluded_sectors = vcat(excluded_sectors, global_params["excluded_sectors"]["energy"])
        excluded_products = vcat(excluded_products, global_params["excluded_products"]["energy"])
    end
	# These are the indexes used for the Macro model
	user_defined_sector_ndxs = findall(.!in(excluded_sectors).(sector_codes))
	user_defined_product_ndxs = findall(.!in(excluded_products).(product_codes))
	global_params["sector-indexes"] = sort(setdiff(user_defined_sector_ndxs, zero_domprod_ndxs))
	global_params["product-indexes"] = sort(setdiff(user_defined_product_ndxs, zero_prod_ndxs))

    global_params["energy-sector-indexes"] = []
    global_params["energy-product-indexes"] = []
    if !include_energy_sectors
        user_defined_energy_sector_ndxs = findall(in(global_params["excluded_sectors"]["energy"]).(sector_codes))
        user_defined_energy_product_ndxs = findall(in(global_params["excluded_products"]["energy"]).(sector_codes))
        global_params["energy-sector-indexes"] = sort(setdiff(user_defined_energy_sector_ndxs, zero_domprod_ndxs))
        global_params["energy-product-indexes"] = sort(setdiff(user_defined_energy_product_ndxs, zero_prod_ndxs))
    end

	global_params["terr-adj-sector-indexes"] = findall(in(global_params["excluded_sectors"]["territorial_adjustment"]).(sector_codes))
	global_params["terr-adj-product-indexes"] = findall(in(global_params["excluded_products"]["territorial_adjustment"]).(sector_codes))

	global_params["included_sector_names"] = all_sectors[:name][global_params["sector-indexes"]]
	global_params["included_sector_codes"] = all_sectors[:code][global_params["sector-indexes"]]
	global_params["included_product_names"] = all_products[:name][global_params["product-indexes"]]
	global_params["included_product_codes"] = all_products[:code][global_params["product-indexes"]]

	# These indices are with respect to the list of included product codes, not the full list
	global_params["non-tradeable-range"] = findall(in(global_params["non_tradeable_products"]).(global_params["included_product_codes"]))

	global_params["LEAP_sector_names"] = [x["name"] for x in global_params["LEAP-sectors"]]
	LEAP_sector_codes = [x["codes"] for x in global_params["LEAP-sectors"]]
	global_params["LEAP_sector_indices"] = Array{Any}(undef,length(LEAP_sector_codes))
	for i in 1:length(LEAP_sector_codes)
		global_params["LEAP_sector_indices"][i] = findall([in(x, LEAP_sector_codes[i]) for x in global_params["included_sector_codes"]])
	end

    return global_params
end

"""
    get_var_params(param_file::String)

Pull in user-specified parameters from different CSV input files with filenames specified in the YAML configuration file.
"""
function get_var_params(param_file::String)
    # Return an ExogParams struct
    retval = ExogParams(
                    Array{Float64}(undef, 0), # πw
                    Array{Float64}(undef, 0), # world_grs
					Array{Float64}(undef, 0), # working_age_grs
					Array{Float64}(undef, 0), # αKV
					Array{Float64}(undef, 0), # βKV
					Array{Float64}(undef, 0), # xr
                    Array{Float64}(undef, 0), # δ
                    [], # export_elast_demand
                    []) # wage_elast_demand

    params = parse_input_file(param_file)
    base_year = params["years"]["start"]
    final_year = params["years"]["end"]
    sec_ndxs = params["sector-indexes"]
    prod_ndxs = params["product-indexes"]

    # Read in info from CSV files
    sector_info = CSV.read(joinpath("inputs",params["files"]["sector_info"]), DataFrame)
    product_info = CSV.read(joinpath("inputs",params["files"]["product_info"]), DataFrame)
    time_series = CSV.read(joinpath("inputs",params["files"]["time_series"]), DataFrame)
    prod_codes = product_info[prod_ndxs,:code]

    #--------------------------------------------------------------------------------------
    # Sector-specific
    #--------------------------------------------------------------------------------------
    # Use same value for each year, but it can vary by sector
    retval.δ = sector_info[sec_ndxs,:depr_rate]

    #--------------------------------------------------------------------------------------
    # Product-specific
    #--------------------------------------------------------------------------------------
    #--- Demand model parameters
    # Note, these are used later to make time-varying parameters
    # Export demand
    export_elast_decay = params["export_elast_demand"]["decay"]
    # Get initial value
    export_elast_demand0 = product_info[prod_ndxs,:export_elast_demand0]
    # Constrain to grow with global economy in the long run (with elasticity 1.0)
    asympt_export_elast = min.(convert(Array,export_elast_demand0), 1.0 * ones(length(export_elast_demand0)))

    # Household demand
    wage_elast_decay = params["wage_elast_demand"]["decay"]
    # Get initial value
    wage_elast_demand0 = product_info[prod_ndxs,:wage_elast_demand0]
    # For asymptotic value, first bring below 1.0
    asympt_wage_elast = min.(convert(Array,wage_elast_demand0), 1.0 * ones(length(wage_elast_demand0)))
    # Next, for Engel's Law, set long-run elasticity for food below 1.0
    engel_products = params["wage_elast_demand"]["engel_prods"]
    engel_ndxs = findall(x -> x in engel_products, prod_codes)
    asympt_wage_elast[engel_ndxs] = params["wage_elast_demand"]["engel_asympt_elast"] * ones(length(engel_products))

    #--------------------------------------------------------------------------------------
    # Time series
    #--------------------------------------------------------------------------------------
    # World inflation rate (apply to world prices only, others induced)
    world_infl_temp = time_series[!,:world_infl_rate]
    world_infl_default = params["global-params"]["infl_default"]
    world_infl_start = floor(Int64, time_series[1,:year])
    world_infl_end = world_infl_start + length(world_infl_temp) - 1

    # World GDP growth rate (including historical, for calibration)
    world_grs_temp = time_series[!,:world_gr]
    world_grs_default = params["global-params"]["gr_default"]
    world_grs_start = floor(Int64, time_series[1,:year])
    world_grs_end = world_grs_start + length(world_grs_temp) - 1

	# Working age growth rate (No default -- must provide all values for specified time period)
    working_age_grs_temp = time_series[!,:working_age_gr]
    working_age_grs_start = floor(Int64, time_series[1,:year])
    working_age_grs_end = working_age_grs_start + length(working_age_grs_temp) - 1

	# Exchange rate (No default -- must provide all values for specified time period)
    xr_temp = time_series[!,:exchange_rate]
    xr_start = floor(Int64, time_series[1,:year])
    xr_end = xr_start + length(xr_temp) - 1

    # Kaldor-Verdoorn parameters
    αKV_temp = time_series[!,:KV_coeff]
    αKV_default = params["labor-prod-fcn"]["KV_coeff_default"]
    αKV_start = floor(Int64, time_series[1,:year])
    αKV_end = αKV_start + length(αKV_temp) - 1

	βKV_temp = time_series[!,:KV_intercept]
    βKV_default = params["labor-prod-fcn"]["KV_intercept_default"]
    βKV_start = floor(Int64, time_series[1,:year])
    βKV_end = βKV_start + length(βKV_temp) - 1
    #--------------------------------------------------------------------------------------
    # Fill in by looping over years
    #--------------------------------------------------------------------------------------
    for year in base_year:final_year
        # creates world inflation rate array for use in the model, using model years
        if year in world_infl_start:world_infl_end
            push!(retval.πw, world_infl_temp[year - world_infl_start + 1])
        else
            push!(retval.πw, world_infl_default)
        end

        # creates world growth rate array for use in the model, using model years
        if year in world_grs_start:world_grs_end
            push!(retval.world_grs, world_grs_temp[year - world_grs_start + 1])
        else
            push!(retval.world_grs, world_grs_default)
        end

		# creates working age growth rate array for use in the model, using model years
        if year in working_age_grs_start:working_age_grs_end
            push!(retval.working_age_grs, working_age_grs_temp[year - working_age_grs_start + 1])
        else
			error_string = "Year is not in range " * string(working_age_grs_start) * ":" * string(working_age_grs_end) * " for working age growth rates"
            throw(DomainError(year, error_string))
        end

		# creates exchange rate array for use in the model, using model years
        if year in xr_start:xr_end
            push!(retval.xr, xr_temp[year - xr_start + 1])
        else
			error_string = "Year is not in range " * string(xr_start) * ":" * string(xr_end) * " for exchange rates"
            throw(DomainError(year, error_string))
        end

		# fill in Kaldor-Verdoorn coefficient and intercept
        if year in αKV_start:αKV_end
            push!(retval.αKV, αKV_temp[year - αKV_start + 1])
        else
            push!(retval.αKV, αKV_default)
        end
		if year in βKV_start:βKV_end
            push!(retval.βKV, βKV_temp[year - βKV_start + 1])
        else
            push!(retval.βKV, βKV_default)
        end

        deltat = max(0, year - base_year + 1) # Year past the end of calibration period

        # Bring high values down, but don't bring low values up
        export_elast_demand_temp = export_elast_demand0 - (1 - exp(-export_elast_decay * deltat)) * (export_elast_demand0 - asympt_export_elast)
        push!(retval.export_elast_demand, export_elast_demand_temp)

        wage_elast_demand_temp = wage_elast_demand0 - (1 - exp(-wage_elast_decay * deltat)) * (wage_elast_demand0 - asympt_wage_elast)
        push!(retval.wage_elast_demand, wage_elast_demand_temp)
    end

    return retval

end

"""
    energy_nonenergy_link_measure(param_file::String)

Calculate a measure of how important the demand for non-energy goods from the energy sector is.

The measure is a ratio of the sum of the elements of two Leontief inverses: one for the whole matrix of
technical coefficients and the other for the matrix with A_{NE} sub-block excluded.

The sum of the elements of a Leontief inverse can be interpreted as the change in total output
from an increase in final demand that is the same across all sectors.
"""
function energy_nonenergy_link_measure(param_file::String)
	params = parse_input_file(param_file)

	full_sector_ndxs = sort(vcat(params["sector-indexes"],params["energy-sector-indexes"]))
	full_product_ndxs = sort(vcat(params["product-indexes"],params["energy-product-indexes"]))
	ns = length(full_sector_ndxs)

	SUT_df = CSV.read(joinpath("inputs",params["files"]["SUT"]), header=false, DataFrame)

	#--------------------------------
	# Compute S matrix
	#--------------------------------
	supply_table = excel_range_to_mat(SUT_df, params["SUT_ranges"]["supply_table"])[full_product_ndxs,full_sector_ndxs]
	V = transpose(supply_table)
	qs = vec(sum(V, dims=1))
	g = vec(sum(V, dims=2))
	S = V * Diagonal(1.0 ./ (qs .+ ϵ))

	#--------------------------------
	# Compute D matrix
	#--------------------------------
	use_table = excel_range_to_mat(SUT_df, params["SUT_ranges"]["use_table"])[full_product_ndxs,full_sector_ndxs]
	D = use_table * Diagonal(1.0 ./ (g .+ ϵ))

	#--------------------------------
	# Compute A matrix and Leontief matrix
	#--------------------------------
	A = S * D
	L = inv(LinearAlgebra.I - A)

	#--------------------------------
	# Set A_NE set to zero and compute the Leontief matrix
	#--------------------------------
	A_reduced = A
	for col in 1:ns
		if !in(full_sector_ndxs[col], params["energy-sector-indexes"])
			continue
		end
		for row in 1:ns
			if in(full_sector_ndxs[row], params["sector-indexes"])
				A_reduced[row,col] = 0.0
			end
		end
	end
	L_reduced = inv(LinearAlgebra.I - A_reduced)

	#--------------------------------
	# Calculate and return measure
	#--------------------------------
	return 1.0 - sum(L_reduced)/sum(L)

end

"""
    supplyusedata(param_file::String)

Pull in supply-use data from CSV input files.
"""
function supplyusedata(param_file::String)
    retval = IOdata(Array{Float64}(undef, 0, 0), # D
                    Array{Float64}(undef, 0, 0), # S
                    Array{Float64}(undef, 0, 0), # Vnorm
                    Array{Float64}(undef, 0), # F
                    Array{Float64}(undef, 0), # I
                    Array{Float64}(undef, 0), # X
                    Array{Float64}(undef, 0), # M
                    Array{Float64}(undef, 0), # W
                    Array{Float64}(undef, 0), # g
                    Array{Float64}(undef, 0), # marg_pos_ratio
                    Array{Float64}(undef, 0), # marg_neg_share
                    Array{Float64}(undef, 0), # m_frac
                    Array{Float64}(undef, 0), # μ
                    Array{Float64}(undef, 0)) # τd
    # Reads in key parameters from the YAML config file
    params = parse_input_file(param_file)

	sector_ndxs = params["sector-indexes"]
	product_ndxs = params["product-indexes"]

	ns = length(sector_ndxs)
	np = length(product_ndxs)

	terr_adj_sector_ndx = params["terr-adj-sector-indexes"]
	terr_adj_product_ndx = params["terr-adj-product-indexes"]

    SUT_df = CSV.read(joinpath("inputs",params["files"]["SUT"]), header=false, DataFrame)

	#--------------------------------
	# Supply table
	#--------------------------------
	supply_table = excel_range_to_mat(SUT_df, params["SUT_ranges"]["supply_table"])[product_ndxs,sector_ndxs]
	V = transpose(supply_table)
	qs = vec(sum(V, dims=1))
	retval.S = V * Diagonal(1 ./ (qs .+ ϵ))
    retval.g = vec(sum(V, dims=2))
	retval.Vnorm = Diagonal(1 ./ retval.g) * V

	# Get total supply column
	tot_supply = excel_range_to_mat(SUT_df, params["SUT_ranges"]["tot_supply"])[product_ndxs]

	#--------------------------------
	# Use table
	#--------------------------------
	use_table = excel_range_to_mat(SUT_df, params["SUT_ranges"]["use_table"])[product_ndxs,sector_ndxs]
	retval.D = use_table * Diagonal(1.0 ./ retval.g)

	#--------------------------------
	# Margins
	#--------------------------------
	margins = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["margins"])[product_ndxs,:], dims=2))
	margins_stat_adj = sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["margins"])[terr_adj_product_ndx,:])
	# Calculate margin ratios
    margins_pos = max.(margins, zeros(np))
    margins_pos = margins_pos * (1.0 + margins_stat_adj/sum(margins_pos))
    margins_neg = max.(-margins, zeros(np))
    margins_neg = margins_neg * (1.0 - margins_stat_adj/sum(margins_neg))
	retval.marg_pos_ratio = margins_pos ./ tot_supply
	retval.marg_neg_share = margins_neg / sum(margins_neg)

	#--------------------------------
	# Taxes
	#--------------------------------
	taxes = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["taxes"])[product_ndxs,:], dims=2))
	retval.τd = taxes ./ tot_supply

	#--------------------------------
	# Imports
	#--------------------------------
	M = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["imports"])[product_ndxs,:], dims=2))
	M_stat_adj = sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["imports"])[terr_adj_product_ndx,:])
    M = M * (1.0 + M_stat_adj/(sum(M) + ϵ))
    M[params["non-tradeable-range"]] .= 0.0 # If it is declared non-tradeable, set imports to zero
	retval.M = M

	#--------------------------------
	# Exports -- will be adjusted later
	#--------------------------------
	X = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["exports"])[product_ndxs,:], dims=2))
	X_stat_adj = sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["exports"])[terr_adj_product_ndx,:])
    X = X * (1.0 + X_stat_adj/(sum(X) + ϵ))
    X[params["non-tradeable-range"]] .= 0.0 # If it is declared non-tradeable, set exports to zero
	retval.X = X

	#--------------------------------
	# Final demand -- will be adjusted later
	#--------------------------------
	F = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["final_demand"])[product_ndxs,:], dims=2))
	F_stat_adj = sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["final_demand"])[terr_adj_product_ndx,:])
    F = F * (1.0 + F_stat_adj/(sum(F) + ϵ))
	retval.F = F

	#--------------------------------
	# Investment -- will be adjusted later
	#--------------------------------
	I = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["investment"])[product_ndxs,:], dims=2))
	I_stat_adj = sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["investment"])[terr_adj_product_ndx,:])
    I = I * (1.0 + I_stat_adj/(sum(I) + ϵ))
	retval.I = I

	#--------------------------------
	# Allocate imports (by a common fraction for each product) to final demand and intermediate demand
	#--------------------------------
	tot_int_sup = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["tot_intermediate_supply"])[product_ndxs,:], dims=2))
    # The fraction m_frac applies to imports excluding imported investment goods
    retval.m_frac = retval.M ./ (tot_int_sup + retval.F + retval.I)
	#--------------------------------
	# Correct demands
	#--------------------------------
	stock_change = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["stock_change"])[product_ndxs,:], dims=2))
    # Correct for stock changes and tax leakage
    corr = 1 .+ (stock_change - taxes) ./ (retval.F + retval.I + retval.X .+ ϵ)
    retval.F = corr .* retval.F
    retval.X = corr .* retval.X
    retval.I = corr .* retval.I

	#--------------------------------
	# Wages
	#--------------------------------
	retval.W = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["wages"])[:,sector_ndxs], dims=1))

	#--------------------------------
	# Profit margin
	#--------------------------------
	tot_int_dmd = vec(sum(excel_range_to_mat(SUT_df, params["SUT_ranges"]["tot_intermediate_demand"])[:,sector_ndxs], dims=1))
	profit = max.(retval.g - tot_int_dmd - retval.W,zeros(ns))
	retval.μ = retval.g ./ (retval.g - profit .+ ϵ)

    if params["report-diagnostics"]
        qd = vec(sum(use_table, dims=2)) # Assign to a variable for clarity
        # Values by product
        write_vector_to_csv(joinpath(params["diagnostics_path"],"domestic_production.csv"), qs, "domestic production", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"imports.csv"),  M, "imports", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"margins.csv"),  margins, "margins", params["included_product_codes"])
		write_vector_to_csv(joinpath(params["diagnostics_path"],"imported_fraction.csv"),  retval.m_frac, "imported fraction", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"tot_intermediate_supply_non-energy_sectors.csv"), qd, "intermediate supply from non-energy sectors", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"tot_intermediate_supply_all_sectors.csv"),  tot_int_sup, "intermediate supply", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"exports.csv"),  retval.X, "exports", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"final_demand.csv"),  retval.F, "final demand", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"investment.csv"),  retval.I, "investment", params["included_product_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"tax_rate.csv"), retval.τd, "taxes on products", params["included_product_codes"])
        # Values by sector
        write_vector_to_csv(joinpath(params["diagnostics_path"],"sector_output.csv"), retval.g, "output", params["included_sector_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"tot_intermediate_demand_all_products.csv"),  tot_int_dmd, "intermediate demand", params["included_sector_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"wages.csv"),  retval.W, "wages", params["included_sector_codes"])
        write_vector_to_csv(joinpath(params["diagnostics_path"],"profit_margins.csv"), retval.μ, "profit margins", params["included_sector_codes"])
        # Matrices
        write_matrix_to_csv(joinpath(params["diagnostics_path"],"supply_fractions.csv"), retval.S, params["included_sector_codes"], params["included_product_codes"])
        write_matrix_to_csv(joinpath(params["diagnostics_path"],"demand_coefficients.csv"), retval.D, params["included_product_codes"], params["included_sector_codes"])
		# Write out nonenergy-energy link metric
        open(joinpath(params["diagnostics_path"],"nonenergy_energy_link_measure.txt"), "w") do fhndl
			R = 100 .* energy_nonenergy_link_measure(param_file)
			A_NE_metric_string = @sprintf("This value should be small: %.2f%%.", R)
			println(fhndl, "Measure of the significance to the economy of the supply of non-energy goods and services to the energy sector:")
		    println(fhndl, A_NE_metric_string)
		end
    end

    return retval, np, ns

end

"""
    prices_init(np::Int64, io::IOdata)

Initialize price structure.
"""
# All domestic prices are initialized to one
function prices_init(np::Int64, io::IOdata)
    return PriceData(ones(np), #pb
                     ones(np), #pd
                     ones(np), #pp
                     ones(np) + io.τd, #pw in domestic currency (converted using xr in Macro.jl)
                     1,1)
end

end

ext: jl

max_stars_repo_name: sei-international/LEAPMacro

max_issues_repo_name: sei-international/LEAPMacro

max_forks_repo_name: sei-international/LEAPMacro


text:
#struct PyTorchModel end

mutable struct PyTorchModel
    model::PyObject
    grad::Bool
end
PyTorchModel(m) = PyTorchModel(m, true)
convertmodel(m::PyTorchModel; device=nothing) = m
convertmodel(p::PyObject; device=nothing) = PyTorchModel(p.to(device))

struct PyTorchParams
    params::PyObject
end

function (m::PyTorchModel)(args...; kwargs...)
    if m.grad
        return m.model(args..., kwargs...)
    else
        @pywith torch.no_grad() begin
            return m.model(args..., kwargs...)
        end
    end
end

function Flux.params(m::PyTorchModel)
    return PyTorchParams(pybuiltin("list")(m.model.parameters()))
end

function Flux.trainmode!(m::PyTorchModel, mode = true)
    m.grad = mode
    m.model = m.model.train(mode)
end
function Flux.testmode!(m::PyTorchModel, mode = true)
    Flux.trainmode!(m, !mode)
end


text:
HERE = @__DIR__
DATA_PATH = joinpath(HERE, "../../test/data/campaspe/")

# Ensure dependent data and packages are available
using Statistics, DataFrames, CSV
using Distributed, BlackBoxOptim

using ModelParameters
using Graphs, MetaGraphs
using YAML
using Streamfall

import Random

Random.seed!(101)
network = YAML.load_file(joinpath(DATA_PATH, "campaspe_network.yml"))
sn = create_network("Example Network", network)

climate_data = CSV.File(joinpath(DATA_PATH, "climate/climate_historic.csv"),
                        comment="#",
                        dateformat="YYYY-mm-dd") |> DataFrame

hist_dam_levels = CSV.File(joinpath(DATA_PATH, "dam/historic_levels_for_fit.csv"),
                           dateformat="YYYY-mm-dd") |> DataFrame
hist_dam_releases = CSV.File(joinpath(DATA_PATH, "dam/historic_releases.csv"),
                             dateformat="YYYY-mm-dd") |> DataFrame

climate_data, hist_dam_levels, hist_dam_releases = Streamfall.align_time_frame(climate_data,
                                                                               hist_dam_levels,
                                                                               hist_dam_releases)

hist_data = Dict(
    "406000" => hist_dam_levels[:, "Dam Level [mAHD]"],
    "406000_extractions" => hist_dam_releases
)

climate = Climate(climate_data, "_rain", "_evap")

"""Calibrate current node."""
function obj_func(params, climate, sn, v_id, calib_data::Dict)

    this_node = sn[v_id]
    update_params!(this_node, params...)

    releases = calib_data["$(this_node.name)_extractions"]
    Streamfall.run_node!(sn, v_id, climate; extraction=releases)

    n_data = this_node.outflow
    h_data = calib_data[this_node.name]

    # Calculate score (NNSE; 0 to 1)
    NNSE = Streamfall.NNSE(h_data, n_data)

    # Switch fitness direction as we want to minimize
    score = 1.0 - NNSE

    # Other metrics are also available.
    # see `metrics.jl`
    # RMSE = Streamfall.RMSE(h_data, node_data)
    # score = RMSE

    # reset to clear stored values
    reset!(sn)

    return score
end

"""Example objective function when performance of current node is
dependent on the next node.
"""
function obj_func(params, climate, sn, v_id, next_vid, calib_data::Dict)

    this_node = sn[v_id]
    update_params!(this_node, params...)

    # Run next node which will run this node
    next_node = sn[next_vid]
    releases = calib_data["$(next_node.name)_extractions"]
    Streamfall.run_node!(sn, next_vid, climate; extraction=releases)

    # Alias data as necessary
    if next_node.name == "406000"
        n_data = next_node.level
        h_data = calib_data[next_node.name]
    elseif this_node.name == "406000"
        n_data = this_node.level
        h_data = calib_data[this_node.name]
    else
        n_data = this_node.outflow
        h_data = calib_data[this_node.name]
    end

    NNSE = Streamfall.NNSE(h_data, n_data)
    score = 1.0 - NNSE

    reset!(sn)

    return score
end

function alt_obj_func(params, climate, sn, v_id, next_vid, calib_data::Dict)
    this_node = sn[v_id]
    update_params!(this_node, params...)

    next_node = sn[next_vid]

    # Run next node (which will also run this node)
    releases = calib_data["$(next_node.name)_extractions"]
    Streamfall.run_node!(sn, next_vid, climate; extraction=releases)

    # Alias data as necessary
    if next_node.name == "406000"
        n_data = next_node.level
        h_data = calib_data[next_node.name]
    elseif this_node.name == "406000"
        n_data = this_node.level
        h_data = calib_data[this_node.name]
    else
        n_data = this_node.outflow
        h_data = calib_data[this_node.name]
    end

    split_NmKGE = Streamfall.naive_split_metric(h_data, n_data; n_members=365, metric=Streamfall.NmKGE, comb_method=mean)
    score = 1.0 - split_NmKGE

    reset!(sn)

    return score
end


text:
struct Entry
    plottype::PlotFunc
    mappings::Arguments
    attributes::Dict{Symbol, Any}
end

Entry(plottype::PlotFunc=Any, mappings=arguments(); attributes...) =
    Entry(plottype, mappings, Dict{Symbol, Any}(attributes))

Entry(mappings::Arguments; attributes...) = Entry(Any, mappings; attributes...)

function separate!(entry::Entry)
    discrete, _ = separate!(entry.mappings)
    return discrete => entry
end

function recombine!(discrete, entry::Entry)
    return Entry(
        entry.plottype,
        recombine!(discrete, entry.mappings),
        entry.attributes
    )
end

struct Entries
    entries::Vector{Entry}
    scales::Arguments
    labels::Arguments
end

Entries() = Entries(Entry[], arguments(), arguments())

function compute_axes_grid(fig, e::Entries; axis=NamedTuple())

    rowcol = (:row, :col)

    layout_scale, scales... = map((:layout, rowcol...)) do sym
        return get(e.scales, sym, nothing)
    end

    grid_size = map(scales, (first, last)) do scale, f
        isnothing(scale) || return maximum(scale.plot)
        isnothing(layout_scale) || return maximum(f, layout_scale.plot)
        return 1
    end

    axes_grid = map(CartesianIndices(grid_size)) do c
        type = get(axis, :type, Axis)
        options = Base.structdiff(axis, (; type))
        ax = type(fig[Tuple(c)...]; options...)
        return AxisEntries(ax, Entry[], e.scales, e.labels)
    end

    for entry in e.entries
        rows, cols = map(rowcol, scales, (first, last)) do sym, scale, f
            v = get(entry.mappings, sym, nothing)
            layout_v = get(entry.mappings, :layout, nothing)
            # without layout info, plot on all axes
            # all values in `v` and `layout_v` are equal
            isnothing(v) || return rescale(v[1:1], scale)
            isnothing(layout_v) || return map(f, rescale(layout_v[1:1], layout_scale))
            return 1:f(grid_size)
        end
        for i in rows, j in cols
            ae = axes_grid[i, j]
            push!(ae.entries, entry)
        end
    end

    return axes_grid

end

function AbstractPlotting.plot(entries::Entries; axis=NamedTuple(), figure=NamedTuple())
    fig = Figure(; figure...)
    grid = plot!(fig, entries; axis)
    return FigureGrid(fig, grid)
end

function AbstractPlotting.plot!(fig, entries::Entries; axis=NamedTuple())
    axes_grid = compute_axes_grid(fig, entries; axis)
    foreach(plot!, axes_grid)
    return axes_grid
end

"""
    AxisEntries(axis::Union{Axis, Nothing}, entries::Vector{Entry}, labels, scales)

Define all ingredients to make plots on an axis.
Each scale can be either a `CategoricalScale` (for discrete collections), such as
`CategoricalScale(["a", "b"], ["red", "blue"])`, or a function,
such as `log10`. Other scales may be supported in the future.
"""
struct AxisEntries
    axis::Union{Axis, Axis3}
    entries::Vector{Entry}
    scales::Arguments
    labels::Arguments
end

AbstractPlotting.Axis(ae::AxisEntries) = ae.axis
Entries(ae::AxisEntries) = Entries(ae.entries, ae.labels, ae.scales)

function prefix(i::Int, sym::Symbol)
    var = (:x, :y, :z)[i]
    return Symbol(var, sym)
end

# Slightly complex machinery to recombine stacked barplots
function mustbemerged(e::Entry)
    isbarplot = e.plottype <: BarPlot
    hasstack = :stack in keys(e.mappings.named) || :stack in keys(e.attributes)
    return isbarplot && hasstack
end

# Combine both entries as a unique entry with longer data
function stack!(e1::Entry, e2::Entry)
    p1, p2 = e1.plottype, e2.plottype
    m1, m2 = e1.mappings, e2.mappings
    a1, a2 = e1.attributes, e2.attributes
    l1, l2 = length(m1[1]), length(m2[1])
    assert_equal(p1, p2)
    for (k, v) in pairs(a1)
        assert_equal(v, a2[k])
    end
    mergewith!(m1, m2) do v1, v2
        long1 = size(v1) == () ? fill(v1[], l1) : v1
        long2 = size(v2) == () ? fill(v2[], l2) : v2
        return vcat(long1, long2)
    end
    return e1
end

function combine(entries::AbstractVector{Entry})
    combinedentries = Entry[]
    for entry in entries
        idx = findfirst(mustbemerged, combinedentries)
        if !isnothing(idx) && mustbemerged(entry)
            stack!(combinedentries[idx], entry)
        else
            push!(combinedentries, entry)
        end
    end
    return combinedentries
end

function AbstractPlotting.plot!(ae::AxisEntries)
    axis, entries, labels, scales = ae.axis, ae.entries, ae.labels, ae.scales
    for entry in combine(entries)
        plottype, mappings, attributes = entry.plottype, entry.mappings, entry.attributes
        trace = map(unwrap∘rescale, mappings, scales)
        positional, named = trace.positional, trace.named
        merge!(named, attributes)

        # Remove layout info
        for sym in [:col, :row, :layout]
            pop!(named, sym, nothing)
        end

        # Implement defaults
        for (key, val) in pairs(default_styles())
            get!(named, key, val)
        end

        # Set dodging and colorrange information
        dodge, colorscale = get(scales, :dodge, nothing), get(scales, :color, nothing)
        isa(dodge, CategoricalScale) && (named[:n_dodge] = maximum(dodge.plot))
        isa(colorscale, ContinuousScale) && (named[:colorrange] = colorscale.extrema)

        # Implement alpha transparency
        alpha = pop!(named, :alpha, nothing)
        color = get(named, :color, nothing)
        color isa Color && alpha isa Number && (named[:color] = (color, alpha))

        plot!(plottype, axis, positional...; named...)
    end
    # TODO: support log colorscale
    ndims = isaxis2d(ae) ? 2 : 3
    for i in 1:ndims
        label, scale = get(labels, i, nothing), get(scales, i, nothing)
        any(isnothing, (label, scale)) && continue
        axislabel, axisticks, axisscale = prefix.(i, (:label, :ticks, :scale))
        getproperty(axis, axisticks)[] = ticks(scale)
        getproperty(axis, axislabel)[] = string(label)
    end
    return axis
end

struct FigureGrid
    figure::Figure
    grid::Matrix{AxisEntries}
end

Base.show(io::IO, fg::FigureGrid) = show(io, fg.figure)
Base.show(io::IO, m::MIME, fg::FigureGrid) = show(io, m, fg.figure)
Base.show(io::IO, ::MIME"text/plain", fg::FigureGrid) = print(io, "FigureGrid()")

Base.showable(mime::MIME{M}, fg::FigureGrid) where {M} = showable(mime, fg.figure)

Base.display(fg::FigureGrid) = display(fg.figure)

function FileIO.save(filename::String, fg::FigureGrid; kwargs...)
    return FileIO.save(FileIO.query(filename), fg; kwargs...)
end

function FileIO.save(file::FileIO.Formatted, fg::FigureGrid; kwargs...)
    return FileIO.save(file, fg.figure; kwargs...)
end

to_tuple(fg) = (fg.figure, fg.grid)

Base.iterate(fg::FigureGrid) = iterate(to_tuple(fg))
Base.iterate(fg::FigureGrid, i) = iterate(to_tuple(fg), i)


text:
module Utilities

using Compat, Base.Meta

export Head, @H_str, issymbol, isexpr, parsefile

# Symbolic dispatch. #

immutable Head{S} end

if VERSION < v"0.4.0-dev"
    macro H_str(text)
        Expr(:(::), Expr(:call, :Union, [Head{symbol(part)} for part in split(text, ", ")]...))
    end
else
    macro H_str(text)
        Expr(:(::), Expr(:curly, :Union, [Head{symbol(part)} for part in split(text, ", ")]...))
    end
end

# File parsing and caching. #

"Internal version of the cache structure."
const CACHE_VER = v"0.0.3"

"Path to Docile's main cache folder."
const CACHE_DIR = normpath(joinpath(dirname(@__FILE__), "..", "cache"))

"The current versioned cache subdirectory set by ``CACHE_VER``"
const CACHE_CUR = joinpath(CACHE_DIR, string(CACHE_VER))

"""
Check cache directory for up to date version directory. Remove older version
directories if any are found. Should be called on module initialisation.
"""
function __init_cache__(cache = CACHE_DIR, current = CACHE_CUR)
    if !isdir(current)
        message("upgrading cache to $(CACHE_VER).")
        mkdir(current)
    end
    for each in readdir(cache)
        handle = joinpath(cache, each)
        if isdir(handle) && handle ≠ current
            rm(handle; recursive = true)
        end
    end
end

"""
Returns the cache path for a given file ``file``.
"""
path_id(file::AbstractString) = joinpath(
    CACHE_CUR,
    string(hash(file)),
    string(hash(mtime(file), hash(VERSION)))
    )

"""
Retrieve the ``Expr`` object from a Julia source file ``file``.

Caches the expression based on it's file name and ``mtime`` value.
"""
function parsefile(file::AbstractString)
    __init_cache__()
    cached = path_id(file)
    if isfile(cached)
        open(deserialize, cached)
    else
        # Read and parse ``file``.
        expr = try
            text = readall(file)
            out  = parse("begin $(text)\n end")
            isexpr(out, :incomplete) ? parse(text) : out
        catch
            warn("could not parse file: '$(file)'.")
            Expr(:block)
        end
        # Remove older versions and cache current ``expr``.
        dir = dirname(cached)
        isdir(dir) && rm(dir; recursive = true)
        mkdir(dir)
        open(cached, "w") do f
            serialize(f, expr)
        end
        expr
    end
end

# Misc. #

"""
Print a 'Docile'-formatted message to ``STDOUT``.
"""
message(msg::AbstractString) = print_with_color(:magenta, "Docile: ", msg, "\n")

"""
Is the module where a function/method is defined the same as ``mod``?
"""
samemodule(mod, def::Method)    = lsdfield(def, :module) == mod
samemodule(mod, func::Function) = lsdfield(func, :module) == mod
samemodule(mod, other)          = false

issymbol(::Symbol) = true
issymbol(::Any)    = false

"""
Path to Julia's base source code.
"""
basepath() = abspath(joinpath(JULIA_HOME, "..", "share", "julia", "base"))

"""
Convert a path to absolute. Relative paths are guessed to be from Julia ``/base``.
"""
expandpath(path) = normpath(isabspath(path) ? path : joinpath(basepath(), path))

# Compat.

lsdfield(x :: Function, f)    = isgeneric(x) ? lsdfield(methods(x), f) : lsdfield(x.code, f)
lsdfield(x :: Method, f)      = lsdfield(x.func, f)
lsdfield(x :: MethodTable, f) = lsdfield(x.defs, f)

lsdfield(x :: LambdaStaticData, f) = getfield(x, f)

end


text:
#using DelimitedFiles # To be able to load .txt files

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

"""
    conversion_factor_arcsec_to_sma

Multiplicative conversion factor from arcseconds to mpc.
"""
const conversion_factor_arcsec_to_sma = d_BH * 2*pi/(360*60*60)

"""
    arcsec_to_mpc(a::Float64)

Converts a semi-major axis given in arcseconds to an semi-major axis given in mpc.
"""
function arcsec_to_mpc(a::Float64)
    return a * conversion_factor_arcsec_to_sma
end

"""
    get_dataCluster()

Recovers the S-cluster data from INPUTCLUSTER.
"""
function get_dataCluster()
    dataCluster, headerCluster = readdlm(INPUTCLUSTER,header=true) # Reading the data file with the header
    return dataCluster # Output the entire data
end

"""
    dataCluster

Table of S-cluster data
"""
const dataCluster = get_dataCluster() # Reading the input data from the Cluster

"""
    nbCluster

Number of S-cluster stars used in the code
"""
const nbCluster = size(dataCluster)[1] # Number of stars in cluster

#################################################
# Getting the relevant physical quantities of the S-cluster
##################################################

"""
    tabaCluster

Table of S-cluster SMAs
"""
const tabaCluster = SVector{nbCluster,Float64}([arcsec_to_mpc(dataCluster[iCluster,2]) for iCluster=1:nbCluster]) # The sma a is stored in column 2 (mpc)

"""
    tabjCluster

Table of S-cluster reduced angular momenta.
"""
const tabjCluster = SVector{nbCluster,Float64}([ecc_to_j(dataCluster[iCluster,3]) for iCluster=1:nbCluster]) # The eccentricity e is stored in column 3

"""
    tabTCluster

Table of S-cluster main-sequence ages.
"""
const tabTCluster = SVector{nbCluster,Float64}([dataCluster[iCluster,4] for iCluster=1:nbCluster]) # The main-sequence age is stored in column 4 (in Myr)

"""
    tabdjCluster

Table of S-cluster SMA uncertainties.
"""
const tabdjCluster = SVector{nbCluster,Float64}([uncertainty_ecc_to_j(dataCluster[iCluster,3],dataCluster[iCluster,5]) for iCluster=1:nbCluster]) # The eccentricity uncertainty de is stored in column 5

"""
    tabdTpCluster

Table of S-cluster upper main-sequence age uncertainties.
"""
const tabdTpCluster = SVector{nbCluster,Float64}([dataCluster[iCluster,6] for iCluster=1:nbCluster]) # The age up-uncertainty dT_p is stored in column 6

"""
    tabdTmCluster

Table of S-cluster lower main-sequence age uncertainties.
"""
const tabdTmCluster = SVector{nbCluster,Float64}([dataCluster[iCluster,7] for iCluster=1:nbCluster]) # The age down-uncertainty dT_m is stored in column 7


text:
module ReinforcementLearningEnvironments

using ReinforcementLearningBase
using Random
using GR
using Requires
using IntervalSets
using Base.Threads:@spawn
using Markdown

const RLEnvs = ReinforcementLearningEnvironments
export RLEnvs

include("base.jl")
include("environments/environments.jl")
include("converters.jl")

# dynamic loading environments
function __init__()
    @require ArcadeLearningEnvironment = "b7f77d8d-088d-5e02-8ac0-89aab2acc977" include(
        "environments/3rd_party/atari.jl",
    )
    @require PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" include(
        "environments/3rd_party/gym.jl",
    )
    @require OpenSpiel = "ceb70bd2-fe3f-44f0-b81f-41608acaf2f2" include(
        "environments/3rd_party/open_spiel.jl",
    )
    @require SnakeGames = "34dccd9f-48d6-4445-aa0f-8c2e373b5429" include(
        "environments/3rd_party/snake.jl",
    )
    @require OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" include(
        "environments/3rd_party/AcrobotEnv.jl",
    )

end

end # module


text:
using LinearAlgebra
using jInv.Mesh
using Test

# setup 3D Mesh
nc = [5, 7, 8]
x0 = [0, 0, 0]
domain = [x0[1], 1, x0[2], 1, x0[3], 1]
domain2D = [x0[1], 1, x0[2], 1]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt = getTensorMesh3D(h1,h2,h3,x0)
Mr = getRegularMesh(domain,nc)

Mr2D = getRegularMesh(domain2D, [nc[1], nc[2]])

Meshes = (Mt, Mr, Mr2D);

for k=1:length(Meshes)

    M = Meshes[k]
    print("\ttesting boundary nodes for $(typeof(M))...")

    if M.dim==2

        X = getNodalGrid(M)[:,1];
        Y = getNodalGrid(M)[:,2];
        ib, iin = getBoundaryNodes(M);

        # Check if there are 0s in X and Y at boundary nodes
        @test length(X[ib]).-length(findall( x->(x < eps(Float32)), abs.(X[ib])))>0
        @test length(Y[ib]).-length(findall( x->(x < eps(Float32)), abs.(Y[ib])))>0

        # Check if there are 1s in X and Y at boundary nodes
        @test length(X[ib]).-length(findall( x->(x < eps(Float32)), abs.(X[ib].-1)))>0
        @test length(Y[ib]).-length(findall( x->(x < eps(Float32)), abs.(Y[ib].-1)))>0

        # Check that there aren't 0s in X, Y at inner nodes
        @test length(X[iin]).-length(findall(x->x!=0, X[iin]))==0
        @test length(Y[iin]).-length(findall(x->x!=0, Y[iin]))==0

        # Check that there aren't 1s in X, Y at inner nodes
        @test length(X[iin]).-length(findall(x->x!=1,X[iin]))==0
        @test length(Y[iin]).-length(findall(x->x!=1,Y[iin]))==0

    elseif M.dim==3

        X = getNodalGrid(M)[:,1];
        Y = getNodalGrid(M)[:,2];
        Z = getNodalGrid(M)[:,3];
        ib, iin = getBoundaryNodes(M);

        # Check if there are 0s in X, Y, and Z at boundary nodes
        @test length(X[ib])-length(findall( x->(x < eps(Float32)), abs.(X[ib])))>0
        @test length(Y[ib])-length(findall( x->(x < eps(Float32)), abs.(Y[ib])))>0
        @test length(Y[ib])-length(findall( x->(x < eps(Float32)), abs.(Z[ib])))>0

        # Check if there are 1s in X, Y, and Z at boundary nodes
        @test length(X[ib])-length(findall( x->(x < eps(Float32)), abs.(X[ib].-1)))>0
        @test length(Y[ib])-length(findall( x->(x < eps(Float32)), abs.(Y[ib].-1)))>0
        @test length(Y[ib])-length(findall( x->(x < eps(Float32)), abs.(Z[ib].-1)))>0

        # Check that there aren't 0s in X,Y,Z at inner nodes
        @test length(X[iin]).-length(findall(x->x!=0,X[iin]))==0
        @test length(Y[iin]).-length(findall(x->x!=0,Y[iin]))==0
        @test length(Z[iin]).-length(findall(x->x!=0,Z[iin]))==0

        # Check that there aren't 1s in X,Y,Z at inner nodes
        @test length(X[iin]).-length(findall(x->x!=1, X[iin]))==0
        @test length(Y[iin]).-length(findall(x->x!=1, Y[iin]))==0
        @test length(Z[iin]).-length(findall(x->x!=1, Z[iin]))==0

    end

    println("passed!\n")
end


text:
function ModelNet40(;mode="point_cloud")

    if mode == "point_cloud"
        return ModelNet40PCloud
    else
        error("Selected mode: $(mode) is not supported (Currently supported mode are {\"point_cloud\"}).")
    end
end

struct MN40DataPoint <: AbstractDataPoint
    idx::Int
    data::Union{PointCloud}
    ground_truth::UInt8
end

struct ModelNet40PCloud <: AbstractDataset
    root::String
    path::String #contains the path to dataset
    train::Bool
    transform #TODO Add type-assertion accordingly
    npoints::Int
    sampling #TODO Add type-assertion accordingly to include two possible option {"top", "uniform"}
    datapaths::Array
    length::Int
    classes_to_idx::Dict{String, UInt8}
    classes::Array{String,1}
end

function MN40_extract(path, npoints)
    pset = Array{Float32}(undef, npoints, 3)
    stream = open(path, "r")
    for i in 1:npoints
        pset[i, :] = map((x->parse(Float32, x)), split(readline(stream, keep=false), ",")[1:3])
    end
    return pset
end

function ModelNet40PCloud(;root::String=default_root, train::Bool=true, npoints::Int=1024, transform=nothing, sampling=nothing)
    _path = normpath(dataset("ModelNet40PCloud", root))
    train ? _split="train" : _split="test"
    cat_file = joinpath(_path, "modelnet40_shape_names.txt")

    classes = []
    classes_to_idx = []
    for (i, line) in enumerate(readlines(cat_file))
        push!(classes, line)
        push!(classes_to_idx, (line, convert(UInt8,i)))
    end
    classes_to_idx = Dict{String, UInt8}(classes_to_idx)

    shapeids = [line for line in readlines(joinpath(_path, "modelnet40_$(_split).txt"))]
    shape_names = [join(split(shapeids[i], "_")[1:end-1], "_") for i in 1:length(shapeids)]
    datapaths = [(shape_names[i], joinpath(_path, shape_names[i], (shapeids[i])*".txt")) for i in 1:length(shapeids)]
    _length = length(datapaths)
    ModelNet40PCloud(root, _path, train, transform, npoints, sampling, datapaths, _length, classes_to_idx, classes)
end

function Base.getindex(v::ModelNet40PCloud, idx::Int)
    cls = v.classes_to_idx[v.datapaths[idx][1]]
    pset = MN40_extract(v.datapaths[idx][2], v.npoints)
    return MN40DataPoint(idx, PointCloud(pset), cls)
end

Base.size(v::ModelNet40PCloud) = (v.length,)
Base.length(v::ModelNet40PCloud) = v.length


text:
using Plots, Colors

rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

function update_ukw_plot(plt, ukw, right_idx; connecting=false, connected=false)
    order = 0
    width = 7.5

    connected && (connecting = true)
    for i=(right_idx, ukw.mapping[right_idx])
        letter = Char(i+64)

        if connected && i == ukw.mapping[right_idx]
            plot!(plt, rectangle(1,1, order*width-1, i-1), color=:yellow)
            plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))
        end
        if right_idx == i || (right_idx !== nothing && i == ukw.mapping[right_idx] && connected)
            plot!(plt, rectangle(1,1, order*width+5, i-1), color=:yellow)
            plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))
        end

        # left lower/upper -
        if connecting && i == ukw.mapping[right_idx]
            plot!(plt, [order*width,order*width+0.3], [i-0.6, i-0.6], color=:yellow)
            plot!(plt, [order*width,order*width+0.3], [i-0.3, i-0.3], color=:yellow)
        end

        # right lower -
        if connecting && right_idx == i
            plot!(plt, [order*width+4.7,order*width+5], [i-0.6, i-0.6], color=:yellow)
        end

        # right upper -
        if connecting && i == ukw.mapping[right_idx]
            plot!(plt, [order*width+4.7,order*width+5], [i-0.3, i-0.3], color=:yellow)
        end

        y_pos_from = i-0.6
        y_pos_to = ukw.mapping[i]-0.6
        # right ro left
        if right_idx == i && (connecting || connected)
            plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:yellow)
        end

        # left to right
        if right_idx !== nothing && ukw.mapping[right_idx] == i && connected
            plot!(plt, [order*width+0.3, order*width+4.7], [y_pos_from+0.3, y_pos_from+0.3], color=:yellow)
        end
    end
    return plt
end

function plot_ukw(plt, ukw)
    order = 0
    width = 7.5

    plot!(plt, [order*width+2.5], [27], annotation=(order*width+2.5, 27, "UKW $(ukw.name)"))
    plot!(plt, rectangle(5,26, order*width, 0), color=:gray)

    for i=1:26
        letter = Char(i+64)
        plot!(plt, rectangle(1,1, order*width-1, i-1), color=:darkgray)
        plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))

        plot!(plt, rectangle(1,1, order*width+5, i-1), color=:darkgray)
        plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))

        # left lower/upper -
        plot!(plt, [order*width,order*width+0.3], [i-0.6, i-0.6], color=:black)
        plot!(plt, [order*width,order*width+0.3], [i-0.3, i-0.3], color=:black)

        # right lower -
        plot!(plt, [order*width+4.7,order*width+5], [i-0.6, i-0.6], color=:black)

        # right upper -
        plot!(plt, [order*width+4.7,order*width+5], [i-0.3, i-0.3], color=:black)

        y_pos_from = i-0.6
        y_pos_to = ukw.mapping[i]-0.6
        # right ro left
        plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:black)

        # left to right
        plot!(plt, [order*width+0.3, order*width+4.7], [y_pos_from+0.3, y_pos_from+0.3], color=RGB(0.3, 0.3, 0.3))
    end
    return plt
end

function update_plugboard_plot(plt, plugboard; right_idx=nothing, left_idx=nothing, connecting=false, connected=false)
    @assert right_idx !== nothing || left_idx !== nothing
    order = 4
    width = 7.5

    connected_to = 0
    if right_idx !== nothing
        connected_to = plugboard[right_idx]
        plot_range = (right_idx, connected_to)
    end
    finished = false
    if left_idx !== nothing
        right_idx    = plugboard[left_idx]
        connected_to = left_idx
        finished = true
        plot_range = (right_idx, connected_to)
    end

    connected && (connecting = true)
    for i in (connected_to, right_idx)
        letter = Char(i+64)

        # left box and letter
        if (connected || finished) && connected_to == i
            plot!(plt, rectangle(1,1, order*width-1, i-1), color=:yellow)
            plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))
        end
        # right box and letter
        if right_idx == i && finished && connected
            plot!(plt, rectangle(1,1, order*width+5, i-1), color=:green)
            plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))
        elseif right_idx == i && !finished
            plot!(plt, rectangle(1,1, order*width+5, i-1), color=:yellow)
            plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))
        end

        y_pos_from = i-0.5
        y_pos_to = plugboard[i]-0.5
        if connecting && right_idx == i
            plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:yellow)
            plot!(plt, [order*width+4.7,order*width+5], [i-0.5, i-0.5], color=:yellow)
        end

        if connecting && connected_to == i
            plot!(plt, [order*width,order*width+0.3], [i-0.5, i-0.5], color=:yellow)
        end
    end
    return plt
end

function plot_plugboard(plt, plugboard)
    order = 4
    width = 7.5

    plot!(plt, [order*width+2.5], [27], annotation=(order*width+2.5, 27, "Steckbrett"))
    plot!(plt, rectangle(5,26, order*width, 0), color=:gray)

    for i=1:26
        letter = Char(i+64)
        plot!(plt, rectangle(1,1, order*width-1, i-1), color=:darkgray)
        plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))

        plot!(plt, rectangle(1,1, order*width+5, i-1), color=:darkgray)
        plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))

        y_pos_from = i-0.5
        y_pos_to = plugboard[i]-0.5
        plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:black)
        plot!(plt, [order*width+4.7,order*width+5], [i-0.5, i-0.5], color=:black)

        plot!(plt, [order*width,order*width+0.3], [i-0.5, i-0.5], color=:black)
    end
    return plt
end

function update_rotor_plot(plt::Plots.Plot, rotor::Rotor; right_idx = nothing, left_idx = nothing, connecting=false, connected=false)
    @assert right_idx !== nothing || left_idx !== nothing

    order = rotor.order
    width = 7.5

    connected_to = 0
    backwards = false
    if right_idx !== nothing
        left_idx = index_connected_to(rotor, right_idx)
        plot_range = (right_idx, left_idx)
    else # left_idx
        right_idx = index_connected_to(rotor, left_idx; backward=true)
        plot_range = (right_idx, left_idx)
        backwards = true
    end
    letter_to_idx = Dict{Char, Int}()
    for i=1:26
        letter = Char((25+i+rotor.position-1)%26 +65)
        letter_to_idx[letter] = i
    end

    for i = (left_idx, right_idx)
        letter = Char((25+i+rotor.position-1)%26 +65)
        letter_to_idx[letter] = i
        # left box
        if left_idx == i && (backwards || connected)
            plot!(plt, rectangle(1,1, order*width-1, i-1), color=:yellow)
            plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))
        end
        # right box
        if right_idx == i && (!backwards || connected)
            plot!(plt, rectangle(1,1, order*width+5, i-1), color=:yellow)
            plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))
        end

        # left side small -
        if left_idx == i
            plot!(plt, [order*width,order*width+0.3], [i-0.5, i-0.5], color=:yellow)
        end

        # left side small -
        if right_idx == i
            plot!(plt, [order*width+4.7,order*width+5], [i-0.5, i-0.5], color=:yellow)
        end
    end

    # plot connections
    for i=1:26
        y_pos_from = letter_to_idx[Char(i+64)]-0.5
        y_pos_to   = letter_to_idx[Char(rotor.mapping[i]+64)]-0.5

        if right_idx == letter_to_idx[Char(i+64)]
            plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:yellow)
        end
    end
    return plt
end

function plot_rotor(plt::Plots.Plot, rotor::Rotor)
    order = rotor.order
    width = 7.5

    plot!(plt, [order*width+2.5], [27], annotation=(order*width+2.5, 27, "Rotor $(rotor.name)"))
    plot!(plt, rectangle(5,26, order*width, 0), color=:gray)
    letter_to_idx = Dict{Char, Int}()
    for i=1:26
        letter = Char((25+i+rotor.position-1)%26 +65)
        letter_to_idx[letter] = i
    end
    for i=1:26
        letter = Char((25+i+rotor.position-1)%26 +65)
        letter_to_idx[letter] = i
        # left box
        plot!(plt, rectangle(1,1, order*width-1, i-1), color=:darkgray)
        plot!(plt, [order*width-1], [i-1], annotation=(order*width-0.5, i-0.5, letter))

        # right box
        plot!(plt, rectangle(1,1, order*width+5, i-1), color=:darkgray)
        plot!(plt, [order*width+5], [i-1], annotation=(order*width+5.5, i-0.5, letter))

        # left side small -
        plot!(plt, [order*width,order*width+0.3], [i-0.5, i-0.5], color=:black)

        # left side small -
        plot!(plt, [order*width+4.7,order*width+5], [i-0.5, i-0.5], color=:black)
    end
    # plot connections
    for i=1:26
        y_pos_from = letter_to_idx[Char(i+64)]-0.5
        y_pos_to   = letter_to_idx[Char(rotor.mapping[i]+64)]-0.5

        plot!(plt, [order*width+4.7, order*width+0.3], [y_pos_from, y_pos_to], color=:black)
    end

    return plt
end

function get_enigma_plot(enigma::EnigmaMachine)
    plt = plot(; size=(1000,800), grid=false, xaxis=false, yaxis=false, aspectratio=:equal, legend=false)

    plot_plugboard(plt, enigma.plugboard)
    for rotor in enigma.rotors
        plot_rotor(plt, rotor)
    end
    plot_ukw(plt, enigma.ukw)
    return plt
end

function get_enigma_decode_plots!(enigma::EnigmaMachine, letter::Char; pngs=false)
    plts = Plots.Plot[]
    plt = get_enigma_plot(enigma)
    push!(plts, deepcopy(plt))
    step_rotors!(enigma)
    plt = get_enigma_plot(enigma)
    push!(plts, deepcopy(plt))

    letter_idx = Int(letter-64)
    update_plugboard_plot(plt, enigma.plugboard; right_idx=letter_idx)
    push!(plts, deepcopy(plt))
    update_plugboard_plot(plt, enigma.plugboard; right_idx=letter_idx, connecting=true)
    push!(plts, deepcopy(plt))
    update_plugboard_plot(plt, enigma.plugboard; right_idx=letter_idx, connected=true)
    push!(plts, deepcopy(plt))
    pngs && png(plt, "visualizations/plugboard_in")

    plugboard_out_idx = enigma.plugboard[letter_idx]
    right_idx = plugboard_out_idx
    for r=3:-1:1
        update_rotor_plot(plt, enigma.rotors[r]; right_idx=right_idx)
        push!(plts, deepcopy(plt))
        update_rotor_plot(plt, enigma.rotors[r]; right_idx=right_idx, connecting=true)
        push!(plts, deepcopy(plt))
        update_rotor_plot(plt, enigma.rotors[r]; right_idx=right_idx, connected=true)
        push!(plts, deepcopy(plt))
        right_idx = index_connected_to(enigma.rotors[r], right_idx)
        pngs && png(plt, "visualizations/rotor_$(r)_in")
    end
    update_ukw_plot(plt, enigma.ukw, right_idx)
    push!(plts, deepcopy(plt))
    update_ukw_plot(plt, enigma.ukw, right_idx, connecting=true)
    push!(plts, deepcopy(plt))
    update_ukw_plot(plt, enigma.ukw, right_idx, connected=true)
    push!(plts, deepcopy(plt))
    left_idx = enigma.ukw.mapping[right_idx]
    pngs && png(plt, "visualizations/ukw")
    for r=1:3
        update_rotor_plot(plt, enigma.rotors[r]; left_idx=left_idx)
        push!(plts, deepcopy(plt))
        update_rotor_plot(plt, enigma.rotors[r]; left_idx=left_idx, connecting=true)
        push!(plts, deepcopy(plt))
        update_rotor_plot(plt, enigma.rotors[r]; left_idx=left_idx, connected=true)
        push!(plts, deepcopy(plt))
        left_idx = index_connected_to(enigma.rotors[r], left_idx; backward=true)
        pngs && png(plt, "visualizations/rotor_$(r)_bw")
    end
    update_plugboard_plot(plt, enigma.plugboard; left_idx=left_idx)
    push!(plts, deepcopy(plt))
    update_plugboard_plot(plt, enigma.plugboard; left_idx=left_idx, connecting=true)
    push!(plts, deepcopy(plt))
    update_plugboard_plot(plt, enigma.plugboard; left_idx=left_idx, connected=true)
    push!(plts, deepcopy(plt))
    pngs && png(plt_plug, "visualizations/plugboard_out")

    return plts
end

function animate_plots(plts::Vector{Plots.Plot}, fname::String; end_extra=5)
    anim = Animation()
    for plt in plts
        plot(plt)
        frame(anim)
    end
    for i=1:end_extra
        plot(plts[end])
        frame(anim)
    end
    gif(anim, "$fname.gif", fps=4)
end

export get_enigma_plot, get_enigma_decode_plots!, animate_plots


text:
# Automatically generated by gen_leap_seconds.jl, do not edit!

const LS_EPOCHS = [
    41317,
    41499,
    41683,
    42048,
    42413,
    42778,
    43144,
    43509,
    43874,
    44239,
    44786,
    45151,
    45516,
    46247,
    47161,
    47892,
    48257,
    48804,
    49169,
    49534,
    50083,
    50630,
    51179,
    53736,
    54832,
    56109,
    57204,
    57754,
]

const LEAP_SECONDS = [
    10.0,
    11.0,
    12.0,
    13.0,
    14.0,
    15.0,
    16.0,
    17.0,
    18.0,
    19.0,
    20.0,
    21.0,
    22.0,
    23.0,
    24.0,
    25.0,
    26.0,
    27.0,
    28.0,
    29.0,
    30.0,
    31.0,
    32.0,
    33.0,
    34.0,
    35.0,
    36.0,
    37.0,
]


text:
import IndexedTables: pkeynames, excludecols, pkeys, reindex
import Dagger: dsort_chunks

export reindex, rechunk

using StatsBase

function reindex(t::DDataset, by=pkeynames(t), select=excludecols(t, by); kwargs...)
    @noinline function _rechunk(c)
        reindex(c, by, select; kwargs...)
    end
    fromchunks(delayedmap(_rechunk, t.chunks))
end

"""
`rechunk(t::Union{DNDSparse, DNDSparse}[, by[, select]]; <options>)`

Reindex and sort a distributed dataset by keys selected by `by`.

Optionally `select` specifies which non-indexed fields are kept. By default this is all fields not mentioned in `by` for Table and the value columns for NDSparse.

# Options:
- `chunks` -- how to distribute the data. This can be:
    1. An integer -- number of chunks to create
    2. An vector of `k` integers -- number of elements in each of the `k` chunks. `sum(k)` must be same as `length(t)`
    3. The distribution of another array. i.e. `vec.subdomains` where `vec` is a distributed array.
- `merge::Function` -- a function which merges two sub-table or sub-ndsparse into one NDSparse. They may have overlaps in their indices.
- `splitters::AbstractVector` -- specify keys to split by. To create `n` chunks you would need to pass `n-1` splitters and also the `chunks=n` option.
- `chunks_sorted::Bool` -- are the chunks sorted locally? If true, this skips sorting or re-indexing them.
- `affinities::Vector{<:Integer}` -- which processes (Int pid) should each output chunk be created on. If unspecified all workers are used.
- `closed::Bool` -- if true, the same key will not be present in multiple chunks (although sorted). `true` by default.
- `nsamples::Integer` -- number of keys to randomly sample from each chunk to estimate splitters in the sorting process. (See [samplesort](https://en.wikipedia.org/wiki/Samplesort)). Defaults to 2000.
- `batchsize::Integer` -- how many chunks at a time from the input should be loaded into memory at any given time. This will essentially sort in batches of `batchsize` chunks.
"""
function rechunk(dt::DDataset,
                 by=pkeynames(dt),
                 select=dt isa DNextTable ? excludecols(dt, by) : valuenames(dt);
                 merge=_merge,
                 splitters=nothing,
                 chunks_presorted=false,
                 affinities=workers(),
                 chunks=nworkers(),
                 closed=true,
                 sortchunks=true,
                 nsamples=2000,
                 batchsize=nworkers())

    if sortchunks
        perm = sortperm(dt.domains, by=first)
        cs = dt.chunks[perm]
    else
        cs = dt.chunks
    end

    function sortandsample(data, nsamples, presorted)
        r = sample(1:length(data), min(length(data), nsamples),
                   replace=false, ordered=true)

        sorted = !presorted ? reindex(data, by, select) : data
        chunk = !presorted ? tochunk(sorted) : nothing

        (chunk, pkeys(sorted)[r])
    end

    dsort_chunks(cs, chunks, nsamples,
                 batchsize=batchsize,
                 sortandsample=sortandsample,
                 affinities=affinities,
                 splitters=splitters,
                 chunks_presorted=chunks_presorted,
                 merge=merge,
                 by=pkeys,
                 sub=subtable) |> fromchunks

end

### Permutedims

function Base.permutedims(t::DNDSparse{K,V}, p::AbstractVector) where {K,V}
    if !(length(p) == ndims(t) && isperm(p))
        throw(ArgumentError("argument to permutedims must be a valid permutation"))
    end

    permuteintv(intv,d) = Interval(first(intv)[d],
                                   last(intv)[d])
    idxs = map(t.domains) do dmn
        IndexSpace(permuteintv(dmn.interval, p),
                   permuteintv(dmn.boundingrect, p),
                   dmn.nrows,
                  )
    end

    chunks = map(delayed(c -> permutedims(c, p)), t.chunks)
    t1 = DNDSparse{eltype(idxs[1]), V}(idxs, chunks)

    cache_thunks(rechunk(t1))
end


text:
module EchogramProcessing

export meanfilter, medianfilter, meanfilterdb, removeholes

using Images
using ImageFiltering
using EchogramUtils
using Statistics

function meanfilter(A, height, width)
    kernel  = Array{Float64}(undef, height, width)
    kernel .= 1 / (width * height)
    kernel = centered(kernel)
    imfilter(A, kernel, Fill(0))
end

function medianfilter(A, height, width)
    mapwindow(median, A, (height, width))
end

"""
    meanfilterdb(db, height, width)

Remember, the Average of the Log is Not the Log of the Average,
so averaging dB requires conversion to and from the linear domain.
"""
function meanfilterdb(db, height, width)
    a = db2pow.(db)
    b = meanfilter(a, height, width)
    pow2db.(b)
end

"""
    removeholes(A; min=9)

Remove holes in binary array A.

e.g.

```
    removeholes([1 1 1; 1 0 1; 1 1 1])

    3×3 Array{Int64,2}:
     1  1  1
     1  1  1
     1  1  1
```
"""
function removeholes(A; min=9)
    c = copy(A)
    A= 1 .- A # invert

    # Label background regions
    labels = label_components(A)
    for i in 1:maximum(labels)
        b = (labels .== i)
        # Remove small regions
        if sum(b) <= min
            c[b] .= 1
        end
    end
    return c
end

end # module

max_stars_repo_name: EchoJulia/EchogramProcessing

max_issues_repo_name: EchoJulia/EchogramProcessing

max_forks_repo_name: EchoJulia/EchogramProcessing


text:
function agregar_sencillo!(array_velocidad, array_flujo, array_densidad, i, t, T)

    v = ( array_densidad[i, t] != 0 ? array_flujo[i, t]/array_densidad[i, t] : 0 )
    array_velocidad[i, t] += v/T
end

function agregar_Frecuencias!(array2, v, i, t, j, x_seccion, T; tipo::Int64 = 2)

  if v.velocidad > 0
    if j-v.velocidad < x_seccion
      if v.tipo == tipo
        array2[i, t] += 1/T
      end
    end
  end
end

function agregar!(i, j, t, x_seccion, v, array, array1, array2, T)

    if v.velocidad > 0
        if j-v.velocidad < x_seccion
            array[i, t] += 1/T

            if v.tipo == 1
                array1[i, t] += 1/T
            end
            if v.tipo == 2
                array2[i, t] += 1/T
            end
        end
    end
end

function Medir_Frecuencias!(carretera, t, Secciones, T, flujo_local2, densidad_local2, flujo_local1, Diagrama_Transicion, j_in, p)

    for (i, Seccion) in enumerate(Secciones)
      for j = Seccion+4:-1:Seccion
            agregar_Frecuencias!(flujo_local2, carretera[j], i, t, j, Seccion, T)
            agregar_Frecuencias!(densidad_local2, carretera[j], i, t, j, Seccion, T*carretera[j].velocidad)
            agregar_Frecuencias!(flujo_local1, carretera[j], i, t, j, Seccion, T, tipo=1)

            if carretera[j].velocidad == 0 && (carretera[j].tipo == 1 || carretera[j].tipo == 2 )
                    Diagrama_Transicion[i, j_in, p] = 0
            end
      end
    end
end

#################################### Mediciones ######################################

function Medir!(C, t, Secciones, T, flujo_local, densidad_local ,
               velocidad_local_promedio, flujo_local1, densidad_local1 ,velocidad_local_promedio1,
               flujo_local2, densidad_local2 , velocidad_local_promedio2)

    #T = 1.  T es el periodo de tiempo por el cual se mide

    for (i, x_seccion) in enumerate(Secciones)

        # Medicion del flujo local
        for j = x_seccion+4:-1:x_seccion
            agregar!(i, j, t, x_seccion, C.carretera[j], flujo_local, flujo_local1, flujo_local2, T)
        end

        # Medicion de la densidad local
        for j = x_seccion+4:-1:x_seccion
            agregar!(i, j, t, x_seccion, C.carretera[j], densidad_local, densidad_local1, densidad_local2,
                    T*C.carretera[j].velocidad)
        end

        # Medicion de la velocidad promedio espacial

        agregar_sencillo!( velocidad_local_promedio, flujo_local, densidad_local, i, t, T )
        agregar_sencillo!( velocidad_local_promedio1, flujo_local1, densidad_local1, i, t, T )
        agregar_sencillo!( velocidad_local_promedio2, flujo_local2, densidad_local2, i, t, T )
    end
end

function ET!(X, T, carretera, t)

    for v in carretera
        if v.tipo == 1 || v.tipo == 2
            push!(X[v.num], v.posicion)
            push!(T[v.num], t)
        end
    end
end

function ET_tipo!(X, T, carretera, t, tipo::Int64)

    for v in carretera
        if v.tipo == tipo
            push!(X[v.num], v.posicion)
            push!(T[v.num], t)
        end
    end
end

function tiempos_salida!(Tiempos_salida::Array{Int64, 1}, Tiempos_salida1::Array{Int64, 1}, Tiempos_salida2::Array{Int64, 1},
                        T0::Array{Any,1}, T1::Array{Any,1}, T2::Array{Any,1},
                        X::Array{Any,1}, X1::Array{Any,1}, X2::Array{Any,1},
                        N::Int64 = 3862)

    for i in 1:length(X)
        if length(X[i]) != 0 && X[i][end] > N-6
            push!(Tiempos_salida, T0[i][end]-T0[i][1])
        end
        if length(X1[i]) != 0 && X1[i][end] > N-6
            push!(Tiempos_salida1, T1[i][end]-T1[i][1])
        end
        if length(X2[i]) != 0 && X2[i][end] > N-6
            push!(Tiempos_salida2, T2[i][end]-T2[i][1])
        end
    end
end

function tiempos_iniciales!(Tiempos_iniciales1::Array{Int64,1}, Tiempos_iniciales2::Array{Int64, 1},
                            T1::Array{Any, 1}, T2::Array{Any, 1}, X1::Array{Any, 1}, X2::Array{Any, 1}, N::Int64 = 3862)

  for i in 1:length(T1)
      if length(T1[i]) != 0 && X1[i][end] > N-6
          push!(Tiempos_iniciales1, T1[i][1])
      end
  end
  for i = 1:length(T2)
      if length(T2[i]) != 0 && X2[i][end] > N-6
          push!(Tiempos_iniciales2, T2[i][1])
      end
  end
end

function Tiempos_promedio(Tf::Int64, Tiempos_salida, Tiempos_iniciales)

    Ts_promedio = Any[]

    h = int64(ceil(Tf/900))

    for i = 1:h
        push!(Ts_promedio, Float64[])
    end

    for (l, lim) in enumerate(linspace(900, Tf, h))
        for i = 1:length(Tiempos_salida)
            if Tiempos_iniciales[i] <= lim
                push!(Ts_promedio[l], Tiempos_salida[i])
                Tiempos_iniciales[i] = Tf+1
            end
        end
    end

    Tiempos_salida_promedio = Float64[0.]
    Desviaciones_salida_promedio = Float64[0.]

    for i = 1:length(Ts_promedio)
        push!(Tiempos_salida_promedio, mean(Ts_promedio[i]))
        push!(Desviaciones_salida_promedio, std(Ts_promedio[i]))
    end

    Ts_promedio = 0
    return Tiempos_salida_promedio, Desviaciones_salida_promedio
end

function Medidas_DF!(Flujos_secciones_promedio::Array{Float64, 3}, D_Flujos_secciones_promedio::Array{Float64, 3}, densidad::Array{Float64, 3},
                     flujo::Array{Float64, 3}, secciones::Array{Int64, 1})

    for k = 1:size(Flujos_secciones_promedio)[3]
        for (i, s) in enumerate(secciones[1:end-1])
            for j =s:secciones[i+1]
                copia_d = copy([de for de in sub(densidad, j, :, k)][1:end-1])
                copia_f = copy([fl for fl in sub(flujo, j, :, k)][1:end-1])
                flujos_promedio, Desviaciones_flujos_promedio = Flujos_promedio(copia_f, copia_d)
                for l = 1:size(Flujos_secciones_promedio)[2]
                    Flujos_secciones_promedio[i, l, k] += flujos_promedio[l]/(secciones[i+1]-s+1)
                    D_Flujos_secciones_promedio[i, l, k] += Desviaciones_flujos_promedio[l]/(secciones[i+1]-s+1)
                end
                copia_d = copia_f = 0
            end
        end
    end
end

function Medidas_ET!(Medidas_secciones_promedio::Array{Float64, 3}, D_Medidas_secciones_promedio::Array{Float64, 3}, Medidas::Array{Float64, 3},
                     Tf::Int64, T::Int64, secciones::Array{Int64, 1})

    for k = 1:size(Medidas_secciones_promedio)[3]
        for (i, s) in enumerate(secciones[1:end-1])
            for j =s:secciones[i+1]
                 M_promedio_T, D_M_promedio_T = Flujos_promedio_T(Tf, T,
                                                    Float64[m for m in sub(Medidas, j, :, k)])
                 for l = 1:length(M_promedio_T)
                    Medidas_secciones_promedio[i, l, k] += M_promedio_T[l]/(secciones[i+1]-s+1)
                    D_Medidas_secciones_promedio[i, l, k] += D_M_promedio_T[l]/(secciones[i+1]-s+1)
                 end
                 M_promedio_T = D_M_promedio_T = 0
            end
        end
    end
end

function Flujos_promedio_T(Tf::Int64, T::Int64, Flujos::Array{Float64, 1})

    F_promedio_T = Float64[0.]
    D_F_promedio_T = Float64[0.]

    h = int(ceil(Tf/900))
    m = 900/T

    for i = 0:h-1
        push!(F_promedio_T, mean(Flujos[i*m+1:m*(i+1)]))
        push!(D_F_promedio_T, std(Flujos[i*m+1:m*(i+1)]))

    end

    return F_promedio_T, D_F_promedio_T
end

function Flujos_promedio(Flujos, Densidades)

    F_promedio = Any[]

    h = 20

    for i = 1:h
        push!(F_promedio, Float64[])
    end

    for (l, lim) in enumerate(linspace(0.05, 1, h))
        for i = 1:length(Flujos)
            if Densidades[i] <= lim
                push!(F_promedio[l], Flujos[i])
                Densidades[i] = 1.1
            end
        end
    end

    flujos_promedio = Float64[0.]
    Desviaciones_flujos_promedio = Float64[0.]

    for i = 1:length(F_promedio)
        push!(flujos_promedio, mean(F_promedio[i]))
        push!(Desviaciones_flujos_promedio, std(F_promedio[i]))
    end

    F_promedio = 0
    return flujos_promedio, Desviaciones_flujos_promedio
end


text:
parse(lines) = split.(lines, "")

# read_data() = parse.(Int, readlines("day11.txt"))
read_data() = parse(readlines("day11.txt"))
# read_data() = parse(readlines("day11_sample1.txt"))

is_occupied(data,x,y) = data[x][y] == "#"
Base.isempty(data,x,y) = data[x][y] == "L"

rows(data) = length(data)
cols(data) = length(data[1])

function occupied_nearby(data, x, y)
    let rows = rows(data), cols = cols(data)
        n  = ((x - 1) >= 1 && is_occupied(data, x - 1, y))
        nw = ((x - 1) >= 1 && (y - 1) >= 1 && is_occupied(data, x - 1, y - 1))
        w  = ((y - 1) >= 1 && is_occupied(data, x, y - 1))
        sw = ((y - 1) >= 1 && (x + 1) <= rows && is_occupied(data, x + 1, y - 1))
        s  = ((x + 1) <= rows && is_occupied(data, x + 1, y))
        se = ((x + 1) <= rows && (y + 1) <= cols && is_occupied(data, x + 1, y + 1))
        e  = ((y + 1) <= cols && is_occupied(data, x, y + 1))
        ne = ((x - 1) >= 1 && (y + 1) <= cols && is_occupied(data, x - 1, y + 1))
        return n + nw + w + sw + s + se + e + ne
    end
end

function printdata(data)
    for r in 1:rows(data)
        for c in 1:cols(data)
            print(data[r][c])
        end
        println()
    end
end

function part1()
    data = read_data()
    h = UInt(0)
    cnt = 1
    while hash(data) != h
        h = hash(data)
        rule1 = [occupied_nearby(data,r,c) == 0 && isempty(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        rule2 = [occupied_nearby(data,r,c) >= 4 && is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        for r in 1:rows(data), c in 1:cols(data)
            if isempty(data,r,c) && rule1[r,c]
                data[r][c] = "#"
            elseif is_occupied(data,r,c) && rule2[r,c]
                data[r][c] = "L"
            end
        end
        cnt += 1
    end
    count(is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data))
end

# Part 2

function visible_occupied_seats(data, x, y)
    function trace(x, y, dx, dy, first=false)
        (x < 1 || y < 1 || x > rows(data) || y > cols(data)) && return 0
        if !first && isempty(data, x, y)
            # @info "found empty seat $x $y"
            return 0   # empty seat obstructs occupied seat
        elseif !first && is_occupied(data, x, y)
            # @info "found occupied seat $x $y"
            return 1
        else
            return trace(x + dx, y + dy, dx, dy, false)  # keep looking
        end
    end
    # trace(x, y, -1, 0, true)
    sum(trace(x, y, dx, dy, true) for dx in [-1, 0, 1], dy in [-1, 0, 1] if !(dx == 0 && dy == 0))
end

function part2()
    data = read_data()
    h = UInt(0)
    cnt = 1
    while hash(data) != h
        h = hash(data)
        rule1 = [visible_occupied_seats(data,r,c) == 0 && isempty(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        rule2 = [visible_occupied_seats(data,r,c) >= 5 && is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        for r in 1:rows(data), c in 1:cols(data)
            if isempty(data,r,c) && rule1[r,c]
                data[r][c] = "#"
            elseif is_occupied(data,r,c) && rule2[r,c]
                data[r][c] = "L"
            end
        end
        # printdata(data)
        cnt += 1
    end
    count(is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data))
end

# Animation

using Plots, ColorSchemes

mapper(x) = x == "L" ? 0.5 : (x == "#" ? 1.0 : 0)

create_matrix(data) = [mapper(data[x][y]) for x in 1:rows(data), y in 1:cols(data)]

make_heatmap(A, i, scale = 5) = heatmap(A, size = size(A).*scale,
    legend=false, ticks=:none, c=:winter, title="AoC Day 11 - Iteration $i")

function part1_animated()
    data = read_data()
    h = UInt(0)
    cnt = 1
    anim = Animation()
    while hash(data) != h
        h = hash(data)
        make_heatmap(create_matrix(data), cnt, 5)
        frame(anim)
        rule1 = [occupied_nearby(data,r,c) == 0 && isempty(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        rule2 = [occupied_nearby(data,r,c) >= 4 && is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data)]
        for r in 1:rows(data), c in 1:cols(data)
            if isempty(data,r,c) && rule1[r,c]
                data[r][c] = "#"
            elseif is_occupied(data,r,c) && rule2[r,c]
                data[r][c] = "L"
            end
        end
        cnt += 1
    end
    gif(anim, "day11_anim.gif", fps = 10)
    count(is_occupied(data,r,c) for r in 1:rows(data), c in 1:cols(data))
end

# Stencil idea

# using OffsetArrays

# function parse_with_pad(file)
#     lines = readlines(file)
#     rows = length(lines)
#     cols = length(lines[1])
#     A = fill('.', (rows+2, cols+2))
#     for r in 1:rows
#         for c in 1:cols
#             A[r+1,c+1] = lines[r][c]
#         end
#     end
#     return OffsetArray(A, 0:rows+1, 0:cols+1)
# end

# function occupied_nearby_using_stencil(data, x, y)
#     stencil = [(-1, -1), (1, 1), (-1, 1), (1, -1), (0, 1), (1, 0), (0, -1), (-1, 0)]
#     delta = CartesianIndex.(stencil)
#     origin = CartesianIndex(x, y)
#     return sum(data[origin + d] == '#' for d in delta)
# end


text:
@testset "testing Abstract Storage With DC - PF" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.Storage, PSI.AbstractStorageForm, PM.DCPlosslessForm, sys5b_storage, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 96
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24
end

@testset "testing Abstract Storage With AC - PF" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.Storage, PSI.AbstractStorageForm, PM.StandardACPForm, sys5b_storage, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 120
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24
end

@testset "testing Basic Storage With DC - PF" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.Storage, PSI.BookKeepingModel, PM.DCPlosslessForm, sys5b_storage, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 96
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24
    end

@testset "testing Basic Storage With AC - PF" begin
    ps_model = PSI._canonical_model_init(bus_numbers, nothing, PM.AbstractPowerFormulation, time_range)
    PSI.construct_device!(ps_model, PSY.Storage, PSI.BookKeepingModel, PM.StandardACPForm, sys5b_storage, time_range);
    @test JuMP.num_variables(ps_model.JuMPmodel) == 120
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.LessThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.GreaterThan{Float64}) == 48
    @test JuMP.num_constraints(ps_model.JuMPmodel,GenericAffExpr{Float64,VariableRef},MOI.EqualTo{Float64}) == 24
 end