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

ext: jl


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


text:
#################### Safe searches ####################

"""
    find_by_binary(tree::GeneralNode, bin::String)

Alias for `find_binary`.
"""
find_by_binary(tree::GeneralNode, bin::String) = find_binary(tree, bin)

"""
    find_by_name(tree::GeneralNode, name::AbstractString)

Alias for `find_name`.
"""
find_by_name(tree::GeneralNode, name::AbstractString) = find_name(tree, name)

"""
    find_by_num(tree::GeneralNode, num::Int64)

Alias for `find_num`.
"""
find_by_num(tree::GeneralNode, num::Int64) = find_num(tree, num)

#################### Unsafe searches ####################

"""
    find_num(root::T, num::Int64)  where T<:AbstractNode

Find a node by its number. The function assumes that the node is present in the
tree.

Do not use this function if you are unsure whether the node is in the tree at all.

Returns reference to Node.

* `root` : root Node of tree to be searched.

* `num` : number of desired Node.
"""
function find_num(root::T, num::I)::T  where {T<:AbstractNode, I<:Integer}
    store = T[]
    found = find_num(root, num, store)
    if length(store) == 0
        throw(ArgumentError("Node not found"))
    else
        return store[1]
    end
end

"""
    find_num(root::T, num::Int64, rn::Vector{T})::Bool  where T<:AbstractNode

Do a post order traversal to find the node corresponding to the `num`.

Returns true if node is found, false otherwise. Desired Node is pushed to rn.

* `root` : root Node of tree to be searched.

* `num` : number of desired Node.

* `rn` : Vector of Nodes; desired Node is pushed to this vector when found.
"""
function find_num(root::T, num::I, rn::Vector{T})::Bool  where {T<:AbstractNode, I<:Integer}
    # if the current node is the correct one store it in rn
    if root.num === num
        push!(rn, root)
        found = true
    else
        found = false
    end

    if !found
        # if the node is not yet found continue
        for child in root.children
            found = find_num(child, num,  rn)
            found && break
        end
    end # if
    return found
end

"""
    find_binary(root::T, bin::String)::T where T<:AbstractNode

Find a node by its binary representation. The function assumes that the node is
present in the tree.

Do not use this function if you are unsure whether the node is in the tree at all.

Returns a reference to the desired Node.

* `root` : root Node of tree to search.

* `bin` : binary representation of desired Node as a String.
"""
function find_binary(root::T, bin::String)::T where T<:AbstractNode
    rv = root
    for i in split(bin[1:end],",")[2:end]
        rv = rv.children[parse(Int64, i)+1]
    end
    rv
end

"""
    find_name(root::T, name::S)::T  where {T<:AbstractNode, S<:AbstractString}

Find a node by its name. Returns reference to Node.

* `root` : root Node of tree to be searched.

* `name` : name of desired Node.
"""
function find_name(root::T, name::S)::T  where {T<:AbstractNode, S<:AbstractString}
    store = T[]
    found = find_name(root, name, store)
    if length(store) == 0
        throw(ArgumentError("Node not found"))
    else
        return store[1]
    end
end

"""
    find_name(root::T, name::S, rn::Vector{T})::Bool where {T<:AbstractNode, S<:AbstractString}

Do a post order traversal to find the node corresponding to the `name`.

Returns true if node is found, false otherwise. Desired Node is pushed to rn.

* `root` : root Node of tree to be searched.

* `name` : name of desired Node.

* `rn` : Vector of Nodes; desired Node is pushed to this vector when found.
"""
function find_name(root::T, name::S, rn::Vector{T})::Bool where {T<:AbstractNode, S<:AbstractString}
    # if the current node is the correct one store it in rn
    if root.name === name
        push!(rn, root)
        found = true
    else
        found = false
    end

    if !found
        # if the node is not yet found continue
        for child in root.children
            found = find_name(child, name,  rn)
            found && break
        end
    end # if
    return found
end

"""
    find_root(node::T)::T where T <: GeneralNode

Finds the root of tree indicated by Node.

Returns reference to root Node of the tree.

* `node` : Node in Tree of interest.
"""
function find_root(node::T)::T where T <: GeneralNode
    while node.root == false
        node = get_mother(node)
    end # while
    return node
end


text:
using LinearAlgebra
using SparseArrays
import Distributions
using Statistics

export GainCache, GainPrior

