text:
""
function constraint_power_balance_shunt_dispatch(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = ref(pm, nw, :bus, i)
    bus_arcs = ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
    bus_arcs_sw = ref(pm, nw, :bus_arcs_sw, i)
    bus_gens = ref(pm, nw, :bus_gens, i)
    bus_loads = ref(pm, nw, :bus_loads, i)
    bus_storage = ref(pm, nw, :bus_storage, i)

    bus_shunts_const = ref(pm, :bus_shunts_const, i)
    bus_shunts_var = ref(pm, :bus_shunts_var, i)

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

    bus_gs_const = Dict(k => ref(pm, :shunt, k, "gs") for k in bus_shunts_const)
    bus_bs_const = Dict(k => ref(pm, :shunt, k, "bs") for k in bus_shunts_const)

    constraint_power_balance_shunt_dispatch(pm, nw, i, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
end

""
function constraint_power_balance_shunt_dispatch_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    bus = ref(pm, nw, :bus, i)
    bus_arcs = ref(pm, nw, :bus_arcs, i)
    bus_arcs_dc = ref(pm, nw, :bus_arcs_dc, i)
    bus_arcs_sw = ref(pm, nw, :bus_arcs_sw, i)
    bus_gens = ref(pm, nw, :bus_gens, i)
    bus_loads = ref(pm, nw, :bus_loads, i)
    bus_storage = ref(pm, nw, :bus_storage, i)

    bus_shunts_const = ref(pm, :bus_shunts_const, i)
    bus_shunts_var = ref(pm, :bus_shunts_var, i)

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

    bus_gs_const = Dict(k => ref(pm, :shunt, k, "gs") for k in bus_shunts_const)
    bus_bs_const = Dict(k => ref(pm, :shunt, k, "bs") for k in bus_shunts_const)

    constraint_power_balance_shunt_dispatch_soft(pm, nw, i, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
end

""
function constraint_ohms_yt_from_goc(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    branch = ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    g_fr = branch["g_fr"]
    b_fr = branch["b_fr"]
    tm = branch["tap"]

    if branch["transformer"]
        constraint_ohms_yt_from_goc(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    else
        PowerModels.constraint_ohms_yt_from(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    end
end

""
function constraint_thermal_limit_from_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    branch = ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)

    if haskey(branch, "rate_a")
        constraint_thermal_limit_from_soft(pm, nw, f_idx, branch["rate_a"])
    end
end

""
function constraint_thermal_limit_to_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    branch = ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    t_idx = (i, t_bus, f_bus)

    if haskey(branch, "rate_a")
        constraint_thermal_limit_to_soft(pm, nw, t_idx, branch["rate_a"])
    end
end

""
function constraint_gen_active_deviation(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    gen = ref(pm, nw, :gen, i)

    constraint_gen_active_deviation(pm, nw, i, gen["pg"])
end

""
function constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    cut = ref(pm, :branch_flow_cuts, i)
    branch = ref(pm, nw, :branch, cut.branch_id)

    if haskey(branch, "rate_c")
        constraint_branch_contingency_ptdf_thermal_limit_from_soft(pm, nw, i, cut.bus_injection, branch["rate_c"])
    end
end

""
function constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    cut = ref(pm, :branch_flow_cuts, i)
    branch = ref(pm, nw, :branch, cut.branch_id)

    if haskey(branch, "rate_c")
        constraint_branch_contingency_ptdf_thermal_limit_to_soft(pm, nw, i, cut.bus_injection, branch["rate_c"])
    end
end

# revisit after issue #14 is closed
# ""
# function constraint_gen_contingency_ptdf_thermal_limit_from_soft(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
#     cut = ref(pm, :gen_flow_cuts, i)
#     branch = ref(pm, nw, :branch, cut.branch_id)
#     gen = ref(pm, :gen, cut.gen_id)
#     gen_bus = ref(pm, :bus, gen["gen_bus"])
#     gen_set = Set(i for i in ref(pm, :area_gens)[gen_bus["area"]] if i != cut.gen_id)

#     gen_alpha = Dict(i => ref(pm, :gen, i)["alpha"] for i in gen_set)
#     gen_bus = Dict(i => ref(pm, :gen, i)["gen_bus"] for i in gen_set)

#     if haskey(branch, "rate_c")
#         constraint_gen_contingency_ptdf_thermal_limit_from_soft(pm, nw, i, cut.bus_injection, branch["rate_c"], gen_set, gen_alpha, gen_bus)
#     end
# end


text:
using Plots, Test
import RecipesPipeline

@testset "Contours" begin
    @testset "check_contour_levels" begin
        @test Plots.check_contour_levels(2) === nothing
        @test Plots.check_contour_levels(-1.0:0.2:10.0) === nothing
        @test Plots.check_contour_levels([-100, -2, -1, 0, 1, 2, 100]) === nothing
        @test_throws ArgumentError Plots.check_contour_levels(1.0)
        @test_throws ArgumentError Plots.check_contour_levels((1, 2, 3))
        @test_throws ArgumentError Plots.check_contour_levels(-3)
    end

    @testset "RecipesPipeline.preprocess_attributes!" begin
        function equal_after_pipeline(kw)
            kw′ = deepcopy(kw)
            RecipesPipeline.preprocess_attributes!(kw′)
            kw == kw′
        end

        @test equal_after_pipeline(KW(:levels => 1))
        @test equal_after_pipeline(KW(:levels => 1:10))
        @test equal_after_pipeline(KW(:levels => [1.0, 3.0, 5.0]))
        @test_throws ArgumentError RecipesPipeline.preprocess_attributes!(
            KW(:levels => 1.0),
        )
        @test_throws ArgumentError RecipesPipeline.preprocess_attributes!(
            KW(:levels => (1, 2, 3)),
        )
        @test_throws ArgumentError RecipesPipeline.preprocess_attributes!(KW(:levels => -3))
    end

    @testset "contour[f]" begin
        x = (-2π):0.1:(2π)
        y = (-π):0.1:π
        z = cos.(y) .* sin.(x')

        @testset "Incorrect input" begin
            @test_throws ArgumentError contour(x, y, z, levels = 1.0)
            @test_throws ArgumentError contour(x, y, z, levels = (1, 2, 3))
            @test_throws ArgumentError contour(x, y, z, levels = -3)
        end

        @testset "Default number" begin
            @test contour(x, y, z)[1][1].plotattributes[:levels] ==
                  Plots._series_defaults[:levels]
        end

        @testset "Number" begin
            @testset "$n contours" for n in (2, 5, 100)
                p = contour(x, y, z, levels = n)
                attr = p[1][1].plotattributes
                @test attr[:seriestype] == :contour
                @test attr[:levels] == n
            end
        end

        @testset "Range" begin
            levels = -1:0.5:1
            @test contour(x, y, z, levels = levels)[1][1].plotattributes[:levels] == levels
        end

        @testset "Set of levels" begin
            levels = [-1, 0.25, 0, 0.25, 1]
            @test contour(x, y, z, levels = levels)[1][1].plotattributes[:levels] == levels
        end
    end
end

ext: jl


text:
# helpers to make addition/multiplication of distance and regularizer possible

import Base.+
function +(a::Function, b::SparseMatrixCSC)
	f(x) = a(x) + b*x
	return f
end


text:
export QuCircuit

module QuCircuit

import QMTK: forward, backward
using Compat, StaticArrays, QMTK

end


text:
using Test
using DataFrames
using Random
using LinearAlgebra

using LinRegOutliers

@testset "Apply function to columns" begin
    ones_vector = ones(Float64, 10)
    zeros_vector = zeros(Float64, 10)
    mat = hcat(ones_vector, zeros_vector)

    @test applyColumns(sum, mat) == [10.0, 0.0]
    @test applyColumns(mean, mat) == [1.0, 0.0]
end

@testset "Ordinary Least Squares" begin
    tol = 0.0001
    #  The model is exatly y = 5 + 5x
    var1 = Float64[1, 2, 3, 4, 5]
    X = hcat(ones(5), var1)
    betas = [5.0, 5.0]
    y = X * betas

    # OLS
    olsreg = ols(X, y)
    @test isapprox(coef(olsreg), betas, atol=tol)
    @test isapprox(residuals(olsreg), zeros(Float64, 5), atol=tol)
    @test isapprox(predict(olsreg), y, atol=tol)
end

@testset "Weighted Least Squares" begin
    tol = 0.0001
    n = 7
    #  The model is exatly y = 5 + 5x
    var1 = Float64[1, 2, 3, 4, 5, 6, 7]
    X = hcat(ones(n), var1)
    betas = [5.0, 5.0]
    y = X * betas
    y[n - 1] = 5000.0
    y[n] = 5000.0
    wts = [1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0]

    # WLS
    olsreg = wls(X, y, wts)
    @test isapprox(coef(olsreg), betas, atol=tol)
    @test isapprox(residuals(olsreg)[1:(n - 2)], zeros(Float64, n - 2), atol=tol)
    @test isapprox(predict(olsreg)[1:(n - 2)], y[1:(n - 2)], atol=tol)
end

@testset "Basis - createRegressionSetting, designMatrix, responseVector" begin
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5],
        y=[2.0, 4, 6, 8, 10]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)

    @test designMatrix(setting) == hcat(ones(5), dataset[:, "x"])
    @test responseVector(setting) == dataset[:, "y"]
    @test setting.formula == @formula(y ~ x)
end

@testset "dffit - single case" begin
    # Since this regression setting is deterministic,
    # omission of a single observation has not an effect on
    # the predicted response.
    myeps = 10.0^(-6.0)
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5],
        y=[2.0, 4, 6, 8, 10]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    n, _ = size(dataset)
    for i in 1:n
        mydiff = dffit(setting, i)
        @test abs(mydiff) < myeps
    end
end

@testset "dffit - for all observations" begin
    # Since this regression setting is deterministic,
    # omission of a single observation has not an effect on
    # the predicted response.
    myeps = 10.0^(-6.0)
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5],
        y=[2.0, 4, 6, 8, 10]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    stats = dffit(setting)
    for element in stats
        @test abs(element) < myeps
    end
end

@testset "hat matrix" begin
    X = [1.0 1.0; 1.0 2.0; 1.0 3.0; 1.0 4.0; 1.0 5.0];
    hats_real = X * inv(X' * X) * X'
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5],
        y=[2.0, 4, 6, 8, 10]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    hats_calculated = hatmatrix(setting)

    @test hats_real == hats_calculated
end

@testset "studentized residuals" begin
    dataset = DataFrame(
        x=Float64[1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
        y=Float64[2.1, 3.9, 6.2, 7.8, 11.0, 18.0]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    resi = studentizedResiduals(setting)
    @test abs(resi[1]) < 1.5
    @test abs(resi[2]) < 1.5
    @test abs(resi[3]) < 1.5
    @test abs(resi[4]) < 1.5
    @test abs(resi[5]) < 1.5
    @test abs(resi[6]) > 1.5
end

@testset "adjusted residuals" begin
    eps = 0.0001
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5],
        y=[2.0, 4, 6, 8, 10]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    resi = adjustedResiduals(setting)
    for element in resi
        @test abs(element) < eps
    end
end

@testset "Cook's distance - Phone data" begin
    eps = 0.00001
    setting = createRegressionSetting(@formula(calls ~ year), phones)
    knowncooks = [0.005344774190779771, 0.0017088194691033181, 0.00016624914057961155,
                    3.16444525831206e-5, 0.0005395058666404081, 0.0014375008774859539,
                    0.0024828140956511258, 0.0036279720445167277, 0.004357605989540906,
                    0.005288503758364767, 0.006313578057565415, 0.0076561205696857254,
                    0.009568574875389256, 0.009970039008782357, 0.02610396373381051,
                    0.029272523880917646, 0.05091236198400663, 0.08176555044049343,
                    0.14380266904640235, 0.26721539425047447, 0.051205153558783356,
                    0.13401084683481085, 0.16860324592350226, 0.2172819114905912]
    cookdists = cooks(setting)
    @test map((x, y) -> abs(x - y) < eps, cookdists, knowncooks) == trues(24)
end

@testset "Jacknifed standard error of regression" begin
    eps = 0.00001
    dataset = DataFrame(
        x=[1.0, 2, 3, 4, 5000],
        y=[2.0, 4, 6, 8, 1000]
    )
    setting = createRegressionSetting(@formula(y ~ x), dataset)
    @test jacknifedS(setting, 1) != 0
    @test jacknifedS(setting, 2) != 0
    @test jacknifedS(setting, 3) != 0
    @test jacknifedS(setting, 4) != 0
    @test jacknifedS(setting, 5) < eps
end

@testset "Mahalanobis squared distances" begin
    myeps = 10.0^(-6.0)
    x = [1.0, 2.0, 3.0, 4.0, 5.0]
    y = reverse(x)
    y[1] = y[1] * 10.0
    x[1] = x[1] * 10.0
    datamat = DataFrame(x=x, y=y)
    dmat = mahalanobisSquaredMatrix(datamat)
    d = diag(dmat)
    @test abs(d[1] - 3.2) < myeps
    @test abs(d[2] - 2.0) < myeps
    @test abs(d[3] - 0.4) < myeps
    @test abs(d[4] - 0.4) < myeps
    @test abs(d[5] - 2.0) < myeps
end

@testset "Find nonzero minimum" begin
    arr = [0.0, 2.0, 0.0, 0.0, 0.0, 9.0, 1.0]
    val = find_minimum_nonzero(arr)
    @test val == 1.0
end

@testset "Hadi & Simonoff 1993 - initial subset" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    y[n - 1] = y[n - 1] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    subset = hs93initialset(reg)
    @test !(49 in subset)
    @test !(50 in subset)
end

@testset "Hadi & Simonoff 1993 - basic subset" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    y[n - 1] = y[n - 1] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    initialsubset = hs93initialset(reg)
    basicsubset = hs93basicsubset(reg, initialsubset)
    @test !(49 in basicsubset)
    @test !(50 in basicsubset)
end

@testset "Hadi & Simonoff 1993 - Algorithm" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    outset = hs93(reg)
    @test 50 in outset["outliers"]
end

@testset "Hadi & Simonoff 1993 - with Hadi & Simonoff random data" begin
    reg = createRegressionSetting(@formula(y ~ x1 + x2), hs93randomdata)
    outset = hs93(reg)
    @test 1 in outset["outliers"]
    @test 2 in outset["outliers"]
    @test 3 in outset["outliers"]
end

@testset "Kianifard & Swallow 1989 - Algorithm" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    outset = ks89(reg, alpha=0.1)
    @test 15 in outset
    @test 16 in outset
    @test 17 in outset
    @test 18 in outset
    @test 19 in outset
    @test 20 in outset

    df2 = stackloss
    reg2 = createRegressionSetting(@formula(stackloss ~ airflow + watertemp + acidcond), stackloss)
    outset2 = ks89(reg2)
    @test 4 in outset2
    @test 21 in outset2
end

@testset "Sebert & Montgomery & Rollier 1998 - Algorithm" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    y[n - 1] = y[n - 1] * 2.0
    y[n - 2] = y[n - 2] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    outset = smr98(reg)
    @test 49 in outset
    @test 50 in outset
end

@testset "LMS - Algorithm - Random data" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    y[n - 1] = y[n - 1] * 2.0
    y[n - 2] = y[n - 2] * 2.0
    x[n - 3] = x[n - 3] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    outset = lms(reg)["outliers"]
    @test 48 in outset
    @test 49 in outset
    @test 50 in outset
end

@testset "LMS - Algorithm - Phone data" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    outset = lms(reg)["outliers"]
    @test 15 in outset
    @test 16 in outset
    @test 17 in outset
    @test 18 in outset
    @test 19 in outset
    @test 20 in outset
    @test 21 in outset
end

@testset "LTS - Algorithm - Random data" begin
    # Create simple data
    rng = MersenneTwister(12345)
    n = 50
    x = collect(1:n)
    e = randn(rng, n) .* 2.0
    y = 5 .+ 5 .* x .+ e
    y[n] = y[n] * 2.0
    y[n - 1] = y[n - 1] * 2.0
    y[n - 2] = y[n - 2] * 2.0
    x[n - 3] = x[n - 3] * 2.0
    df = DataFrame(x=x, y=y)
    reg = createRegressionSetting(@formula(y ~ x), df)
    outset = lms(reg)["outliers"]
    @test 48 in outset
    @test 49 in outset
    @test 50 in outset
end

@testset "LTS - Algorithm - Phone data" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    outset = lts(reg)["outliers"]
    @test 15 in outset
    @test 16 in outset
    @test 17 in outset
    @test 18 in outset
    @test 19 in outset
    @test 20 in outset
    @test 21 in outset
end

@testset "MVE - Algorithm - Phone data" begin
    df = phones
    outset = mve(df)["outliers"]
    @test 15 in outset
    @test 16 in outset
    @test 17 in outset
    @test 18 in outset
    @test 19 in outset
    @test 20 in outset
    @test 21 in outset
end

@testset "MCD - Algorithm - Phone data" begin
    df = phones
    outset = mcd(df)["outliers"]
    @test 15 in outset
    @test 16 in outset
    @test 17 in outset
    @test 18 in outset
    @test 19 in outset
    @test 20 in outset
    @test 21 in outset
end

@testset "MVE & LTS Plot - Algorithm - Phone data" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = mveltsplot(reg, showplot=false)
    regulars = result["regular.points"]
    @test 1 in regulars
    @test 2 in regulars
    @test 3 in regulars
    @test 4 in regulars
    @test 5 in regulars
    @test 6 in regulars
    @test 7 in regulars
    @test 8 in regulars
    @test 9 in regulars
    @test 10 in regulars
    @test 11 in regulars
    @test 12 in regulars
    @test 13 in regulars
end

@testset "BCH - Algorithm - Hawkins & Bradu & Kass data" begin
    df = hbk
    reg = createRegressionSetting(@formula(y ~ x1 + x2 + x3), df)
    result = bch(reg)
    regulars = result["basic.subset"]
    for i in 15:75
        @test i in regulars
    end
end

@testset "PY95 - Algorithm - Hawkins & Bradu & Kass data" begin
    df = hbk
    reg = createRegressionSetting(@formula(y ~ x1 + x2 + x3), df)
    result = py95(reg)
    outliers = result["outliers"]
    for i in 1:14
        @test i in outliers
    end
end

@testset "satman2013 - Algorithm - Hawkins & Bradu & Kass data" begin
    df = hbk
    reg = createRegressionSetting(@formula(y ~ x1 + x2 + x3), df)
    result = satman2013(reg)
    outliers = result["outliers"]
    for i in 1:14
        @test i in outliers
    end
end

@testset "satman2015 - Algorithm - Phone data" begin
    @test dominates([1,2,4], [1,2,3])
    @test dominates([1,2,3], [0,1,2])
    @test !dominates([1,2,3], [1,2,3])
    @test !dominates([1,2,3], [1,2,4])

    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = satman2015(reg)
    outliers = result["outliers"]
    for i in 15:20
        @test i in outliers
    end
end

@testset "asm2000 - Algorithm - Phone data" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = asm2000(reg)
    outliers = result["outliers"]
    for i in 15:20
        @test i in outliers
    end
end

@testset "ccf - Algorithm - Phone data" begin
    eps = 0.0001
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = ccf(reg)
    outliers = result["outliers"]
    for i in 15:20
        @test i in outliers
    end

    @test all(isapprox.(result["betas"], [-63.4816, 1.30406], atol=eps, rtol=0.0))
end

@testset "LAD - Algorithm" begin
    eps = 0.0001
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = lad(reg)
    betas = result["betas"]
    @test betas[1] < 0
    @test betas[2] > 0

    df2 = DataFrame(
        x=[1,2,3,4,5,6,7,8,9,10],
        y=[2,4,6,8,10,12,14,16,18,1000]
    )
    reg2 = createRegressionSetting(@formula(y ~ x), df2)
    result2 = lad(reg2)
    betas2 = result2["betas"]
    @test abs(betas2[1] - 0) < eps
    @test abs(betas2[2] - 2) < eps
end

@testset "LAD with (X, y)" begin
    eps = 0.0001
    df2 = DataFrame(
        x=Float64[1,2,3,4,5,6,7,8,9,10],
        y=Float64[2,4,6,8,10,12,14,16,18,1000]
    )
    n = length(df2[!,"x"])
    X = hcat(ones(Float64, n), df2[!,"x"])
    y = df2[!, "y"]
    result2 = lad(X, y)
    betas2 = result2["betas"]
    @test abs(betas2[1] - 0) < eps
    @test abs(betas2[2] - 2) < eps
end

@testset "LTA - Algorithm - Phone data" begin
    eps = 0.0001
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = lta(reg, exact=true)
    betas = result["betas"]
    @test abs(betas[1] - -55.5) < eps
    @test abs(betas[2] -  1.15) < eps
end

@testset "Hadi 1992 - Algorithm - with several case" begin
    phones_matrix = hcat(phones[:, "calls"], phones[:, "year"])
    result = hadi1992(phones_matrix)
    outlier_indices = result["outliers"]
    for i in 15:20
        @test i in outlier_indices
    end

    hbk_matrix = hcat(hbk[:, "x1"], hbk[:, "x2"], hbk[:, "x3"])
    hbk_n, hbk_p = size(hbk_matrix)
    result = hadi1992(hbk_matrix)
    outlier_indices = result["outliers"]
    for i in 15:hbk_n
        @test !(i in outlier_indices)
    end

end

@testset "Hadi 1994 - Algorithm - with several case" begin
    phones_matrix = hcat(phones[:, "calls"], phones[:, "year"])
    result = hadi1994(phones_matrix)
    outlier_indices = result["outliers"]
    for i in 15:20
        @test i in outlier_indices
    end

    hbk_matrix = hcat(hbk[:, "x1"], hbk[:, "x2"], hbk[:, "x3"])
    hbk_n, hbk_p = size(hbk_matrix)
    result = hadi1994(hbk_matrix)
    outlier_indices = result["outliers"]
    for i in 15:hbk_n
        @test !(i in outlier_indices)
    end

end

@testset "gwcga - Algorithm - phone data" begin
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = gwcga(reg)
    clean_indices = result["clean.subset"]
    for i in 15:20
        @test !(i in clean_indices)
    end
end

@testset "Compact Genetic Algorithm" begin
    Random.seed!(12345)
    function fcost(bits)
        return sum(bits)
    end
    result = cga(chsize=10, costfunction=fcost, popsize=100)
    for element in result
        @test element == 0
    end
end

@testset "Floating-point Genetic Algorithm" begin
    Random.seed!(12345)
    fcost(genes) = (genes[1] - 3.14159265)^2 + (genes[2] - 2.71828)^2
    epsilon = 10.0^(-6.0)
    popsize = 30
    chsize = 2
    mins = [-100.0, -100.0]
    maxs = [100.0, 100.0]
    result = ga(popsize, chsize, fcost, mins, maxs, 0.90, 0.05, 1, 100)
    best::RealChromosome = result[1]

    @test length(result) == 30
    @test abs(best.genes[1] - 3.14159265) < epsilon
    @test abs(best.genes[2] - 2.71828) < epsilon
end

@testset "Satman(2012) (Csteps and GA based LTS) Algorithm - Phones data" begin
    Random.seed!(12345)
    epsilon = 10.0^(-3.0)
    df = phones
    reg = createRegressionSetting(@formula(calls ~ year), df)
    result = galts(reg)
    betas = result["betas"]
    @test abs(betas[1] - (-56.5219)) < epsilon
    @test abs(betas[2] - 1.16488) < epsilon
end

@testset "RANSAC (1981) Algorithm - Paper example" begin
    Random.seed!(12345)
    df = DataFrame(y=[0,1,2,3,3,4,10], x=[0,1,2,2,3,4,2])
    reg = createRegressionSetting(@formula(y ~ x), df)
    result = ransac(reg, t=0.8, w=0.85)
    @test result == [7]
end

@testset "covratio - phone data" begin
    eps = 0.00001
    knownvals = [1.2945913799871505, 1.2700457384289985, 1.247094858041991,
                1.2255082894941416, 1.2056307474474428, 1.1872216590773372,
                1.1702669047202972, 1.1544784969566473, 1.1415354222939595,
                1.129459077701347, 1.1188675334881035, 1.1092977882351849,
                1.1003495479103549, 1.1024131168300708, 1.0438195753346975,
                1.049139979063326, 1.0015616430250294, 0.9497339678816025,
                0.8530116434438035, 0.6899325070071758, 1.1297343653996992,
                1.0347237388802213, 1.039659281049157, 1.0393513929604028
    ]
    setting = createRegressionSetting(@formula(calls ~ year), phones)
    n = length(knownvals)
    for i in 1:n
        calculated = covratio(setting, i)
        @test abs(knownvals[i] - calculated) < eps
    end
end

@testset "dfbeta - phone data" begin
    eps = 0.00001
    reg = createRegressionSetting(@formula(calls ~ year), phones)
    n, p = size(phones)
    knownvalues = [9.6439157 -0.14686166; 5.3459460 -0.08092134;
                    1.6258961 -0.02443345; -0.6866294  0.01022725;
                    -2.7169197  0.04002009; -4.1910124  0.06085238;
                    -5.1029254  0.07267870; -5.5535000  0.07697356;
                    -5.2512176  0.06983887; -4.6721589  0.05791933;
                    -3.6868718  0.03945523; -2.3254391  0.01478033;
                    -0.5673087 -0.01652355; 1.4653937 -0.04958099;
                    -5.4474441  0.12867978; -8.6540162  0.18101030;
                    -14.6631749  0.28835085; -22.0168113  0.41708081;
                    -32.9443226  0.60863505; -49.0851269  0.89065754;
                    22.9820710 -0.41140246; 39.1639294 -0.69370540;
                    45.7655562 -0.80379984; 53.6862082 -0.93638735;]
    for i in 1:n
        for j in 1:p
            dfbetaresult = dfbeta(reg, i)
            @test abs(dfbetaresult[j] - knownvalues[i,j]) < eps
        end
    end
end

@testset "Atkinson 1994 - Algorithm" begin
    df = stackloss
    reg = createRegressionSetting(@formula(stackloss ~ airflow + watertemp + acidcond), stackloss)
    result = atkinson94(reg)
    @test result["outliers"] == [1, 3, 4, 21]
end

@testset "Hadi Measure" begin
    eps = 0.0001
    setting = createRegressionSetting(@formula(calls ~ year), phones)
    knowncooks = [0.19101337, 0.16141894, 0.13677220, 0.11673486, 0.10058688, 0.08815274, 0.07913586,
                 0.07353075, 0.06965672, 0.06906446, 0.07108612, 0.07605141, 0.08428584, 0.08612449,
                0.15005170, 0.15622736, 0.22082452, 0.29817740, 0.44310105, 0.72642163, 0.19942889,
                0.33018728, 0.36907744, 0.41937743 ]
    hm = hadimeasure(setting)["measure"]
    @test map((x, y) -> abs(x - y) < eps, hm, knowncooks) == trues(24)
end

@testset "Imon 2005 - Algorithm" begin
    Random.seed!(12345)
    reg = createRegressionSetting(@formula(y ~ x1 + x2 + x3), hbk)
    result = imon2005(reg)
    @test result["outliers"] == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
end

@testset "BACON 2000 - Algorithm" begin
    df = stackloss
    reg = createRegressionSetting(@formula(stackloss ~ airflow + watertemp + acidcond), stackloss)
    result = bacon(reg, m=12)
    @test result == [1, 3, 4, 21]
end

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
# # 2D decaying turbulence
#
# A simulation of decaying two-dimensional turbulence.

using Pkg; Pkg.instantiate()

using FourierFlows, Printf, Random, JLD2

using Random: seed!
using FFTW: rfft, irfft
using LinearAlgebra: ldiv!
using FourierFlows: parsevalsum

import GeophysicalFlows.TwoDNavierStokes
import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy
import GeophysicalFlows: peakedisotropicspectrum

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

α = 1                      # aspect ratio parameter α = Ly / Lx
ny, Ly  = 2*2*256, 2π            # grid resolution and domain length
nx, Lx  = Int(2*2*256 / α), Ly / α

## Then we pick the time-stepper parameters
	 dt =0.5*1e-3  # timestep
nsubs = 1000    # number of steps between each plot

 ν = 0.2*1e-4

initial_data_wavenumber = 2.0
initial_data_bandwidth = 0.5

k₀ = initial_data_wavenumber

 tfinal = 600.0 #/ (ν * k₀^2)
 nsteps = Int(round(tfinal / dt))

 grid = TwoDGrid(dev, nx, Lx, ny, Ly)
 x, y = gridpoints(grid)

 # No Forcing
 calcF!(args...) = nothing

# Initial condition

# seed!(1234)

E0 = 0.1

function random_initial_condition(grid, initial_data_wavenumber, initial_data_bandwidth; energy=0.1)
	K = @. sqrt(grid.Krsq)
	forcing_spectrum = @. exp(-(K - initial_data_wavenumber)^2 / (2*initial_data_bandwidth^2))

	energy_before_normalization = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)

	forcing_spectrum *= energy / energy_before_normalization

	ξ = exp.(2π * im * rand(eltype(grid), size(K)))

	qh = @. sqrt(forcing_spectrum) * ξ
	q  = irfft(qh, grid.nx)

	return q
end

ζ₀ = random_initial_condition(grid, initial_data_wavenumber, initial_data_bandwidth; energy=E0)

# confirm that energy is E0
# prob = TwoDNavierStokes.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, ν=ν, dt=dt, calcF=calcF!, stepper="ETDRK4")
# TwoDNavierStokes.set_ζ!(prob, ζ₀)
# dx, dy = Lx/nx, Ly/ny
# println(sum(@. (prob.vars.u^2 + prob.vars.v^2)/2 * dx * dy / (Lx * Ly)))

# plot initial condition
# using Plots
# p =heatmap(grid.x, grid.y, ζ₀', aspectratio=1)
# display(p)

include("simulation.jl")

max_stars_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_issues_repo_name: navidcy/TwoDTurbulence-EulerLimit

max_forks_repo_name: navidcy/TwoDTurbulence-EulerLimit


text:
using Dates
using JLD2

#Functions connected to the used Algorithm, but not directly to other domains:
#-------------------------------------------------------------------------------
#Return the L2Reg norm of a vector:
@everywhere L2Reg(x) = mean(x .* x)

#Return the Generalized Advantage Estimation:
@everywhere function td_lambda!(params::Params,traj::Trajectory)
    l_traj = traj.len
    γ      = params.γ_disc
    λ      = params.λ_critic

    Val = traj.val_t
    r_t = traj.r_t

    #Calculate the local TD-error:
    delta_V_i = zeros(l_traj)
    for i_t = 1:l_traj-1
        delta_V_i[i_t] = -Val[i_t] + r_t[i_t] + γ*Val[i_t+1]
    end

    #Discount the TD-Error as the Advantage
    R_t = vec(discount(delta_V_i,λ*γ))

    #Write result in trajectory
    traj.R_t = R_t
end

#Discount the Rewards
@everywhere function discount(rewards, fac)
    discounted = zeros(Float64, length(rewards), 1)
    discounted[end] = rewards[end]

    for i=(length(rewards)-1):-1:1
        discounted[i] = rewards[i] + fac * discounted[i+1]
    end
    return discounted
end

#Better: Do it after the episode was built:
function loss_a3c_mod(θ::Any,params::Params,traj::Trajectory)
    #Take params from struct:
    l_traj = params.l_traj
    γ      = params.γ_disc
    λ      = params.λ_actor

    #Init
    Err_actor  = 0
    Err_critic = 0
    s_t = traj.s_t[1,:]
    val_old    = 0;
    A_t        = 0;

    #Hope: improved peroformance by going backwards...
    for i_t = l_traj-1:-1:1
        #Take from Struct
        s_t = traj.s_t[i_t,:]
        a_t = traj.a_t[i_t]
        r_t = traj.r_t[i_t]

        out_t, val_new = mlp_nlayer(θ,s_t,params)          #Evaluate the Network
        d_val          = r_t + γ*val_old - val_new         #Get local TD-error
        val_old        = val_new                           #Update old error
        A_t            = d_val + (γ*λ)*A_t                 #GAE update
        Err_actor     += logpdf(a_t,out_t[1],params) * A_t #Actor correction
        Err_critic    += 0.5 * (A_t * A_t)                 #Critic correction
    end

    #Return the final Error_Sum
    return Err_critic / l_traj
end

#Function to normalize the inputs of the NN => Chan 1982
@everywhere function update_variance!(rms1::RunMeanStd,rms2::RunMeanStd)
    #Update the mean
    δ = rms1.mean - rms2.mean
    new_mean = (rms1.count .* rms1.mean .- rms2.count .* rms2.mean) ./ (rms1.count + rms2.count)

    #Update the Variance
    m_1 = rms1.var .* (rms1.count - 1)
    m_2 = rms2.var .* (rms2.count - 1)
    M = m_1 .+ m_2 .+ (δ.*δ) .* rms1.count .* rms2.count ./ (rms1.count .+ rms2.count)
    new_count    = rms1.count + rms2.count
    new_variance = M ./ (new_count -1)
    new_std      = sqrt.(new_variance)

    #write new data into rms1:
    rms1.count = min(new_count,1e8) #limit to avoid static behaviour
    rms1.mean  = new_mean
    rms1.var   = new_variance
    rms1.std   = new_std
end

#Function to simply normalize the reward according to the current RMS
@everywhere function normalize_rms(rms::RunMeanStd, s_t)
    #Apply the obersvation-filter:
    return clamp.( (s_t .- rms.mean) ./ (rms.std .+ 1e-2), -5.0, 5.0)
end

#Clip_function
@everywhere function clip(x, lo::Real, hi::Real)
    if x < lo
        return lo
    elseif x > hi
        return hi
    else
        return x
    end
end

#https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm => Two-pass algorithm
#It avoids catastrophic cancellation
@everywhere function two_pass_variance!(all_traj::Traj_Data,rms::RunMeanStd, params::Params)
    #1. Init
    l_state  = length(params.s_init)
    n_traj   = params.n_traj
    x_mean   = zeros(l_state)
    x_sum    = zeros(l_state)
    x_sum_2  = zeros(l_state)
    n        = 0

    # Get the mean:
    for i_traj = 1:n_traj
        l_traj = all_traj.traj[i_traj].len
        x_sum += vec(sum(all_traj.traj[i_traj].s_t2[1:l_traj,:],dims = 1))
        n     += all_traj.traj[i_traj].len
    end
    x_mean = x_sum ./ n

    # Squared error from mean summed:
    for i_traj = 1:n_traj
        l_traj = all_traj.traj[i_traj].len
        for i_t = 1: l_traj
            s_t = all_traj.traj[i_traj].s_t2[i_t,:]
            x_sum_2 += (s_t .- x_mean) .* (s_t .- x_mean)
        end
    end

    #Calculate the final variance:
    variance = x_sum_2 ./ (n - 1)

    #Write data in struct
    rms.count = n
    rms.mean  = x_mean
    rms.var   = variance
    rms.std   = sqrt.(variance)
end

#Build minibatches from attained data:
#Idea Construct Batch_data as partial Trajectories!
@everywhere function build_minibatches(all_traj::Traj_Data,params::Params)
    batch_size = params.batch_size

    timesteps = 0
    #Build current number of timesteps for this iteration
    for i_traj = 1:params.n_traj
        timesteps += all_traj.traj[i_traj].len
    end
    all_traj.timesteps = timesteps

    #Determine number of minibatches:
    number_batches = Int(floor(timesteps / batch_size))

    #Now Initialize new Traj_Data containing only minibatches
    #---------------------------------------------------
    mini_batch = Traj_Data()
    mini_batch.traj = Array{Trajectory,1}(undef,number_batches)
    mini_batch.timesteps = number_batches*batch_size
    mini_batch.n_traj = number_batches

    #Create instance of trajectory:
    traj = init_traj(params,batch_size)

    #Fill the big struct with the small one:
    for i_traj = 1:number_batches
        mini_batch.traj[i_traj] = copy_traj(traj,params)
    end
    #----------------------------------------------------

    #Now fill the minibatch data struct with the data contained from the Trajectories:
    i_traj = 1
    i_tt = 1
    org_len = all_traj.traj[1].len

    for i_batch = 1:number_batches
        mini_batch.traj[i_batch].len = batch_size
        for i_tb = 1:batch_size

            #Copy the data
            copy_to_mini_batch!(mini_batch, all_traj, i_traj, i_batch, i_tt, i_tb)

            #Check if u reached the end of the current trajectory:
            if i_tt == org_len
                i_traj += 1 #switch to new trajectory
                i_tt = 0    #reset the timestep of the trajectory

                #Reset original Trajectory length
                if i_traj <= params.n_traj #avoid undefined step
                    org_len = all_traj.traj[i_traj].len
                end
            end

            #Increase timesteps of original trajectory:
            i_tt += 1
        end
    end

    #Return the mini_batch data struct
    return mini_batch
end

#Copy the data locally to the Minibatch
@everywhere function copy_to_mini_batch!(mini_batch::Traj_Data, all_traj::Traj_Data, i_traj::Int, i_batch::Int, i_tt::Int, i_tb::Int)
    #Copy states
    mini_batch.traj[i_batch].s_t[i_tb,:] = all_traj.traj[i_traj].s_t[i_tt,:]

    #Copy new states
    mini_batch.traj[i_batch].s_t2[i_tb,:] = all_traj.traj[i_traj].s_t2[i_tt,:]

    #Copy actions
    mini_batch.traj[i_batch].a_t[i_tb,:] = all_traj.traj[i_traj].a_t[i_tt,:]

    #Copy probability of actions:
    mini_batch.traj[i_batch].p_a_t[i_tb,:] = all_traj.traj[i_traj].p_a_t[i_tt,:]

    #Copy reward
    mini_batch.traj[i_batch].r_t[i_tb] = all_traj.traj[i_traj].r_t[i_tt]

    #Copy Value Result
    mini_batch.traj[i_batch].val_t[i_tb,:] = all_traj.traj[i_traj].val_t[i_tt,:]

    #Copy Output of the Network
    mini_batch.traj[i_batch].out_t[i_tb,:] = all_traj.traj[i_traj].out_t[i_tt,:]

    #Copy local Advantage
    mini_batch.traj[i_batch].A_t[i_tb] = all_traj.traj[i_traj].A_t[i_tt]

    #Copy local discounted Reward
    mini_batch.traj[i_batch].R_t[i_tb] = all_traj.traj[i_traj].R_t[i_tt]
end

#Save the data for analysis
function save_data(res::Array{Float64,1}, prec::Array{Float64,1}, θa::Any, θv::Any, all_traj::Traj_Data, i_iter::Int)
    filename = string("Results/Actor_Critic_",Dates.today(),"iteartion_",i_iter,".jld2")
    save(filename, "res", res, "prec", prec, "θa", θa, "θv", θv)
end

#This function pumps the new data into the big replay ring-buffer
function pump_buffer!(buffer::Traj_Data,all_traj::Traj_Data, params::Params)
    current_buffer_size = buffer.n_traj
    new_data_size       = all_traj.n_traj
    max_size            = params.buffer_size

    #If buffer not full simply fill:
    if current_buffer_size + new_data_size < max_size
        buffer.n_traj = current_buffer_size + new_data_size
        append!(buffer.traj,all_traj.traj)
        buffer.n_traj = current_buffer_size + new_data_size

    #If buffer not full, but would run full
    elseif current_buffer_size < max_size
        #First calculate the required data to completely fill the buffer
        app_size = max_size - current_buffer_size
        #simply fill to have concluded the buffer
        append!(buffer.traj,all_traj.traj[1:app_size])
        #now get the size required for shifting
        shift_size = new_data_size - app_size
        #now shift the buffer to delete old samples
        buffer.traj[1:end-shift_size] = buffer.traj[shift_size+1:end]
        #fill the "free" space with the still unused free samples
        buffer.traj[end-shift_size+1:end] = all_traj.traj[app_size+1:end]
        #set new buffer size
        buffer.n_traj = max_size

        print("--Buffer is now full!: ", buffer.n_traj,"--")

    #If the buffer is full and the new minibatch is not larger shift
    elseif new_data_size < max_size
        shift_size = new_data_size
        #now shift the buffer to delete old samples
        buffer.traj[1:end-shift_size] = buffer.traj[shift_size+1:end]
        #Now fill the "free" space with the new samples
        buffer.traj[end-shift_size+1:end] = all_traj.traj
        #set new buffer size
        buffer.n_traj = max_size

    #If the new buffer is larger than the oberall buffer just fill partly
    else
        #fill the buffer with the partly new data
        buffer.traj = all_traj.traj[1:max_size]
        #set new buffer size
        buffer.n_traj = max_size
    end
    #set timestepsize of buffer
    buffer.timesteps = buffer.n_traj * params.batch_size
end

"""
Simply trashy function to get the baseline
"""
function simple_baseline!(all_traj::Traj_Data, params::Params)
    n_traj = all_traj.n_traj
    γ = params.γ_disc
    max_len = 0
    #Get R_t for all timesteps
    for i_traj = 1:n_traj
        len = all_traj.traj[i_traj].len
        #update max_len
        if len > max_len
            max_len = len
        end
        R_to = all_traj.traj[i_traj].r_t[1]
        all_traj.traj[i_traj].R_t[1] = R_to
        for i_t = 2:len
            r_t = all_traj.traj[i_traj].r_t[i_t]
            all_traj.traj[i_traj].R_t[i_t] = γ*R_to + r_t
        end
    end
    #Get mean R_t = b_t
    for i_t = 1:max_len
        i_count   = 0
        R_all     = 0
        for i_traj = 1:all_traj.n_traj
            len = all_traj.traj[i_traj].len
            #Count if traj is long enough
            if i_t <= len
                i_count += 1
                R_all += all_traj.traj[i_traj].R_t[i_t]
            end
        end
        b_t = R_all / i_count
        for i_traj = 1:all_traj.n_traj
            len = all_traj.traj[i_traj].len
            #Count if traj is long enough
            if i_t <= len
                R_loc = all_traj.traj[i_traj].R_t[i_t]
                A_t   = R_loc - b_t
            end
        end
    end
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using SpectralElements
using Test, SafeTestsets

@testset "SpectralElements.jl" begin

    @time @safetestset "Field" begin include("Field.jl") end
    @time @safetestset "Operators" begin include("Operators.jl") end

#   @time @safetestset "Examples" begin include("examples.jl") end
end


text:
@test factorial(7,3) == 7*6*5*4
@test_throws DomainError factorial(3,7)
@test_throws DomainError factorial(-3,-7)
@test_throws DomainError factorial(-7,-3)
#JuliaLang/julia#9943
@test factorial(big(100), (80)) == 1303995018204712451095685346159820800000
#JuliaLang/julia#9950
@test_throws OverflowError factorial(1000,80)

# derangement
@test derangement(4) == subfactorial(4) == 9
@test derangement(24) == parse(BigInt,"228250211305338670494289")

# partialderangement
@test partialderangement(7, 3) == 315
@test_throws DomainError partialderangement(8, 9)
@test_throws DomainError partialderangement(-8, 0)

# doublefactorial
@test doublefactorial(70) == parse(BigInt,"355044260642859198243475901411974413130137600000000")
@test_throws DomainError doublefactorial(-1)

# hyperfactorial
@test hyperfactorial(8) == parse(BigInt,"55696437941726556979200000")
@test hyperfactorial(0) == parse(BigInt,"1")
@test hyperfactorial(1) == parse(BigInt,"1")
@test hyperfactorial(2) == parse(BigInt,"4")

# multifactorial
@test multifactorial(40,2) == doublefactorial(40)
@test_throws DomainError multifactorial(-1,1)

# multinomial
@test multinomial(1,4,4,2) == 34650

# primorial
@test primorial(17) == 510510
@test_throws DomainError primorial(-1)


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
A MirrorProxParameters struct specifies the parameters for solving the saddle
point formulation of an problem using mirror prox.

Quadratic Programming Problem (see quadratic_programming.jl):
minimize 1/2 * x' * objective_matrix * x + objective_vector' * x
         + objective_constant

s.t. constraint_matrix[1:num_equalities, :] * x =
     right_hand_side[1:num_equalities]

     constraint_matrix[(num_equalities + 1):end, :] * x >=
     right_hand_side[(num_equalities + 1):end, :]

     variable_lower_bound <= x <= variable_upper_bound

Equivalent saddle point problem:
    min_x max_y phi(x, y)
where
    phi(x, y) = objective_constant + objective_vector' x + right_hand_side' y
                + 1/2 * x' * objective_matrix * x
                - y' constraint_matrix x
    variable_lower_bound <= x <= variable_upper_bound
    y[(num_equalities + 1):end] >= 0

Solve using saddle point mirror prox:
https://blogs.princeton.edu/imabandit/2013/04/23/orf523-mirror-prox/
using mirror map (primal_weight/2) ||x||_X^2 + 1/(2 primal_weight) ||y||_Y^2
where X and Y are diagonal matrices. If `diagonal_scaling`, X[i,i] is set as the
l2 norm of the i-th column of constraint_matrix, and Y[j,j] is set as the l2
norm of the j-th column of constraint_matrix. The intuition behind this setup
is for mirror prox with constant step-size 1, we would like to make sure
[X,-A';-A Y] (A is the constraint_matrix) is positive semi-definite,
which is satisfied automatically for problem with the above X and Y. If not
`diagonal_scaling`, both X and Y are set as identity matrix.

The parameter primal_weight is adjusted dynamically so the two terms contribute
roughly equally.
"""
struct MirrorProxParameters
  """
  Number of L_infinity Ruiz rescaling iterations to apply to the constraint
  matrix. Zero disables this rescaling pass.
  """
  l_inf_ruiz_iterations::Int

  """
  If true, applies L2 norm rescaling after the Ruiz rescaling.
  """
  l2_norm_rescaling::Bool

  """
  If not `nothing`, runs Pock-Chambolle rescaling with the given alpha exponent
  parameter.
  """
  pock_chambolle_alpha::Union{Float64,Nothing}

  """
  Used to bias the computation of the primal/dual balancing parameter
  primal_weight. Must be positive. A value of 1 balances primal and dual
  equally. Somewhat smaller values, e.g. 0.4, seem to yield better dual
  solutions at the cost of increased primal infeasibility.
  """
  primal_importance::Float64

  """
  Use weighted norms of rows/columns as the Bregman divergence so that it
  rescales each primal and dual variable. The allowed values are "l1", "l2", and
  "off".
  If `diagonal_scaling`=l1, we set X[i,i]=||K[:,i]||_1 + ||Q[:,i]||_1
  and Y[j,j]=||K[j,:]||_1;
  If `diagonal_scaling`=l2, we set X[i,i]=sqrt(||K[:,i]||_2^2+||Q[:,i]||_2^2)
  and Y[j,j]=||K[j,:]||_2;
  If `diagonal_scaling`=off, we set X[i,i]=1 and Y[j,j]=1;
  """
  diagonal_scaling::String

  """
  If >= 4 a line of debugging info is printed during some iterations. If >= 2
  some info is printed about the final solution.
  """
  verbosity::Int64

  """
  Whether to record an IterationStats object.
  """
  record_iteration_stats::Bool

  """
  Check for termination with this frequency (in iterations).
  """
  termination_evaluation_frequency::Int32

  """
  The termination criteria for the algorithm.
  """
  termination_criteria::TerminationCriteria

  """
  Parameters that control when the algorithm restarts and whether it resets to
  the average or the current iterate. Also, controls the primal weight updates.
  """
  restart_params::RestartParameters
end
# TODO: Add validate function to check parameters are correct.

struct MirrorProxVector
  data::Vector{Float64}
  primal_size::Int64
end

mutable struct MirrorProxProblem
  problem::QuadraticProgrammingProblem

  """
  [0 A^T;
   -A 0]
  where A = problem.constraint_matrix
  """
  combo_matrix::SparseMatrixCSC{Float64,Int64}

  params::MirrorProxParameters

  """
  Balances the contribution of primal and dual distances to the mirror map.
  """
  primal_weight::Float64

  """
  The weights of the weighted l2 norm (as the Bregman divergence) used in the
  mirror map.
  """
  mirror_map_scaling::MirrorProxVector
end

# This overloads the "-" operator.
function Base.:-(vec1::MirrorProxVector, vec2::MirrorProxVector)
  @assert vec1.primal_size == vec2.primal_size
  return MirrorProxVector(vec1.data - vec2.data, vec1.primal_size)
end

function Base.:+(vec1::MirrorProxVector, vec2::MirrorProxVector)
  @assert vec1.primal_size == vec2.primal_size
  return MirrorProxVector(vec1.data + vec2.data, vec1.primal_size)
end

function Base.:*(scalar::Float64, vec::MirrorProxVector)
  return MirrorProxVector(scalar * vec.data, vec.primal_size)
end

function Base.:/(vec::MirrorProxVector, scalar::Float64)
  return MirrorProxVector(vec.data / scalar, vec.primal_size)
end

function dot(vec1::MirrorProxVector, vec2::MirrorProxVector)
  @assert vec1.primal_size == vec2.primal_size
  return vec1.data' * vec2.data
end

""" Like .* but that's not entirely trivial to implement for a custom type."""
function pointwise_mult(vec1::MirrorProxVector, vec2::MirrorProxVector)
  @assert vec1.primal_size == vec2.primal_size
  return MirrorProxVector(vec1.data .* vec2.data, vec1.primal_size)
end

""" Like ./ but that's not entirely trivial to implement for a custom type."""
function pointwise_div(vec1::MirrorProxVector, vec2::MirrorProxVector)
  @assert vec1.primal_size == vec2.primal_size
  return MirrorProxVector(vec1.data ./ vec2.data, vec1.primal_size)
end

function make_sp_vec(
  primal::AbstractVector{Float64},
  dual::AbstractVector{Float64},
)
  return MirrorProxVector(vcat(primal, dual), length(primal))
end

function zero_sp_vec(problem::QuadraticProgrammingProblem)
  return make_sp_vec(
    zeros(length(problem.variable_lower_bound)),
    zeros(length(problem.right_hand_side)),
  )
end

function const_sp_vec(
  primal_val::Float64,
  dual_val::Float64,
  problem::QuadraticProgrammingProblem,
)
  return make_sp_vec(
    primal_val * ones(length(problem.variable_lower_bound)),
    dual_val * ones(length(problem.right_hand_side)),
  )
end

function primal_part(vec::MirrorProxVector)
  return view(vec.data, 1:vec.primal_size)
end

function dual_part(vec::MirrorProxVector)
  return view(vec.data, (vec.primal_size+1):length(vec.data))
end

function combo_vec(vec::MirrorProxVector)
  return vec.data
end

function mirror_map(vec::MirrorProxVector, problem::MirrorProxProblem)
  return 0.5 * dot(vec, pointwise_mult(problem.mirror_map_scaling, vec))
end

"""The norm that mirror_map is strongly convex with respect to.
This isn't used but is here for documentation."""
function mirror_map_norm(vec::MirrorProxVector, problem::MirrorProxProblem)
  # This is equivalent to:
  # return sqrt(2 * mirror_map(vec, problem))
  return sqrt(dot(vec, pointwise_mult(problem.mirror_map_scaling, vec)))
end

"""The dual of mirror_map_norm().
This isn't used but is here for documentation."""
function dual_norm(vec::MirrorProxVector, problem::MirrorProxProblem)
  return sqrt(dot(vec, pointwise_div(vec, problem.mirror_map_scaling)))
end

""" Bregman divergence of the mirror map. """
function bregman(
  vec1::MirrorProxVector,
  vec2::MirrorProxVector,
  p::MirrorProxProblem,
)
  # This is specific to this mirror map of course.
  return mirror_map(vec1 - vec2, p)
end

function add_to_solution_weighted_average(
  solution_weighted_average::SolutionWeightedAverage,
  current_solution::MirrorProxVector,
  weight::Float64,
)
  add_to_solution_weighted_average(
    solution_weighted_average,
    primal_part(current_solution),
    dual_part(current_solution),
    weight,
  )
end

""" Returns a 4-element array with entries that sum to the saddle point
objective phi."""
function phi_breakdown(vec::MirrorProxVector, problem::MirrorProxProblem)
  problem = problem.problem
  return [
    problem.objective_constant,
    problem.objective_vector' * primal_part(vec),
    problem.right_hand_side' * dual_part(vec),
    -dual_part(vec)' * problem.constraint_matrix * primal_part(vec),
  ]
end

"""
Returns (gradient of phi with respect to x, -gradient of phi with respect to y).
That's c - A'y, Ax - b.
This is used in place of the subgradient in saddle point mirror prox. It can be
thought of as the concatenation of the loss vectors for the two players of the
zero-sum game.
"""
function pseudo_gradient(solution::MirrorProxVector, p::MirrorProxProblem)
  problem = p.problem
  # TODO: cache this matrix product and reuse it for
  # corrected_dual_obj().
  result = p.combo_matrix' * combo_vec(solution)
  # result is now (-A' y, A x).
  for indx in 1:solution.primal_size
    result[indx] += problem.objective_vector[indx]
  end
  result[1:solution.primal_size] +=
    problem.objective_matrix * primal_part(solution)
  for indx in 1:length(problem.right_hand_side)
    result[solution.primal_size+indx] -= problem.right_hand_side[indx]
  end

  return MirrorProxVector(result, solution.primal_size)
end

function corrected_dual_obj(solution::MirrorProxVector, p::MirrorProxProblem)
  return corrected_dual_obj(
    p.problem,
    primal_part(solution),
    dual_part(solution),
  )
end

"""
Returns the proximal step
argmin_z step_size * gradient' * z + bregman(z, initial_solution)
where z ranges over values satisfying the bounds on the primal and dual
variables.

For the particular mirror map and constraint set we're using this minimization
problem is:
min_z step_size gradient' z + 0.5(z - initial_solution) D (z - initial_solution)
where D is a diagonal matrix with problem.mirror_map_scaling on the diagonal.
This is equivalent (up to additive constant) to
min_z 0.5 (z - initial_solution + step_size D^-1 gradient)' D
(z - initial_solution + step_size D^-1 gradient)
This separates into separate trivial optimization problems for each variable.
"""
function proximal_step(
  p::MirrorProxProblem,
  initial_solution::MirrorProxVector,
  gradient::MirrorProxVector,
  step_size::Float64,
)
  problem = p.problem
  solution =
    initial_solution - step_size * pointwise_div(gradient, p.mirror_map_scaling)
  primal = primal_part(solution)
  dual = dual_part(solution)
  project_primal!(primal, problem)
  project_dual!(dual, problem)
  return solution
end

""" Returns a tuple with stats on how many variables are at the bounds and how
many of the matrix non-zeros would remain after eliminating the variables that
are at their bounds. We say a variable is "active" if it is not at its bounds.
"""
function active_variable_stats(
  solution::MirrorProxVector,
  problem::MirrorProxProblem,
)
  problem = problem.problem
  num_prim_lb = 0
  num_prim_active = 0
  num_prim_ub = 0
  primal = primal_part(solution)
  for idx in 1:length(primal)
    if primal[idx] >= problem.variable_upper_bound[idx]
      num_prim_ub += 1
    else
      if primal[idx] <= problem.variable_lower_bound[idx]
        num_prim_lb += 1
      else
        num_prim_active += 1
      end
    end
  end
  num_dual_lb = 0
  num_dual_active = 0
  dual = dual_part(solution)
  for idx in inequality_range(problem)
    if dual[idx] <= 0
      num_dual_lb += 1
    else
      num_dual_active += 1
    end
  end

  # Number of nonzeros where both the row and column variables are active.
  num_active_nonzeros = 0

  rows = SparseArrays.rowvals(problem.constraint_matrix)
  m, n = size(problem.constraint_matrix)
  for col in 1:n
    if primal[col] == problem.variable_upper_bound[col] ||
       primal[col] == problem.variable_lower_bound[col]
      continue
    end
    for j in SparseArrays.nzrange(problem.constraint_matrix, col)
      row = rows[j]
      if row in equality_range(problem) || dual[row] > 0
        num_active_nonzeros += 1
      end
    end
  end

  return (
    num_prim_lb,
    num_prim_active,
    num_prim_ub,
    num_dual_lb,
    num_dual_active,
    num_active_nonzeros,
  )
end

"""
logs info to STDOUT giving detailed information on the current iteration.
"""
function mirror_prox_specific_log(
  params::MirrorProxParameters,
  p::MirrorProxProblem,
  iteration::Int64,
  accepted_iterations::Int64,
  acceptable_nonlinearity::Bool,
  step_size::Float64,
  inv_stepsize_required::Float64,
  current_solution::MirrorProxVector,
)
  phi_components = phi_breakdown(current_solution, p)
  primal_solution = primal_part(current_solution)
  dual_solution = dual_part(current_solution)
  primal_residual = compute_primal_residual(p.problem, primal_solution)
  Printf.@printf(
    "    (%5d): inv_step_size=%9g %s req'd=%9g ",
    accepted_iterations,
    1.0 / step_size,
    (acceptable_nonlinearity ? ">=" : " <"),
    inv_stepsize_required
  )
  # Provides information on *scaled norms*. Useful for diagnosing poor
  # variable scaling.
  Printf.@printf(
    "l2=(%.2e, %.2e) ",
    norm(primal_solution, 2),
    norm(dual_solution, 2)
  )
  Printf.@printf(
    "l_inf=(%.2e, %.2e)\n",
    norm(primal_solution, Inf),
    norm(dual_solution, Inf)
  )
  Printf.@printf(
    "    primal_weight=%9g phi=%12g=",
    p.primal_weight,
    sum(phi_components)
  )
  println(phi_components)
  if params.verbosity >= 6
    stats = active_variable_stats(current_solution, p)
    Printf.@printf(
      "    primal vars [@lb: %d active: %d @ub: %d] ",
      stats[1],
      stats[2],
      stats[3]
    )
    Printf.@printf(
      "dual eq vars: %d dual ineq vars [@lb: %d active: %d] ",
      length(equality_range(p.problem)),
      stats[4],
      stats[5]
    )
    Printf.@printf(
      "active mat nz: %6.3f%%\n",
      100.0 * stats[6] / SparseArrays.nnz(p.problem.constraint_matrix)
    )
  end
end

function primal_obj(solution::MirrorProxVector, p::MirrorProxProblem)
  problem = p.problem
  return problem.objective_constant +
         problem.objective_vector' * primal_part(solution) +
         0.5 *
         primal_part(solution)' *
         problem.objective_matrix *
         primal_part(solution)
end

"""
Creates a MirrorProxProblem including defining static problem scalings
based on the QuadraticProgrammingProblem.
"""
function initialize_saddle_point_problem(
  problem::QuadraticProgrammingProblem,
  params::MirrorProxParameters,
  combo_mat::SparseMatrixCSC{Float64,Int64},
)
  # With mirror map (primal_weight/2) ||x||_X^2 + 1/(2 primal_weight) ||y||_Y^2
  # the theory holds when
  # [primal_weight X, 0, -Q, -A'
  #  0, Y / primal_weight, A, 0
  #  -Q, A', primal_weight X, 0
  #  -A, 0, 0, Y / primal_weight]
  # is positive semi-definte. This is automatically satisfied when we set
  # X[i,i]=||K[:,i]||_1 + ||Q[:,i]||_1 and Y[j,j]=||K[j,:]||_1. In practice,
  # using ||.||_2 instead of ||.||_1 to define X and Y often improves
  # the convergence speed for mirror prox.
  if params.diagonal_scaling == "l2"
    # If `diagonal_scaling`=l2, we set
    # X[i,i]=sqrt(||K[:,i]||_2^2+||Q[:,i]||_2^2) and Y[j,j]=||K[j,:]||_2.
    primal_scalings =
      sqrt.(
        sum(problem.constraint_matrix .^ 2, dims = 1) +
        sum(problem.objective_matrix .^ 2, dims = 1),
      )
    dual_scalings = sqrt.(sum(problem.constraint_matrix .^ 2, dims = 2))
    p = MirrorProxProblem(
      problem,
      combo_mat,
      params,
      params.primal_importance,
      make_sp_vec(vec(primal_scalings), vec(dual_scalings)),
    )
  elseif params.diagonal_scaling == "l1"
    # If `diagonal_scaling`=l1, we set X[i,i]=||K[:,i]||_1+||Q[:,i]||_1 and
    # Y[j,j]=||K[j,:]||_1.
    primal_scalings = vec(
      sum(abs.(problem.constraint_matrix), dims = 1) +
      sum(abs.(problem.objective_matrix), dims = 1),
    )
    dual_scalings = vec(sum(abs.(problem.constraint_matrix), dims = 2))
    p = MirrorProxProblem(
      problem,
      combo_mat,
      params,
      1.0, # primal_weight
      make_sp_vec(vec(primal_scalings), vec(dual_scalings)),
    )
  else
    # If `diagonal_scaling`=off, we set X[i,i]=1 and Y[j,j]=1.
    p = MirrorProxProblem(
      problem,
      combo_mat,
      params,
      1.0, # primal_weight
      const_sp_vec(
        params.primal_importance,
        1.0 / params.primal_importance,
        problem,
      ),
    )
  end

  return p
end

"""
Change the norms that we use to measure the primal and dual variables to reflect
the change in the primal weight
"""
function update_mirror_prox_norms(p::MirrorProxProblem, primal_weight::Float64)
  if p.primal_weight != primal_weight
    problem = p.problem
    primal_size = length(problem.variable_lower_bound)
    dual_size = length(problem.right_hand_side)
    primal_weight_ratio = primal_weight / p.primal_weight
    p.primal_weight = primal_weight

    for i in 1:primal_size
      primal_part(p.mirror_map_scaling)[i] *= primal_weight_ratio
    end
    for i in 1:dual_size
      dual_part(p.mirror_map_scaling)[i] /= primal_weight_ratio
    end
  end
end

"""
Logging while the algorithm is running.
"""
function mirror_prox_display(
  termination_reason::Union{TerminationReason,Bool},
  iteration::Int64,
  params::MirrorProxParameters,
  current_iteration_stats::IterationStats,
  p::MirrorProxProblem,
  accepted_iterations::Int64,
  acceptable_nonlinearity::Bool,
  step_size::Float64,
  inv_stepsize_required::Float64,
  current_solution::MirrorProxVector,
)
  if print_to_screen_this_iteration(
    termination_reason,
    iteration,
    params.verbosity,
    params.termination_evaluation_frequency,
  )
    display_iteration_stats(current_iteration_stats, params.verbosity)
    if iteration > 1 && params.verbosity >= 6
      mirror_prox_specific_log(
        params,
        p,
        iteration,
        accepted_iterations,
        acceptable_nonlinearity,
        step_size,
        inv_stepsize_required,
        current_solution,
      )
    end
  end
end

"""
`optimize(params::MirrorProxParameters,
          original_problem::QuadraticProgrammingProblem)`

Solves a quadratic program using saddle point mirror prox.

# Arguments
- `params::MirrorProxParameters`: parameters.
- `original_problem::QuadraticProgrammingProblem`: the QP to solve.

# Returns
A SaddlePointOutput struct containing the solution found.
"""
function optimize(
  params::MirrorProxParameters,
  original_problem::QuadraticProgrammingProblem,
)
  validate(original_problem)
  # TODO: Split this code into functions to make easier to read.
  qp_cache = cached_quadratic_program_info(original_problem)
  scaled_problem = rescale_problem(
    params.l_inf_ruiz_iterations,
    params.l2_norm_rescaling,
    params.pock_chambolle_alpha,
    params.verbosity,
    original_problem,
  )
  problem = scaled_problem.scaled_qp

  primal_size = length(problem.variable_lower_bound)
  dual_size = length(problem.right_hand_side)
  combo_mat::SparseMatrixCSC{Float64,Int64} = [
    SparseArrays.spzeros(primal_size, primal_size) problem.constraint_matrix'
    (-problem.constraint_matrix) SparseArrays.spzeros(dual_size, dual_size)
  ]
  validate_termination_criteria(params.termination_criteria)
  current_solution = zero_sp_vec(problem)
  # The theory suggests using the inverse of the max singular value as the step
  # size. We could approximate the max singular value using power iteration but
  # that is not only slower but produces step sizes that are too conservative in
  # practice. The step size adjustment algorithm we use decreases the step
  # size faster than it increases it so it's better to err on the side of a
  # bigger initial step size. Instead we use the max absolute value entry in a
  # matrix, which is a lower bound on the max singular value. This seems to work
  # well, although I usually test with Ruiz rescaling pre-processing so this is
  # effectively just initializing the step size to 1.
  step_size = 1.0 / norm(problem.constraint_matrix, Inf)
  solution_weighted_avg = initialize_solution_weighted_average(
    length(primal_part(current_solution)),
    length(dual_part(current_solution)),
  )
  memory_solution = current_solution
  best_dual = -Inf
  accepted_iterations = 0
  iterations_completed = 0
  acceptable_nonlinearity = false
  inv_stepsize_required = NaN
  KKT_PASSES_PER_ITERATION = 2.0

  # Idealized number of KKT passes each time the termination criteria and
  # restart scheme is run. One of these comes from evaluating the gradient at
  # the average solution and evaluating the gradient at the current solution.
  # The number in the current implementation is four.
  KKT_PASSES_PER_TERMINATION_EVALUATION = 2.0

  p = initialize_saddle_point_problem(problem, params, combo_mat)
  primal_weight = select_initial_primal_weight(
    problem,
    primal_part(p.mirror_map_scaling),
    dual_part(p.mirror_map_scaling),
    params.primal_importance,
    params.verbosity,
  )
  update_mirror_prox_norms(p, primal_weight)

  primal_weight_update_smoothing =
    params.restart_params.primal_weight_update_smoothing

  iteration_stats = IterationStats[]
  start_time = time()
  # Basic algorithm refers to the primal and dual steps, and excludes restart
  # schemes and termination evaluation.
  time_spent_doing_basic_algorithm = 0.0

  cumulative_kkt_passes = 0.0

  # This variable is used in the adaptive restart scheme.
  last_restart_info = create_last_restart_info(
    problem,
    primal_part(current_solution),
    dual_part(current_solution),
  ) # TODO: Push running average inside of this, and maybe a
  # counter on number of iterations since a primal weight update (this could be
  # used to give more fined grain control on primal weight update frequency).

  # For termination criteria:
  termination_criteria = params.termination_criteria
  iteration_limit = params.termination_criteria.iteration_limit
  termination_evaluation_frequency = params.termination_evaluation_frequency

  # This flag represents whether a numerical error occurred during the algorithm
  # if it is set to true it will trigger the algorithm to terminate.
  numerical_error = false

  display_iteration_stats_heading(params.verbosity)

  iteration = 0
  while true
    iteration += 1

    current_gradient = pseudo_gradient(current_solution, p)

    # Evaluate the iteration stats at frequency
    # termination_evaluation_frequency, when the iteration_limit is reached,
    # or if a numerical error occurs at the previous iteration.
    if mod(iteration - 1, termination_evaluation_frequency) == 0 ||
       iteration == iteration_limit + 1 ||
       iteration <= 10 ||
       numerical_error
      cumulative_kkt_passes += KKT_PASSES_PER_TERMINATION_EVALUATION

      # Compute the average solution since the last restart point.
      if numerical_error || solution_weighted_avg.sum_solutions_count == 0
        avg_primal_solution = primal_part(current_solution)
        avg_dual_solution = dual_part(current_solution)
      else
        avg_primal_solution, avg_dual_solution =
          compute_average(solution_weighted_avg)
      end
      avg_solution = make_sp_vec(avg_primal_solution, avg_dual_solution)

      # Evalute the iteration stats.
      current_iteration_stats = evaluate_unscaled_iteration_stats(
        scaled_problem,
        qp_cache,
        params.termination_criteria,
        params.record_iteration_stats,
        avg_primal_solution,
        avg_dual_solution,
        iteration,
        time() - start_time,
        cumulative_kkt_passes,
        termination_criteria.eps_optimal_absolute,
        termination_criteria.eps_optimal_relative,
        step_size,
        p.primal_weight,
        POINT_TYPE_AVERAGE_ITERATE,
      )
      method_specific_stats = current_iteration_stats.method_specific_stats
      method_specific_stats["time_spent_doing_basic_algorithm"] =
        time_spent_doing_basic_algorithm
      update_objective_bound_estimates(
        current_iteration_stats.method_specific_stats,
        problem,
        avg_primal_solution,
        avg_dual_solution,
        primal_part(p.mirror_map_scaling),
        dual_part(p.mirror_map_scaling),
      )
      if params.record_iteration_stats
        push!(iteration_stats, current_iteration_stats)
      end

      # Check the termination criteria.
      termination_reason = check_termination_criteria(
        termination_criteria,
        qp_cache,
        current_iteration_stats,
      )
      if numerical_error && termination_reason == false
        termination_reason = TERMINATION_REASON_NUMERICAL_ERROR
      end

      mirror_prox_display(
        termination_reason,
        iteration,
        params,
        current_iteration_stats,
        p,
        accepted_iterations,
        acceptable_nonlinearity,
        step_size,
        inv_stepsize_required,
        avg_solution,
      )

      if termination_reason != false
        # ** Terminate the algorithm **
        # This is the only place the algorithm can terminate. Please keep it
        # this way.
        generic_final_log(
          problem,
          primal_part(avg_solution),
          dual_part(avg_solution),
          current_iteration_stats,
          params.verbosity,
          iteration,
          termination_reason,
        )
        return unscaled_saddle_point_output(
          scaled_problem,
          avg_primal_solution,
          avg_dual_solution,
          termination_reason,
          iterations_completed,
          iteration_stats,
        )
      end

      current_iteration_stats.restart_used = run_restart_scheme(
        problem,
        solution_weighted_avg,
        primal_part(current_solution),
        dual_part(current_solution),
        last_restart_info,
        iterations_completed,
        primal_part(p.mirror_map_scaling),
        dual_part(p.mirror_map_scaling),
        primal_weight,
        params.verbosity,
        params.restart_params,
      )

      # Update primal_weight
      if current_iteration_stats.restart_used != RESTART_CHOICE_NO_RESTART
        primal_weight = compute_new_primal_weight(
          last_restart_info,
          primal_weight,
          primal_weight_update_smoothing,
          params.verbosity,
        )
        update_mirror_prox_norms(p, primal_weight)
      end
    end

    time_spent_doing_basic_algorithm_checkpoint = time()
    current_gradient = pseudo_gradient(current_solution, p)
    # This line can be changed arbitrarily. The analysis will still work except
    # that if the guesses aren't good the steps may be rejected even with a
    # small step size. See Lemma 1 in
    # http://papers.nips.cc/paper/5147-optimization-learning-and-games-with-predictable-sequences.pdf
    guess_gradient = current_gradient
    test_point = proximal_step(p, current_solution, guess_gradient, step_size)
    test_gradient = pseudo_gradient(test_point, p)
    candidate_solution =
      proximal_step(p, current_solution, test_gradient, step_size)

    nonlinearity =
      dot(test_gradient - guess_gradient, test_point - candidate_solution)
    movement =
      bregman(candidate_solution, test_point, p) +
      bregman(test_point, current_solution, p)
    # If we didn't move we must have the problem solved to working precision.
    # Stop now to avoid problems such as divide by zero.
    if movement == 0
      # The algorithm will terminate at the beginning of the next iteration
      numerical_error = true
      continue
    end
    # This must be <= 1/step_size for the mirror prox proof to work.
    # It is at most the max singular value of the constraint matrix.
    inv_stepsize_required = nonlinearity / movement
    acceptable_nonlinearity::Bool = inv_stepsize_required <= 1.0 / step_size

    if acceptable_nonlinearity
      current_solution = candidate_solution

      # Use weighted average, if this was 1.0 it would be standard average.
      weight = step_size
      add_to_solution_weighted_average(
        solution_weighted_avg,
        test_point,
        weight,
      )
      accepted_iterations += 1
    end
    exponent_one = 0.3
    exponent_two = 0.6
    # Our step sizes are a factor 1 - iteration^(-exponent_one) smaller than
    # they could be as margin to reduce rejected steps.
    first_term =
      (1 - (iteration + 1)^(-exponent_one)) / abs(inv_stepsize_required)
    second_term = (1 + (iteration + 1)^(-exponent_two)) * step_size
    # From the first term when we have to reject a step the step_size decreases
    # by a factor of at least 1 - iteration^(-exponent_one).
    # From the second term we increase the step_size by a factor of at most
    # 1 + iteration^(-exponent_two).
    # Therefore if more than order iteration^(exponent_one - exponent_two)
    # fraction of the iterations have a rejected step we overall decrease the
    # step_size. When the step_size is below the inverse of the max singular
    # value we stop having rejected steps.
    step_size = min(first_term, second_term)
    iterations_completed += 1

    time_spent_doing_basic_algorithm +=
      time() - time_spent_doing_basic_algorithm_checkpoint

    cumulative_kkt_passes += KKT_PASSES_PER_ITERATION
  end
end


text:
# pkg> activate .
using BKMaxflow

# construct a n x n image
n = 128
# preallocate
g = create_graph(CImpl{Cdouble}, 17000, 200000)
weights, neighbors = create_graph(JuliaImpl{Cdouble,Int}, n*n+2)  # all pixels + 2 terminals

# data term
imageDims = (n,n)
pixelRange = CartesianIndices(imageDims)
pixelFirst, pixelEnd = first(pixelRange), last(pixelRange)
linearIdx = LinearIndices(imageDims)
for ii in pixelRange
    i = linearIdx[ii]
    source = rand()
    sink = rand()
    # julia
    add_edge!(weights, neighbors, 1, i+2, source, 0.0)
    add_edge!(weights, neighbors, 2, i+2, 0.0, sink)
    # c wrapper
    add_node(CImpl{Cdouble}, g)
    add_tweights(CImpl, g, i-1, source, sink)
end

# smooth term
costs = rand(200000)
idx = 0
for ii in pixelRange
    i = linearIdx[ii]
    neighborFirst = max(pixelFirst, ii-pixelFirst)
    neighborEnd = min(pixelEnd, ii+pixelFirst)
    neighborRange = CartesianIndices((neighborFirst[1]:neighborEnd[1], neighborFirst[2]:neighborEnd[2]))
    for jj in neighborRange
        if ii < jj
            j = linearIdx[jj]
            idx += 2
            # julia
            add_edge!(weights, neighbors, i+2, j+2, costs[idx-1], costs[idx])
            # c wrapper
            add_edge(CImpl, g, i-1, j-1, costs[idx-1], costs[idx])
        end
    end
end

flow, label = boykov_kolmogorov(1, 2, neighbors, weights)

flow = maxflow(CImpl{Cdouble}, g)

delete_graph(CImpl{Cdouble}, g)

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
# ---
# title: 1186. Maximum Subarray Sum with One Deletion
# id: problem1186
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/maximum-subarray-sum-with-one-deletion/description/>
# hidden: true
# ---
#
# Given an array of integers, return the maximum sum for a **non-empty**
# subarray (contiguous elements) with at most one element deletion. In other
# words, you want to choose a subarray and optionally delete one element from it
# so that there is still at least one element left and the sum of the remaining
# elements is maximum possible.
#
# Note that the subarray needs to be **non-empty** after deleting one element.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,-2,0,3]
#     Output: 4
#     Explanation: Because we can choose [1, -2, 0, 3] and drop -2, thus the subarray [1, 0, 3] becomes the maximum value.
#
# **Example 2:**
#
#
#
#     Input: arr = [1,-2,-2,3]
#     Output: 3
#     Explanation: We just choose [3] and it's the maximum sum.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [-1,-1,-1,-1]
#     Output: -1
#     Explanation:  The final subarray needs to be non-empty. You can't choose [-1] and delete -1 from it, then get an empty subarray to make the sum equals to 0.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 105`
#   * `-104 <= arr[i] <= 104`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function is_orthogonal(S, tol=1e-6)
    n = LinearAlgebra.checksquare(S)
    for i in 1:n
        for j in 1:n
            if i != j
                s = LinearAlgebra.dot(S[:, i], S[:, j])
                if abs(s) > tol
                    return false
                end
            end
        end
    end
    return true
end

function _projcoef(A, i, v, w = A[i, :])
    # `LinearAlgebra.dot(v, w)` would not work as it is computing the adjoint
    # of `v`.
    # If `v = α * w + ...` then
    # `dot(w, v) = conj(w) (α * w + ...) = α`
    # However,
    # `dot(v, w) = conj(α * w + ...) w = conj(α)`
    return LinearAlgebra.dot(w, v) / LinearAlgebra.dot(w, w)
end

struct _OrthogonalMatrix end
struct _RowEchelonMatrix end

function __linsolve(A::AbstractMatrix{T}, b::Vector{T}, ::_OrthogonalMatrix) where {T}
    return map(1:size(A, 1)) do i
        return _projcoef(A, i, b)
    end
end

function __linsolve(A::Matrix{T}, b::Vector{T}, ::_RowEchelonMatrix) where {T}
    j = 0
    return map(1:size(A, 1)) do i
        while j < size(A, 2)
            j += 1
            if isone(A[i, j]) && all(k -> k == i || iszero(A[k, j]), 1:size(A, 1))
                return b[j]
            end
        end
        error("Not in row_echelon_form, cannot find for `$i`th entry.")
    end
end

using SparseArrays
function _linsolve(A::AbstractMatrix{T}, b::Vector{T}, form) where {T}
    x = __linsolve(A, b, form)
    @assert transpose(A) * x ≈ b
    return x
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:light,ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.1
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # DFT: The LDA kernel
# ## I. Theory
#
# Previously we described the DFT Fock matrix as
# $$F^{DFT}_{\mu\nu} = H_{\mu\nu} + 2J[D]_{\mu\nu} - \zeta K[D]_{\mu\nu} + V^{\rm{xc}}_{\mu\nu}$$
# upon examination it is revealed that the only quantity that we cannot yet compute is $V^{\rm{xc}}$.
#
# Here we will explore the local density approximation (LDA) functionals where $V^{\rm{xc}} = f[\rho(\hat{r})]$. For these functionals the only required bit of information is the density at the grid point. As we discussed the grid last chapter we will now focus exactly on how to obtain the density on the grid.
#
# Before we begin we should first recall that the Fock matrix is the derivative of the energy with respect to atomic orbitals. Therefore, the $V^{\rm{xc}}$ matrix is not the XC energy, but the derivate of that energy, which can expressed as $\frac{\partial e_{\rm{xc}}}{\partial\rho}$.

# +
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Julia Arrays

mol = psi4.geometry("""
He
symmetry c1
""")
psi4.set_options(Dict("BASIS" => "cc-pvdz",
                      "DFT_SPHERICAL_POINTS" => 6,
                      "DFT_RADIAL_POINTS" => 5))

psi4.core.set_output_file("output.dat", false)

svwn_w, wfn = psi4.energy("SVWN", return_wfn=true)
Vpot = wfn.V_potential()
# -

# ## 2. Density on a Grid
# The density on the grid can be expressed as
# $$\rho(\hat{r}) = \sum\limits_{\mu\nu} D_{\mu\nu}\;\phi_\mu(\hat{r})\phi_\nu(\hat{r})$$
#
# Recall that we compute DFT quanties on a grid, so $\hat{r}$ will run over a grid instead of all space. Using this we can build collocation matrices that map between atomic orbital and grid space $$\phi_\mu(\hat{r}) \rightarrow \phi_\mu^p$$
# where our $p$ index will be the index of individual grid points. Our full expression becomes:
#
# $$\rho_p = \phi_\mu^p D_{\mu\nu} \phi_\nu^p$$
#
# To compute these quantities let us first remember that the DFT grid is blocked loosely over atoms. It should now be apparent to why we do this, consider the $\phi_\mu^p$ objects. The total size of this object would be `nbf` $\times$ `npoints`. To put this in perspective a moderate size molecule could have 1e4 basis functions and 1e8 grid points, so about 8 terabytes of data! As this object is very sparse it is much more convenient to store the grid and compute $\phi\mu^p$ matrices on the fly.
#
# We then need object to compute $\phi_\mu^p$.

# +
# Grab a "points function" to compute the Phi matrices
points_func = Vpot.properties()[1]

# Grab a block and obtain its local mapping
block = Vpot.get_block(1)
npoints = block.npoints()
lpos = np.array(block.functions_local_to_global()) .+ 1
println("Local basis function mapping")
display(lpos)

# Copmute phi, note the number of points and function per phi changes.
phi = np.array(points_func.basis_values()["PHI"])[1:npoints, 1:size(lpos,1)]
println("\nPhi Matrix")
display(phi)
# -

# ## 3. Evaluating the kernel
#
# After building the density on the grid we can then compute the exchange-correlation $f_{xc}$ at every gridpoint. This then need to be reintegrated back to atomic orbital space which can be accomplished like so:
#
# $$V^{\rm{xc}}_{pq}[D_{pq}] = \phi_\mu^a\;\phi_\nu^a\;\; w^a\;f^a_{\rm{xc}}{(\phi_\mu^p D_{\mu\nu} \phi_\nu^p)}$$
#
# Where $w^a$ is our combined Truetler and Lebedev weight at every point.
#
# Unlike SCF theory where the SCF energy can be computed as the sum of the Fock and Density matrices the energy for XC kernels must be computed in grid space. Fortunately, the energy is simply defined as:
#
# $$e_{\rm{xc}} = w^a f^a_{\rm{xc}}$$
#
# We can now put all the pieces together to compute $e_{\rm{xc}}$ and $\frac{\partial E_{\rm{xc}}}{\partial\rho}= V^{\rm{xc}}$.

# +
D = np.array(wfn.Da())

V = zero(D)

rho = []
points_func = Vpot.properties()[1]
superfunc = Vpot.functional()

xc_e, V = let xc_e = 0.0, rho = rho
   # Loop over the blocks
   for b in 1:Vpot.nblocks()

       # Obtain block information
       block = Vpot.get_block(b-1)
       points_func.compute_points(block)
       npoints = block.npoints()
       lpos = np.array(block.functions_local_to_global()) .+ 1

       # Obtain the grid weight
       w = np.array(block.w())

       # Compute ϕᴾμ!
       phi = np.array(points_func.basis_values()["PHI"])[1:npoints, 1:size(lpos,1)]

       # Build a local slice of D
       lD = D[lpos,lpos]

       # Copmute ρ
       # ρᴾ = ϕᴾμ Dμν ϕᴾν
       rho = 2vec(sum( (phi * lD) .* phi, dims=2))

       inp = Dict()
       inp["RHO_A"] = psi4.core.Vector.from_array(rho)

       # Compute the kernel
       ret = superfunc.compute_functional(inp, -1)

       # Compute the XC energy
       vk = np.array(ret["V"])[1:npoints]
       xc_e = w .* vk

       # Compute the XC derivative.
       v_rho_a = np.array(ret["V_RHO_A"])[1:npoints]
       # Vab = ϕᴾ_b Vᴾ Wᴾ ϕᴾ_a
       Vtmp = phi' * (v_rho_a .* w .* phi)

       # Add the temporary back to the larger array by indexing, ensure it is symmetric
       V[lpos, lpos] += 0.5(Vtmp + Vtmp')
   end
   xc_e, V
end

println("XC Energy: ")
display(xc_e)
println("V matrix:")
display(V)

println("\nMatches Psi4 V: ", np.allclose(V, wfn.Va()))
# -

# Refs:
# - Johnson, B. G.; Fisch M. J.; *J. Chem. Phys.*, **1994**, *100*, 7429


text:
@testset "Knapsack" begin

    @testset "Knapsack: DpZeroOnePack!" begin
        dp = zeros(Int, 30)
        @test dp_zero_one_pack!(20, [1,3,11], [2,5,30], dp) == 37
        dp = zeros(Int, 30)
        @test dp_zero_one_pack!(10, [1,3,11], [20,5,80], dp) == 25
    end

    @testset "Knapsack: DpCompletePack!" begin
        dp = zeros(Int, 30)
        @test dp_complete_pack!(20, [1,2,9], [1,3,20], dp) == 43
        dp = zeros(Int, 30)
        @test dp_complete_pack!(10, [1,3,11], [20,5,80], dp) == 200
    end

    @testset "Knapsack: GaZeroOnePack!" begin
        dp = zeros(Int, 30)
        @test ga_zero_one_pack!(20, [1,3,11], [2,5,30]) == 37
        dp = zeros(Int, 30)
        @test ga_zero_one_pack!(10, [1,3,11], [20,5,80]) == 25
    end

    @testset "Knapsack: GaCompletePack!" begin
        dp = zeros(Int, 30)
        @test ga_complete_pack!(20, [1,2,9], [1,3,20]) == 43
        dp = zeros(Int, 30)
        @test ga_complete_pack!(10, [1,3,11], [20,5,80]) == 200
    end

end


text:
export TabularCFRPolicy

struct InfoStateNode
    strategy::Vector{Float64}
    cumulative_regret::Vector{Float64}
    cumulative_strategy::Vector{Float64}
end

InfoStateNode(n) = InfoStateNode(fill(1 / n, n), zeros(n), zeros(n))

function init_info_state_nodes(env::AbstractEnv)
    nodes = Dict{String,InfoStateNode}()
    walk(env) do x
        if !get_terminal(x) && get_current_player(x) != get_chance_player(x)
            get!(nodes, get_state(x), InfoStateNode(length(get_legal_actions(x))))
        end
    end
    nodes
end

"""
    TabularCFRPolicy

See more details: [An Introduction to Counterfactual Regret Minimization](http://modelai.gettysburg.edu/2013/cfr/cfr.pdf)
"""
struct TabularCFRPolicy{S,T,R<:AbstractRNG} <: AbstractPolicy
    nodes::Dict{S,InfoStateNode}
    behavior_policy::QBasedPolicy{TabularLearner{S,T},WeightedExplorer{true,R}}
end

(p::TabularCFRPolicy)(env::AbstractEnv) = p.behavior_policy(env)

RLBase.get_prob(p::TabularCFRPolicy, env::AbstractEnv) = get_prob(p.behavior_policy, env)

"""
    TabularCFRPolicy(;n_iter::Int, env::AbstractEnv)
"""
function TabularCFRPolicy(;
    n_iter::Int,
    env::AbstractEnv,
    rng = Random.GLOBAL_RNG,
    is_reset_neg_regrets = false,
    is_linear_averaging = false,
)
    @assert NumAgentStyle(env) isa MultiAgent
    @assert DynamicStyle(env) === SEQUENTIAL
    @assert RewardStyle(env) === TERMINAL_REWARD
    @assert ChanceStyle(env) === EXPLICIT_STOCHASTIC
    @assert DefaultStateStyle(env) === Information{String}()

    nodes = init_info_state_nodes(env)

    for i in 1:n_iter
        for p in get_players(env)
            if p != get_chance_player(env)
                init_reach_prob =
                    Dict(x => 1.0 for x in get_players(env) if x != get_chance_player(env))
                cfr!(nodes, env, p, init_reach_prob, 1.0, is_linear_averaging ? i : 1)
                update_strategy!(nodes)

                if is_reset_neg_regrets
                    for node in values(nodes)
                        node.cumulative_regret .= max.(node.cumulative_regret, 0)
                    end
                end
            end
        end
    end

    behavior_policy = QBasedPolicy(;
        learner = TabularLearner{String}(),
        explorer = WeightedExplorer(; is_normalized = true, rng = rng),
    )

    for (k, v) in nodes
        s = sum(v.cumulative_strategy)
        if s != 0
            update!(behavior_policy, k => v.cumulative_strategy ./ s)
        end
    end

    TabularCFRPolicy(nodes, behavior_policy)
end

function cfr!(nodes, env, player, reach_probs, chance_player_reach_prob, ratio)
    if get_terminal(env)
        get_reward(env, player)
    else
        if get_current_player(env) == get_chance_player(env)
            v = 0.0
            for a::ActionProbPair in get_legal_actions(env)
                v +=
                    a.prob * cfr!(
                        nodes,
                        child(env, a),
                        player,
                        reach_probs,
                        chance_player_reach_prob * a.prob,
                        ratio,
                    )
            end
            v
        else
            v = 0.0
            node = nodes[get_state(env)]
            legal_actions = get_legal_actions(env)
            U =
                player == get_current_player(env) ?
                Vector{Float64}(undef, length(legal_actions)) : nothing

            for (i, action) in enumerate(legal_actions)
                prob = node.strategy[i]
                new_reach_probs = copy(reach_probs)
                new_reach_probs[get_current_player(env)] *= prob

                u = cfr!(
                    nodes,
                    child(env, action),
                    player,
                    new_reach_probs,
                    chance_player_reach_prob,
                    ratio,
                )
                isnothing(U) || (U[i] = u)
                v += prob * u
            end

            if player == get_current_player(env)
                reach_prob = reach_probs[player]
                counterfactual_reach_prob = reduce(
                    *,
                    (
                        reach_probs[p]
                        for
                        p in get_players(env) if p != player && p != get_chance_player(env)
                    );
                    init = chance_player_reach_prob,
                )
                node.cumulative_regret .+= counterfactual_reach_prob .* (U .- v)
                node.cumulative_strategy .+= ratio .* reach_prob .* node.strategy
            end
            v
        end
    end
end

regret_matching!(node::InfoStateNode) =
    regret_matching!(node.strategy, node.cumulative_regret)

function regret_matching!(strategy, cumulative_regret)
    s = mapreduce(x -> max(0, x), +, cumulative_regret)
    if s > 0
        strategy .= max.(0.0, cumulative_regret) ./ s
    else
        fill!(strategy, 1 / length(strategy))
    end
end

function update_strategy!(nodes)
    for node in values(nodes)
        regret_matching!(node)
    end
end


text:
include("Graph.jl")
include("NodeCentrality.jl")
include("EdgeCentrality.jl")

function NormalizedBetweennessCentrality(G)
    NC = ((G.n > 10000) || (G.m > 10000)) ? ApproxBetweennessCentrality(G) : BetweennessCentrality(G)
    EC = ((G.n > 10000) || (G.m > 10000)) ? ApproxEdgeBetweenness(G) : EdgeBetweenness(G)
    NC ./= sum(NC)
    EC ./= sum(EC)
    return NC, EC
end

function selectNodes(G, selectedNodes, selectedEdges, selectMode)
    NC, EC = NormalizedBetweennessCentrality(G)
    nl = zeros(Int, G.n)
    el = zeros(Int, G.m)
    foreach(i -> nl[i] = i, 1 : G.n)
    foreach(i -> el[i] = i, 1 : G.m)
    sort!(nl, by = x->-NC[x])
    sort!(el, by = x->-EC[x])

    Selected = zeros(Bool, G.n)
    if ('n' in selectMode)
        foreach(i -> Selected[nl[i]] = true, 1 : min(selectedNodes, G.n))
    end
    if ('e' in selectMode)
        for i = 1 : min(selectedEdges, G.m)
            Selected[G.E[el[i]][2]] = true
            Selected[G.E[el[i]][3]] = true
        end
    end

    return NC, EC, Selected
end

function SearchWay(g, status)
    n = size(g, 1)
    h = copy(status)
    Q = zeros(Int, n)
    front = 1
    rear = 0

    pushQ(x) = begin
        rear += 1
        Q[rear] = x
    end
    popQ() = begin
        tmp = Q[front]
        front += 1
        return tmp
    end

    for i = 1 : n
        (h[i] == 1) ? pushQ(i) : nothing
    end

    prior = zeros(Int, n)

    while front <= rear
        u = popQ()
        for v in g[u]
            if h[v] == 0
                h[v] = 1
                prior[v] = u
                pushQ(v)
            elseif h[v] == 2
                nodeList = Array{Int, 1}()
                push!(nodeList, v)
                push!(nodeList, u)
                p = prior[u]
                while p != 0
                    push!(nodeList, p)
                    p = prior[p]
                end
                return nodeList
            end
        end
    end

    return nothing
end

function KeepOne(G, Selected, NC, EC)
    g = getAdjacentList(G)
    F = zeros(Int, G.n)
    foreach(i -> F[i] = i, 1 : G.n)

    find(x) = begin
        if F[x] != x
            F[x] = find(F[x])
        end
        return F[x]
    end
    SetUnion(x, y) = F[find(x)] = find(y)

    for (ID, u, v, w) in G.E
        (Selected[u] && Selected[v]) ? SetUnion(u, v) : nothing
    end
    status = zeros(Int, G.n)

    while true
        clist = Array{Int, 1}()
        for i = 1 : G.n
            if Selected[i] && (F[i] == i)
                push!(clist, i)
            end
        end
        Bst = Array{Int, 1}()
        bstlen = G.n
        fill!(status, 0)
        for u in clist
            for i = 1 : G.n
                if Selected[i]
                    status[i] = (find(i) == u) ? 1 : 2
                end
            end
            tmp = SearchWay(g, status)
            if (tmp != nothing) && (size(tmp, 1) < bstlen)
                bstlen = size(tmp, 1)
                Bst = copy(tmp)
            end
        end
        if bstlen == G.n
            break
        end
        for i = 1 : bstlen
            Selected[Bst[i]] = true
            (i > 1) ? SetUnion(Bst[i], Bst[i - 1]) : nothing
        end
    end

    NodeC = Dict{Int, Float64}()
    NCmax = NC[argmax(NC)]
    for i = 1 : G.n
        if Selected[i]
            NodeC[G.V[i]] = NC[i] / NCmax
        end
    end

    return SubGraph(G, Selected, EC), NodeC, NCmax
end

function SearchAll(g, st, selected)
    n = size(g, 1)
    visited = zeros(Bool, n)
    visited[st] = true
    Q = zeros(Int, n)
    front = 1
    rear = 0

    pushQ(x) = begin
        rear += 1
        Q[rear] = x
    end
    popQ() = begin
        tmp = Q[front]
        front += 1
        return tmp
    end

    pushQ(st)

    prior = zeros(Int, n)

    while front <= rear
        u = popQ()
        for v in g[u]
            if !visited[v]
                visited[v] = true
                prior[v] = u
                pushQ(v)
            end
        end
    end

    needNodes = zeros(Bool, n)
    for v = 1 : n
        if selected[v]
            x = v
            while (x > 0)
                needNodes[x] = true
                x = prior[x]
            end
        end
    end

    return needNodes
end

function KeepAll(G, Selected, NC, EC)
    g = getAdjacentList(G)
    Needed = copy(Selected)
    for i = 1 : G.n
        if Selected[i]
            needNodes = SearchAll(g, i, Selected)
            foreach(j -> Needed[j] = Needed[j] || needNodes[j], 1 : G.n)
        end
    end

    NodeC = Dict{Int, Float64}()
    NCmax = NC[argmax(NC)]
    for i = 1 : G.n
        if Needed[i]
            NodeC[G.V[i]] = NC[i] / NCmax
        end
    end

    return SubGraph(G, Needed, EC), NodeC, NCmax
end

function Compress(G, selectedNodes, selectedEdges, selectMode, addMode)
    NC, EC, Selected = selectNodes(G, selectedNodes, selectedEdges, selectMode)
    return (addMode == "KeepOne") ? KeepOne(G, Selected, NC, EC) : KeepAll(G, Selected, NC, EC)
end

function CentralityInformation(G, NodeC, maxNodeC)
    nci = 0.0
    eci = 0.0
    foreach(i -> nci += (NodeC[G.V[i]]*maxNodeC), 1 : G.n)
    foreach(i -> eci += G.E[i][4], 1 : G.m)
    return nci, eci, (nci+eci)/2
end


text:
# Major changes from original:
#   range(A) -> range(0,length=A) range syntax
#   (lat(i,j) for i ... for j...) -> [lat(i,j) for i ... for j...], comprehension syntax
#   pyplot.show() -> pyplot.gcf() or import PyPlot and use that
#   onsite potential function must be defined as pythong function using py"""..."""
#       with L, L_well, t interpolated, and passed as py"...",
#   A.append(...) -> append!(A,...)

function make_system(a=1, t=1.0, W=10, L=30, L_well=10)
    # Start with an empty tight-binding system and a single square lattice.
    # `a` is the lattice constant (by default set to 1 for simplicity.
    lat = kwant.lattice.square(a)

    syst = kwant.Builder()

    #### Define the scattering region. ####
    # Potential profile
    # note that passing functions is tricky -- haven't yet figured out how to do it
    # easiest way is to define python functions directly and pass them
    py"""
    def potential(site, pot):
        (x, y) = site.pos
        if ($L - $L_well) / 2 < x < ($L + $L_well) / 2:
            return pot
        else:
            return 0

    def onsite(site, pot):
            return 4 * $t + potential(site, pot)
    """

    syst[[lat(x, y) for x in range(0,length=L) for y in range(0,length=W)]] = py"onsite"
    syst[lat.neighbors()] = -t

    #### Define and attach the leads. ####
    lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
    lead[[lat(1, j) for j in range(0,length=W)]] = 4 * t ## note [lat...] not (lat...) for comprehension
    lead[lat.neighbors()] = -t
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    return syst
end

function plot_conductance(syst, energy, welldepths)

    # Compute conductance
    data = []
    for welldepth in welldepths
        smatrix = kwant.smatrix(syst, energy, params=Dict(:pot=>-welldepth)) # note Dict(:pot=>) not dict(pot=)
        append!(data,smatrix.transmission(1, 0)) ## note append! not data.append
    end

    pyplot.figure()
    pyplot.plot(welldepths, data)
    pyplot.xlabel("well depth [t]")
    pyplot.ylabel("conductance [e^2/h]")
    pyplot.gcf() ## note gcf not show
end