text:
# top type - unifying a type with top yeilds the type
struct Top end

# get the type from a named tuple, given a name
get_type(NT, k) = hasfield(NT, k) ? fieldtype(NT, k) : Nothing

# unify two types to a single type
unify(a, b) = unify(b, a)
unify(a::Type{T}, b::Type{S}) where {T,S} = Base.promote_typejoin(T, S)
unify(a::Type{T}, b::Type{S}) where {T,S<:T} = T
unify(a::Type{Top}, b::Type{T}) where {T} = T

function unify(
    a::Type{NamedTuple{A,T}},
    b::Type{NamedTuple{B,S}},
) where {A,T<:Tuple,B,S<:Tuple}
    c = Dict()
    for (k, v) in zip(A, fieldtypes(a))
        c[k] = unify(v, get_type(b, k))
    end

    for (k, v) in zip(B, fieldtypes(b))
        if !haskey(c, k)
            c[k] = unify(v, Nothing)
        end
    end

    return NamedTuple{tuple(keys(c)...),Tuple{values(c)...}}
end

unify(a::Type{Vector{T}}, b::Type{Vector{S}}) where {T,S} = Vector{unify(T, S)}

# parse json into a type
function build_type(o::JSON3.Object)
    d = Dict()
    for (k, v) in o
        d[k] = build_type(v)
    end

    return NamedTuple{tuple(keys(d)...),Tuple{values(d)...}}
end

function build_type(a::JSON3.Array)
    t = Set([])
    nt = Top
    for item in a
        it = build_type(item)
        if it <: NamedTuple
            nt = unify(nt, it)
        else
            push!(t, build_type(item))
        end
    end

    return Vector{foldl(unify, t; init = Union{nt})}
end

build_type(x::T) where {T} = T


text:
module QuantizedStateSystems

using LinearAlgebra

include("qss1.jl")

export qss1

end

ext: jl

max_stars_repo_name: hurak/QuantizedStateSystems

max_issues_repo_name: hurak/QuantizedStateSystems

max_forks_repo_name: hurak/QuantizedStateSystems


text:
import Interpolations:
        CubicSplineInterpolation,
        Line, OnGrid, Throw

"""
Cubic spline interpolation of `x` and `y` on a regular grid `inpt_grid_x`. Returns the
interpolated `y` values.
"""
function cubicspline_interpolation(x::AbstractRange, y, inpt_grid_x;
        bc = Line(OnGrid()),
        extrapolation_bc = NaN,
        sort = true)

    # Knots must be sorted in ascending order.
    if sort
        sort_inds = sortperm(x)
        # Interpolation object
        itp = CubicSplineInterpolation(x[sort_inds], y[sort_inds],
                bc = bc,
                extrapolation_bc = extrapolation_bc)
    else
        itp = CubicSplineInterpolation(x, y,
                bc = bc,
                extrapolation_bc = extrapolation_bc)
    end

    # Evaluate the interpolation object at the x-axis grid we provide. These will
    # be our interpolated y-axis values
    itp(inpt_grid_x)
end

export cubicspline_interpolation


text:
"""
    function get_binding_lines(args...; kwargs...)

Finds the binding lines of a power system which has been solved by OPF, and reorders them in
accordance with their assigned index.

#Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_lines_tol::Float64 = 1e-1`:         Tolerance to identify any binding lines

"""
function get_binding_lines(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_lines_tol::Float64 = 1e-1
    )
    #Get all lines
    lines = get_components(Line, sys)

    #Gets lines with index
    all_lines, all_lines_sorted = assign_lines_index(lines)

    #Get components of the binding lines
    dual_lines = get_duals(res)[:network_flow]
    numberof_bind_lines = 0
    bind_lines_numbers = Array{Int64}(undef, 0)
    bind_lines_names = Array{String}(undef, 0)
    line_count = 0
    for line in lines
        line_count = line_count + 1
        line_name=line.name
        if !isapprox(dual_lines[1,line_name], 0.0; atol = dual_lines_tol)
            #if dual is different from 0.0 the line is binding
            numberof_bind_lines = numberof_bind_lines + 1
            resize!(bind_lines_numbers, numberof_bind_lines)
            resize!(bind_lines_names, numberof_bind_lines)
            bind_lines_numbers[numberof_bind_lines] = all_lines[line_count, 2]
            bind_lines_names[numberof_bind_lines] = line.name
        end
    end
    bind_lines = Array{Any}(undef,(numberof_bind_lines, 2))
    bind_lines[:, 1] = bind_lines_names
    bind_lines[:, 2] = bind_lines_numbers
    bind_lines = bind_lines[sortperm(bind_lines[:, 1]), :]
    return bind_lines
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# Implements node data types and associated helper functions
export
    SumNode,
    ProductNode,
    CategoricalDistribution,
    IndicatorFunction,
    GaussianDistribution

"""
Node Data Structures

Implement a labeled sparse matrix.
"""
abstract type Node end

" Sum node data type "
struct SumNode <: Node
    children::Vector{UInt}
    weights::Vector{Float64}
    # SumNode() = new(Vector{UInt}(),Vector{Float64}())
    # SumNode(children::Vector{<:Integer},weights::Vector{Float64}) = new(children,weights)
end
" Product node data type "
struct ProductNode <: Node
    children::Vector{UInt}
    #ProductNode() = new(Vector{UInt}())
    #ProductNode(children::Vector{<:Integer}) = new(children)
    #ProductNode(children) = new(children)
end
" Abstract leaf node type "
abstract type LeafNode <: Node end
"""
Indicator Function Node. Tolerance sets a maximum discrepancy when evaluating the node at a given value. Its default value is 1e-6.
"""
struct IndicatorFunction <: LeafNode
    scope::UInt
    value::Float64
    tolerance::Float64
    IndicatorFunction(scope::Integer,value::Float64) = new(scope,value,1e-6)
    IndicatorFunction(scope::Integer,value::Integer) = new(scope,Float64(value),1e-6)
end
"""
Univariate Categorical Distribution Node
"""
struct CategoricalDistribution <: LeafNode
    scope::UInt
    values::Vector{Float64}
end
"""
Univariate Gaussian Distribution Node
"""
mutable struct GaussianDistribution <: LeafNode
    scope::UInt
    mean::Float64
    variance::Float64
end

# LeafNode = Union{IndicatorFunction,CategoricalDistribution,GaussianDistribution}

"""
    IndicatorFunction(x::Vector{<:Real})::Float64

Evaluates indicator function at given configuration x.
"""
function (n::IndicatorFunction)(x::AbstractVector{<:Real})::Float64
    return isnan(x[n.scope]) ? 1.0 : n.value ≈ x[n.scope] ? 1.0 : 0.0
end
"""
Evaluates categorical distribution at given configuration
"""
function (n::CategoricalDistribution)(x::AbstractVector{<:Real})::Float64
    return isnan(x[n.scope]) ? 1.0 : n.values[Int(x[n.scope])]
end
"""
Evaluates Gaussian distribution at given configuration
"""
function (n::GaussianDistribution)(x::AbstractVector{<:Real})::Float64
    return isnan(x[n.scope]) ? 1.0 : exp(-(x[n.scope]-n.mean)^2/(2*n.variance))/sqrt(2*π*n.variance)
end

"Is this a leaf node?"
@inline isleaf(n::Node) = isa(n,LeafNode)
"Is this a sum node?"
@inline issum(n::Node) = isa(n,SumNode)
"Is this a product node?"
@inline isprod(n::Node) = isa(n,ProductNode)

"""
    logpdf(node,value)

Evaluates leaf `node` at the given `value` in log domain.
"""
@inline logpdf(n::IndicatorFunction, value::Integer) = isnan(value) ? 0.0 : value == Int(n.value) ? 0.0 : -Inf
@inline logpdf(n::IndicatorFunction, value::Float64) = isnan(value) ? 0.0 : abs(value - n.value) < n.tolerance  ? 0.0 : -Inf
@inline logpdf(n::CategoricalDistribution, value::Integer) = log(n.values[value])
@inline logpdf(n::CategoricalDistribution, value::Float64) = isnan(value) ? 0.0 : logpdf(n,Int(value))
@inline logpdf(n::GaussianDistribution, value::Float64)::Float64 = isnan(value) ? 0.0 : (-(value-n.mean)^2/(2*n.variance)) - log(2*π*n.variance)/2
"""
    maximum(node)

Returns the maximum value of the distribution
"""
@inline Base.maximum(n::IndicatorFunction) = 1.0
@inline Base.maximum(n::CategoricalDistribution) = Base.maximum(n.values)
@inline Base.maximum(n::GaussianDistribution) =  1/sqrt(2*π*n.variance)
"""
    argmax(node)

Returns the value at which the distribution is maximum
"""
@inline Base.argmax(n::IndicatorFunction) = n.value
@inline Base.argmax(n::CategoricalDistribution) = Base.argmax(n.values)
@inline Base.argmax(n::GaussianDistribution) =  n.mean

"""
    scope(node)

Returns the scope of a leaf node
"""
scope(n::LeafNode) = n.scope


text:
module Juriba
    include("Constants.jl")
    include("Particles.jl")
    include("Radiation.jl")
    include("radiation/BBody.jl")
    include("particles/Powerlaw.jl")
    include("particles/Thermal.jl")
    include("particles/Bknpower.jl")

end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
module Bi01KP

using DataStructures

include("problem.jl") #Defines problems and reduce_problem!()
include("mono_problem.jl") #Defines mono_problems
include("solution.jl") #Defines solutions
include("mono_solution.jl") #Defines mono_solutions
include("combo.jl") #API to call Combo to solve mono-objective problems
include("first_phase.jl") #first_phase(pb::problem) -> XSE1m
include("vertex.jl") #Defines vertices
include("arc.jl") #Defines arcs, paths(arclists) and partitions
include("second_phase.jl")
include("graph.jl")
# include("BoundSetReduction.jl")

const libcomboPath = joinpath(@__DIR__,"..","..","deps","libcombo.so")

function solve(pb::problem, output=true)
    reduce_problem!(pb, output)
    XSE1m = first_phase(pb)
    if length(XSE1m) == 1
        X = XSE1m
    else
        X = second_phase(XSE1m, output)
    end
    z1, z2 = map(obj_1, X), map(obj_2, X)
    w = map(weight, X)
    x = map(full_variable_vector, X)
    z1, z2, w, x
end

function parseKP(fname::AbstractString)
    if fname[end-2:end] == "KNP"
        f = open(fname)
        while !(readline(f)[1]=='n') end
        readline(f)
        x = readdlm(f)
        p1 = convert(Vector{Int},x[1:findfirst(x[:,3], "")-1 , 3])
        p2 = convert(Vector{Int},x[1:findfirst(x[:,4], "")-1 , 4])
        w = convert(Vector{Int},x[1:findfirst(x[:,2], "")-1 , 2])
        c = Int(x[findfirst(x[:,1], "W"), 2])
    elseif any(contains.(fname, ["TA","TB","TC","TD"]))
        f = open(fname)
        x = readdlm(f, skipstart=4)
        p1 = convert(Vector{Int},x[1:findfirst(x[:,3], "")-1 , 3])
        p2 = convert(Vector{Int},x[1:findfirst(x[:,4], "")-1 , 4])
        w = convert(Vector{Int},x[1:findfirst(x[:,2], "")-1 , 2])
        c = Int(x[findfirst(x[:,1], "W"), 2])
    else
        f = open(fname)
        x = convert(Array{Int},readdlm(f, skipstart=11))[:,1]
        close(f)
        n = (length(x)-1) ÷ 3
        c = x[end]
        p1 = x[1:n]
        p2 = x[n+1:2n]
        w = x[2n+1:3n]
    end
    return problem(p1,p2,w,c)
end

end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
module WriteVTK

# Documentation on VTK XML formats:
# - https://vtk.org/Wiki/VTK_XML_Formats
# - http://vtk.org/VTK/img/file-formats.pdf

export MeshCell, VTKPolyhedron
export vtk_grid, vtk_save, vtk_point_data, vtk_cell_data
export vtk_multiblock, multiblock_add_block
export paraview_collection, collection_add_timestep, paraview_collection_load
export vtk_write_array
export VTKPointData, VTKCellData, VTKFieldData
export PolyData
export pvtk_grid

import CodecZlib: ZlibCompressorStream
import TranscodingStreams

using LightXML
using FillArrays: Zeros

using Base64: base64encode

import Base: close, isopen, show

# Cell type definitions as in vtkCellType.h
include("VTKCellTypes.jl")
using .VTKCellTypes
export VTKCellTypes, VTKCellType

## Constants ##
const DEFAULT_COMPRESSION_LEVEL = 6
const IS_LITTLE_ENDIAN = ENDIAN_BOM == 0x04030201
const HeaderType = UInt64  # should be UInt32 or UInt64

## Types ##
abstract type VTKFile end

_compression_level(x::Bool) = x ? DEFAULT_COMPRESSION_LEVEL : 0
_compression_level(x) = Int(x)

struct DatasetFile <: VTKFile
    xdoc::XMLDocument
    path::String
    grid_type::String
    Npts::Int           # Number of grid points.
    Ncls::Int           # Number of cells.
    compression_level::Int  # Compression level for zlib (if 0, compression is disabled)
    appended::Bool      # Data is appended? (otherwise it's written inline, base64-encoded)
    ascii::Bool         # if true, inline data is written in ASCII format (only used if `appended` is false)
    buf::IOBuffer       # Buffer with appended data.
    function DatasetFile(xdoc, path, grid_type, Npts, Ncls;
                         compress=true, append=true, ascii=false)
        clevel = _compression_level(compress)
        if !(0 ≤ clevel ≤ 9)
            throw(ArgumentError("unexpected value of `compress` argument: $compress.\n" *
                                "It must be a `Bool` or a value between 0 and 9."))
        end
        ascii = !append && ascii
        if ascii
            clevel = 0  # no compression when writing ASCII
        end
        buf = IOBuffer()
        if !append  # in this case we don't need a buffer
            close(buf)
        end
        new(xdoc, path, grid_type, Npts, Ncls, clevel, append, ascii, buf)
    end
end

DatasetFile(dtype, xdoc::XMLDocument, fname::AbstractString, args...; kwargs...) =
    DatasetFile(xdoc, add_extension(fname, dtype), xml_name(dtype), args...; kwargs...)

function show(io::IO, vtk::DatasetFile)
    open_str = isopen(vtk) ? "open" : "closed"
    print(io, "VTK file '$(vtk.path)' ($(vtk.grid_type) file, $open_str)")
end

struct VTKBlock
    xelm::XMLElement
    blocks::Vector{Union{VTKFile,VTKBlock}}
    # Constructor.
    VTKBlock(xelm) = new(xelm, Union{VTKFile,VTKBlock}[])
end

struct MultiblockFile <: VTKFile
    xdoc::XMLDocument
    path::String
    blocks::Vector{Union{VTKFile,VTKBlock}}
    MultiblockFile(xdoc, path) = new(xdoc, path, Union{VTKFile,VTKBlock}[])
end

struct CollectionFile <: VTKFile
    xdoc::XMLDocument
    path::String
    timeSteps::Vector{String}
    CollectionFile(xdoc, path) = new(xdoc, path, VTKFile[])
end

"""
    close(vtk::VTKFile)

Write and close VTK file.
"""
close(vtk::VTKFile) = free(vtk.xdoc)

"""
    isopen(vtk::VTKFile)

Check if VTK file is still being written.
"""
isopen(vtk::VTKFile) = (vtk.xdoc.ptr != C_NULL)

# Add a default extension to the filename, unless the user have already given
# the correct one.
function add_extension(filename, default_extension::AbstractString) :: String
    path, ext = splitext(filename)
    if ext == default_extension
        return filename
    end
    if !isempty(ext)
        @warn("detected extension '$(ext)' does not correspond to " *
              "dataset type.\nAppending '$(default_extension)' to filename.")
    end
    filename * default_extension
end

add_extension(filename, dtype) = add_extension(filename, file_extension(dtype))

# Common functions and types.
include("write_data.jl")
include("save_files.jl")
include("gridtypes/types.jl")
include("gridtypes/common.jl")

# Multiblock-specific functions and types.
include("gridtypes/multiblock.jl")
include("gridtypes/ParaviewCollection.jl")

# Grid-specific functions and types.
include("gridtypes/structured/imagedata.jl")
include("gridtypes/structured/rectilinear.jl")
include("gridtypes/structured/structured.jl")

include("gridtypes/unstructured/unstructured.jl")
include("gridtypes/unstructured/polyhedron.jl")
include("gridtypes/unstructured/polydata.jl")

include("gridtypes/array.jl")

# Parallel DataSet Formats
include("gridtypes/pvtk_grid.jl")

# This allows using do-block syntax for generation of VTK files.
for func in (:vtk_grid, :vtk_multiblock, :paraview_collection,
             :paraview_collection_load)
    @eval begin
        """
            $($func)(f::Function, args...; kwargs...)

        Create VTK file and apply `f` to it.
        The file is automatically closed by the end of the call.

        This allows to use the do-block syntax for creating VTK files:

        ```julia
        saved_files = $($func)(args...; kwargs...) do vtk
            # do stuff with the `vtk` file handler
        end
        ```
        """
        function ($func)(f::Function, args...; kwargs...)
            vtk = ($func)(args...; kwargs...)
            local outfiles
            try
                f(vtk)
            finally
                outfiles = vtk_save(vtk)
            end
            outfiles :: Vector{String}
        end
    end
end

end


text:
#=
Takes one step in the Markov chain according to the given state's transition
vector.
Inputs:
- transition_vector: An array of length p where the ith element denotes the
                     probability of transitioning into state i from the current
                     state (which is implicitly given).
Outputs:
- An integer in {1,2,..., p} denoting the index of the next state which occurs
with the probability specified in transition_vector.
=#
function one_step(transition_vector)
    cumulative_vector = cumsum(transition_vector)
    state_index = 1
    U = rand()
    while(U > cumulative_vector[state_index])
        state_index += 1
    end

    return state_index
end

#=
Runs a Markov chain a fixed number of times and returns the final state.
Inputs:
- states: An array of length p where each element represents a state of the
          chain. Assumes the ordering is such that the probability of
          transitioning from state i to state j is given by P[i,j].
          Examples:
          [0,1,2,3,4,5]
          ["red", "green", "blue"]
          [0,"blue",3,2]
- P: A (p x p) matrix denoting the transition matrix of the Markov chain.
     Assumes the rows of P sum to 1 and P[i,j] >= 0 for all i,j.
- S0: The initial state of the Markov chain. Assumed to be an element of the
      state variable.
- steps: A nonnegative integer specifying the number of steps for the Markov
         chain to go through.
Outputs:
- An element of the states variable denoting the state the chain is in after
  the specified number of steps.
=#
function run_fixed_simulation(states, P, S0, steps)

    # Define a dictionary mapping state to their index in P
    state_dict = Dict(states[i] => i for i = 1:length(states))

    # Define the current state of the chain
    Sn = S0

    # Run through the specified number of steps
    for i = 1:steps
        state_index = get(state_dict, Sn, -1)
        next_state_index = one_step(P[state_index,:])
        Sn = states[next_state_index]
    end

    return Sn
end

#=
Runs a Markov chain until a state is reached.
Inputs:
- states: Same as defined in ``run_fixed_simulation`` function
- P: Same as defined in ``run_fixed_simulation`` function
- S0: Same as defined in ``run_fixed_simulation`` function
- Sf: An array of states defining a possible final state. Assumed to be a
      subset of the states variable.
Outputs:
- An integer representing the number of steps it took to get to get to state
  contained in the Sf array.
=#
function run_simulation(states, P, S0, Sf)

    # Define a dictionary mapping state to their index in P
    state_dict = Dict(states[i] => i for i = 1:length(states))

    # Define the current state of the chain
    Sn = S0

    # Keep tracks of number of steps taken
    steps = 0

    # Run through the Markov changing
    while !(Sn in Sf) && steps < 1000
        state_index = get(state_dict, Sn, -1)
        next_state_index = one_step(P[state_index,:])
        Sn = states[next_state_index]
        steps += 1
    end

    return steps
end

# Define gambler's ruin Markov chain with N = 5 and p = .3

# gambler's ruin state space (χ)
states = [0,1,2,3,4,5]

# Transition matrix for gambler's ruin
P = [1 0 0 0 0 0;
    .7 0 .3 0 0 0;
    0 .7 0 .3 0 0;
    0 0 .7 0 .3 0;
    0 0 0 .7 0 .3;
    0 0 0 0 0 1]

# Estimate probability of having $5 after 6 turns when starting with $2
num_simulations = 100000
mc_est = mean([run_fixed_simulation(states, P, 2, 6) == 5
 for i = 1:num_simulations])

# Define masters program Markov chain
states = [1,2,"G","D"]

# Transition matrix for masters Markov chain
P = [.4 .5 0 .1;
     0 .3 .6 .1;
     0 0 1 0;
     0 0 0 1]

mc_est = mean([run_simulation(states, P, 1, ["G", "D"]) for i = 1:10000])

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
module MacroModel
using JuMP, GLPK, DelimitedFiles, LinearAlgebra, DataFrames, CSV, Logging, Printf, Suppressor

export runleapmacromodel

include("./IOlib.jl")
include("./LEAPfunctions.jl")
using .IOlib, .LEAPfunctions

"""
    output_var(params, values, filename, index, rowlabel, mode)

Report output values by writing to specified CSV files.
"""
function output_var(params, values, filename, index, rowlabel, mode)
	if isa(values, Array)
		usevalue = join(values, ',')
	else
		usevalue = values
	end
	open(joinpath(params["results_path"], string(filename, "_", index,".csv")), mode) do io
		write(io, string(rowlabel, ',', usevalue, "\r\n"))
	end
end

"""
    to_quoted_string_vec(svec)

Convert a vector of strings by wrapping each string in quotes (for putting into CSV files).
This function converts `svec` in-place, so should only be run once. Double quotation marks
are removed to avoid that problem.
"""
function to_quoted_string_vec(svec)
	for i in 1:length(svec)
		svec[i] = replace(string('\"', svec[i], '\"'), "\"\"" => "\"")
	end
end

"""
    calc_sraffa_inverse(np::Int64, ns::Int64, tradeables::Array{Bool,1}, io::IOdata)

Calculate the Sraffa inverse matrix, which is used to calculate domestic prices
"""
function calc_sraffa_inverse(np::Int64, ns::Int64, tradeables::Array{Bool,1}, io::IOdata)
	sraffa_matrix = Array{Float64}(undef, np, np)
	sraffa_matrix .= 0.0
	# Domestic prices
	for i in findall(.!tradeables)
		for j in 1:np
			sraffa_matrix[i,j] = ((1 + io.τd[j])/(1 + io.τd[i])) * sum(io.μ[r] * io.S[r,i] * io.D[j,r] for r in 1:ns)
		end
	end
    inv(LinearAlgebra.I - sraffa_matrix)
end

"""
    calc_sraffa_RHS(t::Int64, np::Int64, ns::Int64, ω::Array{Float64,1}, tradeables::Array{Bool,1}, prices::PriceData, io::IOdata, exog::ExogParams)

Calculate the Sraffa right-hand side expression, which is multiplied by the Sraffa inverse matrix to give domestic prices
"""
function calc_sraffa_RHS(t::Int64, np::Int64, ns::Int64, ω::Array{Float64,1}, tradeables::Array{Bool,1}, prices::PriceData, io::IOdata, exog::ExogParams)
	sraffa_wage_vector = Array{Float64}(undef, np)
	sraffa_wage_vector .= 0.0 # Initialize to zero
	for i in findall(.!tradeables)
		sraffa_wage_vector[i] = (1/(1 + io.τd[i])) * sum(io.μ[r] * io.S[r,i] * ω[r] for r in 1:ns)
	end
	# Initialize prices.pd
	sraffa_RHS = prices.Pg * sraffa_wage_vector
	for i in findall(tradeables)
		sraffa_RHS[i] = exog.xr[t] * prices.pw[i]
	end
	return sraffa_RHS