abstract type RIMEModel <: AbstractModel end

basemodel(m::RIMEModel) = m.model
flux(m::RIMEModel) = flux(vismodel(m))
radialextent(m::RIMEModel) = radialextent(vismodel(m))

struct DesignMatrix{X,M<:AbstractMatrix{X},T,S} <: AbstractMatrix{X}
    matrix::M
    times::T
    stations::S
end

Base.getindex(m::DesignMatrix, i::Int) = getindex(m.matrix, i)
Base.size(m::DesignMatrix) = size(m.matrix)
Base.IndexStyle(::Type{<:DesignMatrix{X,M}}) where {X,M} = Base.IndexStyle(M)
Base.getindex(m::DesignMatrix, I::Vararg{Int,N}) where {N} = getindex(m.matrix, I)
Base.setindex!(m::DesignMatrix, v, i::Int) = setindex!(m.matrix, v, i)
Base.setindex!(m::DesignMatrix, v, i::Vararg{Int, N}) where {N} = setindex!(m.matrix, v, i)

Base.similar(m::DesignMatrix, ::Type{S}, dims::Dims) where {S} = DesignMatrix(similar(m.matrix, S, dims), m.times, m.stations)

function LinearAlgebra.mul!(y::AbstractArray, M::DesignMatrix, x::AbstractArray)
    LinearAlgebra.mul!(y, M.matrix, x)
end

struct GainPrior{S,T,G} <: Distributions.ContinuousMultivariateDistribution
    stations::S
    times::T
    dist::G
end

function GainPrior(dists, st::ScanTable)
    @argcheck Set(keys(dists)) == Set(stations(st))

    gtimes, gstat = gain_stations(st)

    gprior = Dists.product_distribution([getproperty(dists, g) for g in gstat])
    return GainPrior(gstat, gtimes, gprior)
end

HypercubeTransform.bijector(d::GainPrior) = HypercubeTransform.bijector(d.dist)
HypercubeTransform.asflat(d::GainPrior) = asflat(d.dist)
HypercubeTransform.ascube(d::GainPrior) = ascube(d.dist)

Distributions.sampler(d::GainPrior) = Distributions.sampler(d.dist)
Base.length(d::GainPrior) = length(d.dist)
Base.eltype(d::GainPrior) = eltype(d.dist)
function Distributions._rand!(rng::AbstractRNG, d::GainPrior, x::AbstractVector)
    Distributions._rand!(rng, d.dist, x)
end

function Distributions._logpdf(d::GainPrior, x::AbstractArray)
    return Distributions._logpdf(d.dist, x)
end

Statistics.mean(d::GainPrior) = mean(d.dist)
Statistics.var(d::GainPrior) = var(d.dist)
Distributions.entropy(d::GainPrior) = entropy(d.dist)
Statistics.cov(d::GainPrior) = cov(d.dist)

struct GainCache{D1<:DesignMatrix, D2<:DesignMatrix, T, S}
    m1::D1
    m2::D2
    times::T
    stations::S
end

function GainCache(st::ScanTable)
    gtime, gstat = gain_stations(st)
    m1, m2 = gain_design(st)
    return GainCache(m1, m2, gtime, gstat)
end

"""
    `GainModel`
Construct the corruption model for the telescope gains
"""
struct GainModel{C, G, M} <: RIMEModel
    cache::C
    gains::G
    model::M
end

function intensitymap!(img::IntensityMap, model::GainModel, args...)
    return intensitymap!(img, model.model, args...)
end

function visibilities(model::GainModel, u::AbstractArray, v::AbstractArray)
    vis = visibilities(model.model, u, v)
    return corrupt(vis, model.cache, model.gains)
end

function amplitudes(model::GainModel, u::AbstractArray, v::AbstractArray)
    amp = amplitudes(model.model, u, v)
    return abs.(corrupt(amp, model.cache, model.gains))
end

function closure_phases(model::GainModel, args::Vararg{<:AbstractArray, N}) where {N}
    return closure_phases(model.model, args...)
end

function logclosure_amplitudes(model::GainModel, args::Vararg{<:AbstractArray, N}) where {N}
    return logclosure_amplitudes(model.model, args...)
end

function corrupt(vis::AbstractArray, g::GainCache, gains::AbstractArray)
    g1 = g.m1*gains
    g2 = g.m2*gains
    return @. g1*vis*conj(g2)
