text:
@everywhere using JuMP

@everywhere using PowerModels

@everywhere using InfrastructureModels

@everywhere using Memento

@everywhere const LOGGER = getlogger("GOC")

@everywhere const vm_eq_tol = 1e-4

@everywhere import Statistics: mean

@everywhere function build_pm_model(goc_data) # Function to build a model based on the input 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"]) # Dic

    branch_lookup = Dict()
    for (i,branch) in network["branch"]
        if !branch["transformer"] # Checks if the branch is not the 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 macro allows to optionally specify error message, instead of just printing the failed expression
            @assert branch["source_id"][6] == 0
        end
        branch_lookup[branch_id] = branch
    end

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

    @assert network["per_unit"] # Assert macro allows to optionally specify error message, instead of just printing the failed expression
    mva_base = network["baseMVA"] # Assigns baseMVA from network as mva_base

    dispatch_tbl_lookup = Dict() # Creates a dictionary named dispatch_tbl_lookup
    for dispatch_tbl in goc_data.cost["disptbl"]
        dispatch_tbl_lookup[dispatch_tbl["ctbl"]] = dispatch_tbl # Assigns dispatch_tbl to dictionary named dispatch_tbl_lookup
    end

    cost_tbl_lookup = Dict() # Creates a dictionary named cost_tbl_lookup
    for cost_tbl in goc_data.cost["ctbl"]
        cost_tbl_lookup[cost_tbl["ltbl"]] = cost_tbl # Assigns cost_tbl to dictionary named cost_tbl_lookup
    end

    gen_cost_models = Dict() # Defines an empty dictionary for generation cost models
    for gen_dispatch in goc_data.cost["gen"] # Looks for the cost of gen in goc_data
        gen_id = (gen_dispatch["bus"], strip(gen_dispatch["genid"])) # Removes the leading and trailing characters from genid
        dispatch_tbl = dispatch_tbl_lookup[gen_dispatch["disptbl"]] # Calls the function named dispatch_tbl_lookup for gen dispatch table
        cost_tbl = cost_tbl_lookup[dispatch_tbl["ctbl"]] # Calls the function named cost_tbl_lookup for dispatch table cost

        gen_cost_models[gen_id] = cost_tbl # Assigns the cost_tbl to the specific gen cost model id
    end

    if length(gen_cost_models) != length(network["gen"]) # Checks if the length of both generator network and generator cost model is same
        error(LOGGER, "cost model data missing, network has $(length(network["gen"])) generators, the cost model has $(length(gen_cost_models)) generators") # Mismatch error
    end

    for (gen_id, cost_model) in gen_cost_models
        pm_gen = gen_lookup[gen_id] # Calls gen_lookup for gen_id & stores the result in pm_gen
        pm_gen["model"] = 1
        pm_gen["model_lable"] = cost_model["label"] # Assigns label for cost model to model lable for pm_gen
        pm_gen["ncost"] = length(cost_model["points"]) # Finds the length of points in cost_model

        #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) # Inserts point.y at the end of point_list
        end
        pm_gen["cost"] = point_list # Assigns point_list to cost in pm_gen
    end

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

    if length(goc_data.response) != length(network["gen"]) # Checks if all the generators are active or not
        error(LOGGER, "generator response model data missing, network has $(length(network["gen"])) generators, the response model has $(length(goc_data.response)) generators") # Error to mention the out of service generators
    end

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

        pm_gen = gen_lookup[gen_id]

        pm_gen["alpha"] = gen_response["r"]
    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] # bmin = bmin + nt[n_name]*shunt[b_name]
                else
                    bmax += shunt[n_name]*shunt[b_name] # bmax = bmax + nt[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" # Checks if there is branch contingency
            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")) # Inserts idx, label, and type at the end of branch_ids

        elseif cont["component"] == "generator" # Checks if there is generator contingency
            gen_id = (cont["i"], cont["id"])
            pm_gen = gen_lookup[gen_id]
            push!(generator_ids, (idx=pm_gen["index"], label=cont["label"], type="gen")) # Inserts idx, label, and type at the end of generator_ids

        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 #####

    PowerModels.correct_cost_functions!(network)

    # FYI, this breaks output API
    #PowerModels.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

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"] # Checks if branch is other than 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() # Creates a dictionary for dispatch_tbl_lookup
    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"]) # Checks if the length of generator networks is not equal to generator cost model
        error(LOGGER, "cost model data missing, network has $(length(network["gen"])) generators, the cost model has $(length(gen_cost_models)) generators") # Mentions there is error due to mismatch
    end

    for (gen_id, cost_model) in gen_cost_models
        pm_gen = gen_lookup[gen_id]
        pm_gen["model"] = 1
        pm_gen["model_lable"] = 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

    PowerModels.correct_cost_functions!(network)

    return network
end

@everywhere function read_solution1(network; output_dir="", state_file="solution1.txt") # Function to read solution 1 of base case
    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) # Calls the functin named build_pm_solution to initiate the parsing of solution1.txt
end

@everywhere function build_pm_solution(network, goc_sol_file::String) # build_pm_solution function
    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

@everywhere 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

@everywhere function gens_by_bus(network) # Function to see which are generation buses
    bus_gens = Dict(i => Any[] for (i,bus) in network["bus"])
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            push!(bus_gens["$(gen["gen_bus"])"], gen)
        end
    end
    return bus_gens
end

