text:
"given delta, computes the total power response"
function comp_pg_response_total(network, gens::Set{Int}; delta=network["delta"])
    total_pg = 0.0

    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            if gen["index"] in gens
                pg = gen["pg_base"] + delta*gen["alpha"]
                pg = min(gen["pmax"], pg)
                pg = max(gen["pmin"], pg)
                total_pg += pg
            else
                total_pg += gen["pg_base"]
            end
        end
    end

    return total_pg
end

apply_pg_response!(network, pg_delta::Real) = apply_pg_response!(network, network["response_gens"], pg_delta::Real)

function apply_pg_response!(network, gens::Set{Int}, pg_delta::Real)
    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)")
    status = 0

    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"] && gen["index"] in gens
                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)")
            status = 1
            break
        end

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

        for (i,gen) in network["gen"]
            if gen["gen_status"] != 0 && gen["index"] in gens
                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=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

    alpha_final = 0.0
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0 && !gen["pg_fixed"] && gen["index"] in gens
            alpha_final += gen["alpha"]
        end
    end
    if isapprox(alpha_final, 0.0)
        warn(_LOGGER, "no remaining alpha for generator response (final alpha value $(alpha_final))")
        debug(_LOGGER, "delta $(delta_est)")
    end

    network["delta"] = delta_est
    return status
end

"fixes solution degeneracy issues when qg is a free variable, as is the case in PowerFlow"
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

gen_default = Dict(
#    "gen_bus" => 0,
#    "index" => 0,
    "pg" => 0.0,
    "model" => 1,
    "pg_base" => 0.0,
    "qg" => 0.0,
    "pg_fixed" => false,
    "pmax" => 0.94379,
    "pg_start" => 0.0,
    "qg_start" => 0.0,
    "alpha" => 0.0,
    "cost" => [0.0, 0.0, 99.99, 0.0],
    "gen_status" => 1,
    "qmax" =>  10.0,
    "qmin" => -10.0,
    "qg_fixed" => false,
    "pmin" => 0.0,
    "ncost" => 2,
    "virtual" => true
)

"""
A power flow solver inspired by the ARPA-e GOC Challenge 1 contingency-stage
specification but designed to be faster on large network cases.
Instead conducting PV/PQ bus switching address disjunctive constraints, this
solver simply adds extra reactive capability on buses where the voltage
bounds cannot be enforced and takes an associated power balance penalty.
"""
function run_fixpoint_pf_bqv!(network, pg_lost, solver; iteration_limit=typemax(Int64))
    time_start = time()

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

    #info(_LOGGER, "pg lost: $(pg_lost)")
    #info(_LOGGER, "delta: $(network["delta"])")
    #info(_LOGGER, "pre-solve time: $(time() - time_start)")

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

    gen_idx_max = maximum(gen["index"] for (i,gen) in network["gen"])
    gen_idx = trunc(Int, 10^ceil(log10(gen_idx_max)))

    bus_gens = gens_by_bus(network)

    iteration = 1
    vm_fixed = true
    while vm_fixed && iteration < iteration_limit
        info(_LOGGER, "pf soft fixpoint iteration: $iteration")

        time_start = time()
        result = run_pf_bqv_acr(network, solver, solution_processors=[_PM.sol_data_model!])
        info(_LOGGER, "solve pf time: $(time() - time_start)")

        if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
            correct_qg!(network, result["solution"], bus_gens=bus_gens)
            _PM.update_data!(network, result["solution"])
        else
            warn(_LOGGER, "solve issue with run_pf_bqv_acr, $(result["termination_status"])")
            break
        end

        vm_fixed = false
        for (i,bus) in network["bus"]
            if bus["vm"] < bus["vmin"] - vm_bound_tol || bus["vm"] > bus["vmax"] + vm_bound_tol
                active_gens = 0
                if length(bus_gens[i]) > 0
                    active_gens = sum(gen["gen_status"] != 0 for gen in bus_gens[i])
                end
                @assert(active_gens == 0)

                bus["vm_fixed"] = true
                current_vm = bus["vm"]
                if bus["vm"] < bus["vmin"]
                    bus["vm_base"] = bus["vmin"]
                    bus["vm"] = bus["vmin"]
                end
                if bus["vm"] > bus["vmax"]
                    bus["vm_base"] = bus["vmax"]
                    bus["vm"] = bus["vmax"]
                end

                warn(_LOGGER, "bus $(i) voltage out of bounds $(current_vm) -> $(bus["vm"]), adding virtual generator $(gen_idx)")
                gen_virtual = deepcopy(gen_default)
                gen_virtual["index"] = gen_idx
                gen_virtual["gen_bus"] = bus["index"]
                push!(bus_gens[i], gen_virtual)

                @assert(!haskey(network["gen"], "$(gen_idx)"))
                network["gen"]["$(gen_idx)"] = gen_virtual
                gen_idx += 1

                vm_fixed = true
            end
        end
        iteration += 1
    end

    if iteration >= iteration_limit
        warn(_LOGGER, "hit iteration limit")
    end

    for (i,gen) in collect(network["gen"])
        if haskey(gen, "virtual") && gen["virtual"]
            delete!(network["gen"], i)
            delete!(result["solution"]["gen"], i)
        end
    end

    return result
end