end

function gain_stations(st::ScanTable)
    gainstat = Symbol[]
    times = eltype(st.times)[]
    for i in 1:length(st)
        s = stations(st[i])
        append!(gainstat, s)
        append!(times, fill(st[i].time, length(s)))
    end
    return times, gainstat
end

function gain_design(st::ScanTable)

    # Construct the indices that will tag where gains are
    rowInd1 = Int[]
    colInd1 = Int[]
    rowInd2 = Int[]
    colInd2 = Int[]
    times = st.obs[:time]
    bl = st.obs[:baseline]
    gaintime, gainstat = gain_stations(st)
    gts = collect(zip(gaintime, gainstat))
    for i in 1:length(times)
        t = times[i]
        s1, s2 = bl[i]
        # now get the columns that corresponds to the gain
        c1 = findall(x->((x[1]==t)&&(x[2]==s1)), gts)
        c2 = findall(x->((x[1]==t)&&(x[2]==s2)), gts)
        append!(colInd1, c1)
        append!(rowInd1, fill(i, length(c1)))
        append!(colInd2, c2)
        append!(rowInd2, fill(i, length(c2)))
    end
    z = fill(1, length(rowInd1))
    m1 = sparse(rowInd1, colInd1, z, length(times), length(gaintime))
    m2 = sparse(rowInd2, colInd2, z, length(times), length(gaintime))
    return DesignMatrix(m1, gaintime, gainstat), DesignMatrix(m2, gaintime, gainstat)
end


text:
function _FDEsolver(pos_args, opt_args, JF::Function, par...)

    # extract arguments from pos_args and opt_args structure fields
    β = pos_args.β

    # check compatibility size of the problem with number of fractional orders
    β_length = length(pos_args.β)
    problem_size = size(pos_args.y0, 1)

    if β_length > 1

        # println("Oh! Good! ODE system!")

    else

        β = pos_args.β * ones(problem_size, 1)
        β_length = problem_size

    end

    # Storage of initial conditions
    ic = initial_conditions(pos_args.tSpan[1], pos_args.y0, Int64.(map(ceil, β)), zeros(β_length, Int64.(ceil(maximum(β)))))

    for i in 1:β_length

        for j in 0:ic.m_β[i] - 1

            ic.m_β_factorial[i, j + 1] = factorial(j)

        end

    end

    # Storage of information on the problem
    Probl = JProblem(ic, pos_args.F, problem_size, par, β, β_length, JF)

    # Time discretization
    N = Int64.(cld(pos_args.tSpan[2] - pos_args.tSpan[1], opt_args.h))
    t = pos_args.tSpan[1] .+ collect(0:N) .* opt_args.h

    # Check compatibility size of the problem with size of the vector field
    f_temp = f_value(pos_args.F(t[1], pos_args.y0[:, 1], par...), Probl.problem_size)

    # Number of points in which to evaluate weights and solution
    r = 16
    Nr::Int64 = ceil((N + 1) / r) * r
    Qr::Int64 = ceil(log2(Nr / r)) - 1
    NNr::Int64 = 2^(Qr + 1) * r

    # Preallocation of some variables
    y = zeros(Probl.problem_size, N + 1)
    fy = zeros(Probl.problem_size, N + 1)
    zn = zeros(Probl.problem_size, NNr + 1)

    # Evaluation of coefficients of the PECE method
    nvett = 0:NNr + 1
    an = zeros(Probl.β_length, NNr + 1)
    a0 = zeros(Probl.β_length, NNr + 1)

    for i_β in 1:Probl.β_length

        find_β = findall(β[i_β] == β[1:i_β - 1])

        if !isempty(find_β) # it is for speeding up the computations; we can use multilpe distpach

            an[i_β, :] = an[find_β[1], :]
            a0[i_β, :] = a0[find_β[1], :]

        else

            nβ = nvett .^ β[i_β]
            nβ1 = nβ .* nvett
            an[i_β,:] = [ 1; (nβ1[1:end - 2] - 2 * nβ1[2:end - 1] + nβ1[3:end])]
            a0[i_β, :] = [0; (nβ1[1:end - 2] - nβ[2:end - 1] .* (nvett[2:end - 1] .- β[i_β] .- 1))]

        end

    end

    METH = JMethod(an, a0, opt_args.h .^ β ./ Γ(β .+ 1), opt_args.h .^ β ./ Γ(β .+ 2), opt_args.nc, opt_args.tol, r, opt_args.StopIt, opt_args.itmax)

    # Evaluation of FFT of coefficients of the PECE method
    if Qr >= 0

        index_fft = Int64.(zeros(2, Qr + 1)) # I have tried index_fft::Int64 = zeros(2,Qr+1) and I got an error for converting Type!

        for l in 1:Qr + 1

            if l == 1

                index_fft[1, l] = 1
                index_fft[2, l] = r * 2

            else

                index_fft[1, l] = index_fft[2, l - 1] + 1
                index_fft[2, l] = index_fft[2, l - 1] + 2^l * r

            end

        end

        an_fft =ComplexF64.(zeros(Probl.β_length, index_fft[2, Qr + 1]))

        for l in 1:Qr + 1

            coef_end = 2^l * r

            for i_β in 1:Probl.β_length

                find_β = findall(β[i_β] == β[1:i_β - 1])

                if !isempty(find_β)

                    an_fft[i_β, index_fft[1, l]:index_fft[2, l]] = an_fft[find_β[1], index_fft[1, l]:index_fft[2, l]]

                else

                    an_fft[i_β, index_fft[1, l]:index_fft[2, l]] = fft(METH.an[i_β, 1:coef_end])

                end

            end

        end

        METH_fft = JMethod_fft(an_fft, Int64.(index_fft))

    end

    # Initializing solution and proces of computation
    y[:, 1] = pos_args.y0[:, 1]
    fy[:, 1] = f_temp
    y, fy = JTriangolo(1, r - 1, t, y, fy, zn, N, METH, Probl)

    # Main process of computation by means of the FFT algorithm
    ff = zeros(1, 2^(Qr + 2))
    ff[1:2] = [0, 2]
    card_ff = 2

    nx0 = 0
    ny0 = 0

    for qr in 0:Qr

        L = 2^qr
        y, fy = JDisegnaBlocchi(L, ff, r, Nr, nx0 + L * r, ny0, t, y, fy, zn, N, METH, METH_fft, Probl)
        ff[1:2 * card_ff] = [ff[1:card_ff]; ff[1:card_ff]]
        card_ff = 2 * card_ff
        ff[card_ff] = 4 * L

    end

    # Evaluation solution in T when T is not in the mesh
    if pos_args.tSpan[2] < t[N + 1]

        c = [pos_args.tSpan[2] - t(N)] / opt_args.h
        t[N + 1] = pos_args.tSpan[2]
        y[:, N + 1] = (1 - c) * y[:, N] + c * y[:, N + 1]

    end

    t = t[1:N + 1]
    y = y[:, 1:N + 1]'

    return t, y

end


text:
## File containing different utility functions ##

## Jittering object adapting to the type of the GP ##
struct Jittering end;

const jitt = Jittering()

@inline Base.Float64(::Jittering) = Float64(1e-4)
@inline Base.Float32(::Jittering) = Float32(1e-3)
@inline Base.Float16(::Jittering) = Float16(1e-2)
@inline Base.convert(::Type{Float64}, ::Jittering) = Float64(1e-4)
@inline Base.convert(::Type{Float32}, ::Jittering) = Float32(1e-3)
@inline Base.convert(::Type{Float16}, ::Jittering) = Float16(1e-2)

# Compute expectation given the nodes
function expectation(f, μ::Real, σ²::Real)
    x = pred_nodes * sqrt(max(σ², zero(σ²))) .+ μ
    return dot(pred_weights, f.(x))
end

## delta function `(i,j)`, equal `1` if `i == j`, `0` else ##
@inline function δ(T, i::Integer, j::Integer)
    return ifelse(i == j, one(T), zero(T))
end
δ(i::Integer, j::Integer) = δ(Float64, i, j)

## Hadamard product between two arrays of same size ##
@inline function hadamard(A::AbstractArray{<:Real}, B::AbstractArray{<:Real})
    return A .* B
end

## Add on place the transpose of aZero Mean Prior matrix ##
@inline function add_transpose!(A::AbstractMatrix{<:Real})
    return A .+= A'
end

invquad(a::Cholesky, x::AbstractVector) = sum(abs2, a.L \ x)

