text:
rng = MersenneTwister(912301)

@testset "API" begin
    # * Single-run tests for TestFunctions, Easom
    @test begin
        assert_results = []
        ground_truth = @SVector [π, π]
        for k = 1:10
            val = optimize(Easom(), zeros(2), [-100.0, 100], PSO(); rng=rng)
            result = isapprox(val.x, ground_truth)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Change the keyword arguments
    @test begin
        assert_results = []
        ground_truth = @SVector [π, π]
        args = Dict(:w => 0.8, :c1 => 1.5, :c2 => 1.0)
        for k = 1:10
            val = optimize(Easom(), zeros(2), [-100.0, 100], PSO(); rng=rng, args...)
            result = isapprox(val.x, ground_truth)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Test random RNG
    @test begin
        assert_results = []
        ground_truth = @SVector [π, π]
        for k = 1:10
            val = optimize(Easom(), zeros(2), [-100.0, 100], PSO(); iters=15_000)
            result = isapprox(val.x, ground_truth)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Test fixed RNG, SimulatedAnnealing
    @test begin
        assert_results = []
        ground_truth = 0.0
        for k = 1:10
            val = optimize(Rosenbrock(), zeros(2), [-5.0, 5.0], SimulatedAnnealing();
                rng=rng, iters=20_000)
            result = isapprox(val.min, ground_truth, atol=1e-2)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Test fixed RNG, GeneralSimulatedAnnealing
    @test begin
        assert_results = []
        ground_truth = 0.0
        for k = 1:10
            val = optimize(
                Rosenbrock(),
                zeros(2),
                [-5.0, 5.0],
                GeneralSimulatedAnnealing();
                rng=rng,
                iters=20_000
            )
            result = isapprox(val.min, ground_truth, atol=1e-2)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Test randomly seeded RNG, SimulatedAnnealing
    @test begin
        assert_results = []
        ground_truth = 0.0
        for k = 1:10
            val = optimize(Rosenbrock(), zeros(2), [-5.0, 5.0], SimulatedAnnealing();
                iters=20_000)
            result = isapprox(val.min, ground_truth, atol=1e-2)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Test randomly seeded RNG, GeneralSimulatedAnnealing
    @test begin
        assert_results = []
        ground_truth = 0.0
        for k = 1:10
            val = optimize(
                Rosenbrock(),
                zeros(2),
                [-5.0, 5.0],
                GeneralSimulatedAnnealing();
                iters=20_000
            )
            result = isapprox(val.min, ground_truth, atol=1e-2)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end
end


text:
@testset "CoreGeometry" begin

@test 1==1

end # testset

ext: jl


text:
module SnoopCompileCore

using Core: MethodInstance, CodeInfo

# @snoopi and @snoopc are exported from their files of definition

include("snoopc.jl")

if VERSION >= v"1.2.0-DEV.573"
    include("snoopi.jl")
end

if VERSION >= v"1.6.0-DEV.154"
    include("snoopr.jl")
end

end


text:
"""
1-D Gradient

    gradient(a::Array, dx::Float64=1)
    :param f: Vector
    :param dx: stepsize

    :return g: gradient
"""
function gradient(f::Array{Float64, 1}, dx::Float64=1.)
    n = length(f)
    g = zeros(n)
    h = dx*(1:n)
    # Take forward differences on left and right edges
    g[1] = (f[2] - f[1])/(h[2]-h[1])
    g[n] = (f[n] - f[n-1])/(h[end]-h[end-1])

    # Take centered differences on interior points
    h = h[3:n]-h[1:(n-2)]
    g[2:(n-1)] = (f[3:n]-f[1:(n-2)])/h[1]

    return g
end

"""
2-D Gradient

    gradient2(a::Array, dx::Float64=1, dy::Float=1)
    :param a: matrix
    :param dx: stepsize
    :param dy: stepsize

    :return dxdu: derivatives along x
    :return dydv: derivatives along y
"""
function gradient2(a::Array, dx::Float64=1., dy::Float64=1.)
    m,n = size(a);
    dxdu = zeros(m,n);
    dydv = zeros(m,n);

    for i=1:m
        dxdu[i,:] = gradient(vec(a[i,:]), dx)
    end

    for i=1:n
        dydv[:,i] = gradient(a[:,i], dy)
    end

    return dxdu, dydv
end

"""
Creates Rectangular Grid in 2-D space

    meshgrid(a::LinRange,b::LinRange)
"""
function meshgrid(a::LinRange,b::LinRange)
    grid_a = [i for i in a, j in b]';
    grid_b = [j for i in a, j in b]';

    return grid_a, grid_b
end

"""
Linear interpolation

    approx(xd, yd, xi)
    :param xd: x samples
    :param yd: response samples
    :param xi: new x samples
"""
function approx(xd, yd, xi)
    nd = length(xd);
    ni = length(xi);

    yi = zeros(ni);
    for i in 1:ni
        if (xi[i] <= xd[1])
            t = (xi[i]-xd[1]) / (xd[2] - xd[1]);
            yi[i] = (1.0 - t) * yd[1] + t * yd[2];
        elseif (xd[nd] <= xi[i])
            t = (xi[i] - xd[nd-1]) / (xd[nd] - xd[nd-1]);
            yi[i] = (1.0 - t) * yd[nd-1] + t * yd[nd];
        else
            for k in 2:nd
                if (xd[k-1] <= xi[i] && xi[i] <= xd[k])
                    t = (xi[i] - xd[k-1]) / (xd[k] - xd[k-1]);
                    yi[i] = (1.0 - t) * yd[k-1] + t * yd[k];
                    break
                end
            end
        end
    end

    return yi
end