@everywhere function run_fixpoint_pf_v2_2!(network, pg_lost, model_constructor, solver; iteration_limit=typemax(Int64)) # Core algorithm
    time_start = time()

    delta = apply_pg_response!(network, pg_lost)
    #delta = 0.0

    debug(LOGGER, "pg lost: $(pg_lost)")
    debug(LOGGER, "delta: $(network["delta"])")
    debug(LOGGER, "pre-solve time: $(time() - time_start)")

    base_solution = extract_solution(network)
    base_solution["delta"] = delta
    final_result = Dict(
        "termination_status" => LOCALLY_SOLVED,
        "solution" => base_solution
    )

    bus_gens = gens_by_bus(network)

    #network["delta_start"] = network["delta"]
    for (i,gen) in network["gen"]
        gen["qg_fixed"] = false # Set to false to allow variations in qg_fixed
        gen["pg_start"] = gen["pg"]
        if isapprox(gen["qmin"],gen["qmax"])
            gen["qg_fixed"] = true
            gen["qg"] = gen["qmin"]
        end
        gen["qg_start"] = gen["qg"]
    end

    for (i,bus) in network["bus"]
        active_gens = [gen for gen in bus_gens[i] if !gen["qg_fixed"]]
        if length(active_gens) == 0
            bus["vm_fixed"] = false
        else
            bus["vm_fixed"] = true
        end
        #bus["vm_start"] = bus["vm"]
        #bus["va_start"] = bus["va"]
        bus["vr_start"] = bus["vm"]*cos(bus["va"])
        bus["vi_start"] = bus["vm"]*sin(bus["va"])
    end

    time_start = time()
    result = run_fixed_pf_nbf_rect2(network, model_constructor, solver)
    debug(LOGGER, "pf solve time: $(time() - time_start)")
    if result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED
        correct_qg!(network, result["solution"], bus_gens=bus_gens)
        PowerModels.update_data!(network, result["solution"])
        final_result = result
    else
        warn(LOGGER, "contingency pf solver FAILED with status $(result["termination_status"]) on iteration 0") # Warning if solver is failed to resolve the contingency
        return final_result
    end

    pg_switched = true
    qg_switched = true
    vm_switched = true

    iteration = 1
    deltas = [result["solution"]["delta"]]
    while (pg_switched || qg_switched || vm_switched) && iteration <= iteration_limit
        debug(LOGGER, "obj: $(result["objective"])")
        debug(LOGGER, "delta: $(result["solution"]["delta"])")
        pg_switched = false
        qg_switched = false
        vm_switched = false

        for (i,gen) in network["gen"]
            pg = gen["pg_base"] + network["delta"]*gen["alpha"] # To find the real output power of a responding generator

            if gen["pg_fixed"]
                if !isapprox(gen["pmax"], gen["pmin"]) && pg < gen["pmax"] && pg > gen["pmin"]
                    gen["pg"] = pg
                    gen["pg_fixed"] = false
                    pg_switched = true
                    #info(LOGGER, "unfix pg on gen $(i)")
                end
            else
                if pg >= gen["pmax"]
                    gen["pg"] = gen["pmax"]
                    gen["pg_fixed"] = true
                    pg_switched = true
                    #info(LOGGER, "fix pg to ub on gen $(i)")
                elseif gen["pg"] <= gen["pmin"]
                    gen["pg"] = gen["pmin"]
                    gen["pg_fixed"] = true
                    pg_switched = true
                    #info(LOGGER, "fix pg to lb on gen $(i)")
                end
            end
        end

        for (i,bus) in network["bus"]
            if length(bus_gens[i]) > 0
                qg = sum(gen["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(qmin,qmax)
                    @assert !bus["vm_fixed"]
                    for gen in bus_gens[i]
                        @assert gen["qg_fixed"]
                        @assert isapprox(gen["qg"],gen["qmin"])
                    end
                elseif bus["vm_fixed"]
                    if qg >= qmax
                        bus["vm_fixed"] = false
                        vm_switched = true
                        for gen in bus_gens[i]
                            gen["qg"] = gen["qmax"]
                            gen["qg_fixed"] = true
                        end
                        #info(LOGGER, "fix qg to ub on bus $(i)")
                    end

                    if qg <= qmin
                        bus["vm_fixed"] = false
                        vm_switched = true
                        for gen in bus_gens[i]
                            gen["qg"] = gen["qmin"]
                            gen["qg_fixed"] = true
                        end
                        #info(LOGGER, "fix qg to lb on bus $(i)")
                    end
                else
                    if qg < qmax && qg > qmin
                        bus["vm_fixed"] = true

                        vm_switched = true
                        for gen in bus_gens[i]
                            gen["qg_fixed"] = false
                            gen["qg_start"] = gen["qg"]
                        end
                        #info(LOGGER, "fix vm to on bus $(i)")
                    end
                    if qg >= qmax && bus["vm"] > bus["vm_base"]
                        bus["vm_fixed"] = true
                        vm_switched = true
                        for gen in bus_gens[i]
                            gen["qg_fixed"] = false
                        end
                        #info(LOGGER, "fix vm to on bus $(i)")
                    end
                    if qg <= qmin && bus["vm"] < bus["vm_base"]
                        bus["vm_fixed"] = true
                        vm_switched = true
                        for gen in bus_gens[i]
                            gen["qg_fixed"] = false
                        end
                        #info(LOGGER, "fix vm to on bus $(i)")
                    end
                end
            end
        end

        for (i,gen) in network["gen"]
            gen["pg_start"] = gen["pg"]
            gen["qg_start"] = gen["qg"]
        end

        for (i,bus) in network["bus"]
            bus["vm_start"] = bus["vm"]
            bus["va_start"] = bus["va"]
        end

        if pg_switched || qg_switched || vm_switched # || means OR
            debug(LOGGER, "bus or gen swtiched: $iteration")
            time_start = time()
            result = run_fixed_pf_nbf_rect2(network, model_constructor, solver)
            debug(LOGGER, "pf solve time: $(time() - time_start)")
            if result["termination_status"] == LOCALLY_SOLVED || result["termination_status"] == ALMOST_LOCALLY_SOLVED
                correct_qg!(network, result["solution"], bus_gens=bus_gens)
                PowerModels.update_data!(network, result["solution"])
                final_result = result
            else
                warn(LOGGER, "contingency pf solver FAILED with status $(result["termination_status"]) on iteration 0")
                return final_result
            end

            push!(deltas, result["solution"]["delta"])
            iteration += 1
            if iteration >= iteration_limit
                warn(LOGGER, "hit iteration limit")
            end
            if length(deltas) > 3 && isapprox(deltas[end-2], deltas[end])
                warn(LOGGER, "cycle detected, stopping")
                break
            end
        end
    end

    return final_result
end

@everywhere function apply_pg_response!(network, pg_delta) # Function to compute good starting guess for delta
    for (i,gen) in network["gen"]
        gen["pg_fixed"] = false
    end

    pg_total = 0.0
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            pg_total += gen["pg"]
        end
    end

    pg_target = pg_total + pg_delta
    #info(LOGGER, "total gen:  $(pg_total)")
    #info(LOGGER, "target gen: $(pg_target)")

    delta_est = 0.0
    while !isapprox(pg_total, pg_target)
        alpha_total = 0.0
        for (i,gen) in network["gen"]
            if gen["gen_status"] != 0 && !gen["pg_fixed"]
                alpha_total += gen["alpha"]
            end
        end
        #info(LOGGER, "alpha total: $(alpha_total)")

        if isapprox(alpha_total, 0.0) && !isapprox(pg_total, pg_target)
            warn(LOGGER, "insufficient generator response to meet demand, remaining pg $(pg_total - pg_target), remaining alpha $(alpha_total)")
            break
        end

        delta_est += pg_delta/alpha_total
        #info(LOGGER, "detla: $(delta_est)")

        for (i,gen) in network["gen"]
            if gen["gen_status"] != 0
                pg_cont = gen["pg_base"] + delta_est*gen["alpha"]

                if pg_cont <= gen["pmin"]
                    gen["pg"] = gen["pmin"]
                    if !gen["pg_fixed"]
                        gen["pg_fixed"] = true
                        #info(LOGGER, "gen $(i) hit lb $(gen["pmin"]) with target value of $(pg_cont)")
                    end
                elseif pg_cont >= gen["pmax"]
                    gen["pg"] = gen["pmax"]
                    if !gen["pg_fixed"]
                        gen["pg_fixed"] = true
                        #info(LOGGER, "gen $(i) hit ub $(gen["pmax"]) with target value of $(pg_cont)")
                    end
                else
                    gen["pg"] = pg_cont
                end
            end
        end

        pg_total = 0.0
        for (i,gen) in network["gen"]
            if gen["gen_status"] != 0
                pg_total += gen["pg"]
            end
        end
        #pg_comp = comp_pg_response_total(network, delta_est)
        #info(LOGGER, "detla: $(delta_est)")
        #info(LOGGER, "total gen comp $(pg_comp) - gen inc $(pg_total)")
        #info(LOGGER, "total gen $(pg_total) - target gen $(pg_target)")

        pg_delta = pg_target - pg_total
    end

    network["delta"] = delta_est
    return delta_est
end

@everywhere function extract_solution(network; branch_flow=false)
    sol = Dict{String,Any}()

    sol["bus"] = Dict{String,Any}()
    for (i,bus) in network["bus"]
        bus_dict = Dict{String,Any}()
        bus_dict["va"] = get(bus, "va", 0.0)
        bus_dict["vm"] = get(bus, "vm", 1.0)
        sol["bus"][i] = bus_dict
    end

    sol["shunt"] = Dict{String,Any}()
    for (i,shunt) in network["shunt"]
        shunt_dict = Dict{String,Any}()
        shunt_dict["gs"] = get(shunt, "gs", 0.0)
        shunt_dict["bs"] = get(shunt, "bs", 0.0)
        sol["shunt"][i] = shunt_dict
    end

    sol["gen"] = Dict{String,Any}()
    for (i,gen) in network["gen"]
        gen_dict = Dict{String,Any}()
        gen_dict["pg"] = get(gen, "pg", 0.0)
        gen_dict["qg"] = get(gen, "qg", 0.0)
        sol["gen"][i] = gen_dict
    end

    if branch_flow
        sol["branch"] = Dict{String,Any}()
        for (i,branch) in network["branch"]
            branch_dict = Dict{String,Any}()
            branch_dict["pf"] = get(branch, "pf", 0.0)
            branch_dict["qf"] = get(branch, "qf", 0.0)
            branch_dict["pt"] = get(branch, "pt", 0.0)
            branch_dict["qt"] = get(branch, "qt", 0.0)
            sol["branch"][i] = branch_dict
        end
    end

    return sol
end

""

@everywhere function run_fixed_pf_nbf_rect2(file, model_constructor, solver; kwargs...) # Calls function to check out-of-bound values
    return run_model(file, model_constructor, solver, post_fixed_pf_nbf_rect2; solution_builder = solution_second_stage!, kwargs...)
end

""

@everywhere function post_fixed_pf_nbf_rect2(pm::GenericPowerModel) # Function to look for out-of-bound calues after initial pf.
    start_time = time()
    PowerModels.variable_voltage(pm, bounded=false)
    PowerModels.variable_active_generation(pm, bounded=false)
    PowerModels.variable_reactive_generation(pm, bounded=false)
    #PowerModels.variable_branch_flow(pm, bounded=false)
    #PowerModels.variable_dcline_flow(pm, bounded=false)

    #PowerModels.variable_branch_flow(pm, bounded=false)

    # TODO set bounds bounds on alpha and total gen capacity
    var(pm)[:delta] = @variable(pm.model, delta, base_name="delta", start=0.0)
    #var(pm)[:delta] = @variable(pm.model, delta, base_name="delta", start=ref(pm, :delta_start))
    #Memento.info(LOGGER, "post variable time: $(time() - start_time)")

    start_time = time()

    vr = var(pm, :vr)
    vi = var(pm, :vi)

    PowerModels.constraint_model_voltage(pm)

    for (i,bus) in ref(pm, :bus)
        if bus["vm_fixed"]
            @constraint(pm.model, vr[i]^2 + vi[i]^2 == bus["vm_base"]^2)
        end
    end

    for i in ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end
    #Memento.info(LOGGER, "misc constraints time: $(time() - start_time)")

    start_time = time()
    p = Dict{Tuple{Int64,Int64,Int64},GenericQuadExpr{Float64,VariableRef}}()
    q = Dict{Tuple{Int64,Int64,Int64},GenericQuadExpr{Float64,VariableRef}}()
    for (i,branch) in ref(pm, :branch)
        #PowerModels.constraint_ohms_yt_from(pm, i)
        #PowerModels.constraint_ohms_yt_to(pm, i)

        f_bus_id = branch["f_bus"]
        t_bus_id = branch["t_bus"]
        f_idx = (i, f_bus_id, t_bus_id)
        t_idx = (i, t_bus_id, f_bus_id)

        f_bus = ref(pm, :bus, f_bus_id)
        t_bus = ref(pm, :bus, t_bus_id)

        #g, b = PowerModels.calc_branch_y(branch)
        #tr, ti = PowerModels.calc_branch_t(branch)
        g = branch["g"]
        b = branch["b"]
        tr = branch["tr"]
        ti = branch["ti"]

        g_fr = branch["g_fr"]
        b_fr = branch["b_fr"]
        g_to = branch["g_to"]
        b_to = branch["b_to"]
        tm = branch["tap"]

        #p_fr  = var(pm, :p, f_idx)
        #q_fr  = var(pm, :q, f_idx)
        #p_to  = var(pm, :p, t_idx)
        #q_to  = var(pm, :q, t_idx)
        vr_fr = vr[f_bus_id] #var(pm, :vr, f_bus_id)
        vr_to = vr[t_bus_id] #var(pm, :vr, t_bus_id)
        vi_fr = vi[f_bus_id] #var(pm, :vi, f_bus_id)
        vi_to = vi[t_bus_id] #var(pm, :vi, t_bus_id)

        p[f_idx] = (g+g_fr)/tm^2*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
        q[f_idx] = -(b+b_fr)/tm^2*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
        p[t_idx] = (g+g_to)*(vr_to^2 + vi_to^2) + (-g*tr-b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr+g*ti)/tm^2*(-(vi_fr*vr_to - vr_fr*vi_to))
        q[t_idx] = -(b+b_to)*(vr_to^2 + vi_to^2) - (-b*tr+g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr-b*ti)/tm^2*(-(vi_fr*vr_to - vr_fr*vi_to))
    end
    #Memento.info(LOGGER, "flow expr time: $(time() - start_time)")

    start_time = time()
    pg = Dict{Int,Any}()
    qg = Dict{Int,Any}()
    for (i,gen) in ref(pm, :gen)
        if gen["pg_fixed"]
            pg[i] = gen["pg"]
            @constraint(pm.model, var(pm, :pg, i) == gen["pg"])
        else
            pg[i] = gen["pg_base"] + gen["alpha"]*delta
            @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"] + gen["alpha"]*delta)
        end

        if gen["qg_fixed"]
            qg[i] = gen["qg"]
            @constraint(pm.model, var(pm, :qg, i) == gen["qg"])
        else
            qg[i] = var(pm, :qg, i)
        end
    end
    #Memento.info(LOGGER, "gen expr time: $(time() - start_time)")

    start_time = time()
    for (i,bus) in ref(pm, :bus)
        #PowerModels.constraint_kcl_shunt(pm, i)

        bus_arcs = ref(pm, :bus_arcs, i)
        bus_arcs_dc = ref(pm, :bus_arcs_dc, i)
        bus_gens = ref(pm, :bus_gens, i)
        bus_loads = ref(pm, :bus_loads, i)
        bus_shunts = ref(pm, :bus_shunts, i)

        bus_pd = Dict(k => ref(pm, :load, k, "pd") for k in bus_loads)
        bus_qd = Dict(k => ref(pm, :load, k, "qd") for k in bus_loads)

        bus_gs = Dict(k => ref(pm, :shunt, k, "gs") for k in bus_shunts)
        bus_bs = Dict(k => ref(pm, :shunt, k, "bs") for k in bus_shunts)

        #p = var(pm, :p)
        #q = var(pm, :q)
        #pg = var(pm, :pg)
        #qg = var(pm, :qg)

        @constraint(pm.model, sum(p[a] for a in bus_arcs) == sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs))*(vr[i]^2 + vi[i]^2))
        @constraint(pm.model, sum(q[a] for a in bus_arcs) == sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs))*(vr[i]^2 + vi[i]^2))
    end
    #Memento.info(LOGGER, "power balance constraint time: $(time() - start_time)")