## Return the trace of A*B' ##
@inline function trace_ABt(A::AbstractMatrix{<:Real}, B::AbstractMatrix{<:Real})
    return dot(A, B)
end

## Return the diagonal of A*B' ##
@inline function diag_ABt(A::AbstractMatrix, B::AbstractMatrix)
    return vec(sum(A .* B; dims=2))
end

## Return the multiplication of Diagonal(v)*B ##
function diagv_B(v::AbstractVector, B::AbstractMatrix)
    return v .* B
end

## Return the multiplication of κᵀ*Diagonal(θ)*κ
@inline function κdiagθκ(κ::AbstractMatrix{T}, θ::AbstractVector{T}) where {T<:Real}
    return transpose(θ .* κ) * κ
end

## Return the multiplication of ρ*κᵀ*Diagonal(θ)*κ
@inline function ρκdiagθκ(ρ::T, κ::AbstractMatrix{T}, θ::AbstractVector{T}) where {T<:Real}
    return transpose((ρ * θ) .* κ) * κ
end

## Return the addition of a diagonal to a symmetric matrix ##
function opt_add_diag_mat(v::AbstractVector{T}, B::AbstractMatrix{T}) where {T<:Real}
    A = copy(B)
    @inbounds for i in 1:size(A, 1)
        A[i, i] += v[i]
    end
    return A
end

## Compute exp(μ)/cosh(c) safely if there is an overflow ##
function safe_expcosh(μ::Real, c::Real)
    return isfinite(exp(μ) / cosh(c)) ? exp(μ) / cosh(c) : 2 * logistic(2.0 * max(μ, c))
end

## Return a safe version of log(cosh(c)) ##
function logcosh(c::Real)
    return log(exp(-2.0 * c) + 1.0) + c - logtwo
end

function symcat(S::Symmetric, v::AbstractVector, vv::Real)
    S = vcat(S, v')
    S = hcat(S, vcat(v, vv))
    return Symmetric(S)
end

export make_grid
function make_grid(range1, range2)
    return hcat([i for i in range1, j in range2][:], [j for i in range1, j in range2][:])
end

Base.:*(x::Real, C::Cholesky) = Cholesky(sqrt(x) * C.factors, C.uplo, C.info)

Base.:*(C::Cholesky, x::AbstractVecOrMat) = C.L * (transpose(C.L) * x)

Base.:+(C::Cholesky, x::UniformScaling) = cholesky(Symmetric(Array(C) + x))


text:
function solve(case, N, X, A)
    AA = Channel{Tuple{Int,Int}}(N)
    foreach(x -> put!(AA, x), enumerate(A[]))
    dequeued = []
    while !isempty(AA)
        (i, person) = take!(AA)
        person -= X
        if person <= 0
            push!(dequeued, i)
        else
            put!(AA, (i, person))
        end
    end

    println("Case #$case: $(join(dequeued, " "))")
end

function main()
    T = parse(Int, readline(stdin))

    for case in 1:T
        N, X = parse.(Int, stdin |> readline |> split)
        A = parse.(Int, stdin |> readline |> split)
        solve(case, N, X, Ref(A))
    end
end

isinteractive() || main()

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using CQL
using Base.Test

include("examples.jl")
# write your own tests here
@test cqltest01()
@test cqltest02()
@test cqltest03()


text:
##############################################################################
##
## We use SubDataTable's to maintain a reference to a subset of a DataTable
## without making copies.
##
##############################################################################

if VERSION >= v"0.6.0-dev.2643"
    include_string("""
        immutable SubDataTable{T <: AbstractVector{Int}} <: AbstractDataTable
            parent::DataTable
            rows::T # maps from subdt row indexes to parent row indexes

            function SubDataTable{T}(parent::DataTable, rows::T) where {T <: AbstractVector{Int}}
                if length(rows) > 0
                    rmin, rmax = extrema(rows)
                    if rmin < 1 || rmax > size(parent, 1)
                        throw(BoundsError())
                    end
                end
                new(parent, rows)
            end
        end
    """)
else
    @eval begin
        immutable SubDataTable{T <: AbstractVector{Int}} <: AbstractDataTable
            parent::DataTable
            rows::T # maps from subdt row indexes to parent row indexes

            function SubDataTable(parent::DataTable, rows::T)
                if length(rows) > 0
                    rmin, rmax = extrema(rows)
                    if rmin < 1 || rmax > size(parent, 1)
                        throw(BoundsError())
                    end
                end
                new(parent, rows)
            end
        end
    end
end

"""
A view of row subsets of an AbstractDataTable

A `SubDataTable` is meant to be constructed with `view`.  A
SubDataTable is used frequently in split/apply sorts of operations.

```julia
view(d::AbstractDataTable, rows)
```

### Arguments

* `d` : an AbstractDataTable
* `rows` : any indexing type for rows, typically an Int,
  AbstractVector{Int}, AbstractVector{Bool}, or a Range

### Notes

A `SubDataTable` is an AbstractDataTable, so expect that most
DataTable functions should work. Such methods include `describe`,
`dump`, `nrow`, `size`, `by`, `stack`, and `join`. Indexing is just
like a DataTable; copies are returned.

To subset along columns, use standard column indexing as that creates
a view to the columns by default. To subset along rows and columns,
use column-based indexing with `view`.

### Examples

```julia
dt = DataTable(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
sdt1 = view(dt, 1:6)
sdt2 = view(dt, dt[:a] .> 1)
sdt3 = view(dt[[1,3]], dt[:a] .> 1)  # row and column subsetting
sdt4 = groupby(dt, :a)[1]  # indexing a GroupedDataTable returns a SubDataTable
sdt5 = view(sdt1, 1:3)
sdt1[:,[:a,:b]]
```

"""
SubDataTable

function SubDataTable{T <: AbstractVector{Int}}(parent::DataTable, rows::T)
    return SubDataTable{T}(parent, rows)
end

function SubDataTable(parent::DataTable, row::Integer)
    return SubDataTable(parent, [Int(row)])
end

function SubDataTable{S <: Integer}(parent::DataTable, rows::AbstractVector{S})
    return SubDataTable(parent, convert(Vector{Int}, rows))
end

function SubDataTable(parent::DataTable, rows::AbstractVector{Bool})
    return SubDataTable(parent, find(rows))
end

function SubDataTable{T<:Integer}(sdt::SubDataTable, rowinds::Union{T, AbstractVector{T}})
    return SubDataTable(sdt.parent, sdt.rows[rowinds])
end

function Base.view{T<:Nullable}(adt::AbstractDataTable, rowinds::AbstractVector{T})
    # Vector{<:Nullable} need to be checked for nulls and the values lifted
    any(isnull, rowinds) && throw(NullException())
    return SubDataTable(adt, get.(rowinds))
end

function Base.view(adt::AbstractDataTable, rowinds::NullableVector)
    # convert for NullableVectors will throw NullException if nulls present
    return SubDataTable(adt, convert(Vector, rowinds))
end

function Base.view(adt::AbstractDataTable, rowinds::Any)
    return SubDataTable(adt, rowinds)
end

function Base.view(adt::AbstractDataTable, rowinds::Any, colinds::AbstractVector)
    return SubDataTable(adt[colinds], rowinds)
end

function Base.view(adt::AbstractDataTable, rowinds::Any, colinds::Any)
    return SubDataTable(adt[[colinds]], rowinds)
end

##############################################################################
##
## AbstractDataTable interface
##
##############################################################################

index(sdt::SubDataTable) = index(sdt.parent)

# TODO: Remove these
nrow(sdt::SubDataTable) = ncol(sdt) > 0 ? length(sdt.rows)::Int : 0
ncol(sdt::SubDataTable) = length(index(sdt))

function Base.getindex(sdt::SubDataTable, colinds::Any)
    return sdt.parent[sdt.rows, colinds]
end

function Base.getindex(sdt::SubDataTable, rowinds::Any, colinds::Any)
    return sdt.parent[sdt.rows[rowinds], colinds]
end

function Base.setindex!(sdt::SubDataTable, val::Any, colinds::Any)
    sdt.parent[sdt.rows, colinds] = val
    return sdt
end

function Base.setindex!(sdt::SubDataTable, val::Any, rowinds::Any, colinds::Any)
    sdt.parent[sdt.rows[rowinds], colinds] = val
    return sdt
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.map(f::Function, sdt::SubDataTable) = f(sdt) # TODO: deprecate

without(sdt::SubDataTable, c) = view(without(sdt.parent, c), sdt.rows)


text:
using DifferentialEquations
using LinearAlgebra
using Plots

struct Wave
    f
    g
    c
end

# Exercise

# Visualization

function plot_wave(y, file_name; fps = 60, kwargs...)
    anim = @animate for (i, y_row) in enumerate(eachrow(y))
        plot(
            y_row;
            title = "t = $(i-1)Δt",
            xlabel = "x",
            ylabel = "y(t, x)",
            legend = false,
            linewidth = 2,
            kwargs...
        )
    end
    gif(anim, file_name; fps, show_msg = false)

    return nothing
end

# Exercise

# DiferentialEquations

f(u,p,t) = 0.98*u

u0 = 1.0
tspan = (0.0, 1.0)

prob = ODEProblem(f, u0, tspan)

sol = solve(prob)

plot(sol; label="")

sol(0.8)

# Exercise

# Lorenz system

function lorenz(u, p, t)
    σ, ρ, β = p
    x_t = σ*(u[2]-u[1])
    y_t = u[1]*(ρ-u[3]) - u[2]
    z_t = u[1]*u[2] - β*u[3]
    return [x_t; y_t; z_t]
end

u0 = [1.0; 0.0; 0.0]
p = [10; 28; 8/3]

tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan, p)

sol = solve(prob)

plot(sol)

plt1 = plot(sol, vars=(1,2,3), label="")

plot(sol, vars=(1,2,3), denseplot=false; label="")

traj = hcat(sol.u...)
plot(traj[1,:], traj[2,:], traj[3,:]; label="")

# Exercise

# Solution comparison

hcat(sol(tspan[2]), sol0(tspan[2]))

# Bisection setup

function bisection(f, a, b; tol=1e-6)
    fa = f(a)
    fb = f(b)
    fa == 0 && return a
    fb == 0 && return b
    fa*fb > 0 && error("Wrong initial values for bisection")
    while b-a > tol
        c = (a+b)/2
        fc = f(c)
        fc == 0 && return c
        if fa*fc > 0
            a = c
            fa = fc
        else
            b = c
            fb = fc
        end
    end
    return (a+b)/2
end

# PMSM structure

struct PMSM{T<:Real}
    ρ::T
    ω::T
    A::Matrix{T}
    invA::Matrix{T}

    function PMSM(ρ, ω)
        A = -ρ*[1 0; 0 1] -ω*[0 -1; 1 0]
        return new{eltype(A)}(ρ, ω, A, inv(A))
    end
end

function expA(p::PMSM, t)
    ρ, ω = p.ρ, p.ω
    return exp(-ρ*t)*[cos(ω*t) sin(ω*t); -sin(ω*t) cos(ω*t)]
end

ρ = 0.1
ω = 2
x0 = [0;-0.5]
q = [1;0]

ps = PMSM(ρ, ω)

# Exercise

# Exercise

# Plotting trajectory

ts = 0:0.0001:10

xs1 = trajectory_fin_diff(ps, x0, ts, q)
xs2 = trajectory_exact(ps, x0, ts, q)

plot(xs1[1,:], xs1[2,:], label="Finite differences")
plot!(xs2[1,:], xs2[2,:], label="True value")

# Exercise

# Plotting trajectory

function trajectory_control(p::PMSM, x0, ts, q, U_max, p0)
    xs = zeros(length(x0), length(ts))

    for (i, t) in enumerate(ts)
        eAt  = expA(p, t)
        emAt = expA(p, -t)
        xs[:, i] = eAt*(x0 + p.invA * (I - emAt)*q + U_max/ρ*(exp(ρ*t) - 1)*p0)
    end
    return xs
end

p0 = ps.ρ/(U_max*(exp(ps.ρ*τ)-1))*(expA(ps, -τ)*x_t - x0 - ps.invA*(I-expA(ps, -τ))*q)
p0 /= norm(p0)

ts = range(0, τ; length=100)

traj = trajectory_control(ps, x0, ts, q, U_max, p0)

plot(traj[1,:], traj[2,:], label="Optimal trajectory")
scatter!([x0[1]], [x0[2]], label="Starting point")
scatter!([x_t[1]], [x_t[2]], label="Target point")

# BONUS

ts = 0:0.01:10

plt = plot()
for α = 0:π/4:2*π
    trj = trajectory_control(ps, x0, ts, q, U_max, [sin(α); cos(α)])
    plot!(plt, trj[1,:], trj[2,:], label="")
end
display(plt)

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse