text:
function get_generator_data(data::OPFData, unknown)
    error("Unknown device $typeof(unknown)")
end

function get_generator_data(data::OPFData, ::KA.CPU)
    ngen = length(data.generators)

    pgmin = Array{Float64}(undef, ngen)
    pgmax = Array{Float64}(undef, ngen)
    qgmin = Array{Float64}(undef, ngen)
    qgmax = Array{Float64}(undef, ngen)
    c2 = Array{Float64}(undef, ngen)
    c1 = Array{Float64}(undef, ngen)
    c0 = Array{Float64}(undef, ngen)

    Pmin = Float64[data.generators[g].Pmin for g in 1:ngen]
    Pmax = Float64[data.generators[g].Pmax for g in 1:ngen]
    Qmin = Float64[data.generators[g].Qmin for g in 1:ngen]
    Qmax = Float64[data.generators[g].Qmax for g in 1:ngen]
    coeff0 = Float64[data.generators[g].coeff[3] for g in 1:ngen]
    coeff1 = Float64[data.generators[g].coeff[2] for g in 1:ngen]
    coeff2 = Float64[data.generators[g].coeff[1] for g in 1:ngen]
    copyto!(pgmin, Pmin)
    copyto!(pgmax, Pmax)
    copyto!(qgmin, Qmin)
    copyto!(qgmax, Qmax)
    copyto!(c0, coeff0)
    copyto!(c1, coeff1)
    copyto!(c2, coeff2)

    return pgmin,pgmax,qgmin,qgmax,c2,c1,c0
end

function get_generator_data(data::OPFData, ::CUDADevice)
    ngen = length(data.generators)

    pgmin = CuArray{Float64}(undef, ngen)
    pgmax = CuArray{Float64}(undef, ngen)
    qgmin = CuArray{Float64}(undef, ngen)
    qgmax = CuArray{Float64}(undef, ngen)
    c2 = CuArray{Float64}(undef, ngen)
    c1 = CuArray{Float64}(undef, ngen)
    c0 = CuArray{Float64}(undef, ngen)

    Pmin = Float64[data.generators[g].Pmin for g in 1:ngen]
    Pmax = Float64[data.generators[g].Pmax for g in 1:ngen]
    Qmin = Float64[data.generators[g].Qmin for g in 1:ngen]
    Qmax = Float64[data.generators[g].Qmax for g in 1:ngen]
    coeff0 = Float64[data.generators[g].coeff[3] for g in 1:ngen]
    coeff1 = Float64[data.generators[g].coeff[2] for g in 1:ngen]
    coeff2 = Float64[data.generators[g].coeff[1] for g in 1:ngen]
    copyto!(pgmin, Pmin)
    copyto!(pgmax, Pmax)
    copyto!(qgmin, Qmin)
    copyto!(qgmax, Qmax)
    copyto!(c0, coeff0)
    copyto!(c1, coeff1)
    copyto!(c2, coeff2)

    return pgmin,pgmax,qgmin,qgmax,c2,c1,c0
end

function get_generator_data(data::OPFData, ::ROCDevice)
    ngen = length(data.generators)

    pgmin = ROCArray{Float64}(undef, ngen)
    pgmax = ROCArray{Float64}(undef, ngen)
    qgmin = ROCArray{Float64}(undef, ngen)
    qgmax = ROCArray{Float64}(undef, ngen)
    c2 = ROCArray{Float64}(undef, ngen)
    c1 = ROCArray{Float64}(undef, ngen)
    c0 = ROCArray{Float64}(undef, ngen)

    Pmin = Float64[data.generators[g].Pmin for g in 1:ngen]
    Pmax = Float64[data.generators[g].Pmax for g in 1:ngen]
    Qmin = Float64[data.generators[g].Qmin for g in 1:ngen]
    Qmax = Float64[data.generators[g].Qmax for g in 1:ngen]
    coeff0 = Float64[data.generators[g].coeff[3] for g in 1:ngen]
    coeff1 = Float64[data.generators[g].coeff[2] for g in 1:ngen]
    coeff2 = Float64[data.generators[g].coeff[1] for g in 1:ngen]
    copyto!(pgmin, Pmin)
    copyto!(pgmax, Pmax)
    copyto!(qgmin, Qmin)
    copyto!(qgmax, Qmax)
    copyto!(c0, coeff0)
    copyto!(c1, coeff1)
    copyto!(c2, coeff2)

    return pgmin,pgmax,qgmin,qgmax,c2,c1,c0
end

function get_bus_data(data::OPFData, unknown)
    error("Unknown device $typeof(unknown)")
end

function get_bus_data(data::OPFData, ::KA.CPU)
    nbus = length(data.buses)

    FrIdx = [l for b=1:nbus for l in data.FromLines[b]]
    ToIdx = [l for b=1:nbus for l in data.ToLines[b]]
    GenIdx = [g for b=1:nbus for g in data.BusGenerators[b]]
    FrStart = accumulate(+, vcat([1], [length(data.FromLines[b]) for b=1:nbus]))
    ToStart = accumulate(+, vcat([1], [length(data.ToLines[b]) for b=1:nbus]))
    GenStart = accumulate(+, vcat([1], [length(data.BusGenerators[b]) for b=1:nbus]))

    Pd = Float64[data.buses[i].Pd for i=1:nbus]
    Qd = Float64[data.buses[i].Qd for i=1:nbus]

    return FrStart,FrIdx,ToStart,ToIdx,GenStart,GenIdx,Pd,Qd
end

function get_bus_data(data::OPFData, ::CUDADevice)
    nbus = length(data.buses)

    FrIdx = [l for b=1:nbus for l in data.FromLines[b]]
    ToIdx = [l for b=1:nbus for l in data.ToLines[b]]
    GenIdx = [g for b=1:nbus for g in data.BusGenerators[b]]
    FrStart = accumulate(+, vcat([1], [length(data.FromLines[b]) for b=1:nbus]))
    ToStart = accumulate(+, vcat([1], [length(data.ToLines[b]) for b=1:nbus]))
    GenStart = accumulate(+, vcat([1], [length(data.BusGenerators[b]) for b=1:nbus]))

    Pd = Float64[data.buses[i].Pd for i=1:nbus]
    Qd = Float64[data.buses[i].Qd for i=1:nbus]

    cuFrIdx = CuArray{Int}(undef, length(FrIdx))
    cuToIdx = CuArray{Int}(undef, length(ToIdx))
    cuGenIdx = CuArray{Int}(undef, length(GenIdx))
    cuFrStart = CuArray{Int}(undef, length(FrStart))
    cuToStart = CuArray{Int}(undef, length(ToStart))
    cuGenStart = CuArray{Int}(undef, length(GenStart))
    cuPd = CuArray{Float64}(undef, nbus)
    cuQd = CuArray{Float64}(undef, nbus)

    copyto!(cuFrIdx, FrIdx)
    copyto!(cuToIdx, ToIdx)
    copyto!(cuGenIdx, GenIdx)
    copyto!(cuFrStart, FrStart)
    copyto!(cuToStart, ToStart)
    copyto!(cuGenStart, GenStart)
    copyto!(cuPd, Pd)
    copyto!(cuQd, Qd)

    return cuFrStart,cuFrIdx,cuToStart,cuToIdx,cuGenStart,cuGenIdx,cuPd,cuQd
end

function get_bus_data(data::OPFData, ::ROCDevice)
    nbus = length(data.buses)

    FrIdx = [l for b=1:nbus for l in data.FromLines[b]]
    ToIdx = [l for b=1:nbus for l in data.ToLines[b]]
    GenIdx = [g for b=1:nbus for g in data.BusGenerators[b]]
    FrStart = accumulate(+, vcat([1], [length(data.FromLines[b]) for b=1:nbus]))
    ToStart = accumulate(+, vcat([1], [length(data.ToLines[b]) for b=1:nbus]))
    GenStart = accumulate(+, vcat([1], [length(data.BusGenerators[b]) for b=1:nbus]))

    Pd = Float64[data.buses[i].Pd for i=1:nbus]
    Qd = Float64[data.buses[i].Qd for i=1:nbus]

    cuFrIdx = ROCArray{Int}(undef, length(FrIdx))
    cuToIdx = ROCArray{Int}(undef, length(ToIdx))
    cuGenIdx = ROCArray{Int}(undef, length(GenIdx))
    cuFrStart = ROCArray{Int}(undef, length(FrStart))
    cuToStart = ROCArray{Int}(undef, length(ToStart))
    cuGenStart = ROCArray{Int}(undef, length(GenStart))
    cuPd = ROCArray{Float64}(undef, nbus)
    cuQd = ROCArray{Float64}(undef, nbus)

    copyto!(cuFrIdx, FrIdx)
    copyto!(cuToIdx, ToIdx)
    copyto!(cuGenIdx, GenIdx)
    copyto!(cuFrStart, FrStart)
    copyto!(cuToStart, ToStart)
    copyto!(cuGenStart, GenStart)
    copyto!(cuPd, Pd)
    copyto!(cuQd, Qd)

    return cuFrStart,cuFrIdx,cuToStart,cuToIdx,cuGenStart,cuGenIdx,cuPd,cuQd
end

function get_branch_data(data::OPFData, unknown)
    error("Unknown device $typeof(unknown)")
end

function get_branch_data(data::OPFData, device::KA.CPU)
    buses = data.buses
    lines = data.lines
    BusIdx = data.BusIdx
    nline = length(data.lines)
    ybus = Ybus{Array{Float64}}(computeAdmitances(data.lines, data.buses, data.baseMVA, device; VI=Array{Int}, VD=Array{Float64})...)
    frBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].from]].Vmin^2, buses[BusIdx[lines[l].from]].Vmax^2) ]
    toBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].to]].Vmin^2, buses[BusIdx[lines[l].to]].Vmax^2) ]

    return ybus.YshR, ybus.YshI, ybus.YffR, ybus.YffI, ybus.YftR, ybus.YftI,
            ybus.YttR, ybus.YttI, ybus.YtfR, ybus.YtfI, frBound, toBound
end

function get_branch_data(data::OPFData, device::CUDADevice)
    buses = data.buses
    lines = data.lines
    BusIdx = data.BusIdx
    nline = length(data.lines)
    ybus = Ybus{Array{Float64}}(computeAdmitances(data.lines, data.buses, data.baseMVA, device; VI=Array{Int}, VD=Array{Float64})...)
    frBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].from]].Vmin^2, buses[BusIdx[lines[l].from]].Vmax^2) ]
    toBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].to]].Vmin^2, buses[BusIdx[lines[l].to]].Vmax^2) ]

    cuYshR = CuArray{Float64}(undef, length(ybus.YshR))
    cuYshI = CuArray{Float64}(undef, length(ybus.YshI))
    cuYffR = CuArray{Float64}(undef, nline)
    cuYffI = CuArray{Float64}(undef, nline)
    cuYftR = CuArray{Float64}(undef, nline)
    cuYftI = CuArray{Float64}(undef, nline)
    cuYttR = CuArray{Float64}(undef, nline)
    cuYttI = CuArray{Float64}(undef, nline)
    cuYtfR = CuArray{Float64}(undef, nline)
    cuYtfI = CuArray{Float64}(undef, nline)
    cuFrBound = CuArray{Float64}(undef, 2*nline)
    cuToBound = CuArray{Float64}(undef, 2*nline)
    copyto!(cuYshR, ybus.YshR)
    copyto!(cuYshI, ybus.YshI)
    copyto!(cuYffR, ybus.YffR)
    copyto!(cuYffI, ybus.YffI)
    copyto!(cuYftR, ybus.YftR)
    copyto!(cuYftI, ybus.YftI)
    copyto!(cuYttR, ybus.YttR)
    copyto!(cuYttI, ybus.YttI)
    copyto!(cuYtfR, ybus.YtfR)
    copyto!(cuYtfI, ybus.YtfI)
    copyto!(cuFrBound, frBound)
    copyto!(cuToBound, toBound)

    return cuYshR, cuYshI, cuYffR, cuYffI, cuYftR, cuYftI,
            cuYttR, cuYttI, cuYtfR, cuYtfI, cuFrBound, cuToBound
end

function get_branch_data(data::OPFData, device::ROCDevice)
    buses = data.buses
    lines = data.lines
    BusIdx = data.BusIdx
    nline = length(data.lines)
    ybus = Ybus{Array{Float64}}(computeAdmitances(data.lines, data.buses, data.baseMVA, device; VI=Array{Int}, VD=Array{Float64})...)
    frBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].from]].Vmin^2, buses[BusIdx[lines[l].from]].Vmax^2) ]
    toBound = [ x for l=1:nline for x in (buses[BusIdx[lines[l].to]].Vmin^2, buses[BusIdx[lines[l].to]].Vmax^2) ]

    cuYshR = ROCArray{Float64}(undef, length(ybus.YshR))
    cuYshI = ROCArray{Float64}(undef, length(ybus.YshI))
    cuYffR = ROCArray{Float64}(undef, nline)
    cuYffI = ROCArray{Float64}(undef, nline)
    cuYftR = ROCArray{Float64}(undef, nline)
    cuYftI = ROCArray{Float64}(undef, nline)
    cuYttR = ROCArray{Float64}(undef, nline)
    cuYttI = ROCArray{Float64}(undef, nline)
    cuYtfR = ROCArray{Float64}(undef, nline)
    cuYtfI = ROCArray{Float64}(undef, nline)
    cuFrBound = ROCArray{Float64}(undef, 2*nline)
    cuToBound = ROCArray{Float64}(undef, 2*nline)
    copyto!(cuYshR, ybus.YshR)
    copyto!(cuYshI, ybus.YshI)
    copyto!(cuYffR, ybus.YffR)
    copyto!(cuYffI, ybus.YffI)
    copyto!(cuYftR, ybus.YftR)
    copyto!(cuYftI, ybus.YftI)
    copyto!(cuYttR, ybus.YttR)
    copyto!(cuYttI, ybus.YttI)
    copyto!(cuYtfR, ybus.YtfR)
    copyto!(cuYtfI, ybus.YtfI)
    copyto!(cuFrBound, frBound)
    copyto!(cuToBound, toBound)

    return cuYshR, cuYshI, cuYffR, cuYffI, cuYftR, cuYftI,
            cuYttR, cuYttI, cuYtfR, cuYtfI, cuFrBound, cuToBound
end

function init_solution!(env::AdmmEnv, sol::SolutionOneLevel, ybus::Ybus, rho_pq, rho_va)
    data, model = env.data, env.model

    lines = data.lines
    buses = data.buses
    BusIdx = data.BusIdx
    ngen = length(data.generators)
    nline = length(data.lines)

    YffR = ybus.YffR; YffI = ybus.YffI
    YttR = ybus.YttR; YttI = ybus.YttI
    YftR = ybus.YftR; YftI = ybus.YftI
    YtfR = ybus.YtfR; YtfI = ybus.YtfI
    AT = typeof(sol.v_curr)
    for g=1:ngen
        pg_idx = model.gen_mod.gen_start + 2*(g-1)
        sol.v_curr[pg_idx:pg_idx+1] .= AT([
            0.5*(data.generators[g].Pmin + data.generators[g].Pmax),
            0.5*(data.generators[g].Qmin + data.generators[g].Qmax),
        ])
    end

    sol.rho .= rho_pq
    fill!(sol.u_curr, 0.0)
    fill!(sol.v_curr, 0.0)

    AT = typeof(sol.u_curr)
    for l=1:nline
        wij0 = (buses[BusIdx[lines[l].from]].Vmax^2 + buses[BusIdx[lines[l].from]].Vmin^2) / 2
        wji0 = (buses[BusIdx[lines[l].to]].Vmax^2 + buses[BusIdx[lines[l].to]].Vmin^2) / 2
        wR0 = sqrt(wij0 * wji0)

        pij_idx = model.line_start + 8*(l-1)

        sol.u_curr[pij_idx:pij_idx+3] .= AT([
            YffR[l] * wij0 + YftR[l] * wR0,
            -YffI[l] * wij0 - YftI[l] * wR0,
            YttR[l] * wji0 + YtfR[l] * wR0,
            -YttI[l] * wji0 - YtfI[l] * wR0,
        ])
        #=
        u_curr[pij_idx+4] = wij0
        u_curr[pij_idx+5] = wji0
        u_curr[pij_idx+6] = 0.0
        u_curr[pij_idx+7] = 0.0
        =#
        # wRIij[2*(l-1)+1] = wR0
        # wRIij[2*l] = 0.0
        sol.v_curr[pij_idx+4:pij_idx+7] .= AT([
            wij0,
            wji0,
            0.0,
            0.0,
        ])

        sol.rho[pij_idx+4:pij_idx+7] .= rho_va
    end

    sol.l_curr .= 0
    return
end

@kernel function copy_data_kernel(n, dest, @Const(src))
    I_ = @index(Group, Linear)
    J_ = @index(Local, Linear)
    tx = J_ + (@groupsize()[1] * (I_ - 1))

    if tx <= n
        dest[tx] = src[tx]
    end
end

@kernel function update_multiplier_kernel(n::Int, l_curr,
                                  u_curr,
                                  v_curr,
                                  rho)
    I_ = @index(Group, Linear)
    J_ = @index(Local, Linear)
    tx = J_ + (@groupsize()[1] * (I_ - 1))

    if tx <= n
        l_curr[tx] += rho[tx] * (u_curr[tx] - v_curr[tx])
    end
end

@kernel function primal_residual_kernel(n::Int, rp,
                                u_curr,
                                v_curr)
    I_ = @index(Group, Linear)
    J_ = @index(Local, Linear)
    tx = J_ + (@groupsize()[1] * (I_ - 1))

    if tx <= n
        rp[tx] = u_curr[tx] - v_curr[tx]
    end
end

@kernel function dual_residual_kernel(n::Int, rd,
                              v_prev,
                              v_curr,
                              rho)
    I_ = @index(Group, Linear)
    J_ = @index(Local, Linear)
    tx = J_ + (@groupsize()[1] * (I_ - 1))

    if tx <= n
        rd[tx] = -rho[tx] * (v_curr[tx] - v_prev[tx])
    end
end

function check_linelimit_violation(data::OPFData, u_; device=KA.CPU())
    lines = data.lines
    nline = length(data.lines)
    line_start = 2*length(data.generators) + 1

    rateA_nviols = 0
    rateA_maxviol = 0.0
    rateC_nviols = 0
    rateC_maxviol = 0.0
    u = u_ |> Array

    for l=1:nline
        pij_idx = line_start + 8*(l-1)
        ij_val = 0.0
        ji_val = 0.0
        ij_val = u[pij_idx]^2 + u[pij_idx+1]^2
        ji_val = u[pij_idx+2]^2 + u[pij_idx+3]^2

        limit = (lines[l].rateA / data.baseMVA)^2
        if limit > 0 && limit < 1e10
            if ij_val > limit || ji_val > limit
                rateA_nviols += 1
                rateA_maxviol = max(rateA_maxviol, max(ij_val - limit, ji_val - limit))
            end
        end

        limit = (lines[l].rateC / data.baseMVA)^2
        if limit > 0 && limit < 1e10
            if ij_val > limit || ji_val > limit
                rateC_nviols += 1
                rateC_maxviol = max(rateC_maxviol, max(ij_val - limit, ji_val - limit))
            end
        end
    end
    rateA_maxviol = sqrt(rateA_maxviol)
    rateC_maxviol = sqrt(rateC_maxviol)

    return rateA_nviols, rateA_maxviol, rateC_nviols, rateC_maxviol
end

"""
    admm_restart!

This function restarts the ADMM with a given `env::AdmmEnv` containing solutions and all the other parameters.
"""
admm_restart!(env::AdmmEnv; options...) = admm_solve!(env, env.solution; options...)

function admm_solve!(env::AdmmEnv, sol::SolutionOneLevel; iterlim=800, scale=1e-4)

    data, par, mod = env.data, env.params, env.model

    shift_lines = 0
    shmem_size = sizeof(Float64)*(14*mod.n+3*mod.n^2) + sizeof(Int)*(4*mod.n)

    nblk_gen = div(mod.gen_mod.ngen, 32, RoundUp)
    nblk_br = mod.nline
    nblk_bus = div(mod.nbus, 32, RoundUp)

    it = 0
    time_gen = time_br = time_bus = 0.0
    primres = 0.0
    dualres = 0.0
    eps_pri = 0.0
    eps_dual = 0.0
    auglag_it = 0
    tron_it = 0
    gpu_primres = 0.0
    gpu_dualres = 0.0
    gpu_eps_pri = 0.0
    gpu_eps_dual = 0.0

    @time begin
    while it < iterlim
        it += 1

        if isa(env.device, KA.CPU)
            sol.u_prev .= sol.u_curr
            sol.v_prev .= sol.v_curr
            sol.l_prev .= sol.l_curr

            tcpu = generator_kernel(mod.gen_mod, data.baseMVA, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho, env.device)
            time_gen += tcpu.time

            if env.use_polar
                tcpu = @timed auglag_it, tron_it = polar_kernel_cpu(mod.n, mod.nline, mod.line_start, scale,
                                                                    sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                    shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                    mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            else
                tcpu = @timed auglag_it, tron_it = auglag_kernel_cpu(mod.n, mod.nline, it, par.max_auglag, mod.line_start, par.mu_max,
                                                                     sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                     shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                     mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            end
            time_br += tcpu.time

            tcpu = @timed bus_kernel_cpu(data.baseMVA, mod.nbus, mod.gen_mod.gen_start, mod.line_start,
                                         mod.FrStart, mod.FrIdx, mod.ToStart, mod.ToIdx, mod.GenStart,
                                         mod.GenIdx, mod.Pd, mod.Qd, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho, mod.YshR, mod.YshI)
            time_bus += tcpu.time

            sol.l_curr .+= sol.rho .* (sol.u_curr .- sol.v_curr)
            sol.rd .= -sol.rho .* (sol.v_curr .- sol.v_prev)
            sol.rp .= sol.u_curr .- sol.v_curr
            #sol.rp_old .= sol.u_prev .- sol.v_prev

            primres = norm(sol.rp)
            dualres = norm(sol.rd)

            eps_pri = sqrt(length(sol.l_curr))*par.ABSTOL + par.RELTOL*max(norm(sol.u_curr), norm(-sol.v_curr))
            eps_dual = sqrt(length(sol.u_curr))*par.ABSTOL + par.RELTOL*norm(sol.l_curr)

            (par.verbose > 1) && @printf("[CPU] %10d  %.6e  %.6e  %.6e  %.6e  %6.2f  %6.2f\n",
                    it, primres, dualres, eps_pri, eps_dual, auglag_it, tron_it)

            if primres <= eps_pri && dualres <= eps_dual
                break
            end
        else
            wait(copy_data_kernel(env.device, 64, mod.nvar)(mod.nvar, sol.u_prev, sol.u_curr, dependencies=Event(env.device)))
            wait(copy_data_kernel(env.device, 64, mod.nvar)(mod.nvar, sol.v_prev, sol.v_curr, dependencies=Event(env.device)))
            wait(copy_data_kernel(env.device, 64, mod.nvar)(mod.nvar, sol.l_prev, sol.l_curr, dependencies=Event(env.device)))

            tgpu = generator_kernel(mod.gen_mod, data.baseMVA, sol.u_curr, sol.v_curr, sol.l_curr, sol.rho, env.device)

            # time_gen += tgpu.time
            # if env.use_polar
                wait(polar_kernel(env.device, 32, mod.nline*32)(mod.n, mod.nline, mod.line_start, scale,
                                                                sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
                                                                shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
                                                                mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound,
                                                                dependencies=Event(env.device)
                                                                )
                )
            # else
            #     tgpu = CUDA.@timed @cuda threads=32 blocks=nblk_br shmem=shmem_size auglag_kernel(mod.n, it, par.max_auglag, mod.line_start, scale, par.mu_max,
            #                                                                                sol.u_curr, sol.v_curr, sol.l_curr, sol.rho,
            #                                                                                shift_lines, env.membuf, mod.YffR, mod.YffI, mod.YftR, mod.YftI,
            #                                                                                mod.YttR, mod.YttI, mod.YtfR, mod.YtfI, mod.FrBound, mod.ToBound)
            # end
            # time_br += tgpu.time
            wait(bus_kernel(env.device, 32, mod.nbus)(data.baseMVA, mod.nbus, mod.gen_mod.gen_start, mod.line_start,
                                                                           mod.FrStart, mod.FrIdx, mod.ToStart, mod.ToIdx, mod.GenStart,
                                                                           mod.GenIdx, mod.Pd, mod.Qd, sol.u_curr, sol.v_curr, sol.l_curr,
                                                                           sol.rho, mod.YshR, mod.YshI,
                                                                           dependencies=Event(env.device)
                                                                           )
            )
            # time_bus += tgpu.time

            wait(update_multiplier_kernel(env.device, 32, mod.nvar)(mod.nvar, sol.l_curr,
                    sol.u_curr, sol.v_curr, sol.rho,
                    dependencies=Event(env.device)
                )
            )
            wait(primal_residual_kernel(env.device, 32, mod.nvar)(mod.nvar, sol.rp, sol.u_curr, sol.v_curr,
                    dependencies=Event(env.device)
                )
            )
            wait(dual_residual_kernel(env.device, 32, mod.nvar)(mod.nvar, sol.rd, sol.v_prev, sol.v_curr, sol.rho,
                    dependencies=Event(env.device)
                )
            )

            gpu_primres = norm(sol.rp)
            gpu_dualres = norm(sol.rd)

            gpu_eps_pri = sqrt(length(sol.l_curr))*par.ABSTOL + par.RELTOL*max(norm(sol.u_curr), norm(sol.v_curr))
            gpu_eps_dual = sqrt(length(sol.u_curr))*par.ABSTOL + par.RELTOL*norm(sol.l_curr)

            (par.verbose > 1) && @printf("[GPU] %10d  %.6e  %.6e  %.6e  %.6e\n", it, gpu_primres, gpu_dualres, gpu_eps_pri, gpu_eps_dual)
            if Base.mod(it, 1000) == 0
                (par.verbose > 0) && @printf("[GPU] %10d  %.6e  %.6e  %.6e  %.6e\n", it, gpu_primres, gpu_dualres, gpu_eps_pri, gpu_eps_dual)
            end

            if gpu_primres <= gpu_eps_pri && gpu_dualres <= gpu_eps_dual
                break
            end
        end
    end
    end

    u_curr = zeros(mod.nvar)
    copyto!(u_curr, sol.u_curr)
    objval = sum(data.generators[g].coeff[data.generators[g].n-2]*(data.baseMVA*u_curr[mod.gen_mod.gen_start+2*(g-1)])^2 +
                 data.generators[g].coeff[data.generators[g].n-1]*(data.baseMVA*u_curr[mod.gen_mod.gen_start+2*(g-1)]) +
                 data.generators[g].coeff[data.generators[g].n]
                 for g in 1:mod.gen_mod.ngen)::Float64
    sol.objval = objval

    if it < iterlim
        sol.status = HAS_CONVERGED
    else
        sol.status = MAXIMUM_ITERATIONS
    end

    if par.verbose > 0
        rateA_nviols, rateA_maxviol, rateC_nviols, rateC_maxviol = check_linelimit_violation(data, u_curr)
        if isa(env.device, KA.GPU)
            @printf("[GPU] %10d  %.6e  %.6e  %.6e  %.6e\n", it, gpu_primres, gpu_dualres, gpu_eps_pri, gpu_eps_dual)
        end
        if isa(env.device, KA.CPU)
            @printf("[CPU] %10d  %.6e  %.6e  %.6e  %.6e  %6.2f  %6.2f\n", it, primres, dualres, eps_pri, eps_dual, auglag_it, tron_it)
        end
        status = sol.status == HAS_CONVERGED ? "converged" : "not converged"
        @printf(" ** Line limit violations\n")
        @printf("RateA number of violations = %d (%d)\n", rateA_nviols, mod.nline)
        @printf("RateA maximum violation    = %.2f\n", rateA_maxviol)
        @printf("RateC number of violations = %d (%d)\n", rateC_nviols, mod.nline)
        @printf("RateC maximum violation    = %.2f\n", rateC_maxviol)

        @printf(" ** Time\n")
        @printf("Generator = %.2f\n", time_gen)
        @printf("Branch    = %.2f\n", time_br)
        @printf("Bus       = %.2f\n", time_bus)
        @printf("Total     = %.2f\n", time_gen + time_br + time_bus)

        @printf("Objective value = %.6e\n", objval)
        println("Solver status: $status")
        println("Iterations: $it")
    end
    return
end

function admm_rect_gpu(case::String; iterlim=800, rho_pq=400.0, rho_va=40000.0, scale=1e-4,
                       device=KA.CPU(), use_polar=true, gpu_no=0, verbose=1)
    env = AdmmEnv(case, device, rho_pq, rho_va; use_polar=use_polar, gpu_no=gpu_no, verbose=verbose,)
    admm_restart!(env, iterlim=iterlim, scale=scale)
    return env
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7118f3c-209e-11eb-1f79-bdc2dd17aea0
begin
	using Revise
	using PlutoMustache
	using TropicalTensors
	using Random
	using LightGraphs
	using Compose, Viznet
	using SimpleTensorNetworks
	using HDF5
	Compose.set_default_graphic_size(14cm, 14cm)
end

# ╔═╡ 73cff25e-5434-11eb-0147-21ac9dbcb8b5
function rand_3regular_tn(::Type{T}, n; D=2) where T
	g = LightGraphs.random_regular_graph(n, 3)
	labels = 1:ne(g)
	arrays = [rand(T, fill(2, 3)...) for i=1:n]
	labels = [Int[] for i=1:n]
	for (k, e) in enumerate(edges(g))
		push!(labels[e.src], k)
		push!(labels[e.dst], k)
	end
	tensors = LabeledTensor.(arrays, labels)
	TensorNetwork(tensors)
end

# ╔═╡ 6b01ff60-217a-11eb-2a5c-13eb2716c727
md"number of nodes = $(@bind nnodes Slider(2:2:40, default=16))"

# ╔═╡ 2d5dd844-2136-11eb-176f-e53e2ac97d82
tnr = rand_3regular_tn(Float64, nnodes; D=2)

# ╔═╡ 8d7cb70c-5441-11eb-35cb-e551b9166fd6
md"contract the graph by a random order"

# ╔═╡ 1aabc318-2137-11eb-33b7-cfcec36e631f
@bind step Button(label="step")

# ╔═╡ 72d7474c-2137-11eb-030b-eb894e4b22dc
@regstate step

# ╔═╡ 09ea2b64-2137-11eb-187a-fdf0f5b882ca
@when step begin
	length(tnr.tensors) > 1 && SimpleTensorNetworks.contract_label!(tnr, rand(vcat([t.labels for t in tnr.tensors]...)))
	tnr
end

# ╔═╡ Cell order:
# ╠═c7118f3c-209e-11eb-1f79-bdc2dd17aea0
# ╠═73cff25e-5434-11eb-0147-21ac9dbcb8b5
# ╟─6b01ff60-217a-11eb-2a5c-13eb2716c727
# ╠═2d5dd844-2136-11eb-176f-e53e2ac97d82
# ╟─8d7cb70c-5441-11eb-35cb-e551b9166fd6
# ╟─1aabc318-2137-11eb-33b7-cfcec36e631f
# ╟─72d7474c-2137-11eb-030b-eb894e4b22dc
# ╠═09ea2b64-2137-11eb-187a-fdf0f5b882ca

ext: jl


text:
notary = Notary()

guardiankey = newsigner(notary)
registratorkey = newsigner(notary)

member1key = newsigner(notary)
member2key = newsigner(notary)
member3key = newsigner(notary)

guardian = Guardian(id(guardiankey))
registrator = Registrator(id(registratorkey))
member1 = Member(id(member1key))
member2 = Member(id(member2key))
member3 = Member(id(member3key))

transactions = TransactionVector()
tlog = BraidChain(transactions, guardian, notary)

push!(tlog, Transaction(member1, guardiankey))
Test.@test_throws AssertionError push!(tlog, Transaction(member2, registratorkey))
push!(tlog, Transaction(registrator, guardiankey))
push!(tlog, Transaction(member2, registratorkey))
Test.@test_throws AssertionError push!(tlog, Transaction(member2, registratorkey))

pseudonyms(tlog, 3)

### Braiders

braider = Braider(ID(1332))
push!(tlog, Transaction(braider, guardiankey))

m = MetaData(ID(2323), braider.id, Hash(tlog))

pseudonym1 = newsigner(notary)
pseudonym2 = newsigner(notary)
pseudonym3 = newsigner(notary)

# generated by users
newpseudonyms = ID[id(pseudonym1), id(pseudonym2), id(pseudonym3)]

thisbraid = SynchronicBraid(m, newpseudonyms)

contract = Contract(Dict[Dict(sign(thisbraid, member1key)), Dict(sign(thisbraid, member2key)), Dict(sign(thisbraid, member3key))])

Test.@test_throws AssertionError push!(tlog, Transaction(thisbraid, contract)) ### Should throw an error

push!(tlog, Transaction(member3, registratorkey))
push!(tlog, Transaction(thisbraid, contract)) ### Should throw an error

### New BraidChain

bcnew = BraidChain(tlog.transactions, guardian, notary)
@test pseudonyms(tlog) == pseudonyms(bcnew)


text:
plume_z_params = Dict(
    "A" => (δ=107.7, β=-1.7172, γ=0.2770),
    "B" => (δ=0.1355, β=0.8752, γ=0.0136),
    "C" => (δ=0.09623, β=0.9477, γ=-0.0020),
    "D" => (δ=0.04134, β=1.1737, γ=-0.0316),
    "E" => (δ=0.02275, β=1.3010, γ=-0.0450),
    "F" => (δ=0.01122, β=1.4024, γ=-0.0540)
)

puff_z_params = Dict(
    "A" => (δ=0.60, β=0.75),
    "B" => (δ=0.53, β=0.73),
    "C" => (δ=0.34, β=0.71),
    "D" => (δ=0.15, β=0.70),
    "E" => (δ=0.10, β=0.65),
    "F" => (δ=0.05, β=0.61),
)

plume_y_params = Dict(
    "A" => (δ=0.423, β=0.9, tₐ=18.4),
    "B" => (δ=0.313, β=0.9, tₐ=18.4),
    "C" => (δ=0.210, β=0.9, tₐ=18.4),
    "D" => (δ=0.136, β=0.9, tₐ=18.3),
    "E" => (δ=0.102, β=0.9, tₐ=11.4),
    "F" => (δ=0.0674, β=0.9, tₐ=4.6)
)

puff_y_params = Dict(
    "A" => (δ=0.18, β=0.92),
    "B" => (δ=0.14, β=0.92),
    "C" => (δ=0.10, β=0.92),
    "D" => (δ=0.06, β=0.92),
    "E" => (δ=0.04, β=0.92),
    "F" => (δ=0.02, β=0.89),
)

"""
    crosswind_dispersion(stability_class::String; <keyword arguments>)
returns the crosswind dispersion function σy(x) for a given Pasquill-Gifford
stability class, `x` is assumed to be in meters and `σy` is in meters

# Arguments
`plume` determines if the dispersion is for a plume or instantaneous release
(puff), by default a plume is assumed.
`avg_time` the averaging time in seconds, falls back to an "instantaneous"
averaging time

# References
plume dispersion correlations are from
Spicer, T.O., and J.A. Havens, *User's Guide for the DEGADIS 2.1 Dense Gas
Dispersion Model*, EPA-450/4-89-019, November 1989, pp 45-46

puff dispersion correlations are from:
Slade, D.H., Diffusion from Instantaneous Sources. *Meteorology and Atomic
Energy*, TID-24190, USAEC, 1968, pp 163-175
"""
function crosswind_dispersion(stability_class::String; plume=true, avg_time=600)
    if stability_class ∈ Set(["A","B","C","D","E","F"])
        if plume
            δ, β, tₐ = plume_y_params[stability_class]
            δ = δ*(max(avg_time, tₐ)/600)^0.2
        else
            δ, β = puff_y_params[stability_class]
        end
    else
        err = "$stability_class is not a valid Pasquill-Gifford stability class"
        error(err)
    end

    return x -> δ*x^β
end

"""
vertical_dispersion(stability_class::String, plume=true)
returns the crosswind dispersion function σz(x) for a given Pasquill-Gifford
stability class
`x` is assumed to be in meters and `σz` is in meters
`plume` determines if the dispersion is for a plume or instantaneous release
(puff), by default a plume is assumed.

# References
plume dispersion correlations are from
Seinfeld, J.H., *Atmospheric Chemistry and Physics of Air Pollution*, John
Wiley and Sons, New York, 1986

puff dispersion correlations are from:
Slade, D.H., Diffusion from Instantaneous Sources. *Meteorology and
Atomic Energy*, TID-24190, USAEC, 1968, pp163-175
"""
function vertical_dispersion(stability_class::String; plume=true)
    if stability_class ∈ Set(["A","B","C","D","E","F"])
        if plume
            δ, β, γ = plume_z_params[stability_class]
            return x -> δ*(x^β)*exp(γ*log(x)^2)
        else
            δ, β = puff_z_params[stability_class]
            return x -> δ*x^β
        end
    else
        err = "$stability_class is not a valid Pasquill-Gifford stability class"
        error(err)
    end
end


text:
include("./superres.jl")

using Plots
pyplot();

let ntests=15000,seed=5,Ωₘ=100
    let ℓ=2,R=8,clusterFun=Ζ,amplitudeFun = imag1s
        let pertFun=𝔉ᵣ(Ωₘ,seed), ni = 6
            let ϵR = 10.0.^LinRange(-4,3,ntests), ΩR = LinRange(20,40,ntests), ΔR = LinRange(1e-4,2e-3,ntests)#10.0.^LinRange(-5,-1,ntests)
#            let ϵR = 10.0.^LinRange(-2,1,ntests), ΩR = NrangeLog(1.2,2,ntests), ΔR = 10.0.^LinRange(-5,-1,ntests)
                let res=mapreduce(transpose,vcat,
                        [runSingle_ind(ℓ,R,clusterFun, pertFun, amplitudeFun, ni, ϵR,ΩR,ΔR) for _ in 1:ntests])
                    let clr = map(cc,res[:,2]), sh = map(ss,res[:,2]),idxsSucc = findall(res[:,2].==1),idxsFail=findall(res[:,2].==0)
                        #p=scatter(log10.(res[:,1]),log10.(inv.(res[:,3])),
                        #    mc=clr,msc=:white,ms=7,shape=sh,markerstrokewidth=0.4,markeralpha=0.5,
                        #    xlabel="log SRF",ylabel="log 1/ϵ",lab="",color_palette=:magma)
                        p=scatter(log10.(res[idxsSucc,1]),log10.(inv.(res[idxsSucc,3])),
                            mc=:blue,msc=:white,ms=7,shape=:dtriangle,markerstrokewidth=0.1,markeralpha=0.1,
                            xlabel="log SRF",ylabel="log 1/ϵ",lab="Success",color_palette=:magma)
                        scatter!(log10.(res[idxsFail,1]),log10.(inv.(res[idxsFail,3])),
                            mc=:red,msc=:white,ms=7,shape=:circle,markerstrokewidth=0.1,markeralpha=1.0,
                            xlabel="log SRF",ylabel="log 1/ϵ",lab="Fail",color_palette=:magma)
                        #plot!([t->(2*ℓ-2)*t+5,t->(2*ℓ-1)*t+5],w=3,l=:dash,lc=[:black :blue],lab=["SRF^(2-2ℓ)" "SRF^(1-2ℓ)"])
                        title!("p=$(ℓ),d=$(R), j=$(ni), \$ {\\bf S_1} \$")
                        #ylims!(-.3,0.3)
                        xlims!(1.5,2.5)
                        savefig(p,"./figures/fig6.pdf")
                        #display(sum(res[:,2].==0))
                    end
                end
            end
        end
    end
end

max_stars_repo_name: batenkov-team/iai-code

max_issues_repo_name: batenkov-team/iai-code

max_forks_repo_name: batenkov-team/iai-code


text:
let
bounds = 1100

fabric = zeros(Int, bounds, bounds)
claims = zeros(Int, bounds, bounds)

 function unClaim(claimID)
   for x in 1:bounds, y in 1:bounds
     if claims[x, y] == claimID
       claims[x, y] = -1
     end
   end
 end

 open("d3data.txt") do file
    for line in eachline(file)
      line = replace(line, ":" => "")

      tokens = split(line, " ")
      claimID = parse(Int, replace(tokens[1], "#" => ""))
      pos = split(tokens[3], ",")
      size = split(tokens[4], "x")

      for x in parse(Int, pos[1]): parse(Int, pos[1]) + parse(Int, size[1])-1
        for y in parse(Int, pos[2]): parse(Int, pos[2]) + parse(Int, size[2])-1

          claim = claims[x+1, y+1]
          if claim == 0
            claims[x+1, y+1] = claimID
          elseif claim > 0
            unClaim(claim)
            unClaim(claimID)
            claimID = -1
          end
          fabric[x+1, y+1] = fabric[x+1, y+1] +1
        end
      end
    end
  end
  ctr = 0

  for x in 1:bounds, y in 1:bounds
    if fabric[x, y] > 1
      ctr += 1
    end
    if claims[x, y] > 0
      println(claims[x, y])
      break
    end
  end

  println(ctr)
end


text:
#= Create ODEProblem from the selected Model and Simulation options, and also the selected
    vehicle configurations
=#
import JPLEphemeris: SPK, Dates, position, velocity, state
include("..\\Vehicle\\Spacecraft.jl")
include("..\\Environment\\Environment.jl")
include("State.jl")

#Create EOM from the vehicle configurations
function EOM(vehicle::Type{<:AbstractSpacecraft}, env::Environment)

    # Gravity (Environment-based)
    #

    # VehicleShapeAerodynamics (Vehicle properties, combined with environmental properties)

    # VehicleActuators (Vehicle properties, combined with environmental properties)

    eomfunc = function (y::SimState, t::Float64)

    end #eomfunc
    return eomfunc
end #EOM


text:
datadir = joinpath(@__DIR__, "data")
examples = joinpath(@__DIR__, "../examples")

@testset "Assembler tests" begin
    @testset "Symbol table" begin
        syms = symboltable(joinpath(datadir, "labels.lmc"))
        target = Dict("zero" => 0, "one" => 1, "two" => 2, "three" => 3)
        for (k, v) in target
            @test haskey(syms, k)
            @test syms[k] == v
        end
    end

    @testset "Instruction set" begin
        instructions = [100, 200, 300, 500, 600,
                        700, 800, 901, 902, 000]
        filepath = joinpath(datadir, "instructionset.lmc")
        program = assemble(filepath)
        for (i, expected) in enumerate(instructions)
            @test program[i] == expected
        end
    end

	@testset "Regression test" begin
		files = map(filter(endswith(".machine"), readdir(examples))) do filename
			splitext(filename)[1]
		end
		for file in files
			machinefile = joinpath(examples, string(file, ".machine"))
			srcfile 	= joinpath(examples, string(file, ".lmc"))
			@test isfile(machinefile)
			@test isfile(srcfile)

			target = load(machinefile)
			memory = assemble(srcfile)
			@test memory == target
		end
	end

    @testset "Assemble mnemonic" begin
        using LittleManComputer: assemble_mnemonic

        labels = Dict("foo" => 42, "bar" => 7)

        @test assemble_mnemonic(["ADD", "12"]) == 112
        @test assemble_mnemonic(["SUB", "4"]) == 204
        @test assemble_mnemonic(["ADD", "foo"], labels) == 142
        @test assemble_mnemonic(["SUB", "bar"], labels) == 207

        @test assemble_mnemonic(["BRZ", "bar"], labels) == 707
        @test assemble_mnemonic(["BRZ", "foo"], labels) == 742
        @test assemble_mnemonic(["BRZ", "55"]) == 755
        @test assemble_mnemonic(["LDA", "37"]) == 537
        @test assemble_mnemonic(["LDA", "foo"], labels) == 542
    end

end


text:
struct AddressTable
    readers::Threads.Atomic{UInt}
    entries::Vector{Union{Actor, Nothing}}
    rev_entries::IdDict{Actor, Set{UInt32}}
end

AddressTable() = AddressTable(Threads.Atomic{UInt}(), [], IdDict{Actor, Set{UInt32}}())

struct AddressBook
    write_lock::ReentrantLock
    flips::Threads.Atomic{UInt}
    live::NTuple{2, AddressTable}
end

AddressBook() = AddressBook(ReentrantLock(), Threads.Atomic{UInt}(),
                            (AddressTable(), AddressTable()))

atomic_inc!(atom::Threads.Atomic{T}) where T = Threads.atomic_add!(atom, T(1))
atomic_dec!(atom::Threads.Atomic{T}) where T = Threads.atomic_sub!(atom, T(1))

function read_lock(fn::Function, book::AddressBook)
    retry = 100
    flips = 0
    flop = 0

    while retry > 0
        flips = book.flips[]
        flop = flips % 2 + 1
        atomic_inc!(book.live[flop].readers)
        flips == book.flips[] && break

        atomic_dec!(book.live[flop].readers)
        retry -= 1
    end

    retry > 0 || error("Could not do address lookup due to contention")

    try
        fn(book.live[flop])
    finally
        atomic_dec!(book.live[flop].readers)
    end
end

Base.getindex(book::AddressBook, id::Id) = getindex(book, id.inner)
function Base.getindex(book::AddressBook, id::Integer)
    a = read_lock(book) do table
        table.entries[id]
    end

    a === nothing && throw(KeyError(id))

    a
end

Base.getindex(book::AddressBook, a::Actor) = read_lock(book) do table
    collect(table.rev_entries[a])
end

Base.isempty(book::AddressBook) = read_lock(book) do table
    isempty(table.rev_entries)
end

Base.lastindex(book::AddressBook) = read_lock(book) do table
    lastindex(table.entries)
end

Base.iterate(book::AddressBook, state=1) = try
    (book[state], state + 1)
catch ex
    ex isa BoundsError || rethrow()

    nothing
end

Base.length(book::AddressBook) = read_lock(book) do table
    length(table.entries)
end

function write_lock(fn::Function, book::AddressBook)
    retries = 11
    retry = 0

    lock(book.write_lock)
    try
        table = book.live[2 - book.flips[] % 2]

        @assert table.readers[] == 0 "Previous flip should have drained readers"

        fn(table)

        table = book.live[2 - (atomic_inc!(book.flips) + 1) % 2]

        while table.readers[] > 0
            @assert retry < retries "Waited a long time to drain readers: $(table.readers[])"
            sleep(2^retry * 0.0001)
            retry += 1
        end

        fn(table)
    finally
        unlock(book.write_lock)
    end
end

Base.push!(book::AddressBook, a::Actor{S, M}) where {S, M} = write_lock(book) do table
    i = lastindex(table.entries) + 1

    @assert i < typemax(UInt32) "Can't add actor, no local addresses left"

    push!(table.entries, a)

    # For some reason get!() was behaving like get()
    rev_entries = get(Set, table.rev_entries, a)
    push!(rev_entries, UInt32(i))
    table.rev_entries[a] = rev_entries

    Id{S, M}(UInt32(i))
end

Base.setindex!(book::AddressBook, ::Nothing, id::Id) = write_lock(book) do table
    old_a = table.entries[id.inner]

    rentries = delete!(table.rev_entries[old_a], id.inner)
    isempty(rentries) && delete!(table.rev_entries, old_a)

    table.entries[id.inner] = nothing
end

Base.setindex!(book::AddressBook, a::Actor, id::Id) = write_lock(book) do table
    old_a = table.entries[id.inner]

    rentries = delete!(table.rev_entries[old_a], id.inner)
    isempty(rentries) && delete!(table.rev_entries, old_a)

    table.entries[id.inner] = a

    rev_entries = get(Set, table.rev_entries, a)
    push!(rev_entries, id.inner)
    table.rev_entries[a] = rev_entries

    id
end


text:
"""
`vietorisrips(g::AbstractSparseMatrix{U, T}, dim_max = 2) where {U, T}`

Return a filtered Vietoris Rips complex (up to dimension `dim_max`) from a sparse distance matrix `g`. Only the upper triangular part of `g` is used. Zero values of the matrix are ignored: to obtain the rips filtration up to a threshold `t`, simply remove all entries of `g` greater than `t`. Discretizing values of entries in `g` may also improve performance.

The filtered complex is a tuple (where the first element corresponds to 0-simplices, the second element to 1-simplices and so on). The filtration on n-simplices is represented (slightly improperly) as a `U`-valued n-`Cochain`.

## Examples

```julia
julia> using SparseArrays

julia> M = sparse([1, 2, 1], [2, 3, 3], [0.3, 0.2, 0.4], 3, 3)
3×3 SparseMatrixCSC{Float64,Int64} with 3 stored entries:
  [1, 2]  =  0.3
  [1, 3]  =  0.4
  [2, 3]  =  0.2

julia> cplx = vietorisrips(M, 2)
(Float64-valued 0-Cochain, Float64-valued 1-Cochain, Float64-valued 2-Cochain)

julia> keys(cplx[2])
3-element StructArrays.StructArray{Tuple{Int64,Int64},1,NamedTuple{(:x1, :x2),Tuple{Array{Int64,1},Array{Int64,1}}}}:
 (2, 1)
 (3, 1)
 (3, 2)

julia> values(cplx[2])
3-element Array{Float64,1}:
 0.3
 0.4
 0.2
```
"""
function vietorisrips(g::AbstractSparseMatrix, dim_max = 2)
    lt = sparse_uppertriangular(g)
    vietorisrips_uppertriangular(lt, dim_max)
end

function sparse_uppertriangular(g::AbstractSparseMatrix)
    size1, size2 = size(g)
    @assert size1 == size2
    row, col, vals = findnz(g)
    mask = row .< col
    sparse(row[mask], col[mask], vals[mask], size1, size2)
end

function vietorisrips_uppertriangular(g::SparseMatrixCSC{U, T}, dim_max = 2) where {U, T}
    v = ntuple(dim_max+1) do npts
        Cochain(StructVector{NTuple{npts, T}}(undef, 0), U[])
    end
    compute_simplices!(v, g, dim_max)
    return v
end

function compute_simplices!(v, g::SparseMatrixCSC{U, T}, dim_max) where {U, T}
    us = one(T):T(size(g, 1))
    for u in us
        ptr = g.colptr[u]:g.colptr[u+1]-1
        nu = view(g.rowval, ptr)
        ws = view(g.nzval, ptr)
        add_cofaces!(v, g, (u,), zero(U), StructVector((nu, ws)), dim_max)
    end
end

function add_cofaces!(v, g, sim, sim_weight, weighted_neighbors, dim_max)
    cocycle = v[length(sim)]
    push!(keys(cocycle), sim)
    push!(values(cocycle), sim_weight)

    (length(sim) > dim_max) && return

    for (n, w) in weighted_neighbors
        ptr = g.colptr[n]:g.colptr[n+1]-1
        new_weighted_neighbors = neighbors_intersection(g.rowval[ptr], g.nzval[ptr], weighted_neighbors)
        add_cofaces!(v, g, (sim..., n), max(w, sim_weight), new_weighted_neighbors, dim_max)
    end
end

function neighbors_intersection(neighbors, weights, weighted_neighbors)
    new_weighted_neighbors = similar(weighted_neighbors, 0)
    for (nb, nb_weight) in weighted_neighbors
        idx = searchsortedfirst(neighbors, nb)
        idx in axes(neighbors, 1) || continue
        neighbors[idx] == nb || continue
        push!(new_weighted_neighbors, (nb, max(nb_weight, weights[idx])))
    end
    new_weighted_neighbors
end


text:
#module MPSReader

using JuMP
using MathOptInterface
using LinearAlgebra

#export readmps, mpstomodel

function mpstomodel(filename::String, fixed::Bool = true)
    varBin, varInt, varFloat, bounds, objsense, c, c0, Aeq, beq, Aineq, bineq = mpstomatrices(filename, fixed)

    nVar = length(varBin) + length(varInt) + length(varFloat)

    myModel = Model()
    @variable(myModel, x[1:nVar])
    if !isapprox(c0, 0)
        if objsense == :min
            @objective(myModel, Min, dot(c, x) + c0)
        else
            @objective(myModel, Max, dot(c, x) + c0)
        end
    else
        if objsense == :min
            @objective(myModel, Min, dot(c, x))
        else
            @objective(myModel, Max, dot(c, x))
        end
    end
    # Get type to model and bounds/fixed
    for i in varBin
        JuMP.set_binary(x[i])
        if bounds[i][1] == bounds[i][2]
            JuMP.fix(x[i], bounds[i][1])
        end
    end
    for i in varInt
        JuMP.set_integer(x[i])
        if bounds[i][1] == bounds[i][2]
            JuMP.fix(x[i], bounds[i][1])
        else
            if isfinite(bounds[i][1])
                JuMP.set_lower_bound(x[i], bounds[i][1])
            end
            if isfinite(bounds[i][2])
                JuMP.set_upper_bound(x[i], bounds[i][2])
            end
        end
    end
    for i in varFloat
        if bounds[i][1] == bounds[i][2]
            JuMP.fix(x[i], bounds[i][1])
        else
            if isfinite(bounds[i][1])
                JuMP.set_lower_bound(x[i], bounds[i][1])
            end
            if isfinite(bounds[i][2])
                JuMP.set_upper_bound(x[i], bounds[i][2])
            end
        end
    end

    # Add constraints
    m = length(bineq)
    if objsense == :min
        @constraint(myModel, [i = 1:m], dot(Aineq[i,:], x) >= bineq[i])
    else
        @constraint(myModel, [i = 1:m], dot(Aineq[i,:], x) <= bineq[i])
    end
    m = length(beq)
    @constraint(myModel, [i = 1:m], dot(Aeq[i,:], x) == beq[i])

    return myModel, x
end

function mpstomatrices(filename::String, fixed::Bool = true)
    varTypes, bounds, objsense, c, c0, A, b, conTypes = readmps(filename, fixed)

    # Get precise var type
    varBin = Vector{Int64}()
    varInt = Vector{Int64}()
    varFloat = Vector{Int64}()

    nVar = length(varTypes)
    nCon = length(conTypes)

    for i in 1:nVar
        isFloat = varTypes[i]
        if isFloat
            push!(varFloat, i)
        else
            if isapprox(0, bounds[i][1]) && isapprox(1, bounds[i][2])
                push!(varBin, i)
            else
                push!(varInt, i)
            end
        end
    end

    # Aineq x >= bineq in min
    ismin = 1
    # Aineq x <= bineq in max
    if objsense == :max
        ismin = -1
    end

    # Split A and b in equality and inequality
    Aeq = Matrix{Float64}(undef, 0, nVar)
    beq = Vector{Float64}()
    Aineq = Matrix{Float64}(undef, 0, nVar)
    bineq = Vector{Float64}()

    for i in 1:nCon
        if conTypes[i] == 2
            Aeq = vcat(Aeq, A[i,:]')
            push!(beq, b[i])
        elseif conTypes[i] == 1
            Aineq = vcat(Aineq, -ismin * A[i,:]')
            push!(bineq, -ismin * b[i])
        elseif conTypes[i] == 3
            Aineq = vcat(Aineq, ismin * A[i,:]')
            push!(bineq, ismin * b[i])
        else
            error("Should never happen")
        end
    end

    return varBin, varInt, varFloat, bounds, objsense, c, c0, Aeq, beq, Aineq, bineq
end

function readmps(filename::String, fixed::Bool = true)
    sections = Dict{String, Bool}([
        "NAME" => false,
        "ROWS" => false,
        "COLUMNS" => false,
        "RHS" => false,
        "BOUNDS" => false,
        "RANGES" => false,
        "SOS" => false,
        "ENDATA" => false
    ])
    if !fixed
        sections["OBJSENSE"] = false
        sections["OBJNAME"] = false
    end

    varNames = Dict{String, Int}()
    conNames = Dict{String, Int}()
    varTypes = BitArray{1}() # true : float   false : int
    conTypes = Vector{Int}() # 1 : <=, 2 : ==, 3 : >=

    # Default values
    name = ""
    objectiveName = ""
    objsense = :min
    nVar = 0
    nCon = 0
    A = Matrix{Float64}(undef,0,0)
    b = Vector{Float64}()
    c = Vector{Float64}()
    c0 = 0.0
    bounds = Vector{Tuple{Float64, Float64}}()

    open(filename) do mps
        line = ""
        while !eof(mps) && !sections["ENDATA"]
            line = readline(mps)
            # Read the section name and avoid empty lines and comments
            while strip(line) == "" || line[1] == '*'
                line = readline(mps)
            end
            words = split(line)
            section = words[1]

            # Errors in sections names and order
            if !haskey(sections, section)
                error("Section ", section, " is not a valide section")
            elseif sections[section]
                error("Section ", section, " appears twice")
            elseif section == "COLUMNS" && !sections["ROWS"]
                error("ROWS must come before COLUMNS")
            elseif section == "RHS" && (!sections["COLUMNS"] || !sections["NAME"])
                error("NAME and COLUMNS must come before RHS")
            elseif section == "BOUNDS" && !sections["COLUMNS"]
                error("COLUMNS must come before BOUNDS")
            elseif section == "RANGES" && !sections["RHS"]
                error("RHS must come before RANGES")
            end

            # Read info
            if section == "NAME"
                name = words[2]
            elseif section == "OBJSENSE"
                objsense = readObjsense(mps)
            elseif section == "OBJNAME"
                objectiveName = readObjname(mps)
            elseif section == "ROWS"
                objectiveName, conNames, conTypes = readRows(mps, objectiveName, fixed)
                nCon = length(conTypes)
            elseif section == "COLUMNS"
                varNames, c, A, varTypes = readColumns(mps, objectiveName, conNames, fixed)
                nVar = length(keys(varNames))
                bounds = fill((0.0, Inf), nVar)
            elseif section == "RHS"
                b, c0 = readRhs(mps, objectiveName, conNames, fixed)
            elseif section == "RANGES"
                A, b, conTypes = readRanges(mps, A, b, conNames, conTypes, fixed)
                nCon = length(conTypes)
            elseif section == "BOUNDS"
                bounds, varTypes = readBounds(mps, varNames, varTypes, fixed)
            elseif section == "SOS"
                error("SOS section reader not ready yet")
            end
            sections[section] = true
        end
        if !sections["ENDATA"]
            error("No ENDATA section in the file")
        end
    end

    return varTypes, bounds, objsense, c, c0, A, b, conTypes
end

function readBounds(mps, varNames, varTypesinit, fixed)
    varTypes = copy(varTypesinit)
    pos = position(mps)
    bounds = fill((0.0, Inf), length(keys(varNames)))

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end
    while line[1] == ' '
        pos = position(mps)
        words = Vector{String}()

        if fixed
            lenLine = length(line)
            words = strip.([
                line[2:min(lenLine, 3)],
                line[5:min(lenLine, 12)],
                line[15:min(lenLine, 22)],
                line[25:min(lenLine, 36)]
            ])
        else
            words = split(line)
        end

        boundType = words[1]
        if !fixed
            boundType = uppercase(boundType)
        end
        name = words[3]

        if boundType == "FR"
            # Free variable
            bounds[varNames[name]] = (-Inf, Inf)
        elseif boundType == "MI"
            # Lower bound -inf
            if isfinite(bounds[varNames[name]][2])
                bounds[varNames[name]] = (-Inf, bounds[varNames[name]][2])
            else
                bounds[varNames[name]] = (-Inf, 0.0)
            end
        elseif boundType == "PL"
            # Default value
            bounds[varNames[name]] = (bounds[varNames[name]][1], Inf)
        elseif boundType == "BV"
            # Binary variable
            bounds[varNames[name]] = (0.0, 1.0)
            varTypes[varNames[name]] = false
        elseif boundType == "SC"
            # Semi-continuous variable
            error("SC bound not ready yet")
        elseif boundType == "LO"
            # Lower bound
            boundVal = parse(Float64, words[4])
            bounds[varNames[name]] = (boundVal, bounds[varNames[name]][2])
        elseif boundType == "UP"
            # Upper bound
            boundVal = parse(Float64, words[4])
            bounds[varNames[name]] = (bounds[varNames[name]][1], boundVal)
        elseif boundType == "FX"
            # Fixed variable
            boundVal = parse(Float64, words[4])
            bounds[varNames[name]] = (boundVal, boundVal)
        elseif boundType == "LI"
            # Integer variable / Lower bound
            boundVal = parse(Float64, words[4])
            bounds[varNames[name]] = (boundVal, bounds[varNames[name]][2])
            varTypes[varNames[name]] = false
        elseif boundType == "UI"
            # Integer variable / Upper bound
            boundVal = parse(Float64, words[4])
            bounds[varNames[name]] = (bounds[varNames[name]][1], boundVal)
            varTypes[varNames[name]] = false
        else
            error("Unknown bound type ", boundType)
        end

        line = readline(mps)
        while strip(line) == "" || line[1] == '*'
            pos = position(mps)
            line = readline(mps)
        end
    end

    seek(mps, pos)
    return bounds, varTypes
end

function readRanges(mps, Ainit, binit, conNames, conTypesinit, fixed)
    A = copy(Ainit)
    b = copy(binit)
    conTypes = copy(conTypesinit)
    pos = position(mps)

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end
    while line[1] == ' '
        pos = position(mps)
        words = Vector{String}()

        if fixed
            lenLine = length(line)
            words = strip.([
                line[5:min(lenLine, 12)],
                line[15:min(lenLine, 22)],
                line[25:min(lenLine, 36)],
                line[40:min(lenLine, 47)],
                line[50:min(lenLine, 61)]
            ])
            while words[end] == ""
                pop!(words)
            end
        else
            words = split(line)
        end
        lenWords = length(words)

        for i in 2:2:lenWords
            name = words[i]
            if name != ""
                if conTypes[conNames[name]] == 2
                    error("Range on an equality constraint")
                elseif conTypes[conNames[name]] == 1
                    push!(conTypes, 3)
                    b0 = -abs(parse(Float64, words[i+1]))
                elseif conTypes[conNames[name]] == 3
                    push!(conTypes, 1)
                    b0 = abs(parse(Float64, words[i+1]))
                else
                    error("Should never happen")
                end
                A = vcat(A, A[conNames[name],:]')
                push!(b, b0)
            end
        end

        line = readline(mps)
        while strip(line) == "" || line[1] == '*'
            pos = position(mps)
            line = readline(mps)
        end
    end

    seek(mps, pos)
    return A, b, conTypes
end

function readRhs(mps, objectiveName, conNames, fixed)
    pos = position(mps)
    b = zeros(Float64, length(keys(conNames)))
    c0 = 0.0

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end
    while line[1] == ' '
        pos = position(mps)
        words = Vector{String}()

        if fixed
            lenLine = length(line)
            words = strip.([
                line[5:min(lenLine, 12)],
                line[15:min(lenLine, 22)],
                line[25:min(lenLine, 36)],
                line[40:min(lenLine, 47)],
                line[50:min(lenLine, 61)]
            ])
            while words[end] == ""
                pop!(words)
            end
        else
            words = split(line)
        end
        lenWords = length(words)

        for i in 2:2:lenWords
            name = words[i]
            if name == objectiveName
                if c0 != 0
                    error("Multiple presence of ", objectiveName, " in RHS")
                else
                    c0 = -parse(Float64, words[i+1])
                end
            elseif name != ""
                if b[conNames[name]] != 0
                    error("Multiple presence of ", name, " in RHS")
                else
                    b[conNames[name]] = parse(Float64, words[i+1])
                end
            end
        end

        line = readline(mps)
        while strip(line) == "" || line[1] == '*'
            pos = position(mps)
            line = readline(mps)
        end
    end

    seek(mps, pos)
    return b, c0
end

function readColumns(mps, objectiveName, conNames, fixed)
    pos = position(mps)
    varNames = Dict{String, Int}()
    varTypes = BitArray{1}() # true : float   false : int
    nCon = length(keys(conNames))
    A = zeros(Float64, nCon, 0)
    initVar = zeros(Float64, nCon, 1)
    c = Vector{Float64}()
    nVar = 1
    isFloat = true

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end
    while line[1] == ' '
        pos = position(mps)
        words = Vector{String}()

        if fixed
            lenLine = length(line)
            words = strip.([
                line[5:min(lenLine, 12)],
                line[15:min(lenLine, 22)],
                line[25:min(lenLine, 36)],
                line[40:min(lenLine, 47)],
                line[50:min(lenLine, 61)]
            ])
            while words[end] == ""
                pop!(words)
            end
        else
            words = split(line)
        end
        lenWords = length(words)

        # MARKER for INT
        if words[2] == "'MARKER'"
            ind = 3
            if fixed
                ind = 4
            end
            if words[ind] == "'INTORG'"
                isFloat = false
            elseif words[ind] == "'INTEND'"
                isFloat = true
            else
                error("Unknown MARKER in MPS file")
            end
        else
            name = words[1]
            if !haskey(varNames, name)
                varNames[name] = nVar
                nVar += 1
                A = [A initVar]
                push!(varTypes, isFloat)
                push!(c, 0.0)
            end

            for i in 2:2:lenWords
                if words[i] == objectiveName
                    # Case add to objective
                    c[varNames[name]] = parse(Float64, words[i+1])
                else
                    # Case add to A
                    A[conNames[words[i]], varNames[name]] = parse(Float64, words[i+1])
                end
            end
        end

        line = readline(mps)
        while strip(line) == "" || line[1] == '*'
            pos = position(mps)
            line = readline(mps)
        end
    end

    seek(mps, pos)
    return varNames, c, A, varTypes
end

function readRows(mps, objectiveName::String, fixed::Bool)
    pos = position(mps)
    objective = false
    conNames = Dict{String, Int}()
    conTypes = Vector{Int}() # 1 : <=, 2 : ==, 3 : >=
    nCon = 1

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end
    while line[1] == ' '
        pos = position(mps)
        if fixed
            typeTmp = strip(line[2:3])
            name = strip(line[5:min(length(line), 12)])
        else
            typeTmp, name = split(line)
            typeTmp = uppercase(typeTmp)
        end
        if typeTmp == "N"
            if objective
                error("MultiObjectives")
            else
                objective = true
                objectiveName = name
            end
        else
            if typeTmp == "G"
                push!(conTypes, 3)
            elseif typeTmp == "L"
                push!(conTypes, 1)
            elseif typeTmp == "E"
                push!(conTypes, 2)
            else
                error("Error in type of the row")
            end
            conNames[name] = nCon
            nCon += 1
        end

        line = readline(mps)
        while strip(line) == "" || line[1] == '*'
            pos = position(mps)
            line = readline(mps)
        end
    end

    seek(mps, pos)
    return objectiveName, conNames, conTypes
end

function readObjname(mps)
    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end

    seek(mps, position(mps))
    return strip(line)
end

function readObjsense(mps)
    objectiveSense = ["MIN", "MAX"]

    line = readline(mps)
    while strip(line) == "" || line[1] == '*'
        line = readline(mps)
    end

    if !(strip(line) in objectiveSense)
        error("No valid information in OBJSENSE section")
    end

    seek(mps, position(mps))
    return strip(line) == "MIN" ? :min : :max
end

#end # module

max_stars_repo_name: remi-garcia/MPSReader

max_issues_repo_name: remi-garcia/MPSReader

max_forks_repo_name: remi-garcia/MPSReader


text:
function render(filename::String = tempname();
                filetype=:PNG,
                animation=false, transparent=false)
    # https://blender.stackexchange.com/a/149715
    render_context.film_transparent = transparent
    transparent && (render_context.image_settings.color_mode = "RGBA")

    render_context.use_file_extension = true
    render_context.filepath = filename
    render_context.image_settings.file_format = string(filetype)
    render_ops.render(animation=animation,
                      write_still=!animation)
end


text:
using Documenter, GridLayoutBase

makedocs(;
    modules=[GridLayoutBase],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jkrumbiegel/GridLayoutBase.jl/blob/{commit}{path}#L{line}",
    sitename="GridLayoutBase.jl",
    authors="Julius Krumbiegel",
    assets=String[],
)

deploydocs(;
    repo="github.com/jkrumbiegel/GridLayoutBase.jl",
)


text:
export iauRm2v
"""
Express an r-matrix as an r-vector.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   r        double[3][3]    rotation matrix

Returned:
   w        double[3]       rotation vector (Note 1)

Notes:

   1. A rotation matrix describes a rotation through some angle about
   some arbitrary axis called the Euler axis.  The "rotation vector"
   returned by this function has the same direction as the Euler axis,
   and its magnitude is the angle in radians.  (The magnitude and
   direction can be separated by means of the function iauPn.)

   2. If r is null, so is the result.  If r is not a rotation matrix
   the result is undefined;  r must be proper (i.e. have a positive
   determinant) and real orthogonal (inverse = transpose).

   3. The reference frame rotates clockwise as seen looking along
   the rotation vector from the origin.

This revision:  2015 January 30

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauRm2v(double r[3][3], double w[3])

function iauRm2v(r::AbstractMatrix{<:Real})
   w = zeros(Float64, 3)

   ccall((:iauRm2v, libsofa_c), Cvoid,
         (Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Array{Float64, 2}, r'),
          w)

   return SVector{3}(w)
end


text:
# This file was generated, do not modify it. # hide
using PyPlot

figure(figsize=(8, 6))
for (m, c) in zip((mask1, mask2, mask3), ("red", "green", "blue"))
    plot(W[m, 1], W[m, 2], ls="none", marker=".", markersize=10, color=c)
end

xlabel("PCA-1", fontsize=13)
ylabel("PCA-2", fontsize=13)
legend(["Group 1", "Group 2", "Group 3"], fontsize=13)

savefig(joinpath(@OUTPUT, "ISL-lab-10-cluster.svg")) # hide


text:
"""
    ConvergentCrossMappingTest(timeseries_lengths;
        dim::Int = 3, τ::Int = 1, libsize::Int = 10,
        replace::Bool = false, n_reps::Int = 100, surr_func::Function = randomshuffle,
        which_is_surr::Symbol = :none, exclusion_radius::Int = 0,
        tree_type = NearestNeighbors.KDTree, distance_metric = Distances.Euclidean(),
        correspondence_measure = StatsBase.cor,
        η::Int = 0)

The parameters for a convergent cross mapping [1] test.

## Mandatory keyword arguments

- **`timeseries_lengths`**: The time series lengths over which to cross map and check convergence.

## Optional keyword arguments

- **`dim`**: The dimension of the state space reconstruction (delay embedding)
    constructed from the `response` series. Default is `dim = 3`.
- **`τ`**: The embedding lag for the delay embedding constructed from `response`.
    Default is `τ = 1`.
- **`η`**: The prediction lag to use when predicting scalar values of `driver`
    fromthe delay embedding of `response`.
    `η > 0` are forward lags (causal; `driver`'s past influences `response`'s future),
    and `η < 0` are backwards lags (non-causal; `driver`'s' future influences
    `response`'s past). Adjust the prediction lag if you
    want to performed lagged ccm
    [(Ye et al., 2015)](https://www.nature.com/articles/srep14750).
    Default is `η = 0`, as in
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079).
    *Note: The sign of the lag `η` is organized to conform with the conventions in
    [TransferEntropy.jl](), and is opposite to the convention used in the
    [`rEDM`](https://cran.r-project.org/web/packages/rEDM/index.html) package
    ([Ye et al., 2016](https://cran.r-project.org/web/packages/rEDM/index.html)).*
- **`libsize`**: Among how many delay embedding points should we sample time indices
    and look for nearest neighbours at each cross mapping realization (of which there
    are `n_reps`)?
- **`n_reps`**: The number of times we draw a library of `libsize` points from the
    delay embedding of `response` and try to predict `driver` values. Equivalently,
    how many times do we cross map for this value of `libsize`?
    Default is `n_reps = 100`.
- **`replace`**: Sample delay embedding points with replacement? Default is `replace = true`.
- **`exclusion_radius`**: How many temporal neighbors of the delay embedding
    point `response_embedding(t)` to exclude when searching for neighbors to
    determine weights for predicting the scalar point `driver(t + η)`.
    Default is `exclusion_radius = 0`.
- **`which_is_surr`**: Which data series should be replaced by a surrogate
    realization of the type given by `surr_type`? Must be one of the
    following: `:response`, `:driver`, `:none`, `:both`.
    Default is `:none`.
- **`surr_func`**: A valid surrogate function from TimeseriesSurrogates.jl.
- **`tree_type`**: The type of tree to build when looking for nearest neighbors.
    Must be a tree type from NearestNeighbors.jl. For now, this is either
    `BruteTree`, `KDTree` or `BallTree`.
- **`distance_metric`**: An instance of a `Metric` from Distances.jl. `BallTree` and `BruteTree` work with any `Metric`.
    `KDTree` only works with the axis aligned metrics `Euclidean`, `Chebyshev`,
    `Minkowski` and `Cityblock`. Default is `metric = Euclidean()` *(note the instantiation of the metric)*.
- **`correspondence_measure`**: The function that computes the correspondence
    between actual values of `driver` and predicted values. Can be any
    function returning a similarity measure between two vectors of values.
    Default is `correspondence_measure = StatsBase.cor`, which returns values on ``[-1, 1]``.
    In this case, any negative values are usually filtered out (interpreted as zero coupling) and
    a value of ``1`` means perfect prediction.
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
    also proposes to use the root mean square deviation, for which a value of ``0`` would
    be perfect prediction.
- **`summarise`**: Should cross map skills be summarised for each time series length? Default is
    `summarise = false`.
- **`average_measure`**: Either `:median` or `:mean`. Default is `:median`
- **`uncertainty_measure`**: Either `:quantile` or `:std`.
- **`quantiles`**: Compute uncertainty over quantile(s) if `uncertainty_measure` is `:quantile`.
    Default is `[0.327, 0.673]`, roughly corresponding to `1s` for normally distributed data.

## References

1. Sugihara, George, et al. "Detecting causality in complex ecosystems." Science (2012): 1227079.
    [http://science.sciencemag.org/content/early/2012/09/19/science.1227079](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)

2. Ye, Hao, et al. "Distinguishing time-delayed causal interactions using convergent cross mapping." Scientific
    Reports 5 (2015): 14750. [https://www.nature.com/articles/srep14750](https://www.nature.com/articles/srep14750)

3. Ye, H., et al. "rEDM: Applications of empirical dynamic modeling from time series." R Package Version
    0.4 7 (2016). [https://cran.r-project.org/web/packages/rEDM/index.html](https://cran.r-project.org/web/packages/rEDM/index.html)
"""
Base.@kwdef mutable struct ConvergentCrossMappingTest{N, NL} <: DistanceBasedCausalityTest{N}
    """
    The dimension of the state space reconstruction (delay embedding) constructed
    from the response series. Default is `dim = 3`.
    """
    dim::Int = 3

    """ The embedding lag for the delay embedding constructed from response. Default is τ = 1. """
    τ::Int = 1

    """
    The prediction lag to use when predicting scalar values of `driver`
    from the delay embedding of `response`.
    `η > 0` are forward lags (causal; `driver`'s past influences `response`'s future),
    and `η < 0` are backwards lags (non-causal; `driver`'s' future influences
    `response`'s past). Adjust the prediction lag if you
    want to performed lagged ccm
    [(Ye et al., 2015)](https://www.nature.com/articles/srep14750).
    Default is `η = 0`, as in
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079).
    *Note: The sign of the lag `η` is organized to conform with the conventions in
    [TransferEntropy.jl](), and is opposite to the convention used in the
    [`rEDM`](https://cran.r-project.org/web/packages/rEDM/index.html) package
    ([Ye et al., 2016](https://cran.r-project.org/web/packages/rEDM/index.html)).*
    """
    η::Int = 0

    """
    Among how many delay embedding points should we sample time indices
    and look for nearest neighbours at each cross mapping realization (of which there
    are `n_reps`)?
    """
    libsize::Int = 10

    """ Sample delay embedding points with replacement? Default is `replace = true`. """
    replace::Bool = true

    """
    The number of times we draw a library of `libsize` points from the
    delay embedding of `response` and try to predict `driver` values. Equivalently,
    how many times do we cross map for this value of `libsize`?
    Default is `n_reps = 100`.
    """
    n_reps::Int = 100

    """ A valid surrogate function from TimeseriesSurrogates.jl. """
    surr_func::Function = randomshuffle

    """
    Which data series should be replaced by a surrogate
    realization of the type given by `surr_type`? Must be one of the
    following: `:response`, `:driver`, `:none`, `:both`.
    Default is `:none`.
    """
    which_is_surr::Symbol = :none

    """
    How many temporal neighbors of the delay embedding
    point `response_embedding(t)` to exclude when searching for neighbors to
    determine weights for predicting the scalar point `driver(t + η)`.
    Default is `exclusion_radius = 0`.
    """
    exclusion_radius::Int = 0

    """
    The type of tree to build when looking for nearest neighbors.
    Must be a tree type from NearestNeighbors.jl. For now, this is either
    `BruteTree`, `KDTree` or `BallTree`.
    """
    tree_type = NearestNeighbors.KDTree

    """
    An instance of a `Metric` from Distances.jl. `BallTree` and `BruteTree` work with any `Metric`.
    `KDTree` only works with the axis aligned metrics `Euclidean`, `Chebyshev`,
    `Minkowski` and `Cityblock`. Default is `metric = Euclidean()`
    *(note the instantiation of the metric)*.
    """
    distance_metric = Distances.Euclidean()

    """
    The function that computes the correspondence between actual values of `driver` and predicted
    values. Can be any function returning a similarity measure between two vectors of values.
    Default is `correspondence_measure = StatsBase.cor`, which returns values on ``[-1, 1]``.
    In this case, any negative values are usually filtered out (interpreted as zero coupling) and
    a value of ``1`` means perfect prediction.
    [Sugihara et al. (2012)](http://science.sciencemag.org/content/early/2012/09/19/science.1227079)
    also proposes to use the root mean square deviation, for which a value of ``0`` would
    be perfect prediction.
    """
    correspondence_measure = StatsBase.cor

    """ Should cross map skills be summarised for each time series length? Default is `summarise = false`. """
    summarise::Bool = false

    """ Either `:median` or `:mean`. Default is `:median`. """
    average_measure = :median

    """ Either `:quantile` or `:std`. """
    uncertainty_measure = :quantile

    """
    Compute uncertainty over quantile(s) if `uncertainty_measure` is `:quantile`.
    Default is ``, roughly corresponding to `1s` for normally distributed data.
    """
    quantiles = [0.327, 0.673]

    """ The time series lengths over which to cross map and check convergence """
    timeseries_lengths

    function ConvergentCrossMappingTest(dim::Int, τ::Int, η::Int, libsize::Int, replace::Bool,
        n_reps::Int, surr_func::Function, which_is_surr::Symbol,
        exclusion_radius::Int, tree_type, distance_metric, correspondence_measure,
        summarise::Bool, average_measure, uncertainty_measure, quantiles,
        timeseries_lengths)

        # The number of returned elements. We're performing the test for `n_reps`
        # different randomly selected point libraries, so `N = n_reps`.
        N = n_reps

        # The number of different time series lengths
        NL = length(timeseries_lengths)

        new{N, NL}(dim, τ, η, libsize, replace, n_reps, surr_func, which_is_surr,
            exclusion_radius, tree_type, distance_metric, correspondence_measure,
            summarise, average_measure, uncertainty_measure, quantiles,
            timeseries_lengths)
    end
end

function causality(source::AbstractVector{T}, target::AbstractVector{T}, p::ConvergentCrossMappingTest) where {T<:Real}
     convergentcrossmap(source, target, p.timeseries_lengths;
        dim = p.dim,
        τ = p.τ,
        ν = p.η,
        libsize = p.libsize,
        n_reps = p.n_reps,
        replace = p.replace,
        exclusion_radius = p.exclusion_radius,
        jitter = 1e-5 * max(StatsBase.std(source), StatsBase.std(target)),
        which_is_surr = p.which_is_surr,
        surr_func = p.surr_func,
        tree_type = p.tree_type,
        distance_metric = p.distance_metric,
        correspondence_measure = p.correspondence_measure,
        summarise = p.summarise,
        quantiles = p.quantiles,
        uncertainty_measure = p.uncertainty_measure,
        average_measure = p.average_measure)
end

export ConvergentCrossMappingTest


text:
# Examples presented in class - Lecture 1
# Example 4 - Classification

using JuMP, Ipopt

# Generating random multivariate data.
data11 = rand(Normal(8, 4), 100)
data12 = rand(Normal(12, 2), 100)
data1 = [data11 data12]

data01 = rand(Normal(2, 6), 100)
data02 = rand(Normal(5, 2), 100)
data0 = [data01 data02]

scatter(data1[:,1], data1[:,2], label="Pos. obs.")
scatter!(data0[:,1], data0[:,2], label="Neg. obs.")

function classify(data1, data0, δ)
    m = Model(IpoptOptimizer())

    @variable(m, u[1:100] >= 0)
    @variable(m, v[1:100] >= 0)
    @variable(m, a[1:2])
    @variable(m, b)

    @constraint(m, error1[i=1:100], sum(a[j]*data1[i,j] for j = 1:2) - b - u[i] <= -1)
    @constraint(m, error0[i=1:100], sum(a[j]*data0[i,j] for j = 1:2) - b + v[i] >= 1)

# This parameter trades off accuracy and robustness for the classifier
    @objective(m, Min, sum(u[i] for i=1:100) + sum(v[i] for i=1:100)
        + δ*(sum(a[j]^2 for j=1:2)))

    solve(m);

    return getvalue(a), getvalue(b), getvalue(u), getvalue(v)
end

a, b, u1, v1 = classify(data1, data0, 0.01);

# Calculating the points on plane defined by a and b
x1 = linspace(-10,20,100);
x2 = (b - a[1]*x1)/ a[2];
# boundaries of the slab
x2u = (b + 1 - a[1]*x1)/ a[2];
x2d = (b - 1 - a[1]*x1)/ a[2];
plot!(x1, [x2u x2 x2d],
    xlabel = "x_1",
    ylabel = "x_2",
    label= [ "", "Classifier", ""] ,
    color=[:green :green :green],
    linestyle=[ :dot :solid :dot],
    )

# another classifier, with δ=50
a, b, u2, v2 = classify(data1, data0, 100);
x2 = (b - a[1]*x1)/ a[2];
x2u = (b + 1 - a[1]*x1)/ a[2];
x2d = (b - 1 - a[1]*x1)/ a[2];
plot!(x1, [x2u x2 x2d],
    xlabel = "x_1",
    ylabel = "x_2",
    legend = false,
    color=[:purple :purple :purple],
    linestyle=[:dot :solid :dot],
    )


text:
using Requires

struct D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion} <: Dirac_operator{Dim}   where  T <: AbstractGaugefields
    U::Array{T,1}
    wilsonoperator::Union{Wilson_Dirac_operator{Dim,T,wilsonfermion},Wilson_Dirac_operator_faster{Dim,T,wilsonfermion} }
    mass::Float64
    _temporary_fermi::Array{fermion,1}
    L5::Int64
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    _temporary_fermion_forCG::Vector{fermion}
    boundarycondition::Vector{Int8}
    b::Float64 #coefficient for MobiusDomainwall
    c::Float64 #coefficient for MobiusDomainwall

end

function D5DW_MobiusDomainwall_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters,mass,b,c) where  {NC,Dim}
    @assert haskey(parameters,"L5") "parameters should have the keyword L5"
    L5 = parameters["L5"]
    if L5 != x.L5
        @warn "L5 in Dirac operator and fermion fields is not same. Now L5 = $L5 and x.L5 = $(x.L5)"
        #@assert L5 == x.L5 "L5 in Dirac operator and fermion fields should be same. Now L5 = $L5 and x.L5 = $(x.L5)"
    end

    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    #boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    T = eltype(U)

    #Parameters for Wilson operator
        r = 1
        M = check_parameters(parameters,"M",-1)
        #@assert haskey(parameters,"M") "parameters should have the keyword M"
        #M = parameters["M"]
        κ_wilson = 1/(2*Dim*r+2M)
        parameters_wilson = Dict()
        parameters_wilson["κ"] = κ_wilson
        parameters_wilson["numtempvec"] = 4
        parameters_wilson["numtempvec_CG"] = 1
        parameters_wilson["boundarycondition"] = boundarycondition
        x_wilson  = x.w[1]
        fasterversion = check_parameters(parameters,"faster version",false)
        if fasterversion
            wilsonoperator = Wilson_Dirac_operator_faster(U,x_wilson,parameters_wilson)
        else
            wilsonoperator = Wilson_Dirac_operator(U,x_wilson,parameters_wilson)
        end
    #--------------------------------

    xtype = typeof(x)
    num = 4
    _temporary_fermi = Array{xtype,1}(undef,num)
    for i=1:num
        _temporary_fermi[i] = similar(x)
    end

    numcg = 7
    _temporary_fermion_forCG= Array{xtype,1}(undef,numcg)
    for i=1:numcg
        _temporary_fermion_forCG[i] = similar(x)
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")
    #println("xtype ",xtype)
    #println("D ", typeof(wilsonoperator))

    return D5DW_MobiusDomainwall_operator{Dim,T,xtype,typeof(x_wilson)}(
        U,
        wilsonoperator,
        mass,
        _temporary_fermi,
        L5,
        eps_CG,
        MaxCGstep,
        verbose_level,
        method_CG,
        verbose_print,
        _temporary_fermion_forCG,
        boundarycondition,
        b,
        c
        )
end

function (D::D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion})(U) where {Dim,T,fermion,wilsonfermion}
    return D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion}(
        U,
        D.wilsonoperator(U),
        D.mass,
        D._temporary_fermi,
        D.L5,
        D.eps_CG,
        D.MaxCGstep,
        D.verbose_level,
        D.method_CG,
        D.verbose_print,
        D._temporary_fermion_forCG,
        D.boundarycondition,
        D.b,
        D.c
    )
end

struct Adjoint_D5DW_MobiusDomainwall_operator{T} <: Adjoint_Dirac_operator
    parent::T
end

struct MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion} <: Dirac_operator{Dim}  where T <: AbstractGaugefields
    U::Array{T,1}
    D5DW::D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion}
    D5DW_PV::D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion}
    mass::Float64
    eps_CG::Float64
    MaxCGstep::Int64
    verbose_level::Int8
    method_CG::String
    verbose_print::Verbose_print
    boundarycondition::Vector{Int8}
    b::Float64 #coefficient for MobiusDomainwall
    c::Float64 #coefficient for MobiusDomainwall

end

function MobiusDomainwall_Dirac_operator(U::Array{<:AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
    @assert haskey(parameters,"mass") "parameters should have the keyword mass"
    mass = parameters["mass"]

    b = check_parameters(parameters,"b",2)
    c = check_parameters(parameters,"c",1)

    if b==1 && c ==1
        println_verbose_level1(U[1],"Shamir kernel (standard DW) is used")
    elseif b==2 && c== 0
        println_verbose_level1(U[1],"Borici/Wilson kernel (truncated overlap) is used")
    elseif b == 2&& c == 1
        println_verbose_level1(U[1],"scaled Shamir kernel (Mobius DW) is used")
    end

    D5DW = D5DW_MobiusDomainwall_operator(U,x,parameters,mass,b,c)
    D5DW_PV = D5DW_MobiusDomainwall_operator(U,x,parameters,1,b,c)

    #boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])

    if Dim == 4
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,1,1,-1])
    elseif Dim == 2
        boundarycondition = check_parameters(parameters,"boundarycondition",[1,-1])
    else
        error("Dim should be 2 or 4!")
    end

    eps_CG = check_parameters(parameters,"eps_CG",default_eps_CG)
    #println("eps_CG = ",eps_CG)
    MaxCGstep = check_parameters(parameters,"MaxCGstep",default_MaxCGstep)

    verbose_level = check_parameters(parameters,"verbose_level",2)
    verbose_print = Verbose_print(verbose_level)

    method_CG = check_parameters(parameters,"method_CG","bicg")

    return MobiusDomainwall_Dirac_operator{Dim,eltype(U),typeof(x),typeof(x.w[1])}(
            U,
            D5DW,D5DW_PV,mass,
            eps_CG,
            MaxCGstep,
            verbose_level,
            method_CG,
            verbose_print,
            boundarycondition,
            b,
            c
    )
end

function (D::MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion})(U) where {Dim,T,fermion,wilsonfermion}
    return MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}(
        U,
        D.D5DW(U),D.D5DW_PV(U),
        D.mass,
        D.eps_CG,
        D.MaxCGstep,
        D.verbose_level,
        D.method_CG,
        D.verbose_print,
        D.boundarycondition,
        D.b,
        D.c
    )
end

function get_temporaryvectors(A::T) where T <: MobiusDomainwall_Dirac_operator
    return A.D5DW._temporary_fermi
end

function get_temporaryvectors_forCG(A::T) where T <: MobiusDomainwall_Dirac_operator
    return A.D5DW._temporary_fermion_forCG
end

struct Adjoint_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion} <: Adjoint_Dirac_operator
    parent::MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}
end

function Base.adjoint(A::D5DW_MobiusDomainwall_operator)
    Adjoint_D5DW_MobiusDomainwall_operator(A)
end

function Base.adjoint(A::MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
    Adjoint_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion}(A)
end

struct DdagD_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion} <: DdagD_operator
    dirac::Domainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}
    function DdagD_MobiusDomainwall_operator(U::Array{<: AbstractGaugefields{NC,Dim},1},x,parameters) where  {NC,Dim}
        return new{Dim,eltype(U),typeof(x),typeof(x.w[1])}(MobiusDomainwall_Dirac_operator(U,x,parameters))
    end

    function DdagD_MobiusDomainwall_operator(D::MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        return new{Dim,T,fermion,wilsonfermion}(D)
    end
end

struct MobiusD5DWdagD5DW_Wilson_operator{T} <: DdagD_operator
    dirac::T
    function MobiusD5DWdagD5DW_Wilson_operator(U::Array{T,1},x,parameters,mass,b,c) where  T <: AbstractGaugefields
        return new{T}(D5DW_MobiusDomainwall_operator(U,x,parameters,mass,b,c))
    end

    function MobiusD5DWdagD5DW_Wilson_operator(D::MobiusDomainwall_Dirac_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        dtype = typeof(D.D5DW)
        return new{dtype}(D.D5DW)
    end

    function MobiusD5DWdagD5DW_Wilson_operator(D::D5DW_MobiusDomainwall_operator{Dim,T,fermion,wilsonfermion}) where {Dim,T,fermion,wilsonfermion}
        dtype = typeof(D)
        return new{dtype}(D)
    end
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: MobiusD5DWdagD5DW_Wilson_operator, T3 <:AbstractFermionfields}
    temp = A.dirac._temporary_fermi[1]
    mul!(temp,A.dirac,x)
    mul!(y,A.dirac',temp)

    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: D5DW_MobiusDomainwall_operator, T3 <:AbstractFermionfields}
    D5DWx!(y,A.U,x,A.mass,A.wilsonoperator,A.L5,A.b,A.c,A._temporary_fermi[3],A._temporary_fermi[4])
    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_D5DW_MobiusDomainwall_operator, T3 <:AbstractFermionfields}
    D5DWdagx!(y,A.parent.U,x,A.parent.mass,A.parent.wilsonoperator,A.parent.L5,A.parent.b,A.parent.c,
    A.parent._temporary_fermi[3],A.parent._temporary_fermi[4])
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: Adjoint_MobiusDomainwall_operator, T3 <:AbstractFermionfields}

    #A = D5DW(m)*D5DW(m=1)^(-1)
    #=
    A^+ = [D5DW(m)*D5DW(m=1)^(-1)]^+
        = D5DW(m=1)^(-1)^+ D5DW(m)^+
        = [D5DW(m=1)^+]^(-1) D5DW(m)^+
    y = A^+*x = [D5DW(m=1)^+]^(-1) D5DW(m)^+*x
    =#
    mul!(A.parent.D5DW_PV._temporary_fermi[1],A.parent.D5DW',x)
    bicg(y,A.parent.D5DW_PV',A.parent.D5DW_PV._temporary_fermi[1],
        eps=A.parent.eps_CG,maxsteps = A.parent.MaxCGstep,verbose = A.parent.verbose_print)
        #maxsteps = 10000)
        #verbose = Verbose_3())

    return
end

function LinearAlgebra.mul!(y::T1,A::T2,x::T3) where {T1 <:AbstractFermionfields,T2 <: MobiusDomainwall_Dirac_operator, T3 <:AbstractFermionfields}
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #y = A*x = D5DW(m)*D5DW(m=1))^(-1)*x
    #println("x ",x.w[1][1,1,1,1])
    bicg(A.D5DW_PV._temporary_fermi[1],A.D5DW_PV,x,
        eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
    mul!(y,A.D5DW,A.D5DW_PV._temporary_fermi[1])
    #println("y ",y.w[1][1,1,1,1])

    #bicg(A.D5DW._temporary_fermi[1],A.D5DW,y,
    #eps=A.eps_CG,maxsteps = A.MaxCGstep,verbose = A.verbose_print)
    #mul!(A.D5DW._temporary_fermi[2],A.D5DW_PV,A.D5DW._temporary_fermi[1])
    #println("tmp2 ",A.D5DW._temporary_fermi[2].w[1][1,1,1,1])

    #error("Do not use Domainwall_operator directory. Use D5DW_Domainwall_operator M = D5DW(m)*D5DW(-1)^{-1}")
    #D5DWx!(y,A.U,x,A.m,A.wilsonoperator._temporary_fermi)
    return
end

include("./MobiusDomainwallFermion_5d.jl")

function Initialize_MobiusDomainwallFermion(u::AbstractGaugefields{NC,Dim},L5;nowing=false) where {NC,Dim}
    _,_,NN... = size(u)
    return Initialize_MobiusDomainwallFermion(L5,NC,NN...,nowing = nowing)
end

function Initialize_MobiusDomainwallFermion(L5,NC,NN...;nowing = false)
    Dim = length(NN)
    if Dim == 4
        if nowing
            fermion = MobiusDomainwallFermion_5D(L5,NC,NN...,nowing = nowing)
        else
            fermion = MobiusDomainwallFermion_5D_wing(L5,NC,NN...)
        end
        #fermion = DomainwallFermion_5D_wing(L5,NC,NN...)
        #fermion = WilsonFermion_4D_wing{NC}(NN...)
        #fermion = WilsonFermion_4D_wing(NC,NN...)
    elseif Dim == 2
        fermion = MobiusDomainwallFermion_3D_wing(L5,NC,NN...)
        #fermion = WilsonFermion_2D_wing{NC}(NN...)
    else
        error("Dimension $Dim is not supported")
    end
    return fermion
end

function bicg(x,A::MobiusDomainwall_Dirac_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #A^-1 = D5DW(m=1)*DsDW(m)^-1
    #x = A^-1*b = D5DW(m=1)*DsDW(m)^-1*b
    #println("b ",b.w[1][1,1,1,1])
    bicg(A.D5DW._temporary_fermi[1],A.D5DW,b;eps=eps,maxsteps = maxsteps,verbose = verbose)

    #mul!(x,A.D5DW,A.D5DW._temporary_fermi[1])
    #println("x ",x.w[1][1,1,1,1])

    mul!(x,A.D5DW_PV,A.D5DW._temporary_fermi[1])
end

function bicg(x,A::Adjoint_MobiusDomainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) where T #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1)^(-1)
    #A' = (D5DW(m=1)^+)^(-1) D5DW(m)^+
    #A'^-1 = D5DW(m)^+^-1 D5DW(m=1)^+
    #x = A'^-1*b =D5DW(m)^+^-1 D5DW(m=1)^+*b
    mul!(A.parent.D5DW._temporary_fermi[1],A.parent.D5DW_PV',b)
    bicg(x,A.parent.D5DW',A.parent.D5DW._temporary_fermi[1];eps=eps,maxsteps = maxsteps,verbose = verbose)

end

function bicgstab(x,A::MobiusDomainwall_Dirac_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1))^(-1)
    #A^-1 = D5DW(m=1)*DsDW(m)^-1
    #x = A^-1*b = D5DW(m=1)*DsDW(m)^-1*b

    bicgstab(A.D5DW._temporary_fermi[1],A.D5DW,b;eps=eps,maxsteps = maxsteps,verbose = verbose)
    mul!(x,A.D5DW_PV,A.D5DW._temporary_fermi[1])
end

function bicgstab(x,A::Adjoint_MobiusDomainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2)) where T #A*x = b -> x = A^-1*b
    #A = D5DW(m)*D5DW(m=1)^(-1)
    #A' = (D5DW(m=1)^+)^(-1) D5DW(m)^+
    #A'^-1 = D5DW(m)^+^-1 D5DW(m=1)^+
    #x = A'^-1*b =D5DW(m)^+^-1 D5DW(m=1)^+*b
    mul!(A.parent.D5DW._temporary_fermi[1],A.parent.D5DW_PV',b)
    bicgstab(x,A.parent.D5DW',A.parent.D5DW._temporary_fermi[1];eps=eps,maxsteps = maxsteps,verbose = verbose)

end

function cg(x,A::DdagD_MobiusDomainwall_operator,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_print(2))
    #=
    A^-1 = ( (D5DW(m)*D5DW(m=1))^(-1))^+ D5DW(m)*D5DW(m=1))^(-1) )^-1
      = ( D5DW(m=1)^+)^(-1) D5DW(m)^+ D5DW(m)*D5DW(m=1))^(-1) )^-1
      = D5DW(m=1) (  D5DW(m)^+ D5DW(m) )^(-1) )^-1 D5DW(m=1)^+
    x = A^-1*b = D5DW(m=1) (  D5DW(m)^+ D5DW(m) )^(-1)  D5DW(m=1)^+*b
    =#
    mul!(A.dirac.D5DW_PV._temporary_fermi[1],A.dirac.D5DW_PV',b) #D5DW(m=1)^+*b

    temp = A.dirac.D5DW_PV._temporary_fermi[1]
    cg(A.dirac.D5DW._temporary_fermi[1],A.DdagD,temp;eps=eps,maxsteps = maxsteps,verbose = verbose) #(  D5DW(m)^+ D5DW(m) )^(-1)  D5DW(m=1)^+*b
    mul!(x,A.dirac.D5DW_PV,A.dirac.D5DW._temporary_fermi[1])
end


text:
## merge substitutions
#include("reso.jl")

# this will be used in proving with graph
# not yet implemented

"""
this merge may be a unfiy
"""
function merge(vars::Vlist, t1::Symbol, t2::Symbol)
 if t1==t2; return t1 end
 if isvar(t1,vars) && isvar(t2,vars)
   return t2
 elseif isvar(t1,vars)
   return t2
 elseif isvar(t2,vars)
   return t1
 else
   throw(ICMP(t1,t2,:merge))
 end
end

"""
merge unifys 2 to 1 for substituions
"""
function merge(vars::Vlist, sub1::Tlist, sub2::Tlist)
 if length(sub1) != length(sub2); throw(ICMP(sub1,sub2,:merge)) end
 s1=apply(vars,sub1,sub2) # is s1 == sub1??
 ss=apply(vars,sub2,s1)
 s2=apply(vars,sub1,ss)
 narg=[]
 for i in 1:length(sub1)
   mi = merge(vars, ss[i], s2[i])
   if ndims(sub1) == 1
     narg = vcat(narg,[mi])
   else
     narg = hcat(narg,[mi])
   end
 end
 narg
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
@testset "StraightLineHomotopy" begin
    @polyvar x y z
    f = x + y + 2z
    g = x + z
    h = y + z

    @test StraightLineHomotopy(f, g) isa StraightLineHomotopy{Float64}
    @test StraightLineHomotopy([f, f], [g, h]) isa StraightLineHomotopy{Float64}
    @test StraightLineHomotopy{Float64}([f, g], [f, h]) isa StraightLineHomotopy{Float64}
    @test StraightLineHomotopy{Float64}(f, g) isa StraightLineHomotopy{Float64}
    @test StraightLineHomotopy([FP.Polynomial(f)], [FP.Polynomial{Complex128}(f)]) isa StraightLineHomotopy{Complex128}
    @test_throws AssertionError StraightLineHomotopy{Complex128}([f], [g])
    @test_throws AssertionError StraightLineHomotopy{Complex128}([f], [f, f])

    @test StraightLineHomotopy{BigFloat}(FP.Polynomial{Float64}.([f, f]), convert(Vector{FP.Polynomial{Float64}}, [g, h])) isa StraightLineHomotopy{BigFloat}

    H = StraightLineHomotopy{Float64}([f, g], [f, h])
    K = StraightLineHomotopy([f, g], [f, h])
    @test promote_type(typeof(H), typeof(K)) == StraightLineHomotopy{Float64}
    @test convert(StraightLineHomotopy{Float64}, K) isa StraightLineHomotopy{Float64}

    w = [1.0, 2.0, 2.0]
    @test evaluate(H, w, 1.0) == [7, 3]
    @test H(w, 1.0) == [7, 3]
    u = zeros(2)
    evaluate!(u, H, w, 1.0)
    @test u == [7, 3]

    cfg = PolynomialHomotopyConfig(H)

    @test evaluate(H, w, 1.0, cfg) == [7, 3]
    u = zeros(2)
    evaluate!(u, H, w, 1.0, cfg)
    @test u == [7, 3]

    @test promote_type(StraightLineHomotopy{Float64}, Complex128) == StraightLineHomotopy{Complex128}
    @test promote_type(StraightLineHomotopy{Float64}, StraightLineHomotopy{Complex128}) == StraightLineHomotopy{Complex128}

    @test PolynomialHomotopyConfig(H, rand(Complex128, 3)) isa PolynomialHomotopyConfig{Complex128}

    @test jacobian(H, rand(3), 0.0, cfg) == [1 1 2; 0 1 1]
    @test jacobian(H, rand(3), 1.0, cfg) == [1 1 2; 1 0 1]
    u = zeros(2, 3)
    @test jacobian!(u, H, rand(3), 0.0, cfg) == [1 1 2; 0 1 1]
    @test u == [1 1 2; 0 1 1]

    r = JacobianDiffResult(cfg)
    w = rand(3)
    jacobian!(r, H, w, 0.24, cfg)
    @test value(r) == H(w, 0.24)
    @test jacobian(r) == jacobian(H, w, 0.24, cfg)

    @test dt(H, [1.0, 2, 2], 1.0, cfg) == [0, -1]
    u = zeros(2)
    dt!(u, H, [1, 2, 2.0], 1.0, cfg)
    @test u == [0, -1]
    @test string(H) == "Homotopies.StraightLineHomotopy{Float64} with 2 polynomials in 3 variables.\n"

    r = DtDiffResult(cfg)
    w = rand(3)
    dt!(r, H, w, 0.24, cfg)
    @test value(r) == H(w, 0.24)
    @test dt(r) == dt(H, w, 0.24, cfg)

    N = weylnorm(H)
    @test N(0.0)==sqrt(8)
    @test N(0.5)==sqrt(7.5)

    H = StraightLineHomotopy(x^2+y+z, z^2+2+x+y^2)
    @test ishomogenous(H) == false
    @test ishomogenized(H) == false
    HH = homogenize(H)
    @test ishomogenous(HH)
    @test ishomogenized(HH)
    @test H != HH
    @test HH == homogenize(HH)
    @test isequal(HH, homogenize(HH))
    @test dehomogenize(HH) == H

    @test length(H) == 1
    @test nvariables(H) == 3
end