text:
abstract type GOCInput <: AbstractInput end

##TODO
function get_baseMVA(bus::String)
    return 100
end

"""
    read_input(input_type::T, instance_path::String) where T<:Type{GOCInput}

Read instance in `instance_path` depending on `input_type`.\n
Return a structure OPFProblems.\n
For `GOCInput`, read files pscopf_data.mat, contingency.csv and generator.csv, store them in arrays and use information to fulfill a DataSource structure for the basecase scenario. Finally take into account contingencies.

"""
function read_input(input_type::T, instance_path::String) where T<:Type{GOCInput}
    #read
    power_data,generator_data,contingency_data = getdataGOC(instance_path)
    ##info in bus
    bus,index = read_data_bus(power_data)
    ##info in generator
    gendata_dict = generator_data_to_dict(generator_data,index)
    ##info in generator
    bus = add_generator_data!(power_data,gendata_dict,bus,index)
    ##info in branch
    link = read_branch_data(power_data,index)
    #info basecase
    ds = DataSource(bus,link)
    node_linksin, node_linksout = SortedDict{String, SortedSet{Link}}(), SortedDict{String, SortedSet{Link}}()
    node_vars = SortedDict{String, SortedDict{String, Variable}}()
    link_vars = SortedDict{Link, SortedDict{String, Variable}}()
    gs = GridStructure("BaseCase", node_linksin, node_linksout)
    node_formulations = SortedDict{String, SortedDict{String, Symbol}}()
    link_formulations = SortedDict{Link, SortedDict{String, Symbol}}()
    mp = MathematicalProgramming(node_formulations, link_formulations, node_vars,link_vars)
    ##read scenarios
    OPFproblems = scenarios_data(ds, gs, mp, contingency_data,index)
    return OPFproblems
end

#########################
## Utils functions
#########################
"""
    getdataGOC(instance_path)

Read 3 files (pscopf_data.mat,generator.csv,contingency.csv) in `instance_path` and returns 3 arrays of data (power_data,generator_data,contingency_data)
Return an error if one the three files does not exist in path_folder.
"""
function getdataGOC(instance_path) ##String instance+scenario
    power_data = getpowerdata(instance_path)
    generator_data = getgeneratordata(instance_path)
    contingency_data = getcontingencydata(instance_path)
    return power_data,generator_data,contingency_data
end

function getpowerdata(instance_path)
    path_folder = instance_path
    if !("pscopf_data.mat" ∈ readdir(path_folder))
        error("missing file : pscopf_data.mat must be in $path_folder")
    end
    try
        return read(matopen(joinpath(path_folder,"pscopf_data.mat")))
    catch
        warn("Not possible to read matopen(pscopf_data.mat)")
        return
    end
end

function getgeneratordata(instance_path)
    path_folder = instance_path
    if !("generator.csv" ∈ readdir(path_folder))
        error("missing file : generator.csv must be in $path_folder")
    end
    try
        return readdlm(joinpath(path_folder,"generator.csv"),',',skipstart=1)
    catch
        warn("Not possible to readdlm generator.csv")
        return
    end

end

function getcontingencydata(instance_path)
    path_folder = instance_path
    if !("contingency.csv" ∈ readdir(path_folder))
        error("missing file : contingency.csv must be in $path_folder")
    end
    try
        return readdlm(joinpath(path_folder,"contingency.csv"),',',skipstart=1)
    catch
        warn("Not possible to readdlm contingency.csv")
    end
end

#########################
function get_bus_data(bus_data, listofkeys, id)
    listofvalues = []
    for key in listofkeys
        if !haskey(bus_data, key)
            error("bus_data dictionary does not have a key $key")
        else
            append!(listofvalues, bus_data[key][id])
        end
    end
    return listofvalues
end

function get_bus_index(power_data)
    nb_bus = Int(power_data["totalbus"])
    all_busdata_key = filter(x->ismatch(r"all_busdata", x), collect(keys(power_data)))[1]
    raw_bus_data = power_data[all_busdata_key]
    nb_lines = size(raw_bus_data,1)
    if nb_lines!=nb_bus
        error("number of lines of all_busdata not equal to number of bus in power_data at key totalbus")
    end
    index = SortedDict{Int64,Int64}()
    for line in 1:nb_lines
        id_bus = Int64(raw_bus_data[line,1])
        index[id_bus] = id_bus
    end
    return index
end

"""
    read_data_bus(power_data)

Reads data in `power_data` to construct the DataSource field bus without generator data (only VOLTM, load and shunt data)
"""
function read_data_bus(power_data)
    nb_bus = Int(power_data["totalbus"])
    index = get_bus_index(power_data)
    bus_data = power_data["bus"]
    # all_busName_key = filter(x->ismatch(r"all_busName", x), collect(keys(power_data)))[1]
    # busname_data = power_data[all_busName_key]
    all_busdata_key = filter(x->ismatch(r"all_busdata", x), collect(keys(power_data)))[1]
    raw_bus_data = power_data[all_busdata_key]
    nb_lines = size(raw_bus_data,1)
    if nb_lines!=length(bus_data["BaseKV"])
        error("sizes of all_busName and bus not equal")
    end
    bus = SortedDict{String, SortedDict{String,Any}}()
    # bus = SortedDict(bus_name(line) => SortedDict{String,Any}() for line in 1:nb_bus)
    for line in 1:nb_lines
        id_bus = Int(raw_bus_data[line,1])
        # index[id_bus] = line
        busname = bus_name(index[id_bus])
        baseKV = bus_data["BaseKV"][line]
        baseMVA = 1.0
        voltage_magnitude_min , voltage_magnitude_max  = get_bus_data(bus_data, ["VoltageMagnitudeMin", "VoltageMagnitudeMax"], line)
        if !haskey(bus, busname)
            bus[busname] = SortedDict{String,Any}()
        end
        bus[busname][volt_name()] = GOCVolt(id_bus, baseKV, baseMVA, voltage_magnitude_min, voltage_magnitude_max)
    end

    all_loadID_key = filter(x->ismatch(r"all_loadID", x), collect(keys(power_data)))[1]
    load_id = power_data[all_loadID_key]
    all_loaddata_key = filter(x->ismatch(r"all_loaddata", x), collect(keys(power_data)))[1]
    load_data = power_data[all_loaddata_key]

    for line in 1:size(load_data,1)
        id_bus = Int(load_data[line,1])
        busname = bus_name(index[id_bus])
        load = load_data[line,5] + im * load_data[line,6]
        id_load = id_elem(load_id,line)
        loadname = load_name(id_load)
        bus[busname][loadname] = GOCLoad(id_bus,loadname,load)
    end

    all_fixedshuntID_key = filter(x->ismatch(r"all_fixedshuntID", x), collect(keys(power_data)))[1]
    shunt_id = power_data[all_fixedshuntID_key]
    all_fixedshuntdata_key = filter(x->ismatch(r"all_fixedshuntdata", x), collect(keys(power_data)))[1]
    shunt_data = power_data[all_fixedshuntdata_key]
    for line in 1:size(shunt_data,1)
        id_bus = Int(shunt_data[line,1])
        busname = bus_name(index[id_bus])
        shunt = shunt_data[line,3] + im * shunt_data[line,4]
        id_shunt = id_elem(shunt_id,line)
        shuntname = shunt_name(id_shunt)
        bus[busname][shuntname] = GOCShunt(id_bus,shuntname,shunt)
    end

    return bus,index
end

function id_elem(array_id,line)
    if typeof(array_id) == String
        temp_string = array_id
    else
        temp_string = array_id[line]
        if typeof(temp_string)==Float64
            temp_string = Int(temp_string)
        end
    end
    return temp_string
    #return matchall(r"\S+" ,temp_string)[1]
end

#########################
function get_generator_data(generator_data, listofkeys, id)
    listofvalues = []
    for key in listofkeys
        if !haskey(generator_data, key)
            error("generator_data dictionary does not have a key $key")
        else
            append!(listofvalues, generator_data[key][id])
        end
    end
    return listofvalues
end

function generator_data_to_dict(generator_data,index) ##conversion array to dict
    generator_data_dict = SortedDict{String, SortedDict{String,SortedDict{Int64,Float64}}}()
    for line in 1:size(generator_data,1)
        bus_id = Int(generator_data[line,1])
        busname = bus_name(index[bus_id])
        gen_id = generator_data[line,2]
        gen_id = remove_simple_quotes_if_present(gen_id)
        if typeof(gen_id)==Float64
            gen_id = Int(gen_id)
        end
        genname = generator_name(gen_id)
        if !haskey(generator_data_dict, busname)
            generator_data_dict[busname] = SortedDict{String,SortedDict{Int64,Float64}}()
        end
        if !haskey(generator_data_dict[busname],genname)
            generator_data_dict[busname][genname] = SortedDict{Int64,Float64}()
        end
        term_id = Int(generator_data[line,3])
        value = generator_data[line,4]
        if term_id!=9
            generator_data_dict[busname][genname][term_id] = value
        end
    end
    return generator_data_dict
end

function remove_simple_quotes_if_present(id)
    if typeof(id)==SubString{String}
        if id[end-1]==' '
            return id[2:(end-2)]
        else
            return id[2:(end-1)]
        end
    else
        if id[end]==' '
            return id[1:(end-1)]
        else
            return id
        end
    end
end

"""
    add_generator_data!(power_data,generator_data_dict,bus,index)

Add generator data contained in `power_data` in DataSource field `bus` depending on `generator_data_dict` and `index`.\n
Return `bus`.

"""
function add_generator_data!(power_data,generator_data_dict,bus,index)
    nb_bus = Int(power_data["totalbus"])
    generator_data = power_data["generator"]
    for gen in 1:length(generator_data["BusNum"])
        bus_id = Int(get_generator_data(generator_data,["BusNum"],gen)[1])
        # all_busName_key = filter(x->ismatch(r"all_busName", x), collect(keys(power_data)))[1]
        # busname_data = power_data[all_busName_key]
        busname = bus_name(index[bus_id])
        gen_id, Pmin, Qmin, Pmax, Qmax, participation_factor =  get_generator_data(generator_data, ["UnitID","RealPowerMin","ReactivePowerMin","RealPowerMax","ReactivePowerMax", "ParticipationFactor"], gen)
        power_min = Pmin + im*Qmin
        power_max = Pmax + im*Qmax
        # cost_degrees = generator_data["RealPowerCostExponent"][gen,:]
        # cost_coeffs = generator_data["RealPowerCostCoefficient"][gen,:]
        gen_id2 = gen_id
        gen_id2 = remove_simple_quotes_if_present(gen_id2)
        if typeof(gen_id2)==Float64
            gen_id2 = Int(gen_id2)
        end
        generatorname = generator_name(gen_id2)
        dict_obj_coeffs = SortedDict{Int64,Float64}()
        for (degree, value) in generator_data_dict[busname][generatorname]
            if (degree ∈ [0,1,2])
                dict_obj_coeffs[degree] = value
            end
        end
        if length(dict_obj_coeffs)<3
            #warn("A 2-order polynom cost with less than 3 values for generator $generatorname at bus $busname, default value 0")
            if !haskey(dict_obj_coeffs,0) dict_obj_coeffs[0] = 0 end
            if !haskey(dict_obj_coeffs,1) dict_obj_coeffs[1] = 0 end
            if !haskey(dict_obj_coeffs,2) dict_obj_coeffs[2] = 0 end
       end
        bus[busname][generatorname] = GOCGenerator(bus_id,gen_id,power_min,power_max,participation_factor,dict_obj_coeffs)
    end
    return bus
end

#########################
function get_branch_data(branch_data, listofkeys, id)
    listofvalues = []
    for key in listofkeys
        if !haskey(branch_data, key)
            error("branch_data dictionary does not have a key $key")
        else
            append!(listofvalues, branch_data[key][id])
        end
    end
    return listofvalues
end