end

"""
    ModelCalculations(file::String, I_en::Array, run::Int64)

Implement the Macro model. This is the main function.
"""
function ModelCalculations(file::String, I_en::Array, run::Int64)

    #------------status
    @info "Loading data..."
    #------------status

    params = IOlib.parse_input_file(file)

	# Clean up folders if requested
	if run == 0
		if params["clear-folders"]["results"] & isdir(params["results_path"])
			for f in readdir(params["results_path"])
				rm(joinpath(params["results_path"], f))
			end
		end
		if params["clear-folders"]["calibration"] & isdir(params["calibration_path"])
			for f in readdir(params["calibration_path"])
				rm(joinpath(params["calibration_path"], f))
			end
		end
		if params["clear-folders"]["diagnostics"] & isdir(params["diagnostics_path"])
			for f in readdir(params["diagnostics_path"])
				rm(joinpath(params["diagnostics_path"], f))
			end
			if !params["report-diagnostics"]
				rm(params["diagnostics_path"])
			end
		end

		# Ensure that they exist
		mkpath(params["results_path"])
		mkpath(params["calibration_path"])
		if params["report-diagnostics"]
			mkpath(params["diagnostics_path"])
		end

	end

	io, np, ns = IOlib.supplyusedata(file)
	prices = IOlib.prices_init(np, io)

	#############################################################################
	#
    # Assign selected configuration file parameters to variables
	#
	#############################################################################

	# Years
    base_year = params["years"]["start"]
    final_year = params["years"]["end"]
	# Investment function
    neutral_growth = params["investment-fcn"]["init_neutral_growth"]
    util_sens_scalar = params["investment-fcn"]["util_sens"]
    profit_sens = params["investment-fcn"]["profit_sens"]
    intrate_sens = params["investment-fcn"]["intrate_sens"]
    growth_adj = params["investment-fcn"]["growth_adj"]
	# Linear program objective function
    wu = params["objective-fcn"]["category_weights"]["utilization"]
    wf = params["objective-fcn"]["category_weights"]["final_demand_cov"]
    wx = params["objective-fcn"]["category_weights"]["exports_cov"]
    wm = params["objective-fcn"]["category_weights"]["imports_cov"]
	ϕu = params["objective-fcn"]["product_sector_weight_factors"]["utilization"]
	ϕf = params["objective-fcn"]["product_sector_weight_factors"]["final_demand_cov"]
	ϕx = params["objective-fcn"]["product_sector_weight_factors"]["exports_cov"]
	# Taylor function
	tf_γ0 = neutral_growth
	tf_min_γ0 = params["taylor-fcn"]["neutral_growth_band"][1]
	tf_max_γ0 = params["taylor-fcn"]["neutral_growth_band"][2]
    tf_i_targ = params["taylor-fcn"]["target_intrate"]
    tf_gr_resp = params["taylor-fcn"]["gr_resp"]
    tf_π_targ = params["taylor-fcn"]["target_infl"]
	tf_π_targ_use_πw = isnothing(tf_π_targ)
    tf_infl_resp = params["taylor-fcn"]["infl_resp"]
	# Wage rate function
	infl_passthrough = params["wage-fcn"]["infl_passthrough"]
	lab_constr_coeff = params["wage-fcn"]["lab_constr_coeff"]
	LEAP_indices = params["LEAP_sector_indices"]

	#----------------------------------
    # Calculate variables based on parameters
    #----------------------------------
	# Number of time steps
	ntime = 1 + (final_year - base_year)
	# Initial autonomous growth rate
	γ_0 = neutral_growth * ones(ns)
	# Set up a "tradeables" filter
    tradeables = [true for i in 1:np] # Initialize to true
    tradeables[params["non-tradeable-range"]] .= false

	#----------------------------------
    # Read in exogenous time-varying parameters and sector-specific parameters
    #----------------------------------
	exog = IOlib.get_var_params(file)
	# Convert pw into global currency
	prices.pw = prices.pw/exog.xr[1]

	#############################################################################
	#
    # Initialize variables (some of them adjusted after a later calibration step)
	#
	#############################################################################

	#----------------------------------
    # For calibration, apply a user-specified adjustment to exports, final demands, and potential output
    #----------------------------------
	Xmax = (1.0 + params["calib"]["max_export_adj_factor"]) * io.X
    Fmax = (1.0 + params["calib"]["max_hh_dmd_adj_factor"]) * max.(io.F,zeros(np))
	z = (1.0 + params["calib"]["pot_output_adj_factor"]) * io.g

	#----------------------------------
    # Investment
    #----------------------------------
	# Total investment
    I_total = sum(io.I)
	# Track non-energy investment separate from energy investment (which is supplied from LEAP)
    I_ne = I_total - I_en[1]
    # Share of supply of investment goods
    θ = io.I / sum(io.I)
	# Investment function parameter -- make a vector
    util_sens = util_sens_scalar * ones(ns)
	# Initial values: These will change endogenously over time
	γ = γ_0
	i_bank = tf_i_targ

	#----------------------------------
    # Wages
    #----------------------------------
    W = io.W
	ω = io.W ./ io.g
    replace!(ω, NaN=>0)
	wage_ratio = 1.0

	#----------------------------------
    # Domestic prices using a Sraffian price system
    #----------------------------------
	sraffa_inverse = calc_sraffa_inverse(np, ns, tradeables, io)
    prices.pd = sraffa_inverse * calc_sraffa_RHS(1, np, ns, ω, tradeables, prices, io, exog)
	prices.pb = prices.pd

	#----------------------------------
    # Capital productivity of new investment
    #----------------------------------
	# Intermediate variables
	profit_per_output = io.Vnorm * prices.pd - (prices.Pg .* ω +  transpose(io.D) * ((ones(np) + io.τd) .* prices.pd))
	price_of_capital = dot(θ,(ones(np) + io.τd) .* prices.pd)
	I_nextper = (1 + neutral_growth) * (1 + params["calib"]["nextper_inv_adj_factor"]) * I_ne
	profit_share_rel_capprice = (1/price_of_capital) * profit_per_output
	init_profit = sum((γ_0 + exog.δ) .* io.g .* profit_share_rel_capprice)
    init_profit_rate = init_profit/I_nextper
	# Capital productivity
    capital_output_ratio = profit_share_rel_capprice / init_profit_rate
	replace!(capital_output_ratio, NaN=>0) # replace NaN with 0
	# Initialize for the investment function
	targ_profit_rate = init_profit_rate

	#############################################################################
	#
    # Define linear goal program
	#
	#############################################################################

    #------------status
    @info "Preparing model..."
    #------------status

    # model
    mdl = Model(GLPK.Optimizer)

	#----------------------------------
	# Product and sector weights for the objective function
	#----------------------------------
	u_sector_wt = ϕu * io.g / sum(io.g) + (1.0 - ϕu) * ones(ns)/ns
	f_sector_wt = ϕf * io.F / sum(io.F) + (1.0 - ϕf) * ones(np)/np
	x_sector_wt = ϕx * io.X / sum(io.X) + (1.0 - ϕx) * ones(np)/np

	#----------------------------------
	# A filter for products that are domestically supplied
	#----------------------------------
	no_dom_production = vec(sum(io.S, dims=1)) .== 0

	#----------------------------------
    # Dynamic parameters
    #----------------------------------
	param_z = z
	param_Xmax = Xmax
	param_Fmax = Fmax
	param_I_tot = I_total
	param_prev_I_tot = 2 * I_total # Allow extra slack -- this sets a scale
	param_pw = prices.pw
	param_pd = prices.pd
	param_Pg = prices.Pg
	param_pb = prices.pb
	param_Mref = 2 * io.M # Allow for some extra slack -- this just sets a scale
	param_mfrac = io.m_frac
    #----------------------------------
    # Variables
    #----------------------------------
    # For objective function
    @variable(mdl, 0 <= ugap[i in 1:ns] <= 1)
    @variable(mdl, 0 <= fgap[i in 1:np] <= 1)
    @variable(mdl, 0 <= xgap[i in 1:np] <= 1)
    @variable(mdl, 0 <= ψ_pos[i in 1:np] <= 1) # Excess over target imports (as a relative amount)
    @variable(mdl, 0 <= ψ_neg[i in 1:np] <= 1) # Deficit below target imports (as a relative amount)
    # Supply and demand quantities
    @variable(mdl, qs[i in 1:np] >= 0)
    @variable(mdl, qd[i in 1:np] >= 0)
    # Final demand
    @variable(mdl, F[i in 1:np] >= 0)
    @variable(mdl, 0 <= fshare[i in 1:np] <= 1)
    @variable(mdl, M[i in 1:np] >= 0)
    @variable(mdl, 0 <= M_share[i in 1:np] <= 1)
    @variable(mdl, X[i in 1:np] >= 0)
    @variable(mdl, 0 <= xshare[i in 1:np] <= 1)
	@variable(mdl, I_tot) # Fix to equal param_I_tot
	@variable(mdl, I_supply[i in 1:np]) # Supply of investment goods by product
    # Output
    @variable(mdl, 0 <= u[i in 1:ns] <= 1)
    # Margins
    @variable(mdl, margins_pos[i in 1:np] >= 0)
    @variable(mdl, margins_neg[i in 1:np] >= 0)

    #----------------------------------
    # Objective function (goal program)
    #----------------------------------
    @objective(mdl, Min, wu * sum(u_sector_wt[i] * ugap[i] for i in 1:ns) +
					     wf * sum(f_sector_wt[i] * fgap[i] for i in 1:np) +
						 wx * sum(x_sector_wt[i] * xgap[i] for i in 1:np) +
						 wm * (sum(ψ_pos[i] + ψ_neg[i] for i in 1:np)))

    #----------------------------------
    # Constraints
    #----------------------------------
    # Sector constraints
	@constraint(mdl, eq_util[i = 1:ns], u[i] + ugap[i] == 1.0)
	# Fundamental input-output equation
	@constraint(mdl, eq_io[i = 1:ns], sum(io.S[i,j] * param_pb[j] * qs[j] for j in 1:np) - param_Pg * param_z[i] * u[i] == 0)
    # Product constraints
	# Variables for objective function (goal program)
	@constraint(mdl, eq_xshare[i = 1:np], xshare[i] + xgap[i] == 1.0)
	@constraint(mdl, eq_X[i = 1:np], X[i] == xshare[i] * param_Xmax[i])
	@constraint(mdl, eq_fshare[i = 1:np], fshare[i] + fgap[i] == 1.0)
	@constraint(mdl, eq_F[i = 1:np], F[i] == fshare[i] * param_Fmax[i])
	# Intermediate demand
	@constraint(mdl, eq_intdmd[i = 1:np], qd[i] - sum(io.D[i,j] * u[j] * param_z[j] for j in 1:ns) == 0)
	# Investment
	fix(I_tot, param_I_tot)
	@constraint(mdl, eq_inv_supply[i = 1:np], I_supply[i] == θ[i] * I_tot)
	# Total supply
	@constraint(mdl, eq_totsupply[i = 1:np], qs[i] == qd[i] - margins_pos[i] + margins_neg[i] + X[i] + F[i] + I_supply[i] - M[i])
	@constraint(mdl, eq_no_dom_prod[i = 1:np], qs[i] * no_dom_production[i] == 0)
	# Imports
	@constraint(mdl, eq_M[i = 1:np], M[i] == param_mfrac[i] * (qd[i] + F[i] + I_supply[i]) + param_Mref[i] * (ψ_pos[i] - ψ_neg[i]))
	# Margins
	@constraint(mdl, eq_margpos[i = 1:np], margins_pos[i] == io.marg_pos_ratio[i] * (qs[i] + M[i]))
	@constraint(mdl, eq_margneg[i = 1:np], margins_neg[i] == io.marg_neg_share[i] * sum(margins_pos[j] for j in 1:np))

    #------------------------------------------
    # Calibration run
    #------------------------------------------
	if params["report-diagnostics"]
		open(joinpath(params["diagnostics_path"], string("model_", run, "_calib_1", ".txt")), "w") do f
			print(f, mdl)
		end
	end
    optim_output = @capture_out begin
		optimize!(mdl)
	end
	if !isempty(optim_output)
		@info optim_output
	end
    status = primal_status(mdl)
    @info "Calibrating: $status"

	# Sector variables
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("capacity_utilization_",run,".csv")), value.(u), "capacity utilization", params["included_sector_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("sector_output_",run,".csv")), value.(u) .* z, "sector output", params["included_sector_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("wage_share_",run,".csv")), ω, "wage share", params["included_sector_codes"])
	IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("capital_output_ratio_",run,".csv")), capital_output_ratio, "capital-output ratio", params["included_sector_codes"])

	# Product variables
	IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("exports_",run,".csv")), value.(X), "exports", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("final_demand_",run,".csv")), value.(F), "final demand", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("imports_",run,".csv")), value.(M), "imports", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("domestic_production_",run,".csv")), value.(qs), "domestic production", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("tot_intermediate_supply_non-energy_sectors_",run,".csv")), value.(qd), "intermediate supply from non-energy sectors", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("basic_prices_",run,".csv")), param_pb, "basic prices", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("margins_neg_",run,".csv")), value.(margins_neg), "negative margins", params["included_product_codes"])
    IOlib.write_vector_to_csv(joinpath(params["calibration_path"], string("margins_pos_",run,".csv")), value.(margins_pos), "positive margins", params["included_product_codes"])

    # First run is calibration -- now set Xmax and Fmax based on solution and run again
    Xmax = max.(Xmax, value.(X))
    Fmax = max.(Fmax, value.(F))
	param_Xmax = Xmax
	param_Fmax = Fmax
	for i in 1:np
		set_normalized_coefficient(eq_X[i], xshare[i], -param_Xmax[i])
		set_normalized_coefficient(eq_F[i], fshare[i], -param_Fmax[i])
	end
	if params["report-diagnostics"]
		open(joinpath(params["diagnostics_path"], string("model_", run, "_calib_2", ".txt")), "w") do f
			print(f, mdl)
		end
	end
	optim_output = @capture_out begin
		optimize!(mdl)
	end
	if !isempty(optim_output)
		@info optim_output
	end
    status = primal_status(mdl)

	#--------------------------------
	# Initialize variables for calculating inflation, growth rates, etc.
	#--------------------------------
    pd_prev = prices.pd/(1 + params["global-params"]["infl_default"])
    pb_prev = param_pb/(1 + params["global-params"]["infl_default"])
    prev_GDP = sum(param_pb[i] * (value.(qs) - value.(qd))[i] for i in 1:np)/(1 + neutral_growth)

	prev_GDP_deflator = 1
    prev_GDP_gr = neutral_growth
    prev_πg = params["global-params"]["infl_default"]

	lab_force_index = 1

    year = base_year

	#--------------------------------
	# Initialize array of indices to pass to LEAP
	#--------------------------------
	labels = ["Year"]
	if haskey(params, "GDP-branch")
		labels = vcat(labels, params["GDP-branch"]["name"])
	end
	if haskey(params, "Employment-branch")
		labels = vcat(labels, params["Employment-branch"]["name"])
	end
	labels = vcat(labels, params["LEAP_sector_names"])
    indices = Array{Float64}(undef, ntime, length(labels))

	#------------------------------------------
    # Initialize files for writing results
    #------------------------------------------
	# Get sector and product labels
	sector_names = params["included_sector_names"]
	product_names = params["included_product_names"]
	# Quote the sector and product names (for putting into CSV files)
	to_quoted_string_vec(sector_names)
	to_quoted_string_vec(product_names)

	# Create files for sector variables
	output_var(params, sector_names, "sector_output", run, "", "w")
	output_var(params, sector_names, "potential_sector_output", run, "", "w")
	output_var(params, sector_names, "capacity_utilization", run, "", "w")
	output_var(params, sector_names, "real_value_added", run, "", "w")
	output_var(params, sector_names, "profit_rate", run, "", "w")
	output_var(params, sector_names, "autonomous_investment_rate", run, "", "w")
	# Create files for product variables
	output_var(params, product_names, "final_demand", run, "", "w")
	output_var(params, product_names, "imports", run, "", "w")
	output_var(params, product_names, "exports", run, "", "w")
	output_var(params, product_names, "basic_prices", run, "", "w")
	# Create a file to hold scalar variables
	scalar_var_list = ["GDP gr", "curr acct surplus to GDP ratio", "curr acct surplus", "real GDP",
					   "GDP deflator", "labor productivity gr", "labor force gr", "wage rate gr",
					   "wage per effective worker gr", "real investment", "central bank rate"]
	to_quoted_string_vec(scalar_var_list)
	output_var(params, scalar_var_list, "collected_variables", run, "", "w")

	#############################################################################
	#
    # Loop over years
	#
	#############################################################################
    @info "Running from $base_year to $final_year:"

    previous_failed = false
    for t in 1:ntime
		#--------------------------------
		# Output
		#--------------------------------
        g = value.(u) .* z
		g_share = pb_prev .* value.(qs)/sum(pb_prev .* value.(qs))

		#--------------------------------
		# Price indices and deflators
		#--------------------------------
        # This is used here if the previous calculation failed, but is also reported below
        va_at_prev_prices_1 = prices.Pg * g
        va_at_prev_prices_2 = sum(pd_prev[j] * io.D[j,:] for j in 1:np) .* g
        value_added_at_prev_prices = va_at_prev_prices_1 - va_at_prev_prices_2
        if previous_failed
            πg = prev_πg
            pd_prev = pd_prev * (1 + πg)
            pb_prev = pb_prev * (1 + πg)
            va1 = (1 + πg) * prices.Pg * g
            va2 = sum(value.(param_pd)[j] * io.D[j,:] for j in 1:np) .* g
            value_added = va1 - va2
            πGDP = sum(value_added)/sum(value_added_at_prev_prices) - 1
        else
			πd = (param_pd - pd_prev) ./ pd_prev
	        πb = (param_pb - pb_prev) ./ pb_prev
            πg = sum(g_share .* πb)
            prev_πg = πg
			val_findmd = (ones(np) + io.τd) .* pd_prev .* (value.(X) + value.(F) + value.(I_supply) - value.(M))
	        val_findmd_share = val_findmd/sum(val_findmd)
            πGDP = sum(val_findmd_share .* πd)
            pd_prev = param_pd
            pb_prev = param_pb
        end
        GDP_deflator = prev_GDP_deflator
        prev_GDP_deflator = (1 + πGDP) * prev_GDP_deflator

		if tf_π_targ_use_πw
			tf_π_targ = exog.πw[t]
		end

		#--------------------------------
		# GDP
		#--------------------------------
        GDP = sum(param_pb[i] * (value.(qs) - value.(qd))[i] for i in 1:np)/GDP_deflator
        if previous_failed
            GDP_gr = prev_GDP_gr
            prev_GDP = prev_GDP * (1 + GDP_gr)
        else
            GDP_gr = GDP/prev_GDP - 1.0
            prev_GDP = GDP
            prev_GDP_gr = GDP_gr
        end

		#--------------------------------
		# Wages and the labor force
		#--------------------------------
		λ_gr = exog.αKV[t] * GDP_gr + exog.βKV[t] # Kaldor-Verdoorn equation for labor productivity
		L_gr = GDP_gr - λ_gr # Labor force growth rate

		lab_force_index *= 1 + L_gr

		w_gr = infl_passthrough * πGDP + λ_gr * (1.0 + lab_constr_coeff * (L_gr - exog.working_age_grs[t]))
		ω_gr = w_gr - λ_gr - πg
		ω = (1.0 + ω_gr) * ω

		#--------------------------------
		# Profits
		#--------------------------------
		# First, update the Vnorm matrix
		io.Vnorm = Diagonal(1 ./ g) * io.S * Diagonal(value.(qs))
		profit_per_output = io.Vnorm * prices.pd - (prices.Pg .* ω +  transpose(io.D) * ((ones(np) + io.τd) .* prices.pd))
		pK = dot(θ,(ones(np) + io.τd) .* prices.pd)
		profit_rate = profit_per_output ./ (pK * capital_output_ratio)

		#--------------------------------
		# Investment
		#--------------------------------
		# Investment function
		γ_u = util_sens .* (value.(u) .- 1)
		γ_r = profit_sens * (profit_rate .- targ_profit_rate)
		γ_i = -intrate_sens * (i_bank - tf_i_targ) * ones(ns)
        γ = max.(γ_0 + γ_u + γ_r + γ_i, -exog.δ)
		# Non-energy investment, by sector and total
        I_ne_disag = z .* (γ + exog.δ) .* capital_output_ratio
        I_ne = sum(I_ne_disag)
		# Combine with energy investment
        I_total = I_ne + I_en[t]
		# Update autonomous investment term
        γ_0 = γ_0 + growth_adj * (γ - γ_0)

		#--------------------------------
		# Export demand
		#--------------------------------
        Xmax = Xmax .* (1 + exog.world_grs[t]).^exog.export_elast_demand[t]

		#--------------------------------
		# Final domestic consumption demand
		#--------------------------------
		# Wage ratio drives increase in quantity, so use real wage increase, deflate by Pg
		W_curr = sum(W)
        W = ((1 + w_gr)/(1 + λ_gr)) * W .* (1 .+ γ)
        wage_ratio = (1/(1 + πGDP)) * sum(W)/W_curr
        Fmax = Fmax .* wage_ratio.^exog.wage_elast_demand[t]

		#--------------------------------
		# Calculate indices to pass to LEAP
		#--------------------------------
        indices[t,1] = year
		curr_index = 1
		if haskey(params, "GDP-branch")
			curr_index += 1
			indices[t,curr_index] = GDP
		end
		if haskey(params, "Employment-branch")
			curr_index += 1
			indices[t,curr_index] = lab_force_index
		end

 		for i in 1:length(LEAP_indices)
			indices[t, curr_index + i] = 0
			for j in LEAP_indices[i]
				indices[t, curr_index + i] += g[j]
			end
		end

		#--------------------------------
		# Reporting
		#--------------------------------
		# These variables are only reported, not used
		CA_surplus = sum(prices.pw .* (value.(X) - value.(M)))/GDP_deflator
		CA_to_GDP_ratio = CA_surplus/GDP

		# Sector variables
		output_var(params, g, "sector_output", run, year, "a")
		output_var(params, z, "potential_sector_output", run, year, "a")
		output_var(params, value.(u), "capacity_utilization", run, year, "a")
		output_var(params, value_added_at_prev_prices/prev_GDP_deflator, "real_value_added", run, year, "a")
		output_var(params, profit_rate, "profit_rate", run, year, "a")
		output_var(params, γ_0, "autonomous_investment_rate", run, year, "a")
	    # Product variables
		output_var(params, value.(F), "final_demand", run, year, "a")
		output_var(params, value.(M), "imports", run, year, "a")
		output_var(params, value.(X), "exports", run, year, "a")
		output_var(params, param_pb, "basic_prices", run, year, "a")
		# Scalar variables
		scalar_var_vals = [GDP_gr, CA_to_GDP_ratio, CA_surplus, GDP, GDP_deflator, λ_gr, L_gr,
							w_gr, ω_gr, value.(I_tot), i_bank]
		output_var(params, scalar_var_vals, "collected_variables", run, year, "a")

		#--------------------------------
		# Update Taylor rule
		#--------------------------------
        i_bank = tf_i_targ + tf_gr_resp * (GDP_gr - tf_γ0) + tf_infl_resp * (πGDP - tf_π_targ)
		# Apply adaptive expectations, then bound
		tf_γ0 = min(max(tf_γ0 + growth_adj * (GDP_gr - tf_γ0), tf_min_γ0), tf_max_γ0)

		#--------------------------------
		# Update potential output
		#--------------------------------
        z = (1 .+ γ) .* z

		#--------------------------------
		# Update prices
		#--------------------------------
		prices.Pg = (1 + πg) * prices.Pg
		prices.pw = prices.pw * (1 + exog.πw[t])
		prices.pd = sraffa_inverse * calc_sraffa_RHS(t, np, ns, ω, tradeables, prices, io, exog)
		prices.pb = prices.pd

		#--------------------------------
		# Update the linear goal program
		#--------------------------------
		param_Mref = 2 * value.(M) # Allow for some extra slack -- this just sets a scale
		param_Pg = prices.Pg
		param_Xmax = Xmax
		param_Fmax = Fmax
		param_I_tot = I_total
		param_prev_I_tot = 2 * I_total # Allow extra slack -- this sets a scale
		param_pw = prices.pw
		param_pd = prices.pd
		param_pb = prices.pb
		param_z = z
		param_mfrac = (value.(M) + IOlib.ϵ * param_mfrac) ./ (value.(qd) + value.(F) + value.(I_supply) .+ IOlib.ϵ)

		for i in 1:ns
			set_normalized_coefficient(eq_io[i], u[i], -param_Pg * param_z[i])
			for j in 1:np
				set_normalized_coefficient(eq_io[i], qs[j], io.S[i,j] * param_pb[j])
			end
		end
		for i in 1:np
			set_normalized_coefficient(eq_X[i], xshare[i], -param_Xmax[i])
			set_normalized_coefficient(eq_F[i], fshare[i], -param_Fmax[i])
			set_normalized_coefficient(eq_M[i], ψ_pos[i], -param_Mref[i])
			set_normalized_coefficient(eq_M[i], ψ_neg[i], param_Mref[i])
			set_normalized_coefficient(eq_M[i], qd[i], -param_mfrac[i])
			set_normalized_coefficient(eq_M[i], F[i], -param_mfrac[i])
			set_normalized_coefficient(eq_M[i], I_supply[i], -param_mfrac[i])
			for j in 1:ns
				set_normalized_coefficient(eq_intdmd[i], u[j], -io.D[i,j] * param_z[j])
			end
		end
		fix(I_tot, param_I_tot)

		if params["report-diagnostics"]
			open(joinpath(params["diagnostics_path"], string("model_", run, "_", year, ".txt")), "w") do f
				print(f, mdl)
			end
		end
		optim_output = @capture_out begin
			optimize!(mdl)
		end
		if !isempty(optim_output)
			@info optim_output
		end
        status = primal_status(mdl)
        @info "Simulating for $year: $status"
        previous_failed = status != MOI.FEASIBLE_POINT
        if previous_failed
			finndx = length(LEAP_indices) + 2 # Adds column for year and for GDP
            indices[t,2:finndx] = fill(NaN, (finndx - 2) + 1)
        end

		year += 1

    end

    # Make indices into indices
    indices_0 = indices[1,:]
    for t = 1:ntime
        indices[t,2:end] = indices[t,2:end] ./ indices_0[2:end]
    end
    open(joinpath(params["results_path"], string("indices_",run,".csv")), "w") do io
               writedlm(io, reshape(labels, 1, :), ',')
               writedlm(io, indices, ',')
           end;

	# Release the model
	mdl = nothing

    return indices
end # ModelCalculations

"""
    resultcomparison(params, run::Int64)

Compare the GDP and value added indices generated from different runs and calculate the maximum difference.
Indices are specified in the YAML configuration file.
"""
function resultcomparison(params, run::Int64)
    # Obtains data from indices files for current run and previous run
    file1 = joinpath(params["results_path"], string("indices_",run,".csv"))
    file2 = joinpath(params["results_path"], string("indices_",run-1,".csv"))

    file1_dat = CSV.read(file1, header=1, missingstring=["0"], DataFrame)
    file2_dat = CSV.read(file2, header=1, missingstring=["0"], DataFrame)

    max_diff = 0

    # Compares run data
    for i = 2:size(file1_dat, 2)
        diff = abs.((file1_dat[:,i]-file2_dat[:,i])./file2_dat[:,i]*100)
        deleteat!(diff, findall(x->isnan(x), diff))
        if maximum(diff) > max_diff
            max_diff = maximum(diff)
        end
    end

    return max_diff
end

"""
    runleapmacromodel(file::String, logile::IOStream, include_energy_sectors::Bool = false)

Iteratively run the Macro model and LEAP until convergence.
"""
function runleapmacromodel(file::String, logfile::IOStream, include_energy_sectors::Bool = false)

    ## get base_year and final_year, and force fresh start with global_params
    params = IOlib.parse_input_file(file, force = true, include_energy_sectors = include_energy_sectors)
    base_year = params["years"]["start"]
    final_year = params["years"]["end"]
    ntime = 1 + (final_year - base_year)

    # set model run parameters and initial values
    I_en = zeros(ntime, 1)
    run_leap = params["model"]["run_leap"]
    if run_leap
        max_runs = params["model"]["max_runs"]
        ## checks that user has LEAP installed
        if ismissing(LEAPfunctions.connecttoleap())
            @error "Cannot connect to LEAP. Please check that LEAP is installed, or set 'run_leap: false' in the configuration file."
            return
        end
    else
        max_runs = 0
    end
    max_tolerance = params["model"]["max_tolerance"]

    for run = 0:max_runs
        ## Run Macro model
        #------------status
		print("Macro model run ($run)...")
        @info "Macro model run ($run)..."
        #------------status
        indices = ModelCalculations(file, I_en, run)

        ## Compare run results
        if run >= 1
            tolerance = resultcomparison(params, run)
            if tolerance <= max_tolerance
                @info @sprintf("Convergence in run number %d at %.2f%% ≤ %.2f%% target...", run, tolerance, max_tolerance)
                return
            else
                @info @sprintf("Results did not converge: %.2f%% > %.2f%% target...", tolerance, max_tolerance)
                if run == max_runs
                    return
                end
            end
        end

        ## Send output to LEAP
        if run_leap
			if params["model"]["hide_leap"]
				#------------status
				@info "Hiding LEAP to improve performance..."
				#------------status
				LEAPfunctions.visible(false)
			end

            #------------status
            @info "Sending Macro output to LEAP..."
            #------------status
            LEAPfunctions.outputtoleap(file, indices, run)
            ## Run LEAP model
			try
				#------------status
				@info "Running LEAP model..."
				flush(logfile)
				#------------status
				LEAPfunctions.calculateleap(params["LEAP-info"]["result_scenario"])

				## Obtain energy investment data from LEAP
				#------------status
				@info "Obtaining LEAP results..."
				flush(logfile)
				#------------status
				I_en = LEAPfunctions.energyinvestment(file, run)
			finally
				if params["model"]["hide_leap"]
					#------------status
					@info "Restoring LEAP..."
					#------------status
					LEAPfunctions.visible(true)
				end
				flush(logfile)
			end
        end
		#------------status
		println("completed")
        #------------status

    end

end # runleapmacromodel

end

max_stars_repo_name: sei-international/LEAPMacro

max_issues_repo_name: sei-international/LEAPMacro

max_forks_repo_name: sei-international/LEAPMacro


text:
using CompScienceMeshes
using BEAST
using LinearAlgebra
using FastBEAST

CM = CompScienceMeshes
function test_beast_laplace_singlelayer(h)
    Γ = CM.meshsphere(1, h)

    X = lagrangecxd0(Γ)

    𝒱 = Helmholtz3D.singlelayer(wavenumber=0.0)

    hmat = hassemble(𝒱,X,X, treeoptions = BoxTreeOptions(nmin=100))

    return  hmat
end

hmat = test_beast_laplace_singlelayer(0.1)

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
#
# Details on the Google Sheets API can be found at:
# https://developers.google.com/sheets/api/guides/concepts
#
# The quickstart guide is also useful
# https://developers.google.com/sheets/api/quickstart/python
#

"""
A package for working with Google Sheets.
"""
module GoogleSheets

using PyCall
using RateLimiter
using JSON
import MacroTools
import DataFrames: DataFrame, nrow, ncol, names
import ColorTypes: Colorant, RGBA, red, green, blue, alpha
using Colors

include("enums.jl")
include("types.jl")
include("json.jl")
include("client.jl")
include("spreadsheet.jl")
include("sheet.jl")
include("meta.jl")
include("io.jl")
include("format.jl")

#TODO: add chart -> https://developers.google.com/sheets/api/reference/rest/v4/spreadsheets/charts
#TODO: add filterview -> batch_update -> https://developers.google.com/sheets/api/guides/filters

end # module


text:
function f(x; y=0)
  x * y
end

println(f(10))
println(f(10, y=10))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
"""

```
sphericalshellvolume(i,step)
```

Computes the volume of the spherical shell defined within [(i-1)*step,i*step].

"""
function sphericalshellvolume(i,step)
  rmin = (i-1)*step
  return (4*pi/3)*( (rmin+step)^3 - rmin^3 )
end


text:
module Sudoku

using Random

export Board
export solve
export solve!

"""
	`Board(array)`

A struct that represents a Sudoku board. A 0 means the space is unknown and spaces with 1 through 9 are tiles with those numbers.

A 9×9 `array` can be used to construct.

	`Board(string)`

A board can be constructed with a `string` that is 81 in length an of only digit characters.
The string is read to fill the rows first from left to right, starting with the left most
character in the string first.

	ex. `Board("123456789890000000...")` = [[1 2 3 4 5 6 7 8 9]; [8 9 0 0 0 0 0 0 0 ]; ...]

Board properties...
`tiles` = A 9×9 array of Int64
"""
struct Board
	tiles::Array{Int64,2}
end

"""
	`Board()`

	`Board(; state=:unkown)` creates an board with all zeros (thus a board with no known spaces.)
"""
function Board(; state=:unknown)
	(state == :unknown) && return Board(zeros(Int64, 9, 9))

	if state == :random_solved
		board = Board(; state=:unknown)
		solve!(board;shuffle=true)
		return board
	end

	error("The there is not state option, $state.")
end

function Board(tile_string::AbstractString)
	length(tile_string) != 81 && error("The tile string, $tile_string, must of length 81.")

	for char in tile_string
		isdigit(char) || error("Only digits in your tile string, $tile_string.")
	end

	board = Board()

	for i in 1:81
		element = parse(Int64,tile_string[i])
		set_element!(board,(i-1)÷9 + 1, ((i-1) % 9) + 1, element)
	end

	return board
end

function get_element(board::Board, row_num::Int64,col_num::Int64)
	return board.tiles[row_num,col_num]
end

function get_row(board::Board, row_num::Int64)
	return board.tiles[row_num,:]
end

function get_column(board::Board, col_num::Int64)
	return board.tiles[:,col_num]
end

to_block_index(index::Int64) = index + 1 - ((index - 1) % 3)

function get_block(board::Board, row_num::Int64,col_num::Int64)
	block_row = to_block_index(row_num)
	block_col = to_block_index(col_num)
	return (board.tiles[block_row-1:block_row+1,block_col-1:block_col+1])[:]
end

function is_valid(board::Board)

	for i∈1:9,k ∈ 1:9
			count(==(k),get_row(board,i)) > 1 && return false
	end

	for i∈1:9,k ∈ 1:9
			count(==(k),get_column(board,i)) > 1 && return false
	end

	for i∈2:3:8,j∈2:3:8,k ∈ 1:9
		count(==(k),get_block(board,i,j)) > 1 && return false
	end

	return true
end

function is_solved(board::Board)
	return 0 ∉ board.tiles
end

function set_element!(board::Board, row_num::Int64, col_num::Int64, element::Int64)
	board.tiles[row_num,col_num] = element
end

function possible_moves(board::Board,row_num::Int64,col_num::Int64; shuffle = false)
	element = get_element(board,row_num, col_num)

	element == 0 || return Int64[]

	moves = [i for i ∈ 1:9 if i∉get_row(board, row_num)
													 && i∉get_column(board, col_num)
													 && i∉get_block(board, row_num, col_num)]

	shuffle && Random.shuffle!(moves)

	return moves
end

"""
	`solve!(board::Board)`

	Changes solve to a solved board. This function returns nothing.
	The algorithim is a depth search.

	WARNING: If there are multiple solutions, this algorithm will not find both of them.
"""
function solve!(board::Board; shuffle = false)
	index = 1
	backtrack_dict = Dict()

	unknown_tiles = [(i,j) for i∈1:9,j∈1:9 if get_element(board, i, j) == 0][:]

	while !is_solved(board) && index <= length(unknown_tiles) && index > 0

		row, column = unknown_tiles[index]

		if !haskey(backtrack_dict, (row, column))
			push!(backtrack_dict, (row, column) => possible_moves(board, row, column; shuffle = shuffle))
		end

		if isempty(backtrack_dict[(row,column)])
			delete!(backtrack_dict, (row,column))
			set_element!(board,row, column, 0)
			index -= 1
			continue
		end

		next_move = popfirst!(backtrack_dict[(row, column)])
		set_element!(board, row, column, next_move)
		index += 1

	end
end

"""
	`solve(tile_string::AbstractString)`

	`solve(tile_string)` returns a solved board, starting from the board given as a string.
"""
function solve(tile_string::AbstractString; shuffle = false)
	board = Board(tile_string)
	solve!(board; shuffle = shuffle)
	return board
end

"""
	`solve(board::Board)`

	`solve(board)` returns a solved board, starting from a generic board.
"""
function solve(board::Board; shuffle= false)
	solved_board = board
	solve!(solved_board; shuffle = shuffle)
	return solved_board
end

end # module


text:
using Test
using Spandex

@testset "norm!" begin
    local a = from_csr(
        10,
        10,
        [1, 2, 3, 6, 7, 9, 10, 11, 12, 16, 19],
        [1, 2, 1, 2, 3, 4, 4, 5, 6, 7, 8, 6, 7, 8, 9, 8, 9, 10],
        [
            0.32502429653806530,
            0.68646720167973219,
            0.20664713231230356,
            0.80911717430065078,
            1.0850649558194219,
            0.37484683221445758,
            0.29211785639326571,
            0.56725008142048217,
            0.30569332800228499,
            0.12993461529471456,
            0.48493422490889138,
            0.22043327522151054,
            0.056305455607475616,
            0.33290016729588462,
            0.70740209966807233,
            0.45713455373057488,
            0.73539821745905609,
            1.5818072373358281,
        ],
    )

    local norm = Spandex.norm!(a)
    local e = [
        1.7540504746792074,
        1.2069522923516747,
        0.96000198563023442,
        1.6333267599804464,
        1.3277390136750080,
        1.8086603584984466,
        2.7741987233555063,
        1.4360136957244634,
        1.1889588605089314,
        0.79510268651896931,
    ]

    norm = round.(norm, digits = 16)
    e = round.(e, digits = 16)

    @test e == norm
end


text:
module EscapeAnalysis

let
    README = normpath(dirname(@__DIR__), "README.md")
    include_dependency(README)
    @doc read(README, String) EscapeAnalysis
end

export
    analyze_escapes,
    @analyze_escapes

import Core:
    CodeInfo,
    CodeInstance,
    MethodInstance,
    Const,
    Argument,
    SSAValue,
    PiNode,
    PhiNode,
    UpsilonNode,
    PhiCNode,
    ReturnNode,
    GotoNode,
    GotoIfNot,
    SimpleVector

const CC = Core.Compiler

import .CC:
    AbstractInterpreter,
    NativeInterpreter,
    WorldView,
    WorldRange,
    InferenceParams,
    OptimizationParams,
    get_world_counter,
    get_inference_cache,
    lock_mi_inference,
    unlock_mi_inference,
    add_remark!,
    may_optimize,
    may_compress,
    may_discard_trees,
    verbose_stmt_info,
    code_cache,
    get_inference_cache,
    OptimizationState,
    IRCode,
    optimize,
    widenconst,
    argextype,
    singleton_type,
    IR_FLAG_EFFECT_FREE,
    is_meta_expr_head

import Base.Meta:
    isexpr

import Base:
    destructure_callex

using InteractiveUtils

let __init_hooks__ = []
    global __init__() = foreach(f->f(), __init_hooks__)
    global register_init_hook!(@nospecialize(f)) = push!(__init_hooks__, f)
end

mutable struct EscapeAnalyzer{State} <: AbstractInterpreter
    native::NativeInterpreter
    ir::IRCode
    state::State
    linfo::MethodInstance
    EscapeAnalyzer(native::NativeInterpreter) = new{EscapeState}(native)
end

CC.InferenceParams(interp::EscapeAnalyzer)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::EscapeAnalyzer) = OptimizationParams(interp.native)
CC.get_world_counter(interp::EscapeAnalyzer)  = get_world_counter(interp.native)

CC.lock_mi_inference(::EscapeAnalyzer,   ::MethodInstance) = nothing
CC.unlock_mi_inference(::EscapeAnalyzer, ::MethodInstance) = nothing

CC.add_remark!(interp::EscapeAnalyzer, sv, s) = add_remark!(interp.native, sv, s)

CC.may_optimize(interp::EscapeAnalyzer)      = may_optimize(interp.native)
CC.may_compress(interp::EscapeAnalyzer)      = may_compress(interp.native)
CC.may_discard_trees(interp::EscapeAnalyzer) = may_discard_trees(interp.native)
CC.verbose_stmt_info(interp::EscapeAnalyzer) = verbose_stmt_info(interp.native)

CC.get_inference_cache(interp::EscapeAnalyzer) = get_inference_cache(interp.native)

const GLOBAL_CODE_CACHE = IdDict{MethodInstance,CodeInstance}()
__clear_code_cache!() = empty!(GLOBAL_CODE_CACHE)

function CC.code_cache(interp::EscapeAnalyzer)
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(GlobalCache(), worlds)
end

struct GlobalCache end

CC.haskey(wvc::WorldView{GlobalCache}, mi::MethodInstance) = haskey(GLOBAL_CODE_CACHE, mi)

CC.get(wvc::WorldView{GlobalCache}, mi::MethodInstance, default) = get(GLOBAL_CODE_CACHE, mi, default)

CC.getindex(wvc::WorldView{GlobalCache}, mi::MethodInstance) = getindex(GLOBAL_CODE_CACHE, mi)

function CC.setindex!(wvc::WorldView{GlobalCache}, ci::CodeInstance, mi::MethodInstance)
    GLOBAL_CODE_CACHE[mi] = ci
    add_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_cache!)
        end
    end
    return nothing
end

function invalidate_cache!(replaced, max_world, depth = 0)
    delete!(GLOBAL_CODE_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !haskey(GLOBAL_CODE_CACHE, mi)
                continue # otherwise fall into infinite loop
            end
            invalidate_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# analysis
# ========

"""
    x::EscapeLattice

A lattice for escape information, which holds the following properties:
- `x.Analyzed::Bool`: not formally part of the lattice, indicates this statement has not been analyzed at all
- `x.ReturnEscape::BitSet`: keeps SSA numbers of return statements where it can be returned to the caller
  * `isempty(x.ReturnEscape)` means it never escapes to the caller
  * otherwise it indicates it will escape to the caller via return (possibly as a field),
    where `0 ∈ x.ReturnEscape` has the special meaning that it's visible to the caller
    simply because it's passed as call argument
- `x.ThrownEscape::Bool`: indicates it may escape to somewhere through an exception (possibly as a field)
- `x.GlobalEscape::Bool`: indicates it may escape to a global space an exception (possibly as a field)
- `x.ArgEscape::Int` (not implemented yet): indicates it will escape to the caller through `setfield!` on argument(s)
  * `-1` : no escape
  * `0` : unknown or multiple
  * `n` : through argument N

These attributes can be combined to create a partial lattice that has a finite height, given
that input program has a finite number of statements, which is assured by Julia's semantics.

There are utility constructors to create common `EscapeLattice`s, e.g.,
- `NoEscape()`: the bottom element of this lattice, meaning it won't escape to anywhere
- `AllEscape()`: the topmost element of this lattice, meaning it will escape to everywhere

The escape analysis will transition these elements from the bottom to the top,
in the same direction as Julia's native type inference routine.
An abstract state will be initialized with the bottom(-like) elements:
- the call arguments are initialized as `ArgumentReturnEscape()`, because they're visible from a caller immediately
- the other states are initialized as `NotAnalyzed()`, which is a special lattice element that
  is slightly lower than `NoEscape`, but at the same time doesn't represent any meaning
  other than it's not analyzed yet (thus it's not formally part of the lattice).
"""
struct EscapeLattice
    Analyzed::Bool
    ReturnEscape::BitSet
    ThrownEscape::Bool
    GlobalEscape::Bool
    # TODO: ArgEscape::Int
end

function Base.:(==)(x::EscapeLattice, y::EscapeLattice)
    return x.Analyzed === y.Analyzed &&
           x.ReturnEscape == y.ReturnEscape &&
           x.ThrownEscape === y.ThrownEscape &&
           x.GlobalEscape === y.GlobalEscape
end

# lattice constructors
# precompute default values in order to eliminate computations at callsites
const NO_RETURN = BitSet()
const ARGUMENT_RETURN = BitSet(0)
NotAnalyzed() = EscapeLattice(false, NO_RETURN, false, false) # not formally part of the lattice
NoEscape() = EscapeLattice(true, NO_RETURN, false, false)
ReturnEscape(pcs::BitSet) = EscapeLattice(true, pcs, false, false)
ReturnEscape(pc::Int) = ReturnEscape(BitSet(pc))
ArgumentReturnEscape() = ReturnEscape(ARGUMENT_RETURN)
ThrownEscape() = EscapeLattice(true, NO_RETURN, true, false)
GlobalEscape() = EscapeLattice(true, NO_RETURN, false, true)
let
    all_return = BitSet(0:100_000)
    global AllReturnEscape() = ReturnEscape(all_return) # used for `show`
    global AllEscape() = EscapeLattice(true, all_return, true, true)
end

# Convenience names for some ⊑ queries
export
    has_not_analyzed,
    has_no_escape,
    has_return_escape,
    has_thrown_escape,
    has_global_escape,
    has_all_escape,
    can_elide_finalizer
has_not_analyzed(x::EscapeLattice) = x == NotAnalyzed()
has_no_escape(x::EscapeLattice) = x ⊑ NoEscape()
has_return_escape(x::EscapeLattice) = !isempty(x.ReturnEscape)
has_return_escape(x::EscapeLattice, pc::Int) = pc in x.ReturnEscape
has_thrown_escape(x::EscapeLattice) = x.ThrownEscape
has_global_escape(x::EscapeLattice) = x.GlobalEscape
has_all_escape(x::EscapeLattice) = AllEscape() == x

"""
    can_elide_finalizer(x::EscapeLattice, pc::Int) -> Bool

Queries the validity of the finalizer elision optimization at the `return` site of statement `pc`,
which inserts `finalize` call when the lifetime of interested object ends.
Note that we don't need to take `x.ThrownEscape` into account because it would have never
been thrown when the program execution reaches the `return` site.
"""
function can_elide_finalizer(x::EscapeLattice, pc::Int)
    x.GlobalEscape && return false
    0 in x.ReturnEscape && return false
    return pc ∉ x.ReturnEscape
end

function ⊑(x::EscapeLattice, y::EscapeLattice)
    if x.Analyzed ≤ y.Analyzed &&
       x.ReturnEscape ⊆ y.ReturnEscape &&
       x.ThrownEscape ≤ y.ThrownEscape &&
       x.GlobalEscape ≤ y.GlobalEscape
       return true
    end
    return false
end
⋤(x::EscapeLattice, y::EscapeLattice) = ⊑(x, y) && !⊑(y, x)

function ⊔(x::EscapeLattice, y::EscapeLattice)
    return EscapeLattice(
        x.Analyzed | y.Analyzed,
        x.ReturnEscape ∪ y.ReturnEscape,
        x.ThrownEscape | y.ThrownEscape,
        x.GlobalEscape | y.GlobalEscape,
        )
end

function ⊓(x::EscapeLattice, y::EscapeLattice)
    return EscapeLattice(
        x.Analyzed & y.Analyzed,
        x.ReturnEscape ∩ y.ReturnEscape,
        x.ThrownEscape & y.ThrownEscape,
        x.GlobalEscape & y.GlobalEscape,
        )
end

# TODO setup a more effient struct for cache
# which can discard escape information on SSS values and arguments that don't join dispatch signature

"""
    state::EscapeState

Extended lattice that maps arguments and SSA values to escape information represented as `EscapeLattice`:
- `state.arguments::Vector{EscapeLattice}`: escape information about "arguments" – note that
  "argument" can include both call arguments and slots appearing in analysis frame
- `ssavalues::Vector{EscapeLattice}`: escape information about each SSA value
"""
struct EscapeState
    arguments::Vector{EscapeLattice}
    ssavalues::Vector{EscapeLattice}
end
function EscapeState(nslots::Int, nargs::Int, nstmts::Int)
    arguments = EscapeLattice[
        1 ≤ i ≤ nargs ? ArgumentReturnEscape() : NotAnalyzed() for i in 1:nslots]
    ssavalues = EscapeLattice[NotAnalyzed() for _ in 1:nstmts]
    return EscapeState(arguments, ssavalues)
end

# we preserve `IRCode` as well just for debugging purpose
const GLOBAL_ESCAPE_CACHE = IdDict{MethodInstance,Tuple{EscapeState,IRCode}}()
__clear_escape_cache!() = empty!(GLOBAL_ESCAPE_CACHE)

const Change  = Pair{Union{Argument,SSAValue},EscapeLattice}
const Changes = Vector{Change}

"""
    find_escapes(ir::IRCode, nargs::Int) -> EscapeState

Escape analysis implementation is based on the data-flow algorithm described in the paper [^MM02].
The analysis works on the lattice of [`EscapeLattice`](@ref) and transitions lattice elements
from the bottom to the top in a _backward_ way, i.e. data flows from usage cites to definitions,
until every lattice gets converged to a fixed point by maintaining a (conceptual) working set
that contains program counters corresponding to remaining SSA statements to be analyzed.
Note that the analysis only manages a single global state, with some flow-sensitivity
encoded as property of `EscapeLattice`.

[^MM02]: A Graph-Free approach to Data-Flow Analysis.
         Markas Mohnen, 2002, April.
         <https://api.semanticscholar.org/CorpusID:28519618>
"""
function find_escapes(ir::IRCode, nargs::Int)
    (; stmts, sptypes, argtypes) = ir
    nstmts = length(stmts)

    # only manage a single state, some flow-sensitivity is encoded as `EscapeLattice` properties
    state = EscapeState(length(ir.argtypes), nargs, nstmts)
    changes = Changes() # stashes changes that happen at current statement

    while true
        local anyupdate = false

        for pc in nstmts:-1:1
            stmt = stmts.inst[pc]

            # we escape statements with the `ThrownEscape` property using the effect-freeness
            # information computed by the inliner
            is_effect_free = stmts.flag[pc] & IR_FLAG_EFFECT_FREE ≠ 0

            # collect escape information
            if isa(stmt, Expr)
                head = stmt.head
                if head === :call
                    has_changes = escape_call!(stmt.args, pc, state, ir, changes)
                    if !is_effect_free
                        add_changes!(stmt.args, ir, ThrownEscape(), changes)
                    else
                        has_changes || continue
                    end
                elseif head === :invoke
                    escape_invoke!(stmt.args, pc, state, ir, changes)
                elseif head === :new
                    info = state.ssavalues[pc]
                    if info == NotAnalyzed()
                        info = NoEscape()
                        add_change!(SSAValue(pc), ir, info, changes) # we will be interested in if this allocation escapes or not
                    end
                    add_changes!(stmt.args[2:end], ir, info, changes)
                elseif head === :splatnew
                    info = state.ssavalues[pc]
                    if info == NotAnalyzed()
                        info = NoEscape()
                        add_change!(SSAValue(pc), ir, info, changes) # we will be interested in if this allocation escapes or not
                    end
                    # splatnew passes field values using a single tuple (args[2])
                    add_change!(stmt.args[2], ir, info, changes)
                elseif head === :(=)
                    lhs, rhs = stmt.args
                    if isa(lhs, GlobalRef) # global store
                        add_change!(rhs, ir, GlobalEscape(), changes)
                    end
                elseif head === :foreigncall
                    # for foreigncall we simply escape every argument (args[6:length(args[3])])
                    # and its name (args[1])
                    # TODO: we can apply a similar strategy like builtin calls to specialize some foreigncalls
                    foreigncall_nargs = length((stmt.args[3])::SimpleVector)
                    name = stmt.args[1]
                    # if normalize(name) === :jl_gc_add_finalizer_th
                    #     continue # XXX assume this finalizer call is valid for finalizer elision
                    # end
                    add_change!(name, ir, ThrownEscape(), changes)
                    add_changes!(stmt.args[6:5+foreigncall_nargs], ir, ThrownEscape(), changes)
                elseif head === :throw_undef_if_not # XXX when is this expression inserted ?
                    add_change!(stmt.args[1], ir, ThrownEscape(), changes)
                elseif is_meta_expr_head(head)
                    # meta expressions doesn't account for any usages
                    continue
                elseif head === :static_parameter
                    # :static_parameter refers any of static parameters, but since they exist
                    # statically, we're really not interested in their escapes
                    continue
                elseif head === :copyast
                    # copyast simply copies a surface syntax AST, and should never use any of arguments or SSA values
                    continue
                elseif head === :undefcheck
                    # undefcheck is temporarily inserted by compiler
                    # it will be processd be later pass so it won't change any of escape states
                    continue
                elseif head === :the_exception
                    # we don't propagate escape information on exceptions via this expression, but rather
                    # use a dedicated lattice property `ThrownEscape`
                    continue
                elseif head === :isdefined
                    # just returns `Bool`, nothing accounts for any usages
                    continue
                elseif head === :enter || head === :leave || head === :pop_exception
                    # these exception frame managements doesn't account for any usages
                    # we can just ignore escape information from
                    continue
                elseif head === :gc_preserve_begin || head === :gc_preserve_end
                    # `GC.@preserve` may "use" arbitrary values, but we can just ignore the escape information
                    # imposed on `GC.@preserve` expressions since they're supposed to never be used elsewhere
                    continue
                else
                    add_changes!(stmt.args, ir, AllEscape(), changes)
                end
            elseif isa(stmt, GlobalRef) # global load
                add_change!(SSAValue(pc), ir, GlobalEscape(), changes)
            elseif isa(stmt, PiNode)
                if isdefined(stmt, :val)
                    info = state.ssavalues[pc]
                    add_change!(stmt.val, ir, info, changes)
                end
            elseif isa(stmt, PhiNode)
                info = state.ssavalues[pc]
                values = stmt.values
                for i in 1:length(values)
                    if isassigned(values, i)
                        add_change!(values[i], ir, info, changes)
                    end
                end
            elseif isa(stmt, PhiCNode)
                info = state.ssavalues[pc]
                values = stmt.values
                for i in 1:length(values)
                    if isassigned(values, i)
                        add_change!(values[i], ir, info, changes)
                    end
                end
            elseif isa(stmt, UpsilonNode)
                if isdefined(stmt, :val)
                    info = state.ssavalues[pc]
                    add_change!(stmt.val, ir, info, changes)
                end
            elseif isa(stmt, ReturnNode)
                if isdefined(stmt, :val)
                    add_change!(stmt.val, ir, ReturnEscape(pc), changes)
                end
            else
                @assert stmt isa GotoNode || stmt isa GotoIfNot || stmt isa GlobalRef || isnothing(stmt) # TODO remove me
                continue
            end

            isempty(changes) && continue

            anyupdate |= propagate_changes!(state, changes)

            empty!(changes)
        end

        anyupdate || break
    end

    return state
end

# propagate changes, and check convergence
function propagate_changes!(state::EscapeState, changes::Changes)
    local anychanged = false

    for (x, info) in changes
        if isa(x, Argument)
            old = state.arguments[x.n]
            new = old ⊔ info
            if old ≠ new
                state.arguments[x.n] = new
                anychanged |= true
            end
        else
            x = x::SSAValue
            old = state.ssavalues[x.id]
            new = old ⊔ info
            if old ≠ new
                state.ssavalues[x.id] = new
                anychanged |= true
            end
        end
    end

    return anychanged
end

# function normalize(@nospecialize(x))
#     if isa(x, QuoteNode)
#         return x.value
#     else
#         return x
#     end
# end

function add_changes!(args::Vector{Any}, ir::IRCode, info::EscapeLattice, changes::Changes)
    for x in args
        add_change!(x, ir, info, changes)
    end
end

function add_change!(@nospecialize(x), ir::IRCode, info::EscapeLattice, changes::Changes)
    if isa(x, Argument) || isa(x, SSAValue)
        if !isbitstype(widenconst(argextype(x, ir, ir.sptypes, ir.argtypes)))
            push!(changes, Change(x, info))
        end
    end
end

function escape_invoke!(args::Vector{Any}, pc::Int,
                        state::EscapeState, ir::IRCode, changes::Changes)
    linfo = first(args)::MethodInstance
    cache = get(GLOBAL_ESCAPE_CACHE, linfo, nothing)
    args = args[2:end]
    if isnothing(cache)
        add_changes!(args, ir, AllEscape(), changes)
    else
        (linfostate, _ #=ir::IRCode=#) = cache
        retinfo = state.ssavalues[pc] # escape information imposed on the call statement
        method = linfo.def::Method
        nargs = Int(method.nargs)
        for i in 1:length(args)
            arg = args[i]
            if i ≤ nargs
                arginfo = linfostate.arguments[i]
            else # handle isva signature: COMBAK will this invalid once we encode alias information ?
                arginfo = linfostate.arguments[nargs]
            end
            if isempty(arginfo.ReturnEscape)
                @eval Main (ir = $ir; linfo = $linfo)
                error("invalid escape lattice element returned from inter-procedural context: inspect `Main.ir` and `Main.linfo`")
            end
            info = from_interprocedural(arginfo, retinfo)
            add_change!(arg, ir, info, changes)
        end
    end
end

# reinterpret the escape information imposed on the callee argument (`arginfo`) in the
# context of the caller frame using the escape information imposed on the return value (`retinfo`)
function from_interprocedural(arginfo::EscapeLattice, retinfo::EscapeLattice)
    ar = arginfo.ReturnEscape
    newarginfo = EscapeLattice(true, NO_RETURN, arginfo.ThrownEscape, arginfo.GlobalEscape)
    if ar == ARGUMENT_RETURN
        # if this is simply passed as the call argument, we can discard the `ReturnEscape`
        # information and just propagate the other escape information
        return newarginfo
    else
        # if this can be a return value, we have to merge it with the escape information
        return newarginfo ⊔ retinfo
    end
end

function escape_call!(args::Vector{Any}, pc::Int,
                      state::EscapeState, ir::IRCode, changes::Changes)
    ft = argextype(first(args), ir, ir.sptypes, ir.argtypes)
    f = singleton_type(ft)
    if isa(f, Core.IntrinsicFunction)
        return false # COMBAK we may break soundness here, e.g. `pointerref`
    end
    ishandled = escape_builtin!(f, args, pc, state, ir, changes)::Union{Nothing,Bool}
    isnothing(ishandled) && return false # nothing to propagate
    if !ishandled
        # if this call hasn't been handled by any of pre-defined handlers,
        # we escape this call conservatively
        add_changes!(args[2:end], ir, AllEscape(), changes)
    end
    return true
end

# TODO implement more builtins, make them more accurate
# TODO use `T_IFUNC`-like logic and don't not abuse dispatch ?

escape_builtin!(@nospecialize(f), _...) = return false

escape_builtin!(::typeof(isa), _...) = return nothing
escape_builtin!(::typeof(typeof), _...) = return nothing
escape_builtin!(::typeof(Core.sizeof), _...) = return nothing
escape_builtin!(::typeof(===), _...) = return nothing

function escape_builtin!(::typeof(ifelse), args::Vector{Any}, pc::Int, state::EscapeState, ir::IRCode, changes::Changes)
    length(args) == 4 || return false
    f, cond, th, el = args
    info = state.ssavalues[pc]
    condt = argextype(cond, ir, ir.sptypes, ir.argtypes)
    if isa(condt, Const) && (cond = condt.val; isa(cond, Bool))
        if cond
            add_change!(th, ir, info, changes)
        else
            add_change!(el, ir, info, changes)
        end
    else
        add_change!(th, ir, info, changes)
        add_change!(el, ir, info, changes)
    end
    return true
end

function escape_builtin!(::typeof(typeassert), args::Vector{Any}, pc::Int, state::EscapeState, ir::IRCode, changes::Changes)
    length(args) == 3 || return false
    f, obj, typ = args
    info = state.ssavalues[pc]
    add_change!(obj, ir, info, changes)
    return true
end

function escape_builtin!(::typeof(tuple), args::Vector{Any}, pc::Int, state::EscapeState, ir::IRCode, changes::Changes)
    info = state.ssavalues[pc]
    if info == NotAnalyzed()
        info = NoEscape()
    end
    add_changes!(args[2:end], ir, info, changes)
    return true
end

# TODO don't propagate escape information to the 1st argument, but propagate information to aliased field
function escape_builtin!(::typeof(getfield), args::Vector{Any}, pc::Int, state::EscapeState, ir::IRCode, changes::Changes)
    info = state.ssavalues[pc]
    if info == NotAnalyzed()
        info = NoEscape()
    end
    # only propagate info when the field itself is non-bitstype
    if !isbitstype(widenconst(ir.stmts.type[pc]))
        add_changes!(args[2:end], ir, info, changes)
    end
    return true
end

# entries
# =======

function CC.optimize(interp::EscapeAnalyzer, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    ir = run_passes_with_escape_analysis(interp, opt.src, opt)
    return CC.finish(interp, opt, params, ir, result)
end

# TODO implement finalizer elision optimization
function elide_finalizers(ir::IRCode, state::EscapeState)
    return ir
end

# HACK enable copy and paste from Core.Compiler
function run_passes_with_escape_analysis end
register_init_hook!() do
@eval CC begin
    function $EscapeAnalysis.run_passes_with_escape_analysis(interp::$EscapeAnalyzer, ci::CodeInfo, sv::OptimizationState)
        preserve_coverage = coverage_enabled(sv.mod)
        ir = convert_to_ircode(ci, copy_exprargs(ci.code), preserve_coverage, sv)
        ir = slot2reg(ir, ci, sv)
        #@Base.show ("after_construct", ir)
        # TODO: Domsorting can produce an updated domtree - no need to recompute here
        @timeit "compact 1" ir = compact!(ir)
        @timeit "Inlining" ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
        #@timeit "verify 2" verify_ir(ir)
        ir = compact!(ir)
        svdef = sv.linfo.def
        nargs = isa(svdef, Method) ? Int(svdef.nargs) : 0
        @timeit "collect escape information" state = $find_escapes(ir, nargs)
        cacheir = copy(ir)
        # cache
        $setindex!($GLOBAL_ESCAPE_CACHE, (state, cacheir), sv.linfo)
        # for entry result
        interp.ir = cacheir
        interp.state = state
        interp.linfo = sv.linfo
        @timeit "finalizer elision" ir = $elide_finalizers(ir, state)
        #@Base.show ("before_sroa", ir)
        @timeit "SROA" ir = getfield_elim_pass!(ir)
        #@Base.show ir.new_nodes
        #@Base.show ("after_sroa", ir)
        ir = adce_pass!(ir)
        #@Base.show ("after_adce", ir)
        @timeit "type lift" ir = type_lift_pass!(ir)
        @timeit "compact 3" ir = compact!(ir)
        #@Base.show ir
        if JLOptions().debug_level == 2
            @timeit "verify 3" (verify_ir(ir); verify_linetable(ir.linetable))
        end
        return ir
    end
end
end # register_init_hook!() do

macro analyze_escapes(ex0...)
    return InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :analyze_escapes, ex0)
end

function analyze_escapes(@nospecialize(f), @nospecialize(types=Tuple{});
                         world = get_world_counter(),
                         interp = Core.Compiler.NativeInterpreter(world))
    interp = EscapeAnalyzer(interp)
    results = code_typed(f, types; optimize=true, world, interp)
    isone(length(results)) || throw(ArgumentError("`analyze_escapes` only supports single analysis result"))
    return EscapeResult(interp.ir, interp.state, interp.linfo)
end

# utilities
# =========

# in order to run a whole analysis from ground zero (e.g. for benchmarking, etc.)
__clear_caches!() = (__clear_code_cache!(); __clear_escape_cache!())

function get_name_color(x::EscapeLattice, symbol::Bool = false)
    getname(x) = string(nameof(x))
    if x == NotAnalyzed()
        name, color = (getname(NotAnalyzed), '◌'), :plain
    elseif x == NoEscape()
        name, color = (getname(NoEscape), '✓'), :green
    elseif NoEscape() ⋤ x ⊑ AllReturnEscape()
        pcs = sprint(show, collect(x.ReturnEscape); context=:limit=>true)
        name1 = string(getname(ReturnEscape), '(', pcs, ')')
        name = name1, '↑'
        color = :cyan
    elseif NoEscape() ⋤ x ⊑ ThrownEscape()
        name, color = (getname(ThrownEscape), '↓'), :yellow
    elseif NoEscape() ⋤ x ⊑ GlobalEscape()
        name, color = (getname(GlobalEscape), 'G'), :red
    elseif x == AllEscape()
        name, color = (getname(AllEscape), 'X'), :red
    else
        name, color = (nothing, '*'), :red
    end
    return (symbol ? last(name) : first(name), color)
end

function Base.show(io::IO, x::EscapeLattice)
    name, color = get_name_color(x)
    if isnothing(name)
        Base.@invoke show(io::IO, x::Any)
    else
        printstyled(io, name; color)
    end
end
function Base.show(io::IO, ::MIME"application/prs.juno.inline", x::EscapeLattice)
    name, color = get_name_color(x)
    if isnothing(name)
        return x # use fancy tree-view
    else
        printstyled(io, name; color)
    end
end

struct EscapeResult
    ir::IRCode
    state::EscapeState
    linfo::Union{Nothing,MethodInstance}
    EscapeResult(ir::IRCode, state::EscapeState, linfo::Union{Nothing,MethodInstance} = nothing) =
        new(ir, state, linfo)
end
Base.show(io::IO, result::EscapeResult) = print_with_info(io, result.ir, result.state, result.linfo)
@eval Base.iterate(res::EscapeResult, state=1) =
    return state > $(fieldcount(EscapeResult)) ? nothing : (getfield(res, state), state+1)

# adapted from https://github.com/JuliaDebug/LoweredCodeUtils.jl/blob/4612349432447e868cf9285f647108f43bd0a11c/src/codeedges.jl#L881-L897
function print_with_info(io::IO,
    ir::IRCode, (; arguments, ssavalues)::EscapeState, linfo::Union{Nothing,MethodInstance})
    # print escape information on SSA values
    function preprint(io::IO)
        ft = ir.argtypes[1]
        f = singleton_type(ft)
        if f === nothing
            f = widenconst(ft)
        end
        print(io, f, '(')
        for (i, arg) in enumerate(arguments)
            i == 1 && continue
            c, color = get_name_color(arg, true)
            printstyled(io, '_', i, "::", ir.argtypes[i], ' ', c; color)
            i ≠ length(arguments) && print(io, ", ")
        end
        print(io, ')')
        if !isnothing(linfo)
            def = linfo.def
            printstyled(io, " in ", (isa(def, Module) ? (def,) : (def.module, " at ", def.file, ':', def.line))...; color=:bold)
        end
        println(io)
    end

    # print escape information on SSA values
    # nd = ndigits(length(ssavalues))
    function preprint(io::IO, idx::Int)
        c, color = get_name_color(ssavalues[idx], true)
        # printstyled(io, lpad(idx, nd), ' ', c, ' '; color)
        printstyled(io, c, ' '; color)
    end

    print_with_info(preprint, (args...)->nothing, io, ir)
end

function print_with_info(preprint, postprint, io::IO, ir::IRCode)
    io = IOContext(io, :displaysize=>displaysize(io))
    used = Base.IRShow.stmts_used(io, ir)
    # line_info_preprinter = Base.IRShow.lineinfo_disabled
    line_info_preprinter = function (io::IO, indent::String, idx::Int)
        r = Base.IRShow.inline_linfo_printer(ir)(io, indent, idx)
        idx ≠ 0 && preprint(io, idx)
        return r
    end
    line_info_postprinter = Base.IRShow.default_expr_type_printer
    preprint(io)
    bb_idx_prev = bb_idx = 1
    for idx = 1:length(ir.stmts)
        preprint(io, idx)
        bb_idx = Base.IRShow.show_ir_stmt(io, ir, idx, line_info_preprinter, line_info_postprinter, used, ir.cfg, bb_idx)
        postprint(io, idx, bb_idx != bb_idx_prev)
        bb_idx_prev = bb_idx
    end
    max_bb_idx_size = ndigits(length(ir.cfg.blocks))
    line_info_preprinter(io, " "^(max_bb_idx_size + 2), 0)
    postprint(io)
    return nothing
end

end # module EscapeAnalysis


text:
using CUDA
using ExaPF
using KernelAbstractions
using Test
using Printf

import ExaPF: PowerSystem, LinearSolvers, TimerOutputs

# For debugging in REPL use the following lines
# empty!(ARGS)
# push!(ARGS, "BICGSTAB")
# push!(ARGS, "CUDADevice")
# push!(ARGS, "caseGO30R-025.raw")

# We do need the time in ms, and not with time units all over the place
function ExaPF.TimerOutputs.prettytime(t)
    value = t / 1e6 # "ms"

    if round(value) >= 100
        str = string(@sprintf("%.0f", value))
    elseif round(value * 10) >= 100
        str = string(@sprintf("%.1f", value))
    elseif round(value * 100) >= 100
        str = string(@sprintf("%.2f", value))
    else
        str = string(@sprintf("%.3f", value))
    end
    return lpad(str, 6, " ")
end

function printtimer(timers, key::String)
   prettytime(timers[key].accumulated_data.time)
end

linsolver = eval(Meta.parse("LinearSolvers.$(ARGS[1])"))
device = eval(Meta.parse("$(ARGS[2])()"))
datafile = joinpath(dirname(@__FILE__), ARGS[3])
if endswith(datafile, ".m")
    pf = PowerSystem.PowerNetwork(datafile, 1)
else
    pf = PowerSystem.PowerNetwork(datafile)
end
# Parameters
tolerance = 1e-6
polar = PolarForm(pf, device)
jac = ExaPF._state_jacobian(polar)
@show size(jac)
@show npartitions = ceil(Int64,(size(jac,1)/64))
precond = ExaPF.LinearSolvers.BlockJacobiPreconditioner(jac, npartitions, device)
# Retrieve initial state of network
x0 = ExaPF.initial(polar, State())
uk = ExaPF.initial(polar, Control())
p = ExaPF.initial(polar, Parameters())

algo = linsolver(precond)
xk = copy(x0)
nlp = ExaPF.ReducedSpaceEvaluator(polar, xk, uk, p;
                                    ε_tol=tolerance, linear_solver=algo)
convergence = ExaPF.update!(nlp, uk; verbose_level=ExaPF.VERBOSE_LEVEL_HIGH)
nlp.x .= x0
convergence = ExaPF.update!(nlp, uk; verbose_level=ExaPF.VERBOSE_LEVEL_HIGH)
nlp.x .= x0
ExaPF.reset_timer!(ExaPF.TIMER)
convergence = ExaPF.update!(nlp, uk; verbose_level=ExaPF.VERBOSE_LEVEL_HIGH)

# Make sure we are converged
@assert(convergence.has_converged)

# Output
prettytime = ExaPF.TimerOutputs.prettytime
timers = ExaPF.TIMER.inner_timers
inner_timer = timers["Newton"]
println("$(ARGS[1]), $(ARGS[2]), $(ARGS[3]),",
        printtimer(timers, "Newton"),",",
        printtimer(inner_timer, "Jacobian"),",",
        printtimer(inner_timer, "Linear Solver"))