text:
using CapacityExpansion
using Clp
exact_scenarios=Dict{String,OptResult}()
for (state, years) in [["GER_1", [2016]],["CA_1", [2016]]]
       # laod data
       ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
       cep_data = load_cep_data_provided(state)
       ## CLUSTERING ##
       ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
       ts_full_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=30)
       ## OPTIMIZATION ##
       optimizer=Clp.Optimizer
       optimizer_config=Dict{Symbol,Any}(:LogLevel => 0)
       exact_scenarios["$state-$years-co2"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="co2",limit_emission=Dict{String,Number}("CO2/electricity"=>1000),optimizer_config=optimizer_config)
       exact_scenarios["$state-$years-slack"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="slack",lost_load_cost=Dict{String,Number}("electricity"=>1e6), lost_emission_cost=Dict{String,Number}("CO2"=>700),optimizer_config=optimizer_config)
       exact_scenarios["$state-$years-ex"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="ex",infrastructure=Dict{String,Array}("existing"=>["all"]),optimizer_config=optimizer_config)
       exact_scenarios["$state-$years-simple"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="simple storage",storage_type="simple",conversion=true,optimizer_config=optimizer_config)
       exact_scenarios["$state-$years-seasonal"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="seasonal storage",storage_type="seasonal",conversion=true,optimizer_config=optimizer_config)
       design_result=run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="des&op",optimizer_config=optimizer_config)
       exact_scenarios["$state-$years-des&op"] = run_opt(ts_full_data.clust_data,cep_data,design_result.config,get_cep_design_variables(design_result),optimizer;lost_load_cost=Dict{String,Number}("electricity"=>1e6), lost_emission_cost=Dict{String,Number}("CO2"=>700))
end
for (state, years) in [["GER_18", [2016]],["CA_14", [2016]]]
       # laod data
       ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
       cep_data = load_cep_data_provided(state)
       ## CLUSTERING ##
       ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
       ## OPTIMIZATION ##
       optimizer=Clp.Optimizer
       optimizer_config=Dict{Symbol,Any}(:LogLevel => 0)
       exact_scenarios["$state-$years-trans"]=run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="trans",transmission=true,optimizer_config=optimizer_config)
end

@save normpath(joinpath(dirname(@__FILE__),"cep_exact_data.jld2")) exact_scenarios


text:
module BitTwiddlingConvenienceFunctions

using Static

intlog2(N::I) where {I <: Integer} = (8sizeof(I) - one(I) - leading_zeros(N)) % I
intlog2(::Type{T}) where {T} = intlog2(sizeof(T))
nextpow2(W::T) where {T<:Base.BitInteger} = (one(T) << (T(8sizeof(W)) - leading_zeros((W - one(T)))))
prevpow2(W::T) where {T<:Base.BitInteger} = (one(T) << (((T(8sizeof(W))) - one(T)) - leading_zeros(W)))
@generated nextpow2(::StaticInt{N}) where {N} = Expr(:call, Expr(:curly, :StaticInt, nextpow2(N)))
@generated prevpow2(::StaticInt{N}) where {N} = Expr(:call, Expr(:curly, :StaticInt, prevpow2(N)))
@generated intlog2(::StaticInt{N}) where {N} = Expr(:call, Expr(:curly, :StaticInt, intlog2(N)))

end

ext: jl


text:
"""Functions for computing radio models."""

abstract type AbstractRadioModel{T<:Real} <: AbstractForegroundModel end

"""
    Radio_Sehgal2009{T}(model parameters...)

Define CIB model parameters. Defaults are from Viero et al. 2013.

```@example
model = CIBModel{Float32}(a_0=0.4)
```
"""
Base.@kwdef struct Radio_Sehgal2009{T<:Real} <: AbstractRadioModel{T}
    nside::Int64    = 4096
    min_redshift::T = 0.0
    max_redshift::T = 5.0
    min_mass::T     = 1e13
    box_size::T     = 40000

    # these coefficients are shared for type I and II
    a_0::T   = -1.0
    a_1_dist::Distributions.Uniform{T} = Distributions.Uniform(T(-0.12), T(0.07))
    a_2_dist::Distributions.Uniform{T} = Distributions.Uniform(T(-0.34), T(0.99))
    a_3_dist::Distributions.Uniform{T} = Distributions.Uniform(T(-0.75), T(-0.23))

    I_R_int::T = 10^(-2.6)
    I_γ::T = 6.0
    I_N_0::T = 1.0
    I_M_0::T = 4e13
    I_α::T = 0.1
    I_L_b::T = 10^(24.0)
    I_m::T = -1.55
    I_n::T = 0.0
    I_δ::T = 3.0
    I_z_p::T = 0.8

    II_R_int::T = 10^(-2.8)
    II_γ::T = 8.0
    II_N_0::T = 0.015
    II_M_0::T = 3e15
    II_α::T = 0.1
    II_L_b::T = 10^(27.5)
    II_m::T = -1.6
    II_n::T = -0.65
    II_z_p::T = 1.3
    II_σ_l::T = 0.4
    II_σ_h::T = 0.73

    I_L_min::T = 2.5e23
    II_L_min::T = 4e23

    # physical constants
    phys_h::T     = 6.62606957e-27      # erg.s
    phys_c::T     = 3e+10               # cm/s
    phys_k::T     = 1.3806488e-16       # erg/K
    phys_Msun::T  = 2e33                # g
    phys_Mpc::T   = 3.086e24            # cm
end

pdf_norm1(x, μ, σ) = exp( -(x-μ)^2 / (2 * σ^2) );

function FR_I_redshift_evolution(z, model)
    if z < model.I_z_p
        return (1+z)^model.I_δ
    else
        return (1+model.I_z_p)^model.I_δ
    end
end

function FR_II_redshift_evolution(z, model)
    norm = pdf_norm1(0, model.II_z_p, model.II_σ_l)
    if z < model.II_z_p
        return pdf_norm1(z, model.II_z_p, model.II_σ_l) / norm
    else
        return pdf_norm1(z, model.II_z_p, model.II_σ_h) / norm
    end
end

"""
Populate halos with radio sources according to the HOD in Sehgal et al. 2009.

The optional rng parameter provides an array of random number generators, one
for each thread.
"""
function hod_sehgal(
        halo_mass, redshift, model::Radio_Sehgal2009{T}) where T

    N_halos = size(halo_mass,1)
    nsources_I = Array{Int32}(undef, N_halos)
    nsources_II = Array{Int32}(undef, N_halos)

    # compute poisson mean, then draw it
    Threads.@threads for i = 1:N_halos
        I_HON = model.I_N_0 * (halo_mass[i] / model.I_M_0)^model.I_α
        I_HON *= FR_I_redshift_evolution(redshift[i], model)
        nsources_I[i] = rand(Distributions.Poisson(Float64.(I_HON)))
        II_HON = model.II_N_0 * (halo_mass[i] / model.II_M_0)^model.II_α
        II_HON *= FR_II_redshift_evolution(redshift[i], model)
        nsources_II[i] = rand(Distributions.Poisson(Float64.(II_HON)))
    end

    return nsources_I, nsources_II
end

function sehgal_LFn0(rand_unit::T, m, Lb, Lmin) where T <: Real
    crossover = log( Lb/Lmin )
    x = rand_unit * (-(1/m)+log(Lb/Lmin))
    if x < crossover
        return exp(x) * Lmin
    else
        return Lb * (1+m * x + m * log(Lmin/Lb))^(1/m)
    end
end

function sehgal_LF(rand_unit::T, m, n, Lb, Lmin) where T
    crossover = (T(1) - Lb^(-n) * Lmin^n)/n
    x = (rand_unit * (m-Lb^(-n) * Lmin^n * m-n)/(m * n))
    if x < crossover
        return (Lmin^n + Lb^n * n * x)^(1/n)
    else
        return Lb * (1 + (m * (-1 + (Lmin/Lb)^n + n * x))/n)^(1/m)
    end
end

"""
Fills the result array with draws from the luminosity function.
"""
function sehgal_LF!(result::Array{T,1}, m, n, Lb::T, Lmin,
    rand_array::Array{T,1}) where T

    if n ≈ 0.0
        Threads.@threads for i = 1:size(result,1)
            result[i] = sehgal_LFn0(
                rand_array[i],
                Float64(m), Float64(Lb), Float64(Lmin))
        end
    else
        Threads.@threads for i = 1:size(result,1)
            result[i] = sehgal_LF(
                rand_array[i],
                Float64(m), Float64(n), Float64(Lb), Float64(Lmin))
        end
    end
end

B(cosθ, β) = ( (1-β*cosθ)^(-2) + (1+β*cosθ)^(-2) ) / 2

function get_core_lobe_lum(L_beam::T, ν_Hz, R_int, γ, a_0, a_1, a_2, a_3, cosθ) where T
    β = sqrt(1-γ^(-2))
    R_obs = R_int * B(cosθ, β)
    L_int = L_beam * (1 + R_int) / (1 + R_obs)
    L_l_int = L_int / (1 + R_int)
    L_c_beam = R_obs * L_l_int

    # now put in the frequency dependence. a_0 = 0.0
    lν = log10(ν_Hz / T(1e9) ) # normed to 151 MHz
    f_core =  T(10) ^ ( a_1 * lν + a_2 * lν^2 + a_3 * lν^3 ) # a_0 +

    lν_norm = T(log10(151e6 / 1e9 )) # normed to 151 MHz
    norm_core =  10 ^ ( a_1 * lν_norm + a_2 * lν_norm^2 + a_3 * lν_norm^3 )

    f_lobe = (ν_Hz / T(151e6) ).^T(-0.8)
    return T(L_c_beam * f_core * 0.3f0 / norm_core), T(L_l_int * f_lobe)
end

function draw_coeff!(a_coeff, model::Radio_Sehgal2009)
    a1_task = Threads.@spawn rand!(model.a_1_dist, @view a_coeff[1,:])
    a2_task = Threads.@spawn rand!(model.a_2_dist, @view a_coeff[2,:])
    a3_task = Threads.@spawn rand!(model.a_3_dist, @view a_coeff[3,:])

    wait(a1_task)
    wait(a2_task)
    wait(a3_task)
end

"""
Produce a source catalog from a model and halo catalog.
"""
function generate_sources(
        # model parameters
        model::AbstractRadioModel{T}, cosmo::Cosmology.FlatLCDM{T},
        # halo arrays
        halo_pos_inp::AbstractArray{TH,2}, halo_mass_inp::AbstractArray{TH,1};
        verbose=true) where {T, TH}

    # make sure halo inputs are the model type
    halo_pos = convert(Array{T,2}, halo_pos_inp)
    halo_mass = convert(Array{T,1}, halo_mass_inp)

    res = Resolution(model.nside)

    verbose && println("Culling halos below mass $(model.min_mass).")
    mass_cut = halo_mass .> model.min_mass
    halo_pos = halo_pos[:, mass_cut]
    halo_mass = halo_mass[mass_cut]
    N_halos = size(halo_mass, 1)

    verbose && println("Allocating for $(N_halos) halos.")
    dist, redshift, hp_ind = get_basic_halo_properties(halo_pos, model, cosmo, res)
    θ, ϕ = get_angles(halo_pos)

    verbose && println("Populating HOD.")
    nsources_I, nsources_II = hod_sehgal(
        halo_mass, redshift, model)

     # total number of sources
    total_n_I, total_n_II = sum(nsources_I), sum(nsources_II)

    # the last element of cumsat_I is the total number of type I sources
    verbose && println("Drawing spectral coefficients.")
    a_coeff_I = Array{T, 2}(undef, 3, total_n_I)
    a_coeff_II = Array{T, 2}(undef, 3, total_n_II)
    draw_coeff!(a_coeff_I, model)
    draw_coeff!(a_coeff_II, model)

    verbose && println("Drawing from luminosity function.")
    # generate random numbers for LF
    rand_buffer_I = Array{T, 1}(undef, total_n_I)
    rand_buffer_II = Array{T, 1}(undef, total_n_II)
    threaded_rand!(rand_buffer_I)
    threaded_rand!(rand_buffer_II)

    # draw luminosities for both populations
    L_I_151 = Array{T, 1}(undef, total_n_I)
    L_II_151 = Array{T, 1}(undef, total_n_II)
    sehgal_LF!(L_I_151, model.I_m, model.I_n, model.I_L_b,
        model.I_L_min, rand_buffer_I)
    sehgal_LF!(L_II_151, model.II_m, model.II_n, model.II_L_b,
        model.II_L_min, rand_buffer_II)

    verbose && println("Drawing for impact parameter.")
    # reuse rand buffers for the impact parameters
    threaded_rand!(rand_buffer_I)
    threaded_rand!(rand_buffer_II)
    cosθ_I = rand_buffer_I
    cosθ_II = rand_buffer_II

    return (
        redshift=redshift, halo_mass=halo_mass, halo_pos=halo_pos,
        halo_hp_ind=hp_ind, θ=θ, ϕ=ϕ,
        dist=dist, nsources_I=nsources_I, nsources_II=nsources_II,
        a_coeff_I=a_coeff_I, a_coeff_II=a_coeff_II,
        L_I_151=L_I_151, L_II_151=L_II_151,
        cosθ_I=cosθ_I, cosθ_II=cosθ_II,
        total_n_I=total_n_I, total_n_II=total_n_II
    )
end

function paint!(result_map::HealpixMap{T_map,RingOrder},
                nu_obs::T, model::AbstractRadioModel, sources;
                return_fluxes=false) where {T_map, T}

    pixel_array = result_map.pixels
    fill!(pixel_array, zero(T))  # prepare the frequency map

    flux_to_Jy = ustrip(u"Jy", 1u"W/Hz*Mpc^-2")

    source_offset_I = generate_subhalo_offsets(sources.nsources_I)
    source_offset_II = generate_subhalo_offsets(sources.nsources_II)
    N_halo = size(sources.halo_mass, 1)
    total_n_sat_I = size(sources.L_I_151, 1)
    total_n_sat_II = size(sources.L_II_151, 1)

    flux_I = Array{T, 1}(undef, total_n_sat_I)
    flux_II = Array{T, 1}(undef, total_n_sat_II)
    redshift_I = Array{T, 1}(undef, total_n_sat_I)
    redshift_II = Array{T, 1}(undef, total_n_sat_II)
    θ_I = Array{T, 1}(undef, total_n_sat_I)
    θ_II = Array{T, 1}(undef, total_n_sat_II)
    ϕ_I = Array{T, 1}(undef, total_n_sat_I)
    ϕ_II = Array{T, 1}(undef, total_n_sat_II)

    Threads.@threads for i_halo in 1:N_halo
        nu = (1 + sources.redshift[i_halo]) * nu_obs
        hp_ind = sources.halo_hp_ind[i_halo]

        # Generate FR I fluxes
        for source_index in 1:sources.nsources_I[i_halo]
            # index of satellite in satellite arrays
            i_sat = source_offset_I[i_halo] + source_index
            L_c, L_l = get_core_lobe_lum(
                sources.L_I_151[i_sat], nu, model.I_R_int, model.I_γ,
                model.a_0,
                sources.a_coeff_I[1,i_sat], sources.a_coeff_I[2,i_sat],
                sources.a_coeff_I[3,i_sat], sources.cosθ_I[i_sat])

            flux_I[i_sat] = l2f(L_c+L_l, sources.dist[i_halo],
                sources.redshift[i_halo]) * flux_to_Jy
            redshift_I[i_sat] = sources.redshift[i_halo]
            θ_I[i_sat] = sources.θ[i_halo]
            ϕ_I[i_sat] = sources.ϕ[i_halo]
            pixel_array[hp_ind] += flux_I[i_sat]
        end

        # Generate FR II fluxes
        for source_index in 1:sources.nsources_II[i_halo]
            # index of satellite in satellite arrays
            i_sat = source_offset_II[i_halo] + source_index
            L_c, L_l = get_core_lobe_lum(
                sources.L_II_151[i_sat], nu, model.II_R_int, model.II_γ,
                model.a_0,
                sources.a_coeff_II[1,i_sat], sources.a_coeff_II[2,i_sat],
                sources.a_coeff_II[3,i_sat], sources.cosθ_II[i_sat])

            flux_II[i_sat] = l2f(L_c+L_l, sources.dist[i_halo],
                sources.redshift[i_halo]) * flux_to_Jy
            redshift_II[i_sat] = sources.redshift[i_halo]
            θ_II[i_sat] = sources.θ[i_halo]
            ϕ_II[i_sat] = sources.ϕ[i_halo]
            pixel_array[hp_ind] += flux_II[i_sat]
        end
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(result_map.resolution.nside)
    pixel_array .*= per_pixel_steradian

    # maps are in Jansky per steradian, fluxes are in Jansky
    if return_fluxes
        return flux_I, redshift_I, θ_I, ϕ_I, flux_II, redshift_II, θ_II, ϕ_II
    end

end

max_stars_repo_name: xzackli/XGPaint

max_issues_repo_name: xzackli/XGPaint

max_forks_repo_name: xzackli/XGPaint


text:
"""
    Reprint(fn) -- apply the lambda function when printed
"""
mutable struct Reprint
    content::Function
end

Base.print(io::IO, r::Reprint) = r.content(io)

"""
    Render(data) - printed object shows its text/html
"""
struct Render{T}
    content::T
end

Base.print(io::IO, r::Render) =
    show(io, MIME"text/html"(), r.content)

"""
    Bypass(data) - printed object passes though EscapeProxy unescaped
"""

mutable struct Bypass{T}
    content::T
end

Base.print(io::IO, x::Bypass) = print(io, x.content)

"""
    EscapeProxy(io) - wrap an `io` to perform HTML escaping

This is a transparent proxy that performs HTML escaping so that objects
that are printed are properly converted into valid HTML values. As a
special case, objects wrapped with `Bypass` are not escaped, and
bypass the proxy.

# Examples
```julia-repl
julia> ep = EscapeProxy(stdout);
julia> print(ep, "A&B")
A&B
julia> print(ep, Bypass("<tag/>"))
<tag/>
```
"""
struct EscapeProxy{T<:IO} <: IO
    io::T
end

EscapeProxy(io::EscapeProxy) = io

Base.print(ep::EscapeProxy, h::Reprint) = h.content(ep)
Base.print(ep::EscapeProxy, w::Render) =
    show(ep.io, MIME"text/html"(), w.content)
Base.print(ep::EscapeProxy, x::Bypass) = print(ep.io, x.content)

function Base.write(ep::EscapeProxy, octet::UInt8)
    if octet == Int('&')
        write(ep.io, "&")
    elseif octet == Int('<')
        write(ep.io, "<")
    elseif octet == Int('"')
        write(ep.io, """)
    elseif octet == Int('\'')
        write(ep.io, "'")
    else
        write(ep.io, octet)
    end
end

function Base.unsafe_write(ep::EscapeProxy, input::Ptr{UInt8}, nbytes::UInt)
    written = 0
    last = cursor = input
    final = input + nbytes
    while cursor < final
        ch = unsafe_load(cursor)
        if ch == Int('&')
            written += unsafe_write(ep.io, last, cursor - last)
            written += unsafe_write(ep.io, pointer("&"), 5)
            cursor += 1
            last = cursor
            continue
        end
        if ch == Int('<')
            written += unsafe_write(ep.io, last, cursor - last)
            written += unsafe_write(ep.io, pointer("<"), 4)
            cursor += 1
            last = cursor
            continue
        end
        if ch == Int('\'')
            written += unsafe_write(ep.io, last, cursor - last)
            written += unsafe_write(ep.io, pointer("'"), 6)
            cursor += 1
            last = cursor
            continue
        end
        if ch == Int('"')
            written += unsafe_write(ep.io, last, cursor - last)
            written += unsafe_write(ep.io, pointer("""), 6)
            cursor += 1
            last = cursor
            continue
        end
        cursor += 1
    end
    if last < final
        written += unsafe_write(ep.io, last, final - last)
    end
    return written
end


text:
## General purpose functions

#User is expected to pass functions
#x_k(p0, alpha, mhatk)  : provides a solution with data mhatk
#c_ik(x)
#These are often collected into arrays of functions
# xs			: [x_k for k = 1:K ]
# c_k  			: [c_ik for i = 1:d ]
# cs			: (d, K) Matrix with elements c_ik

#General purpose functions
###
function shrink(phat_k, p0, alpha, Nhat_k)
	if Nhat_k <= 0
		return p0
	end
	@. alpha * p0 / (Nhat_k + alpha) + Nhat_k * phat_k / (Nhat_k + alpha)
end

#need to be careful because some problem smay have zero data
function get_GM_anchor(mhats)
	Nhats = vec(sum(mhats, dims=1))
	non_zero_indx = Nhats .> 0
	vec(mean(mhats[:, non_zero_indx] ./ Nhats[non_zero_indx]', dims=2))
end

function sim_path(p_k, N::Int)
	@assert N > 0 "Something weird"
    mhats_k = zeros(length(p_k))
    for j in rand(Categorical(p_k), N)
        mhats_k[j] += 1
    end
    mhats_k
end

#ps is a matrix of pk
function sim_path(ps, Nhats)
    K = size(ps, 2)
    out = zeros(size(ps))
    for k = 1:K
        out[:, k] = sim_path(ps[:, k], Nhats[k])
    end
    out
end

#empirical bayes moment-matched estimates.
#returns p0, alpha0
#VG this method does not seem to consistenty yield alpha in [0, 1]
function eb_mm_estimates(mhats)
	Nhats = sum(mhats, dims=1)
	p0 = mean(mhats ./ Nhats, dims=2)
	K = size(mhats, 2)

	C0 = sum(mhats.^2)/K
	sq_norm_p0 = norm(p0)^2
	Nbar = mean(Nhats)
	var_n = mean( @.(Nhats^2 - Nhats) )

	alpha0 = C0 - sq_norm_p0 * var_n - Nbar
	@assert alpha0 > 0 "Moment matching fails"
	alpha0 /= (1 - sq_norm_p0) / var_n
	alpha0 = 1 / alpha0 - 1

	return p0, alpha0
end

#solves an approximation to the Apriori MSE equation for alpha
#for large K this is like AlphaOR for MSE
#for now, super lazy and just search a grid.
#returns AlphaEstimate, minimizingAlphaIndex
function mse_estimates(mhats, supps, p0, alpha_grid)
	K = size(mhats, 2)
	Nhats = sum(mhats, dims=1)
	phats = mhats ./ Nhats
	mu0s = vec(p0' * supps)

	#summary statistics
	muhats = Vector{Float64}(undef, K)
	sigmas = Vector{Float64}(undef, K)
	mse0 = Vector{Float64}(undef, K)
	for k = 1:K
		muhats[k] = dot(phats[:, k], supps[:, k])
		sigmas[k] = dot(phats[:, k], supps[:, k].^2) - muhats[k]^2
		mse0[k] = dot(phats[:, k], @.((supps[:, k] - mu0s[k])^2) )
	end

	#worker function
	function mse(a)
		out = 0.
		for k = 1:K
			out += (Nhats[k] - a^2) / (Nhats[k] + a)^2 * sigmas[k]
			out += ( a / (a + Nhats[k]) )^2 * mse0[k]
		end
		out/K
	end
	out = map(mse, alpha_grid)
	alpha_grid[argmin(out)], argmin(out)
end

function z_k(x_k, c_k, p0, alpha, mhat_k, ps_k, lam_k, lamavg)
	x = x_k(p0, alpha, mhat_k)
	cs = map(c_ik -> c_ik(x), c_k)
	lam_k/lamavg * dot(ps_k, cs)
end

#actually returns zLOO_k * N * lamavg
function zLOO_k_unsc(x_k, c_k, p0, alpha, mhat_k)
	out = 0.
	mhatloo = mhat_k[:] #copy data
	#only compute for terms where mhat_k[i] > 0
	#compute a base solution to resuse as necessary
	x_base = x_k(p0, alpha, mhat_k)
	for i = 1:length(mhat_k)
		#correct previous toggle
		if i > 1 && mhat_k[i - 1] > 0
			mhatloo[i - 1] += 1
		end

		if mhat_k[i] > 0
			mhatloo[i] -= 1
		    x = x_k(p0, alpha, mhatloo)
		else
			x = x_base
		end
		out += mhat_k[i] * c_k[i](x)
	end
	out
end

#average true performance
function zbar(xs, cs, p0, alpha, mhats, ps, lams)
	lamavg = mean(lams)
	K = size(cs, 2)
	out = 0.
	#VG maybe change this to a more numerically stable way to do avg?
	for k = 1:K
		out += z_k(xs[k], cs[:, k], p0, alpha, mhats[:, k], ps[:, k], lams[k], lamavg)
	end
	out/K
end

#full-info for scaling/comparison
#uses the fact that data scale doesn't matter.
zstar(xs, cs, ps, lams) = zbar(xs, cs, ps[:, 1], 0., ps, ps, lams)

function zLOObar_unsc(xs, cs, p0, alpha, mhats)
	K = size(cs, 2)
	out = 0.
	for k = 1:K
		out += zLOO_k_unsc(xs[k], cs[:, k], p0, alpha, mhats[:, k])
	end
	out /K
end

#returns alphaOR, minimizingAlphaIndex, curveInAlpha
function oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
	alphaOR = 0.
	jstar = -1
	best_val = Inf
	out = zeros(length(alpha_grid))
	for (j, alpha) in enumerate(alpha_grid)
		out[j] = zbar(xs, cs, p0, alpha, mhats, ps, lams)
		if out[j] < best_val
			jstar = j
			alphaOR = alpha
			best_val = out[j]
		end
	end
	return alphaOR, jstar, out
end

#returns alphaLOO, jstar, and cuveInAlpha_unsc
function loo_alpha(xs, cs, mhats, p0, alpha_grid)
	alphaLOO = 0.
	jstar = -1
	best_val = Inf
	out = zeros(length(alpha_grid))
	for (j, alpha) in enumerate(alpha_grid)
		out[j] = zLOObar_unsc(xs, cs, p0, alpha, mhats)
		if out[j] < best_val
			jstar = j
			alphaLOO = alpha
			best_val = out[j]
		end
	end
	return alphaLOO, jstar, out
end

###
#### Deprecated
#Convenience functions to generate a sequence of newsvendor problems
#  with common support, potentially different service levels
#  assumes supp is ordered vector of Reals
# returns cs, xs
# function genNewsvendors(supp, ss, K)
# 	@assert issorted(supp) "Supp must be an ordered vector"

# 	#Generic computation of the sth quantile
# 	function x_k(p0, alpha, mhat_k, s)
# 	    const Nhat_k = sum(mhat_k)
# 	    palpha = JS.shrink(mhat_k./Nhat_k, p0, alpha, Nhat_k)
# 	    indx = quantile(Categorical(palpha), s)
# 	    supp[indx]
# 	end

# 	function c_ik(i, x, s)
# 	    supp[i] > x ? s/(1-s) * (supp[i] - x) : (x - supp[i])
# 	end

# 	xs  = [(p0, alpha, mhat_k)-> x_k(p0, alpha, mhat_k, ss[k]) for k = 1:K]
# 	cs = [x->c_ik(i, x, ss[k]) for i = 1:length(supp), k = 1:K]
# 	cs, xs
# end

###
# Ideally should have a broadcast implementation that combines
# previous two...
function genNewsvendorsDiffSupp(supps, s, K)
	#Generic computation of the sth quantile
	function x_k(p0, alpha, mhat_k, k, s)
	    Nhat_k = sum(mhat_k)
	    palpha = JS.shrink(mhat_k./Nhat_k, p0, alpha, Nhat_k)
	    if !isapprox(sum(palpha), 1) || sum(mhat_k .< 0) > 0
	    	println("Nhat_k:\t", Nhat_k)
	    	println("alpha:\t", alpha)
	    	println("mhat_k:\n", mhat_k)
	    	println("phat_k:\n", palpha)
	    end

	    indx = quantile(Categorical(palpha), s)
	    supps[indx, k]
	end

	function c_ik(i, k, x, s)
	    supps[i, k] > x ? s/(1-s) * (supps[i, k] - x) : (x - supps[i, k])
	end

	xs  = [(p0, alpha, mhat_k)-> x_k(p0, alpha, mhat_k, k, s) for k = 1:K]
	cs = [x->c_ik(i, k, x, s) for i = 1:size(supps, 1), k = 1:K]
	cs, xs
end

#### Helper function for binning data
#dat_vec is vector of actual realizations
#computes a grid of size d
#returns
	#d vector of lefthand binpoints
	#d vector of counts
	#raw histogram object for safety
#assumes NA's already filtered from calcualtion
function bin_data(dat_vec, d; TOL = .001)
	#Explicitly make bins to prevent histogram from rounding
	l = minimum(dat_vec)
	u = maximum(dat_vec)
	bin_width = (u - l)/d
    bin_edges = collect(range(l, stop=u, length=d + 1))
    bin_edges[end] += TOL * bin_width  #inflate the last one a little...
    bin_edges[1] -= TOL * bin_width
    dat_hist = fit(Histogram, dat_vec, bin_edges, closed =:left)
    bin_edges[1:d], dat_hist.weights, dat_hist
end

max_stars_repo_name: CausalML/DataPooling

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
"""
	msfm()

	Calculate the shortest distance from a list of
	points to all other pixels in an image, using the
	Multistencil Fast Marching Method (MSFM).

	This method gives more accurate
	distances by using second order derivatives and cross neighbours.

	T=msfm2d(F, SourcePoints, UseSecond, UseCross)

	# inputs,
	  F: The speed image. The speed function must always be larger
				than zero (min value 1e-8), otherwise some regions will
				never be reached because the time will go to infinity.
	  SourcePoints : A list of starting points [2 x N] (distance zero)
	  UseSecond : Boolean Set to true if not only first but also second
	               order derivatives are used (default)
	  UseCross : Boolean Set to true if also cross neighbours
	               are used (default)
	# outputs
	  T : Image with distance from SourcePoints to all pixels
"""
function msfm(
    speedimage::AbstractArray{T,2},
    SourcePoints::AbstractArray{T},
    usesecond::Bool=true,
    usecross::Bool=true,
    Ed=false::Bool,
) where {T}
    distanceimage = zeros(T, size(speedimage))

    Frozen = falses(size(speedimage)) # Pixels which are processed and have a final distance are frozen

    msfm!(distanceimage, Frozen, speedimage, SourcePoints, usesecond, usecross, Ed)
    return distanceimage
end

function msfm!(
    distanceimage,
    Frozen,
    speedimage::AbstractArray{T,2},
    SourcePoints::AbstractArray{T},
    usesecond::Bool=true,
    usecross::Bool=true,
    Ed=false::Bool,
) where {T}
    fill!(distanceimage, T(-Inf))
    # Augmented Fast Marching (For skeletonize)

    # Euclidian distance image
    if Ed
        Y = zeros(T, size(speedimage))
    end

    # Store neighbours of the (segmented) region
    neg_free = length(speedimage)
    neg_pos = 0
    neg_list_1 = zeros(T, neg_free)
    neg_list_x = fill(1, neg_free)
    neg_list_y = fill(1, neg_free)
    if Ed
        neg_list_4 = zeros(T, neg_free)
    end

    # (There are 3 pixel classes:
    #   - frozen (processed)
    #   - narrow band (boundary) (in list to check for the next pixel with smallest distance)
    #   - far (not yet used)

    # Neighbours
    ne = [
        -1 0
        1 0
        0 -1
        0 1
    ]
    Tt = T(0)
    Ty = T(0)

    # set all starting points to distance zero and frozen
    for z in 1:size(SourcePoints, 2)
        # starting point
        x = round(Int, SourcePoints[1, z])
        y = round(Int, SourcePoints[2, z])
        # Set starting point to frozen and distance to zero
        Frozen[x, y] = true
        distanceimage[x, y] = 0
    end

    # Add all neighbours of the starting points to narrow list
    @inbounds for z in 1:size(SourcePoints, 2)
        # starting point
        x = round(Int, SourcePoints[1, z])
        y = round(Int, SourcePoints[2, z])

        for k in 1:4
            # Location of neighbour
            i = x + ne[k, 1]
            j = y + ne[k, 2]
            # Check if current neighbour is not yet frozen and inside the
            # picture
            if (
                (i > 0) &&
                (j > 0) &&
                (i <= size(speedimage, 1)) &&
                (j <= size(speedimage, 2)) &&
                (~Frozen[i, j])
            )
                Tt = 1 / max(speedimage[i, j], eps(T))
                Ty = T(1)
                # Update distance in neighbour list or add to neighbour list
                if distanceimage[i, j] > 0
                    if neg_list_1[distanceimage[i, j]] > Tt
                        neg_list_1[distanceimage[i, j]] = Tt
                    end
                    if Ed
                        neg_list_4[distanceimage[i, j]] = min(
                            Ty, neg_list_4[distanceimage[i, j]]
                        )
                    end
                else
                    neg_pos = neg_pos + 1
                    neg_list_1[neg_pos] = Tt
                    neg_list_x[neg_pos] = i
                    neg_list_y[neg_pos] = j
                    if Ed
                        neg_list_4[neg_pos] = Ty
                    end
                    distanceimage[i, j] = neg_pos
                end
            end
        end
    end

    # Reuse variables:
    Tpatch = zeros(T, 5, 5)
    Order = zeros(Int8, 1, 4)
    Tm = zeros(T, 1, 4)
    Tm2 = zeros(T, 1, 4)
    Coeff = zeros(T, 3)
    TT = zeros(T, 2)
    TT2 = zeros(T, 2)

    # Loop through all pixels of the image
    for itt in 1:length(speedimage)
        # Get the pixel from narrow list (boundary list) with smallest
        # distance value and set it to current pixel location

        if neg_pos == 0
            break
        end
        (t, index) = findmin(view(neg_list_1, 1:neg_pos))
        x = neg_list_x[index]
        y = neg_list_y[index]
        Frozen[x, y] = true
        distanceimage[x, y] = neg_list_1[index]

        if Ed
            Y[x, y] = neg_list_4[index]
        end

        # Remove min value by replacing it with the last value in the array
        if index < neg_pos
            neg_list_1[index] = neg_list_1[neg_pos]
            neg_list_x[index] = neg_list_x[neg_pos]
            neg_list_y[index] = neg_list_y[neg_pos]
            if Ed
                neg_list_4[index] = neg_list_4[neg_pos]
            end
            x2 = neg_list_x[index]
            y2 = neg_list_y[index]
            distanceimage[x2, y2] = index
        end
        neg_pos = neg_pos - 1
        # Loop through all 4 neighbours of current pixel
        for k in 1:4
            # Location of neighbour
            i = x + ne[k, 1]
            j = y + ne[k, 2]

            # Check if current neighbour is not yet frozen and inside the
            # picture
            if (i > 0) &&
               (j > 0) &&
               (i <= size(speedimage, 1)) &&
               (j <= size(speedimage, 2)) &&
               (~Frozen[i, j])
                Tt = calculatedistance!(
                    distanceimage,
                    Tpatch,
                    Order,
                    Coeff,
                    Tm,
                    Tm2,
                    TT,
                    TT2,
                    speedimage[i, j],
                    size(speedimage),
                    i,
                    j,
                    usesecond,
                    usecross,
                    Frozen,
                )
                if Ed
                    Ty = calculatedistance!(
                        Y,
                        Tpatch,
                        Order,
                        Coeff,
                        Tm,
                        Tm2,
                        TT,
                        TT2,
                        speedimage[i, j],
                        size(speedimage),
                        i,
                        j,
                        usesecond,
                        usecross,
                        Frozen,
                    )
                end

                # Update distance in neighbour list or add to neighbour list
                if distanceimage[i, j] > 0
                    neg_list_1[round(Int, distanceimage[i, j])] = min(
                        Tt, neg_list_1[round(Int, distanceimage[i, j])]
                    )
                    if Ed
                        neg_list_4[round(Int, distanceimage[i, j])] = min(
                            Ty, neg_list_4[round(Int, distanceimage[i, j])]
                        )
                    end
                else
                    neg_pos = neg_pos + 1
                    neg_list_1[neg_pos] = Tt
                    neg_list_x[neg_pos] = i
                    neg_list_y[neg_pos] = j
                    if Ed
                        neg_list_4[neg_pos] = Ty
                    end
                    distanceimage[i, j] = neg_pos
                end
            end
        end
    end
    return distanceimage
end # End function msfm

function calculatedistance(
    TI::AbstractArray{T}, Fij, sizeF, i, j, usesecond, usecross, Frozen
) where {T}
    Tpatch = zeros(T, 5, 5)
    Order = zeros(Int8, 1, 4)
    Tm = zeros(T, 1, 4)
    Tm2 = zeros(T, 1, 4)
    Coeff = zeros(T, 3)
    TT = zeros(T, 2)
    TT2 = zeros(T, 2)
    return calculatedistance!(
        TI,
        Tpatch,
        Order,
        Coeff,
        Tm,
        Tm2,
        TT,
        TT2,
        Fij,
        sizeF,
        i,
        j,
        usesecond,
        usecross,
        Frozen,
    )
end

function calculatedistance!(
    TI::AbstractArray{T},
    Tpatch,
    Order,
    Coeff,
    Tm,
    Tm2,
    TT,
    TT2,
    Fij,
    sizeF,
    i,
    j,
    usesecond,
    usecross,
    Frozen,
) where {T}
    fill!(Tpatch, T(Inf))
    @inbounds for nx in -2:2
        for ny in -2:2
            i_n = i + nx
            j_n = j + ny
            if i_n > 0 && j_n > 0 && i_n <= sizeF[1] && j_n <= sizeF[2] && Frozen[i_n, j_n]
                Tpatch[nx + 3, ny + 3] = TI[i_n, j_n]
            end
        end
    end

    # The values in order is 0 if no neighbours in that direction
    # 1 if 1e order derivatives is used and 2 if second order
    # derivatives are used
    fill!(Order, 0)

    # Make 1e order derivatives in x and y direction
    fill!(Tm, 0)
    Tm[1] = min(Tpatch[2, 3], Tpatch[4, 3])
    if isfinite(Tm[1])
        Order[1] = 1
    end

    Tm[2] = min(Tpatch[3, 2], Tpatch[3, 4])
    if isfinite(Tm[2])
        Order[2] = 1
    end

    # Make 1e order derivatives in cross directions
    if usecross
        Tm[3] = min(Tpatch[2, 2], Tpatch[4, 4])
        if isfinite(Tm[3])
            Order[3] = 1
        end
        Tm[4] = min(Tpatch[2, 4], Tpatch[4, 2])
        if isfinite(Tm[4])
            Order[4] = 1
        end
    end

    # Make 2e order derivatives
    if usesecond
        fill!(Tm2, 0)
        # pixels with a pixeldistance 2 from the center must be
        # lower in value otherwise use other side or first order
        ch1 = (Tpatch[1, 3] < Tpatch[2, 3]) && isfinite(Tpatch[2, 3])
        ch2 = (Tpatch[5, 3] < Tpatch[4, 3]) && isfinite(Tpatch[4, 3])

        if ch1 && ch2
            Tm2[1] = min(
                (4 * Tpatch[2, 3] - Tpatch[1, 3]) / 3, (4 * Tpatch[4, 3] - Tpatch[5, 3]) / 3
            )
            Order[1] = 2
        elseif ch1
            Tm2[1] = (4 * Tpatch[2, 3] - Tpatch[1, 3]) / 3
            Order[1] = 2
        elseif ch2
            Tm2[1] = (4 * Tpatch[4, 3] - Tpatch[5, 3]) / 3
            Order[1] = 2
        end

        ch1 = (Tpatch[3, 1] < Tpatch[3, 2]) && isfinite(Tpatch[3, 2])
        ch2 = (Tpatch[3, 5] < Tpatch[3, 4]) && isfinite(Tpatch[3, 4])

        if ch1 && ch2
            Tm2[2] = min(
                (4 * Tpatch[3, 2] - Tpatch[3, 1]) / 3, (4 * Tpatch[3, 4] - Tpatch[3, 5]) / 3
            )
            Order[2] = 2
        elseif ch1
            Tm2[2] = (4 * Tpatch[3, 2] - Tpatch[3, 1]) / 3
            Order[2] = 2
        elseif ch2
            Tm2[2] = (4 * Tpatch[3, 4] - Tpatch[3, 5]) / 3
            Order[2] = 2
        end

        if usecross
            ch1 = (Tpatch[1, 1] < Tpatch[2, 2]) && isfinite(Tpatch[2, 2])
            ch2 = (Tpatch[5, 5] < Tpatch[4, 4]) && isfinite(Tpatch[4, 4])
            if ch1 && ch2
                Tm2[3] = min(
                    (4 * Tpatch[2, 2] - Tpatch[1, 1]) / 3,
                    (4 * Tpatch[4, 4] - Tpatch[5, 5]) / 3,
                )
                Order[3] = 2
            elseif ch1
                Tm2[3] = (4 * Tpatch[2, 2] - Tpatch[1, 1]) / 3
                Order[3] = 2
            elseif ch2
                Tm2[3] = (4 * Tpatch[4, 4] - Tpatch[5, 5]) / 3
                Order[3] = 2
            end

            ch1 = (Tpatch[1, 5] < Tpatch[2, 4]) && isfinite(Tpatch[2, 4])
            ch2 = (Tpatch[5, 1] < Tpatch[4, 2]) && isfinite(Tpatch[4, 2])
            if ch1 && ch2
                Tm2[4] = min(
                    (4 * Tpatch[2, 4] - Tpatch[1, 5]) / 3,
                    (4 * Tpatch[4, 2] - Tpatch[5, 1]) / 3,
                )
                Order[4] = 2
            elseif ch1
                Tm2[4] = (4 * Tpatch[2, 4] - Tpatch[1, 5]) / 3
                Order[4] = 2
            elseif ch2
                Tm2[4] = (4 * Tpatch[4, 2] - Tpatch[5, 1]) / 3
                Order[4] = 2
            end
        end
    else
        fill!(Tm2, 0)
    end

    # Calculate the distance using x and y direction
    Coeff[1] = 0
    Coeff[2] = 0
    Coeff[3] = -1 / (max(Fij^2, eps(T)))
    for t in 1:2
        if Order[t] == 1
            Coeff[1] += 1
            Coeff[2] -= 2 * Tm[t]
            Coeff[3] += Tm[t]^2
        elseif Order[t] == 2
            Coeff[1] += 2.25
            Coeff[2] -= 2.25 * 2 * Tm2[t]
            Coeff[3] += 2.25 * Tm2[t]^2
        end
    end

    roots!(TT, Coeff)
    # Calculate the distance using the cross directions
    if usecross
        # Coeff[1] += 0
        # Coeff[2] += 0
        Coeff[3] += -1 / (max(Fij^2, eps(T)))
        for t in 3:4
            if Order[t] == 1
                Coeff[1] += 0.5
                Coeff[2] -= Tm[t]
                Coeff[3] += 0.5 * Tm[t]^2
            elseif Order[t] == 2
                Coeff[1] += 0.5 * 2.25
                Coeff[2] -= 2.25 * Tm2[t]
                Coeff[3] += 0.5 * 2.25 * Tm2[t]^2
            end
        end
        roots!(TT2, Coeff)
        # Select minimum distance value of both stensils
        if ~isempty(TT2)
            TT = min(maximum(TT), maximum(TT2))
        end
    else
        TT = maximum(TT)
    end

    # Upwind condition check, current distance must be larger
    # then direct neighbours used in solution
    # DirectNeigbInSol=Tm(isfinite(Tm));
    # if(nnz(DirectNeigbInSol>=Tt)>0) # Will this ever happen?
    #    Tt=min(DirectNeigbInSol)+(1/(max(Fij,eps)));
    # end
    return minimum(TT)
end

function roots!(z, Coeff::AbstractArray{T}) where {T}
    a = Coeff[1]
    b = Coeff[2]
    c = Coeff[3]
    d = max(b^2 - 4.0 * a * c, zero(T))
    if a != 0
        z[1] = (-b - sqrt(d)) / 2a
        z[2] = (-b + sqrt(d)) / 2a
    else
        z[1] = 2c / (-b - sqrt(d))
        z[2] = 2c / (-b + sqrt(d))
    end
end

function roots(Coeff::AbstractArray{T}) where {T}
    z = zeros(T, 2, 1)
    roots!(z, Coeff)
    return z
end


text:
using LinearAlgebra
using Plots
using Trixi

function compute_spectrum(volume_flux, alpha)
  # semidiscretization of the compressible Euler equations
  equations = HartensCompressibleEulerEquations2D(1.4, alpha)

  initial_condition = initial_condition_density_wave

  solver = DGSEM(polydeg=5, surface_flux=volume_flux,
                 volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

  coordinates_min = (-1.0, -1.0)
  coordinates_max = ( 1.0,  1.0)
  mesh = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=2,
                  n_cells_max=30_000)

  semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

  # compute the spectrum
  J = jacobian_ad_forward(semi)
  λ = eigvals(J)
end

function plot_spectrum(volume_flux, alpha, filename)
  λ = compute_spectrum(volume_flux, alpha)

  fontsizes = (xtickfontsize=18, ytickfontsize=18,
               xguidefontsize=20,  yguidefontsize=20, legendfontsize=18)
  kwargs = (xguide="\$\\mathrm{Re}(\\lambda)\$", yguide="\$\\mathrm{Im}(\\lambda)\$",
            legend=:none, size=(700, 500), fontsizes...)
  scatter(real.(λ), imag.(λ); kwargs...)
  savefig(filename)
end

# β = (α + γ) / (1 - γ) = 1 for α = 1 - 2γ (γ = 1.4)
plot_spectrum(flux_sjögreen_yee, -1.8,
              joinpath(dirname(@__DIR__), "figures", "spectrum_beta_1.pdf"))

# β = (α + γ) / (1 - γ) = 2 for α = 2 - 3γ (γ = 1.4)
plot_spectrum(flux_sjögreen_yee, -2.2,
              joinpath(dirname(@__DIR__), "figures", "spectrum_beta_2.pdf"))

# β = (α + γ) / (1 - γ) = -6 for α = 1 (γ = 1.4)
plot_spectrum(flux_sjögreen_yee, 1.0,
              joinpath(dirname(@__DIR__), "figures", "spectrum_beta_m6.pdf"))


text:
function normalizegraph(graph::ModelGraph)
    n = 1
    for node in getnodes(graph)
        m = getmodel(node)
        if m.objSense == :Max
            m.objSense = :Min
            m.obj = -m.obj
            n = -1
        end
    end
    setattribute(graph, :normalized, n)
    return n
end

function fix(var::JuMP.Variable,value::Real)
  setlowerbound(var,value)
  setupperbound(var,value)
end

"""
Checks if n1 is a child node of n2
"""
ischildnode(graph::ModelGraph, n1::ModelNode, n2::ModelNode) = in(n2,in_neighbors(graph,n1))

function savenodeobjective(mf::JuMP.Model)
    g = mf.ext[:Graph]
    numnodes = length(getnodes(g))
    nodeindex = Dict("node$i" => i for i in 1:numnodes)
    nov = mf.ext[:nodeobj] = [AffExpr(0.0) for i in 1:numnodes]
    obj = mf.obj.aff
    for (i,var) in enumerate(obj.vars)
        coeff = obj.coeffs[i]
        varname = mf.colNames[var.col]
        nodename = varname[1:search(varname,'.')-1]
        index = nodeindex[nodename]
        push!(nov[index],coeff,var)
    end
end


text:
##### GOC Data Tools #####

##### Generic Helper Functions #####

function _remove_comment(string)
    return split(string, "/")[1]
end

##### GOC Initialization File Parser (.ini) #####

"""
Given an `.ini` file and a scenario id (i.e. a directory name), parses the
files from ARPA-e GOC Challenge 1 data format.
* `.raw` network model
* `.rop` generator costs
* `.con` contingency set
* `.inl` generator contingency response parameters
"""
function parse_goc_files(ini_file; scenario_id="")
    files, scenario_id = find_goc_files(ini_file, scenario_id=scenario_id)
    return parse_goc_files(files["con"], files["inl"], files["raw"], files["rop"], ini_file=ini_file, scenario_id=scenario_id)
end

function find_goc_files(ini_file; scenario_id="")
    files = Dict(
        "rop" => "x",
        "raw" => "x",
        "con" => "x",
        "inl" => "x"
    )

    if !endswith(ini_file, ".ini")
        warn(_LOGGER, "given init file does not end with .ini, $(ini_file)")
    end

    open(ini_file) do io
        for line in readlines(io)
            line = strip(line)
            #println(line)
            if startswith(line, "[INPUTS]")
                # do nothing
            elseif startswith(line, "ROP")
                files["rop"] = strip(split(line,"=")[2])
            elseif startswith(line, "RAW")
                files["raw"] = strip(split(line,"=")[2])
            elseif startswith(line, "CON")
                files["con"] = strip(split(line,"=")[2])
            elseif startswith(line, "INL")
                files["inl"] = strip(split(line,"=")[2])
            else
                warn(_LOGGER, "unknown input given in ini file: $(line)")
            end
        end
    end

    #println(files)

    ini_dir = dirname(ini_file)

    #println(ini_dir)
    scenario_dirs = [file for file in readdir(ini_dir) if isdir(joinpath(ini_dir, file))]
    scenario_dirs = sort(scenario_dirs)
    #println(scenario_dirs)

    if length(scenario_id) == 0
        scenario_id = scenario_dirs[1]
        info(_LOGGER, "no scenario specified, selected directory \"$(scenario_id)\"")
    else
        if !(scenario_id in scenario_dirs)
            error(_LOGGER, "$(scenario_id) not found in $(scenario_dirs)")
        end
    end

    for (id, path) in files
        if path == "."
            files[id] = ini_dir
        elseif path == "x"
            files[id] = joinpath(ini_dir, scenario_id)
        else
            error(_LOGGER, "unknown file path directive $(path) for file $(id)")
        end
    end

    files["raw"] = joinpath(files["raw"], "case.raw")
    files["rop"] = joinpath(files["rop"], "case.rop")
    files["inl"] = joinpath(files["inl"], "case.inl")
    files["con"] = joinpath(files["con"], "case.con")

    return files, scenario_id
end

function parse_goc_files(con_file, inl_file, raw_file, rop_file; ini_file="", scenario_id="none")
    files = Dict(
        "rop" => rop_file,
        "raw" => raw_file,
        "con" => con_file,
        "inl" => inl_file
    )

    info(_LOGGER, "Parsing Files")
    info(_LOGGER, "  raw: $(files["raw"])")
    info(_LOGGER, "  rop: $(files["rop"])")
    info(_LOGGER, "  inl: $(files["inl"])")
    info(_LOGGER, "  con: $(files["con"])")

    info(_LOGGER, "skipping power models data warnings")
    pm_logger_level = getlevel(getlogger(PowerModels))
    setlevel!(getlogger(PowerModels), "error")
    network_model = _PM.parse_file(files["raw"], import_all=true)
    #network_model = parse_psse(files["raw"], import_all=true)
    #@time network_model = parse_psse(files["raw"], import_all=true)

    setlevel!(getlogger(PowerModels), pm_logger_level)

    gen_cost = parse_rop_file(files["rop"])
    response = parse_inl_file(files["inl"])
    contingencies = parse_con_file(files["con"])

    return (ini_file=ini_file, scenario=scenario_id, network=network_model, cost=gen_cost, response=response, contingencies=contingencies, files=files)
end

function parse_goc_opf_files(ini_file; scenario_id="")
    files = Dict(
        "rop" => "x",
        "raw" => "x",
    )

    if !endswith(ini_file, ".ini")
        warn(_LOGGER, "given init file does not end with .ini, $(ini_file)")
    end

    open(ini_file) do io
        for line in readlines(io)
            line = strip(line)
            #println(line)
            if startswith(line, "[INPUTS]")
                # do nothing
            elseif startswith(line, "ROP")
                files["rop"] = strip(split(line,"=")[2])
            elseif startswith(line, "RAW")
                files["raw"] = strip(split(line,"=")[2])
            else
                warn(_LOGGER, "unknown input given in ini file: $(line)")
            end
        end
    end

    #println(files)

    ini_dir = dirname(ini_file)

    #println(ini_dir)
    scenario_dirs = [file for file in readdir(ini_dir) if isdir(joinpath(ini_dir, file))]
    scenario_dirs = sort(scenario_dirs)
    #println(scenario_dirs)

    if length(scenario_id) == 0
        scenario_id = scenario_dirs[1]
        info(_LOGGER, "no scenario specified, selected directory \"$(scenario_id)\"")
    else
        if !(scenario_id in scenario_dirs)
            error(_LOGGER, "$(scenario_id) not found in $(scenario_dirs)")
        end
    end

    for (id, path) in files
        if path == "."
            files[id] = ini_dir
        elseif path == "x"
            files[id] = joinpath(ini_dir, scenario_id)
        else
            error(_LOGGER, "unknown file path directive $(path) for file $(id)")
        end
    end

    files["raw"] = joinpath(files["raw"], "case.raw")
    files["rop"] = joinpath(files["rop"], "case.rop")

    info(_LOGGER, "Parsing Files")
    info(_LOGGER, "  raw: $(files["raw"])")
    info(_LOGGER, "  rop: $(files["rop"])")

    network_model = _PM.parse_file(files["raw"], import_all=true)
    gen_cost = parse_rop_file(files["rop"])

    return (ini_file=ini_file, scenario=scenario_id, network=network_model, cost=gen_cost, files=files)
end

##### Unit Inertia and Governor Response Data File Parser (.inl) #####

function parse_inl_file(file::String)
    open(file) do io
        return parse_inl_file(io)
    end
end

function parse_inl_file(io::IO)
    inl_list = []
    for line in readlines(io)
        #line = _remove_comment(line)

        if startswith(strip(line), "0")
            debug(_LOGGER, "inl file sentinel found")
            break
        end
        line_parts = split(line, ",")
        @assert length(line_parts) >= 7

        inl_data = Dict(
            "i"    => parse(Int, line_parts[1]),
            "id"   => strip(line_parts[2]),
            "h"    => strip(line_parts[3]),
            "pmax" => strip(line_parts[4]),
            "pmin" => strip(line_parts[5]),
            "r"    => parse(Float64, line_parts[6]),
            "d"    => strip(line_parts[7])
        )

        @assert inl_data["r"] >= 0.0

        #println(inl_data)
        push!(inl_list, inl_data)
    end
    return inl_list
end

##### Generator Cost Data File Parser (.rop) #####

_rop_sections = [
    "mod" => "Modification Code",
    "bus_vm" => "Bus Voltage Attributes",
    "shunt_adj" => "Adjustable Bus Shunts",
    "load" => "Bus Loads",
    "load_adj" => "Adjustable Bus Load Tables",
    "gen" => "Generator Dispatch Units",
    "disptbl" => "Active Power Dispatch Tables",
    "gen_reserve" => "Generator Reserve Units",
    "qg" => "Generation Reactive Capability",
    "branch_x" => "Adjustable Branch Reactance",
    "ctbl" => "Piecewise Linear Cost Curve Tables",
    "pwc" => "Piecewise Quadratic Cost Curve Tables",
    "pec" => "Polynomial & Exponential Cost Curve Tables",
    "reserve" => "Period Reserves",
    "branch_flow" => "Branch Flows",
    "int_flow" => "Interface Flows",
    "lin_const" => "Linear Constraint Dependencies",
    "dc_const" => "Two Terminal DC Line Constraint Dependencies",
]

function parse_rop_file(file::String)
    open(file) do io
        return parse_rop_file(io)
    end
end

function parse_rop_file(io::IO)
    active_section_idx = 1
    active_section = _rop_sections[active_section_idx]

    section_data = Dict()
    section_data[active_section.first] = []

    line_idx = 1
    lines = readlines(io)
    while line_idx < length(lines)
        #line = _remove_comment(lines[line_idx])
        line = lines[line_idx]
        if startswith(strip(line), "0")
            debug(_LOGGER, "finished reading rop section $(active_section.second) with $(length(section_data[active_section.first])) items")
            active_section_idx += 1
            if active_section_idx > length(_rop_sections)
                debug(_LOGGER, "finished reading known rop sections")
                break
            end
            active_section = _rop_sections[active_section_idx]
            section_data[active_section.first] = []
            line_idx += 1
            continue
        end

        if active_section.first == "gen"
            push!(section_data[active_section.first], _parse_rop_gen(line))
        elseif active_section.first == "disptbl"
            push!(section_data[active_section.first], _parse_rop_pg(line))
        elseif active_section.first == "ctbl"
            pwl_line_parts = split(line, ",")
            @assert length(pwl_line_parts) >= 3

            num_pwl_lines = parse(Int, pwl_line_parts[3])
            @assert num_pwl_lines > 0

            pwl_point_lines = lines[line_idx+1:line_idx+num_pwl_lines]
            #pwl_point_lines = remove_comment.(pwl_point_lines)
            push!(section_data[active_section.first], _parse_rop_pwl(pwl_line_parts, pwl_point_lines))
            line_idx += num_pwl_lines
        else
            info(_LOGGER, "skipping data line: $(line)")
        end
        line_idx += 1
    end
    return section_data
end

function _parse_rop_gen(line)
    line_parts = split(line, ",")
    @assert length(line_parts) >= 4

    data = Dict(
        "bus"     => parse(Int, line_parts[1]),
        "genid"   => strip(line_parts[2]),
        "disp"    => strip(line_parts[3]),
        "disptbl" => parse(Int, line_parts[4]),
    )

    @assert data["disptbl"] >= 0

    return data
end

function _parse_rop_pg(line)
    line_parts = split(line, ",")
    @assert length(line_parts) >= 7

    data = Dict(
        "tbl"      => parse(Int, line_parts[1]),
        "pmax"     => strip(line_parts[2]),
        "pmin"     => strip(line_parts[3]),
        "fuelcost" => strip(line_parts[4]),
        "ctyp"     => strip(line_parts[5]),
        "status"   => strip(line_parts[6]),
        "ctbl"     => parse(Int, line_parts[7]),
    )

    @assert data["tbl"] >= 0
    @assert data["ctbl"] >= 0

    return data
end

function _parse_rop_pwl(pwl_parts, point_lines)
    @assert length(pwl_parts) >= 2

    points = []

    for point_line in point_lines
        point_line_parts = split(point_line, ",")
        @assert length(point_line_parts) >= 2
        x = parse(Float64, point_line_parts[1])
        y = parse(Float64, point_line_parts[2])

        push!(points, (x=x, y=y))
    end

    data = Dict(
        "ltbl"   =>  parse(Int, pwl_parts[1]),
        "label"  => strip(pwl_parts[2]),
        "points" => points
    )

    @assert data["ltbl"] >= 0

    return data
end

##### Contingency Description Data File (.con) #####

# OPEN BRANCH FROM BUS *I TO BUS *J CIRCUIT *1CKT
_branch_contigency_structure = [
    1 => "OPEN",
    2 => "BRANCH",
    3 => "FROM",
    4 => "BUS",
    #5 => "I",
    6 => "TO",
    7 => "BUS",
    #8 => "J",
    9 => "CIRCUIT",
    #10 => "CKT"
]

# REMOVE UNIT *ID FROM BUS *I
_generator_contigency_structure = [
    1 => "REMOVE",
    2 => "UNIT",
    #3 => "ID",
    4 => "FROM",
    5 => "BUS",
    #6 => "I"
]

function parse_con_file(file::String)
    open(file) do io
        return parse_con_file(io)
    end
end

function parse_con_file(io::IO)
    con_lists = []

    tokens = []

    for line in readlines(io)
        #line_tokens = split(strip(_remove_comment(line)))
        line_tokens = split(strip(line))
        #println(line_tokens)
        append!(tokens, line_tokens)
    end

    #println(tokens)

    token_idx = 1
    while token_idx <= length(tokens)
        token = tokens[token_idx]
        if token == "END"
            debug(_LOGGER, "end of contingency file found")
            break
        elseif token == "CONTINGENCY"
            # start reading contingencies

            contingency_name = tokens[token_idx+1]
            debug(_LOGGER, "reading contingency $(contingency_name)")

            token_idx += 2
            token = tokens[token_idx]
            remaining_tokens = length(tokens) - token_idx

            if token == "OPEN" # branch contingency case
                # OPEN BRANCH FROM BUS *I TO BUS *J CIRCUIT *1CKT

                @assert remaining_tokens >= 9
                branch_tokens = tokens[token_idx:token_idx+9]
                #println(branch_tokens)

                #if !all(branch_tokens[idx] == val for (idx, val) in _branch_contigency_structure) && !all(branch_tokens[idx] == val for (idx, val) in _branch_contigency_structure_alt)
                if any(branch_tokens[idx] != val for (idx, val) in _branch_contigency_structure)
                    error(_LOGGER, "incorrect branch contingency structure: $(branch_tokens)")
                end

                bus_i = parse(Int, branch_tokens[5])
                @assert bus_i >= 0

                bus_j = parse(Int, branch_tokens[8])
                @assert bus_j >= 0

                ckt = branch_tokens[10]

                branch_contingency = Dict(
                    "label" => contingency_name,
                    "component" => "branch",
                    "action" => "open",
                    "i" => bus_i,
                    "j" => bus_j,
                    "ckt" => ckt,
                )

                push!(con_lists, branch_contingency)

                token_idx += 9
            elseif token == "REMOVE"
                # REMOVE UNIT *ID FROM BUS *I

                @assert remaining_tokens >= 5
                generator_tokens = tokens[token_idx:token_idx+5]
                #println(generator_tokens)

                if any(generator_tokens[idx] != val for (idx, val) in _generator_contigency_structure)
                    error(_LOGGER, "incorrect generator contingency structure: $(generator_tokens)")
                end

                gen_id = generator_tokens[3]

                bus_i = parse(Int, generator_tokens[6])
                @assert bus_i >= 0

                generator_contingency = Dict(
                    "label" => contingency_name,
                    "component" => "generator",
                    "action" => "remove",
                    "id" => gen_id,
                    "i" => bus_i,
                )

                push!(con_lists, generator_contingency)

                token_idx += 5
            elseif token == "END"
                warn(_LOGGER, "no action provided for contingency $(contingency_name)")
                token_idx -= 1
            else
                warn(_LOGGER, "unrecognized token $(token)")
            end

            token_idx += 1
            token = tokens[token_idx]
            if token != "END"
                error(_LOGGER, "expected END token at end of CONTINGENCY, got $(token)")
            end
        else
            warn(_LOGGER, "unrecognized token $(token)")
        end
        token_idx += 1
    end

    return con_lists
end

function parse_solution1_file(file::String)
    open(file) do io
        return parse_solution1_file(io)
    end
end

function parse_solution1_file(io::IO)
    bus_data_list = []
    gen_data_list = []

    lines = readlines(io)

    # skip bus list header section
    idx = 1

    separator_count = 0
    skip_next = false

    while idx <= length(lines)
        line = lines[idx]
        if length(strip(line)) == 0
            warn(_LOGGER, "skipping blank line in solution1 file ($(idx))")
        elseif skip_next
            skip_next = false
        elseif startswith(strip(line), "--")
            separator_count += 1
            skip_next = true
        else
            if separator_count == 1
                parts = split(line, ",")
                @assert length(parts) >= 4
                bus_data = (
                    bus = parse(Int, parts[1]),
                    vm = parse(Float64, parts[2]),
                    va = parse(Float64, parts[3]),
                    bcs = parse(Float64, parts[4])
                )
                push!(bus_data_list, bus_data)
            elseif separator_count == 2
                parts = split(line, ",")
                @assert length(parts) >= 4
                gen_data = (
                    bus = parse(Int, parts[1]),
                    id = strip(strip(parts[2]), ['\'', ' ']),
                    pg = parse(Float64, parts[3]),
                    qg = parse(Float64, parts[4])
                )
                push!(gen_data_list, gen_data)
            else
                warn(_LOGGER, "skipping line in solution1 file ($(idx)): $(line)")
            end
        end
        idx += 1
    end

    return (bus=bus_data_list, gen=gen_data_list)
end

##### Transform GOC Data in PM Data #####

"""
transforms files from ARPA-e GOC Challenge 1 data format in to the PowerModels
data format.  This consists of taking the data from multiple data structures
and putting into a network data dictionary.
"""
function build_pm_model(goc_data)
    scenario = goc_data.scenario
    network = goc_data.network

    ##### General Helpers #####
    gen_lookup = Dict(tuple(gen["source_id"][2], strip(gen["source_id"][3])) => gen for (i,gen) in network["gen"])

    branch_lookup = Dict()
    for (i,branch) in network["branch"]
        if !branch["transformer"]
            branch_id = tuple(branch["source_id"][2], branch["source_id"][3], strip(branch["source_id"][4]))
        else
            branch_id = tuple(branch["source_id"][2], branch["source_id"][3], strip(branch["source_id"][5]))
            @assert branch["source_id"][4] == 0
            @assert branch["source_id"][6] == 0
        end
        branch_lookup[branch_id] = branch
    end

    ##### Link Generator Cost Data #####

    @assert network["per_unit"]
    mva_base = network["baseMVA"]

    dispatch_tbl_lookup = Dict()
    for dispatch_tbl in goc_data.cost["disptbl"]
        dispatch_tbl_lookup[dispatch_tbl["ctbl"]] = dispatch_tbl
    end

    cost_tbl_lookup = Dict()
    for cost_tbl in goc_data.cost["ctbl"]
        cost_tbl_lookup[cost_tbl["ltbl"]] = cost_tbl
    end

    gen_cost_models = Dict()
    for gen_dispatch in goc_data.cost["gen"]
        gen_id = (gen_dispatch["bus"], strip(replace(gen_dispatch["genid"], "'" => "")))
        dispatch_tbl = dispatch_tbl_lookup[gen_dispatch["disptbl"]]
        cost_tbl = cost_tbl_lookup[dispatch_tbl["ctbl"]]

        gen_cost_models[gen_id] = cost_tbl
    end

    if length(gen_cost_models) != length(network["gen"])
        error(_LOGGER, "cost model data missing, network has $(length(network["gen"])) generators, the cost model has $(length(gen_cost_models)) generators")
    end

    for (gen_id, cost_model) in gen_cost_models
        pm_gen = gen_lookup[gen_id]
        pm_gen["model"] = 1
        pm_gen["model_label"] = cost_model["label"]
        pm_gen["ncost"] = length(cost_model["points"])

        #println(cost_model["points"])
        point_list = Float64[]
        for point in cost_model["points"]
            push!(point_list, point.x/mva_base)
            push!(point_list, point.y)
        end
        pm_gen["cost"] = point_list
    end

    ##### Link Generator Participation Data #####

    if length(goc_data.response) != length(network["gen"])
        error(_LOGGER, "generator response model data missing, network has $(length(network["gen"])) generators, the response model has $(length(goc_data.response)) generators")
    end

    for gen_response in goc_data.response
        gen_id = (gen_response["i"], strip(gen_response["id"]))

        pm_gen = gen_lookup[gen_id]

        pm_gen["alpha"] = gen_response["r"]
    end

    ##### Setup Generator Area Group Data #####

    area_gens = Dict{Int,Set{Int}}()
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            gen_bus = network["bus"]["$(gen["gen_bus"])"]
            area = gen_bus["area"]
            if !haskey(area_gens, area)
                area_gens[area] = Set{Int}()
            end
            push!(area_gens[area], gen["index"])
        end
    end
    network["area_gens"] = area_gens

    ##### Flexible Shunt Data #####

    for (i,shunt) in network["shunt"]
        #println(shunt["source_id"])
        if shunt["source_id"][1] == "switched shunt"
            #@assert shunt["source_id"][3] == 0
            @assert shunt["gs"] == 0.0
            shunt["dispatchable"] = true

            bmin = 0.0
            bmax = 0.0
            for (n_name,b_name) in [("n1","b1"),("n2","b2"),("n3","b3"),("n4","b4"),("n5","b5"),("n6","b6"),("n7","b7"),("n8","b8")]
                if shunt[b_name] <= 0.0
                    bmin += shunt[n_name]*shunt[b_name]
                else
                    bmax += shunt[n_name]*shunt[b_name]
                end
            end
            shunt["bmin"] = bmin/mva_base
            shunt["bmax"] = bmax/mva_base
        else
            shunt["dispatchable"] = false
        end
    end

    ##### Add Contingency Lists #####

    generator_ids = []
    branch_ids = []

    for (i,cont) in enumerate(goc_data.contingencies)
        if cont["component"] == "branch"
            branch_id = (cont["i"], cont["j"], cont["ckt"])
            pm_branch = branch_lookup[branch_id]
            push!(branch_ids, (idx=pm_branch["index"], label=cont["label"], type="branch"))

        elseif cont["component"] == "generator"
            gen_id = (cont["i"], cont["id"])
            pm_gen = gen_lookup[gen_id]
            push!(generator_ids, (idx=pm_gen["index"], label=cont["label"], type="gen"))

        else
            error(_LOGGER, "unrecognized contingency component type $(cont["component"]) at contingency $(i)")
        end
    end

    network["branch_contingencies"] = branch_ids
    network["gen_contingencies"] = generator_ids

    network["branch_contingencies_active"] = []
    network["gen_contingencies_active"] = []

    ##### Fix Broken Data #####

    _PM.correct_cost_functions!(network)

    # FYI, this breaks output API
    #_PM.propagate_topology_status!(network)

    for (i,shunt) in network["shunt"]
        # test checks if a "switched shunt" in the orginal data model
        if shunt["dispatchable"]
            if shunt["bs"] < shunt["bmin"]
                warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmin"])")
                shunt["bs"] = shunt["bmin"]
            end
            if shunt["bs"] > shunt["bmax"]
                warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmax"])")
                shunt["bs"] = shunt["bmax"]
            end
        end
    end

    return network
end

"a simpler version of `build_pm_model` that does not require contingency information"
function build_pm_opf_model(goc_data)
    scenario = goc_data.scenario
    network = goc_data.network

    ##### General Helpers #####

    gen_lookup = Dict(tuple(gen["source_id"][2], strip(gen["source_id"][3])) => gen for (i,gen) in network["gen"])

    branch_lookup = Dict()
    for (i,branch) in network["branch"]
        if !branch["transformer"]
            branch_id = tuple(branch["source_id"][2], branch["source_id"][3], strip(branch["source_id"][4]))
        else
            branch_id = tuple(branch["source_id"][2], branch["source_id"][3], strip(branch["source_id"][5]))
            @assert branch["source_id"][4] == 0
            @assert branch["source_id"][6] == 0
        end
        branch_lookup[branch_id] = branch
    end

    ##### Link Generator Cost Data #####

    @assert network["per_unit"]
    mva_base = network["baseMVA"]

    dispatch_tbl_lookup = Dict()
    for dispatch_tbl in goc_data.cost["disptbl"]
        dispatch_tbl_lookup[dispatch_tbl["ctbl"]] = dispatch_tbl
    end

    cost_tbl_lookup = Dict()
    for cost_tbl in goc_data.cost["ctbl"]
        cost_tbl_lookup[cost_tbl["ltbl"]] = cost_tbl
    end

    gen_cost_models = Dict()
    for gen_dispatch in goc_data.cost["gen"]
        gen_id = (gen_dispatch["bus"], strip(gen_dispatch["genid"]))
        dispatch_tbl = dispatch_tbl_lookup[gen_dispatch["disptbl"]]
        cost_tbl = cost_tbl_lookup[dispatch_tbl["ctbl"]]

        gen_cost_models[gen_id] = cost_tbl
    end

    if length(gen_cost_models) != length(network["gen"])
        error(_LOGGER, "cost model data missing, network has $(length(network["gen"])) generators, the cost model has $(length(gen_cost_models)) generators")
    end

    for (gen_id, cost_model) in gen_cost_models
        pm_gen = gen_lookup[gen_id]
        pm_gen["model"] = 1
        pm_gen["model_label"] = cost_model["label"]
        pm_gen["ncost"] = length(cost_model["points"])

        #println(cost_model["points"])
        point_list = Float64[]
        for point in cost_model["points"]
            push!(point_list, point.x/mva_base)
            push!(point_list, point.y)
        end
        pm_gen["cost"] = point_list
    end

    ##### Flexible Shunt Data #####

    for (i,shunt) in network["shunt"]
        if shunt["source_id"][1] == "switched shunt"
            @assert shunt["source_id"][3] == 0
            @assert shunt["gs"] == 0.0
            shunt["dispatchable"] = true

            bmin = 0.0
            bmax = 0.0
            for (n_name,b_name) in [("n1","b1"),("n2","b2"),("n3","b3"),("n4","b4"),("n5","b5"),("n6","b6"),("n7","b7"),("n8","b8")]
                if shunt[b_name] <= 0.0
                    bmin += shunt[n_name]*shunt[b_name]
                else
                    bmax += shunt[n_name]*shunt[b_name]
                end
            end
            shunt["bmin"] = bmin/mva_base
            shunt["bmax"] = bmax/mva_base
        else
            shunt["dispatchable"] = false
        end
    end

    ##### Fix Broken Data #####

    _PM.correct_cost_functions!(network)

    # FYI, this breaks output API
    #_PM.propagate_topology_status!(network)

    for (i,shunt) in network["shunt"]
        # test checks if a "switched shunt" in the orginal data model
        if shunt["dispatchable"]
            if shunt["bs"] < shunt["bmin"]
                warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmin"])")
                shunt["bs"] = shunt["bmin"]
            end
            if shunt["bs"] > shunt["bmax"]
                warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmax"])")
                shunt["bs"] = shunt["bmax"]
            end
        end
    end

    return network
end

##### Working with GOC Solution Files #####

# this Does not appear to be relivent (03/14/2020)
# function build_contingency_solutions(network, solution)
#     contingency_solutions = Dict{String,Any}()

#     solution = deepcopy(solution)
#     solution["delta"] = 0.0

#     for cont in network["gen_contingencies"]
#         contingency_solutions[cont.label] = solution
#     end

#     for cont in network["branch_contingencies"]
#         contingency_solutions[cont.label] = solution
#     end

#     return contingency_solutions
# end

"checks feasibility criteria of network solution, produces an error if a problem is found"
function check_network_solution(network)
    for (i,bus) in network["bus"]
        if bus["bus_type"] != 4
            if bus["vm"] > bus["vmax"] || bus["vm"] < bus["vmin"]
                error(_LOGGER, "vm on $(bus["source_id"]) is not in bounds $(bus["vmin"]) to $(bus["vmax"]), given $(bus["vm"])")
            end
        end
    end

    for (i,shunt) in network["shunt"]
        if shunt["status"] != 0
            if haskey(shunt, "dispatchable")
                if shunt["dispatchable"]
                    @assert shunt["gs"] == 0.0
                    @assert haskey(shunt, "bmin") && haskey(shunt, "bmax")
                    if shunt["bs"] > shunt["bmax"] || shunt["bs"] < shunt["bmin"]
                        error(_LOGGER, "bs on $(shunt["source_id"]) is not in bounds $(shunt["bmin"]) to $(shunt["bmax"]), given $(shunt["bs"])")
                    end
                end
            end
        end
    end

    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            if gen["pg"] > gen["pmax"] || gen["pg"] < gen["pmin"]
                error(_LOGGER, "pg on gen $(gen["source_id"]) not in bounds $(gen["pmin"]) to $(gen["pmax"]), given $(gen["pg"])")
            end

            if gen["qg"] > gen["qmax"] || gen["qg"] < gen["qmin"]
                error(_LOGGER, "pg on gen $(gen["source_id"]) not in bounds $(gen["qmin"]) to $(gen["qmax"]), given $(gen["qg"])")
            end
        end
    end
end

"checks feasibility criteria of network solution, corrects when possible"
function correct_network_solution!(network)

    # default value is required for correctness of max and mean computations when no changes are made
    vm_changes = [0.0]
    for (i,bus) in network["bus"]
        if bus["bus_type"] != 4
            if bus["vm"] > bus["vmax"]
                warn(_LOGGER, "update vm on bus $(i) to be in bounds $(bus["vm"]) -> $(bus["vmax"])")
                push!(vm_changes, bus["vm"] - bus["vmax"])
                bus["vm"] = bus["vmax"]
            end

            if bus["vm"] < bus["vmin"]
                warn(_LOGGER, "update vm on bus $(i) to be in bounds $(bus["vm"]) -> $(bus["vmin"])")
                push!(vm_changes, bus["vmin"] - bus["vm"])
                bus["vm"] = bus["vmin"]
            end
        else
            bus["vm"] = 0.0
            bus["va"] = 0.0
        end
    end

    bs_changes = [0.0]
    for (i,shunt) in network["shunt"]
        if haskey(shunt, "dispatchable")
            if shunt["dispatchable"]
                @assert shunt["gs"] == 0.0
                @assert haskey(shunt, "bmin") && haskey(shunt, "bmax")
                if shunt["bs"] > shunt["bmax"]
                    warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmax"])")
                    push!(bs_changes, shunt["bs"] - shunt["bmax"])
                    shunt["bs"] = shunt["bmax"]
                end
                if shunt["bs"] < shunt["bmin"]
                    warn(_LOGGER, "update bs on shunt $(i) to be in bounds $(shunt["bs"]) -> $(shunt["bmin"])")
                    push!(bs_changes, shunt["bmin"] - shunt["bs"])
                    shunt["bs"] = shunt["bmin"]
                end
            end
        else
            warn(_LOGGER, "shunt $(i) missing dispatchable parameter")
        end
    end

    pg_changes = [0.0]
    qg_changes = [0.0]
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            if gen["pg"] > gen["pmax"]
                warn(_LOGGER, "update pg on gen $(i) to be in bounds $(gen["pg"]) -> $(gen["pmax"])")
                push!(pg_changes, gen["pg"] - gen["pmax"])
                gen["pg"] = gen["pmax"]
            end
            if gen["pg"] < gen["pmin"]
                warn(_LOGGER, "update pg on gen $(i) to be in bounds $(gen["pg"]) -> $(gen["pmin"])")
                push!(pg_changes, gen["pmin"] - gen["pg"])
                gen["pg"] = gen["pmin"]
            end

            if gen["qg"] > gen["qmax"]
                warn(_LOGGER, "update qg on gen $(i) to be in bounds $(gen["qg"]) -> $(gen["qmax"])")
                push!(qg_changes, gen["qg"] - gen["qmax"])
                gen["qg"] = gen["qmax"]
            end
            if gen["qg"] < gen["qmin"]
                warn(_LOGGER, "update qg on gen $(i) to be in bounds $(gen["qg"]) -> $(gen["qmin"])")
                push!(qg_changes, gen["qmin"] - gen["qg"])
                gen["qg"] = gen["qmin"]
            end
        else
            gen["pg"] = 0.0
            gen["qg"] = 0.0
        end
    end

    _summary_changes(network, "base_case", vm_changes, bs_changes, pg_changes, qg_changes)
end

"checks feasibility criteria of contingencies, corrects when possible"
function correct_contingency_solutions!(network, contingency_solutions)
    bus_gens = gens_by_bus(network)

    cont_changes = Int64[]
    cont_vm_changes_max = [0.0]
    cont_bs_changes_max = [0.0]
    cont_pg_changes_max = [0.0]
    cont_qg_changes_max = [0.0]

    for cont_sol in contingency_solutions
        changes = correct_contingency_solution!(network, cont_sol; bus_gens=bus_gens)

        push!(cont_changes, changes.changed)
        push!(cont_vm_changes_max, changes.vm_changes_max)
        push!(cont_bs_changes_max, changes.bs_changes_max)
        push!(cont_pg_changes_max, changes.pg_changes_max)
        push!(cont_qg_changes_max, changes.qg_changes_max)
    end

    println("")

    data = [
        "----",
        "bus",
        "branch",
        "gen_cont",
        "branch_cont",
        "changes_count",
        "vm_max_max",
        "bs_max_max",
        "pg_max_max",
        "qg_max_max",
        "vm_max_mean",
        "bs_max_mean",
        "pg_max_mean",
        "qg_max_mean",
    ]
    println(join(data, ", "))

    data = [
        "DATA_CCS",
        length(network["bus"]),
        length(network["branch"]),
        length(network["gen_contingencies"]),
        length(network["branch_contingencies"]),
        sum(cont_changes),
        maximum(cont_vm_changes_max),
        maximum(cont_bs_changes_max),
        maximum(cont_pg_changes_max),
        maximum(cont_qg_changes_max),
        mean(cont_vm_changes_max),
        mean(cont_bs_changes_max),
        mean(cont_pg_changes_max),
        mean(cont_qg_changes_max),
    ]
    println(join(data, ", "))
end

"checks feasibility criteria of contingencies, corrects when possible"
function correct_contingency_solution!(network, cont_sol; bus_gens = gens_by_bus(network))
    label = cont_sol["label"]
    vm_changes = [0.0]
    for (i,bus) in cont_sol["bus"]
        nw_bus = network["bus"][i]

        if nw_bus["bus_type"] != 4
            if length(bus_gens[i]) > 0
                qg = sum(cont_sol["gen"]["$(gen["index"])"]["qg"] for gen in bus_gens[i])
                qmin = sum(gen["qmin"] for gen in bus_gens[i])
                qmax = sum(gen["qmax"] for gen in bus_gens[i])

                if !isapprox(abs(qmin - qmax), 0.0)
                    if qg >= qmax && bus["vm"] - vm_eq_tol/10 > nw_bus["vm"]
                        #println(qmin, " ", qg, " ", qmax)
                        #warn(_LOGGER, "update qg $(qmin) $(qg) $(qmax)")
                        warn(_LOGGER, "update vm on bus $(i) in contingency $(label) to match set-point $(bus["vm"]) -> $(nw_bus["vm"]) due to qg upper bound and vm direction")
                        push!(vm_changes, abs(bus["vm"] - nw_bus["vm"]))
                        bus["vm"] = nw_bus["vm"]
                    end

                    if qg <= qmin && bus["vm"] + vm_eq_tol/10 < nw_bus["vm"]
                        #println(qmin, " ", qg, " ", qmax)
                        #warn(_LOGGER, "update qg $(qmin) $(qg) $(qmax)")
                        warn(_LOGGER, "update vm on bus $(i) in contingency $(label) to match set-point $(bus["vm"]) -> $(nw_bus["vm"]) due to qg lower bound and vm direction")
                        push!(vm_changes, abs(bus["vm"] - nw_bus["vm"]))
                        bus["vm"] = nw_bus["vm"]
                    end
                end
            end

            if bus["vm"] > nw_bus["vmax"]
                warn(_LOGGER, "update vm on bus $(i) in contingency $(label) to match ub $(bus["vm"]) -> $(nw_bus["vmax"]) due to out of bounds")
                push!(vm_changes, bus["vm"] - nw_bus["vmax"])
                bus["vm"] = nw_bus["vmax"]
            end

            if bus["vm"] < nw_bus["vmin"]
                warn(_LOGGER, "update vm on bus $(i) in contingency $(label) to match lb $(bus["vm"]) -> $(nw_bus["vmin"]) due to out of bounds")
                push!(vm_changes, nw_bus["vmin"] - bus["vm"])
                bus["vm"] = nw_bus["vmin"]
            end
        else
            bus["vm"] = 0.0
            bus["va"] = 0.0
        end
    end

    bs_changes = [0.0]
    if haskey(cont_sol, "shunt")
        for (i,shunt) in cont_sol["shunt"]
            nw_shunt = network["shunt"][i]
            if haskey(nw_shunt, "dispatchable") && nw_shunt["dispatchable"]
                @assert nw_shunt["gs"] == 0.0
                @assert haskey(nw_shunt, "bmin") && haskey(nw_shunt, "bmax")
                if shunt["bs"] > nw_shunt["bmax"]
                    warn(_LOGGER, "update bs on shunt $(i) in contingency $(label) to be in bounds $(shunt["bs"]) -> $(nw_shunt["bmax"])")
                    push!(bs_changes, shunt["bs"] - nw_shunt["bmax"])
                    shunt["bs"] = nw_shunt["bmax"]
                end
                if shunt["bs"] < nw_shunt["bmin"]
                    warn(_LOGGER, "update bs on shunt $(i) in contingency $(label) to be in bounds $(shunt["bs"]) -> $(nw_shunt["bmin"])")
                    push!(bs_changes, nw_shunt["bmin"] - shunt["bs"])
                    shunt["bs"] = nw_shunt["bmin"]
                end
            end
        end
    end

    response_gens = Set{Int}()

    gen_id = -1
    if cont_sol["cont_type"] == "gen"
        gen_id = cont_sol["cont_comp_id"]
        nw_gen = network["gen"]["$(gen_id)"]
        nw_gen_bus = network["bus"]["$(nw_gen["gen_bus"])"]
        response_gens = network["area_gens"][nw_gen_bus["area"]]
    end

    if cont_sol["cont_type"] == "branch"
        branch_id = cont_sol["cont_comp_id"]
        nw_branch = network["branch"]["$(branch_id)"]
        nw_fr_bus = network["bus"]["$(nw_branch["f_bus"])"]
        nw_to_bus = network["bus"]["$(nw_branch["t_bus"])"]
        response_gens = Set()
        if haskey(network["area_gens"], nw_fr_bus["area"])
            response_gens = network["area_gens"][nw_fr_bus["area"]]
        end
        if nw_fr_bus["area"] != nw_to_bus["area"] && haskey(network["area_gens"], nw_to_bus["area"])
            response_gens = union(response_gens, network["area_gens"][nw_to_bus["area"]])
        end
    end

    pg_changes = [0.0]
    qg_changes = [0.0]
    delta = cont_sol["delta"]
    for (i,gen) in cont_sol["gen"]
        nw_gen = network["gen"][i]

        if !(nw_gen["gen_status"] == 0 || (gen_id >= 0 && nw_gen["index"] == gen_id))
            bus_id = nw_gen["gen_bus"]
            nw_bus = network["bus"]["$(bus_id)"]

            if gen["qg"] < nw_gen["qmax"] && gen["qg"] > nw_gen["qmin"]
                bus = cont_sol["bus"]["$(bus_id)"]
                #if !isapprox(bus["vm"], nw_bus["vm"])
                if !isapprox(bus["vm"], nw_bus["vm"], atol=vm_eq_tol/2)
                    #debug(_LOGGER, "bus $(bus_id) : vm_base $(nw_bus["vm"]) - vm $(bus["vm"]) : reactive bounds $(nw_gen["qmin"]) - $(gen["qg"]) - $(nw_gen["qmax"])")
                    warn(_LOGGER, "update vm on bus $(bus_id) in contingency $(label) to match base case $(bus["vm"]) -> $(nw_bus["vm"]) due to within reactive bounds")
                end
                bus["vm"] = nw_bus["vm"]
            end

            pg_calc = nw_gen["pg"]
            if nw_gen["index"] in response_gens
                pg_calc += nw_gen["alpha"]*delta
                pg_calc = max(pg_calc, nw_gen["pmin"])
                pg_calc = min(pg_calc, nw_gen["pmax"])
            end

            if !isapprox(gen["pg"], pg_calc, atol=1e-5)
                warn(_LOGGER, "pg value on gen $(i) $(nw_gen["source_id"]) in contingency $(label) is not consistent with the computed value given:$(gen["pg"]) calc:$(pg_calc)")
            end

            if gen["pg"] > nw_gen["pmax"]
                warn(_LOGGER, "update pg on gen $(i) $(nw_gen["source_id"]) in contingency $(label) to match ub $(gen["pg"]) -> $(nw_gen["pmax"])")
                push!(pg_changes, gen["pg"] - nw_gen["pmax"])
                gen["pg"] = nw_gen["pmax"]
            end

            if gen["pg"] < nw_gen["pmin"]
                warn(_LOGGER, "update pg on gen $(i) $(nw_gen["source_id"]) in contingency $(label) to match lb $(gen["pg"]) -> $(nw_gen["pmin"])")
                push!(pg_changes, nw_gen["pmin"] - gen["pg"])
                gen["pg"] = nw_gen["pmin"]
            end

            if gen["qg"] > nw_gen["qmax"]
                warn(_LOGGER, "update qg on gen $(i) $(nw_gen["source_id"]) in contingency $(label) to match ub $(gen["qg"]) -> $(nw_gen["qmax"])")
                push!(qg_changes, gen["qg"] - nw_gen["qmax"])
                gen["qg"] = nw_gen["qmax"]
            end

            if gen["qg"] < nw_gen["qmin"]
                warn(_LOGGER, "update qg on gen $(i) $(nw_gen["source_id"]) in contingency $(label) to match lb $(gen["qg"]) -> $(nw_gen["qmin"])")
                push!(qg_changes, nw_gen["qmin"] - gen["qg"])
                gen["qg"] = nw_gen["qmin"]
            end
        else
            gen["pg"] = 0.0
            gen["qg"] = 0.0
        end
    end

    # test imbalance on branch conts after flow corrections
    #=
    network_tmp = deepcopy(network)
    cont_type = cont_sol["cont_type"]
    cont_idx = cont_sol["cont_comp_id"]
    network_tmp["branch"]["$(cont_idx)"]["br_status"] = 0
    _PM.update_data!(network_tmp, cont_sol)
    deltas = compute_power_balance_deltas!(network_tmp)
    println("$(label): $(deltas)")
    =#

    cont_changed = length(vm_changes) > 1 || length(bs_changes) > 1 || length(pg_changes) > 1 || length(qg_changes) > 1

    if cont_changed
        _summary_changes(network, label, vm_changes, bs_changes, pg_changes, qg_changes)
    end

    return (changed=Int(cont_changed), vm_changes_max=maximum(vm_changes), bs_changes_max=maximum(bs_changes), pg_changes_max=maximum(pg_changes), qg_changes_max=maximum(qg_changes))
end

function _summary_changes(network, contingency, vm_changes, bs_changes, pg_changes, qg_changes)
    #println("")

    data = [
        "----",
        "contingency",
        "bus",
        "branch",
        "gen_cont",
        "branch_cont",
        "vm_count",
        "bs_count",
        "pg_count",
        "qg_count",
        "vm_max",
        "bs_max",
        "pg_max",
        "qg_max",
        "vm_mean",
        "bs_mean",
        "pg_mean",
        "qg_mean",
    ]
    info(_LOGGER, join(data, ", "))

    data = [
        "DATA_CHANGES",
        contingency,
        length(network["bus"]),
        length(network["branch"]),
        length(network["gen_contingencies"]),
        length(network["branch_contingencies"]),
        length(vm_changes)-1,
        length(bs_changes)-1,
        length(pg_changes)-1,
        length(qg_changes)-1,
        maximum(vm_changes),
        maximum(bs_changes),
        maximum(pg_changes),
        maximum(qg_changes),
        mean(vm_changes),
        mean(bs_changes),
        mean(pg_changes),
        mean(qg_changes),
    ]
    info(_LOGGER, join(data, ", "))
end

function write_solution(goc_data, network, contingencies; solution_1="", solution_2="")
    error(_LOGGER, "the write_solution function has been replaced by write_solution1, write_solution2")
end

function write_solution1(network; output_dir="", solution_file="solution1.txt")
    if length(output_dir) > 0
        solution_path = joinpath(output_dir, solution_file)
    else
        solution_path = solution_file
    end

    open(solution_path, "w") do sol1
        base_mva = network["baseMVA"]

        bus_switched_shunt_b = Dict(i => 0.0 for (i,bus) in network["bus"])
        for (i,shunt) in network["shunt"]
            # test checks if a "switched shunt" in the orginal data model
            if shunt["dispatchable"] && shunt["status"] == 1
                @assert shunt["gs"] == 0.0
                bus_switched_shunt_b["$(shunt["shunt_bus"])"] += shunt["bs"]
            end
        end

        write(sol1, "-- bus section\n")
        write(sol1, "i, v(p.u.), theta(deg), bcs(MVAR at v = 1 p.u.)\n")
        for (i,bus) in network["bus"]
            write(sol1, "$(bus["index"]), $(bus["vm"]), $(rad2deg(bus["va"])), $(base_mva*bus_switched_shunt_b[i])\n")
        end

        write(sol1, "-- generator section\n")
        write(sol1, "i, id, p(MW), q(MVAR)\n")
        for (i,gen) in network["gen"]
            bus_index = gen["source_id"][2]
            gen_id = gen["source_id"][3]
            write(sol1, "$(bus_index), $(gen_id), $(base_mva*gen["pg"]), $(base_mva*gen["qg"])\n")
        end
    end
end

function write_solution2(network, contingencies; output_dir="", solution_file="solution2.txt")
    if length(output_dir) > 0
        solution_path = joinpath(output_dir, solution_file)
    else
        solution_path = solution_file
    end

    open(solution_path, "w") do sol2
        base_mva = network["baseMVA"]

        for cont_solution in contingencies
            write_solution2_contingency(sol2, network, cont_solution)
        end
    end

    return solution_path
end

function write_solution2_contingency(io::IO, network, contingency_solution)
    base_mva = network["baseMVA"]

    bus_switched_shunt_b = Dict(i => 0.0 for (i,bus) in network["bus"])
    if haskey(network, "shunt") && haskey(contingency_solution, "shunt")
        for (i,nw_shunt) in network["shunt"]
            if nw_shunt["dispatchable"] && nw_shunt["status"] == 1
                #@assert nw_shunt["gs"] == 0.0
                shunt = contingency_solution["shunt"][i]
                bus_switched_shunt_b["$(nw_shunt["shunt_bus"])"] += shunt["bs"]
            end
        end
    end

    write(io, "-- contingency\n")
    write(io, "label\n")
    write(io, "$(contingency_solution["label"])\n")

    write(io, "-- bus section\n")
    write(io, "i, v(p.u.), theta(deg), bcs(MVAR at v = 1 p.u.)\n")
    for (i,nw_bus) in network["bus"]
        vm = 0.0
        va = 0.0

        if haskey(contingency_solution["bus"], i)
            bus = contingency_solution["bus"][i]
            vm = bus["vm"]
            va = bus["va"]
        end

        write(io, "$(nw_bus["index"]), $(vm), $(rad2deg(va)), $(base_mva*bus_switched_shunt_b[i])\n")
    end

    write(io, "-- generator section\n")
    write(io, "i, id, p(MW), q(MVAR)\n")
    for (i,nw_gen) in network["gen"]

        nw_gen = network["gen"][i]
        bus_index = nw_gen["source_id"][2]
        gen_id = nw_gen["source_id"][3]

        pg = 0.0
        qg = 0.0

        if haskey(contingency_solution["gen"], i)
            gen = contingency_solution["gen"][i]
            pg = gen["pg"]
            qg = gen["qg"]
        end

        write(io, "$(bus_index), $(gen_id), $(base_mva*pg), $(base_mva*qg)\n")
    end

    write(io, "-- delta section\n")
    write(io, "delta(MW)\n")
    write(io, "$(base_mva*contingency_solution["delta"])\n")
end

function combine_files(files, output_file_name; output_dir="")
    if length(output_dir) > 0
        output_path = joinpath(output_dir, output_file_name)
    else
        output_path = output_file_name
    end

    open(output_path, "w") do output
        for file in files
            open(file, "r") do input
                for line in readlines(input, keep=true)
                    write(output, line)
                end
            end
        end
    end

    return output_path
end

function remove_files(files)
    for file in files
        if isfile(file)
            info(_LOGGER, "deleting: $(file)")
            rm(file)
        else
            info(_LOGGER, "skipping file: $(file)")
        end
    end
end

function remove_solution_files(;output_dir="", solution1_file="solution1.txt", solution2_file="solution2.txt")
    if length(output_dir) > 0
        solution1_path = joinpath(output_dir, solution1_file)
    else
        solution1_path = solution1_file
    end

    if isfile(solution1_path)
        info(_LOGGER, "deleting: $(solution1_path)")
        rm(solution1_path)
    else
        info(_LOGGER, "skipping file: $(solution1_path)")
    end

    if length(output_dir) > 0
        solution2_path = joinpath(output_dir, solution2_file)
    else
        solution2_path = solution2_file
    end

    if isfile(solution2_path)
        info(_LOGGER, "deleting: $(solution2_path)")
        rm(solution2_path)
    else
        info(_LOGGER, "skipping file: $(solution2_path)")
    end
end

function remove_detail_file(;output_dir="", detail_file="detail.csv")
    if length(output_dir) > 0
        detail_file_path = joinpath(output_dir, detail_file)
    else
        detail_file_path = detail_file
    end

    if isfile(detail_file_path)
        info(_LOGGER, "deleting: $(detail_file_path)")
        rm(detail_file_path)
    else
        info(_LOGGER, "skipping file: $(detail_file_path)")
    end
end

function write_file_paths(files; output_dir="", solution1_file="solution1.txt", solution2_file="solution2.txt")
    if length(output_dir) > 0
        solution1_path = joinpath(output_dir, solution1_file)
    else
        solution1_path = solution1_file
    end

    if length(output_dir) > 0
        solution2_path = joinpath(output_dir, solution2_file)
    else
        solution2_path = solution2_file
    end

    files["sol1"] = solution1_path
    files["sol2"] = solution2_path
    open("files.json", "w") do input
        JSON.print(input, files)
    end
end

function read_solution1(network; output_dir="", state_file="solution1.txt")
    if length(output_dir) > 0
        solution1_path = joinpath(output_dir, state_file)
    else
        solution1_path = state_file
    end

    return build_pm_solution(network, solution1_path)
end

function build_pm_solution(network, goc_sol_file::String)
    info(_LOGGER, "loading solution file: $(goc_sol_file)")
    goc_sol = parse_solution1_file(goc_sol_file)

    info(_LOGGER, "converting GOC solution to PowerModels solution")
    pm_sol = build_pm_solution(network, goc_sol)

    return pm_sol
end

function build_pm_solution(network, goc_sol)
    bus_lookup = Dict(parse(Int, bus["source_id"][2]) => bus for (i,bus) in network["bus"])
    gen_lookup = Dict((gen["source_id"][2], strip(gen["source_id"][3])) => gen for (i,gen) in network["gen"])
    shunt_lookup = Dict{Int,Any}()
    for (i,shunt) in network["shunt"]
        if shunt["source_id"][1] == "switched shunt"
            @assert shunt["source_id"][3] == 0
            shunt_lookup[shunt["source_id"][2]] = shunt
        end
    end

    base_mva = network["baseMVA"]

    bus_data = Dict{String,Any}()
    shunt_data = Dict{String,Any}()
    for bus_sol in goc_sol.bus
        pm_bus = bus_lookup[bus_sol.bus]
        bus_data["$(pm_bus["index"])"] = Dict(
            "vm" => bus_sol.vm,
            "va" => deg2rad(bus_sol.va)
        )

        if haskey(shunt_lookup, bus_sol.bus)
            pm_shunt = shunt_lookup[bus_sol.bus]
            shunt_data["$(pm_shunt["index"])"] = Dict(
                "gs" => 0.0,
                "bs" => bus_sol.bcs/base_mva
            )
        else
            @assert bus_sol.bcs == 0.0
        end
    end

    gen_data = Dict{String,Any}()
    for gen_sol in goc_sol.gen
        pm_gen = gen_lookup[(gen_sol.bus, gen_sol.id)]
        gen_data["$(pm_gen["index"])"] = Dict(
            "pg" => gen_sol.pg/base_mva,
            "qg" => gen_sol.qg/base_mva
        )
    end

    solution = Dict(
        "per_unit" => true,
        "bus" => bus_data,
        "shunt" => shunt_data,
        "gen" => gen_data
    )

    return solution
end


text:
function pressure(model::EoSModel, v::__VolumeKind, T::Unitful.Temperature, z=SA[1.]; output=u"Pa")
    st = standarize(model,v,T,z)
    _v,_T,_z = state_to_vt(model,st)
    res = pressure(model, _v, _T,_z)*u"Pa"
    return uconvert(output, res)
end

for (fn,unit) in Iterators.zip(
    [:enthalpy,
    :entropy,
    :internal_energy,
    :gibbs_free_energy,
    :helmholtz_free_energy,
    :isochoric_heat_capacity,
    :isobaric_heat_capacity,
    :isothermal_compressibility,
    :isentropic_compressibility,
    :speed_of_sound,
    :isobaric_expansivity,
    :joule_thomson_coefficient],
    [u"J",
    u"J/K",
    u"J",
    u"J",
    u"J",
    u"J/K",
    u"J/K",
    u"1/Pa",
    u"1/Pa",
    u"m/s",
    u"1/K",
    u"K/Pa"])
    VT_fn = Symbol(:VT_,fn)
    @eval begin
        function $fn(model::EoSModel, v::__VolumeKind, T::Unitful.Temperature, z=SA[1.]; output=$unit)
            st = standarize(model,v,T,z)
            _v,_T,_z = state_to_vt(model,st)
            res = $VT_fn(model, _v, _T,_z)*$unit
            return uconvert(output, res)
        end

        function $fn(model::EoSModel, p::Unitful.Pressure, T::Unitful.Temperature, z=SA[1.]; phase="unknown", output=$unit)
            st = standarize(model,p,T,z)
            _p,_T,_z = state_to_pt(model,st)
            res = $fn(model, _p, _T, _z; phase=phase)*($unit)
            return uconvert(output, res)
        end
    end
end
function volume(model::EoSModel, p::Unitful.Pressure, T::Unitful.Temperature, z=SA[1.]; phase="unknown", output=u"m^3")
    st = standarize(model,p,T,z)
    _p,_T,_z = state_to_pt(model,st)
    res = volume(model, _p, _T, _z; phase=phase)*u"m^3"
    return uconvert(output, res)
end

#compressibility_factor
function compressibility_factor(model::EoSModel, v::__VolumeKind, T::Unitful.Temperature, z=SA[1.])
    st = standarize(model,v,T,z)
    _v,_T,_z = state_to_vt(model,st)
    res = VT_compressibility_factor(model, _v, _T,_z)
    return res
end

function compressibility_factor(model::EoSModel, p::Unitful.Pressure, T::Unitful.Temperature, z=SA[1.]; phase="unknown")
    st = standarize(model,p,T,z)
    _p,_T,_z = state_to_pt(model,st)
    res = compressibility_factor(model, _p, _T, _z; phase=phase)
    return res
end

#second_virial_coefficient
function second_virial_coefficient(model::EoSModel, T::Unitful.Temperature, z=SA[1.];output=u"m^3")
    st = standarize(model,-1,T,z)
    _,_T,_z = state_to_pt(model,st)
    res = second_virial_coefficient(model, _T,_z)*u"m^3"
    return uconvert(output, res)
end

#pip
function pip(model::EoSModel, v::__VolumeKind, T::Unitful.Temperature, z=SA[1.])
    st = standarize(model,v,T,z)
    _v,_T,_z = state_to_vt(model,st)
    res = pip(model, _v, _T,_z)
    return res
end

#inversion_temperature
function inversion_temperature(model::EoSModel, p::Unitful.Pressure, z=SA[1.]; output=u"K")
    st = standarize(model,p,-1,z)
    _p,_,_z = state_to_pt(model,st)
    res = inversion_temperature(model, _p, _z)*u"K"
    return uconvert(output, res)
end

#enthalpy_vap
function enthalpy_vap(model::EoSModel, T::Unitful.Temperature; output=u"J")
    st = standarize(model,-1,T,SA[1.0])
    _,_T,_ = state_to_pt(model,st)
    res = enthalpy_vap(model, _T)*u"J"
    return uconvert(output,res)
end

function saturation_pressure(model::EoSModel, T::Unitful.Temperature; output=[u"Pa", u"m^3", u"m^3"])
    st = standarize(model,-1,T,SA[1.0])
    _,_T,_ = state_to_pt(model,st)
    (P_sat, v_l, v_v) = saturation_pressure(model,_T)
    _P_sat = uconvert(output[1],P_sat*u"Pa")
    _v_l = uconvert(output[2],v_l*u"m^3")
    _v_v = uconvert(output[3],v_v*u"m^3")
    return (_P_sat,_v_l,_v_v)
end

#molar density
function molar_density(model::EoSModel, p::Unitful.Pressure, T::Unitful.Temperature, z=SA[1.]; phase="unknown",threaded=true, output=u"mol/m^3")
    st = standarize(model,p,T,z)
    _p,_T,_z = state_to_pt(model,st)

    res = molar_density(model,_p,_T,_z;phase=phase,threaded=threaded) *u"mol/m^3"
    return uconvert(output, res)
end

#mass density
function mass_density(model::EoSModel, p::Unitful.Pressure, T::Unitful.Temperature, z=SA[1.]; phase="unknown",threaded=true, output=u"kg/m^3")
    st = standarize(model,p,T,z)
    _p,_T,_z = state_to_pt(model,st)

    res = mass_density(model,_p,_T,_z;phase=phase,threaded=threaded) *u"kg/m^3"
    return uconvert(output, res)
end