"""
Trapezoidal Integration

    trapz(x, y, dim=1)
    :param x: vector of time samples
    :param y: array of response samples
    :param dim: dimension along which to integrate
"""
function trapz(x::Array{Float64, 1}, y::Array{Float64}, dim::Integer=1)
    perm = [dim:max(ndims(y),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        m = 1;
    else
        m = size(y,1);
    end

    if m == 1
        M = length(y);
        out = sum(diff(x).*(y[1:M-1] + y[2:M])/2.0);
    else
        out = transpose(diff(x)) * (reshape(y,Val(2))[1:m-1,:] + reshape(y,Val(2))[2:m,:])/2.0;
        siz = size(y); siz = collect(siz); siz[1] = 1;
        out = reshape(out, tuple(siz...));
        out = permutedims(out, invperm(perm));
        ind = findall(collect(size(out)).==1);
        out = dropdims(out,dims=ind[1]);
        if length(out) == 1;
            out = out[1];
        end
    end

    return out
end

"""
Cumulative Trapezoidal Integration

    cumtrapz(x, y, dim=1)
    :param x: vector describing time samples
    :param y: array describing response
    :param dim: dimension to integrate over
"""
function cumtrapz(x::Array{Float64, 1}, y::Array{Float64}, dim::Integer=1)
    perm = [dim:max(length(size(y)),dim); 1:dim-1];
    y = permutedims(y, perm);
    if ndims(y) == 1
        n = 1;
        m = length(y);
    else
        m, n = size(y);
    end

    if n == 1
        dt = diff(x)/2.0;
        z = [0; cumsum(dt.*(y[1:(m-1)] + y[2:m]))];
    else
        dt = repeat(diff(x)/2.0,1,n);
        z = [zeros(1,n); cumsum(dt.*(y[1:(m-1), :] + y[2:m, :]),dims=1)];
        z = permutedims(z, invperm(perm));
    end

    return z
end

"""
Cumulative Trapezoidal Integration using midpoint

    cumtrapzmid(x, y, c)
    :param x: time samples
    :param y: resposne samples
    :param c: midpoint
    :param mid: midpoint location
"""
function cumtrapzmid(x, y, c, mid)
    a = length(x);

    # case < mid
    fn = zeros(a);
    tmpx = x[(mid-1):-1:1];
    tmpy = y[(mid-1):-1:1];
    tmp = c .+ cumtrapz(tmpx, tmpy);
    fn[1:(mid-1)] = reverse(tmp);

    # case >= mid
    fn[mid:a] = c .+ cumtrapz(x[mid:a],y[mid:a]);

    return fn

end

"""
Multivariate Normal random number generation

    mvnrand(mu, C, n)
    :param mu: mean vector
    :param C: covariance matrix
    :param n: number of samples
"""
function mvnrand(mu, C, n)
    tmp = cholesky(C, :U, Val{true});
    R = tmp[:U];
    R = Array(R);
    R = R[:, tmp.piv];
    retval = randn(n, size(R,1)) * R;
    retval += transpose(repeat(mu, 1, n));
    return transpose(retval)
end

"""
Linear interpolation when response contains flat regions

    interp1_flat(x, y, xx)
    :param x: time samples
    :param y: response samples
    :param xx: new time samples
"""
function interp1_flat(x,y,xx)
    flat = find(diff(x).<=0);
    n = length(flat);

    if n==0
        tmp = interpolate((x,), y, Gridded(Linear()))
        yy = tmp[xx];
    else
        yy = zeros(size(xx));
        i1 = 1;
        if flat[1] == 1
            i2 = 1;
            j = xx.==x[i2];
            yy[j] = minimum(y[i2:i2+1]);
        else
            i2 = flat[1];
            j = (xx.>=x[i1]) & (xx.<=x[i2]);
            tmp = interpolate((x[i1:i2],), y[i1:i2], Gridded(Linear()))
            yy[j] = tmp[xx[j]];
            i1 = copy(i2);
        end
        for k = 2:n
            i2 = flat[k];
            if i2 > i1+1
                j = (xx.>=x[i1]) & (xx.<=x[i2]);
                yi = interpolate((x[i1+1:i2],), y[i1+1:i2], Gridded(Linear()))
                yy[j] = tmp[xx[j]];
            end
            j = xx.==x[i2];
            yy[j] = minimum(y[i2:i2+1]);
            i1 = copy(i2);
        end
        i2 = length(x);
        j = (xx.>=x[i1]) & (xx.<=x[i2]);
        if i1+1 == i2
            yy[j] = y[i2];
        else
            tmp = interpolate((x[i1+1:i2],), y[i1+1:i2], Gridded(Linear()))
            yy[j] = tmp[xx[j]];
        end
    end

    return yy
end

function basis_fourier(f_domain::Vector, numBasis::Int64, fourier_p::Int64)
    result = Array{Float64}(undef,length(f_domain),2*numBasis)
    for i in 1:2*numBasis
        j = ceil(i/2)
        if ((i % 2) == 1)
            result[:,i] = sqrt(2) * sin.(2*j*pi*f_domain/fourier_p)
        end
        if ((i % 2) == 0)
            result[:,i] = sqrt(2) * cos.(2*j*pi*f_domain/fourier_p)
        end
    end

    out = basis(f_domain,result)
    return out
end

function f_predictfunction(f::func, at; deriv=0, method="linear")
    if method=="linear"
        if (deriv==0)
            result = approx(f.x, f.y, at)
        end
        if (deriv == 1)
            fmod = approx(f.x, f.y, at)
            diffy1 = [0; diff(fmod)]
            diffy2 = [diff(fmod); 0]
            diffx1 = [0; diff(at)]
            diffx2 = [diff(at); 0]

            result = (diffy2 + diffy1) ./ (diffx2 + diffx1)
        end
    else
        error("Method not implemented")
    end

    out = func(at, result)
    return out
end

function f_basistofunction(f_domain, coef, cbasis::basis; coefconst=0)
    if (size(cbasis.matrix,2)<length(coef))
        error("coefficients exceeds basis functions")
    end

    y1 = cbasis.matrix[:,(1:length(coef))] * coef .+ coefconst
    result = func(cbasis.x, y1)
    y1 = f_predictfunction(result, f_domain)

    return y1
end


text:
module ARLToolkit

using Requires

export AIS, CTD, Logs, Bathymetry, BLAS, GPXFile
export UnetLogs, bathy

include("utils.jl")
include("AIS.jl")
include("CTD.jl")
include("Logs.jl")
include("Bathymetry.jl")
include("BLAS.jl")

import .Logs: UnetLogs
import .Bathymetry: bathy

function __init__()
  @eval Base.show(io::IO, x::LLA) = prettyprint(io, x)
  @eval Base.print(io::IO, x::ZonedDateTime) = prettyprint(io, x)
  @require GPX="b55ef746-885f-40a4-ab22-c8118be08013" include("GPX.jl")
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    include("plot-recipes.jl")
    include("TimingDiagrams.jl")
  end
  @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("makie-recipes.jl")
end

end # module


text:
if !any([dirname(path) == dirname(@__FILE__) for path in LOAD_PATH])
    push!(LOAD_PATH, dirname(@__FILE__))
end

ENV["LINES"] = 15

using RBT: Entity, ToyDatabase, Schema, Class, Arrow, Instance, Iso, Opt, Seq

CSV_PATH = joinpath(dirname(@__FILE__), "../data/hr/hr.csv")

TREASURER_DEPT = "TREASURER"
TREASURER_HEAD = "CITY TREASURER"
ACC_HEAD = "DIR OF ACCOUNTING"
ACCS = ["AUDITOR IV", "ACCOUNTANT IV", "ACCOUNTANT III", "STAFF ASST", "ACCOUNTANT I"]

schema = Schema(
    Class(
        :department,
        Arrow(:name, String, runique=true),
        Arrow(
            :employee,
            lunique=false, ltotal=false, runique=true, rtotal=true,
            select=(:name, :surname), inverse=:department),
        select=(:name,)),
    Class(
        :employee,
        Arrow(:name, String),
        Arrow(:surname, String),
        Arrow(:position, String),
        Arrow(:salary, Int),
        Arrow(:department, select=:name),
        Arrow(
            :managed_by, :employee,
            ltotal=false,
            select=(:name, :surname, :position)),
        Arrow(
            :manages, :employee,
            lunique=false, ltotal=false, runique=true,
            select=(:name, :surname, :position),
            inverse=:managed_by),
        select=(:name, :surname, :department, :position, :salary)))

D = Entity{:department}
E = Entity{:employee}

dept = Vector{D}()
empl = Vector{E}()

dept_name = Vector{Iso{String}}()
dept_empl = Vector{Seq{E}}()

empl_name = Vector{Iso{String}}()
empl_surname = Vector{Iso{String}}()
empl_position = Vector{Iso{String}}()
empl_salary = Vector{Iso{Int}}()
empl_dept = Vector{Iso{D}}()
empl_managed_by = Vector{Opt{E}}()
empl_manages = Vector{Seq{E}}()

tr_head = nothing
trs = []
acc_head = nothing
accs = []

csv = readcsv(CSV_PATH, header=true)[1]
did_by_name = Dict{String, D}()
for i = 1:size(csv, 1)
    department_name, name, surname, position, salary = csv[i, :]
    if !haskey(did_by_name, department_name)
        did = D(length(dept)+1)
        push!(dept, did)
        push!(dept_name, Iso{String}(department_name))
        push!(dept_empl, Seq{E}(E[]))
        did_by_name[department_name] = did
    end
    did = did_by_name[department_name]
    eid = E(length(empl)+1)
    push!(empl, eid)
    push!(dept_empl[did.id].data, eid)
    push!(empl_name, Iso{String}(name))
    push!(empl_surname, Iso{String}(surname))
    push!(empl_position, Iso{String}(position))
    push!(empl_salary, Iso{Int}(salary))
    push!(empl_dept, Iso{D}(did))
    if department_name == TREASURER_DEPT
        if position == TREASURER_HEAD
            tr_head = eid
        elseif position == ACC_HEAD
            acc_head = eid
            push!(trs, eid)
        elseif position in ACCS
            push!(accs, eid)
        else
            push!(trs, eid)
        end
    end
    push!(empl_manages, Seq{E}(E[]))
    push!(empl_managed_by, Opt{E}())
end

empl_manages[tr_head.id] = Seq{E}(trs)
empl_manages[acc_head.id] = Seq{E}(accs)
for tr in trs
    empl_managed_by[tr.id] = Opt{E}(tr_head)
end
for acc in accs
    empl_managed_by[acc.id] = Opt{E}(acc_head)
end

instance = Instance(
    Dict(
        :department => dept,
        :employee => empl),
    Dict(
        (:department, :name) => dept_name,
        (:department, :employee) => dept_empl,
        (:employee, :name) => empl_name,
        (:employee, :surname) => empl_surname,
        (:employee, :position) => empl_position,
        (:employee, :salary) => empl_salary,
        (:employee, :department) => empl_dept,
        (:employee, :managed_by) => empl_managed_by,
        (:employee, :manages) => empl_manages))

citydb = ToyDatabase(schema, instance)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    show_states_initial_value(sim::Simulation)

Function to print initial states.

# Arguments

- `sim::Simulation` : Simulation object that contains the initial condition
"""
function show_states_initial_value(sim::Simulation)
    inputs = sim.inputs
    bus_size = get_bus_count(inputs)
    system = get_system(sim)
    global_state_map = make_global_state_map(inputs)
    println("Voltage Variables")
    println("====================")
    buses_sorted =
        sort(collect(PSY.get_components(PSY.Bus, system)); by = x -> PSY.get_number(x))
    for bus in buses_sorted
        name = PSY.get_name(bus)
        println(name)
        println("====================")
        bus_n = PSY.get_number(bus)
        bus_ix = PSID.get_lookup(sim.inputs)[bus_n]
        V_R = sim.x0_init[bus_ix]
        V_I = sim.x0_init[bus_ix + bus_size]
        Vm = sqrt(V_R^2 + V_I^2)
        θ = angle(V_R + V_I * 1im)
        print("Vm ", round(Vm, digits = 4), "\n")
        print("θ ", round(θ, digits = 4), "\n")
        println("====================")
    end
    println("====================")
    for device in get_dynamic_injectors(inputs)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        println("Differential States")
        println(name)
        println("====================")
        global_index = global_state_map[name]
        for s in states
            print(s, " ", round(sim.x0_init[global_index[s]], digits = 4), "\n")
        end
        println("====================")
    end
    dyn_branches = get_dynamic_branches(inputs)
    if !isempty(dyn_branches)
        println("====================")
        println("Line Current States")
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            printed_name = "Line " * name
            println("====================")
            println(printed_name)
            global_index = global_state_map[name]
            x0_br = Dict{Symbol, Float64}()
            for (i, s) in enumerate(states)
                print(s, " ", round(sim.x0_init[global_index[s]], digits = 5), "\n")
            end
            println("====================")
        end
    end
    return
end


text:
module SlicePass
using Test
using Cassette
import Cassette: @context, canoverdub, overdub, fallback, similarcontext
print("   running SliceCtx test...")
 before_time = time()
 using Core: CodeInfo, SlotNumber, SSAValue
 @context SliceCtx
 function Cassette.execute(ctx::SliceCtx, callback, f, args...)
    if canoverdub(ctx, f, args...)
        _ctx = similarcontext(ctx, metadata = callback)
        return overdub(_ctx, f, args...) # return result, callback
    else
        return fallback(ctx, f, args...), callback
    end
end
 const global_test_cache = Any[]
 push_to_global_test_cache!(x) = push!(global_test_cache, x)
 function Cassette.execute(ctx::SliceCtx, callback, ::typeof(push_to_global_test_cache!), x)
    return nothing, () -> (callback(); push_to_global_test_cache!(x))
end
 # handle Core._apply calls; Cassette might do this for you in a future update
function Cassette.execute(ctx::SliceCtx, callback, ::typeof(Core._apply), f, args...)
    return Core._apply(Cassette.execute, (ctx,), (callback,), (f,), args...)
end

function varname(ir::CodeInfo, sym::Symbol)
    s = string(sym)[2:end]
    i = parse(Int,s)
    varname = ir.slotnames[i]
    return varname
end

 function sliceprintln(::Type{<:SliceCtx}, ::Type{S}, ir::CodeInfo) where {S}
    slotnames = ir.slotnames
    vn = slotnames[end]
    varnames = Symbol[]
     funccalls = Any[]
     s = ""
     i = 1
    for expr in ir.code
        isa(expr, Expr) || continue
        if expr.head == :(=)
            @show expr
            try
                sym = Symbol(expr.args[1])
                vn = varname(ir, sym)
                push!(varnames, vn)
                if expr.args[2].head == :(call)
                    fname = expr.args[2].args[1]
                    push!(funccalls, (varname, fname))
                else
                    @show expr.args[2]
                    @show expr.args[2]
                end
            catch ex
                @show ex
                # @warn "could not find slotname for $(expr.args[1])"
                # @show slotnames
                # @show expr.args[1], s, i
            end
        end
    end
     if length(varnames) > 0
         @show varnames
     end
     # TODO: do something with functions without explicit assignment.
     # if length(varnames) == 0
     #     @show ir
     # end

     if length(funccalls) > 0
         @show funccalls
     end

    callbackslotname = gensym("callback")
    push!(ir.slotnames, callbackslotname)
    push!(ir.slotflags, 0x00)
    callbackslot = SlotNumber(length(ir.slotnames))
    getmetadata = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), Expr(:contextslot), QuoteNode(:metadata))
     # insert the initial `callbackslot` assignment into the IR.
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                 (stmt, i) -> i == 1 ? 2 : nothing,
                                 (stmt, i) -> [Expr(:(=), callbackslot, getmetadata), stmt])
     # replace all calls of the form `f(args...)` with `callback(f, args...)`, taking care to
    # properly destructure the returned `(result, callback)` into the appropriate statements
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                 (stmt, i) -> begin
                                    i > 1 || return nothing # don't slice the callback assignment
                                    stmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                    return Base.Meta.isexpr(stmt, :call) ? 3 : nothing
                                 end,
                                 (stmt, i) -> begin
                                     items = Any[]
                                     callstmt = Base.Meta.isexpr(stmt, :(=)) ? stmt.args[2] : stmt
                                     push!(items, Expr(:call, callbackslot, callstmt.args...))
                                     push!(items, Expr(:(=), callbackslot, Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), SSAValue(i), 2)))
                                     result = Expr(:call, Expr(:nooverdub, GlobalRef(Core, :getfield)), SSAValue(i), 1)
                                     if Base.Meta.isexpr(stmt, :(=))
                                         result = Expr(:(=), stmt.args[1], result)
                                     end
                                     push!(items, result)
                                     return items
                                 end)
     # replace return statements of the form `return x` with `return (x, callback)`
    Cassette.insert_statements!(ir.code, ir.codelocs,
                                  (stmt, i) -> Base.Meta.isexpr(stmt, :return) ? 2 : nothing,
                                  (stmt, i) -> begin
                                      return [
                                          Expr(:call, Expr(:nooverdub, GlobalRef(Core, :tuple)), stmt.args[1], callbackslot)
                                          Expr(:return, SSAValue(i))
                                      ]
                                  end)
    return ir
end
 const sliceprintlnpass = Cassette.@pass sliceprintln
 a = rand(3)
b = rand(3)
 function add(a, b)
    push_to_global_test_cache!(a)
    push_to_global_test_cache!(b)
    c = a + b
    push_to_global_test_cache!(c)
    return c
end
ctx = SliceCtx(pass=sliceprintlnpass, metadata = () -> nothing)
result, callback = Cassette.overdub(ctx, add, a, b)
@test result == a + b
@test isempty(global_test_cache)
@show callback()
@test global_test_cache == [a, b, result]
 println("done (took ", time() - before_time, " seconds)")
end


text:
function add_vec!(w::DenseVector, x::DenseVector, f::Function)
	if length(w) != length(x)
		throw(DimensionMismatch())
	end
	for i in 1:length(x)
		w[i] += f(x[i])
	end
end

function add_vec!(w::DenseVector, X::SparseMatrixCSC, f::Function)
	n, l = size(X)
	# Assumes that X is a column vector
	if length(w) != n || l != 1
		throw(DimensionMismatch())
	end
	I, J, V = findnz(X)
	for k in 1:length(I)
		w[I[k]] += f(V[k])
	end
end

function inner_prod{T<:Real}(x::DenseVector{T}, y::DenseVector{T})
	dot(x, y)
end

function inner_prod{T<:Real}(X::SparseMatrixCSC{T}, y::DenseVector{T})
	inner_product(y, X)
end

function inner_prod{T<:Real}(x::DenseVector{T}, Y::SparseMatrixCSC{T})
	n, l = size(Y)
	# Assumes that X is a column vector
	if length(x) != n || l != 1
		throw(DimensionMismatch())
	end
	if n == 0
		return zero(eltype(x))*zero(eltype(Y))
	end
	I, J, V = findnz(Y)
	if size(I) == 0
		return zero(eltype(x))*zero(eltype(Y))
	end
	s = x[I[1]] * V[1]
	@inbounds for k in 2:length(I)
		s += x[I[k]] * V[k]
	end
	s
end


text:
module Services

include("eutils.jl")

end


text:
load("core.jl")

arm1 = BernoulliArm(0.7)
draw(arm1)

arm2 = BernoulliArm(0.2)
draw(arm2)

arms = [arm1, arm2]

algo = EpsilonGreedy(0.7, 2)

initialize(algo, arms)

horizon = 100

for t = 1:horizon
  chosen_arm = select_arm(algo)
  reward = draw(arms[chosen_arm])
  update(algo, chosen_arm, reward)
end

initialize(algo, arms)

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
"""
    Phases
"""
module Phases

abstract type Phase end
abstract type AbstractTrainingPhase <: Phase end

struct TrainingPhase <: AbstractTrainingPhase end

struct ValidationPhase <: Phase end
struct TestPhase <: Phase end
struct InitializationPhase <: Phase end
struct CleanupPhase <: Phase end

export
    Phase,
    AbstractTrainingPhase,
    TrainingPhase,
    ValidationPhase,
    TestPhase,
    InitializationPhase,
    CleanupPhase
end # module

using .Phases


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libclangex_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libclangex")
JLLWrappers.@generate_main_file("libclangex", UUID("82994860-12e5-56b1-9971-c34819be692e"))
end  # module libclangex_jll


text:
cd("..")
using Pkg
Pkg.activate(".")

using MLKernelsMissing

kpca = fit(kPCA, GaussianKernel(1.0), convert(Array{Union{Float64, Missing}}, rand(10, 100)))

transform(kpca, convert(Array{Union{Float64, Missing}}, rand(10, 1)))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CUDNN_CUDA110_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CUDNN_CUDA110")
JLLWrappers.@generate_main_file("CUDNN_CUDA110", UUID("5fa5d4a9-0408-52e0-9638-7667eddd2fce"))
end  # module CUDNN_CUDA110_jll


text:
# ============================================================
# Automatic Differentiation Methods
# ============================================================

"""
$(SIGNATURES)

Calculate derivatives of Maxwellian
"""
function ∂maxwellian(u::Real, ρ, U, λ)
    Mu = u -> maxwellian(u, ρ, U, λ)
    return ForwardDiff.derivative(Mu, u)
end

"""
$(SIGNATURES)
"""
∂maxwellian(u::Real, prim::AV) = ∂maxwellian(u, prim[1], prim[2], prim[end])

function ∂maxwellian(u::AV, ρ, U, λ)
    ∂maxwellian.(u::AV, ρ, U, λ)
end

∂maxwellian(u::AV, prim::AV) = ∂maxwellian(u, prim[1], prim[2], prim[end])


text:
module GWAlpha

using DelimitedFiles
using Statistics
using RCall

include("sync_processing_module.jl")
include("relatedness_module.jl")
include("LMM_module.jl")
include("GWAlpha_ML_module.jl")
include("significance_testing_module.jl")

"""
# __________________________________________________________________
# Genomic prediction and genome-wide association using Pool-seq data

`PoolGPAS(;filename_sync::String, filename_phen::String, maf::Float64=0.001, depth::Int64=1, model::String=["GWAlpha", "MIXED", "GLMNET"][1], filename_random_covariate=nothing, random_covariate::String=["FST", "RELATEDNESS"][1], varcomp_est::String=["ML", "REML"][1], glmnet_alpha::Float64=collect(range(0.0,1.0,step=0.01,))[1], fpr::Float64=0.01, plot::Bool=false)`

Build genomic prediction models and perform genome-wide association on quantitative traits (GPAS) by inferring additive allelic effects using pool sequencing (Pool-seq) data.

# Inputs
1. *filename_sync* [String]: filename of the genotype data in [synchronized pileup format (.sync)](https://sourceforge.net/p/popoolation2/wiki/Manual/)
2. *filename_phen* [String]: filename of the phenotype data in one of two formats:
- **.py** extension for iterative maximum likelihood estimation i.e. `MODEL="GWAlpha"`, e.g.:
```
	Pheno_name='Phenotype Name';
	sig=0.06724693662723039;		# standard deviation
	MIN=0.0;						# minimum phenotype value
	MAX=0.424591738712776;			# maximum phenotype value
	perc=[0.2,0.4,0.6,0.8];			# cummulative pool sizes percentiles excluding the last pool
	q=[0.16,0.20,0.23,0.27,0.42];	# phenotype values corresponding to each percentile
```
- **.csv** extension for comma-separated headerless pool sizes and corresponding mean phenotypic values, e.g.:
```
	200.0,0.11988952929875112
	200.0,0.18030259365994225
	200.0,0.21548030739673382
	200.0,0.24966378482228616
	200.0,0.31328530259365983
```
3. *maf* [Float64]: minimum allele frequency threshold (default=0.001)
4. *depth* [Int64]: minimum sequencing depth threshold (default=1)
5. *model* [String]: model to use (default="GWAlpha")
- **"GWAlpha"** - iterative maximum likelihood estimation
- **"MIXED"** - linear mixed model
- **"GLMNET"** - [elastic-net penalization](https://web.stanford.edu/~hastie/Papers/glmnet.pdf)
6. *filename_random_covariate* [String]: filename of a precomputed headerless square symmetric matrix of pool relatedness (default=nothing)
7. *random_covariate* [String]: type of relatedness matrix to compute if `filename_random_covariate==nothing` (default="FST")
- **"FST"** - pairwise estimation of fixation indices using Pool-seq data using [Weir and Cockerham, 1984 method](https://www.jstor.org/stable/2408641?seq=1) (additionally [Hivert et al, 2018 method](https://www.biorxiv.org/content/biorxiv/early/2018/03/20/282400.full.pdf) is also available: see `?GWAlpha.relatedness_module.Fst_pairwise`)
- **"RELATEDNESS"** - simple standardized relatedness matrix `XX'/p`, where `X` is the allele frequency matrix (Pool-seq data) and `p` is the number of columns of `X`
8. *glmnet_alpha* [Float64]: elastic-net penalty (default=0.00 or ridge regression penalty)
9. *varcomp_est* [String]: variance component estimation method
- **"ML"** - maximum likelihood estimation
- **"REML"** - restricted maximum likelihood estimation
10. *glmnet_alpha* [Float64]: elastic-net penalty (from 0.00 for ridge regression up to 1.00 for LASSO)
11. *fpr* [Float64]: false positive rate for computing the Bonferroni threshold in significance testing (default=0.01)
12. *plot* [Bool]: generate Manhattan and quantile-quantile (QQ) plots, and save in a portable network graphics (.png) format (default=false)

# Outputs
1. Additive allelic effects array (header: CHROM, POS, ALLELE, FREQ, ALPHA, PVALUES, LOD) written into a comma-separated (.csv) file
- "GWAlpha": `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "-GWAlpha-OUTPUT.csv")`
- "MIXED": `string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, varcomp_est, "_", random_covariate, "-OUTPUT.csv")`
- "GLMNET": `string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, "_ALPHA", glmnet_alpha, "-OUTPUT.csv")`
2. Random covariate effects vector (headerless: RANDOM_EFFECTS) written into a comma-separated (.csv) file
- "GWAlpha": nothing
- "MIXED": `string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, "_", random_covariate, "-RANEF-OUTPUT.csv")`
- "GLMNET": nothing
3. Manhattan and QQ plots in .png format
- `string(join(split(filename_output_csv, ".")[1:(end-1)], '.'), ".png")`
4. Parsing, filtering, and relatedness matrix outputs:
- Parsed sync data into a .csv file of allele frequency data (headerless: CHROM, POS, ALLELE, FREQ_POOL1, ..., FREQ_POOLn): `string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_ALLELEFREQ.csv")`
- Filtered sync data into a sync file: filename_sync_filtered = `string(join(split(filename_sync, ".")[1:(end-1)], "."), "_MAF", maf, "_DEPTH", depth, ".sync")`
- Pairwise pool relatedness (square, symmetric and headerless): `string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "_COVARIATE_", random_covariate, ".csv")`

# Examples
```
### Input files:
filename_sync = "test/test.sync"
filename_phen_py = "test/test.py"
filename_phen_csv = "test/test.csv"
### Single-threaded execution:
using GWAlpha
@time OUT_GWAS = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_py, maf=0.001, depth=10, model="GWAlpha", fpr=0.01, plot=true)
@time OUT_MIXEDML_FST = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_csv, maf=0.001, depth=10, model="MIXED", random_covariate="FST", fpr=0.01, plot=true)
@time OUT_MIXEDREML_RELATEDNESS = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_csv, maf=0.001, depth=10, model="MIXED", random_covariate="RELATEDNESS", varcomp_est="REML", fpr=0.01, plot=true)
@time OUT_GLMNET_RIDGE = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_csv, maf=0.001, depth=10, model="GLMNET", glmnet_alpha=0.00, fpr=0.01, plot=true)
@time OUT_GLMNET_LASSO = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_csv, maf=0.001, depth=10, model="GLMNET", glmnet_alpha=1.00, fpr=0.01, plot=true)
### Multi-threaded execution (parallel execution only applicable to model=="GWAlpha"):
using Distributed
Distributed.addprocs(length(Sys.cpu_info())-1)
@everywhere using GWAlpha
@time OUT_GWAS = GWAlpha.PoolGPAS(filename_sync=filename_sync, filename_phen=filename_phen_py, maf=0.001, depth=10, model="GWAlpha", fpr=0.01, plot=true)
```
# Details
The GWAlpha model is defined as α = W(μₐₗₗₑₗₑ-μₐₗₜₑᵣₙₐₜᵢᵥₑ)/σᵧ, where:
- μ is the mean of the beta distribution Beta(θ) where θ={θ₁,θ₂}
- θ is estimated via maximum likelihood L(θ|Q) ∝ πᵢ₌₁₋ₖf(qᵢ|θ)
- Q = {q₁,...,qₖ} is the cumulative sum of allele frequencies across increasing-phenotypic-value-sorted pools where k is the number of pools
- E(allele|θ) = Beta_cdf(yᵢ',θ) - Beta_cdf(yᵢ₋₁',θ), where yᵢ' ∈ Y'
- Y' is the inverse quantile-normalized into phenotype data such that Y' ∈ [0,1]
- W = 2√{E(allele)*(1-E(allele))} is the penalization for low allele frequency

The linear mixed model is defined as y = Xb + Zu + e, where:
- X [n,p] is the centered matrix of allele frequencies
- Z [n,n] is the square symmetric matrix of relatedness
- y [n,1] is the centered vector of phenotypic values
- no intercept is explicitly fitted but implicitly set at the mean phenotypic value as a consequence of centering y
- u ~ N(0, σ²uI)
- e ~ N(0, σ²eI)
- y ~ N(0, V); V = (Z (σ²uI) Z') + (σ²eI)
- variance component (σ²e, σ²u) are estimated via maximum likelihood (ML) or restricted maximum likelihood (REML)
- fixed effects (b) are estimated via least squares: (X' * inverse(V)) * inverse(X * X' * inverse(V)) * y
- random effects (u) are estimated by solving: (σ²uI) * Z' * inverse(V) * (y - (X*b))

Empirical p-values were calculated by modelling the additive allelic effects (α) using a normal distribution with mean and variance estimated using maximum likelihood.

"""
function PoolGPAS(;filename_sync::String, filename_phen::String, maf::Float64=0.001, depth::Int64=1, model::String=["GWAlpha", "MIXED", "GLMNET"][1], filename_random_covariate=nothing, random_covariate::String=["FST", "RELATEDNESS"][1], varcomp_est::String=["ML", "REML"][1], glmnet_alpha::Float64=collect(range(0.0,1.0,step=0.01,))[1], fpr::Float64=0.01, plot::Bool=false)
	# ### TEST
	# filename_sync = "../test/test.sync"
	# filename_phen = ["../test/test.py", "../test/test.csv"][2]
	# maf = 0.001
	# depth = 10
	# model = ["GWAlpha", "MIXED", "GLMNET"][2]
	# filename_random_covariate = nothing
	# random_covariate = ["FST", "RELATEDNESS"][1]
	# varcomp_est = "REML"
	# glmnet_alpha = 0.50
	# fpr = 0.05
	if model == "GWAlpha"
		println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
		println("PERFORMING ITERATIVE GWALPHA")
		println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
		OUT = GWAlpha_ML_module.GWAlpha_ML(filename_sync=filename_sync, filename_phen_py=filename_phen, MAF=maf)
		filename_output_csv = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "-GWAlpha-OUTPUT.csv")
		writedlm(filename_output_csv, reshape(vcat(string.(keys(OUT))...), 1, length(OUT)), ',') ### write header
		io = open(filename_output_csv, "a"); writedlm(io, hcat(OUT...), ','); close(io) ### append the column-binded GWAlpha output
	else
		println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
		println("PERFORMING NON-ITERATIVE GWALPHA")
		println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
		### filter the sync file by maf and depth if the input file has not yet been filtered
		split_fname = split(filename_sync, "_")
		bool_maf = try
						maf == parse(Float64, split(split_fname[.!isnothing.(match.(r"MAF", split_fname))][1], "MAF")[2])
					catch
						false
					end ### filtered by minimum allele frequency?
		bool_depth = try
						depth == parse(Int64, split(split(split_fname[.!isnothing.(match.(r"DEPTH", split_fname))][1], "DEPTH")[2], ".")[1])
					catch
						false
					end ### filtered by minimum sequencing depth?
		if (bool_maf & bool_depth) == false
			println(string("Filtering the genotype data (sync format): \"", filename_sync, "\" by minimum allele frequency: ", maf, " and minimum depth: ", depth, "."))
			sync_processing_module.sync_filter(filename_sync=filename_sync, MAF=maf, DEPTH=depth);
			filename_sync_filtered = string(join(split(filename_sync, ".")[1:(end-1)], "."), "_MAF", maf, "_DEPTH", depth, ".sync")

		else
			filename_sync_filtered = filename_sync
		end
		### parse the filtered sync file if the input file has not yet been parsed
		filename_sync_parsed = string(join(split(filename_sync, ".")[1:(end-1)], '.'), "_ALLELEFREQ.csv")
		if !isfile(filename_sync_parsed)
			println(string("Parsing the filtered genotype data (sync format): \"", filename_sync, "\""))
			sync_processing_module.sync_parse(filename_sync_filtered)
			filename_sync_parsed = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "_ALLELEFREQ.csv")
		else
			filename_sync_parsed = filename_sync_parsed
		end
		println(string("Load the filtered and parsed genotype data (csv format): \"", filename_sync_parsed, "\"."))
		GENO = DelimitedFiles.readdlm(filename_sync_parsed, ',')
		idx_non_zero = sum(GENO[:,4:end],dims=2)[:,1] .!= 0.0
		X = convert(Array{Float64,2}, GENO[idx_non_zero, 4:end]')
		println(string("Loading phenotype data (csv format): ", filename_phen))
		PHENO = DelimitedFiles.readdlm(filename_phen, ',')
		n_pools = size(PHENO,1)
		pool_sizes = convert(Array{Int64,1},PHENO[:,1])
		y = convert(Array{Float64,1},PHENO[:,2])
		println(string("Extracting loci information from \"", filename_sync_parsed, "\"."))
		CHROM = convert(Array{Any,1}, GENO[idx_non_zero, 1])
		POS = convert(Array{Int64,1}, GENO[idx_non_zero, 2])
		ALLELE = convert(Array{String,1}, GENO[idx_non_zero, 3])
		FREQ = convert(Array{Float64,1}, sum(X .* (pool_sizes/sum(pool_sizes)'), dims=1)[1,:])
		if model == "MIXED"
			println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
			println("Performing Mixed Modelling")
			println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
			### load random covariate of relatedness (Fst matrix or standardized relatedness matrix) computed from the filtered sync file
			if isnothing(filename_random_covariate)
				filename_random_covariate = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "_COVARIATE_", random_covariate, ".csv")
				println(string("Building relatedness matrix to be used as the random covariate: ", filename_random_covariate))
				if random_covariate == "FST"
					relatedness_module.Fst_pairwise(filename_sync=filename_sync_filtered, window_size=100000, pool_sizes=pool_sizes) ### default="WeirCock" [Weir and Cockerham, 1984 method](https://www.jstor.org/stable/2408641?seq=1)
				elseif random_covariate == "RELATEDNESS"
					relatedness_module.standardized_relatedness(filename_sync_filtered)
				end
			else
				random_covariate = "USER_DEFINED_COVARIATE"
			end
			println(string("Loading relatedness matrix to be used as the random covariate: ", filename_random_covariate))
			Z = DelimitedFiles.readdlm(filename_random_covariate, ',')
			println(string("Identifying the mixed model to use: ", model))
			METHOD_VAR_EST = string(split(model, "_")[1])
			println(string("Fitting a mixed model using ", METHOD_VAR_EST, " variance estimation with ", random_covariate, " as the random covariate."))
			b0, b_hat, u_hat, Ve, Vu = LMM_module.LMM(X=X, y=y, Z=Z, METHOD_VAR_EST=varcomp_est)
			filename_output_csv = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, varcomp_est, "_", random_covariate, "-OUTPUT.csv")
			filename_ranef_csv = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, varcomp_est, "_", random_covariate, "-RANEF-OUTPUT.csv")
			filename_varcomp_csv = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, varcomp_est, "_", random_covariate, "-VARCOMP_Ve_Vu-OUTPUT.csv")
			writedlm(filename_ranef_csv, u_hat, ',') ### write headerless random effects
			writedlm(filename_varcomp_csv, [Ve, Vu], ',') ### write headerless variances
		else
			println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
			println("Performing Elastic-Net Regression")
			println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
			@rput X
			@rput y
			@rput glmnet_alpha
			R"library(glmnet)";
			R"lambda = glmnet::cv.glmnet(x=X, y=y, alpha=glmnet_alpha, intercept=TRUE)$lambda.min";
			R"glmnet_out = glmnet::glmnet(x=X, y=y, alpha=glmnet_alpha, intercept=TRUE, lambda=lambda)";
			R"b0 = as.numeric(glmnet_out$a0)";
			R"b_hat = as.numeric(glmnet_out$beta)";
			@rget b0
			@rget b_hat
			filename_output_csv = string(join(split(filename_sync_filtered, ".")[1:(end-1)], '.'), "-", model, "_ALPHA", glmnet_alpha, "-OUTPUT.csv")
		end
		### write-out the additive allelic effects
		println("Significance testing")
		PVAL, LOD = significance_testing_module.estimate_pval_lod(b_hat)
		println("Summarizing")
		OUT = (CHROM=vcat(["Intercept"], CHROM), POS=vcat([0], POS), ALLELE=vcat(["Intercept"], ALLELE), FREQ=vcat([0.0], FREQ), ALPHA=vcat([b0], b_hat), PVALUES=vcat([0.0], PVAL), LOD=vcat([0.0], LOD))
		writedlm(filename_output_csv, reshape(vcat(string.(keys(OUT))...), 1, length(OUT)), ',') ### write header
		io = open(filename_output_csv, "a"); writedlm(io, hcat(OUT...), ','); close(io) ### append the column-binded GWAlpha output
	end
	### Manhattan plotting
	if plot
		filename_output_png = string(join(split(filename_output_csv, ".")[1:(end-1)], '.'), ".png")
		significance_testing_module.plot_manhattan(chrom=OUT.CHROM, pos=OUT.POS, pval=OUT.PVALUES, lod=OUT.LOD, fpr=fpr, png_fname=filename_output_png)
	end
	### Output reporting
	println("===============================================================")
	println("Everything went well. Please check the output files:")
	println(filename_output_csv)
	if model == "MIXED"
		println(filename_ranef_csv)
		println(filename_varcomp_csv)
	end
	if plot
		println(filename_output_png)
	end
	println("===============================================================")
	return(0)
end

end

max_stars_repo_name: jeffersonfparil/GWAlpha


text:
function adaptive_sampling(A,b,obj,noisyp,Atest,btest,testx,instance,S)
    Experiment = 1

    I = 100
	TS = 100
    re_solve_counter = 0
    cone = fill(zeros(c,n), 1)
    Pred_Error_Evolution = zeros(Float64, I, 10)
    ind = 0

    A_train = zeros(Float64, c, n, I)
    b_train = zeros(Float64, I, c)
    x_train = zeros(Float64, 30, n, I)

    # Generate random folders to store BARON temp files
    id=randstring(5)
    @eval @everywhere id=$id
    mkdir(string(root,"Results//ex",id))
    for t=1:S
        mkdir(string(root,"Results//ex",id,"//opts",t))
    end

    while(Experiment <= I)
        T = Experiment
        A_train[:,:,Experiment] = deepcopy(A);
        b_train[Experiment, :] = deepcopy(b);
        x = SolveModel(A, obj, b, n)
        J = 30;
        noisy_x = zeros(Float64, J, n)
        for j = 1:J
            for k = 1:n
                noisy_x[j, k] = x[k] + 3*randn(rng, Float64)
            end
        end
        x_train[:,:,Experiment] = deepcopy(noisy_x);

        hat_x = zeros(Float64, I, n)
        hat_z = zeros(Float64, I, c)

        hat_x[T, :], r, z = SolveP1(A_train[:,:,T], b_train[T,:], noisy_x, 1, J)
        istight = []
        istight, z = tightind(A_train[:,:,T], b_train[T,:], hat_x[T, :])
        hat_z[T, :] = z

        if Experiment == 1
            for j = 1:length(istight)
                cone[Experiment][j,:] = A[istight[j], :]/norm(A[istight[j], :])
            end
        else
            B = zeros(Float64, length(istight), n)
            for j = 1:length(istight)
                B[j,:] = A[istight[j], :]/norm(A[istight[j], :])
            end
            push!(cone, B)
        end
        status = SolveFP(cone, T)

        if status == :Infeasible
            re_solve_counter += 1
            z = []
            x, r, z, status = SolveP1m(A_train[:,:,1:T], b_train[1:T,:], x_train[:,:,1:T], T, J, instance, hat_x, hat_z)

            if status != :Optimal
                break
            end
            hat_x[1:T, :] = x
            hat_z[1:T, :] = z

            cone = fill(zeros(c,n), 1)

            for k = 1:T
                istight = []
                istight, z = tightind(A_train[:,:,k], b_train[k,:], hat_x[k, :])
                if k == 1
                    for j = 1:length(istight)
                        cone[k][j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                else
                    B = zeros(Float64, length(istight), n)
                    for j = 1:length(istight)
                        B[j,:] = A_train[istight[j], :, k]/norm(A_train[istight[j], :, k])
                    end
                    push!(cone, B)
                end
            end
        end

        for estimate = 1:10
            hatc, loss = SolveP2(noisyp[estimate, :], cone, T)
            x_pred  = zeros(Float16, TS, n)

            for test = 1:TS
                x_pred[test, 1:n] = SolveModel(Atest[:, :, test], hatc, btest[test, :], n)
            end
            Pred_Error_Evolution[T, estimate] = infnorm(testx, x_pred)
        end

        # Initialize arrays to store samples for heuristic-based adaptive sampling
        global bprime = zeros(Float64, S, c)
        global Aprime = zeros(Float64, c, n, S)
        xprime = zeros(Float64, S, dim)
        global zprime = zeros(Int16, S, c)
        phi = zeros(Float64, S) # To store the objective values of (ASP)
        Boun = zeros(Float64, S) # To store the best bound for (ASP)
        G = zeros(Float64, S)

        for sample = 1:S
            mutrain = zeros(Float64, II, JJ)
            Dtrain = zeros(Float64, JJ, TT)
            for i = 1:II
                for j = 1:JJ
                    mutrain[i,j] = (rand(rng, 50:100)/100)*Mumax[i,j]
                    for t = 1:TT
                        Dtrain[j, t] = (rand(rng, 90:110)/100)*Nominal_Demand[j]
                    end
                end
            end
            Aprime[:,:,sample], bprime[sample, :] = ModelParameters(II,JJ,TT,mutrain,Dtrain)
            hat_c = SampleC(cone) # Sample a random vector from the admissible set
            xprime[sample, :] = SolveModel(Aprime[:,:,sample], hat_c, bprime[sample,:], n)
            temp, zprime[sample, :] = tightind(Aprime[:,:,sample], bprime[sample,:], xprime[sample,:])
        end

        # Solve the adaptive sampling problem for the samples in Aprime, bprime
        Results = zeros(Float64, S, 3)
        Results = pmap(Tr->maxmin(Tr, cone, Aprime[:,:,Tr], zprime[Tr,:]), 1:S, retry_delays=ones(10))
        for i = 1:S
            phi[i] = Results[i][1]
            Boun[i] = Results[i][2]
        end
        ind = findmax(phi)

        # Update A, b to the sample with the highest value of \eta
        b = bprime[ind[2],:]
        A = Aprime[:,:, ind[2]]
        Experiment += 1
    end
    cd(string(root,"Results"))
    df = DataFrame(Pred_Error_Evolution)
    s = @sprintf "PredictionError_dim_%03.d_apaptivesamp_instance_%03.d.csv" dim instance
    CSV.write(s, df)
end

max_stars_repo_name: ddolab/InvLP

max_issues_repo_name: ddolab/InvLP

max_forks_repo_name: ddolab/InvLP


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: iot_events_data
using AWS.Compat
using AWS.UUIDs

"""
    batch_put_message(messages)
    batch_put_message(messages, params::Dict{String,<:Any})

Sends a set of messages to the AWS IoT Events system. Each message payload is transformed
into the input you specify (\"inputName\") and ingested into any detectors that monitor
that input. If multiple messages are sent, the order in which the messages are processed
isn't guaranteed. To guarantee ordering, you must send messages one at a time and wait for
a successful response.

# Arguments
- `messages`: The list of messages to send. Each message has the following format: '{
  \"messageId\": \"string\", \"inputName\": \"string\", \"payload\": \"string\"}'

"""
batch_put_message(messages; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}("messages"=>messages); aws_config=aws_config)
batch_put_message(messages, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/inputs/messages", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("messages"=>messages), params)); aws_config=aws_config)

"""
    batch_update_detector(detectors)
    batch_update_detector(detectors, params::Dict{String,<:Any})

Updates the state, variable values, and timer settings of one or more detectors (instances)
of a specified detector model.

# Arguments
- `detectors`: The list of detectors (instances) to update, along with the values to update.

"""
batch_update_detector(detectors; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}("detectors"=>detectors); aws_config=aws_config)
batch_update_detector(detectors, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("POST", "/detectors", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("detectors"=>detectors), params)); aws_config=aws_config)

"""
    describe_detector(detector_model_name)
    describe_detector(detector_model_name, params::Dict{String,<:Any})

Returns information about the specified detector (instance).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) you
  want information about.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"keyValue"`: A filter used to limit results to detectors (instances) created because of
  the given key ID.
"""
describe_detector(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/"; aws_config=aws_config)
describe_detector(detectorModelName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)/keyValues/", params; aws_config=aws_config)

"""
    list_detectors(detector_model_name)
    list_detectors(detector_model_name, params::Dict{String,<:Any})

Lists detectors (the instances of a detector model).

# Arguments
- `detector_model_name`: The name of the detector model whose detectors (instances) are
  listed.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results to return at one time.
- `"nextToken"`: The token for the next set of results.
- `"stateName"`: A filter that limits results to those detectors (instances) in the given
  state.
"""
list_detectors(detectorModelName; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)"; aws_config=aws_config)
list_detectors(detectorModelName, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = iot_events_data("GET", "/detectors/$(detectorModelName)", params; aws_config=aws_config)


text:
# la variable de la boucle for est locale

function f(x)
    print(x, ".")
    x
end

n = 42
println(n)

for n = f(1) : f(10)
    print(n)
end
println()

println(n)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test

@testset "PEPSAD.jl" begin
    @testset "$filename" for filename in
                             ["peps_test.jl", "models_test.jl", "optimizations_test.jl"]
        println("Running $filename")
        include(filename)
    end
end


text:
module Results

using Printf
include("valuetostate.jl")

function createcsv(nodeSampleResults, columns, resultfilename)
    outfile = open(resultfilename, "w")

    sortedColumns = sort(collect(columns))
    sortedColumnsFull = copy(sortedColumns)
    pushfirst!(sortedColumnsFull, "node")

    write(outfile, join(sortedColumnsFull, "\t"))
    write(outfile, "\n")

    for node in keys(nodeSampleResults)
        if occursin("PSEUDONODE", node)
            continue
        end
        write(outfile, "$node\t")

        index = 0
        lengthColumns = length(columns)
        for column in sortedColumns
            index += 1
            value = @sprintf "%.2f" nodeSampleResults[node][column]
            write(outfile, "$value")
            if lengthColumns == index
                write(outfile, "\n")
            else
                write(outfile, "\t")
            end
        end
    end

    close(outfile)
end

function outputAllResults(nodesampleresults, columns, resultfilename)
    outfile = open(resultfilename, "w")

    write(outfile, join(columns, "\t"))
    write(outfile, "\n")

    for node in keys(nodesampleresults)
        write(outfile, "$node\t")

        index = 0
        lengthColumns = length(columns)
        for column in columns
            index += 1
            if column == "gene"
                continue
            end
            (x, x_bar) = nodesampleresults[node][column]
            write(outfile, "$x:$x_bar")
            if lengthColumns == index
                write(outfile, "\n")
            else
                write(outfile, "\t")
            end
        end
    end

    close(outfile)
end

function getResults(resultfilename)
    result_data = readlines(resultfilename)

    header = popfirst!(result_data)
    headerparts = split(chomp(header), "\t")

    samplenodestate = Dict()

    i = 0
    for column in headerparts
        i = i + 1
        if i != 1
            samplenodestate[column] = Dict()
        end
    end

    counts = Dict("0" => 0, "1" => 0, "2" => 0)

    node = ""
    for line in result_data
        lineparts = split(chomp(line), "\t")
        i = 0
        for column in headerparts
            i = i + 1
            if i == 1
                node = lineparts[1]
            else
                samplenodestate[column][node] = lineparts[i]
            end
        end
    end

    return samplenodestate
end

function getExpected(expectedfile)
    expected_data = readlines(expectedfile)

    header = popfirst!(expected_data)
    headerparts = split(chomp(header), "\t")

    samplenodestate = Dict()

    i = 0
    for column in headerparts
        i = i + 1
        if i > 1
            samplenodestate[column] = Dict()
        end
    end

    counts = Dict("0" => 0, "1" => 0, "2" => 0)

    node = ""
    for line in expected_data
        lineparts = split(chomp(line), "\t")
        i = 0
        for column in headerparts
            i = i + 1
            if i == 1
                node = lineparts[1]
            else
                state = lineparts[i]
                if state != "-999"
                    counts[state] = counts[state] + 1
                    samplenodestate[column][node] = state
                end
            end
        end
    end

    return Dict("counts" => counts, "samplenodestate" => samplenodestate)
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
function parse_contents(contents, filename)
  content_type, content_string = split(contents, ',')
  decoded = base64decode(content_string)
  df = DataFrame()
  try
    if occursin("csv", filename)
      str = String(decoded)
      df =  CSV.read(IOBuffer(str), DataFrame)
    end
  catch e
    print(e)
    return html_div([
        "There was an error processing this file."
    ])
  end
   #shuffle data
   n_row = nrow(df)
   shuffle_data = df[shuffle(1:n_row)[:], :]
   shuffle_data[!,"id"] = 1:n_row
  return shuffle_data
end

function reset_vars()
  global m, ep, st, x_graph, y_graph
  m, ep, st, x_graph, y_graph = nothing, 0, false, [], []
  return 0
end

max_stars_repo_name: DEICIDIA/fluxdash

max_issues_repo_name: DEICIDIA/fluxdash

max_forks_repo_name: DEICIDIA/fluxdash


text:
module ColoredGraphs

using LightGraphs
using MacroTools: @forward

import Base:
    eltype, show, ==, Pair, Tuple, copy, length, issubset, reverse, zero, in, iterate

import LightGraphs:
    _NI, AbstractGraph, AbstractEdge, AbstractEdgeIter,
    src, dst, edgetype, nv, ne, vertices, edges, is_directed,
    has_vertex, has_edge, inneighbors, outneighbors,

    indegree, outdegree, degree, has_self_loops, num_self_loops, insorted

export ColoredGraph, ColoredDiGraph, ColoredEdge

export HyperCube, translational_symm_table

"""
    AbstractColoredEdge{T} <: AbstractEdge{T}

Abstract base type for an edge with a color.
"""
abstract type AbstractColoredEdge{T} <: AbstractEdge{T} end

"""
    AbstractColoredGraph{T} <: AbstractGraph{T}

Abstract base type for colored graphs.
"""
abstract type AbstractColoredGraph{T} <: AbstractGraph{T} end

include("ColoredGraphs/ColoredEdge.jl")
include("ColoredGraphs/ColoredEdgeIter.jl")
include("ColoredGraphs/ColoredGraph.jl")

include("WrappedSimpleGraphs/AbstractWrappedSimpleGraph.jl")
include("WrappedSimpleGraphs/HyperCube.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 066c5b84-c19e-11ec-3f27-115d87575e0d
struct Dual{T}
	val::T
	der::T
end

# ╔═╡ fd3fcf64-b334-4ba7-862f-cdc02b221d6e
a = Dual(3.0,5.0)

# ╔═╡ 393d4637-21a3-4ce3-b720-bca1439c3233
typeof(a.val)

# ╔═╡ ca1251ce-b63f-4075-94c8-c0d184655bdc
Base.:+(f::Dual, g::Dual) = Dual(f.val + g.val, f.der + g.der)

# ╔═╡ 5d59360d-fd53-4996-8400-7d218e4eb29f
Base.:+(f::Dual, α::Number) = Dual(f.val + α, f.der)

# ╔═╡ 32d7e717-6a90-4975-a20a-80c7cabcb85c
Base.:+(α::Number, f::Dual) = f + α

# ╔═╡ 8cdaede1-9825-4404-8898-295c4b0f6e94
Base.:-(f::Dual, g::Dual) = Dual(f.val - g.val, f.der - g.der)

# ╔═╡ 4acc7af2-7eba-4d68-a106-245fb310b9b9
Base.:*(f::Dual, g::Dual) = Dual(f.val * g.val, f.der*g.val + g.der*f.val)

# ╔═╡ 0b9d1885-e113-42f3-b808-90064338f9e2
Base.:^(f::Dual, n::Integer) = Base.power_by_squaring(f,n)

# ╔═╡ 0336e41f-aabd-40a9-af47-dfeed25daa68
begin
	import Base:exp
	exp(f::Dual) = Dual( exp(f.val), exp(f.val)*f.der )
end

# ╔═╡ 9f4518e5-b444-4bd9-8ecf-23450a4d9125
b = Dual(2.0,7.0)

# ╔═╡ e63e723a-e023-4a35-bdaa-be6a3bf4b628
a + b

# ╔═╡ eb39443c-1b71-4137-b777-086700d44a51
 a - b

# ╔═╡ 6b25fab1-aee7-497d-8055-444a4c03aedb
a * b

# ╔═╡ 85588519-64e0-41ba-872a-1942c5a6925d
a^3

# ╔═╡ d1227abf-b2e8-43f8-8efb-aa353455d5ba
a*(a+b)

# ╔═╡ 0609e93e-f7b4-462c-82ba-d674871dc4d5
exp(a)

# ╔═╡ Cell order:
# ╠═066c5b84-c19e-11ec-3f27-115d87575e0d
# ╠═fd3fcf64-b334-4ba7-862f-cdc02b221d6e
# ╠═393d4637-21a3-4ce3-b720-bca1439c3233
# ╠═ca1251ce-b63f-4075-94c8-c0d184655bdc
# ╠═5d59360d-fd53-4996-8400-7d218e4eb29f
# ╠═32d7e717-6a90-4975-a20a-80c7cabcb85c
# ╠═8cdaede1-9825-4404-8898-295c4b0f6e94
# ╠═4acc7af2-7eba-4d68-a106-245fb310b9b9
# ╠═0b9d1885-e113-42f3-b808-90064338f9e2
# ╠═0336e41f-aabd-40a9-af47-dfeed25daa68
# ╠═9f4518e5-b444-4bd9-8ecf-23450a4d9125
# ╠═e63e723a-e023-4a35-bdaa-be6a3bf4b628
# ╠═eb39443c-1b71-4137-b777-086700d44a51
# ╠═6b25fab1-aee7-497d-8055-444a4c03aedb
# ╠═85588519-64e0-41ba-872a-1942c5a6925d
# ╠═d1227abf-b2e8-43f8-8efb-aa353455d5ba
# ╠═0609e93e-f7b4-462c-82ba-d674871dc4d5


text:
#
# For the manifolds that are nonmutating only, we have to introduce a few special cases
#
function get_gradient!(p::GradientProblem{AllocatingEvaluation}, ::AbstractFloat, x)
    X = p.gradient!!(p.M, x)
    return X
end
function get_hessian!(p::HessianProblem{AllocatingEvaluation}, ::AbstractFloat, x, X)
    Y = p.hessian!!(p.M, x, X)
    return Y
end
function linesearch_backtrack(
    M::NONMUTATINGMANIFOLDS,
    F::TF,
    x,
    gradF::T,
    s,
    decrease,
    contract,
    retr::AbstractRetractionMethod=ExponentialRetraction(),
    η::T=-gradF,
    f0=F(x),
) where {TF,T}
    x_new = retract(M, x, s * η, retr)
    fNew = F(x_new)
    while fNew < f0 + decrease * s * inner(M, x, η, gradF) # increase
        x_new = retract(M, x, s * η, retr)
        fNew = F(x_new)
        s = s / contract
    end
    s = s * contract # correct last
    while fNew > f0 + decrease * s * inner(M, x, η, gradF) # decrease
        s = contract * s
        x_new = retract(M, x, s * η, retr)
        fNew = F(x_new)
    end
    return s
end
# modify gradient descent step_solver
function step_solver!(
    p::GradientProblem{T,<:NONMUTATINGMANIFOLDS}, o::GradientDescentOptions, iter
) where {T}
    s, o.gradient = o.direction(p, o, iter)
    o.x = retract(p.M, o.x, -s .* o.gradient, o.retraction_method)
    return o
end


text:
using Zarr, AWS, Downloads, DataFrames, CSV, CFTime, Dates, Statistics, NetCDF

"""
    cmip(institution_id,source_id,variable_id)

Access CMIP6 climate model archive (https://bit.ly/2WiWmoh) via
`AWS.jl` and `Zarr.jl` and compute (1) time mean global map and
(2) time evolving global mean.

This example was partly inspired by @rabernat 's https://bit.ly/2VRMgvl notebook

```
using ClimateModels
(mm,gm,meta)=cmip()
nm=meta["long_name"]*" in "*meta["units"]

using Plots
heatmap(mm["lon"], mm["lat"], transpose(mm["m"]),
        title=nm*" (time mean)")
plot(gm["t"][1:12:end],gm["y"][1:12:end],xlabel="time",ylabel=nm,
     title=meta["institution_id"]*" (global mean, month by month)")
display.([plot!(gm["t"][i:12:end],gm["y"][i:12:end], leg = false) for i in 2:12])

```
"""
function cmip(institution_id="IPSL",source_id="IPSL-CM6A-LR",
    variable_id="tas")

    #choose model and variable
    S=[institution_id, source_id, variable_id]

    #get list of contents for cloud storage unit
    url="https://storage.googleapis.com/cmip6/cmip6-zarr-consolidated-stores.csv"
    ξ = CSV.read(Downloads.download(url),DataFrame)

    # get model grid cell areas
    ii=findall( (ξ[!,:source_id].==S[2]).&(ξ[!,:variable_id].=="areacella") )
    μ=ξ[ii,:]
    cmip6,p = Zarr.storefromstring(μ.zstore[end])
    ζ = zopen(cmip6,path=p)
    Å = ζ["areacella"][:, :]

    # get model solution ensemble list
    i=findall( (ξ[!,:activity_id].=="CMIP").&(ξ[!,:table_id].=="Amon").&
    (ξ[!,:variable_id].==S[3]).&(ξ[!,:experiment_id].=="historical").&
    (ξ[!,:institution_id].==S[1]) )
    μ=ξ[i,:]

    # access one model ensemble member
    cmip6,p = Zarr.storefromstring(μ.zstore[end])
    ζ = zopen(cmip6,path=p)

    meta=Dict("institution_id" => institution_id,"source_id" => source_id,
        "variable_id" => variable_id, "units" => ζ[S[3]].attrs["units"],
        "long_name" => ζ[S[3]].attrs["long_name"])

    # time mean global map
    m = convert(Array{Union{Missing, Float32},3},ζ[S[3]][:,:,:])
    m = dropdims(mean(m,dims=3),dims=3)

    mm=Dict("lon" => ζ["lon"], "lat" => ζ["lat"], "m" => m)

    # time evolving global mean
    t = ζ["time"]
    t = timedecode(t[:], t.attrs["units"], t.attrs["calendar"])

    y = ζ[S[3]][:,:,:]
    y=[sum(y[:, :, i].*Å) for i in 1:length(t)]./sum(Å)

    gm=Dict("t" => t, "y" => y)

    return mm,gm,meta
end

"""
	IPCC_hexagons()

Read hexagons used in IPPC AR6 report. DataFrame contains
acronym, name, region, and the x, y coordinate (integers)
of each losange. There arrangement mimics the distribution
of continents as done in the report and interactive atlas.

```
df=IPCC_hexagons()
```
"""
function IPCC_hexagons()
	fil2=joinpath(IPCC_SPM_path,"reference-regions","hexagon_grid_locations.csv")
	DataFrame(CSV.File(fil2))
end

"""
	IPCC_fig3_example(df)

```
df=IPCC_hexagons()
clv, ttl, colors=IPCC_fig3_example(df)
```
"""
function IPCC_fig3_example(df)
	fil2=joinpath(IPCC_SPM_path,"reference-regions","hexagon_grid_locations.csv")
	df=DataFrame(CSV.File(fil2))

	#T=Union{Makie.LinePattern,Symbol}
	colors=Array{Symbol}(undef,size(df,1))

    #set colors to co where df.acronym.==acr
    function set_col!(df,acr,colors,co)
        k=findall(df.acronym.==acr)[1]
        colors[k]=co
    end

	colors.=:lightcoral
	set_col!(df,"CNA",colors,:yellow)
	set_col!(df,"ENA",colors,:yellow)
	set_col!(df,"SSA",colors,:lightgray)
	set_col!(df,"CAF",colors,:lightgray)

	ttl="a) Synthesis of assessment of observed change in hot extremes \n and confidence in human contribution to the observed changes \n in the world’s regions"

	confidencelevels=fill("⚪?",size(df,1))
	cl=["⚫⚫⚫","⚫⚫","⚫","⚪"]
	clv=Array{Any}(undef,4)
	clv[1] = [ "NWN","NEN","NEU","WCE","EEU","WSB","ESB","RFE","MED","WCA","ECA","TIB","EAS",
				"SAS","SEA","NWS","SES","WSAF","ESAF","NAU","CAU","EAU","SAU"]
	clv[2] = [ "GIC","RAR","WNA","NCA","SCA","CAR","SAH","ARP","PAC","NWF","NSA",
				"WAF","NEAF","SAM","NES","SEAF","SWS"]
	clv[3] = [ "CNA","ENA","MDG","NZ"]
	clv[4] = [ "CAF","SSA"]

    return clv, ttl, colors

end

#(dat_1b,meta_1b)=ClimateModels.IPCC_fig1b_read()
function IPCC_fig1b_read()
    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_01","v20210809")
    fil_1b=joinpath(pth_ipcc,"panel_b","gmst_changes_model_and_obs.csv")
    dat_1b=DataFrame(CSV.File(fil_1b; header=false, skipto=37, footerskip=1))
    rename!(dat_1b, Dict(:Column1 => :year, :Column2 => :HNm, :Column3 => :HN5, :Column4 => :HN95))
    rename!(dat_1b, Dict(:Column5 => :Nm, :Column6 => :N5, :Column7 => :N95, :Column8 => :obs))

	header=readlines(fil_1b)[12:34]

	long_name=String[]
	units=String[]
	comment=String[]
	types=String[]

	for i in findall(occursin.("long_name",header))
			tmp2=split(header[i],",")
			push!(long_name,tmp2[3])
			push!(units,tmp2[4])
	end
	for i in findall(occursin.("comment",header))
			tmp2=split(header[i],",")
			push!(comment,tmp2[3])
	end
	for i in findall(occursin.("type",header))
			tmp2=split(header[i],",")
			push!(types,tmp2[3])
	end

    meta_1b=(long_name,units,comment,types)

    return dat_1b,meta_1b
end

#(dat, dat1, dat2)=IPCC_fig1a_read()
function IPCC_fig1a_read()
    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_01","v20210809")
	files=readdir(joinpath(pth_ipcc,"panel_a"))

	input=joinpath(pth_ipcc,"panel_a",files[1])
	dat=DataFrame(CSV.File(input))
	rename!(dat, Dict(Symbol("5%") => :t5, Symbol("95%") => :t95))

	input=joinpath(pth_ipcc,"panel_a",files[2])
	dat1=sort(DataFrame(CSV.File(input)),:year)

	input=joinpath(pth_ipcc,"panel_a",files[3])
	dat2=DataFrame(CSV.File(input))
	rename!(dat2, Dict(Symbol("5%") => :t5, Symbol("95%") => :t95))

    dat, dat1, dat2
end

#(dat2a,dat2b,dat2c)=IPCC_fig2_read()
function IPCC_fig2_read()
    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_02","v20210809")
    files=readdir(joinpath(pth_ipcc,"panel_a"))

	input=joinpath(pth_ipcc,"panel_a",files[1])
	dat=DataFrame(CSV.File(input))

    files_b=readdir(joinpath(pth_ipcc,"panel_b"))
	input_b=joinpath(pth_ipcc,"panel_b",files_b[1])
	dat_b=DataFrame(CSV.File(input_b))

    files_c=readdir(joinpath(pth_ipcc,"panel_c"))
	input_c=joinpath(pth_ipcc,"panel_c",files_c[1])
	dat_c=DataFrame(CSV.File(input_c))

	dat_c[6,1]="Volatile organic compounds and \n carbon monoxide (NMVOC + CO)"
	dat_c[4,1]="Halogenated gases (CFC + \n HCFC + HFC)"
	dat_c[11,1]="Land-use reflectance and \n irrigation (irrig+albedo)"

    return dat,dat_b,dat_c
end

#dat=IPCC_fig4_read()
function IPCC_fig4a_read()
    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_04","v20210809")

    C=DataFrame(CSV.File(joinpath(pth_ipcc,"panel_a","Carbon_dioxide_Gt_CO2_yr.csv")))
    M=DataFrame(CSV.File(joinpath(pth_ipcc,"panel_a","Methane_Mt_CO2_yr.csv")))
    N=DataFrame(CSV.File(joinpath(pth_ipcc,"panel_a","Nitrous_oxide_Mt_N2O_yr.csv")))
    S=DataFrame(CSV.File(joinpath(pth_ipcc,"panel_a","Sulfur_dioxide_Mt_SO2_yr.csv")))

    (C=C,M=M,N=N,S=S)
end

#dat=IPCC_fig4b_read()
function IPCC_fig4b_read()
    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_04","v20210809")
    fil=joinpath(pth_ipcc,"panel_b","ts_warming_ranges_1850-1900_base_panel_b.csv")
    DataFrame(CSV.File(fil))
end

#dat=IPCC_fig5_read()
function IPCC_fig5_read(fil="Panel_a2_Simulated_temperature_change_at_1C.nc")

    nam=""
    occursin("temperature",fil) ? nam="tas" : nothing
    occursin("precipitation",fil) ? nam="pr" : nothing
    occursin("SM_tot",fil) ? nam="mrso" : nothing

    pth_ipcc=joinpath(IPCC_SPM_path,"spm","spm_05","v20210809")
	lst=readdir(pth_ipcc)
	readme="Readme_for_figure_SPM5.txt"

    lon = Float64.(NetCDF.open(joinpath(pth_ipcc,fil), "lon")[:])
    lat = Float64.(NetCDF.open(joinpath(pth_ipcc,fil), "lat")[:])
    var = Float64.(NetCDF.open(joinpath(pth_ipcc,fil), nam)[:,:,1])

    tmp=(;lon,lat)
    lon=[tmp.lon[i] for i in 1:length(tmp.lon), j in 1:length(tmp.lat)]
    lat=[tmp.lat[j] for i in 1:length(tmp.lon), j in 1:length(tmp.lat)]

    if nam=="tas"
        ttl="Annual mean temperature change (°C) relative to 1850-1900"
        meta=(colorrange=(0,6),cmap=:Reds_9,ttl=ttl)
        var=0.5*floor.(2*var,digits=0)
    elseif nam=="pr"
        ttl="Annual mean precipitation change (%) relative to 1850-1900"
        meta=(colorrange=(-40,40),cmap=:BrBG_10,ttl=ttl)
        var=5.0*floor.(0.2*var,digits=0)
    elseif nam=="mrso"
        ttl="Annual mean total column soil moisture change (standard deviation)"
        meta=(colorrange=(-1.5,1.5),cmap=:BrBG_10,ttl=ttl)
        var=0.25*floor.(4*var,digits=0)
    else
        error("unknown case")
    end

    (lon=lon,lat=lat,var=var,meta=meta)
end


text:
"""
    function dwz_ranking(
        players::Vector{Player},
        ranking_scheme_to_update::String,
        last_tournament_id::String
    )

Update the players rankings using the Deutsche Wertungszahl Algorithm

Arguments
=========
- `players`:
    an array of players in the simulation.
- `ranking_scheme_to_update`:
    the name of the ranking scheme to update for each Player.
- `last_tournament_id`:
    the last tournament that was completed.
"""

function dwz_ranking(
    players::Vector{Player},
    ranking_scheme_to_update::String,
    last_tournament_id::String
)
    updatedRankings = Dict{Int64, Int64}();
    for player in players
        if haskey(player.doubles_results, last_tournament_id)
            player_ranking = player.rankings["Dwz"]
            teammate_ranking = search_dwz_ranking(player.doubles_results.partner_player_id, players)
            achieved_points = 0
            expected_points = 0
            n = length(player.doubles_results[last_tournament_id])
            for match in player.doubles_results[last_tournament_id].matches
                zA = (player_ranking + teammate_ranking) / 2
                zG = search_dwz_ranking(match.opponent1_player_id, players) + search_dwz_ranking(match.opponent2_player_id, players)
                wE = 1/(1+(10^((zG-zA)/400)))
                wA = match.result == "Won" ? 1 : 0
                achieved_points += wA
                expected_points += wE
            end
            if player.age < 20
                J = 5
            elseif player.age <= 25
                J = 10
            else
                J = 15
            end
            E0 = (player_ranking/1000)^4 + J
            a = player.age <= 20 && achieved_points > expected_points ? player_ranking/2000 : 1
            B = player_ranking < 1300 && achieved_points <= expected_points ? exp((1300-player_ranking)/150)-1 : 0
            E = a * E0 + B
            Zn = player_ranking + (800/(E+n))(achieved_points-expected_points)
            updatedRankings[player.id] = Zn
        else
            # If player didn't participate in tournament, their ranking doesn't change
            updatedRankings[player.id] = player.rankings["Dwz"];
        end
    end
    map(player -> player.rankings[ranking_scheme_to_update] = updatedRankings[player.id], players)
    return nothing
end

"""
    function search_dwz_ranking(
        player_id::Int64,
        players::Vector{Player}
    )::Float64

Search for a player's current Dwz Ranking

Arguments
=========
- `player_id`:
    the player whose ranking is under search.
- `players`:
    an array of players in the simulation.
"""

function search_dwz_ranking(
    player_id::Int64,
    players::Vector{Player}
)::Float64
    player_ranking = filter(p -> p.id == player_id, players)[1].rankings["Dwz"]
    return player_ranking
end

max_stars_repo_name: Danial-Hussain/Tennis-Ranking

max_issues_repo_name: Danial-Hussain/Tennis-Ranking

max_forks_repo_name: Danial-Hussain/Tennis-Ranking


text:
using Printf
using LaTeXStrings
using Plots, ColorSchemes
using JLD2, FileIO

import CairoMakie

include("plot-util.jl")

Base.zero(x::Type{Array{Float64,1}}) = Float64[]

function unpack(key)
    elt = split(key,'/')
    return (
       hgt=parse(Float64,elt[1]),
       snp=parse(Float64,elt[2]),
       nit=parse(Int64,elt[3]),
       elt=elt[4],
    )
end

function params(data)
    hgt = Float64[]
    snp = Float64[]

    for key in keys(data)
        param = unpack(key)
        param.hgt ∈ hgt || push!(hgt, param.hgt)
        param.snp ∈ snp || push!(snp, param.snp)
    end

    sort!(hgt)
    sort!(snp)

    return hgt,snp
end

function collectentry(data, type, name, accumulate!, post!)
    hgt, snp = params(data)
    val = [ zero(type) for i in 1:length(hgt), j in 1:length(snp) ]
    num = zeros(Int,length(hgt),length(snp))
    for key in keys(data)
        param = unpack(key)
        param.elt == name || continue

        i = findfirst(hgt .== param.hgt)
        j = findfirst(snp .== param.snp)

        data[key] !== nothing   || continue
        !any(isnan.(data[key])) || continue

        accumulate!(val,i,j,data[key])
        num[i,j] += 1
    end
    return post!(val, num), hgt, snp
end

# specific getters
entropy(data)    = collectentry(data,Float64,"tiles",(arr,i,j,x)->arr[i,j]+=x, (arr,num)->arr./num)
accuracy(data)   = collectentry(data,Array{Float64,1},"costs",(arr,i,j,x)->append!(arr[i,j],x),(arr,num)->arr)
diversity(data)  = collectentry(data,Float64,"dists",(arr,i,j,x)->arr[i,j]+=x, (arr,num)->arr./num)
complexity(data) = collectentry(data,Float64,"nblks",(arr,i,j,x)->arr[i,j]+=x, (arr,num)->arr./num)

function plotgrid(data, x, y; group="", labels=true)
    heatmap(string.(x), [@sprintf("%.1E",5*Y) for Y in y], data';
        xlabel = labels ? "HGT rate / genome / generation" : "",
        ylabel = labels ? "pairwise diversity" : "",
        title  = length(group) > 0 ? "entropy $(group)" : "",
        legend = :none,
        clim   = (0,1.5),
    )
end

function plotcdfs(data, x, y; group="", fontsize=12, kwargs...)
    i = size(data,1)÷2
    c = cgrad(:matter, size(data,2), categorical=true)
    p = plot(;
        legend = :bottomright,
        xscale = :log10,
        xlabel = "breakpoint misplacement (bp)",
        ylabel = "fraction of breakpoints",
        title  = length(group) > 0 ? "accuracy $(group)" : "",
        xtickfontsize  = fontsize,
        ytickfontsize  = fontsize,
        xguidefontsize = fontsize,
        yguidefontsize = fontsize,
        legendfontsize = fontsize,
        kwargs...
    )

    for j in 1:size(data,2)
        cdfplot!(data[i,j].+1;
            linewidth = 1,
            color     = c[j],
            label     = @sprintf("%.1E", 5*y[j]),
         )
    end

    return p
end

function publication(data, x, y)
    i = size(data,1)÷2
    c = cgrad(:matter, size(data,2), categorical=true)

    fig  = CairoMakie.Figure(font="Latin Modern Math", fontsize=26)
    axis = CairoMakie.Axis(fig[1,1],
                xscale=CairoMakie.log10,
                xlabel="breakpoint misplacement (bp)",
                ylabel="fraction of breakpoints",
                xticks=([1, 10, 100, 1000], [L"10^0", L"10^1", L"10^2", L"10^3"]),
                yticks=CairoMakie.LinearTicks(5)
    )
    CairoMakie.hidespines!(axis, :t, :r)

    for j in 1:size(data,2)
        points = sort(reduce(vcat, data[i,j] for i in 1:size(data,1))).+1
        CairoMakie.lines!(axis, points , range(0,1,length(points)),
            color = c[j],
            label = @sprintf("%.1e", 5*y[j])
        )
    end

    CairoMakie.axislegend("Avg. pairwise diversity", position = :rb, nbanks=3)

    return fig
end

const re = r"-([0-9]+).jld2"
function group(path)
    return match(re, path)[1]
end

function main(path, destdir)
    data = load(path)
    grid = plotgrid(entropy(data)...;  group=group(path))
    cdfs = plotcdfs(accuracy(data)...; group=group(path))

    save(base, name) = savefig("$(destdir)/$(base)-$(name).png")
    save(plot, base, name) = savefig(plot, "$(destdir)/$(base)-$(name).png")

    base(name) = basename(name)[1:end-5]

    save(grid,"heatmap",base(path))
    save(cdfs,"cdf",base(path))

    # publication plot
    fig = publication(accuracy(data)...)
    CairoMakie.save("$(destdir)/paper-$(base(path)).png", fig, px_per_unit=2)
#=
    TODO: return one figure with correct layout
    l = @layout[Plots.grid(1,2) a{0.05w}]
    return (
        plot(grid...,
            heatmap((0:0.01:1).*ones(101,1);
                    legend=:none,
                    xticks=:none,
                    yticks=(1:10:101, [@sprintf("%1.1f",1.5*x) for x in 0:0.1:1]),
             ); layout=l
        ),
        plot(cdfs...),
        plot(grid[1], cdfs[1])
    )
=#
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS[1], ARGS[2])
end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
function SHEMS_optimizer_sesu(sh, hp, fh, hw, b, m, bc_violations)
    flows = [:PV_DE, :B_DE, :GR_DE, :PV_B, :PV_GR, :PV_HP, :GR_HP, :B_HP, :HP_FH, :HP_HW];
    p_concr = 2400.0;   #kg/m^3
    c_concr = 1.0;      #kJ/(kg*°C)
    p_water = 997.0;    #kg/m^3
    c_water = 4.184;    #kJ/(kg*°C)

    # Input data__________________________________________________________________________________________________________________________________________________________
    df = CSV.read("data/200124_datafile_all_details_right_timestamp.csv");
    h_last = sh.h_start + m.h_predict -1;                     # optimization horizon

    # all peers have the same demand
    d_e =  df[sh.h_start:h_last,:electkwh];                     # electricity demand 1 year from BeOpt
    d_fh = df[sh.h_start:h_last,:heatingkwh];                  # heating demand 1 year from BeOpt
    d_hw = df[sh.h_start:h_last,:hotwaterkwh];                 # hot water demand 1 year from BeOpt
    g_e =  df[sh.h_start:h_last,:PV_generation];                # electricity generation (from Renewable ninjas)
    t_outside = df[sh.h_start:h_last,:Temperature];            # °C (from Renewable ninjas)

    # Calculate coefficients of performance for every time period (1:h_predict)___________________________________________
    cop_fh = COPcalc(fh, t_outside);
    cop_hw = COPcalc(hw, t_outside);

    # Model start________________________________________________________________________________________________________________________________________________________
    #Define Model und Solver and settings
    if m.solver == "Cbc"
        model = Model(optimizer_with_attributes(() -> Cbc.Optimizer()));
        set_optimizer_attribute(model, "ratioGap", m.mip_gap);
        set_optimizer_attribute(model, "logLevel", m.output_flag);
    elseif m.solver == "Gurobi"
        model = Model(optimizer_with_attributes(() -> Gurobi.Optimizer(GUROBI_ENV)));
        set_optimizer_attribute(model, "MIPGap", m.mip_gap);
        set_optimizer_attribute(model, "Presolve", m.presolve_flag);
        set_optimizer_attribute(model, "OutputFlag", m.output_flag);
    end

    @variables(model, begin
        X[1:m.h_predict, 1:length(flows)] >= 0;    #1:PV_DE, 2:B_DE, 3:GR_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 7:GR_HP, 8:B_HP, 9:HP_FH, 10:HP_HW
        SOC_b[1:m.h_predict]  >= 0;
        B_switch[1:m.h_predict], Bin;
        T_fh[1:m.h_predict] >= 0; V_hw[1:m.h_predict] >= 0;
        Mod_fh[1:m.h_predict] >= 0; Mod_hw[1:m.h_predict]   >= 0;
        T_fh_plus[1:m.h_predict] >= 0; T_fh_minus[1:m.h_predict] >= 0; V_hw_plus[1:m.h_predict] >= 0; V_hw_minus[1:m.h_predict] >= 0;
        HP_switch[1:m.h_predict], Bin;
        Hot[1:m.h_predict], Bin;
    end)

    # Fix start SoCs
    fix.(SOC_b[1], sh.soc_b; force=true);
    fix.(T_fh[1],sh.soc_fh; force=true);
    fix.(V_hw[1],sh.soc_hw; force=true);

    #1:PV_DE, 2:B_DE, 3:GR_DE, 4:PV_B, 5:PV_GR, 6:PV_HP, 7:GR_HP, 8:B_HP, 9:HP_FH, 10:HP_HW
    # Objective function: maximize self-sufficiency _____________________________________________________________________
    @objective(model, Min, sum(X[h,3] + X[h,7] +
        (sh.costfactor *(T_fh_plus[h] +T_fh_minus[h] +V_hw_plus[h] +V_hw_minus[h])) for h=1:m.h_predict));

    # Self-consumption: Comfort violations according to base case
    @constraint(model, SeSu_comfort, sum(T_fh_plus[h] + T_fh_minus[h] + V_hw_plus[h] + V_hw_minus[h] for h in 1:m.h_predict) <= bc_violations);

    # Electricity demand, generation, market clearing___________________________________
    @constraints(model, begin
        [h=1:m.h_predict],     sum(X[h,i] for i in 1:3) == d_e[h];                          # fulfill energy demand
        [h=1:m.h_predict],     sum(X[h,i] for i=[1,4,5,6]) == g_e[h];                       # restricted by PV generation
    end)
    # Battery__________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1],   SOC_b[h+1] == ((1 -b.loss) *SOC_b[h])+
                                    (b.eta *X[h,4]) -sum((1.0/(b.eta)) *X[h,i] for i=[2,8]);  # State of Charge, loss for unique solutions
        [h=1:m.h_predict],     b.soc_min <= SOC_b[h] <= b.soc_max;                          # Limits Battery usable capacity
        [h=1:m.h_predict],     X[h,4] <= B_switch[h]*b.rate_max;                              # Either charging
        [h=1:m.h_predict],     sum(X[h,i] for i=[2,8]) <= (1-B_switch[h])*b.rate_max;         # or discharging
    end)
    # Heat pump_______________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict],     sum(X[h,i] for i in 9:10) == sum(X[h,i] for i in 6:8);       # level power heat pump=and out
        [h=1:m.h_predict],     X[h,9] == Mod_fh[h] *hp.rate_max;                            # heating energy FH max cap
        [h=1:m.h_predict],     X[h,10] == Mod_hw[h] *hp.rate_max;                           # heating energy HW max cap
        [h=1:m.h_predict],     Mod_fh[h] <= 1 -HP_switch[h];                                # switch FH or HW
        [h=1:m.h_predict],     Mod_hw[h] <= HP_switch[h];                                   # switch HW or FH
    end)
    # Floor_heating__________________________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1],   T_fh[h+1] == T_fh[h]+
                                    (60*60)/(p_concr *fh.volume *c_concr) *
                                    ((cop_fh[h] *X[h,9]) -d_fh[h] -
                                    ((1-Hot[h]) *fh.loss) +(Hot[h] *fh.loss));              # SoC floor heating (temperature)
        [h=1:m.h_predict],     T_fh[h] - ((1 -Hot[h]) *m.big) <= t_outside[h];              # force hot binary on if hotter outside than inside
        [h=1:m.h_predict],     t_outside[h] -(Hot[h] *m.big) <= T_fh[h];
        [h=1:m.h_predict],     T_fh[h] <= fh.soc_max +T_fh_plus[h];                         # Limits temperature FH max
        [h=1:m.h_predict],     fh.soc_min -T_fh_minus[h] <= T_fh[h];                        # Limits temperature FH min
    end)
    #Hot water_____________________________________________________________________________________________
    @constraints(model, begin
        [h=1:m.h_predict-1],   V_hw[h+1] == V_hw[h]+
                                    (60*60)/((p_water *hw.t_supply *c_water)/1000) *
                                    ((cop_hw[h] *X[h,10]) -d_hw[h] -hw.loss);               # SoC hot water (volume)
        [h=1:m.h_predict],     V_hw[h] <= hw.soc_max +V_hw_plus[h];                         # Limits volume HW max
        [h=1:m.h_predict],     hw.soc_min - V_hw_minus[h] <= V_hw[h];                       # Limits volume HW min
    end)

    JuMP.optimize!(model);

    # collect returns
    profits = (sh.p_sell .*JuMP.value.(X[1:m.h_control,5])) .-sh.p_buy .*(JuMP.value.(X[1:m.h_control,3]) .+ JuMP.value.(X[1:m.h_control,7]));

    results = hcat(JuMP.value.(T_fh[1:m.h_control]), JuMP.value.(V_hw[1:m.h_control]), JuMP.value.(SOC_b[1:m.h_control]),
                    JuMP.value.(V_hw_plus[1:m.h_control]), JuMP.value.(V_hw_minus[1:m.h_control]),
                    JuMP.value.(T_fh_plus[1:m.h_control]), JuMP.value.(T_fh_minus[1:m.h_control]),
                    profits, cop_fh[1:m.h_control], cop_hw[1:m.h_control],
                    JuMP.value.(X[1:m.h_control, :]),                                               # flow variables
                    Matrix(df[sh.h_start:(sh.h_start + m.h_control -1), 13:15]));                   # month + day + hour

    if m.rolling_flag==1
        return JuMP.value.(SOC_b[m.h_control+1]), JuMP.value.(T_fh[m.h_control+1]), JuMP.value.(V_hw[m.h_control+1]), results
    else
        return results
    end
end


text:
using Pkg
Pkg.activate(@__DIR__)

using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--tests"
        default = false
        arg_type = Bool
    "--benchmarks"
        default = true
        arg_type = Bool
    "--storage"
        range_tester = s -> s in ("Array","CuArray")
        default = "Array"
    "--benchmark_accuracy"
        default = 0.01
        arg_type = Float64
    "-N"
        default = 256
        arg_type = Int
    "-T"
        default = Float32
        range_tester = s -> s in (Float32,Float64)
        eval_arg = true
end
args = parse_args(ARGS, s)

if args["storage"]=="CuArray"
    using CuArrays
end
storage = eval(Symbol(args["storage"]))
macro belapsed(ex)
    if args["storage"] == "CuArray"
        esc(:(BenchmarkTools.@belapsed CuArrays.@sync $ex))
    else
        esc(:(BenchmarkTools.@belapsed $ex))
    end
end
N = args["N"]
T = args["T"]

##

using CMBLensing
using CMBLensing: adapt

##

import BenchmarkTools
using Base.Broadcast: broadcastable
using Crayons
using DelimitedFiles
using LibGit2
using LinearAlgebra
using PrettyTables
using Printf
using Test
using Zygote

##

# benchmark to 1% and check to 5%, so only a "5σ event" should trigger a failure
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = args["benchmark_accuracy"]
rtol = 0.05
≲(a,b; rtol=rtol) = (a < (1+rtol)*b)

##

if args["tests"]
try

##

@testset "Performance Tests" begin

# spin 0
f = adapt(storage,FlatMap(rand(T,N,N)))
g = adapt(storage,FlatMap(rand(T,N,N)))
Ðf = Ð(f)
@test @belapsed(@.      $f + $f) ≲ @belapsed(@.         $f.Ix + $f.Ix)
@test @belapsed(@. $g = $f + $f) ≲ @belapsed(@. $g.Ix = $f.Ix + $f.Ix)
@test @belapsed(∇[1]*$Ðf)        ≲ @belapsed($(broadcastable(typeof(Ðf), ∇[1].diag)) .* $Ðf.Il)

# spin 2
ft = adapt(storage,FlatQUMap(rand(T,N,N),rand(T,N,N)))
gt = adapt(storage,FlatQUMap(rand(T,N,N),rand(T,N,N)))
@test @belapsed(@.       $ft + $ft) ≲ @belapsed(@. (         $ft.Qx + $ft.Qx,          $ft.Ux + $ft.Ux))
@test @belapsed(@. $gt = $ft + $ft) ≲ @belapsed(@. ($gt.Qx = $ft.Qx + $ft.Qx, $gt.Ux = $ft.Ux + $ft.Ux))

end

##

catch
end
end

##

if args["benchmarks"]

Cℓ = camb()
timing = []

for (s,pol) in [(0,:I),(2,:P)]

    @unpack f,f̃,ϕ,ds,ds₀ = load_sim(θpix=3, Nside=N, T=T, pol=pol, storage=storage, Cℓ=Cℓ);
    f°,ϕ° = mix(f,ϕ,ds₀)
    Lϕ = cache(LenseFlow(ϕ),f)

    append!(timing,[
        "Spin-$s Cache L" => @belapsed cache(LenseFlow($ϕ), $f);
        "Spin-$s L"       => @belapsed $Lϕ  * $f;
        "Spin-$s L†"      => @belapsed $Lϕ' * $f;
        "Spin-$s (∇L)†"   => @belapsed gradient(ϕ->norm($Lϕ(ϕ)*$f), $ϕ);
        "Spin-$s lnP"     => @belapsed                 lnP(:mix, $f°, $ϕ°, $ds₀);
        "Spin-$s ∇lnP"    => @belapsed gradient((f,ϕ)->lnP(:mix,  f,   ϕ,  $ds₀), $f°, $ϕ°);
    ])

end

meta = [
    "COMMIT" => string(LibGit2.head_oid(GitRepo(joinpath(dirname(@__FILE__),".."))))[1:7],
    "JULIA_NUM_THREADS" => Base.Threads.nthreads(),
    "FFTW_NUM_THREADS" => CMBLensing.FFTW_NUM_THREADS[],
    (k=>args[k] for k in ["storage","N","T"])...
]

reference_timing = Dict(
    "Spin-0 Cache L" => 25,
    "Spin-0 L"       => 13,
    "Spin-0 L†"      => 13,
    "Spin-0 (∇L)†"   => 85,
    "Spin-0 lnP"     => 65,
    "Spin-0 ∇lnP"    => 240,
    "Spin-2 Cache L" => 25,
    "Spin-2 L"       => 30,
    "Spin-2 L†"      => 30,
    "Spin-2 (∇L)†"   => 140,
    "Spin-2 lnP"     => 110,
    "Spin-2 ∇lnP"    => 380
)

# print benchmarks

pretty_table(Dict(meta), crop=:none)

pretty_table(
    vcat(([k v reference_timing[k]*1e-3] for (k,v) in timing)...),
    ["Operation","Time","Reference"],
    formatter = Dict(
        1 => (v,_) -> v,
        2 => (v,_) -> @sprintf("%.1f ms",1000v),
        3 => (v,_) -> @sprintf("%.0f ms",1000v)
    ),
    highlighters = (
        Highlighter((data, i, j) -> j==2 && data[i,j] > (1+rtol) * data[i,j+1], foreground=:red),
        Highlighter((data, i, j) -> j==2 && data[i,j] < (1-rtol) * data[i,j+1], foreground=:green)
    ),
    crop=:none,
    alignment=[:l,:r,:r],
    hlines=[6]
)

# save benchmarks
filename = "benchmarks/"*join(["$(k)_$(v)" for (k,v) in meta], "__")*".txt"
!ispath("benchmarks") && mkdir("benchmarks")
open(filename,"w") do f
    writedlm.(Ref(f),(meta,timing))
end

end


text:
using Crayons
# using ColorTypes

function print_density(som::Som)
    population_to_array(som) |> print_greyscale_rectangle
end

function print_density(som::Som,train,neighbors=4)
    train_density_to_array(som::Som,train,neighbors) |> print_greyscale_rectangle
end

function print_colored_text(som::Som,train,row_labels)
    neighbors,_ = nearest_neighbors(som,train,1)
    neighbors = map(x->x[1],neighbors)
    labels = map(x-> row_labels[x] ,neighbors)
    name_array = neuron_attribute_to_array(som,labels)
    color_array = train_density_to_array(som,train)
    print_greyscale_text(name_array,color_array)
end

function print_colored_text_uniq(som::Som,train,row_labels)
    included = []
    neighbors,_ = nearest_neighbors(som,train,10)
    neighbors = map(neighbors) do set
        # print(set)
        i = popfirst!(set)
        while ((i in included) && isempty(set))
            i = popfirst!(set)
        end
        push!(included,i)
        # print(i)
        i
    end
    # neighbors = map(x->x[1],neighbors)
    labels = map(x-> row_labels[x] ,neighbors)
    name_array = neuron_attribute_to_array(som,labels)
    color_array = train_density_to_array(som,train)
    print_greyscale_text(name_array,color_array)
end

"""
Attemps to map each training data row onto a unique neuron.
Row labels most be ordered in order of preference (higher up more likely to
be assigned to nearest neighbor neuron)
"""
function one_to_one_mapping(som::Som,train,row_labels)
    mapping = fill("",som.nCodes)
    neighbors, distances = nearest_neighbors(som::Som,train,16)
    for pair in zip(row_labels,distances)

    end
end

#
#

function train_density_to_array(som::Som,train,neighbors=4)
    densities = train_density(som,train,neighbors)
    neuron_attribute_to_array(som::Som,densities)
end

function population_to_array(som::Som)
    neuron_attribute_to_array(som,som.population)
end

function neuron_attribute_to_array(som::Som,attribute::Array{String})
    xdim, ydim = som.xdim, som.ydim
    a = fill("",xdim,ydim)
    _neuron_attribute_to_array(som::Som,attribute,a)
end

function neuron_attribute_to_array(som::Som,attributes)
    xdim, ydim = som.xdim, som.ydim
    a = zeros(xdim,ydim)
    _neuron_attribute_to_array(som::Som,attributes,a)
end

function _neuron_attribute_to_array(som::Som,attributes,a)
    for n in 1:som.nCodes
        x,y= som.indices[n,1:2]
        a[x,y] = attributes[n]
    end
    return a
end

function print_greyscale_rectangle(data,io::IO = stdout)

    f(x)= x*255 |> floor |> x->convert(Int,x)
    data = normalize01(data)
    data = map(f,data)
    width, height = size(data)

    str = "██"
    for i in 1:width
        for j in 1:height
            value = data[i,j]
            print(io, Crayon(; foreground = (value, value, value)), str, Crayon(reset = true))
       end
       print("\n")
    end

end

function print_greyscale_text(text_array,data,io::IO = stdout)

    f(x)= x*255 |> floor |> x->convert(Int,x)
    data = normalize01(data)
    data = map(f,data)
    width, height = size(data)

    for i in 1:width
        for j in 1:height
            str = text_array[i,j]
            value = data[i,j]
            print(io, Crayon(; foreground = (value, value, value)), str*" ", Crayon(reset = true))
       end
       print("\n")
    end

end

function normalize01(data)
    min,max = minimum(data),maximum(data)
    factor = max-min
    f(x)=(x-min)/factor
    map(f,data)
end