end

""

@everywhere function solution_second_stage!(pm::GenericPowerModel, sol::Dict{String,Any})
    #start_time = time()
    PowerModels.add_setpoint_bus_voltage!(sol, pm)
    #Memento.info(LOGGER, "voltage solution time: $(time() - start_time)")

    #start_time = time()
    PowerModels.add_setpoint_generator_power!(sol, pm)
    #Memento.info(LOGGER, "generator solution time: $(time() - start_time)")

    #start_time = time()
    PowerModels.add_setpoint_branch_flow!(sol, pm)
    #Memento.info(LOGGER, "branch solution time: $(time() - start_time)")

    sol["delta"] = JuMP.value(var(pm, :delta))

    #start_time = time()
    PowerModels.add_setpoint_fixed!(sol, pm, "shunt", "bs", default_value = (item) -> item["bs"])
    #Memento.info(LOGGER, "shunt solution time: $(time() - start_time)")

    #PowerModels.print_summary(sol)
end

"fixes solution degeneracy issues when qg is a free variable, as is the case in PowerFlow"

@everywhere function correct_qg!(network, solution; bus_gens=gens_by_bus(network))
    for (i,gens) in bus_gens
        if length(gens) > 1
            gen_ids = [gen["index"] for gen in gens]
            qgs = [solution["gen"]["$(j)"]["qg"] for j in gen_ids]
            if !isapprox(abs(sum(qgs)), sum(abs.(qgs)))
                #info(LOGGER, "$(i) - $(gen_ids) - $(qgs) - output requires correction!")
                qg_total = sum(qgs)

                qg_remaining = sum(qgs)
                qg_assignment = Dict(j => 0.0 for j in gen_ids)
                for (i,gen) in enumerate(gens)
                    gen_qg = qg_remaining
                    gen_qg = max(gen_qg, gen["qmin"])
                    gen_qg = min(gen_qg, gen["qmax"])
                    qg_assignment[gen["index"]] = gen_qg
                    qg_remaining = qg_remaining - gen_qg
                    if i == length(gens) && abs(qg_remaining) > 0.0
                        qg_assignment[gen["index"]] = gen_qg + qg_remaining
                    end
                end
                #info(LOGGER, "$(qg_assignment)")
                for (j,qg) in qg_assignment
                    solution["gen"]["$(j)"]["qg"] = qg
                end

                sol_qg_total = sum(solution["gen"]["$(j)"]["qg"] for j in gen_ids)
                #info(LOGGER, "$(qg_total) - $(sol_qg_total)")
                @assert isapprox(qg_total, sol_qg_total)

                #info(LOGGER, "updated to $([solution["gen"]["$(j)"]["qg"] for j in gen_ids])")
            end
        end
    end
end

"""
assumes there is one reference bus and one connected component and adjusts voltage
angles to be centered around zero at the reference bus.
"""

@everywhere function correct_voltage_angles!(network)
    ref_bus = -1
    for (i,bus) in network["bus"]
        if bus["bus_type"] == 3
            @assert ref_bus == -1
            ref_bus = bus
        end
    end

    if !isapprox(ref_bus["va"], 0.0, atol=1e-8)
        warn(LOGGER, "shifting voltage angles by $(-ref_bus["va"]) to set reference bus to 0.0")
        shift_voltage_anlges!(network, -ref_bus["va"])
    end
end

"shift networks voltage angles by a specified amount"

@everywhere  function shift_voltage_anlges!(network, shift::Number)
    for (i,bus) in network["bus"]
        bus["va"] = bus["va"] + shift
    end
end

"build a static ordering of all contigencies"

@everywhere function contingency_order(pm_network)
    gen_cont_order = sort(pm_network["gen_contingencies"], by=(x) -> x.label)
    branch_cont_order = sort(pm_network["branch_contingencies"], by=(x) -> x.label)

    gen_cont_total = length(gen_cont_order)
    branch_cont_total = length(branch_cont_order)

    gen_rate = 1.0
    branch_rate = 1.0
    steps = 1

    if gen_cont_total == 0 && branch_cont_total == 0
        # defaults are good
    elseif gen_cont_total == 0 && branch_cont_total != 0
        steps = branch_cont_total
    elseif gen_cont_total != 0 < branch_cont_total == 0
        steps = gen_cont_total
    elseif gen_cont_total == branch_cont_total
        steps = branch_cont_total
    elseif gen_cont_total < branch_cont_total
        gen_rate = 1.0
        branch_rate = branch_cont_total/gen_cont_total
        steps = gen_cont_total
    elseif gen_cont_total > branch_cont_total
        gen_rate = gen_cont_total/branch_cont_total
        branch_rate = 1.0
        steps = branch_cont_total
    end

    cont_order = []
    gen_cont_start = 1
    branch_cont_start = 1
    for s in 1:steps
        gen_cont_end = min(gen_cont_total, trunc(Int,ceil(s*gen_rate)))
        #println(gen_cont_start:gen_cont_end)
        for j in gen_cont_start:gen_cont_end
            push!(cont_order, gen_cont_order[j])
        end
        gen_cont_start = gen_cont_end+1

        branch_cont_end = min(branch_cont_total, trunc(Int,ceil(s*branch_rate)))
        #println("$(s) - $(branch_cont_start:branch_cont_end)")
        for j in branch_cont_start:branch_cont_end
            push!(cont_order, branch_cont_order[j])
        end
        branch_cont_start = branch_cont_end+1
    end

    @assert(length(cont_order) == gen_cont_total + branch_cont_total)

    return cont_order
end

"checks feasibility criteria of contingencies, corrects when possible"

@everywhere function correct_contingency_solution!(network, cont_sol; bus_gens = gens_by_bus(network)) # Function for post-processing
    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"] > 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"] < 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]
    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

    gen_id = -1
    if cont_sol["cont_type"] == "gen"
        gen_id = cont_sol["cont_comp_id"]
    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"] + nw_gen["alpha"]*delta
            pg_calc = max(pg_calc, nw_gen["pmin"])
            pg_calc = min(pg_calc, nw_gen["pmax"])

            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
    PowerModels.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

@everywhere 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",
    ]
    println(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),
    ]
    println(join(data, ", "))
end

@everywhere function write_solution2_contingency(io::IO, pm_network, contingency_solution) # Writes solution2 for contingency
    base_mva = pm_network["baseMVA"]

    bus_switched_shunt_b = Dict(i => 0.0 for (i,bus) in pm_network["bus"])
    for (i,nw_shunt) in pm_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

    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,bus) in contingency_solution["bus"]
        nw_bus = pm_network["bus"][i]
        write(io, "$(nw_bus["index"]), $(bus["vm"]), $(rad2deg(bus["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,gen) in contingency_solution["gen"]
        nw_gen = pm_network["gen"][i]
        bus_index = nw_gen["source_id"][2]
        gen_id = nw_gen["source_id"][3]
        write(io, "$(bus_index), $(gen_id), $(base_mva*gen["pg"]), $(base_mva*gen["qg"])\n")
    end

    write(io, "-- delta section\n")
    write(io, "delta(MW)\n")
    write(io, "$(base_mva*contingency_solution["delta"])\n")
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

function compute_power_balance_deltas!(network)
    flows = PowerModels.calc_branch_flow_ac(network)
    PowerModels.update_data!(network, flows)

    balance = PowerModels.calc_power_balance(network)
    PowerModels.update_data!(network, balance)

    p_delta_abs = [abs(bus["p_delta"]) for (i,bus) in network["bus"] if bus["bus_type"] != 4]
    q_delta_abs = [abs(bus["q_delta"]) for (i,bus) in network["bus"] if bus["bus_type"] != 4]

    return (
        p_delta_abs_max = maximum(p_delta_abs),
        p_delta_abs_mean = mean(p_delta_abs),
        q_delta_abs_max = maximum(q_delta_abs),
        q_delta_abs_mean = mean(q_delta_abs),
    )
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


text:
using DataFrames, Documenter, CrossfilterCharts

makedocs(
    # options
    modules = [CrossfilterCharts]
)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/tawheeler/CrossfilterCharts.jl.git",
    julia  = "release",
    osname = "linux"
)

ext: jl


text:
# struct to model a subset

@tm1def mutable struct Subset
  Name::Union{String,Nothing}
  UniqueName::Union{String,Nothing}
  Expression::Union{String,Nothing}
  Attributes::Union{Dict,Nothing}
  Alias::Union{String,Nothing}
  Hierarchy::Union{Hierarchy,Nothing}
  Elements::Union{Vector{Element},Nothing}
end

Subset(name::AbstractString) = Subset(Dict("Name" => name))

namefield(subset::Subset) = subset.Name

# functions for endpoints

@api_default function get_subset(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString,
  subset_name::AbstractString,
  private::Bool = false;
  options...,
)
  subset_type = subset_type_string(private)

  # exact parameters to pass need review
  params = Dict("\$expand" => "Hierarchy(\$select=Dimension,Name),Elements(\$select=Name)")

  result = tm1_get_json(
    api,
    "Dimensions('" *
    dimension_name *
    "')/Hierarchies('" *
    hierarchy_name *
    "')/" *
    subset_type *
    "('" *
    subset_name *
    "')";
    params = params,
    options...,
  )
  Subset(result)
end

@api_default function delete_subset(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString,
  subset_name::AbstractString,
  private::Bool = false;
  options...,
)
  subset_type = subset_type_string(private)

  tm1_delete(
    api,
    "Dimensions('" *
    dimension_name *
    "')/Hierarchies('" *
    hierarchy_name *
    "')/" *
    subset_type *
    "('" *
    subset_name *
    "')";
    options...,
  )
end

function subset_type_string(private::Bool)
  if private
    subset_type = "PrivateSubsets"
  else
    subset_type = "Subsets"
  end
  return subset_type
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot, ColorSchemes
# using Random, Colors # hide
#Random.seed!(123) # hide
#function custom_palette(colormap=:viridis) # hide
#    cmap = get(colorschemes[colormap], LinRange(0,1,256)) # hide
#    ctmp = "0 '#$(hex(cmap[1]))'," # hide
#    for i in 2:256; ctmp = ctmp*"$(i-1) '#$(hex(cmap[i]))'," end; # hide
#    "set palette defined("*ctmp[1:end-1]*")" # hide
#end # hide
# Archimedes spiral
a = 1.5
b = -2.4
t = LinRange(0,5*π,500)
x = (a .+ b*t) .* cos.(t)
y = (a .+ b*t) .* sin.(t)
cmaps = [:cool, :viridis, :plasma, :inferno, :thermal,
    :leonardo, :winter, :spring, :ice]
@gp "set multiplot layout 3,3; set key off;
    unset ytics; unset xtics; unset border" :-
for i in 1:9
    @gp :- i title = "$(cmaps[i])" "set size square" :-
    @gp :- x y t "w l lw 3 dt 1 lc palette" palette(cmaps[i]) :-
    @gp :- "set cbtics out nomirror" :-
end
@gp
save(term="pngcairo size 900,800", output="./code/plt_ex9.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
"""
Bouncing ball
q = (x,z), u = (xdot, zdot),
g(q) = q[2]
W = del g / del q = [0; 1]
M = [m 0; 0 m]
h = [0 -m*g]
"""

using LinearAlgebra
using PyPlot
using TimeStepping

const mass = 1.0
const g = 9.81
qA = [0.0; 1.0]
uA = [1.0; 2.0]

function gap(q::AbstractArray, u::AbstractArray)
    W = zeros(Float32, 2,1)
    W[1,1] = 0.0
    W[2,1] = 1.0
    return [q[2]], W
end

function dynamics(q::AbstractArray, u::AbstractArray)
    return diagm(ones(Float32, 2)*mass), [0; -mass*g]
end

bouncing_ball = Integrator(gap, dynamics, convert.(Float32, qA),
    convert.(Float32, uA); Δt=1e-3, ε=0.5)
integrate(bouncing_ball, 2.0)

fig = figure(1, clear=true)
fig.add_subplot(1,1,1)
fig.axes[1].clear()
fig.axes[1].plot(bouncing_ball.t, getindex.(bouncing_ball.q,2), "k", linewidth=2)

fig.axes[1].set_title("Bouncing Ball (z-axis)", fontsize=18)
fig.axes[1].set_xlabel(L"$t$", fontsize=18)
fig.axes[1].set_ylabel(L"$z(t)$", fontsize=18)

fig.axes[1].autoscale_view()
fig.tight_layout()


text:
# ──────────────────────────────────────────────────────────────
#
# Diag
# Wrapper around Diagonal, extending its functionality.
# NOTE: we do not extend Diagonal to avoid contaminating the base
# NOTE: Julia 0.5 has all the functionality needed here. Can
# 		replace with the generic Diagonal once it comes out.
#       I'm looking forward to getting rid of this forever.
#
# ──────────────────────────────────────────────────────────────

import Base:+,*,-,\,^,sparse

ViewTypes   = Union{SubArray}
VectorTypes = Union{AbstractMatrix, Vector, ViewTypes}
MatrixTypes = Union{AbstractMatrix, Array{Real,2},
                    SparseMatrixCSC{Real,Integer}}

type Diag
    diag::Vector
end

+(A::Diag,B::Diag)             = Diag(A.diag + B.diag)
+(A::Diag,T::UniformScaling)   = Diag(A.diag + T.λ)
+(T::UniformScaling,A::Diag)   = Diag(A.diag + T.λ)
*(T::UniformScaling,A::Diag)   = Diag(A.diag*T.λ)
*(A::Diag,T::UniformScaling)   = Diag(A.diag*T.λ)
*(α::Real,A::Diag)             = Diag(A.diag*α)
*(A::Diag,α::Real)             = Diag(A.diag*α)
*(A::Diag,B::Diag)             = Diag(A.diag.*B.diag)
function +(A::MatrixTypes,B::Diag)
    O = copy(A)
    for i = 1:size(A,2); O[i,i] = A[i,i] + B.diag[i]; end
    return O
end
*(A::Diag,B::AbstractMatrix)   = A.diag.*B
*(A::Diag,b::AbstractVector)   = A.diag.*b
+(B::Diag,A::MatrixTypes)      = A + B
^(A::Diag, n::Integer)         = Diag(A.diag.^n)
Base.full(A::Diag)             = full(Diagonal(A.diag))
Base.inv(A::Diag)              = Diag(1.0./A.diag)
*(A::Diag, b::ViewTypes)       = A.diag.*b
Base.size(A::Diag, k::Integer) = (k == 1 || k == 2) ? size(A.diag,1) : 1
Base.size(A::Diag)             = (size(A.diag,1), size(A.diag,1))
Id(n::Integer)                 = Diag(ones(n))
Base.sparse(A::Diag)           = spdiagm(A.diag)

Base.Ac_mul_B(O1::Diag, O2::Diag) = O1*O2

# returns a pointer to the original matrix, this is consistent with the
# behavior of Symmetric in Base.
Base.ctranspose(A::Diag) = A
Base.Ac_mul_B(A::Diag, x::VectorTypes) = A*x


text:
module Monotones

export relative_entropy
import RTFOOL.Context

function relative_entropy(ps::AbstractVector{Float64}, qs::AbstractVector{Float64})
    if length(ps) != length(qs)
        error("states have different sizes")
    elseif !(sum(ps) ≈ 1.0) || !(sum(qs) ≈ 1.0)
        error("state does not sum to 1.0")
    end
    sum(map((p,q) -> (p > 0.0) ? p*log(p/q) : 0.0, ps, qs))
end
relative_entropy(ctx::Context) = relative_entropy(ctx.system_state, ctx.bath_state)

end


text:
using GLVisualize, GeometryTypes, Reactive, GLAbstraction

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end
# let the visualization rotate later on
rotation = map(rotationmatrix_z, const_lift(*, timesignal, 2f0*pi))
# create some random 3D vectors
vectors3d = rand(Vec3f0, 5,5,5)
# this is not the best way to use 2D sprites, but this will space them on
# a 3D grid and use the rotation from `vectors3d` and the length of them
# to look up the a color from the optional keyword argument `color_map`.
arrows = visualize(('➤', vectors3d), scale=Vec2f0(0.1), model=rotation)

_view(arrows, camera=:perspective)

if !isdefined(:runtests)
    renderloop(window)
end


text:
using Documenter, DocTestPkg

makedocs(
  pages = Any[
  "home" => "index.md",
  "main" => "main.md",
  "Other" => "other.md",
  "Other2" => "other2.md"
  ]
)

deploydocs(
  deps = Deps.pip("mkdocs", "python-markdown-math"),
  repo = "github.com/JaredCrean2/DocTestPkg.jl.git",
  julia = "0.4",
  osname = "linux",
)


text:
# (C) Jarett Revels ForwardDiff.jl, License MIT
function tupexpr(f,N)
    ex = Expr(:tuple, [f(i) for i=1:N]...)
    return quote
        $(Expr(:meta, :inline))
        @inbounds return $ex
    end
end

@inline zero_tuple(::Type{Tuple{}}) = tuple()

@generated function zero_tuple{N,T}(::Type{NTuple{N,T}})
    result = tupexpr(i -> :z, N)
    return quote
        $(Expr(:meta, :inline))
        z = zero($T)
        return $result
    end
end

@generated function one_tuple{N,T}(::Type{NTuple{N,T}})
    result = tupexpr(i -> :z, N)
    return quote
        $(Expr(:meta, :inline))
        z = one($T)
        return $result
    end
end

@inline rand_tuple(::Type{Tuple{}}) = tuple()

@generated function rand_tuple{N,T}(::Type{NTuple{N,T}})
    return tupexpr(i -> :(rand($T)), N)
end

for f in UNARY_FUNCS
    tuple_f_string = symbol(string(f) * "_tuple")
    @eval begin
        @generated function $(tuple_f_string){N}(a::NTuple{N})
            return quote
                $(Expr(:meta, :inline))
                $(tupexpr(i -> :($($f)(a[$i])), N))
            end
        end
    end
end

for f in BINARY_FUNCS
    tuple_f_string = symbol(string(f) * "_tuple")
    @eval begin
        @generated function $(tuple_f_string){N}(a::NTuple{N}, b::NTuple{N})
            return quote
            $(Expr(:meta, :inline))
            $(tupexpr(i -> :($($f)(a[$i], b[$i])), N))
            end
        end

        @generated function $(tuple_f_string){N}(a::NTuple{N}, b::Number)
            return quote
                $(Expr(:meta, :inline))
                $(tupexpr(i -> :($($f)(a[$i], b)), N))
            end
        end

        @generated function $(tuple_f_string){N}(b::Number, a::NTuple{N})
            return quote
                $(Expr(:meta, :inline))
                $(tupexpr(i -> :($($f)(b, a[$i])), N))
            end
        end
    end
end


text:
using GeometryTypes

abstract type TrialState end

mutable struct FixationState <: TrialState
    area::Rectangle
    clock::ClockTimer
    current_pos::Point2f0
end

mutable struct TargetState <: TrialState
    area::Rectangle
    clock::ClockTimer
    current_pos::Point2f0
end

mutable struct ResponseState <: TrialState
    area::Rectangle
    clock::ClockTimer
    current_pos::Point2f0
end

mutable struct TrialEndState <: TrialState
    clock::ClockTimer
end

mutable struct StateTransitions
    states::Vector{TrialState}
    edges::Matrix{Int64}
    current::Int64
end

mutable struct ExperimentRecord
    clock::Clock
    timestamp::Vector{Float64}
    event::Vector{String}
end

ExperimentRecord() = ExperimentRecord(Clock(), Float64[], String[])

function next!(tt::StateTransitions,engaged::Bool, record::ExperimentRecord)
    nn = tt.current
    if engaged
        nn = tt.edges[1,tt.current]
    else
        nn = tt.edges[2, tt.current]
    end
    #log this event
    t1 = get_time(record.clock)
    push!(record.timestamp, t1)
    push!(record.event, "$(tt.current) -> $(nn)")
    #make sure we reset the clock
    reset!(tt.states[nn].clock)
    tt.current = nn
    tt.current
end

check_state(state::TrialState) = error("Not implemented")
next(state::TrialState, engaged::Bool) = error("Not implemented")
get_position(state::TrialState) = state.current_pos
get_area(state::TrialState) = state.area
get_area(state::TrialEndState) = Rectangle(-Inf32, -Inf32, Inf32, Inf32)

function set_position!(state::TrialState, pos::Point2f0)
    state.current_pos = pos
    nothing
end
set_position!(state::TrialEndState,pos::Point2f0) = nothing

function next!(transitions::StateTransitions, record::ExperimentRecord)
    state = transitions.states[transitions.current]
    b = check_state(state)
    #check if we are transitioning and if so to which state
    if b == 1
        next!(transitions, true, record)
    elseif b == 2
        next!(transitions, false, record)
    end
    nothing
end

"""
Check whether we should transition out of `state`. Returns 0 if we should stay in the same state, 1 if we should transition to the next valid state and 2 if we should transition to the next invalid state.
"""
function check_state(state::Union{FixationState,ResponseState})
    b = 0
    if is_inside(state)
        if is_done!(state.clock)
            b = 1
        end
    else
        b = 2
    end
    b
end

function check_state(state::TrialEndState)
    b = 0
    if is_done!(state.clock)
        b = 1
    end
    b
end

function is_inside(state::TrialState)
    pos = get_position(state)
    is_inside(state, pos)
end

function is_inside(state::TrialState, pos::Point2f0)
    area = get_area(state)
    b = (area.x <= pos[1] <= area.x + area.w)
    b &= (area.y <= pos[2] <= area.y + area.h)
    return b
end


text:
#!/usr/bin/env julia
#=
Common struct
del2z <delta.z@aliyun.com>
=#
module Data

include("loader.jl")

end


text:
########### Defines the parameters & plot for Fig 2 a,c,e,g,i ###########
c = 0.5; # set the host contact rate as a constant 0.5
# import required variables/modules common to all simulations
include("../Model_scripts/fig2params.jl");
include("../Model_scripts/Functions.jl");
# now that we have all the requirements, let's run the simulation
run_simulation();
Np = N';
########### Figure 2a ###########
# plot the number of infected/uninfected hosts
labels = vcat("Infected", ["" for i in 1:1:n_parasites]);
plot_population_numbers(
    hcat(labels...),
    "c=0.5\n \nNumber of infected and uninfected hosts",
    100,
    size(Np)[1],
    "$figure_directory/graph_2a.png"
)
########### Figure 2e ###########
# calculate the weighted average mortality
ui_w_avg = calculate_mortality()
# calculate weighted V0
V0_w = calculate_V0()
# plot the average V0
plot_vertical_reproductive_ratio(1, size(Np)[1], "$figure_directory/graph_2e.png")
########### Figure 2c ###########
# calculate the weighted average β
βy_w_avg = calculate_horizontal_transmission()
# calculate weighted H0
H0_w = calculate_H0()
# plot average weighted R0
plot_reproductive_rate(5, size(Np)[1], "$figure_directory/graph_2c.png")
########### Figure 2g ###########
# calculate the weighted virulence (for noise reduction)
vir_w_avg = calculate_average_virulence()
# plot the average horizontal transmission & weighted virulence
plot_horizontal_and_virulence(1, size(Np)[1], "$figure_directory/graph_2g.png")
########### Figure 2i ###########
# calculate & plot the population evenness through time
plot_evenness(1, size(Np)[1], "$figure_directory/graph_2i.png")


text:
@testset "547.friend-circles.jl" begin
    @test find_circle_num([[1, 1, 0], [1, 1, 1], [0, 1, 1]]) == 1
    @test find_circle_num([[1, 1, 0], [1, 1, 0], [0, 0, 1]]) == 2
    @test find_circle_num([[1, 0, 0], [0, 1, 0], [0, 0, 1]]) == 3
    @test find_circle_num([[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
                            0, 0, 1, 0, 0, 0, 0, 0, 0],
                           [1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 1, 0, 0],
                           [0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 1, 0, 0],
                           [1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
                            1, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
                            0, 0, 0, 1, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                            0, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                            1, 0, 0, 0, 0, 0, 0, 0, 0],
                           [0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 1, 0, 0, 0, 0, 0, 0, 0],
                           [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 1, 0, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
                            1, 0, 0, 1, 0, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 1, 0, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            1, 0, 0, 0, 0, 1, 0, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 1, 0, 0],
                           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                            0, 0, 0, 0, 0, 0, 0, 1, 0],
                           [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
                            0, 0, 0, 0, 0, 0, 0, 0, 1]]) == 11
end


text:
function initialize!(p::GridVisualizer,::Type{PlutoVistaType})
    PlutoVista=p.context[:Plotter]
    layout=p.context[:layout]
    @assert(layout==(1,1))
    p.context[:scene]=PlutoVista.PlutoVistaPlot(resolution=p.context[:resolution])
    PlutoVista.backend!(p.context[:scene])
    for I in CartesianIndices(layout)
        ctx=p.subplots[I]
        ctx[:figure]=p.context[:scene]
    end
end

function reveal(p::GridVisualizer,::Type{PlutoVistaType})
    p.context[:scene]
end

function reveal(ctx::SubVisualizer,TP::Type{PlutoVistaType})
    if ctx[:show]||ctx[:reveal]
        reveal(ctx[:GridVisualizer],TP)
    end
end

function gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{1}}, grid)
    PlutoVista=ctx[:Plotter]
    coord=grid[Coordinates]

    cellregions=grid[CellRegions]
    cellnodes=grid[CellNodes]
    coord=grid[Coordinates]
    ncellregions=grid[NumCellRegions]
    bfacenodes=grid[BFaceNodes]
    bfaceregions=grid[BFaceRegions]
    nbfaceregions=grid[NumBFaceRegions]
    ncellregions=grid[NumCellRegions]

    crflag=ones(Bool,ncellregions)
    brflag=ones(Bool,nbfaceregions)

    xmin=minimum(coord)
    xmax=maximum(coord)
    h=(xmax-xmin)/20.0

    #    ax.set_aspect(ctx[:aspect])
#    ax.get_yaxis().set_ticks([])
#    ax.set_ylim(-5*h,xmax-xmin)
    cmap=region_cmap(ncellregions)

    for icell=1:num_cells(grid)
        ireg=cellregions[icell]
        label = crflag[ireg] ? "c$(ireg)" : ""
        crflag[ireg]=false

        x1=coord[1,cellnodes[1,icell]]
        x2=coord[1,cellnodes[2,icell]]

        PlutoVista.plot!(ctx[:figure],[x1,x2],[0,0],clear=false,linewidth=3.0,color=rgbtuple(cmap[cellregions[icell]]),label=label)
        PlutoVista.plot!(ctx[:figure],[x1,x1],[-h,h],clear=false,linewidth=ctx[:linewidth],color=:black)
        PlutoVista.plot!(ctx[:figure],[x2,x2],[-h,h],clear=false,linewidth=ctx[:linewidth],color=:black)
    end

    cmap=bregion_cmap(nbfaceregions)
    for ibface=1:num_bfaces(grid)
        ireg=bfaceregions[ibface]
        if ireg >0
            label = brflag[ireg] ? "b$(ireg)" : ""
            brflag[ireg]=false
            x1=coord[1,bfacenodes[1,ibface]]
            PlutoVista.plot!(ctx[:figure],[x1,x1],[-2*h,2*h],clear=false,linewidth=3.0,color=rgbtuple(cmap[ireg]),label=label,legend=leglocs[ctx[:legend]])
        end
    end
    reveal(ctx,TP)
end

function scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{1}}, grid,func)
    PlutoVista=ctx[:Plotter]
    coord=grid[Coordinates]
    PlutoVista.plot!(ctx[:figure],coord[1,:],func,
                     color=ctx[:color],
                     markertype=ctx[:markershape],
                     markercount=length(func)÷ctx[:markevery],
                     linestyle=ctx[:linestyle],
                     xlabel=ctx[:xlabel],
                     ylabel=ctx[:ylabel],
                     label=ctx[:label],
                     linewidth=ctx[:linewidth],
                     legend=ctx[:legend],
                     xlimits=ctx[:xlimits],
                     ylimits=ctx[:flimits],
                     clear=ctx[:clear]
                     )
    reveal(ctx,TP)
end

gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{2}}, grid)=nothing
function scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{2}}, grid,func)
    PlutoVista=ctx[:Plotter]
    pts=grid[Coordinates]
    tris=grid[CellNodes]

    isolines=ctx[:isolines]
    fmin,fmax=ctx[:flimits]

    if fmin<fmax
        isolines=collect(range(fmin,fmax,length=isolines))
    end
    PlutoVista.backend!(ctx[:figure],backend=ctx[:backend],datadim=2)
    PlutoVista.tricontour!(ctx[:figure],pts,tris,func,
                           colormap=ctx[:colormap],
                           isolines=isolines
                           )
    reveal(ctx,TP)

end

gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{3}}, grid)=nothing
scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{3}}, grid,func)=nothing


text:
const DT = Integer

function find_amicable(limit::DT)
  T = typeof(limit)
  v = Vector{Tuple{Pair{T, Vector{T}}, Pair{T, Vector{T}}}}()
  t = Set{T}()

  for n ∈ one(T):limit
    d₁ = divisors(n)
    p = sum(d₁)
    d₂ = divisors(p)
    s = sum(d₂)

    if s == n && d₁ ≠ d₂ && s ∉ t
      push!(t, s)
      push!(t, p)
      push!(v, (n => d₁, p => d₂))
    end
  end

  (t, v)
end

# (Set([2924, 5564, 6232, 5020, 220, 284, 1184, 2620, 1210, 6368]),
#  [(220 => [1, 2, 4, 5, 10, 11, 20, 22, 44, 55, 110], 284 => [1, 2, 4, 71, 142]), (1184 => [1, 2, 4, 8, 16, 32, 37, 74, 148, 296, 592], 1210 => [1, 2, 5, 10, 11, 22, 55, 110, 121, 242, 605]), (2620 => [1, 2, 4, 5, 10, 20, 131, 262, 524, 655, 1310], 2924 => [1, 2, 4, 17, 34, 43, 68, 86, 172, 731, 1462]), (5020 => [1, 2, 4, 5, 10, 20, 251, 502, 1004, 1255, 2510], 5564 => [1, 2, 4, 13, 26, 52, 107, 214, 428, 1391, 2782]), (6232 => [1, 2, 4, 8, 19, 38, 41, 76, 82, 152, 164, 328, 779, 1558, 3116], 6368 => [1, 2, 4, 8, 16, 32, 199, 398, 796, 1592, 3184])])
#
# julia> s |> sum
# 31626

function is_amicable(n₁::Any, n₂::Any)
  throw(ArgumentError("n₁ and n₂ should be Integers"))
end

function is_amicable(n₁::DT, n₂::DT)::Bool
  d₁ = divisors(n₁)
  d₂ = divisors(n₂)
  n₂ == sum(d₁) && n₁ == sum(d₂)
end

function divisors(n::DT)::Vector{<: DT}
  T = typeof(n)
  ary::Vector{T} = [one(T), ]

  d = T(2)
  while d * d ≤ n
    q, r = divrem(n, d)
    if r == 0
      q == d ? push!(ary, d) : push!(ary, d, q)
    end
    d += 1
  end

  return ary |> sort
end

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
Base.show(io::IO, db::Database) = print(io, "Database($(db.client), \"", db.name, "\")")
Base.getindex(database::Database, collection_name::String) = Collection(database, collection_name)

"""
    command_simple(database::Database, command::BSON) :: BSON
    command_simple(collection::Collection, command::BSON) :: BSON

Executes a `command` given by a JSON string or a BSON instance.

It returns the first document from the result cursor.

# Example

```julia
julia> client = Mongoc.Client() # connects to localhost at port 27017
Client(URI("mongodb://localhost:27017"))

julia> bson_result = Mongoc.command_simple(client[\"admin\"], "{ \"ping\" : 1 }")
BSON("{ "ok" : 1.0 }")
```

# C API

* [`mongoc_database_command_simple`](http://mongoc.org/libmongoc/current/mongoc_database_command_simple.html)
"""
function command_simple(database::Database, command::BSON) :: BSON
    reply = BSON()
    err_ref = Ref{BSONError}()
    ok = mongoc_database_command_simple(database.handle, command.handle, C_NULL, reply.handle,
        err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

"""
    write_command(database::Database, command::BSON;
        options::Union{Nothing, BSON}=nothing) :: BSON

Issue a command with write semantics.

See http://mongoc.org/libmongoc/current/mongoc_database_read_write_command_with_opts.html.
"""
function write_command(database::Database, command::BSON;
        options::Union{Nothing, BSON}=nothing) :: BSON

    options_handle = options == nothing ? C_NULL : options.handle
    reply = BSON()
    err_ref = Ref{BSONError}()
    ok = mongoc_database_write_command_with_opts(database.handle, command.handle,
            options_handle, reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

"""
    read_command(database::Database, command::BSON;
        options::Union{Nothing, BSON}=nothing) :: BSON

Issue a command with read semantics.

See http://mongoc.org/libmongoc/current/mongoc_database_read_command_with_opts.html.
"""
function read_command(database::Database, command::BSON;
        options::Union{Nothing, BSON}=nothing) :: BSON

    options_handle = options == nothing ? C_NULL : options.handle
    reply = BSON()
    err_ref = Ref{BSONError}()
    ok = mongoc_database_read_command_with_opts(database.handle, command.handle,
            C_NULL, options_handle, reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

"""
    add_user(database::Database, username::String, password::String, roles::Union{Nothing, BSON},
        custom_data::Union{Nothing, BSON}=nothing)

This function shall create a new user with access to database.

**Warning:** Do not call this function without TLS.
"""
function add_user(database::Database, username::String, password::String,
        roles::Union{Nothing, BSON}, custom_data::Union{Nothing, BSON}=nothing)

    err_ref = Ref{BSONError}()
    roles_handle = roles == nothing ? C_NULL : roles.handle
    custom_data_handle = custom_data == nothing ? C_NULL : custom_data.handle
    ok = mongoc_database_add_user(database.handle, username, password, roles_handle,
                                  custom_data_handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end

"""
    remove_user(database::Database, username::String)

Removes a user from database.
"""
function remove_user(database::Database, username::String)
    err_ref = Ref{BSONError}()
    ok = mongoc_database_remove_user(database.handle, username, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end

"""
    has_user(database::Database, user_name::String) :: Bool

Checks if `database` has a user named `user_name`.
"""
function has_user(database::Database, user_name::String) :: Bool
    cmd_result = command_simple(database, BSON("""{ "usersInfo": "$user_name" }"""))
    return !isempty(cmd_result["users"])
end

"""
    find_collections(database::Database; options::Union{Nothing, BSON}=nothing) :: Cursor

Queries for collections in a `database`.
"""
function find_collections(database::Database; options::Union{Nothing, BSON}=nothing) :: Cursor
    options_handle = options == nothing ? C_NULL : options.handle
    cursor_handle = mongoc_database_find_collections_with_opts(database.handle, options_handle)
    if cursor_handle == C_NULL
        error("Couldn't execute query.")
    end
    return Cursor(database, cursor_handle)
end

"""
    get_collection_names(database::Database;
        options::Union{Nothing, BSON}=nothing) :: Vector{String}

Helper method to get collection names.

See also [`Mongoc.find_collections`](@ref).
"""
function get_collection_names(database::Database;
                              options::Union{Nothing, BSON}=nothing) :: Vector{String}
    result = Vector{String}()
    for bson_collection in find_collections(database, options=options)
        push!(result, bson_collection["name"])
    end
    return result
end

# docstring for drop is at collection.jl
function drop(database::Database;
        options::Union{Nothing, BSON}=nothing)

    err_ref = Ref{BSONError}()
    opts_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_database_drop_with_opts(database.handle, opts_handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    nothing
end


text:
include("task_func.jl")
function drawAll(r)
    listOfShags = []
    while !isborder(r, West) || !isborder(r, Sud)
        push!(listOfShags, goToSide!(r, West, save = true))
        push!(listOfShags, goToSide!(r, Sud, save = true))
    end
    side = Ost
    putmarker!(r)
    goToSide!(r, Ost, mark = true)
    side = reverseSide(side)
    while !isborder(r, Nord)
        move!(r, Nord)
        putmarker!(r)
        while true
            goToSide!(r, side, mark = true)
            if isRect(r, side)
                aroundBorder!(r, side)
                putmarker!(r)
            else
                break
            end
        end
        side = reverseSide(side)
    end
    goToSide!(r, West)
    goToSide!(r, Sud)
    goToStart!(r, listOfShags)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
module TestNormalizer

using Dates
using Test
using Random
using Statistics
using DataFrames
using TSML

function generatedf()
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval1 = rand(length(gdate))
    gval2 = rand(length(gdate))
    gval3 = rand(length(gdate))
    X = DataFrame(Date=gdate,Value1=gval1,Value2=gval2,Value3=gval3)
    X
end

function test_normalizer()
   Random.seed!(123)
   X = generatedf()
   norm = Normalizer(Dict(:method => :zscore))
   res=fit_transform!(norm,X)
   @test isapprox(mean(res[:,1]),0.0,atol=1e-8)
   @test isapprox(mean(res[:,2]),0.0,atol=1e-8)
   @test isapprox(std(res[:,1]),1.0,atol=1e-8)
   @test isapprox(std(res[:,2]),1.0,atol=1e-8)
   res=fit_transform(norm,X)
   @test isapprox(mean(res[:,1]),0.0,atol=1e-8)
   @test isapprox(mean(res[:,2]),0.0,atol=1e-8)
   @test isapprox(std(res[:,1]),1.0,atol=1e-8)
   @test isapprox(std(res[:,2]),1.0,atol=1e-8)
   norm = Normalizer(Dict(:method => :unitrange))
   res=fit_transform!(norm,X)
   @test isapprox(minimum(res[:,1]),0.0,atol=1e-8)
   @test isapprox(minimum(res[:,2]),0.0,atol=1e-8)
   @test isapprox(maximum(res[:,1]),1.0,atol=1e-8)
   @test isapprox(maximum(res[:,2]),1.0,atol=1e-8)
   res=fit_transform(norm,X)
   @test isapprox(minimum(res[:,1]),0.0,atol=1e-8)
   @test isapprox(minimum(res[:,2]),0.0,atol=1e-8)
   @test isapprox(maximum(res[:,1]),1.0,atol=1e-8)
   @test isapprox(maximum(res[:,2]),1.0,atol=1e-8)
   norm = Normalizer(Dict(:method => :pca))
   res=fit_transform!(norm,X)
   @test isapprox(std(res[:,1]),0.28996,atol=0.1)
   res=fit_transform(norm,X)
   @test isapprox(std(res[:,1]),0.28996,atol=0.1)
   norm = Normalizer(Dict(:method => :fa))
   res = fit_transform!(norm,X)
   @test isapprox(std(res[:,1]),0.81670,atol=0.1)
   res = fit_transform(norm,X)
   @test isapprox(std(res[:,1]),0.81670,atol=0.1)
   norm = Normalizer(Dict(:method => :ppca))
   res = fit_transform!(norm,X)
   @test isapprox(std(res[:,1]),0.11169,atol=0.1)
   res = fit_transform(norm,X)
   @test isapprox(std(res[:,1]),0.11169,atol=0.1)
   norm = Normalizer(Dict(:method => :log))
   res = fit_transform!(norm,X)
   @test isapprox(std(res[:,1]),0.99941,atol=0.1)
   res = fit_transform(norm,X)
   @test isapprox(std(res[:,1]),0.99941,atol=0.1)
   norm = Normalizer(Dict(:method => :sqrt))
   res = fit_transform!(norm,X)
   @test isapprox(std(res[:,1]),0.2355,atol=0.1)
   res = fit_transform(norm,X)
   @test isapprox(std(res[:,1]),0.2355,atol=0.1)
end
@testset "Normalizer: zscore, unitrange, pca, ppca, fa" begin
   test_normalizer()
end

end

max_stars_repo_name: IBM/tsml

max_issues_repo_name: IBM/tsml

max_forks_repo_name: IBM/tsml


text:
struct NewtonTrustRegion
end

struct NewtonTrustRegionCache{Tx} <: AbstractSolverCache
    x::Tx
    xold::Tx
    r::Tx
    r_predict::Tx
    p::Tx
    p_c::Tx
    pi::Tx
    d::Tx
end
function NewtonTrustRegionCache(df)
    x = similar(df.x_f) # Current point
    xold = similar(x) # Old point
    r = similar(df.F)       # Current residual
    r_predict = similar(x)  # predicted residual
    p = similar(x)          # Step
    p_c = similar(x)          # Cauchy point
    pi = similar(x)
    d = similar(x)          # Scaling vector
    NewtonTrustRegionCache(x, xold, r, r_predict, p, p_c, pi, d)
end
macro trustregiontrace(stepnorm)
    esc(quote
        if tracing
            dt = Dict()
            if extended_trace
                dt["x"] = copy(cache.x)
                dt["f(x)"] = copy(value(df))
                dt["g(x)"] = copy(jacobian(df))
                dt["delta"] = delta
                dt["rho"] = rho
            end
            update!(tr,
                    it,
                    maximum(abs, value(df)),
                    $stepnorm,
                    dt,
                    store_trace,
                    show_trace)
        end
    end)
end

function dogleg!(p::AbstractArray{T}, p_c::AbstractArray{T}, p_i,
                 r::AbstractArray{T}, d::AbstractArray{T}, J::AbstractMatrix{T},
                 delta::Real) where T
    try
        copyto!(p_i, J \ vec(r)) # Gauss-Newton step
    catch e
        if isa(e, LAPACKException) || isa(e, SingularException)
            # If Jacobian is singular, compute a least-squares solution to J*x+r=0
            U, S, V = svd(convert(Matrix{T}, J)) # Convert to full matrix because sparse SVD not implemented as of Julia 0.3
            k = sum(S .> eps())
            mrinv = V * Matrix(Diagonal([1 ./ S[1:k]; zeros(eltype(S), length(S)-k)])) * U' # Moore-Penrose generalized inverse of J
            vecpi = vec(p_i)
            mul!(vecpi,mrinv,vec(r))
        else
            throw(e)
        end
    end
    rmul!(p_i, -one(T))

    # Test if Gauss-Newton step is within the region
    if wnorm(d, p_i) <= delta
        copyto!(p, p_i)   # accepts equation 4.13 from N&W for this step
    else
        # For intermediate we will use the output array `p` as a buffer to hold
        # the gradient. To make it easy to remember which variable that array
        # is representing we make g an alias to p and use g when we want the
        # gradient

        # compute g = J'r ./ (d .^ 2)
        g = p
        mul!(vec(g), transpose(J), vec(r))
        g .= g ./ d.^2

        # compute Cauchy point
        p_c .= -wnorm(d, g)^2 / sum(abs2, J*vec(g)) .* g

        if wnorm(d, p_c) >= delta
            # Cauchy point is out of the region, take the largest step along
            # gradient direction
            rmul!(g, -delta/wnorm(d, g))

            # now we want to set p = g, but that is already true, so we're done

        else
            # from this point on we will only need p_i in the term p_i-p_c.
            # so we reuse the vector p_i by computing p_i = p_i - p_c and then
            # just so we aren't confused we name that p_diff
            p_i .-= p_c
            p_diff = p_i

            # Compute the optimal point on dogleg path
            b = 2 * wdot(d, p_c, d, p_diff)
            a = wnorm(d, p_diff)^2
            tau = (-b + sqrt(b^2 - 4a*(wnorm(d, p_c)^2 - delta^2)))/(2a)
            p_c .+= tau .* p_diff
            copyto!(p, p_c)
        end
    end
end

function trust_region_(df::OnceDifferentiable,
                       initial_x::AbstractArray{T},
                       xtol::Real,
                       ftol::Real,
                       iterations::Integer,
                       store_trace::Bool,
                       show_trace::Bool,
                       extended_trace::Bool,
                       factor::Real,
                       autoscale::Bool,
                       cache = NewtonTrustRegionCache(df)) where T
    copyto!(cache.x, initial_x)
    value_jacobian!!(df, cache.x)
    cache.r .= value(df)
    check_isfinite(cache.r)

    it = 0
    x_converged, f_converged, converged = assess_convergence(value(df), ftol)
    delta = convert(real(T), NaN)
    rho = convert(real(T), NaN)
    if converged
        tr = SolverTrace()
        name = "Trust-region with dogleg"
        if autoscale
            name *= " and autoscaling"
        end

        return SolverResults(name,
        #initial_x, reshape(cache.x, size(initial_x)...), norm(cache.r, Inf),
        initial_x, copy(cache.x), norm(cache.r, Inf),
        it, x_converged, xtol, f_converged, ftol, tr,
        first(df.f_calls), first(df.df_calls))
    end

    tr = SolverTrace()
    tracing = store_trace || show_trace || extended_trace
    @trustregiontrace convert(real(T), NaN)
    nn = length(cache.x)
    if autoscale
        for j = 1:nn
            cache.d[j] = norm(view(jacobian(df), :, j))
            if cache.d[j] == zero(cache.d[j])
                cache.d[j] = one(cache.d[j])
            end
        end
    else
        fill!(cache.d, one(real(T)))
    end

    delta = factor * wnorm(cache.d, cache.x)
    if delta == zero(delta)
        delta = factor
    end

    eta = convert(real(T), 1e-4)

    while !converged && it < iterations
        it += 1

        # Compute proposed iteration step
        dogleg!(cache.p, cache.p_c, cache.pi, cache.r, cache.d, jacobian(df), delta)

        copyto!(cache.xold, cache.x)
        cache.x .+= cache.p
        value!(df, cache.x)

        # Ratio of actual to predicted reduction (equation 11.47 in N&W)
        mul!(vec(cache.r_predict), jacobian(df), vec(cache.p))
        cache.r_predict .+= cache.r
        rho = (sum(abs2, cache.r) - sum(abs2, value(df))) / (sum(abs2, cache.r) - sum(abs2, cache.r_predict))

        if rho > eta
            # Successful iteration
            cache.r .= value(df)
            jacobian!(df, cache.x)

            # Update scaling vector
            if autoscale
                for j = 1:nn
                    cache.d[j] = max(convert(real(T), 0.1) * real(cache.d[j]), norm(view(jacobian(df), :, j)))
                end
            end

            x_converged, f_converged, converged = assess_convergence(cache.x, cache.xold, cache.r, xtol, ftol)
        else
            cache.x .-= cache.p
            x_converged, converged = false, false
        end

        @trustregiontrace euclidean(cache.x, cache.xold)

        # Update size of trust region
        if rho < 0.1
            delta = delta/2
        elseif rho >= 0.9
            delta = 2 * wnorm(cache.d, cache.p)
        elseif rho >= 0.5
            delta = max(delta, 2 * wnorm(cache.d, cache.p))
        end
    end

    name = "Trust-region with dogleg"
    if autoscale
        name *= " and autoscaling"
    end
    return SolverResults(name,
                         initial_x, copy(cache.x), norm(cache.r, Inf),
                         it, x_converged, xtol, f_converged, ftol, tr,
                         first(df.f_calls), first(df.df_calls))
end

function trust_region(df::OnceDifferentiable,
                      initial_x::AbstractArray{T},
                      xtol::Real,
                      ftol::Real,
                      iterations::Integer,
                      store_trace::Bool,
                      show_trace::Bool,
                      extended_trace::Bool,
                      factor::Real,
                      autoscale::Bool,
                      cache = NewtonTrustRegionCache(df)) where T
    trust_region_(df, initial_x, convert(real(T), xtol), convert(real(T), ftol), iterations, store_trace, show_trace, extended_trace, convert(real(T), factor), autoscale, cache)
end


text:
"""
	random_density(n)

Returns a random ``N × N`` density matrix, distributed according to the
Hilbert-Schmidt measure.

NOTE: No guarantee about the quality of the algorithm
For more details see [V Al Osipov et al 2010 J. Phys. A: Math. Theor. 43 055302]
(https://doi.org/10.1088/1751-8113/43/5/055302)
"""
function random_density_matrix(n)
    A = randn(n,n) + 1im * randn(n,n)
    rho = A * A'
    return rho / tr(rho)
end

"""
	unitary_op(H, t)

Returns the unitary evolution operator ``U = \\exp[-i H t]``.
"""
function unitary_op(H, t)
   return exp(-1im * Matrix(H) * t)
end

"""
    commutator(A, B)

Returns the commutator [A,B] = A * B - B * A between two matrices
"""
function commutator(a, b)
	a * b - b * a
end

"""
	localized_state(N, i)

Returns the density operator for a state ``|i⟩`` in an ``N``-dimensional Hilbert
space
"""
function localized_state(n, i)
   sparse([i,n],[i,n],[1,0])
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using Random

f(x) = string(first(keys(x))) * " => " * string(first(values(x)))

function test_namedtuple(n)
    S = String[]
    for _ in 1:n
        symb = Symbol(randstring(rand(3:10)))
        str = randstring(rand(3:10))
        nt = NamedTuple{(symb,)}((str,))
        push!(S, f(nt))
    end
    S
end

function test_dict(n)
    S = String[]
    for _ in 1:n
        symb = Symbol(randstring(rand(3:10)))
        str = randstring(rand(3:10))
        dic = Dict(symb => str)
        push!(S, f(dic))
    end
    S
end

@show test_namedtuple(5)
@show test_dict(5)
println()
print("test_namedtuple(2^8):")
@time test_namedtuple(2^8)
print("test_dict(2^8):      ")
@time test_dict(2^8);

# %%

max_stars_repo_path: 0009/namedtuple causes a large amount of specialization.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public