"""
    read_branch_data(power_data, index)

Read branch data contained in `power_data` to fulfill DataSource field `link` depending on `index`. \n
Return dictionary `link`.
"""
function read_branch_data(power_data, index)
    nb_branch = Int(power_data["totalbranch"])
    nb_bus = Int(power_data["totalbus"])
    # all_busName_key = filter(x->ismatch(r"all_busName", x), collect(keys(power_data)))[1]
    # busname_data = power_data[all_busName_key]
    all_lineCKT_key = filter(x->ismatch(r"all_lineCKT", x), collect(keys(power_data)))[1]
    lines_data = power_data[all_lineCKT_key]
    nb_lines = length(lines_data)
    link = SortedDict{Link, SortedDict{String,Any}}()
    branch_data = power_data["branch"]

    #lines without transformer
    for branch in 1:nb_lines
        origin = Int(get_branch_data(branch_data, ["Origin"], branch)[1])
        destination = Int(get_branch_data(branch_data, ["Destination"], branch)[1])
        linkname = Link(bus_name(index[origin]),bus_name(index[destination]))
        br_id = lines_data[branch]
        br_id2 = br_id
        branch_id = remove_simple_quotes_if_present(br_id2)
        if typeof(branch_id)==Float64
            branch_id = Int(branch_id)
        end
        resistance, reactance, susceptance, power_magnitude_max = get_branch_data(branch_data, ["SeriesResistance","SeriesReactance","ChargingSusceptance","PowerMagnitudeMax"], branch)
        if !haskey(link, linkname)
            link[linkname] = SortedDict{String,Any}()
        end
        if resistance == reactance == 0
            println(linkname, " nullimpedance line without transformer")
            name_line = nullimpedance_notransformer_name(branch_id)
            link[linkname][name_line] = GOCNullImpedance_notransformer(origin,destination,br_id,susceptance, power_magnitude_max)
        else
            name_line = linepi_notransformer_name(branch_id)
            link[linkname][name_line] = GOCLineπ_notransformer(origin,destination,br_id,resistance,reactance,susceptance, power_magnitude_max)
        end
    end

    all_transformerCKT_key = filter(x->ismatch(r"all_transformerCKT", x), collect(keys(power_data)))[1]
    transformername_data = power_data[all_transformerCKT_key]

    #lines with transformers
    for branch in (nb_lines+1):nb_branch
        origin = Int(get_branch_data(branch_data, ["Origin"], branch)[1])
        destination = Int(get_branch_data(branch_data, ["Destination"], branch)[1])
        linkname = Link(bus_name(index[origin]),bus_name(index[destination]))
        br_id = transformername_data[branch-nb_lines]
        br_id2 = br_id
        branch_id = remove_simple_quotes_if_present(br_id2)
        if typeof(branch_id)==Float64
            branch_id = Int(branch_id)
        end
        resistance, reactance, susceptance, power_magnitude_max = get_branch_data(branch_data, ["SeriesResistance","SeriesReactance","ChargingSusceptance","PowerMagnitudeMax"], branch)
        transfo_ratio, transfo_phase = get_branch_data(branch_data, ["TapRatio", "PhaseShift"], branch)
        if !haskey(link, linkname)
            link[linkname] = SortedDict{String,Any}()
        end
        if resistance == reactance == 0
            println(linkname, " nullimpedance line with transformer")
            name_line = nullimpedance_withtransformer_name(branch_id)
            link[linkname][name_line] = GOCNullImpedance_withtransformer(origin,destination,br_id,susceptance, transfo_ratio,transfo_phase, power_magnitude_max)
        else
            name_line = linepi_withtransformer_name(branch_id)
            link[linkname][name_line] = GOCLineπ_withtransformer(origin,destination,br_id,resistance,reactance,susceptance, transfo_ratio,transfo_phase, power_magnitude_max)
        end
    end

    return link
end

#########################
"""
    scenarios_data(ds,gs,mp,power_data, contingency_data,index)

Use dictionaries `ds`, `gs` and `mp` to create `Scenario` structures for each contingency contained in `contingency_data`.\n
Return an `OPFproblems` structure : "scenario" => Scenario

"""
function scenarios_data(ds,gs,mp,contingency_data,index)
    output = SortedDict("BaseCase" => Scenario(ds,gs,mp))
    ds = output["BaseCase"].ds
    nb_contingencies = size(contingency_data,1)
    for ct in 1:nb_contingencies
        id_contingency = contingency_data[ct,1]
        scenario_name = scenarioname(id_contingency)
        type_contingency = contingency_data[ct,2]
        #B, T or G (Branch, Transformer, Generator)
        gs_scenario = copy(gs)
        gs_scenario.scenario = scenario_name
        ds_scenario_bus = SortedDict{String, SortedDict{String,Any}}()
        for busname in keys(ds.bus)
            ds_scenario_bus[busname] = SortedDict{String,Any}()
        end
        ds_scenario_link = SortedDict{Link, SortedDict{String,Any}}()
        for linkname in keys(ds.link)
            ds_scenario_link[linkname] = SortedDict{String,Any}()
        end
        ds_scenario = DataSource(ds_scenario_bus,ds_scenario_link)
        if type_contingency == "G"
            println(type_contingency)
            bus_id = contingency_data[ct,3]
            busname = bus_name(index[bus_id])
            gen_num = contingency_data[ct,4]
            CID = contingency_data[ct,5]
            if typeof(CID)==SubString{String}
                gen_ID = CID[2:(end-1)]
            else
                gen_ID = string(CID)
            end
            for (bus,dict) in ds.bus
                if bus==busname
                    for (buslabel,data) in dict
                        if buslabel != generator_name(gen_ID)
                            ds_scenario.bus[bus][buslabel] = data
                        end
                    end
                else
                    for (buslabel,data) in dict
                        ds_scenario.bus[bus][buslabel] = data
                    end
                end
            end
            for (linkname,dict) in ds_scenario.link
                for (linklabel,data) in ds.link[linkname]
                    ds_scenario.link[linkname][linklabel] = ds.link[linkname][linklabel]
                end
            end
        elseif type_contingency=="B" || type_contingency=="T"
            #println("Type of contingency $ct is $type_contingency")
            origin = Int(contingency_data[ct,3])
            destination = Int(contingency_data[ct,4])
            link_to_remove = Link(bus_name(index[origin]),bus_name(index[destination]))
            # println("link_to_remove : ", link_to_remove)
            CID = contingency_data[ct,5]
            if typeof(CID)==SubString{String}
                id_line = CID[2:(end-1)]
            else
                id_line = string(CID)
            end
            for (busname,dict) in ds.bus
                for (buslabel,data) in dict
                    ds_scenario.bus[busname][buslabel] = data
                end
            end
            for (linkname,dict) in ds.link
                if linkname.orig == link_to_remove.orig && linkname.dest == link_to_remove.dest
                    # println("link to remove : ", link_to_remove)
                    for (linklabel,data) in dict
                        if !contains(linklabel,id_line)
                            ds_scenario.link[linkname][linklabel] = data
                        end
                    end
                    if isempty(ds_scenario.link[linkname])
                        delete!(ds_scenario.link,linkname)
                    end
                else
                    for (linklabel,data) in dict
                    ds_scenario.link[linkname][linklabel] = data
                    end
                end
            end
        end
        output[scenario_name] = Scenario(ds_scenario,gs_scenario,mp)
    end

    return output
end


text:
using JuMP
using LightGraphs, SimpleWeightedGraphs
using Random

lazy_dcc_constraints = 0

function directed_cutset_constraints!(
    model,
    graph::SimpleWeightedGraph,
    k::Int)
    global lazy_dcc_constraints = 0
    @constraint(model, sum(variable_by_name(model, "y[1,$j]") for j in 2:nv(graph)) == 1)
    @constraint(model, sum(variable_by_name(model, "y[$j,1]") for j in 2:nv(graph)) >= 1)
    @constraint(model, sum(variable_by_name(model, "y[$j,1]") for j in 2:nv(graph)) <= k-1)
    for i in 2:nv(graph)
        @constraint(model, sum(variable_by_name(model, "y[$i,$j]")
                               for j in 1:nv(graph) if has_edge(graph, i, j))
                    >= sum(variable_by_name(model, "y[$j,$i]")
                           for j in 1:nv(graph) if has_edge(graph, j, i)))

        @constraint(model, sum(variable_by_name(model, "y[$i,$j]")
                               for j in 1:nv(graph) if has_edge(graph, i, j))
                    <= sum(variable_by_name(model, "y[$j,$i]")
                           for j in 1:nv(graph) if has_edge(graph, j, i)) * (k-1))
        # @constraint(model, sum(variable_by_name(model, "y[$i,$j]") for j in 1:nv(graph) if has_edge(graph, i,j)) <= k-1)
        @constraint(model, sum(variable_by_name(model, "y[$j,$i]") for j in 1:nv(graph) if has_edge(graph, i,j)) <= 1)
    end
    MOI.set(model,
            MOI.LazyConstraintCallback(),
            dcc_lazy_clause_generator(model, graph, k))
end

function dcc_lazy_clause_generator(model, graph::SimpleWeightedGraph, k::Int)
    es = collect(SimpleWeightedEdge, edges(graph))
    sort!(es, by=dst)
    sort!(es, by=src, alg=MergeSort)
    y::Dict{Tuple{Int,Int},VariableRef} = Dict{Tuple{Int,Int},VariableRef}()
    for e in es
        i::Int = src(e)
        j::Int = dst(e)
        y[i,j] = variable_by_name(model, "y[$i,$j]")
        y[j,i] = variable_by_name(model, "y[$j,$i]")
    end
    return function(cb_data)
        status = callback_node_status(cb_data, model)
        if status == MOI.CALLBACK_NODE_STATUS_FRACTIONAL
            @debug "Found fractional solution"
            return
        elseif status == MOI.CALLBACK_NODE_STATUS_INTEGER
            @debug "Found integer solution"
            y_val::Dict{Tuple{Int,Int},Int} = Dict{Tuple{Int,Int},Int}()
            for e in es
                i::Int = src(e)
                j::Int = dst(e)
                y_val[i,j] = round(callback_value(cb_data, y[i,j]))
                y_val[j,i] = round(callback_value(cb_data, y[j,i]))
            end
            cutset = find_cutset(graph, y_val, k)
            if isnothing(cutset)
                return
            end
            @debug "Found cutset of size $(length(cutset)): $cutset."
            complement::Set{Int} = Set{Int}(2:nv(graph))
            filter!(v -> !(v in cutset), complement)
            @assert length(complement) + length(cutset) == nv(graph) "Not all nodes in cutset and complement"
            filter!(v -> v != 1, cutset)
            constraint = @build_constraint(sum(y[i,j] for i in cutset for j in complement if has_edge(graph, i, j)) >= sum(y[1, i] for i in cutset)) # Constraint 39
            global lazy_dcc_constraints += 1
            MOI.submit(model, MOI.LazyConstraint(cb_data), constraint)
        end
    end
end

function find_cutset(
    graph::SimpleWeightedGraph,
    y_val::Dict{Tuple{Int,Int},Int64},
    k::Int)::Union{Nothing,Set{Int}}
    visited::Set{Int} = Set{Int}()
    queue::Set{Int} = Set{Int}([1])
    while !isempty(queue)
        @assert !isempty(queue) "queue should not be empty."
        c::Int = pop!(queue)
        push!(visited, c)
        ns = neighbors(graph, c)
        for n::Int in ns
            if y_val[c, n] == 1 && n != 1
                push!(queue, n)
            end
        end
    end
    if length(visited) == k+1
        return # no cutset needed
    end
    @assert length(visited) <= k "Solution should be a forest, but only has one tree."
    return visited

end

ext: jl

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
"""
Wrapper type to reduce the visual noise of Pluto's `Base.show` on `RTree`.
"""
struct RTreeWrapper
    t::RTree{Float64,3}
end

"""
Adapted from https://github.com/alyst/SpatialIndexing.jl/blob/105b420e/src/rtree/show.jl.
"""
function Base.show(io::IO, rtree::RTreeWrapper)
    tree = rtree.t
    print(io, typeof(tree))
    print(io, "(variant=")
    print(io, tree.variant)
    print(io, ", tight_mbrs=")
    print(io, tree.tight_mbrs)
    print(io, ", nearmin_overlap=")
    print(io, tree.nearmin_overlap)
    print(io, ", fill_factor=")
    print(io, tree.fill_factor)
    print(io, ", split_factor=")
    print(io, tree.split_factor)
    print(io, ", reinsert_factor=")
    print(io, tree.reinsert_factor)

    print(io, ", leaf_capacity=")
    print(io, SI.capacity(SI.Leaf, tree))
    print(io, ", branch_capacity=")
    print(io, SI.capacity(SI.Branch, tree))
    println(io, ")")
    print(io, SI.length(tree))
    print(io, " element(s) in ")
    print(io, SI.height(tree))
    print(io, " level(s) (")
    print(io, join(reverse(tree.nnodes_perlevel), ", "))
    println(io, " node(s) per level)")
end

"""
Geometric data structure that stores a tetrahedral mesh with potentials.
"""
struct Mesh
    n::Int32 # number of vertices
    m::Int32 # number of tets
    verts::Matrix{Float64} # shape: [3, n]
    tets::Matrix{Int64} # shape: [4, m]
    potentials::Vector{Float64} # shape: [n]
    com::Vector{Float64} # center of mass
    rtree::RTreeWrapper # bounding box index data structure

    Mesh(verts::Matrix{Float64}, tets::Matrix{Int64}, potentials::Vector{Float64}) = begin
        @assert(size(verts, 1) == 3, "verts should be in R3")
        @assert(size(tets, 1) == 4, "tets should be 4-tuples of indices")
        @assert(
            size(potentials, 1) == size(verts, 2),
            "potentials should be the same size as verts"
        )
        n, m = size(verts, 2), size(tets, 2)
        com = center_of_mass(verts, tets)
        rtree = RTree{Float64,3}(Int64)
        for i = 1:m
            points = verts[:, tets[:, i]]
            bb_min, bb_max = bounding_box(points)
            rect = SI.Rect(
                (bb_min[1], bb_min[2], bb_min[3]),
                (bb_max[1], bb_max[2], bb_max[3]),
            )
            SI.insert!(rtree, rect, i)
        end
        new(n, m, verts, tets, potentials, com, RTreeWrapper(rtree))
    end
end

"""
Compute center of mass, given coords of vertices and tets.
"""
function center_of_mass(verts::Matrix{Float64}, tets::Matrix{Int64})
    # verts [3,n]; tets [4, m]
    tet_centers = Vector{Float64}[]
    vols = Vector{Float64}()
    for i = 1:size(tets, 2)
        tet_vtxs = verts[:, tets[:, i]] # 3 x 4
        push!(tet_centers, mean(eachcol(tet_vtxs)))
        tet_vtxs_4 = [tet_vtxs; ones(Int, 1, 4)]
        append!(vols, abs(det(tet_vtxs_4) / 6))
    end
    sum(tet_centers .* vols) / sum(vols)
end

"""Compute the volume of a mesh."""
function volume(mesh::Mesh)
    values = mesh.verts[:, mesh.tets]
    volume = 0.0
    for tet_values in eachslice(values, dims = 3)
        a, b, c, d = eachcol(tet_values)
        # compute triple product
        volume += abs(dot(cross(b - a, c - a), d - a)) / 6.0
    end
    volume
end

"""Compute the lower and upper bounding box of a collection of points."""
function bounding_box(points::Matrix{Float64})
    min.(eachcol(points)...), max.(eachcol(points)...)
end

const Pose = SMatrix{4,4,Float64}

struct Object
    mesh::Mesh
    pose::Pose
    v::SVector{3,Float64}
    ω::SVector{3,Float64}

    """Construct a new object from a mesh, with optional pose and velocities."""
    Object(
        mesh::Mesh,
        pose::Pose = Pose(I),
        v::SVector{3,Float64} = (@SVector zeros(3)),
        ω::SVector{3,Float64} = (@SVector zeros(3)),
    ) = new(mesh, pose, v, ω)
end

"""Compute the volume of an object."""
function volume(obj::Object)
    volume(obj.mesh)
end

"""Applies transform to an object."""
function transform(obj::Object, matrix::Pose)::Object
    Object(obj.mesh, matrix * obj.pose, obj.v, obj.ω)
end

"""
Applies transform to vertices.
"""
function transform(vertices::Matrix{Float64}, matrix::Pose)
    # concatenate ones to vertices
    ones_arr = ones(1, size(vertices, 2))
    mat = vcat(vertices, ones_arr)
    vertices = matrix * mat
    return vertices[1:3, :]
end

"""
Applies transform to a single vertex.
"""
function transform(vertex::Vector{Float64}, matrix::Pose)
    return (matrix*[vertex; 1.0])[1:3]
end

function translate(obj::Object, disp::SVector{3})::Object
    transform(obj, @SMatrix [
        1 0 0 disp[1]
        0 1 0 disp[2]
        0 0 1 disp[3]
        0 0 0 1
    ])
end

function rotateX(obj::Object, θ::Float64)::Object
    transform(obj, @SMatrix [
        1 0 0 0
        0 cos(θ) -sin(θ) 0
        0 sin(θ) cos(θ) 0
        0 0 0 1
    ])
end

function rotateY(obj::Object, θ::Float64)::Object
    transform(obj, @SMatrix [
        cos(θ) 0 sin(θ) 0
        0 1 0 0
        -sin(θ) 0 cos(θ) 0
        0 0 0 1
    ])
end

function rotateZ(obj::Object, θ::Float64)::Object
    transform(obj, @SMatrix [
        cos(θ) -sin(θ) 0 0
        sin(θ) cos(θ) 0 0
        0 0 1 0
        0 0 0 1
    ])
end

"""
Updates the current object's pose according to an angular velocity.

Uses the exponential map on SO(3), see https://arwilliams.github.io/so3-exp.pdf.
"""
function changeRotation(obj::Object, ω::SVector{3,Float64})::Object
    θ = norm(ω)
    if θ < 1e-12
        return obj
    end

    # Compute the skew matrix.
    ωx = @SMatrix [
        0 -ω[3] ω[2]
        ω[3] 0 -ω[1]
        -ω[2] ω[1] 0
    ]

    # Compute the matrix exponential / exponential map.
    ω_exp = I + (sin(θ) / θ) .* ωx + ((1 - cos(θ)) / θ^2) .* (ωx * ωx)
    transform(obj, Pose([ω_exp zeros(3); zeros(3)' 1]))
end

export Mesh,
    Pose, Object, volume, transform, translate, rotateX, rotateY, rotateZ, changeRotation

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
using Test
import IJulia: InlineIOContext

@testset "Custom Jupyter inline display" begin
    @eval struct TestDataType
        payload
    end

    function Base.show(io::IO, m::MIME"text/plain", data::TestDataType)
        print(io, "TestDataType: ")
        if get(io, :jupyter, false)
            print(io, "Jupyter: ")
        end
        Base.show(io, m, data.payload)
    end

    data = TestDataType("foo")
    buf = IOBuffer()

    Base.show(buf, MIME("text/plain"), data)
    @test String(take!(buf)) == "TestDataType: \"foo\""

    Base.show(InlineIOContext(buf), MIME("text/plain"), data)
    @test String(take!(buf)) == "TestDataType: Jupyter: \"foo\""
end


text:
using Test, SimGBS

my_tests = ["gbs.jl", "genome.jl", "pop.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
import Flux.Tracker.back
using SparseArrays
import Base.*

function multrans(A::Matrix, B::SparseMatrixCSC)
  mA, nA = size(A)
  mB, nB = size(B)
  @boundscheck nA == nB || throw(DimensionMismatch())
  C = zeros(promote_type(eltype(A),eltype(B)), mA, mB)
  @inbounds for jB in 1:nB
      for kB in B.colptr[jB]:(B.colptr[jB+1] - 1)
          iB = B.rowval[kB]
          xB = B.nzval[kB]
          for iA in 1:mA
              C[iA, iB] += A[iA, jB] * xB
          end
      end
  end
  return C
end

function mul(A::Matrix, B::SparseMatrixCSC)
  mA, nA = size(A)
  mB, nB = size(B)
  @boundscheck nA == mB || throw(DimensionMismatch())
  C = zeros(promote_type(eltype(A),eltype(B)), mA, nB)
  @inbounds for jB in 1:nB
      for kB in B.colptr[jB]:(B.colptr[jB+1] - 1)
          iB = B.rowval[kB]
          xB = B.nzval[kB]
          for iA in 1:mA
              C[iA, jB] += A[iA, iB] * xB
          end
      end
  end
  return C
end

"""
  mul(A::Matrix, B::SparseMatrixCSC, mask::Vector{Int})
  mul(A::Matrix, B::SparseMatrixCSC, mask::BitArray)

  multiply transposed matrix `A` with `B` while assuming that only rows of `B` with indexes in `mask` are non-zeros
"""
function multrans(A::Matrix, B::SparseMatrixCSC, mask::BitArray)
  mA, nA = size(A)
  mB, nB = size(B)
  @boundscheck nA == nB || throw(DimensionMismatch())
  C = zeros(promote_type(eltype(A),eltype(B)), mA, mB)
  @inbounds for jB in 1:nB
      for kB in B.colptr[jB]:(B.colptr[jB+1] - 1)
          iB = B.rowval[kB]
          xB = B.nzval[kB]
          !mask[iB] && continue
          for iA in 1:mA
              C[iA, iB] += A[iA, jB] * xB
          end
      end
  end
  return C
end
multrans(A::Matrix, B::SparseMatrixCSC, mask::Vector{Int}) = (m = falses(size(A,2)); m[mask] .= true; multrans(A, B, m))

"""
  mul(A::Matrix, B::SparseMatrixCSC, mask::Vector{Int})
  mul(A::Matrix, B::SparseMatrixCSC, mask::BitArray)

  multiply the matrix `A` with `B` while assuming that only rows of `B`  with indexes in `mask` are non-zeros
"""
function mul(A::Matrix, B::SparseMatrixCSC, mask::BitArray)
  mA, nA = size(A)
  mB, nB = size(B)
  @boundscheck nA == mB || throw(DimensionMismatch())
  C = zeros(promote_type(eltype(A),eltype(B)), mA, nB)
  @inbounds for jB in 1:nB
      for kB in B.colptr[jB]:(B.colptr[jB+1] - 1)
          iB = B.rowval[kB]
          xB = B.nzval[kB]
          !mask[iB] && continue
          for iA in 1:mA
              C[iA, jB] += A[iA, iB] * xB
          end
      end
  end
  return C
end
mul(A::Matrix, B::SparseMatrixCSC, mask::Vector{Int}) = (m = falses(size(A,2)); m[mask] .= true; mul(A, B, m))

a::Flux.Tracker.TrackedMatrix * b::SparseMatrixCSC = Flux.Tracker.track(mul, a, b)
Flux.Tracker.@grad function mul(a::AbstractMatrix, b::SparseMatrixCSC)
  return mul(Flux.data(a),b) , Δ -> (multrans(Δ, b),nothing)
end

mul(a::Flux.Tracker.TrackedMatrix, b::SparseMatrixCSC, mask) = Flux.Tracker.track(mul, a, b, mask)
Flux.Tracker.@grad function mul(a::AbstractMatrix, b::SparseMatrixCSC, mask)
  return mul(Flux.data(a), b, mask) , Δ -> (multrans(Δ, b),nothing, nothing)
end

function Base.reduce(::typeof(hcat), A::AbstractVector{T}) where {T<: SparseMatrixCSC{Tv, Ti} where {Tv, Ti}}
  isempty(A) && return(A)
  length(A) == 1 && return(A[1])
  Tv, Ti = eltype(A[1].nzval), eltype(A[1].colptr)
  coloffset = 0
  Is, Js, Vs = Vector{Vector{Ti}}(undef, length(A)), Vector{Vector{Ti}}(undef, length(A)), Vector{Vector{Tv}}(undef, length(A))
  for (i, a) in enumerate(A)
    (I, J, V) = findnz(a)
    Is[i] =  I
    Js[i] =  J .+ coloffset
    Vs[i] =  V
    coloffset += a.n
  end
  sparse(reduce(vcat, Is), reduce(vcat, Js), reduce(vcat, Vs),  A[1].m, coloffset)
end

# using CuArrays
# using CUSPARSE
# CuArrays.allowscalar(false)
# back(::typeof(*), Δ, a::CuMatrix, b::CudaSparseMatrixCSC) = Flux.Tracker.@back(a, transpose(A_mul_Bt(b,Δ)))
# back(::typeof(mul), Δ, a::CuMatrix, b::CudaSparseMatrixCSC) = Flux.Tracker.@back(a, transpose(A_mul_Bt(b,Δ)))
# mul(A::CuMatrix{T},B::CudaSparseMatrixCSC{T}) where T = transpose(At_mul_Bt(A,B))
# a::Flux.Tracker.TrackedMatrix * b::CudaSparseMatrixCSC = Flux.Tracker.TrackedArray(Flux.Tracker.Call(mul, a, b))

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
# simple program to test the new k-means (not ready yet)

using Base.Test
using Clustering

srand(34568)

m = 3
n = 1000
k = 10

x = rand(m, n)

# non-weighted
r = kmeans(x, k; maxiter=50)
@test isa(r, KmeansResult{Float64})
@test size(r.centers) == (m, k)
@test length(r.assignments) == n
@test all(r.assignments .>= 1) && all(r.assignments .<= k)
@test length(r.costs) == n
@test length(r.counts) == k
@test sum(r.counts) == n
@test r.cweights == @compat(map(Float64, r.counts))
@test_approx_eq sum(r.costs) r.totalcost

# non-weighted (float32)
r = kmeans(@compat(map(Float32, x)), k; maxiter=50)
@test isa(r, KmeansResult{Float32})
@test size(r.centers) == (m, k)
@test length(r.assignments) == n
@test all(r.assignments .>= 1) && all(r.assignments .<= k)
@test length(r.costs) == n
@test length(r.counts) == k
@test sum(r.counts) == n
@test r.cweights == @compat(map(Float64, r.counts))
@test_approx_eq sum(r.costs) r.totalcost

# weighted
w = rand(n)
r = kmeans(x, k; maxiter=50, weights=w)
@test isa(r, KmeansResult{Float64})
@test size(r.centers) == (m, k)
@test length(r.assignments) == n
@test all(r.assignments .>= 1) && all(r.assignments .<= k)
@test length(r.costs) == n
@test length(r.counts) == k
@test sum(r.counts) == n

cw = zeros(k)
for i = 1:n
    cw[r.assignments[i]] += w[i]
end
@test_approx_eq r.cweights cw

@test_approx_eq dot(r.costs, w) r.totalcost


text:
@testset "Internal Methods" begin
    @testset "_support_sum_coeffs" begin
        @test IOMT._support_sum_coeffs([1, 2]) == [1, 1]
    end
end

@testset "Support Sum" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in [0, 1])
    @infinite_parameter(m, y[1:2] in [0, 1])
    @variable(m, inf, Infinite(x, y))

    @test InfiniteOpt._index_type(support_sum(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(support_sum(inf, y)) == MeasureIndex
    @test name(support_sum(inf, x)) == "support_sum"
    @test support_sum(inf, y, label = UserDefined) isa GeneralVariableRef
end

@testset "Macro" begin
    m = InfiniteModel()
    @infinite_parameter(m, x in [0, 1])
    @infinite_parameter(m, y[1:2] in [0, 1])
    @variable(m, inf, Infinite(x, y))

    @test InfiniteOpt._index_type(@support_sum(inf, x)) == MeasureIndex
    @test InfiniteOpt._index_type(@support_sum(inf, y)) == MeasureIndex
    @test_macro_throws ErrorException @support_sum(inf, y, 2)
end


text:
export Diagm

"C[:,:,k]=diagm(A[:,k])"
function Diagm(A::Array{Float64,2})
  C=zeros(size(A,1),size(A,1),size(A,2))
  @simd for k=1:1:size(A,2)
    @inbounds C[:,:,k]=diagm(vec(A[:,k]))
  end
  C
end


text:
function nll(m::NystromModel)
	if !same_hyp(m)
		recompute_kernel_eig!(m)
	end
	refresh_cache!(m)
	sigma2 = exp(2m.lik.log_sigma)
	neig = length(m.nystrom.values)
	LambdaInv = diagm(neig./m.nystrom.values)
	Z = (m.PhiTPhi + sigma2 * LambdaInv)
	yTQy = (m.yTy - dot(m.PhiTy, Z \ m.PhiTy)) / sigma2
	logdetZ = logdet(Z)
	n = size(m.X, 1)
	scale = n - length(m.nystrom.values)
	logdetQ = scale*log(sigma2) + logdetZ + sum(log(m.nystrom.values))
	yTQy / 2 + logdetQ / 2 + scale*log(2pi)
end

function nll(m::FullModel)
	refresh_cache!(m)
	n = size(m.X, 1)
	dot(m.y, m.CinvY) / 2 + logdet(m.C) / 2 + n*log(2pi)
end

function nll(m::FITCModel)
	sigma2 = exp(2m.lik.log_sigma)
	Kmn = kernel(m.k, m.u, m.X)
	Km = kernel(m.k, m.u)
	Knn = kernel(m.k, m.X, m.X)
	Lambda = diagm(diag(Knn - Kmn' * (Km \ Kmn)))
	C = Kmn' * (Km \ Kmn) + Lambda + sigma2 * eye(Lambda)
	n = size(m.y)[1]
	(dot(m.y,(C \ m.y))  .+ logdet(C) .+ n*log(2pi))/2
end

# function nll(m::LaplacianModel)
# 	refresh_cache!(m)
# 	n = size(m.X, 1)
# 	sigma2 = exp(2m.lik.log_sigma)
# 	scale = n - length(m.lambdas)
# 	Ssqrtlambda = Float64[S(m.k, sqrt(lambda)) for lambda in m.lambdas]
# 	LambdaInv = diagm(1./Ssqrtlambda)
# 	Z = (m.PhiTPhi + sigma2 * LambdaInv)
# 	logdetZ = logdet(Z)
# 	logdetQ = scale*log(sigma2) + logdetZ + sum(log(Ssqrtlambda))
# 	yTQy = (m.yTy - dot(m.PhiTy, Z \ m.PhiTy)) / sigma2
# 	yTQy / 2 + logdetQ / 2 + scale*log(2pi)
# end

include("negloglik_laplacian.jl")


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 1, Zadanie 6.
#

fl = Float64

function f(x::fl)
    return fl(sqrt(fl(fl(x^2) + 1)) - 1)
end

function g(x::fl)
    return fl(fl(x^2) / fl(sqrt(fl(fl(x^2) + 1)) + 1))
end

# print out the results
for i in (-1:-1:-16)
    o1 = f(fl(fl(8)^i))
    o2 = g(fl(fl(8)^i))
    println(
        "f(8^{", i, "}) &= ", o1, "\\\\ \n",
        "g(8^{", i, "}) &= ", o2, "\\\\",
        # "equal? ", (o1 == o2), "\n"
    )
end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
using Gadfly

for i in 0:5
if i < 4 println( i , " is less than 4") end
if i >= 4 println(i, " is greater or equal to 4") end
end

i = 0

while i < 5
println(i)
i=i+1
end

output = open("fake_data_trial.txt", "w")
dataout = rand(2,100)

writedlm(output, dataout)
close(output)

datain = readdlm("fake_data_trial.txt")

plot(x = datain[1, :], y = datain[2, :])

plot(x = datain[1, :], y = datain[2, :], Geom.point, Geom.line)

myfirstplot = plot(x = datain[1, :], y = datain[2, :],
Geom.point, Geom.line)

draw(PDF("my_first_plot.pdf", 9inch, 6inch), myfirstplot)

max_stars_repo_name: thehackerwithin/swinburne

max_issues_repo_name: thehackerwithin/swinburne

max_forks_repo_name: thehackerwithin/swinburne


text:
pprint(io, x) = print(io, x)
function pprint(io, x::Expr)
    if x.head == :block
        args = [a for a in x.args if !(a isa LineNumberNode)]
        if length(args) == 1
            pprint(io, args[1])
        else
            for (i, a) in enumerate(args)
                i == 1 || println(io)
                pprint(io, a)
            end
        end
    elseif x.head == :(=)
        @assert length(x.args) == 2
        pprint(io, x.args[1])
        print(io, " = ")
        pprint(io, x.args[2])
    else
        print(io, x)
    end
end

prettyexpr(x) = sprint(pprint, x)


text:
#fill in structs when geoarray/geostats ecosystem stabilize?
struct AbstractLandsatImage end

#TODO: add older landsat types
#struct Landsat5Image <: AbstractLandsatImage end

struct Landsat7ETMImage <: AbstractLandsatImage
    b1
    b2
    b3
    b4
    b5
    b6
    b7
    b8
    row
    path
    date
end

struct Landsat8OLIImage <: AbstractLandsatImage
    b1
    b2
    b3
    b4
    b5
    b6
    b7
    b8
    b9
    b10
    b11
    row
    path
    date
end

struct AbstractSentinelImage end

#other metadata?
struct Sentinel2AImage <: AbstractSentinelImage
    b1
    b2
    b3
    b4
    b5
    b6
    b7
    b8
    b8a
    b9
    b10
    b11
    b12
end

struct Sentinel2BImage <: AbstractSentinelImage
    b1
    b2
    b3
    b4
    b5
    b6
    b7
    b8
    b8a
    b9
    b10
    b11
    b12
end


text:
module SimplePaint

import MiniFB as MFB
import REPL
import SimpleDraw as SD

end


text:
module Utils

export Protocol, Endpoint

@enum Protocol begin
    rest_xml = 1
    json = 2
end

struct Endpoint
    service::String
    protocol::Protocol
    global_endpoint::String
    api_version::String
end

end  # End Module

max_stars_repo_name: mattBrzezinski/AWSCore-Prototype

max_issues_repo_name: mattBrzezinski/AWSCore-Prototype

max_forks_repo_name: mattBrzezinski/AWSCore-Prototype


text:
## 02.jl
# 交互に連結

s2= "タクシー"
s1= "パトカー"
s = join([string(x,y) for (x,y) in zip(s1,s2)],"")
println(s)


text:
using Jecco, Jecco.AdS5_3_1
import Jecco.AdS5_3_1: geta4, getfx2, getfy2, getxi, getB1, getB2, getG, getphi

grid = SpecCartGrid3D(
    x_min            = -5.0,
    x_max            =  5.0,
    x_nodes          =  32,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  32,
    u_outer_min      =  0.1,
    u_outer_max      =  1.005,
    u_outer_domains  =  3,
    u_outer_nodes    =  32,
    u_inner_nodes    =  12,
    fd_order         =  6,
    sigma_diss       =  0.2,
)

potential = Phi8Potential(
    oophiM2 = -1.0,
    oophiQ  = 0.1,
)

evoleq = AffineNull(
    phi0           = 1.0,
    potential      = potential,
    gaugecondition = ConstantAH(u_AH = 1.0),
)

dirname = "/home/mzilhao/dev/julia/Jecco/examples/spinodal2D_model1_en0.9_run01/"

it = 0

ts_bulk  = OpenPMDTimeSeries(dirname, prefix="bulk_")
ts_bdry  = OpenPMDTimeSeries(dirname, prefix="boundary_")
ts_gauge = OpenPMDTimeSeries(dirname, prefix="gauge_")

# atlas of grid configuration and respective SystemPartition
atlas     = Atlas(grid)
systems   = SystemPartition(grid)
Nsys      = length(systems)

# allocate variables
boundary       = Boundary(grid)
boundary_t     = Boundary(grid)
gauge          = Gauge(grid)
gauge_t        = Gauge(grid)
bulkevols      = BulkEvolvedPartition(grid)
bulkevols_t    = BulkEvolvedPartition(grid)
bulkconstrains = BulkConstrainedPartition(grid)
bulkderivs     = BulkDerivPartition(grid)
cache          = AdS5_3_1.HorizonCache(systems[end], evoleq.gaugecondition.fd_order)

charts = atlas.charts

# for the boundary/xi grid
empty   = Cartesian{1}("u", systems[1].ucoord[1], systems[1].ucoord[1], 1)
chart2D = Chart(empty, systems[1].xcoord, systems[1].ycoord)

# restore data
AdS5_3_1.restore!(bulkevols, ts_bulk, it)
AdS5_3_1.restore!(boundary, ts_bdry, it)
AdS5_3_1.restore!(gauge, ts_gauge, it)

t0 = ts_bulk.current_t

tinfo  = Jecco.TimeInfo(it, t0, 0.0, 0.0)

# nested system function
nested = AdS5_3_1.Nested(systems, bulkconstrains, bulkderivs)

AdS5_3_1.compute_boundary_t!(boundary_t, bulkevols[1], boundary, gauge, systems[1], evoleq)

# solve nested system for the constrained variables
nested(bulkevols, boundary, gauge, evoleq)

AdS5_3_1.compute_xi_t!(gauge_t, bulkconstrains[Nsys], bulkevols[Nsys],
                       bulkderivs[Nsys], gauge, cache, systems[Nsys],
                       evoleq.gaugecondition)

for aa in 1:Nsys
    sys           = systems[aa]
    bulkevol_t    = bulkevols_t[aa]
    bulkevol      = bulkevols[aa]
    bulkconstrain = bulkconstrains[aa]

    AdS5_3_1.compute_bulkevolved_t!(bulkevol_t, bulkconstrain, gauge_t, bulkevol,
                                    boundary, gauge, sys, evoleq)
end
AdS5_3_1.sync_bulkevolved!(bulkevols_t, bulkconstrains, gauge_t, systems, evoleq)

# output structures
out_bdry_t  = Jecco.Output(dirname, "boundary_t_", tinfo)
out_gauge_t = Jecco.Output(dirname, "gauge_t_", tinfo)
out_bulk_t  = Jecco.Output(dirname, "bulk_t_", tinfo)

# output fields
boundary_fields = (
    Jecco.Field("a4_t",   geta4(boundary_t),  chart2D),
    Jecco.Field("fx2_t",  getfx2(boundary_t), chart2D),
    Jecco.Field("fy2_t",  getfy2(boundary_t), chart2D),
)

gauge_fields = Jecco.Field("xi_t", getxi(gauge_t), chart2D)

bulkevols_fields = ntuple(i -> (
    Jecco.Field("B1_t c=$i",  getB1(bulkevols_t[i]),  charts[i]),
    Jecco.Field("B2_t c=$i",  getB2(bulkevols_t[i]),  charts[i]),
    Jecco.Field("G_t c=$i",   getG(bulkevols_t[i]),   charts[i]),
    Jecco.Field("phi_t c=$i", getphi(bulkevols_t[i]), charts[i])
), Nsys)

out_bdry_t(boundary_fields)
out_gauge_t(gauge_fields)
out_bulk_t.(bulkevols_fields)


text:
##### model_diagnostics.jl
# Created by Alexander Dowling

###############
##### Exported functions for model diagnostics

function round_to_bounds!(m::Model, epsilon::Float64=1.0E-6, f=STDOUT)
	# This function rounds each element of the solution vector to the closest bound if
	# it is within epsilon.
	# Why is this necessary? Some solvers (e.g., IPOPT) give solutions that are slightly
	# infeasible or violate bounds. This is problematic in a meta-algorithm with calls
	# to other solvers.

	x = get_x(m, :Unsolved)

	up = m.colUpper
	lo = m.colLower

	for i = 1:length(x)

		v = Variable(m, i)

		if(up[i] < lo[i])
			println(f,"Warning: Upper bound (",up[i],") is smaller than lower bound (",lo[i],") for variable ",v)
		end

		x_ = NaN

		if(x[i] > up[i])
			x_ = up[i]
		elseif(x[i] < lo[i])
			x_ = lo[i]
		elseif(up[i] - x[i] < epsilon)
			x_ = up[i]
		elseif(x[i] - lo[i] < epsilon)
			x_ = lo[i]
		end

		if(!isnan(x_))
			setvalue(v, x_)
		end

	end

	return nothing

end

function printVariableDiagnostics(m::Model, epsilon::Float64=1.0E-6, f=STDOUT, status::Symbol=:Unknown)
	x = get_x(m, status)

	println(f,"Uninitialized Variables: ")
	for i = 1:length(x)
		if(isnan(x[i]))
			v = Variable(m, i)
			#println("x[",i,"] = ",v)
			println(f,v)
		end
	end
	println(f," ")

	println(f,"Variable Violating Bounds: ")
	up = m.colUpper
	lo = m.colLower

	for i = 1:length(x)
		viol_up = x[i] - up[i] > epsilon
		viol_lo = x[i] - lo[i] < -epsilon
		if(viol_up || viol_lo)

			v = Variable(m, i)
			print(f,v, " = ",x[i], ", bounds: [ ")

			if(viol_lo)
				print_with_color(:red, f, string(lo[i]))
			else
				print(f, string(lo[i]))
			end

			print(f, " , ")

			if(viol_up)
				print_with_color(:green, f, string(up[i]))
			else
				print(f, string(up[i]))
			end

			println(f, " ]")
		end
	end

	return nothing

end

function printInfeasibleEquations(m2::Model, eqns, f=STDOUT, status::Symbol=:Unknown)

	dd = DegenData(m2, status)
	return printInfeasibleEquations(m2, dd, eqns, f)

end

function printInfeasibleEquations(m2::Model, dd::DegenData, eqns, f=STDOUT)

	println(f,"Checking infeasibility for specified equations: ")

	r = reshape(min(dd.g - dd.gLB, 0) + max(dd.g - dd.gUB, 0),length(dd.g))

	for i = eqns
		print(f, "r[",string(i),"] = ")
		print_with_color(:green,f,string(r[i]))
		print(f,"\n")

		printEquation(m2, dd, i, f)

		println(f, " ")

	end

	return nothing

end

function printInfeasibleEquations(m2::Model, epsilon::Float64, f=STDOUT, status::Symbol=:Unknown)

	dd = DegenData(m2, f, status)
	return printInfeasibleEquations(m2, dd, epsilon, f)
end

function printInfeasibleEquations(m2::Model, dd::DegenData, epsilon::Float64, f=STDOUT)

	println(f,"Infeasible equations: ")

	#=
	println("size(dd.g) = ",size(dd.g))
	println("size(dd.gLB) = ",size(dd.gLB))
	println("size(dd.gUB) = ",size(dd.gUB))
	=#

	r = reshape(min(dd.g - dd.gLB, 0) + max(dd.g - dd.gUB, 0),length(dd.g))

	k = sortperm(abs(r))

	for j = 1:length(k)

		i = k[j]

		if(abs(r[i]) >= epsilon || isnan(r[i]))
			print(f,"r[",string(i),"] = ")
			print_with_color(:green,f,string(r[i]))
			print(f,"\n")

			printEquation(m2, dd, i,f)

			printVariablesInEquation(m2, dd, i, true, true, f)

			println(f," ")

		end

	end

	return nothing

end

function printInactiveEquations(m2::Model, epsilon::Float64=1E-6, f=STDOUT)

	dd = DegenData(m2, f)
	return printInactiveEquations(m2, dd, epsilon, f)

end

function printInactiveEquations(m2::Model, dd::DegenData, epsilon::Float64=1E-6, f=STDOUT)

	println(f,"Inactive Equations: ")

	inactive = ((dd.gUB - dd.g) .> epsilon) & ((dd.g - dd.gLB) .> epsilon)

	for i = 1:length(dd.g)
		if(((dd.gUB[i] - dd.g[i]) .> epsilon) & ((dd.g[i] - dd.gLB[i]) .> epsilon))
			printEquation(m2, dd, i, f)
		end
	end

	return nothing

end

function checkVarBounds(m::Model, f=STDOUT)

	low = m.colLower
	up = m.colUpper
	val = m.colVal

	n = length(val)

	for i = 1:n
		if(up[i] < low[i] || val[i] < low[i] || val[i] > up[i])

			print(f,getname(m,i),":\t\t","Lower: ")

			if(up[i] < low[i] || val[i] < low[i])
				print_with_color(:red,f,string(low[i]))
			else
				print(f,low[i])
			end

			print(f,"\t Upper: ")

			if(up[i] < low[i] || val[i] > up[i])
				print_with_color(:blue,f,string(up[i]))
			else
				print(f,up[i])
			end

			print(f,"\t Value: ")

			if(val[i] < low[i] || val[i] > up[i])
				print_with_color(:green,f,string(val[i]))
			else
				print(f,val[i])
			end

			print(f,"\n")
		end
	end

	return nothing

end

function printBound(m::Model, i::Int64, epsiActive::Float64, f=STDOUT)

	v = Variable(m, i)

	up = m.colUpper[i]
	lo = m.colLower[i]
	val = m.colVal[i]

	if(abs(val-lo) < epsiActive)
		print(f,string(lo)," <= ")
	end

	print(getname(v))

	if(abs(val-up) < epsiActive)
		print(f," <= ",string(up))
	end

	print(f,"\n")

	return nothing

end

function printRows(m::Model, rows::Int64, f=STDOUT; lite::Bool=false)
	return printRows(m, Array([rows]), f, lite)
end

function printRows(m::Model, rows::Array{Int64,1}, f=STDOUT; lite::Bool=false)

	dd = DegenData(m, f)
	return printRows(m, dd, rows, f, lite)

end

function printRows(m::Model, dd::DegenData, rows::Int64, f=STDOUT; lite::Bool=false)
	return printRows(m, dd, Array([rows]), f, lite=lite)
end

function printRows(m::Model, dd::DegenData, rows::Array{Int64,1}, f=STDOUT; lite::Bool=false)

	if(lite && length(rows) == 1)

		r = rows[1]

		if(r <= dd.nConstr)
			printEquation(m, dd, rows[1], f)
		else
			printBound(m, r - dd.nConstr, 1E-6, f)
		end
	else

		for i = 1:length(rows)
			r = rows[i]

			println(f,"**********************************")

			if(r <= dd.nConstr)

				println(f,"Constraint ",string(r),"...")
				printEquation(m, dd, r, f)
				println(f,string(dd.gLB[r])," <= ", string(dd.g[r]), " <= ",string(dd.gUB[r]))

				println(f," ")
				println(f,"Involved Variables:")
				k = find(dd.iR .== r)

				for j in unique(dd.jC[k])
					v = Variable(m,j)
					print(f,getname(v)," \t")
					println(f,getlowerbound(v)," <= ",getvalue(v)," <= ",getupperbound(v))
				end

			else
				v = r - dd.nConstr
				println(f,"Bound for variable ",v)
				printBounds(m, v, 1E-6, f)

			end

			println(f," ")

		end
	end

	return nothing

end

function checkEquationScaling(m::Model, status::Symbol=:Unknown, f=STDOUT)

	dd = DegenData(m, status)

	checkJacobianScaling(m, dd, f)

	return nothing

end


text:
using Test
using Forecast

@testset "summarize" begin

    x = summarize(rand(100))
    @test size(x.quantiles) == (1,8)
    @test size(x.moments) == (1,5)
    @test size(x.format) == (1,2)

    @test x.quantiles isa DataFrame
    @test x.moments isa DataFrame
    @test x.format isa DataFrame

    x = summarize(rand(100,2))
    @test size(x.quantiles) == (2,8)
    @test size(x.moments) == (2,5)
    @test size(x.format) == (2,2)

end


text:
"""
    ExactReach

ExactReach performs exact reachability analysis to compute the output reachable set for a network.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: HPolytope
3. Output: HPolytope

# Return
`ReachabilityResult`

# Method
Exact reachability analysis.

# Property
Sound and complete.

# Reference
[W. Xiang, H.-D. Tran, and T. T. Johnson,
"Reachable Set Computation and Safety Verification for Neural Networks with ReLU Activations,"
*ArXiv Preprint ArXiv:1712.08163*, 2017.](https://arxiv.org/abs/1712.08163)
"""
struct ExactReach end

function solve(solver::ExactReach, problem::Problem)
    reach = forward_network(solver, problem.network, problem.input)
    return check_inclusion(reach, problem.output)
end

function forward_layer(solver::ExactReach, layer::Layer, input::Vector{HPolytope})
    output = Vector{HPolytope}(undef, 0)
    for i in 1:length(input)
        input[i] = linear_transformation(layer, input[i])
        append!(output, forward_partition(layer.activation, input[i]))
    end
    return output
end

function forward_layer(solver::ExactReach, layer::Layer, input::HPolytope)
    input = linear_transformation(layer, input)
    return forward_partition(layer.activation, input)
end

function forward_partition(act::ReLU, input::HPolytope)
    n = dim(input)
    output = Vector{HPolytope}(undef, 0)
    C, d = tosimplehrep(input)
    dh = [d; zeros(n)]
    for h in 0:2^n-1
        P = getP(h, n)
        Ch = [C; I - 2P]
        input_h = HPolytope(Ch, dh)
        if !isempty(input_h)
            push!(output, linear_transformation(Matrix(P), input_h))
        end
    end
    return output
end

function getP(h::Int64, n::Int64)
    str = string(h-1, pad = n, base = 2)
    vec = Vector{Int64}(undef, n)
    for i in 1:n
        vec[i] = ifelse(str[i] == '1', 1, 0)
    end
    return Diagonal(vec)
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/25/2014

if !isdefined(:PACKAGE_PATH)
    PACKAGE_PATH = Pkg.dir("SISLES", "src")
end

push!(LOAD_PATH, "$PACKAGE_PATH/include")
push!(LOAD_PATH, "$PACKAGE_PATH/common")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter/CorrAEMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter/LLAEMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/PilotResponse")
push!(LOAD_PATH, "$PACKAGE_PATH/PilotResponse/SimplePilotResponseImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/DynamicModel")
push!(LOAD_PATH, "$PACKAGE_PATH/DynamicModel/SimpleADMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/WorldModel")
push!(LOAD_PATH, "$PACKAGE_PATH/WorldModel/AirspaceImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Sensor")
push!(LOAD_PATH, "$PACKAGE_PATH/Sensor/SimpleTCASSensorImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/CollisionAvoidanceSystem")
push!(LOAD_PATH, "$PACKAGE_PATH/CollisionAvoidanceSystem/SimpleTCASImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Simulator")
push!(LOAD_PATH, "$PACKAGE_PATH/Simulator/TCASSimulatorImpl")

using Util
using Encounter
using PilotResponse
using DynamicModel
using WorldModel
using Sensor
using CollisionAvoidanceSystem
using Simulator

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using CSV
using PyPlot
using HTTP
using StableRNGs

req = HTTP.get("https://raw.githubusercontent.com/rupakc/UCI-Data-Analysis/master/Airfoil%20Dataset/airfoil_self_noise.dat");

df = CSV.read(req.body, DataFrames.DataFrame; header=[
                   "Frequency","Attack_Angle","Chord+Length",
                   "Free_Velocity","Suction_Side","Scaled_Sound"
                   ]
              );
df[1:5, :] |> pretty

schema(df)

y, X = unpack(df, ==(:Scaled_Sound));

X = MLJ.transform(fit!(machine(Standardizer(), X)), X);

train, test = partition(collect(eachindex(y)), 0.7, shuffle=true, rng=StableRNG(612));

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

coerce!(X, :Frequency=>Continuous)

for model in models(matching(X, y))
       print("Model Name: " , model.name , " , Package: " , model.package_name , "\n")
end

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree

dcrm = machine(DecisionTreeRegressor(), X, y)

fit!(dcrm, rows=train)
pred_dcrm = predict(dcrm, rows=test);

rms(pred_dcrm, y[test])

RandomForestRegressor = @load RandomForestRegressor pkg=DecisionTree
rfr = RandomForestRegressor()

rfr_m = machine(rfr, X, y);

fit!(rfr_m, rows=train);

pred_rfr = predict(rfr_m, rows=test);
rms(pred_rfr, y[test])

r_maxD = range(rfr, :n_trees, lower=9, upper=15)
r_samF = range(rfr, :sampling_fraction, lower=0.6, upper=0.8)
r = [r_maxD, r_samF];

tuning = Grid(resolution=7)
resampling = CV(nfolds=6)

tm = TunedModel(model=rfr, tuning=tuning,
                resampling=resampling, ranges=r, measure=rms)

rfr_tm = machine(tm, X, y);

fit!(rfr_tm, rows=train);

pred_rfr_tm = predict(rfr_tm, rows=test);
rms(pred_rfr_tm, y[test])

fitted_params(rfr_tm).best_model

r = report(rfr_tm)
res = r.plotting

md = res.parameter_values[:,1]
mcw = res.parameter_values[:,2]

figure(figsize=(8,6))
tricontourf(md, mcw, res.measurements)

xlabel("Number of trees", fontsize=14)
ylabel("Sampling fraction", fontsize=14)
xticks(9:1:15, fontsize=12)
yticks(fontsize=12)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using LocalREPLHistory
using Test

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


text:
using TopOpt.TopOptProblems.IO.INP
using JuAFEM

using Test

cube = INP.Parser.import_inp(joinpath(@__DIR__, "testcube.inp"))
dh = cube.dh
grid = dh.grid
cells = grid.cells
nodes = grid.nodes
getdim(::JuAFEM.Cell{N}) where {N} = N

@test dh.field_dims == [3]
@test getdim(cells[1]) == 3 # 3D cells
@test JuAFEM.nfaces(cells[1]) == 4 # Tetrahedron
@test JuAFEM.nnodes(cells[1]) == 10 # Quadratic tetrahedron
@test length(grid.boundary_matrix.nzval) == 16

raw_inp = INP.Parser.extract_inp(joinpath(@__DIR__, "testcube.inp"))
@test raw_inp.celltype == "C3D10"
@test raw_inp.E == 70_000
@test raw_inp.ν == 0.3

@test raw_inp.cellsets["Eall"] == 1:5
@test raw_inp.cellsets["Evolumes"] == 1:5
@test raw_inp.cellsets["SolidMaterialSolid"] == 1:5
force_node = collect(keys(raw_inp.cloads))[1]
@test raw_inp.node_coords[force_node] == (10, 10, 10)
@test raw_inp.cloads[force_node] == [0, -1, 0]

@test raw_inp.facesets["DLOAD_SET_1"] == [(1,3), (5,2)]
@test raw_inp.dloads["DLOAD_SET_1"] == 1

@test raw_inp.nodedbcs["FemConstraintDisplacement"] == [(1, 0), (2, 0), (3, 0)]
for n in raw_inp.nodesets["FemConstraintDisplacement"]
    @test raw_inp.node_coords[n][3] == 0
end


text:
using OpenCL

function compile(input, output)
	device, ctx, queue = cl.create_compute_context()
	kernel = read(input, String)
	p = cl.Program(ctx, source=kernel) |> cl.build!
	bin = cl.info(p, Symbol("binaries"))

	write(output, bin[device])
end

compile("./kernel/kernel.cl", "./bin/kernel.bin")


text:
@safetestset "isiso31661alpha2" begin
    using Validator

    @test Validator.isiso31661alpha2("FR") === true
    @test Validator.isiso31661alpha2("fR") === true
    @test Validator.isiso31661alpha2("GB") === true
    @test Validator.isiso31661alpha2("PT") === true
    @test Validator.isiso31661alpha2("CM") === true
    @test Validator.isiso31661alpha2("JP") === true
    @test Validator.isiso31661alpha2("PM") === true
    @test Validator.isiso31661alpha2("ZW") === true
    @test Validator.isiso31661alpha2("MM") === true
    @test Validator.isiso31661alpha2("cc") === true
    @test Validator.isiso31661alpha2("GG") === true

    @test Validator.isiso31661alpha2("") === false
    @test Validator.isiso31661alpha2("FRA") === false
    @test Validator.isiso31661alpha2("AA") === false
    @test Validator.isiso31661alpha2("PI") === false
    @test Validator.isiso31661alpha2("RP") === false
    @test Validator.isiso31661alpha2("WV") === false
    @test Validator.isiso31661alpha2("WL") === false
    @test Validator.isiso31661alpha2("UK") === false
    @test Validator.isiso31661alpha2("ZZ") === false
end


text:
function _compute_lnn_to_coverage()
    lock(Main.MicroCoverage_lock) do
        lnn_to_tracker_indices = Dict{LineNumberNode, Vector{Int}}()
        for tracker_idx = 1:length(Main.MicroCoverage_tracker_linenumbernodes)
            lnn = Main.MicroCoverage_tracker_linenumbernodes[tracker_idx]
            if !haskey(lnn_to_tracker_indices, lnn)
                lnn_to_tracker_indices[lnn] = Vector{Int}()
            end
            push!(lnn_to_tracker_indices[lnn], tracker_idx)
        end
        lnn_to_coverage = Dict{LineNumberNode, Vector{Bool}}()
        for k in keys(lnn_to_tracker_indices)
            tracker_indices = lnn_to_tracker_indices[k]
            lnn_to_coverage[k] = Main.MicroCoverage_tracker[tracker_indices]
        end
        return lnn_to_coverage
    end
end

function _compute_coverage_file_contents(filename::Symbol;
                                         min_padding::Integer,
                                         max_padding::Integer,
                                         source_code_filename::Symbol = filename)::Vector{String}
    _source_code_filename = string(source_code_filename)
    original_source_string::String = read(_source_code_filename, String)::String
    original_source_lines = split(original_source_string, '\n')
    num_lines = length(original_source_lines)
    lnn_to_coverage = _compute_lnn_to_coverage()
    coverage_lines_prefixes = Vector{String}(undef, num_lines)
    for i = 1:num_lines
        coverage_lines_prefixes[i] = _coverage_boolvector_to_string(get(lnn_to_coverage,
                                                                        LineNumberNode(i,
                                                                        filename), Bool[]))
    end
    coverage_lines = Vector{String}(undef, num_lines)
    for i = 1:num_lines
        this_line_coverage_prefix = coverage_lines_prefixes[i]
        this_line_padding_length = max(min_padding,
                                       max_padding - length(this_line_coverage_prefix))
        this_line_padding = repeat(" ",
                                   max(1, this_line_padding_length))
        coverage_lines[i] = string(this_line_coverage_prefix,
                                   this_line_padding,
                                   original_source_lines[i])
    end
    return coverage_lines
end

function _coverage_boolvector_to_string(line_coverage::Vector{Bool})
    if length(line_coverage) < 1
        return "-"
    end
    return string("[", join(string.(convert(Vector{Int}, line_coverage)), ","), "]")
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 2dad73b0-b0ea-11eb-108d-419e444c4b1d
md"""
# [Day 23: Opening the Turing Lock](https://adventofcode.com/2015/day/23)

Little Jane Marie just got her very first computer for Christmas from some unknown benefactor. It comes with instructions and an example program, but the computer itself seems to be malfunctioning. She's curious what the program does, and would like you to help her run it.

The manual explains that the computer supports two [registers](https://en.wikipedia.org/wiki/Processor_register) and six [instructions](https://en.wikipedia.org/wiki/Instruction_set) (truly, it goes on to remind the reader, a state-of-the-art technology). The registers are named `a` and `b`, can hold any [non-negative integer](https://en.wikipedia.org/wiki/Natural_number), and begin with a value of `0`. The instructions are as follows:

-   `hlf r` sets register `r` to half its current value, then continues with the next instruction.
-   `tpl r` sets register `r` to triple its current value, then continues with the next instruction.
-   `inc r` increments register `r`, adding `1` to it, then continues with the next instruction.
-   `jmp offset` is a jump; it continues with the instruction `offset` away relative to itself.
-   `jie r, offset` is like `jmp`, but only jumps if register `r` is even ("jump if even").
-   `jio r, offset` is like `jmp`, but only jumps if register `r` is `1` ("jump if one", not odd).

All three jump instructions work with an offset relative to that instruction. The offset is always written with a prefix `+` or `-` to indicate the direction of the jump (forward or backward, respectively). For example, `jmp +1` would simply continue with the next instruction, while `jmp +0` would continuously jump back to itself forever.

The program exits when it tries to run an instruction beyond the ones defined.

For example, this program sets `a` to `2`, because the `jio` instruction causes it to skip the `tpl` instruction:

```
inc a
jio a, +2
tpl a
inc a

```

**What is the value in register `b` when the program in your puzzle input is finished executing?**
"""

# ╔═╡ a7b969c5-1aa9-46a7-bc21-7dbdae27128c
puzzle_input = "jio a, +18
inc a
tpl a
inc a
tpl a
tpl a
tpl a
inc a
tpl a
inc a
tpl a
inc a
inc a
tpl a
tpl a
tpl a
inc a
jmp +22
tpl a
inc a
tpl a
inc a
inc a
tpl a
inc a
tpl a
inc a
inc a
tpl a
tpl a
inc a
inc a
tpl a
inc a
inc a
tpl a
inc a
inc a
tpl a
jio a, +8
inc b
jie a, +4
tpl a
inc a
jmp +2
hlf a
jmp -7"

# ╔═╡ Cell order:
# ╟─2dad73b0-b0ea-11eb-108d-419e444c4b1d
# ╟─a7b969c5-1aa9-46a7-bc21-7dbdae27128c


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "bsrsv2" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "bsrsv2!" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_X = CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
            h_Y = to_host(d_X)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_X = CudaArray(rand(elty,n))
            @test_throws DimensionMismatch CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.sv2!('N','U',alpha,d_A,d_X,'O')
        end

        @testset "bsrsv2" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_Y = CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_Y = UpperTriangular(d_A)\d_X
            h_Y = to_host(d_Y)
            @test h_Y ≈ A\X
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.sv2('N','U',alpha,d_A,d_X,'O')
        end
    end
end


text:
using HDF5

struct LoadInstace
    parcels::Any
    crowd_shippers::Any
    distance::Any
    wait::Any
end

function load_instance(file_name, entry_name)
    h5open(file_name) do file
        load(file[entry_name])
    end
end

function make_adjacent_matrix(li::LoadInstace)
    sz = size(li.distance, 1)
    adj = ones(sz, sz)
    for i = 1:sz
        adj[i, i] = 0
    end
    adj
end

function load(data)
    parcels = read(data["parcels"])
    csrs = read(data["csrs"])
    distance = read(data["delay"])
    wait = read(data["wait"])
    li = LoadInstace(parcels, csrs, distance, wait)
    T = make_sparse_csr_array(li)
    A = make_arrive_matrix(li)
    li, T, A
end

function make_sparse_csr_array(li::LoadInstace)
    stations = size(li.distance, 1)
    csr = size(li.crowd_shippers, 1)
    T = [spzeros(Int8, stations, stations) for i = 1:csr]
    for i = 1:csr
        s, d = li.crowd_shippers[i, 2:3]
        T[i][s, d] = 1
    end
    T
end

function make_arrive_matrix(li::LoadInstace)
    c = size(li.crowd_shippers, 1)
    s = size(li.distance, 1)
    A = Inf * ones(s, c)
    for (i, (w, s, d)) in enumerate(eachrow(li.crowd_shippers))
        A[s, i] = w
        A[d, i] = w + li.wait[s] + li.distance[s, d]
    end
    A
end

max_stars_repo_name: FilippoRanza/crowd-shipping-models

max_issues_repo_name: FilippoRanza/crowd-shipping-models

max_forks_repo_name: FilippoRanza/crowd-shipping-models


text:
function getdescriptors(descriptors)

    for i = length(descriptors):-1:1
        if descriptors[i] === nothing
            deleteat!(descriptors, i)
        end
    end

    n = length(descriptors)
    emdescriptors = Array{Any}(nothing, n)
    mldescriptors = nothing;

    n = 0
    for i = 1:length(descriptors)
        if typeof(descriptors[i]) == Potential.PotentialStruct
            n = n + 1
            emdescriptors[n] = descriptors[i]
        else
            mldescriptors = descriptors[i]
        end
    end

    if n>0
        if length(descriptors) > n
            for i = length(descriptors):-1:(n+1)
                deleteat!(emdescriptors, i)
            end
        end
    else
        emdescriptors = []
    end

    return emdescriptors, mldescriptors
end

function maximumcutoff(descriptors)

    emdescriptors, mldescriptor = getdescriptors(descriptors)

    rcutmax = 0.0
    for i = 1:length(emdescriptors)
        rcutmax = max(rcutmax, emdescriptors[i].rcut)
    end
    if mldescriptor !== nothing
        if lowercase(mldescriptor.name) == "snap"
            rcutmax = max(mldescriptor.rcutmax,rcutmax)
        end
    end

    return rcutmax
end

function emlpotential(x, q, t, alist, neighlist, neighnum, eta, kappa, descriptors)

    empotential, mlpotential = getdescriptors(descriptors)
    etot, e, f, v, vtot = empiricalpotential(x, q, t, alist, neighlist, neighnum, eta, kappa, empotential)

    if mlpotential !== nothing
        eatom, fatom, vatom  = snappotential(x, t, alist, neighlist, neighnum, mlpotential)
        e = e + eatom
        f = f + fatom
        v = v + vatom
        etot = sum(e);
        vtot = sum(v,dims=2);
    end

    return etot, e, f, v, vtot
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using JuMP, CSV, DataFrames
import GLPK

function solve(env)
    data = env["data"]

    LT = env["LT"]

    MC = keys(env["data"])

    m = Model(with_optimizer(GLPK.Optimizer, tm_lim = 60000))

    @variables m begin
        assignment[ lt = LT, mc = [m for m in MC], periods = 1:3] >= 0, Bin
        z
    end

    T = 1:3

    for t in T
        for mc in MC
            @constraint(m, sum(assignment[:, mc, t]) == 1)
        end

        for lt in LT
            @constraint(m, sum(assignment[lt, :, t]) <= 1)
        end
    end

    for lt in LT
        for mc in MC
            @constraint(m, sum(assignment[lt, mc, :]) <= 1)
        end
    end

    for mc in MC
        @constraint(m, sum(assignment[data[mc].alreadyAssigned, mc, :]) ==0)
    end

    for mc in MC
        if data[mc].enableAssignmentToCurrentLT != true
            @constraint(m, sum(assignment[data[mc].currentLT, mc, :]) == 0)
        end
    end

    for mc in MC
        @constraint(m, sum(data[mc].preferences[lt] * assignment[lt, mc, t] for lt in LT, t in T) <= z)
    end

    #Make worst person as good as possible
    @objective(m, Min, z)

    optimize!(m)
    println(value(z))

    @constraint(m, z <= objective_value(m))

    #Add constraint so they don't get any worse

    #Make everyone as good as possible
    @objective(m, Min, sum(data[mc].preferences[lt] * assignment[lt, mc, t] for lt in LT, mc in MC, t in T))
    optimize!(m)

    f = open("FileName.csv", "w")

    results = Dict{String, Any}()

    # output = DataFrame(index = [mc for mc in MC], columns = [t for t in T])
    # println(output)
    columns = [Symbol(col) => String[] for col in T]
    df = DataFrame(fill(String, length(MC)), Symbol.(MC), length(T))
    df[2, Symbol.("Jonathan Symons")] = "5"
    println("HERE")
    println(df)
    for mc in MC
        println(mc)
        # println(sum(data[mc].preferences[lt] * value(assignment[lt, mc, t]) for lt in LT, t in T))
        results[mc] = Dict{Int, String}()
        for lt in LT

            for t in T
                if value.(assignment[lt, mc, t]) == 1
                    println(assignment[lt, mc, t])
                    # results[mc][t] = lt
                    # output[mc,t] = lt
                    df[t, Symbol.(mc)] = lt
                end
            end
        end
        # push!(df, results[mc])
        # println(results[mc])
        # df
    end
    println(df)
    # itr = (merge((; key = key), val) for (key, val) in results)
    # itr |> CSV.write("test.csv")

    println("Objective is ", objective_value(m))

    # CSV.write(f, DataFrame(results))
    CSV.write(f, df)
    close(f)

    # println(output)
    # output.to_csv(path_or_buf = "output.csv")

    return m
end


text:
using Test
using HolidayCalendars
using Dates: Date, dayofweek
using RDates

@testset "Holiday Calendars" verbose=true begin
    @testset "United Kingdom → Settlements" begin
        holidays_2004 = Set([
            Date(2004,1,1),
            Date(2004,4,9),
            Date(2004,4,12),
            Date(2004,5,3),
            Date(2004,5,31),
            Date(2004,8,30),
            Date(2004,12,27),
            Date(2004,12,28),
        ])

        cal = RDates.calendar(CALENDARS, ["UK/SETTLEMENT"])
        for d in range(Date(2004,1,1), Date(2004,12,31), rd"1b@WEEKEND"; cal_mgr=CALENDARS)
            @test is_holiday(cal, d) == (d in holidays_2004 || dayofweek(d) > 5)
        end

        holidays_2005 = Set([
            Date(2005,1,3),
            Date(2005,3,25),
            Date(2005,3,28),
            Date(2005,5,2),
            Date(2005,5,30),
            Date(2005,8,29),
            Date(2005,12,26),
            Date(2005,12,27),
        ])

        for d in range(Date(2005,1,1), Date(2005,12,31), rd"1b@WEEKEND"; cal_mgr=CALENDARS)
            @test is_holiday(cal, d) == (d in holidays_2005 || dayofweek(d) > 5)
        end

        holidays_2021 = Set([
            Date(2021,1,1),
            Date(2021,4,2),
            Date(2021,4,5),
            Date(2021,5,3),
            Date(2021,5,31),
            Date(2021,8,30),
            Date(2021,12,27),
            Date(2021,12,28),
        ])

        for d in range(Date(2021,1,1), Date(2021,12,31), rd"1b@WEEKEND"; cal_mgr=CALENDARS)
            @test is_holiday(cal, d) == (d in holidays_2021 || dayofweek(d) > 5)
        end

        holidays_2022 = Set([
            Date(2022,1,3),
            Date(2022,4,15),
            Date(2022,4,18),
            Date(2022,5,2),
            Date(2022,6,2),
            Date(2022,6,3),
            Date(2022,8,29),
            Date(2022,12,26),
            Date(2022,12,27),
        ])

        for d in range(Date(2022,1,1), Date(2022,12,31), rd"1b@WEEKEND"; cal_mgr=CALENDARS)
            @test is_holiday(cal, d) == (d in holidays_2022 || dayofweek(d) > 5)
        end
    end
end


text:
using CSTParser

"""
Represents one API alias.

LHS is the undecorated function name
RHS is the API version specific function name
"""
struct KuberAPIAlias
    lhs_fn_name::String
    rhs_fn_name::String
    args::Vector{Pair{String,String}}
    kwargs::Vector{Pair{String,String}}
end

"""
List of all aliases parsed.
Used to maintain state and write out the aliases file at the end.
"""
const KuberAPIAliasesSet = Dict{String,Vector{KuberAPIAlias}}

"""
Returns one `KuberAPIAlias` instance for the function expression.
Returns `nothing` if it is not required to have an alias for the function.
"""
function emit_alias(fn_expr::CSTParser.EXPR, api_decoration::String)
    fn_name = CSTParser.str_value(CSTParser.get_name(fn_expr))
    undec_fn_name = join(split(fn_name, api_decoration))

    (fn_name == undec_fn_name) && (return nothing)
    fn_sig = CSTParser.get_sig(fn_expr)
    fn_args = filter(x->((x.typ === CSTParser.BinaryOpCall) || (x.typ === CSTParser.IDENTIFIER)), fn_sig.args[3:end])
    fn_kwargs_container = filter(x->(x.typ === CSTParser.Parameters), fn_sig.args[3:end])
    fn_kwargs = isempty(fn_kwargs_container) ? [] : filter(x->(x.typ === CSTParser.Kw), fn_kwargs_container[1].args)

    arg_pairs = map(x->CSTParser.str_value(CSTParser.get_name(x))=>((x.typ === CSTParser.IDENTIFIER) ? "" : CSTParser.str_value(x.args[3])), fn_args)
    kwarg_pairs = map(x->CSTParser.str_value(CSTParser.get_name(x.args[1]))=>CSTParser.str_value(x.args[3]), fn_kwargs)

    KuberAPIAlias(undec_fn_name, fn_name, arg_pairs, kwarg_pairs)
end

"""
Returns all API aliases detected in the specified file.
"""
function kuberapi(file::String)
    api_decoration = ""
    aliases = Vector{KuberAPIAlias}()

    x, ps = CSTParser.parse(ParseState(String(readchomp(file))))

    while !ps.done
        if CSTParser.defines_struct(x)
            structsig = CSTParser.get_sig(x)
            if structsig.typ === CSTParser.BinaryOpCall
                op = structsig.args[2]
                if CSTParser.is_issubt(op)
                    subtyp = CSTParser.str_value(CSTParser.get_name(structsig.args[3]))
                    if subtyp == "SwaggerApi"
                        typename = CSTParser.str_value(CSTParser.get_name(structsig.args[1]))
                        if endswith(typename, "Api")
                            api_decoration = typename[1:(end-3)]
                        end
                    end
                end
            end
        elseif !isempty(api_decoration)
            # Note: We are guaranteed to receive the struct definition before methods
            #       because of the sequence in which Swagger code is generated and also
            #       because of the Julia restriction of type being defined before use
            fn_expr = ((x.typ === CSTParser.MacroCall) && CSTParser.defines_function(x.args[3])) ? x.args[3] : CSTParser.defines_function(x) ? x : nothing
            if fn_expr !== nothing
                alias = emit_alias(fn_expr, api_decoration)
                (alias === nothing) || push!(aliases, alias)
            end
        end
        x, ps = CSTParser.parse(ps)
    end

    aliases
end

function detect_aliases(folder::String)
    aliases_set = KuberAPIAliasesSet()

    for file in filter(x->startswith(x, "api_"), readdir(folder))
        @info("    - " * file)
        for alias in kuberapi(joinpath(folder, file))
            @info("        - " * alias.lhs_fn_name * " => " * alias.rhs_fn_name)
            push!(get!(()->Vector{KuberAPIAlias}(), aliases_set, alias.lhs_fn_name), alias)
        end
    end
    @info("got $(length(aliases_set)) sets with $(sum(map(length, values(aliases_set)))) aliases")
    aliases_set
end

function gen_aliases(folder::String, output::String)
    @info("reading $folder/api_*.jl")
    @info("generating $output")

    aliases_set = detect_aliases(folder)

    open("$output", "w") do faliases
        for (api_name, api_aliases) in aliases_set
            println(faliases, "# ", api_name)
            for api_alias in api_aliases
                if isempty(api_alias.args)
                    lhs_args = rhs_args = ""
                else
                    lhs_args = join(map(x->(isempty(x[2]) ? x[1] : join(x,"::")), api_alias.args), ", ")
                    rhs_args = join(map(x->x[1], api_alias.args), ", ")
                end
                if isempty(api_alias.kwargs)
                    lhs_kwargs = rhs_kwargs = ""
                else
                    lhs_kwargs = "; " * join(map(x->join(x,'='), api_alias.kwargs), ", ")
                    rhs_kwargs = "; " * join(map(x->join((x[1],x[1]),'='), api_alias.kwargs), ", ")
                end
                println(faliases, api_alias.lhs_fn_name, "(", lhs_args, lhs_kwargs, ") = ", api_alias.rhs_fn_name, "(", rhs_args, rhs_kwargs, ")")
            end
            println(faliases, "export ", api_name, "\n")
        end
    end
end

function main()
    DIR = dirname(@__FILE__)
    gen_aliases(joinpath(DIR, "../src/api"), joinpath(DIR, "../src/apialiases.jl"))
end

main()


text:
# Compare stacking using circular and plane wave approximations
using Beamforming, Seis, TauPy, Plots
import Printf

function synthetic_prem_arrivals(β, Δ, phase;
		aperture=20, nsta=20, depth=0, pad=100, freq=0.1, delta=1/20freq)
	lon = aperture.*(rand(nsta) .- 0.5)
	lat = aperture.*(rand(nsta) .- 0.5)
	x, y, z, mean_lon, mean_lat = Beamforming.array_geometry(lon, lat, zeros(nsta))
	# Event
	evt = Seis.Event()
	evt.lon, evt.lat = Beamforming.geodesic_endpoint(mean_lon, mean_lat, β, Δ)
	evt.dep = depth
	# Distances
	Δs = Beamforming.delta.(lon, lat, evt.lon, evt.lat)
	# Get arrival times
	arr = try
		first.(TauPy.travel_time.(evt.dep, Δs, phase))
	catch err
		error("cannot find arrivals for $phase")
	end
	times = [a.time for a in arr]
	t1, t2 = extrema(times) .+ (-pad, pad)
	npts = round(Int, (t2 - t1)/delta) + 1
	# Traces
	s = [Seis.Synth.gaussian1(t1, delta, npts, freq, time) for time in times]
	s.sta.lon = lon
	s.sta.lat = lat
	s.sta.elev = 0
	s.sta.net = "AN"
	s.sta.sta = [Printf.@sprintf("A%01d", i) for i in 1:length(s)]
	s.sta.cha = "BXR"
	s.evt = evt
	s
end

s = synthetic_prem_arrivals(45, 120, "SKS", pad=500, freq=0.2, aperture=10)
arf = array_response(s, 1.5, 0.05, 0.01, 0.2, 0.01, true)
bfs = [beamform(s, 1500, 1600, 1, 4, 1, 4, 0.05, method=:linear, n=3,
	            wavefront=wave) for wave in (:plane, :circle)]
colormap = cgrad([:red,:orange,:yellow,:green,:blue,:white] |> reverse)
pl = plot(
	scatter(s.sta.lon, s.sta.lat, aspect_ratio=:equal, label="",
		framestyle=:box, grid=false, m=(:utriangle,:blue)),
	plot(arf, c=colormap),
	plot.(bfs, phases="SKS", c=colormap)...,
	title=["Array map" "ARF" "Plane" "Circle"],
	size=(700,700))

max_stars_repo_name: anowacki/Beamforming

max_issues_repo_name: anowacki/Beamforming

max_forks_repo_name: anowacki/Beamforming


text:
#!/usr/bin/julia --color=yes

for n in 1:30
    header=
    """
    #!/bin/sh
    #\$ -o \$HOME/corrthing_$(n).out -j y
    #\$ -N corrthing_$(n)

    cd tfg/

    ./exe/annealer -p -t 1e6 -T 1e8 -n 100 -i 1e9 -m 100 0.8 corrthing_$(n).net

    """
    open("job_$n","w") do io
        println(io,header)
    end
end

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg


text:
using CVChannel
using Test

"""
This script verifies the multiplicativity of the depolarizing channel
with itself for small dimensions.
"""

@testset "Verifying Depolarizing Channel Multiplicativity for d=3" begin
    println("Here we verify multiplicativity of the depolarizing channel with itself for d=3.")
    d = 3
    scan_range = [0:0.05:1;]
    is_mult = true
    max_val = 0
    for q in scan_range
        println("Now checking q = ", q)
        depolqChan(X) = depolarizingChannel(X,q)
        depolq_chan = Choi(depolqChan,d,d)

        par_cv, = pptCV(parChoi(depolq_chan, depolq_chan), :dual)

        target = (d * (1-q) + q)^2
        !isapprox(par_cv,target,atol=3e-6) ? is_mult = false : nothing
        max_val < abs(par_cv-target) ? max_val = abs(par_cv-target) : nothing
    end
    @test is_mult
end

"""
This verifies multiplicativity for d=4 but it takes forever
and doesn't matter for running in parallel with other channels
because it is too big.

@testset "Verifying Multiplicativity for d=4" begin
    d = 4
    scan_range = [0.75:0.05:1;]
    is_mult = true
    max_val = 2.264111881089e-6
    for q in scan_range
        println("Now checking q =", q, ".")
        depolqChan(X) = depolarizingChannel(X,q)
        depolq_chan = Choi(depolqChan,d,d)

        par_dims = d*[1,1,1,1]
        par_JN = permuteSubsystems(kron(depolq_chan.JN, depolq_chan.JN), [1,3,2,4], par_dims)
        par_in_dim = d^2
        par_out_dim = d^2
        par_cv, = pptCVDual(par_JN, d^2, d^2)

        target = (d*(1-q)+q)^2
        !isapprox(par_cv,target,atol=3e-6) ? is_mult = false : nothing
        max_val < abs(par_cv-target) ? max_val = abs(par_cv-target) : nothing
    end
    @test is_mult
end

"""

max_forks_repo_name: ChitambarLab/cv-channel


text:
#julia
using Makie, GeometryTypes, Colors, MacroTools
scene = Scene()
vx = -1:0.1:1;
vy = -1:0.1:1;

f(x, y) = (sin(x*10) + cos(y*10)) / 4
psurf = surface(vx, vy, f)

pos = lift_node(psurf[:x], psurf[:y], psurf[:z]) do x, y, z
    vec(Point3f0.(x, y', z .+ 0.5))
end
pscat = scatter(pos)
plines = lines(view(pos, 1:2:length(pos)))
scene

@theme theme = begin
    markersize = to_markersize(0.01)
    strokecolor = to_color(:white)
    strokewidth = to_float(0.01)
end
# this pushes all the values from theme to the plot
# push!(pscat, theme)
# Or update the entire surface node with this
# scene[:scatter] = theme
# Or permananently (to be more precise: just for this session) change the theme for scatter
scene[:theme, :scatter] = theme
scatter(lift_node(x-> x .+ (Point3f0(0, 0, 1),), pos)) # will now use new theme
scene

# Make a completely new theme
function custom_theme(scene)
    @theme theme = begin
        linewidth = to_float(3)
        colormap = to_colormap(:BuGn)
        scatter = begin
            marker = to_spritemarker(Circle)
            markersize = to_float(0.03)
            strokecolor = to_color(:white)
            strokewidth = to_float(0.01)
            glowcolor = to_color(RGBA(0, 0, 0, 0.4))
            glowwidth = to_float(0.1)
        end
    end
    # update theme values
    scene[:theme] = theme
end

# apply it to the scene
custom_theme(scene)

# From now everything will be plotted with new theme
psurf = surface(vx, 1:0.1:2, psurf[:z])


text:
module RamirezReyes_Yang_SpontaneousCyclogenesis

using Distributed: @everywhere, @spawnat, fetch
#using ClusterManagers
using ImageSegmentation
using AvailablePotentialEnergyFramework
using JLD
using HDF5
using NCDatasets
using Statistics

#addprocs(SlurmManager(9),m="cyclic")

# Write your package code here.

include("APE_calculation_from_sam_output.jl")
include("APE_calculation_from_sam_output_nosmoothing.jl")
include("APE_calculation_from_sam_output_addcondensateloading.jl")
include("Composite_creation_from_sam_output_and_diagnostics.jl")
include("Composite_creation_from_sam_output_and_diagnostics_nosmoothing.jl")
include("Composite_creation_from_sam_output_and_diagnostics_in_chunks_9hpa_15days.jl")

export getapeanalysis,getapeanalysis_nosmoothing,get_composites, getapeanalysis_last15days
export get_composites_nosmoothing
export smooth_vars_and_write_to_netcdf!, cyclone_detect_and_save, get_composites_in_chunks
#import animate_tests.jl

end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# contributed by Simon Danisch
# based on the [C++ implementation](https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/binarytrees-gpp-9.html)

using Printf
struct Node
    left::UInt32
    right::UInt32
end
function alloc!(pool, left, right)
    push!(pool, Node(left, right))
    return length(pool)
end
function make(pool, d)
    d == 0 && return 0
    alloc!(pool, make(pool, d - 1), make(pool, d - 1))
end
check(pool, t::Node) = 1 + check(pool, t.left) + check(pool, t.right)
function check(pool, node::Integer)
    node == 0 && return 1
    @inbounds return check(pool, pool[node])
end
function threads_inner(pool, d, min_depth, max_depth)
    niter = 1 << (max_depth - d + min_depth)
    c = 0
    for j = 1:niter
        c += check(pool, make(pool, d))
        empty!(pool)
    end
    return (niter, d, c)
end

function loop_depths(io, d, min_depth, max_depth, ::Val{N}) where N
    threadstore = ntuple(x-> NTuple{3, Int}[], N)
    Threads.@threads for d in min_depth:2:max_depth
        pool = Node[]
        push!(threadstore[Threads.threadid()], threads_inner(pool, d, min_depth, max_depth))
    end
    for results in threadstore, result in results
        @printf(io, "%i\t trees of depth %i\t check: %i\n", result...)
    end
end
function perf_binary_trees(io, N::Int=10)
    min_depth = 4
    max_depth = N
    stretch_depth = max_depth + 1
    pool = Node[]
    # create and check stretch tree
    c = check(pool, make(pool, stretch_depth))
    @printf(io, "stretch tree of depth %i\t check: %i\n", stretch_depth, c)

    long_lived_tree = make(pool, max_depth)

    loop_depths(io, min_depth, min_depth, max_depth, Val(Threads.nthreads()))
    @printf(io, "long lived tree of depth %i\t check: %i\n", max_depth, check(pool, long_lived_tree))
end

perf_binary_trees(stdout, parse(Int, ARGS[1]))


text:
include("utils.jl")

struct ProjManifold <: Manifold end

ManifoldsBase.inner(::ProjManifold, x, w, v) = dot(w, v)
ManifoldsBase.project_tangent!(S::ProjManifold, w, x, v) = (w .= v .- dot(x, v) .* x)
ManifoldsBase.representation_size(::ProjManifold) = (2,3)
ManifoldsBase.manifold_dimension(::ProjManifold) = 5
Manifolds.get_vector(::ProjManifold, x, v, ::ArbitraryOrthonormalBasis) = reverse(v)

@testset "Projected and arbitrary orthonormal basis" begin
    M = ProjManifold()
    x = [sqrt(2)/2 0.0 0.0;
         0.0 sqrt(2)/2 0.0]

    pb = get_basis(M, x, ProjectedOrthonormalBasis(:svd))
    @test number_system(pb) == ℝ
    @test get_basis(M, x, pb) == pb
    N = manifold_dimension(M)
    @test isa(pb, PrecomputedOrthonormalBasis)
    @test length(pb.vectors) == N
    # test orthonormality
    for i in 1:N
        @test norm(M, x, pb.vectors[i]) ≈ 1
        for j in i+1:N
            @test inner(M, x, pb.vectors[i], pb.vectors[j]) ≈ 0 atol = 1e-15
        end
    end
    # check projection idempotency
    for i in 1:N
        @test project_tangent(M, x, pb.vectors[i]) ≈ pb.vectors[i]
    end

    aonb = get_basis(M, x, ArbitraryOrthonormalBasis())
    @test size(aonb.vectors) == (5,)
    @test aonb.vectors[1] ≈ [0, 0, 0, 0, 1]
end

struct NonManifold <: Manifold end
struct NonBasis <: Manifolds.AbstractBasis{ℝ} end

@testset "ManifoldsBase.jl stuff" begin

    @testset "Errors" begin
        m = NonManifold()
        onb = ArbitraryOrthonormalBasis()

        @test_throws ErrorException get_basis(m, [0], onb)
        @test_throws ErrorException get_basis(m, [0], NonBasis())
        @test_throws ErrorException get_coordinates(m, [0], [0], onb)
        @test_throws ErrorException get_vector(m, [0], [0], onb)
        @test_throws ErrorException get_vectors(m, [0], NonBasis())
    end

    M = ManifoldsBase.DefaultManifold(3)
    pts = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    @testset "basis representation" begin
        v1 = log(M, pts[1], pts[2])

        vb = get_coordinates(M, pts[1], v1, ArbitraryOrthonormalBasis())
        @test isa(vb, AbstractVector)
        vbi = get_vector(M, pts[1], vb, ArbitraryOrthonormalBasis())
        @test isapprox(M, pts[1], v1, vbi)

        b = get_basis(M, pts[1], ArbitraryOrthonormalBasis())
        @test isa(b, PrecomputedOrthonormalBasis)
        N = manifold_dimension(M)
        @test length(b.vectors) == N
        # check orthonormality
        for i in 1:N
            @test norm(M, pts[1], b.vectors[i]) ≈ 1
            for j in i+1:N
                @test inner(M, pts[1], b.vectors[i], b.vectors[j]) ≈ 0
            end
        end
        # check that the coefficients correspond to the basis
        for i in 1:N
            @test inner(M, pts[1], v1, b.vectors[i]) ≈ vb[i]
        end

        @test get_coordinates(M, pts[1], v1, b) ≈ get_coordinates(M, pts[1], v1, ArbitraryOrthonormalBasis())
        @test get_vector(M, pts[1], vb, b) ≈ get_vector(M, pts[1], vb, ArbitraryOrthonormalBasis())

    end

    @testset "ArrayManifold basis" begin
        A = ArrayManifold(M)
        aonb = ArbitraryOrthonormalBasis()
        b = get_basis(A, pts[1], aonb)
        @test_throws ErrorException get_vector(A, pts[1], [], aonb)
        @test_throws DimensionMismatch get_coordinates(A, pts[1], [], aonb)
        @test_throws ArgumentError get_basis(A, pts[1], PrecomputedOrthonormalBasis([pts[1]]))
        @test_throws ArgumentError get_basis(A, pts[1], PrecomputedOrthonormalBasis([pts[1], pts[1], pts[1]]))
        @test_throws ArgumentError get_basis(A, pts[1], PrecomputedOrthonormalBasis([2*pts[1], pts[1], pts[1]]))
    end
end


text:
module FrictionModels

using ..NQCModels: NQCModels
using ..AdiabaticModels: AdiabaticModel

export friction, friction!

"""
    AdiabaticFrictionModel <: AdiabaticModel

`AdiabaticFrictionModel`s must implement `potential!`, `derivative!`, and `friction!`

`potential!` and `friction!` should be the same as for the `AdiabaticModel`.

`friction!` must fill an `AbstractMatrix` with `size = (ndofs*natoms, ndofs*natoms)`.
"""
abstract type AdiabaticFrictionModel <: AdiabaticModel end

"""
    ElectronicFrictionProvider

Abstract type for defining models that provide electronic friction only.
Subtypes of this should implement `friction!` and `ndofs`.
"""
abstract type ElectronicFrictionProvider end

"""
    friction!(model::AdiabaticFrictionModel, F, R:AbstractMatrix)

Fill `F` with the electronic friction as a function of the positions `R`.

This need only be implemented for `AdiabaticFrictionModel`s.
"""
function friction! end

"""
    friction(model::Model, R)

Obtain the friction for the current position `R`.

This is an allocating version of `friction!`.
"""
function friction(model::AdiabaticFrictionModel, R)
    F = zero_friction(model, R)
    friction!(model, F, R)
    return F
end

zero_friction(::AdiabaticFrictionModel, R) = zeros(eltype(R), length(R), length(R))

include("composite_friction_model.jl")
export CompositeFrictionModel
include("ase_friction_interface.jl")
export ASEFrictionProvider

include("constant_friction.jl")
export ConstantFriction
include("random_friction.jl")
export RandomFriction

end # module