""
function run_pf_bqv_acr(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACRPowerModel, solver, build_pf_bqv_acr; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_pf_bqv_acr(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power_real(pm, bounded=false)
    _PM.variable_gen_power_imaginary(pm, bounded=false)

    delta = ref(pm, :delta)
    var(pm)[:delta] = @variable(pm.model, base_name="delta", start=0.0)
    sol(pm)[:delta] = var(pm)[:delta]
    #@constraint(pm.model, var(pm, :delta) == 0.0)

    shunt_values = Dict(sid => ref(pm, :shunt, sid)["bs"] for sid in ids(pm, :shunt_var))
    _PM.sol_component_value(pm, nw_id_default, :shunt, :bs, ids(pm, :shunt_var), shunt_values)

    var(pm)[:p_slack] = @variable(pm.model, p_slack, base_name="p_slack", start=0.0)

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

    _PM.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)
        _PM.constraint_theta_ref(pm, i)
    end
    #Memento.info(_LOGGER, "misc constraints time: $(time() - start_time)")

    start_time = time()
    for i in ids(pm, :branch)
        expression_branch_power_ohms_yt_from_goc(pm, i)
        _PM.expression_branch_power_ohms_yt_to(pm, i)
    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)
        pg[i] = gen["pg_base"]
        @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"])

        qg[i] = var(pm, :qg, i)
    end
    #Memento.info(_LOGGER, "gen expr time: $(time() - start_time)")

    start_time = time()
    p = var(pm, :p)
    q = var(pm, :q)
    for (i,bus) in ref(pm, :bus)
        #_PM.constraint_power_balance(pm, i)

        bus_arcs = ref(pm, :bus_arcs, 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)

        @constraint(pm.model, p_slack + 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

"""
A power flow solver conforming to the ARPA-e GOC Challenge 1 contingency-stage
specification.  The solver conducts multiple rounds of PV/PQ bus switching
to heuristically enforce the disjunctive constraints.
"""
function run_fixpoint_pf_pvpq!(network, pg_lost, solver; iteration_limit=typemax(Int64))
    time_start = time()

    network_backup = deepcopy(network)

    response_status = apply_pg_response!(network, pg_lost)

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

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

    bus_gens = gens_by_bus(network)

    active_response_gens = []
    for i in network["response_gens"]
        gen = network["gen"]["$(i)"]
        if gen["gen_status"] != 0
            push!(active_response_gens, gen)
        end
    end

    #network["delta_start"] = network["delta"]
    for (i,gen) in network["gen"]
        gen["qg_fixed"] = false
        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

    pf_fixed_all = all(gen["pg_fixed"] for gen in active_response_gens)

    cont_pf_failed = false
    time_start = time()
    #result = run_fixed_pf_nbf_rect(network, solver)
    if !pf_fixed_all
        result = run_pf_fixed_acr(network, solver, solution_processors=[_PM.sol_data_model!])
    else
        result = run_pf_fixed_bp_slack_acr(network, solver, solution_processors=[_PM.sol_data_model!])
    end
    debug(_LOGGER, "pf solve time: $(time() - time_start)")
    if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
        correct_qg!(network, result["solution"], bus_gens=bus_gens)
        _PM.update_data!(network, result["solution"])
        final_result = result
    else
        warn(_LOGGER, "$(network["cont_label"]) contingency pf solver FAILED with status $(result["termination_status"]) on iteration 0")
        cont_pf_failed = true
    end

    pg_switched = true
    qg_switched = true
    vm_switched = true

    iteration = 1
    deltas = [result["solution"]["delta"]]
    while (pg_switched || qg_switched || vm_switched) && !cont_pf_failed && 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"]
            if gen["index"] in network["response_gens"]
                pg = gen["pg_base"] + network["delta"]*gen["alpha"]
            else
                pg = gen["pg_base"]
            end

            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

                #=
                leads to solver infeasiblity
                for gen in bus_gens[i]
                    if isapprox(gen["qmin"],gen["qmax"])
                        gen["qg_fixed"] = true
                        gen["qg"] = gen["qmin"]
                        gen["qg_start"] = gen["qg"]
                    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
            debug(_LOGGER, "bus or gen switched: $iteration")
            time_start = time()

            pf_fixed_all = all(gen["pg_fixed"] for gen in active_response_gens)

            #result = run_fixed_pf_nbf_rect(network, solver, solution_processors=[_PM.sol_data_model!])
            #result = run_pf_fixed_acr(network, solver, solution_processors=[_PM.sol_data_model!])
            if !pf_fixed_all
                result = run_pf_fixed_acr(network, solver, solution_processors=[_PM.sol_data_model!])
            else
                result = run_pf_fixed_bp_slack_acr(network, solver, solution_processors=[_PM.sol_data_model!])
            end
            debug(_LOGGER, "pf solve time: $(time() - time_start)")
            if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
                correct_qg!(network, result["solution"], bus_gens=bus_gens)
                _PM.update_data!(network, result["solution"])
                final_result = result
            else
                warn(_LOGGER, "$(network["cont_label"]) contingency pf solver FAILED with status $(result["termination_status"]) on iteration 0")
                break
            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

    vm_bound_vio = false
    for (i,bus) in network["bus"]
        if bus["bus_type"] != 4
            bus_sol = final_result["solution"]["bus"][i]
            if bus_sol["vm"] - vm_bound_tol >= bus["vmax"] || bus_sol["vm"] + vm_bound_tol <= bus["vmin"]
                vm_bound_vio = true
                warn(_LOGGER, "$(network["cont_label"]) vm bound out of range on bus $(i): $(bus["vmin"]) - $(bus_sol["vm"]) - $(bus["vmax"])")
            end
        end
    end

    qg_bound_vio = false
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            gen_sol = final_result["solution"]["gen"][i]
            if gen_sol["qg"] - qg_bound_tol >= gen["qmax"] || gen_sol["qg"] + qg_bound_tol <= gen["qmin"]
                qg_bound_vio = true
                warn(_LOGGER, "$(network["cont_label"]) qg bound out of range on gen $(i): $(gen["qmin"]) - $(gen_sol["qg"]) - $(gen["qmax"])")
            end
        end
    end

    if vm_bound_vio || qg_bound_vio
        warn(_LOGGER, "$(network["cont_label"]) running voltage profile correction")
        result = run_fixpoint_opf!(network_backup, pg_lost, solver, iteration_limit=iteration_limit)
        if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
            _PM.update_data!(network, result["solution"])
            final_result = result
        else
            warn(_LOGGER, "$(network["cont_label"]) voltage profile correction solver FAILED with status $(result["termination_status"])")
        end
    end

    return final_result
end

""
function run_pf_fixed_acr(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACRPowerModel, solver, build_pf_fixed_acr; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_pf_fixed_acr(pm::_PM.AbstractPowerModel)
    start_time = time()
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power_real(pm, bounded=false)
    _PM.variable_gen_power_imaginary(pm, bounded=false)
    #_PM.variable_branch_power(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)")
    sol(pm)[:delta] = var(pm)[:delta]

    shunt_values = Dict(sid => ref(pm, :shunt, sid)["bs"] for sid in ids(pm, :shunt_var))
    _PM.sol_component_value(pm, nw_id_default, :shunt, :bs, ids(pm, :shunt_var), shunt_values)

    start_time = time()

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

    _PM.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)
        _PM.constraint_theta_ref(pm, i)
    end
    #Memento.info(_LOGGER, "misc constraints time: $(time() - start_time)")

    start_time = time()
    for i in ids(pm, :branch)
        expression_branch_power_ohms_yt_from_goc(pm, i)
        _PM.expression_branch_power_ohms_yt_to(pm, i)
    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"]
            @constraint(pm.model, var(pm, :pg, i) == gen["pg"])
            pg[i] = gen["pg"]
        else
            if i in ref(pm, :response_gens)
                @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"] + gen["alpha"]*delta)
                pg[i] = gen["pg_base"] + gen["alpha"]*delta
            else
                @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"])
                pg[i] = gen["pg_base"]
            end
        end

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

    start_time = time()
    for (i,bus) in ref(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end
    #Memento.info(_LOGGER, "power balance constraint time: $(time() - start_time)")
end

"a variant of fixed_pf_nbf_rect2 with a distributed active power slack"
function run_pf_fixed_bp_slack_acr(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACRPowerModel, solver, build_pf_fixed_bp_slack_acr; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_pf_fixed_bp_slack_acr(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power_real(pm, bounded=false)
    _PM.variable_gen_power_imaginary(pm, bounded=false)
    #_PM.variable_branch_power(pm, bounded=false)

    delta = ref(pm, :delta)
    var(pm)[:delta] = @variable(pm.model, base_name="delta", start=0.0)
    sol(pm)[:delta] = var(pm)[:delta]
    @constraint(pm.model, var(pm, :delta) == delta)

    shunt_values = Dict(sid => ref(pm, :shunt, sid)["bs"] for sid in ids(pm, :shunt_var))
    _PM.sol_component_value(pm, nw_id_default, :shunt, :bs, ids(pm, :shunt_var), shunt_values)

    active_response_gens = intersect(ids(pm, :gen), ref(pm, :response_gens))
    # for i in active_response_gens
    #     gen = ref(pm, :gen, i)
    #     println("$(i) $(gen["pg_fixed"]) $(gen["pmin"]) $(gen["pg_base"]) $(gen["pg"]) $(gen["pmax"])")
    # end

    if !all(ref(pm, :gen, i)["pg_fixed"] for i in active_response_gens)
        Memento.error(_LOGGER, "fixed_pf_nbf_rect2_ds model requires all response_gens have pg_fixed set to true")
    end
    var(pm)[:p_slack] = @variable(pm.model, p_slack, base_name="p_slack", start=0.0)

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

    _PM.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)
        _PM.constraint_theta_ref(pm, i)
    end
    #Memento.info(_LOGGER, "misc constraints time: $(time() - start_time)")

    start_time = time()
    for i in ids(pm, :branch)
        expression_branch_power_ohms_yt_from_goc(pm, i)
        _PM.expression_branch_power_ohms_yt_to(pm, i)
    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
            if i in ref(pm, :response_gens)
                pg[i] = gen["pg_base"] + gen["alpha"]*delta
                @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"] + gen["alpha"]*delta)
            else
                pg[i] = gen["pg_base"]
                @constraint(pm.model, var(pm, :pg, i) == gen["pg_base"])
            end
        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()
    p = var(pm, :p)
    q = var(pm, :q)
    for (i,bus) in ref(pm, :bus)
        #_PM.constraint_power_balance(pm, i)

        bus_arcs = ref(pm, :bus_arcs, 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, p_slack + 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

function run_fixpoint_opf!(network, pg_lost, solver; iteration_limit=typemax(Int64))
    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" => _PM.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
        gen["pg_start"] = gen["pg"]
        gen["qg_start"] = gen["qg"]
    end

    for (i,bus) in network["bus"]
        if length(bus_gens[i]) == 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_opf_contingency_acp(network, solver, solution_processors=[_PM.sol_data_model!])
    debug(_LOGGER, "pf solve time: $(time() - time_start)")
    if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
        warn(_LOGGER, "$(network["cont_label"]) voltage profile correction objective: $(result["objective"])")
        _PM.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

    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"]

            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 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
            debug(_LOGGER, "bus or gen swtiched: $iteration")
            time_start = time()
            #result = run_fixed_pf_nbf_rect(network, solver, solution_processors=[_PM.sol_data_model!])
            result = run_pf_fixed_acr(network, solver, solution_processors=[_PM.sol_data_model!])
            debug(_LOGGER, "pf solve time: $(time() - time_start)")
            if result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED
                correct_qg!(network, result["solution"], bus_gens=bus_gens)
                _PM.update_data!(network, result["solution"])
                final_result = result
            else
                warn(_LOGGER, "contingency pf solver FAILED with status $(result["termination_status"]) on iteration 0")
                break
            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

# this approach has a minor ipopt runtime performance hit, but extracting branch flow values is much faster

"attempts to correct voltage profile only"
function run_opf_contingency_acp(file, solver; kwargs...)
    return _PM.run_model(file, _PM.ACPPowerModel, solver, build_opf_contingency; ref_extensions=[ref_add_goc!], kwargs...)
end

""
function build_opf_contingency(pm::_PM.AbstractPowerModel)
    _PM.variable_bus_voltage(pm, bounded=false)
    _PM.variable_gen_power_real(pm, bounded=false)
    _PM.variable_gen_power_imaginary(pm, bounded=false)
    #_PM.variable_branch_power(pm, bounded=false)

    variable_shunt_admittance_imaginary(pm)

    qg_vio = @variable(pm.model,
        [i in ids(pm, :gen)], base_name="qg_vio",
        lower_bound = 0.0,
        start = 0.0
    )

    vm_vio = @variable(pm.model,
        [i in ids(pm, :bus)], base_name="vm_vio",
        lower_bound = 0.0,
        start = 0.0
    )

    var(pm)[:delta] = @variable(pm.model, delta, base_name="delta", start = 0.0)
    sol(pm)[:delta] = var(pm)[:delta]
    if haskey(pm.ref, :delta_start)
        set_start_value(delta, ref(pm, :delta_start))
    end

    #@objective(pm.model, Min,
    #    sum( 5e5*qg_vio[i]^2 for (i,gen) in ref(pm, :gen) ) +
    #    sum( 1e7*vm_vio[i]^2 for (i,bus) in ref(pm, :bus) )
    #)

    @objective(pm.model, Min,
        sum( 5e5*qg_vio[i] for (i,gen) in ref(pm, :gen) ) +
        sum( 1e7*vm_vio[i] for (i,bus) in ref(pm, :bus) )
    )

    _PM.constraint_model_voltage(pm)

    vm = Dict{Int,Any}()
    for (i,bus) in ref(pm, :bus)
        if bus["vm_fixed"]
            @constraint(pm.model, var(pm, :vm, i) == bus["vm_base"])
            @constraint(pm.model, vm_vio[i] == 0)
            vm[i] = bus["vm_base"]
        else
            @constraint(pm.model,  vm_vio[i] >=  var(pm, :vm, i) - bus["vmax"])
            @constraint(pm.model,  vm_vio[i] >= -var(pm, :vm, i) + bus["vmin"])
            vm[i] = var(pm, :vm, i)
        end
    end
    var(pm)[:vm] = vm

    for i in ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    start_time = time()
    for i in ids(pm, :branch)
        expression_branch_power_ohms_yt_from_goc(pm, i)
        _PM.expression_branch_power_ohms_yt_to(pm, i)
    end
    #Memento.info(_LOGGER, "flow expr time: $(time() - start_time)")

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

        if gen["qg_fixed"]
            @constraint(pm.model, var(pm, :qg, i) == gen["qg"])
            @constraint(pm.model, qg_vio[i] == 0)
            qg[i] = gen["qg"]
        else
            @constraint(pm.model,  qg_vio[i] >=  var(pm, :qg, i) - gen["qmax"])
            @constraint(pm.model,  qg_vio[i] >= -var(pm, :qg, i) + gen["qmin"])
            qg[i] = var(pm, :qg, i)
        end
    end
    var(pm)[:pg] = pg
    var(pm)[:qg] = qg

    for (i,bus) in ref(pm, :bus)
        constraint_power_balance_shunt_dispatch(pm, i)
    end

end


text:
using KernelAbstractions
using CUDAKernels
using CUDA

nij = 14
nn1 = 1000
nn0 = 100

iet⁺⁺ = rand(nij,nij,nn1,nn0);
ier⁻⁺ = rand(nij,nij,nn1,nn0);

ϖ_λ₀λ₁ = rand(nn1,nn0);
dτ₀= rand(1)[1];
dτ₁= rand(1)[1];
dτ_λ= rand(nn1);
Z⁻⁺_λ₀λ₁= rand(nij,nij);
Z⁺⁺_λ₀λ₁= rand(nn1,nn0);
qp_μN= rand(nij);
wct2 = rand(nij);
ϖ_λ = rand(nn1);

@kernel function get_elem_rt!(ier⁻⁺, iet⁺⁺, ϖ_λ, ϖ_λ₀λ₁, dτ₀, dτ₁, dτ_λ, Z⁻⁺_λ₀λ₁, Z⁺⁺_λ₀λ₁, qp_μN, wct2)
    i, j, n₁, n₀ = @index(Global, NTuple)
    # let n₁ cover the full range of wavelengths, while n₀ only includes wavelengths at intervals
    # that contribute significantly enough to inelastic scattering, so that n₀≪n₁
    if (wct2[j]>1.e-8)
        # dτ₀, dτ₁ are the purely scattering (elastic+inelastic) molecular elemental
        # optical thicknesses at wavelengths λ₀ and λ₁
        # 𝐑⁻⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁻⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ+μⱼ)) ̇(1 - exp{-τ ̇(1/μᵢ + 1/μⱼ)}) ̇𝑤ⱼ
        ier⁻⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * (dτ₀/dτ₁) * Z⁻⁺_λ₀λ₁[i,j] * (qp_μN[j]*dτ₁ / (qp_μN[i]*dτ₀ + qp_μN[j]*dτ₁)) * (1 - exp(-((dτ_λ[n₁] / qp_μN[i]) + (dτ_λ[n₀] / qp_μN[j])))) * (wct2[j])

        if (qp_μN[i] == qp_μN[j])
            # @show i,j
            # 𝐓⁺⁺(μᵢ, μᵢ) = (exp{-τ/μᵢ} + ϖ ̇𝐙⁺⁺(μᵢ, μᵢ) ̇(τ/μᵢ) ̇exp{-τ/μᵢ}) ̇𝑤ᵢ
            if i == j
                if abs(dτ_λ[n₀]-dτ_λ[n₁])>1.e-6
                    iet⁺⁺[i,j,n₁,n₀] = ((exp(-dτ_λ[n₀] / qp_μN[i]) - exp(-dτ_λ[n₁] / qp_μN[i]))/(dτ_λ[n₁]-dτ_λ[n₀])) * ϖ_λ₀λ₁[n₁,n₀] * dτ₀ * Z⁺⁺_λ₀λ₁[i,i] * wct2[i]
                else
                    iet⁺⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * dτ₀ * Z⁺⁺_λ₀λ₁[i,i] * wct2[i] * exp(-dτ_λ[n₀] / qp_μN[j])/ qp_μN[j]
                end
            else
                iet⁺⁺[i,j,n₁,n₀] = 0.0
            end
        else
            #@show  qp_μN[i], qp_μN[j]
            # 𝐓⁺⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁺⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ-μⱼ)) ̇(exp{-τ/μᵢ} - exp{-τ/μⱼ}) ̇𝑤ⱼ
            # (𝑖 ≠ 𝑗)
            iet⁺⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * (dτ₀/dτ₁) * Z⁺⁺_λ₀λ₁[i,j] * (qp_μN[j]*dτ₁ / (qp_μN[i]*dτ₀ - qp_μN[j]*dτ₁)) * (exp(-dτ_λ[n₁] / qp_μN[i]) - exp(-dτ_λ[n₀] / qp_μN[j])) * wct2[j]
        end
    else
        ier⁻⁺[i,j,n₁,n₀] = 0.0
        if i==j
            iet⁺⁺[i,j,n₁,n₀] = 0.0
        else
            iet⁺⁺[i,j,n₁,n₀] = 0.0
        end
    end
end

function get_elem_rt!(ier⁻⁺, iet⁺⁺, ϖ_λ, ϖ_λ₀λ₁, dτ₀, dτ₁, dτ_λ, Z⁻⁺_λ₀λ₁, Z⁺⁺_λ₀λ₁, qp_μN, wct2)
    d1,d2,d3,d4 = size(ier⁻⁺)
    for i in 1:d1, j in 1:d2, n₁ in 1:d3, n₀ in 1:d4
        #i, j, n₁, n₀ = @index(Global, NTuple)
        # let n₁ cover the full range of wavelengths, while n₀ only includes wavelengths at intervals
        # that contribute significantly enough to inelastic scattering, so that n₀≪n₁
        if (wct2[j]>1.e-8)
            # dτ₀, dτ₁ are the purely scattering (elastic+inelastic) molecular elemental
            # optical thicknesses at wavelengths λ₀ and λ₁
            # 𝐑⁻⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁻⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ+μⱼ)) ̇(1 - exp{-τ ̇(1/μᵢ + 1/μⱼ)}) ̇𝑤ⱼ
            ier⁻⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * (dτ₀/dτ₁) * Z⁻⁺_λ₀λ₁[i,j] * (qp_μN[j]*dτ₁ / (qp_μN[i]*dτ₀ + qp_μN[j]*dτ₁)) * (1 - exp(-((dτ_λ[n₁] / qp_μN[i]) + (dτ_λ[n₀] / qp_μN[j])))) * (wct2[j])

            if (qp_μN[i] == qp_μN[j])

                # 𝐓⁺⁺(μᵢ, μᵢ) = (exp{-τ/μᵢ} + ϖ ̇𝐙⁺⁺(μᵢ, μᵢ) ̇(τ/μᵢ) ̇exp{-τ/μᵢ}) ̇𝑤ᵢ
                if i == j
                    if abs(dτ_λ[n₀]-dτ_λ[n₁])>1.e-6
                        iet⁺⁺[i,j,n₁,n₀] = ((exp(-dτ_λ[n₀] / qp_μN[i]) - exp(-dτ_λ[n₁] / qp_μN[i]))/(dτ_λ[n₁]-dτ_λ[n₀])) * ϖ_λ₀λ₁[n₁,n₀] * dτ₀ * Z⁺⁺_λ₀λ₁[i,i] * wct2[i]
                    else
                        iet⁺⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * dτ₀ * Z⁺⁺_λ₀λ₁[i,i] * wct2[i] * exp(-dτ_λ[n₀] / qp_μN[j])/ qp_μN[j]
                    end
                else
                    iet⁺⁺[i,j,n₁,n₀] = 0.0
                end
            else

                # 𝐓⁺⁺(μᵢ, μⱼ) = ϖ ̇𝐙⁺⁺(μᵢ, μⱼ) ̇(μⱼ/(μᵢ-μⱼ)) ̇(exp{-τ/μᵢ} - exp{-τ/μⱼ}) ̇𝑤ⱼ
                # (𝑖 ≠ 𝑗)
                iet⁺⁺[i,j,n₁,n₀] = ϖ_λ₀λ₁[n₁,n₀] * (dτ₀/dτ₁) * Z⁺⁺_λ₀λ₁[i,j] * (qp_μN[j]*dτ₁ / (qp_μN[i]*dτ₀ - qp_μN[j]*dτ₁)) * (exp(-dτ_λ[n₁] / qp_μN[i]) - exp(-dτ_λ[n₀] / qp_μN[j])) * wct2[j]
            end
        else
            ier⁻⁺[i,j,n₁,n₀] = 0.0
            if i==j
                iet⁺⁺[i,j,n₁,n₀] = 0.0
            else
                iet⁺⁺[i,j,n₁,n₀] = 0.0
            end
        end
    end
end

# Test w/o kernel:
get_elem_rt!(ier⁻⁺, iet⁺⁺, ϖ_λ, ϖ_λ₀λ₁,dτ₀, dτ₁, dτ_λ, Z⁻⁺_λ₀λ₁, Z⁺⁺_λ₀λ₁, qp_μN, wct2);
base_ier⁻⁺ = deepcopy(ier⁻⁺);
base_iet⁺⁺ = deepcopy(iet⁺⁺);

#ier⁻⁺ .= 0;
#iet⁺⁺ .= 0;
# Test CPU kernel version:
device = CPU()
kernel! = get_elem_rt!(device)
event = kernel!(ier⁻⁺, iet⁺⁺, ϖ_λ, ϖ_λ₀λ₁,dτ₀, dτ₁, dτ_λ, Z⁻⁺_λ₀λ₁, Z⁺⁺_λ₀λ₁, qp_μN, wct2, ndrange=size(ier⁻⁺));
wait(device, event)

base_ier⁻⁺ ≈ ier⁻⁺
base_iet⁺⁺ ≈ iet⁺⁺

# Test GPU kernel version:
#ier⁻⁺ .= 0;
#iet⁺⁺ .= 0;
if has_cuda()
    c_ϖ_λ₀λ₁ = CuArray(ϖ_λ₀λ₁);
    #dτ₀= rand(1)[1];
    #dτ₁= rand(1)[1];
    c_dτ_λ= CuArray(dτ_λ);
    c_Z⁻⁺_λ₀λ₁= CuArray(Z⁻⁺_λ₀λ₁);
    c_Z⁺⁺_λ₀λ₁= CuArray(Z⁺⁺_λ₀λ₁);
    c_qp_μN = CuArray(qp_μN);
    c_wct2 = CuArray(wct2);
    c_ϖ_λ = CuArray(ϖ_λ);
    c_iet⁺⁺ = CuArray(iet⁺⁺);
    c_ier⁻⁺ = CuArray(ier⁻⁺);

    device = CUDAKernels.CUDADevice()
    kernel! = get_elem_rt!(device)
    event = kernel!(c_ier⁻⁺, c_iet⁺⁺, c_ϖ_λ, c_ϖ_λ₀λ₁,dτ₀, dτ₁, c_dτ_λ, c_Z⁻⁺_λ₀λ₁, c_Z⁺⁺_λ₀λ₁, c_qp_μN, c_wct2, ndrange=size(c_ier⁻⁺));
    wait(device, event)

    cuda_iet⁺⁺ = Array(c_iet⁺⁺);
    cuda_ier⁻⁺ = Array(c_ier⁻⁺);

    cuda_iet⁺⁺ ≈ base_iet⁺⁺
    cuda_ier⁻⁺ ≈ base_ier⁻⁺
end

ext: jl


text:
# Formatted Output (B.17)

export @cuprintf

"""
Print a formatted string in device context on the host standard output:

    @cuprintf("%Fmt", args...)

Note that this is not a fully C-compliant `printf` implementation; see the CUDA
documentation for supported options and inputs.

Also beware that it is an untyped, and unforgiving `printf` implementation. Type widths need
to match, eg. printing a 64-bit Julia integer requires the `%ld` formatting string.
"""
macro cuprintf(fmt::String, args...)
    fmt_val = Val(Symbol(fmt))
    return :(_cuprintf($fmt_val, $(map(esc, args)...)))
end

@generated function _cuprintf(::Val{fmt}, argspec...) where {fmt}
    arg_exprs = [:( argspec[$i] ) for i in 1:length(argspec)]
    arg_types = [argspec...]

    # TODO: needs to adhere to C vararg promotion (short -> int, float -> double, etc.)

    T_void = LLVM.VoidType(JuliaContext())
    T_int32 = LLVM.Int32Type(JuliaContext())
    T_pint8 = LLVM.PointerType(LLVM.Int8Type(JuliaContext()))

    # create functions
    param_types = LLVMType[convert.(LLVMType, arg_types)...]
    llvm_f, _ = create_function(T_int32, param_types)
    mod = LLVM.parent(llvm_f)

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        str = globalstring_ptr!(builder, String(fmt))

        # construct and fill args buffer
        if isempty(argspec)
            buffer = LLVM.PointerNull(T_pint8)
        else
            argtypes = LLVM.StructType("vprintf_args", JuliaContext())
            elements!(argtypes, param_types)

            args = alloca!(builder, argtypes, "args")
            for (i, param) in enumerate(parameters(llvm_f))
                p = struct_gep!(builder, args, i-1)
                store!(builder, param, p)
            end

            buffer = bitcast!(builder, args, T_pint8)
        end

        # invoke vprintf and return
        vprintf_typ = LLVM.FunctionType(T_int32, [T_pint8, T_pint8])
        vprintf = LLVM.Function(mod, "vprintf", vprintf_typ)
        chars = call!(builder, vprintf, [str, buffer])

        ret!(builder, chars)
    end

    arg_tuple = Expr(:tuple, arg_exprs...)
    call_function(llvm_f, Int32, Tuple{arg_types...}, arg_tuple)
end


text:
"""
    This example shows how to use Term's styled logging

It's as simple as calling one function, Term takes care of the rest
so you can use Julia's logging functionality as you would normally
"""

import Term: install_term_logger, uninstall_term_logger, Panel

# insall term's logger as the global logger
install_term_logger()

# log
@info "my log!"

"""
as you can see the output is similar to the default logging system, with just
a couple more details. However you can now add markdown style to your
log messages
"""

@warn "tell us if this was [bold red]undexpected![/bold red]"

"""
but there's more!
Just like in the normal logging system you can pass additional arguments to the logging macro.
Term's logger prints them with additional information, and some style of course.
"""

x = collect(1:2:20)
y = x * x'
name = "the name is {bold blue}Term{/bold blue}"
p1 = Panel("text")

print("\n"^3)
@error "{italic green bold}fancy logs!{/italic green bold}" x y name √9 install_term_logger p1

"""
You can reset Julia's default logger like this:
"""

uninstall_term_logger()
@error "{italic green bold}fancy logs!{/italic green bold}" x y name √9 install_term_logger p1


text:
module Solvers
using DataArrays

    function cg(A, b, x, maxiter, eps)
    # Conjugate Gradients method
    # INPUT: A, b, starting value x, max numer of iterations maxiter, error tolerance eps<1

        i           = 0;
        r           = b - A*x;
        d           = r;
        delta_new   = dot(r,r); #r'*r;
        delta_0     = delta_new;
        while i <= maxiter && delta_new[1] > (eps^2)*delta_0
            q       = A*d;
            alpha   = delta_new/dot(d,q) #(d'*q);
            x       = x + alpha*d;

            # Error
            if mod(i,50)==0    # if stucked somewhere
                r = b - A*x;
            else
                r = r - alpha*q;
            end

            # Update
            delta_old = delta_new;
            delta_new = dot(r,r) #r'*r;
            beta      = delta_new/delta_old;
            d         = r + beta*d;
            i         = i+1;

        end #while
        return x
    end #function

end


text:
using CSG
using Test

@testset "CSG.jl" begin
    # Write your tests here.
end


text:
using Dates

mutable struct ToDo
  id::Int32
  description::String
  completed::Bool
  created::Date
  priority::Int8
end

todo1 = ToDo(1, "Getting groceries", false, Date("2022-04-01", "yyyy-mm-dd"), 5)
todo2 = ToDo(2, "Preparing dinner", false, Date("2022-04-02", "yyyy-mm-dd"), 6)

# if:
if todo2.priority > todo1.priority
  println("Better do todo2 first")
else
  println("Better do todo1 first")
end
# "Better do todo2 first"


text:
# kruskal_wallis.jl
# Kruskal-Wallis rank sum test
#
# Copyright (C) 2014   Christoph Sawade
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export KruskalWallisTest

immutable KruskalWallisTest <: HypothesisTest
    n_i::Vector{Int}         # number of observations in each group
    df::Int                  # degrees of freedom
    R_i::Vector{Float64}     # rank sums
    H::Float64               # test statistic: chi-square statistic
    tie_adjustment::Float64  # adjustment for ties
end

function KruskalWallisTest{T<:Real}(groups::AbstractVector{T}...)
    (H, R_i, tieadj, n_i) = kwstats(groups...)
    if length(groups)<=3 && any(n_i .< 6)
        warn("This test is only asymptotically correct and might be inaccurate for the given group size")
    end
    df = length(groups) - 1
    KruskalWallisTest(n_i, df, R_i, H, tieadj)
end

testname(::KruskalWallisTest) = "Kruskal-Wallis rank sum test (chi-square approximation)"
population_param_of_interest(x::KruskalWallisTest) = ("Location parameters", "all equal", NaN) # parameter of interest: name, value under h0, point estimate

function show_params(io::IO, x::KruskalWallisTest, ident)
    println(io, ident, "number of observation in each group: ", x.n_i)
    println(io, ident, "χ²-statistic:                        ", x.H)
    println(io, ident, "rank sums:                           ", x.R_i)
    println(io, ident, "degrees of freedom:                  ", x.df)
    println(io, ident, "adjustment for ties:                 ", x.tie_adjustment)
end

pvalue(x::KruskalWallisTest) = pvalue(Chisq(x.df), x.H; tail=:right)

## helper

# Get H, rank sums, and tie adjustment for Kruskal-Wallis test
function kwstats{T<:Real}(groups::AbstractVector{T}...)
    n_i = [length(g) for g in groups]
    n = sum(n_i)

    # get ranks and adjustment for ties
    (ranks, tieadj) = tiedrank_adj([groups...;])
    C = 1-tieadj/(n^3 - n)

    # compute rank sums
    R_i = Array(Float64, length(groups))
    n_end = 0
    for i=1:length(groups)
        R_i[i] = sum(ranks[n_end+1:n_end+n_i[i]])
        n_end += n_i[i]
    end

    # compute test statistic and correct for ties
    H = 12 * sum(R_i.^2./n_i) / (n * (n + 1)) - 3 * (n + 1)
    H /= C

    (H, R_i, C, n_i)
end


text:
"""
    GenExpansionParams(h::Function, hj::Function, X::Vector{Interval{T}},
                       P::Vector{Interval{T}},pmid::Vector{T},mc_opts::mc_opts)

Generates relaxation of state variable at pmid in P and outputs results to use
as the parameters in subsequent calculations. Use a direct calculation of the
inverse function for preconditioning. Inputs are:
* `h::Function`: h(z,p) which implicit defines state function
* `hj::Function`: hj(z,p) w.r.t z which implicit defines state function
* `P::Vector{Interval}`: State variable bounds
* `X::Vector{Interval{T}}`: Decision variable bounds
* `pmid::Vector{T}`: Point at which to generate relaxation
* `mc_opts`: Options for generating implicit function relaxation
Outs the following the tuple `sto_out`:
* `sto_out::Vector{SMCg{np,T}}` - McCormick parameter object stack
--------------------------------------------------------------------------------
"""
function GenExpansionParams(h::Function, hj::Function,
                            X::Vector{Interval{T}},
                            P::Vector{Interval{T}},
                            pmid::Vector{T},mc_opts::mc_opts{T}) where {T<:AbstractFloat}

  nxi::Int64 = length(X)
  np::Int64 = length(P)
  szero::SVector{np,T} = @SVector zeros(np)
  sone::SVector{np,T} = @SVector ones(np)
  exp_opt::Vector{Any} = Any[nxi,np,mc_opts.lambda]
  sto_out::Vector{Vector{SMCg{np,T}}} = [[zero(SMCg{np,T}) for i=1:nxi] for j=1:(mc_opts.kmax+1)]

  x_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].lo,szero,szero,@interval(X[i].lo,X[i].hi),false,[∅],[1.0]) for i=1:nxi]
  xa_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].lo,X[i].lo,szero,szero,@interval(X[i].lo,X[i].lo),false,[∅],[1.0]) for i=1:nxi]
  xA_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].hi,szero,szero,@interval(X[i].hi,X[i].hi),false,[∅],[1.0]) for i=1:nxi]
  if mc_opts.z_rnd_all == true
    z_mc::Vector{SMCg{np,T}} = Rnd_Out_Z_All(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_all_eps)
  elseif mc_opts.z_rnd == true
    z_mc = Rnd_Out_Z_Intv(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_eps)
  else
    z_mc = mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc
  end

  p_mc::Vector{SMCg{np,T}} = [SMCg(pmid[i],pmid[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  pref_mc::Vector{SMCg{np,T}} = copy(p_mc)
  aff_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(xA_mc[i].cc,xa_mc[i].cv,szero,szero,@interval(xa_mc[i].cv,xA_mc[i].cc),false,[∅],[1.0]) for i=1:nxi]

  H_mc::Vector{SMCg{np,T}} = h(z_mc,p_mc)
  dH_mc::Union{Vector{SMCg{np,T}},Array{SMCg{np,T},2}} = hj(aff_mc,p_mc)

  Y = (nxi == 1) ? [one(T)/mid(dH_mc[1].Intv)] : mid.(Intv.(dH_mc))
  Precondition!(H_mc,dH_mc,Y,nxi)
  # stores some things for repeated us
  optc = Any[szero,sone]

  # Begins loop to generate parameters
  sto_out[1] = copy(x_mc)
  for k=1:mc_opts.kmax

    aff_mc = [SMCg{np,Float64}(x_mc[i].cc,x_mc[i].cv,szero,szero,
                   @interval(x_mc[i].cv,x_mc[i].cc),false,[∅],[1.0]) for i=1:nxi]
    if mc_opts.aff_rnd_all == true
      aff_mc = Rnd_Out_Z_All(aff_mc,mc_opts.aff_rnd_all_eps)
    elseif mc_opts.aff_rnd == true
      aff_mc = Rnd_Out_Z_Intv(aff_mc,mc_opts.aff_rnd_eps)
    end

    if mc_opts.hhj_rnd_all == true
      H_mc,dH_mc = Rnd_Out_H_All(h(z_mc,p_mc),hj(aff_mc,p_mc),mc_opts.hhj_rnd_all_eps)
    elseif mc_opts.hhj_rnd == true
      H_mc,dH_mc = Rnd_Out_H_Intv(h(z_mc,p_mc),hj(aff_mc,p_mc),mc_opts.hhj_rnd_eps)
    else
      H_mc,dH_mc = h(z_mc,p_mc),hj(aff_mc,p_mc)
    end

    Y = (nxi == 1) ? [one(T)/mid(dH_mc[1].Intv)] : mid.(Intv.(dH_mc))

    # applies preconditioner
    Precondition!(H_mc,dH_mc,Y,nxi)
    # applies parametric iteration
    if (mc_opts.style == "NewtonGS")
      MC_NewtonGS!(z_mc,x_mc,dH_mc,H_mc,nxi)
    elseif (mc_opts.style == "KrawczykCW")
      MC_KrawczykCW!(z_mc,x_mc,dH_mc,H_mc,nxi)
    else
      error("Unsupported Style of Implicit Relaxation")
    end

    # update affine relaxations & correct
    Affine_Exp!(x_mc,p_mc,p_mc,xa_mc,xA_mc,z_mc,exp_opt)

    Correct_Exp!(z_mc,x_mc,X,nxi,np,mc_opts.aff_correct_eps)

    if mc_opts.z_rnd_all == true
      z_mc = Rnd_Out_Z_All(z_mc,mc_opts.z_rnd_all_eps)
    elseif mc_opts.z_rnd == true
      z_mc = Rnd_Out_Z_Intv(z_mc,mc_opts.z_rnd_eps)
    end

    # store relaxation
    sto_out[k+1] = copy(x_mc)
    end
  #return sto_out,sto_z,sto_x
  return sto_out
end

"""
    InGenExpansionParams(h::Function, hj::Function, X::Vector{Interval{T}},
                       P::Vector{Interval{T}},pmid::Vector{T},mc_opts::mc_opts)

Generates relaxation of state variable at pmid in P and outputs results to use
as the parameters in subsequent calculations. Use a sparse LDU factorization with
full pivoting to precondition and performs function evaluations in-place. Inputs are:
* `h::Function`: h!(H,z,p) which implicit defines state function
* `hj::Function`: hj!(H,z,p) w.r.t z which implicit defines state function
* `P::Vector{Interval}`: State variable bounds
* `X::Vector{Interval{T}}`: Decision variable bounds
* `pmid::Vector{T}`: Point at which to generate relaxation
* `mc_opts`: Options for generating implicit function relaxation
Outs the following the tuple `sto_out`:
* `sto_out::Vector{SMCg{np,T}}` - McCormick parameter object stack
--------------------------------------------------------------------------------
"""
function InGenExpansionParams(h!::Function, hj!::Function,
                            X::Vector{Interval{T}},
                            P::Vector{Interval{T}},
                            pmid::Vector{T},mc_opts::mc_opts{T}) where {T<:AbstractFloat}

  nx::Int64 = length(X)
  np::Int64 = length(P)
  szero::SVector{np,T} = @SVector zeros(np)
  sone::SVector{np,T} = @SVector ones(np)
  exp_opt::Vector{Any} = Any[nx,np,mc_opts.lambda]
  sto_out::Vector{Vector{SMCg{np,T}}} = [[zero(SMCg{np,T}) for i=1:nx] for j=1:(mc_opts.kmax+1)]

  x_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].lo,szero,szero,@interval(X[i].lo,X[i].hi),false,[∅],[1.0]) for i=1:nx]
  xa_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].lo,X[i].lo,szero,szero,@interval(X[i].lo,X[i].lo),false,[∅],[1.0]) for i=1:nx]
  xA_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].hi,szero,szero,@interval(X[i].hi,X[i].hi),false,[∅],[1.0]) for i=1:nx]

  if mc_opts.z_rnd_all == true
    z_mc::Vector{SMCg{np,T}} = Rnd_Out_Z_All(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_all_eps)
  elseif mc_opts.z_rnd == true
    z_mc = Rnd_Out_Z_Intv(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_eps)
  else
    z_mc = mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc
  end

  p_mc::Vector{SMCg{np,T}} = [SMCg(pmid[i],pmid[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  pref_mc::Vector{SMCg{np,T}} = copy(p_mc)
  aff_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(xA_mc[i].cc,xa_mc[i].cv,szero,szero,@interval(xa_mc[i].cv,xA_mc[i].cc),false,[∅],[1.0]) for i=1:nx]
  H_mc::Vector{SMCg{np,T}} = zeros(SMCg{np,T},nx)
  dH_mc::SparseMatrixCSC{SMCg{np,T},Int64} = spzeros(SMCg{np,T},nx,nx)
  Y::SparseMatrixCSC{T,Int64} = spzeros(T,nx,nx)

  # initalizes storage object for in-place sparse calculations
  SSto = SparseInSto()
  SSto.Xh = copy(H_mc)
  SSto.Yh = copy(H_mc)
  SSto.Zh = copy(H_mc)
  SSto.Xj = spzeros(SMCg{np,T},nx,nx)
  SSto.Yj = spzeros(SMCg{np,T},nx,nx)
  SSto.Zj = spzeros(SMCg{np,T},nx,nx)
  SSto.nx = copy(nx)

  # stores some things for repeated us
  optc = Any[szero,sone]

  # Begins loop to generate parameters
  sto_out[1] = copy(x_mc)
  for k=1:mc_opts.kmax
    aff_mc = [SMCg{np,Float64}(x_mc[i].cc,x_mc[i].cv,szero,szero,
                   @interval(x_mc[i].cv,x_mc[i].cc),false,[∅],[1.0]) for i=1:nx]
    if mc_opts.aff_rnd_all == true
      aff_mc = Rnd_Out_Z_All(aff_mc,mc_opts.aff_rnd_all_eps)
    elseif mc_opts.aff_rnd == true
      aff_mc = Rnd_Out_Z_Intv(aff_mc,mc_opts.aff_rnd_eps)
    end
    h!(H_mc,z_mc,p_mc)
    hj!(dH_mc,aff_mc,p_mc)
    Sparse_Precondition!(H_mc,dH_mc,mid.(Intv.(dH_mc)),SSto)
    dH_mc = transpose(dH_mc)

    # applies parametric iteration
    if (mc_opts.style == "NewtonGS")
      MCn_NewtonGS!(z_mc,x_mc,dH_mc,H_mc,nx)
    elseif (mc_opts.style == "KrawczykCW")
      MCn_KrawczykCW!(z_mc,x_mc,dH_mc,H_mc,nx)
    else
      error("Unsupported Style of Implicit Relaxation")
    end
    dH_mc = transpose(dH_mc)

    # update affine relaxations & correct
    Affine_Exp!(x_mc,p_mc,p_mc,xa_mc,xA_mc,z_mc,exp_opt)
    Correct_Exp!(z_mc,x_mc,X,nx,np,mc_opts.aff_correct_eps)
    if mc_opts.z_rnd_all == true
      z_mc = Rnd_Out_Z_All(z_mc,mc_opts.z_rnd_all_eps)
    elseif mc_opts.z_rnd == true
      z_mc = Rnd_Out_Z_Intv(z_mc,mc_opts.z_rnd_eps)
    end
    # store relaxation
    sto_out[k+1] = copy(x_mc)
    end
  return sto_out
end

"""
    InGenExpansionParams(h::Function, hj::Function, X::Vector{Interval{T}},
                       P::Vector{Interval{T}},pmid::Vector{T},mc_opts::mc_opts)

Generates relaxation of state variable at pmid in P and outputs results to use
as the parameters in subsequent calculations. Use a sparse LDU factorization with
full pivoting to precondition and performs function evaluations in-place. Inputs are:
* `h::Function`: h!(H,z,p) which implicit defines state function
* `hj::Function`: hj!(H,z,p) w.r.t z which implicit defines state function
* `P::Vector{Interval}`: State variable bounds
* `X::Vector{Interval{T}}`: Decision variable bounds
* `pmid::Vector{T}`: Point at which to generate relaxation
* `mc_opts`: Options for generating implicit function relaxation
Outs the following the tuple `sto_out`:
* `sto_out::Vector{SMCg{np,T}}` - McCormick parameter object stack
--------------------------------------------------------------------------------
"""
function IndGenExpansionParams(h!::Function, hj!::Function,
                            X::Vector{Interval{T}},
                            P::Vector{Interval{T}},
                            pmid::Vector{T},mc_opts::mc_opts{T}) where {T<:AbstractFloat}

  nx::Int64 = length(X)
  np::Int64 = length(P)
  szero::SVector{np,T} = @SVector zeros(np)
  sone::SVector{np,T} = @SVector ones(np)
  exp_opt::Vector{Any} = Any[nx,np,mc_opts.lambda]
  sto_out::Vector{Vector{SMCg{np,T}}} = [[zero(SMCg{np,T}) for i=1:nx] for j=1:(mc_opts.kmax+1)]

  x_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].lo,szero,szero,@interval(X[i].lo,X[i].hi),false,[∅],[1.0]) for i=1:nx]
  xa_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].lo,X[i].lo,szero,szero,@interval(X[i].lo,X[i].lo),false,[∅],[1.0]) for i=1:nx]
  xA_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(X[i].hi,X[i].hi,szero,szero,@interval(X[i].hi,X[i].hi),false,[∅],[1.0]) for i=1:nx]

  if mc_opts.z_rnd_all == true
    z_mc::Vector{SMCg{np,T}} = Rnd_Out_Z_All(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_all_eps)
  elseif mc_opts.z_rnd == true
    z_mc = Rnd_Out_Z_Intv(mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc,mc_opts.z_rnd_eps)
  else
    z_mc = mc_opts.lambda*xa_mc+(1.0-mc_opts.lambda)*xA_mc
  end

  p_mc::Vector{SMCg{np,T}} = [SMCg(pmid[i],pmid[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]
  pref_mc::Vector{SMCg{np,T}} = copy(p_mc)
  aff_mc::Vector{SMCg{np,T}} = [SMCg{np,T}(xA_mc[i].cc,xa_mc[i].cv,szero,szero,@interval(xa_mc[i].cv,xA_mc[i].cc),false,[∅],[1.0]) for i=1:nx]
  H_mc::Vector{SMCg{np,T}} = zeros(SMCg{np,T},nx)
  dH_mc::Array{SMCg{np,T},2} = zeros(SMCg{np,T},nx,nx)
  Y::Array{T,2} = zeros(T,nx,nx)

  # stores some things for repeated us
  optc = Any[szero,sone]

  # Begins loop to generate parameters
  sto_out[1] = copy(x_mc)
  for k=1:mc_opts.kmax
    aff_mc = [SMCg{np,Float64}(x_mc[i].cc,x_mc[i].cv,szero,szero,
                   @interval(x_mc[i].cv,x_mc[i].cc),false,[∅],[1.0]) for i=1:nx]
    if mc_opts.aff_rnd_all == true
      aff_mc = Rnd_Out_Z_All(aff_mc,mc_opts.aff_rnd_all_eps)
    elseif mc_opts.aff_rnd == true
      aff_mc = Rnd_Out_Z_Intv(aff_mc,mc_opts.aff_rnd_eps)
    end

    h!(H_mc,z_mc,p_mc)
    hj!(dH_mc,aff_mc,p_mc)
    Dense_Precondition!(H_mc,dH_mc,mid.(Intv.(dH_mc)),nx)

    # applies parametric iteration
    if (mc_opts.style == "NewtonGS")
      MC_NewtonGS!(z_mc,x_mc,dH_mc,H_mc,nx)
    elseif (mc_opts.style == "KrawczykCW")
      MC_KrawczykCW!(z_mc,x_mc,dH_mc,H_mc,nx)
    else
      error("Unsupported Style of Implicit Relaxation")
    end

    # update affine relaxations & correct
    Affine_Exp!(x_mc,p_mc,p_mc,xa_mc,xA_mc,z_mc,exp_opt)
    Correct_Exp!(z_mc,x_mc,X,nx,np,mc_opts.aff_correct_eps)
    if mc_opts.z_rnd_all == true
      z_mc = Rnd_Out_Z_All(z_mc,mc_opts.z_rnd_all_eps)
    elseif mc_opts.z_rnd == true
      z_mc = Rnd_Out_Z_Intv(z_mc,mc_opts.z_rnd_eps)
    end
    # store relaxation
    sto_out[k+1] = copy(x_mc)
    end
  #return sto_out,sto_z,sto_x
  return sto_out
end


text:
"""
    Fgamma(m, x, small)

Returns the value of Boys integral
```math
F_m(x) = \\int_0^1 t^{2m} e^{-xt^2} dt =
\\frac{1}{2 x^{m + \\frac{1}{2}}} \\gamma \\left( m + \\frac{1}{2}, x \\right),
```
where ``\\gamma(n, x)`` is the lower incomplete gamma function.

Taylor series for Boys integral, of use for small values of ``x``
```math
F_m(x) = \\sum_{k=0}^{\\infty} \\frac{(-x)^k}{k! (2(m + k) + 1)}
```

External links: [DLMF](https://dlmf.nist.gov/8.2.4),
[Wikipedia](https://en.wikipedia.org/wiki/Incomplete_gamma_function)

See also [`gamma_inc(a,x,ind=0)`](@ref SpecialFunctions.gamma_inc)
"""
function Fgamma(m::Int64, x::Float64, SMALL::Float64 = 1e-12)
    if x < SMALL
        # use small x expansion
        bf = 1/(2*m + 1) - x/(2*m + 3)
    else
        bf = gamma(0.5 + m) * gamma_inc(0.5 + m, x)[1] / (2 * x^(0.5 + m))
    end
    return bf
end

function test_fgamma()
    @testset "test_fgamma" begin
        for (x, res) in [
            (0.0, 1),
            (30.0, 0.161802159),
            (60.0, 0.114411404),
            (90.0, 0.0934165203),
            (120.0, 0.08090108),
            (300.0, 0.051166336),
        ]
            @test res ≈ Fgamma(0, x)
        end
    end # testset
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Miscellaneous functions related to the Text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _compute_table_size_data(header_str::Matrix{String}, data_str::Matrix{Vector{String}}, id_cols::Vector{Int}, Δc::Int, # Configurations columns_width::Vector{Int}, crop_subheader::Bool, maximum_columns_width::Vector{Int}, minimum_columns_width::Vector{Int}, noheader::Bool)

Compute the following table size data:

* The width of each column; and
* The maximum number of lines in the cells of each column.
"""
function _compute_table_size_data(
    header_str::Matrix{String},
    data_str::Matrix{Vector{String}},
    id_cols::Vector{Int},
    Δc::Int,
    # Configurations
    columns_width::Vector{Int},
    crop_subheader::Bool,
    maximum_columns_width::Vector{Int},
    minimum_columns_width::Vector{Int},
    noheader::Bool
)
    num_printed_rows, num_printed_cols = size(data_str)
    num_header_printed_rows, ~         = size(header_str)

    # The width of the columns in the table. Notice that the minimum allowed
    # column size must be one.
    cols_width = ones(Int, num_printed_cols)

    # The number of lines in each row.
    num_lines_in_row = ones(Int, num_printed_rows)

    # Regex to remove ANSI escape sequences so that we can compute the printable
    # size of the cell. This is required for Markdown cells.
    r_ansi_escape = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"

    @inbounds for j = 1:num_printed_cols
        largest_cell_width = 0

        # Header
        # ----------------------------------------------------------------------

        if !noheader
            for i = 1:num_header_printed_rows
                !isassigned(header_str, i, j) && continue

                # Remove the escape sequences to obtain only printable
                # characters.
                str_p = replace(header_str[i,j], r_ansi_escape => "")

                header_ij_len   = textwidth(str_p)

                # If the user wants to crop the subheader, then it should not be
                # used to compute the largest cell width of this column.
                if (i == 1) || (!crop_subheader)
                    largest_cell_width = max(largest_cell_width, header_ij_len)
                end
            end
        end

        # Data
        # ----------------------------------------------------------------------

        for i = 1:num_printed_rows
            !isassigned(data_str, i, j) && continue

            num_lines = length(data_str[i,j])

            num_lines_in_row[i] = max(num_lines_in_row[i], num_lines)

            for k = 1:num_lines
                # Remove the escape sequences to obtain only printable
                # characters.
                str_p = replace(data_str[i,j][k], r_ansi_escape => "")
                data_ijk_len = textwidth(str_p)
                largest_cell_width = max(largest_cell_width, data_ijk_len)
            end
        end

        # If we are in data columns, then we need to check the user preferences
        # before updating the table widths. Otherwise, we are in the added
        # columns (row number or row name) and we always need to update the
        # widths.
        if j ≤ Δc
            cols_width[j] = max(cols_width[j], largest_cell_width)
        else
            jc = id_cols[j-Δc]

            # Check if we need to replace the column width.
            if columns_width[jc] ≤ 0
                cols_width[j] = max(cols_width[j], largest_cell_width)

                # Make sure that the maximum column width is respected.
                if maximum_columns_width[jc] > 0 &&
                    maximum_columns_width[jc] < cols_width[j]
                    cols_width[j] = maximum_columns_width[jc]
                end

                # Make sure that the minimum column width is respected.
                if minimum_columns_width[jc] > 0 &&
                    minimum_columns_width[jc] > cols_width[j]
                    cols_width[j] = minimum_columns_width[jc]
                end
            else
                cols_width[j] = columns_width[jc]
            end
        end
    end

    return cols_width, num_lines_in_row
end


text:
using Test

@test foo(1) == "hello"


text:
# Command: Minimize Window
"""
	minimize!(session::Session; window::AbstractString = "current")::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}

The Minimize Window command invokes the window manager-specific "minimize" operation, if any, on the window containing the current top-level browsing context.
This typically hides the window in the system tray.

!!! note

    This command is implemented using the JSON Wire Protocol.
"""
function minimize!(
    session::Session;
    window::AbstractString = "current",
)::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}
    window = "current"
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/window/$window/minimize",
        [("Content-Type" => "application/json")],
    )
    @assert response.status == 200
    output = JSON3.read(response.body).value
    (
        width = output.width,
        height = output.height,
        x = output.x,
        y = output.y,
    )::NamedTuple{(:width, :height, :x, :y),NTuple{4,Int64}}
end


text:
# Model hyperparameters

"""
RAPIDS Docs for Logistic Regression:
    https://docs.rapids.ai/api/cuml/stable/api.html#logistic-regression

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = [repeat([0], 50)..., repeat([1], 50)...]

model = LogisticRegression()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct LogisticRegression <: MMI.Probabilistic
    handle = nothing
    penalty::String = "l2"::(_ in ("none", "l1", "l2", "elasticnet"))
    tol::Float64 = 1e-4::(_ > 0)
    C::Float64 = 1.0::(_ > 0)
    fit_intercept::Bool = true
    class_weight = nothing
    max_iter::Int = 1000::(_ > 0)
    linesearch_max_iter::Int = 50::(_ > 0)
    verbose::Bool = false
    l1_ratio = nothing
    solver::String = "qn"::(_ in ("qn", "lbfgs", "owl"))
end

"""
RAPIDS Docs for the MBSGD Classifier:
    https://docs.rapids.ai/api/cuml/stable/api.html#mini-batch-sgd-classifier

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = [repeat([0], 50)..., repeat([1], 50)...]

model = MBSGDClassifier()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct MBSGDClassifier <: MMI.Probabilistic
    handle = nothing
    loss::String = "squared_loss"::(_ in ("hinge", "log", "squared_loss"))
    penalty::String = "none"::(_ in ("none", "l1", "l2", "elasticnet"))
    alpha::Float64 = 0.0001::(_ >= 0)
    l1_ratio::Float64 = 0.15::(_ >= 0)
    batch_size::Int = 32::(_ > 0)
    fit_intercept::Bool = true
    epochs::Int = 1000::(_ > 0)
    tol::Float64 = 1e-3::(_ > 0)
    shuffle::Bool = true
    eta0::Float64 = 0.001::(_ > 0)
    power_t::Float64 = 0.5::(_ > 0)
    learning_rate::String = "constant"::(_ in ("adaptive", "constant", "invscaling", "optimal"))
    n_iter_no_change::Int = 5::(_ > 0)
    verbose::Bool = false
end

"""
RAPIDS Docs for the Nearest Neighbors Classifier:
    https://docs.rapids.ai/api/cuml/stable/api.html#nearest-neighbors-classification

Example:
```
using RAPIDS
using MLJ

X = rand(100, 5)
y = [repeat([0], 50)..., repeat([1], 50)...]

model = KNeighborsClassifier()
mach = machine(model, X, y)
fit!(mach)
preds = predict(mach, X)
```
"""
MLJModelInterface.@mlj_model mutable struct KNeighborsClassifier <: MMI.Probabilistic
    handle = nothing
    algorithm::String = "brute"::(_ in ("brute",))
    metric::String = "euclidean"
    weights::String = "uniform"::(_ in ("uniform",))
    verbose::Bool = false
end

# Multiple dispatch for initializing models
model_init(mlj_model::LogisticRegression) = cuml.LogisticRegression(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::MBSGDClassifier) = cuml.MBSGDClassifier(; mlj_to_kwargs(mlj_model)...)
model_init(mlj_model::KNeighborsClassifier) = cuml.KNeighborsClassifier(; mlj_to_kwargs(mlj_model)...)

const CUML_CLASSIFICATION = Union{LogisticRegression, MBSGDClassifier, KNeighborsClassifier}

# add metadata
MMI.load_path(::Type{<:LogisticRegression}) = "$PKG.LogisticRegression"
MMI.load_path(::Type{<:MBSGDClassifier}) = "$PKG.MBSGDClassifier"
MMI.load_path(::Type{<:KNeighborsClassifier}) = "$PKG.KNeighborsClassifier"

MMI.input_scitype(::Type{<:CUML_CLASSIFICATION}) = Union{AbstractMatrix, Table(Continuous)}
MMI.target_scitype(::Type{<:CUML_CLASSIFICATION}) = Union{AbstractVector, Table(Continuous)} #should we use count even though we input floats?

MMI.docstring(::Type{<:LogisticRegression}) = "cuML's LogisticRegression: https://docs.rapids.ai/api/cuml/stable/api.html#logistic-regression"
MMI.docstring(::Type{<:MBSGDClassifier}) = "cuML's MBSGDClassifier: https://docs.rapids.ai/api/cuml/stable/api.html#mini-batch-sgd-classifier"
MMI.docstring(::Type{<:KNeighborsClassifier}) = "cuML's KNeighborsClassifier: https://docs.rapids.ai/api/cuml/stable/api.html#nearest-neighbors-classification"

# fit methods
function MMI.fit(mlj_model::CUML_CLASSIFICATION, verbosity, X, y, w=nothing)
    # fit the model
    model = model_init(mlj_model)
    model.fit(prepare_input(X), prepare_input(y))
    fitresult = model

    # save result
    cache = nothing
    report = ()
    return (fitresult, cache, report)
end

# predict methods
function MMI.predict(mlj_model::CUML_CLASSIFICATION, fitresult, Xnew)
    model  = fitresult
    py_preds = model.predict(prepare_input(Xnew))
    preds = pyconvert(Array, py_preds)

    return preds
end

# Classification metadata
MMI.metadata_pkg.((LogisticRegression, MBSGDClassifier, KNeighborsClassifier),
    name = "cuML Classification Methods",
    uuid = "2764e59e-7dd7-4b2d-a28d-ce06411bac13", # see your Project.toml
    url  = "https://github.com/tylerjthomas9/RAPIDS.jl",  # URL to your package repo
    julia = false,          # is it written entirely in Julia?
    license = "MIT",        # your package license
    is_wrapper = true,      # does it wrap around some other package?
)


text:
using PlotlyJS
using Statistics
using SurrogateModelOptim
using BlackBoxOptim
using Parameters
using Printf

options=SurrogateModelOptim.Options(
    iterations=10,
    smooth=false, # Should only be false for smooth deterministic functions (noise free)
    num_interpolants=10, #Preferably even number of added processes
    num_start_samples=4,
    constrained_seed_gens=0,
    rbf_opt_gens=1000,
    create_final_surrogate=true, #Use the results from last iteration to
                                 #re-create the surrogate before using it for plotting
        )

fitness_schaffer1(x) = (sum(abs2, x), sum(abs2, x .- 2.0))

search_range = [(-10.0,10.0) for _ in 1:3]

function multi_objective_smoptimize(f,search_range,options)

    #2 objectives
    n_funcs = 2

    #Load some option values
    @unpack num_start_samples, sampling_plan_opt_gens,
            iterations, trace = options

    lhc_plan = SurrogateModelOptim.scaled_LHC_sampling_plan(search_range,num_start_samples,10000)

    #Evaluate sampling plan
    lhc_samples = [SurrogateModelOptim.f_opt_eval(x->f(x)[i],lhc_plan) for i in 1:n_funcs]

    #Initialize variables to be returned
    sm_interpolant = [nothing for _ in 1:n_funcs]
    infill_type = [Array{Symbol,1}(undef,0) for _ in 1:n_funcs]
    infill_prediction = [Array{Float64,1}(undef,0) for _ in 1:n_funcs]
    optres = [nothing for _ in 1:n_funcs]
    infill_plan = Array{Float64,2}(undef,size(lhc_plan,1),0)
    infill_sample = [Array{Float64,2}(undef,1,0) for _ in 1:n_funcs]
    sm_med = nothing

    #Run the entire optimization iterations number of times
    for i = 1:iterations
        SurrogateModelOptim.print_iteration(trace,i,iterations)

        plan_all = [lhc_plan infill_plan]
        samples_all = [[lhc_samples[j] infill_sample[j]] for j in 1:n_funcs]

        #Create the optimized Radial Basis Function interpolant for each objective function
        sm_fs = Tuple(SurrogateModelOptim.surrogate_model(plan_all, samples_all[j]; options=options)[1] for j in 1:n_funcs)
        sm = x -> Tuple(sm_fs[i](x) for i in 1:n_funcs)

        #We add two points per iteration, one exploring and one exploiting sample
        sm_med = x-> median.(sm(x))
        sm_std = x-> -1 .*std.(sm(x))

        res_med = bboptimize(sm_med; Method=:borg_moea,
                FitnessScheme=ParetoFitnessScheme{n_funcs}(is_minimizing=true),
                SearchRange=search_range, NumDimensions=size(plan_all,2), ϵ=0.05,
                MaxSteps=10_000, TraceInterval=1.0, TraceMode=:silent);

        res_std = bboptimize(sm_std; Method=:borg_moea,
                FitnessScheme=ParetoFitnessScheme{n_funcs}(is_minimizing=true),
                SearchRange=search_range, NumDimensions=size(plan_all,2), ϵ=0.05,
                MaxSteps=10_000, TraceInterval=1.0, TraceMode=:silent);

        #Add two points randomly selected points from the pareto front
        pf_res_med = pareto_frontier(res_med); p_med = pf_res_med[rand(1:length(pf_res_med))]
        pf_res_std = pareto_frontier(res_std); p_std = pf_res_std[rand(1:length(pf_res_std))]

        infill_plan_new = [permutedims(p_med.inner.params') permutedims(p_std.inner.params')]

        #Evaluate the new infill points
        infill_sample_new = [SurrogateModelOptim.f_opt_eval(x->f(x)[j],infill_plan_new,samples_all[j];trace=:silent) for j in 1:n_funcs]

        println("\nEvaluating functions------------")
        [str = printstyled(@sprintf("Objective value: %-15.7g \n",minimum(infill_sample_new[i]))) for i in 1:n_funcs]
        println("--------------------------------")

        #Add infill points
        infill_plan = [infill_plan infill_plan_new]
        infill_sample = [[infill_sample[j] infill_sample_new[j]] for j in 1:n_funcs]
        samples_all = [[lhc_samples[j] infill_sample[j]] for j in 1:n_funcs]
    end

    return sm_med
end
sm = multi_objective_smoptimize(fitness_schaffer1,search_range,options)

function plot_pareto(fun,sm)

    res_bbopt = bboptimize(fun; Method=:borg_moea,
                FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true),
                SearchRange=search_range, NumDimensions=3, ϵ=0.05,
                MaxSteps=50000, TraceInterval=1.0, TraceMode=:silent);
    pf1= pareto_frontier(res_bbopt)
    bb_f1 = [pf1[i].fitness[1] for i in 1:length(pf1)]
    bb_f2 = [pf1[i].fitness[2] for i in 1:length(pf1)]

    res_smopt = bboptimize(sm; Method=:borg_moea,
                FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true),
                SearchRange=search_range, NumDimensions=3, ϵ=0.05,
                MaxSteps=50000, TraceInterval=1.0, TraceMode=:silent);
    pf2 = pareto_frontier(res_smopt)
    sm_f1 = [pf2[i].fitness[1] for i in 1:length(pf2)]
    sm_f2 = [pf2[i].fitness[2] for i in 1:length(pf2)]

    trace1 = scatter(;x=bb_f1, y=bb_f2, mode="markers",name= "BB_optim")
    trace2 = scatter(;x=sm_f1, y=sm_f2, mode="markers",name= "SM_optim")
    layout = Layout(title="Pareto")
    plot([trace1, trace2],layout)
end

# Plot the results
display(plot_pareto(fitness_schaffer1,sm))


text:
using Pkg

Pkg.activate(dirname(@__FILE__))
Pkg.instantiate()

using HTTP
using JSON
using Blink
using WebIO

function draw_queens(url; N=8, W=30)
    try
        resp = HTTP.get(url*"/webqueens/$N")
        json = JSON.parse(String(resp.body))
        allqueens = json["data"]
        w = Window()
        body!(w, dom"div"(
            dom"canvas.#ChessBoard[width=$(N*W), height=$(N*W)]"(),
            dom"script"(
                """
                    function draw_queen(ctx, sx, sy){
                        ctx.fillStyle = "#C0C0C0";
                        ctx.fillText("Q", sx+$W/5, sy+4*$W/5);
                        ctx.fillStyle = "#000000";
                    }

                    function draw_board(ctx, queens){
                        n = queens.length
                        ctx.fillStyle = "#FFFFFF";
                        ctx.fillRect(0, 0, n*$W, n*$W);
                        ctx.fillStyle = "#000000";
                        sx = 0; sy = 0;

                        for (i=n-1; i >= 0; i--){
                            for (j=0; j < n; j++){
                                if ((i+j)%2 == 0)
                                    ctx.fillRect(sx, sy, $W, $W);
                                if (j == queens[i] - 1)
                                    draw_queen(ctx, sx, sy)
                                sx += $W;
                            }
                            sy += $W; sx = 0;
                        }
                    }

                    var c = document.getElementById("ChessBoard");
                    var ctx = c.getContext("2d");
                    ctx.font = "$(4*div(W,5))px sans-serif";
                    draw_board(ctx, $([allqueens...]));
                """
            )
        ))
    catch e
        return "Error occurred : $e"
    end
end

max_stars_repo_path: Chapter 14/client/client.jl

max_issues_repo_path: Chapter 14/client/client.jl

max_forks_repo_path: Chapter 14/client/client.jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DeploymentWhatIf <: SwaggerModel
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    properties::Any # spec type: Union{ Nothing, DeploymentWhatIfProperties } # spec name: properties

    function DeploymentWhatIf(;location=nothing, properties=nothing)
        o = new()
        validate_property(DeploymentWhatIf, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(DeploymentWhatIf, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type DeploymentWhatIf

const _property_map_DeploymentWhatIf = Dict{Symbol,Symbol}(Symbol("location")=>Symbol("location"), Symbol("properties")=>Symbol("properties"))
const _property_types_DeploymentWhatIf = Dict{Symbol,String}(Symbol("location")=>"String", Symbol("properties")=>"DeploymentWhatIfProperties")
Base.propertynames(::Type{ DeploymentWhatIf }) = collect(keys(_property_map_DeploymentWhatIf))
Swagger.property_type(::Type{ DeploymentWhatIf }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DeploymentWhatIf[name]))}
Swagger.field_name(::Type{ DeploymentWhatIf }, property_name::Symbol) =  _property_map_DeploymentWhatIf[property_name]

function check_required(o::DeploymentWhatIf)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ DeploymentWhatIf }, name::Symbol, val)
end


text:
using PyCall

@pyimport tensorflow as tf
@pyimport tensorflow.models.rnn.rnn_cell as tfrnncell
@pyimport tensorflow.python.platform.flags as tf_flags
import TensorFlowBuilder: TFParser

target_path = ARGS[1]
const outdir = joinpath(target_path, "src", "API")

println("Generating source code to $outdir")

if !ispath(outdir)
  println("mkdir $outdir")
  mkdir(outdir)
end

@pyimport tensorflow as tf
for (pyimp, pymod, pyname, jlname) in [
                                       (:(tensorflow), :tf, :tf, :Tf),
                                       (:(tensorflow.python.ops.nn), :tf_nn, :(tf.nn), :TfNn),
                                       (:(tensorflow.python.training.training), :tf_train, :(tf.train), :TfTrain),
                                       (:(tensorflow.models.rnn.rnn_cell), :tf_rnn_cell, :(tfrnncell), :TfRnnCell),
                                       (:(tensorflow.python.platform.flags), :tf_flags, :(tf_flags), :TfFlags),
                                       ]
  @eval m = $pyname
  fname = joinpath(outdir, "$jlname.jl")
  if ispath(fname)
    println("Deleting $fname")
    rm(fname)
  end
    f = open(fname, "w")
    version = "?"
    try
        version = tf.__version__
    catch
    end
  write(f, """
"Generated automatically by TensorFlowBuilder, from TensorFlow Python version $(version)"
#"TensorFlow, the TensorFlow logo and any related marks are trademarks of Google Inc.""
module $jlname
using PyCall
@pyimport tensorflow as tf
@pyimport $pyimp as $pymod
import TensorFlow.CoreTypes: *
using TensorFlow.CoreTypes
""")

  TFParser.tfwritejulia(m, pymod, f)
  write(f, "end\n")
  close(f)
end

println("Loading generated module")

include(joinpath("..", target_path, "src", "TensorFlow.jl"))
using .TensorFlow


text:
using LinearAlgebra, Plots
using KitBase, ProgressMeter
import KitML

cd(@__DIR__)
include("math.jl")

model = KitML.load_model("best_model.h5"; mode = :tf)

begin
    # space
    x0 = -3.5
    x1 = 3.5
    y0 = -3.5
    y1 = 3.5
    nx = 50#100
    ny = 50#100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 5.0
    cfl = 0.7
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # moments
    L = 1
    #ne = (L + 1)^2
    ne = GetBasisSize(L,3)

    phi = zeros(Float32, ne, nx, ny)
    α = zeros(Float32, ne, nx, ny)
    h =
    #m = KitBase.eval_spherharmonic(points, L)
    m = ComputeSphericalBasisAnalytical(points)
end

function is_absorb(x::T, y::T) where {T<:Real}
    cds = Array{Bool}(undef, 11) # conditions

    cds[1] = -2.5<x<-1.5 && 1.5<y<2.5
    cds[2] = -2.5<x<-1.5 && -0.5<y<0.5
    cds[3] = -2.5<x<-1.5 && -2.5<y<-1.5
    cds[4] = -1.5<x<-0.5 && 0.5<y<1.5
    cds[5] = -1.5<x<-0.5 && -1.5<y<-0.5
    cds[6] = -0.5<x<0.5 && -2.5<y<-1.5
    cds[7] = 0.5<x<1.5 && 0.5<y<1.5
    cds[8] = 0.5<x<1.5 && -1.5<y<-0.5
    cds[9] = 1.5<x<2.5 && 1.5<y<2.5
    cds[10] = 1.5<x<2.5 && -0.5<y<0.5
    cds[11] = 1.5<x<2.5 && -2.5<y<-1.5

    if any(cds) == true
        return true
    else
        return false
    end
end

begin
    σs = zeros(Float64, nx, ny)
    σa = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if is_absorb(pspace.x[i, j], pspace.y[i, j])
            σs[i, j] = 0.0
            σa[i, j] = 10.0
        else
            σs[i, j] = 1.0
            σa[i, j] = 0.0
        end
    end
    σt = σs + σa
    σq = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if -0.5<pspace.x[i, j]<0.5 && -0.5<pspace.y[i, j]<0.5
            σq[i, j] = 1.0 / (4.0 * π)
        else
            σq[i, j] = 0.0
        end
    end
end

#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], σs')

#Training loop

global t = 0.0
flux1 = zeros(ne, nx + 1, ny)
flux2 = zeros(ne, nx, ny + 1)

αT = zeros(nx*ny, ne)
phiT = zeros(nx*ny, ne)

    # Start Time loop

#begin
#for metaIter = 1:1
metaIter = 1
#initial conditions

phi = zeros(Float32, ne, nx, ny)
phi_old = zeros(Float32, ne, nx, ny)
for j = 1:nx
    for i = 1:ny
        phi[1, i, j] = 1e-6
    end
end

anim = @animate for iter = 1:1000
    println("Iteration $(iter)")

    #restructuring of the matrix
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :, :][:]
    end
    #Call the neural network
    αT = KitML.neural_closure(model, phiT)
    @inbounds for i = 1:ne
        α[i, :, :] .= reshape(αT[:, i], nx, ny)
    end

    #### Introduce helper variables
    phi_old .= phi
    phi_train = Vector()
    h_train = Vector()
    errorsFound = false
    #####

    # realizability reconstruction
    error = 0
    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            phi[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)
            # print( phi[:, i, j])
            # print("|")
            # println( phi_old[:, i, j] )

            # Check the prediction error
            error = norm( phi[:, i, j]-phi_old[:, i, j],2)^2 # take care of the criterion
            if error > 0.0001
                #println(error)
                errorsFound = true
                #add point to retraining batch. compute alpha and h w.r.t phi_old
                res = KitBase.optimize_closure(α[:, i, j], m, weights, phi_old[:, i, j], KitBase.maxwell_boltzmann_dual)
                α[:, i, j] .= res.minimizer
                phi[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)

                #Setup trianing data
                h = computeEntropyH(α[:, i, j], m, weights,maxwell_boltzmann_primal, KitBase.maxwell_boltzmann_dual_prime)

                h_train =vcat(h_train, h)
                phi_train =vcat(phi_train, phi[:, i, j])

            end

        end
    end

    #retrain network on cells with too much error, if there is any
    if errorsFound
        print("Current number of  errorous datapoints: ")
        print(size(h_train)[1])
        print(" of ")
        print(nx*ny)
        print(" in iter ")
        println(iter)
        # Transform into matrix format
        phiTrainMat = zeros(size(h_train)[1],4)
        hTrainMat = zeros(size(h_train)[1],1)
        for i in 1:size(h_train)[1]
            phiTrainMat[i,1] = phi_train[4*(i-1)+1]
            phiTrainMat[i,2] = phi_train[4*(i-1)+2]
            phiTrainMat[i,3] = phi_train[4*(i-1)+3]
            phiTrainMat[i,4] = phi_train[4*(i-1)+4]
            hTrainMat[i,1]= h_train[i]
        end
        model.fit(phiTrainMat,hTrainMat,epochs =50,batch_size = 32,validation_split=0.0,verbose= 0)
        println("model retrained")
        model.save("best_model.h5")
    end

    # flux
    fη1 = zeros(nq)
    @inbounds for j = 1:ny
        for i = 2:nx
            KitBase.flux_kfvs!(fη1, KitBase.maxwell_boltzmann_dual.(α[:, i-1, j]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 1], dt)

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    @inbounds for i = 1:nx
        for j = 2:ny
            KitBase.flux_kfvs!(fη2, KitBase.maxwell_boltzmann_dual.(α[:, i, j-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 2], dt)

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (weights .* fη2))
            end
        end
    end

    # update
    @inbounds for j = 2:ny-1
        for i = 2:nx-1
            for q = 1:1
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (σs[i, j] * phi[q, i, j] - σt[i, j] * phi[q, i, j]) * dt +
                    σq[i, j] * dt * 100.0
            end

            for q = 2:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (-σt[i, j] * phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
    contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])
end
cd(@__DIR__)
gif(anim, "lattice_mn_neural_reinforcement.gif")

#global t=0
#end
#println(phi[1, :, :])
#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])


text:
export Dup

struct Dup{T} <: Bundle{T}
    x::T
    twin::T
end
function Dup(x::T) where T
   Dup{T}(x, copy(x))
end

@fieldview NiLangCore.value(x::Dup) = x.x
Base.isapprox(a::Dup, b::Dup; kwargs...) = isapprox(a.x, b.x; kwargs...) && isapprox(a.twin, b.twin; kwargs...)


text:
using BinaryBuilder

cuda_version = v"11.0.3"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/11.0_20200825/cudnn-11.0-linux-x64-v8.0.3.33.tgz",
                  "8924bcc4f833734bdd0009050d110ad0c8419d3796010cf7bc515df654f6065a")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/11.0_20200825/cudnn-11.0-linux-ppc64le-v8.0.3.33.tgz",
                  "c2d0519831137b43d0eebe07522edb4ef5d62320e65e5d5fa840a9856f25923d")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/11.0_20200825/cudnn-11.0-windows-x64-v8.0.3.33.zip",
                  "67f64109a677fb43ddbf730b1fbf8f9ba9bf0214d3a9e1027ac546255d4ae3c6")
]

include("../common.jl")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil