text:
## --- Import packages we'll need

    using StatGeochem, Statistics
    using Plots; gr(); default(fmt = :svg)

## --- Read data files

    kernel_sigma = 330 # Myr

    abbott = importdataset("data/Abbott.psv", '|')
    abbott["Age_sigma"] = kernel_sigma * ones(size(abbott["Age"]))

    herzberg = importdataset("data/Herzberg.psv",'|')
    herzberg["Age_sigma"] = kernel_sigma * ones(size(herzberg["Age"]))

    ganne = importdataset("data/Ganne_FeTi.psv",'|')
    ganne["Age_sigma"] = kernel_sigma * ones(size(ganne["Age"]))

    condieHM = importdataset("data/CondieHM.psv",'|')
    condieHM["Age_sigma"] = kernel_sigma * ones(size(condieHM["Age"]))

    condieDM = importdataset("data/CondieDM.psv",'|')
    condieDM["Age_sigma"] = kernel_sigma * ones(size(condieDM["Age"]))

    condieEM = importdataset("data/CondieEM.psv",'|')
    condieEM["Age_sigma"] = kernel_sigma * ones(size(condieEM["Age"]))

    ks2012 = importdataset("data/KellerSchoene.psv",'|')

    kimura = Dict()
    kimura["Age"]  = [68.051, 144.932, 345.154, 454.212, 497.127, 742.077, 802.843, 801.042, 1283.93, 1578.87, 1743.42, 1873.8, 2027.54, 2463.79, 2463.86, 2465.69, 2665.88, 2662.36, 2762.47, 2758.93, 2937.62, 3030.64, 3030.67, 3291.55, 3398.95, 3398.98, 3470.46]
    kimura["T"] = [1389.64, 1396.78, 1368.1, 1368.07, 1379.99, 1406.16, 1371.56, 1348.89, 1629.08, 1550.28, 1663.55, 1438.07, 1412.98, 1443.88, 1558.39, 1626.38, 1541.64, 1635.87, 1625.11, 1687.14, 1533.22, 1609.53, 1660.83, 1417.42, 1642.84, 1684.59, 1627.31]

## --- Plot Herzberg, Kimura and Abbott data and bootstrap-resampled averages

    h1 = plot( xlabel="Age (Ma)", ylabel="T (C)", framestyle=:box,
        xlims=(0,4000), ylims=(1200,1700))

    # Plot dots
    plot!(h1, ganne["Age"], ganne["T_p"], st=:scatter, msalpha=0, label="")
    plot!(h1, abbott["Age"], abbott["T_l"], # yerror=2*abbott["T_l_sigma"],
        st=:scatter, msalpha=0, label="")
    plot!(h1, herzberg["Age"], herzberg["T_p"], st=:scatter, color=:blue, msalpha=0, label="")
    plot!(h1, kimura["Age"], kimura["T"], st=:scatter, color=:darkblue, msalpha=0, label="")

    nresamplings = 10000
    nbins = 10
    tmin = 0
    tmax = 4000

    c = Dict(); m = Dict(); e = Dict();

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    k = invweight_age([ganne["Age"]; herzberg["Age"]])
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the abbott in each resampling
    (c["primelts"], m["primelts"], e["primelts"]) = bin_bsr([ganne["Age"]; herzberg["Age"]],
        [ganne["T_p"]; herzberg["T_p"]], tmin, tmax, nbins, [ganne["Age_sigma"]; herzberg["Age_sigma"]], nresamplings, p)
    plot!(h1, c["primelts"], m["primelts"], yerror=2*e["primelts"], label="", color=:darkblue, mscolor=:auto)

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    k = invweight_age(abbott["Age"])
    p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the abbott in each resampling
    (c["abbott"], m["abbott"], e["abbott"]) = bin_bsr(abbott["Age"],abbott["T_l"], tmin, tmax, nbins, abbott["Age_sigma"],nresamplings,p)
    plot!(h1, c["abbott"], m["abbott"], yerror=2*e["abbott"], label="", color=:darkred, mscolor=:auto)

    # Plot results
    display(h1)
    savefig(h1,"AbbottHerzbergKimuraT.pdf")

## --- Plot Condie data and bootstrap-resampled averages

    # Calculate binned means and uncertainties
    # (c = bincenters, m = mean, el = lower 95% CI, eu = upper 95% CI)
    nbins = 10

    h2 = plot( xlabel="Age (Ma)", ylabel="T (C)", framestyle=:box, xlims=(0,4000), ylims=(1300,1650))
    plot!(h2, condieEM["Age"], condieEM["T_g"], st=:scatter, color=:red, msalpha=0, label="", markersize=2)
    k = invweight_age(condieEM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the abbott in each resampling
    (c["condieEM"], m["condieEM"], e["condieEM"]) = bin_bsr(condieEM["Age"], condieEM["T_g"], tmin, tmax, nbins, condieEM["Age_sigma"], nresamplings, p)
    plot!(h2, c["condieEM"], m["condieEM"], yerror=2*e["condieEM"], label="", color=:red, mscolor=:auto)

    plot!(h2, condieDM["Age"], condieDM["T_g"], st=:scatter, color=:cyan, msalpha=0, label="", markersize=2)
    k = invweight_age(condieDM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the abbott in each resampling
    (c["condieDM"], m["condieDM"], e["condieDM"]) = bin_bsr(condieDM["Age"], condieDM["T_g"], tmin, tmax, nbins, condieDM["Age_sigma"], nresamplings, p)
    plot!(h2, c["condieDM"], m["condieDM"], yerror=2*e["condieDM"], label="", color=:cyan, mscolor=:auto)

    plot!(h2, condieHM["Age"], condieHM["T_g"], st=:scatter, color=:blue, msalpha=0, label="", markersize=2)
    k = invweight_age(condieHM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0) # Keep rougly one-fith of the abbott in each resampling
    (c["condieHM"], m["condieHM"], e["condieHM"]) = bin_bsr(condieHM["Age"], condieHM["T_g"], tmin, tmax, nbins, condieHM["Age_sigma"], nresamplings, p)
    plot!(h2, c["condieHM"], m["condieHM"], yerror=2*e["condieHM"], label="", color=:blue, mscolor=:auto)

    display(h2)
    savefig(h2,"CondieDMEMHM.pdf")

    k = invweight_age([condieHM["Age"]; condieDM["Age"]]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["condie"], m["condie"], e["condie"]) = bin_bsr([condieHM["Age"]; condieDM["Age"]], [condieHM["T_g"]; condieDM["T_g"]], tmin, tmax, nbins, [condieHM["Age_sigma"]; condieDM["Age_sigma"]], nresamplings, p)

## --- Plot relative cooling curves

    h3 = plot(xlabel="Age (Ma)", ylabel="ΔT (C)", framestyle=:box, xlims=(0,4000), ylims=(-30, 230), legend=:topleft, fg_color_legend=:white)
    # plot!(h3, c["condieHM"], m["condieHM"] .- m["condieHM"][1], yerror=2*e["condieHM"], label="condieHM", color=:orange, mscolor=:auto)
    # plot!(h3, c["condieDM"], m["condieDM"] .- m["condieDM"][1], yerror=2*e["condieDM"], label="condieDM", color=:red, mscolor=:auto)
    plot!(h3, c["condie"], m["condie"] .- m["condie"][1], yerror=2*e["condie"], label="condie", color=:orangered, mscolor=:auto)
    plot!(h3, c["abbott"], m["abbott"] .- m["abbott"][1], yerror=2*e["abbott"], label="abbott", color=:purple, mscolor=:auto)
    plot!(h3, ks2012["Age"], ks2012["T_p"] .- ks2012["T_p"][1], yerror=2*ks2012["T_p_sigma"], label="ks2012", color=:darkred, mscolor=:auto)
    plot!(h3, c["primelts"], m["primelts"] .- m["primelts"][1], yerror=2*e["primelts"], label="primelts", color=:blue, mscolor=:auto)
    display(h3)
    savefig(h3,"Trel.pdf")

## --- Calculate boostrap-resampled averages in 100 Myr bins

    nbins = 39
    tmin = 0
    tmax = 3900

    k = invweight_age([ganne["Age"]; herzberg["Age"]]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["primelts"], m["primelts"], e["primelts"]) = bin_bsr([ganne["Age"]; herzberg["Age"]], [ganne["T_p"]; herzberg["T_p"]], tmin, tmax, nbins, [ganne["Age_sigma"]; herzberg["Age_sigma"]], nresamplings, p)

    k = invweight_age(abbott["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["abbott"], m["abbott"], e["abbott"]) = bin_bsr(abbott["Age"],abbott["T_l"], tmin, tmax, nbins, abbott["Age_sigma"],nresamplings,p)

    k = invweight_age(condieDM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["condieDM"], m["condieDM"], e["condieDM"]) = bin_bsr(condieDM["Age"], condieDM["T_g"], tmin, tmax, nbins, condieDM["Age_sigma"], nresamplings, p)

    k = invweight_age(condieHM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["condieHM"], m["condieHM"], e["condieHM"]) = bin_bsr(condieHM["Age"], condieHM["T_g"], tmin, tmax, nbins, condieHM["Age_sigma"], nresamplings, p)

    k = invweight_age(condieDM["Age"]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["condieDM"], m["condieDM"], e["condieDM"]) = bin_bsr(condieDM["Age"], condieDM["T_g"], tmin, tmax, nbins, condieDM["Age_sigma"], nresamplings, p)

    k = invweight_age([condieHM["Age"]; condieDM["Age"]]); p = 1.0 ./ ((k .* median(5.0 ./ k)) .+ 1.0)
    (c["condie"], m["condie"], e["condie"]) = bin_bsr([condieHM["Age"]; condieDM["Age"]], [condieHM["T_g"]; condieDM["T_g"]], tmin, tmax, nbins, [condieHM["Age_sigma"]; condieDM["Age_sigma"]], nresamplings, p)

    Age = c["primelts"]
    T_ave = ones(nbins)
    T_ave_sigma = ones(nbins)
    for i=1:nbins
        T = [m["condieDM"][i], m["condieHM"][i], m["abbott"][i], m["primelts"][i]]
        T_sigma = [e["condieDM"][i], e["condieHM"][i], e["abbott"][i], e["primelts"][i]]
        (T_ave[i], T_ave_sigma[i], mswd) = awmean(T, T_sigma)
    end
    plot(Age,T_ave,yerror=2*T_ave_sigma,mscolor=:auto,label="",framestyle=:box)
    savefig("TrelNon-KS2012.pdf")

## --- Plot weighted-mean Trel

    Age = c["primelts"]
    T_ave = ones(nbins)
    T_ave_sigma = ones(nbins)
    for i=1:nbins
        T = [m["condieDM"][i], m["condieHM"][i], m["abbott"][i], m["primelts"][i], ks2012["T_p"][i]]
        T_sigma = [e["condieDM"][i], e["condieHM"][i], e["abbott"][i], e["primelts"][i], ks2012["T_p_sigma"][i]]
        (T_ave[i], T_ave_sigma[i], mswd) = awmean(T, T_sigma)
    end
    plot(Age,T_ave,yerror=2*T_ave_sigma,mscolor=:auto,label="",framestyle=:box)
    plot!(ks2012["Age"], ks2012["T_p"], yerror=2*ks2012["T_p_sigma"], label="", color=:darkred, mscolor=:auto)
    savefig("TrelCombined.pdf")

## --- Export weighted-mean Trel to csv

    CoolingAverage=Dict()
    CoolingAverage["Age"] = Age
    CoolingAverage["T"] = T_ave
    CoolingAverage["T_sigma"] = T_ave_sigma
    CoolingAverage["Trel"] = T_ave .- T_ave[1]
    CoolingAverage["Trel_sigma"] = T_ave_sigma
    CoolingAverage["elements"] = collect(keys(CoolingAverage))
    exportdataset(CoolingAverage, "data/CoolingAverage.csv", ',')

## --- End of File

max_stars_repo_name: brenhinkeller/SecularCooling

max_issues_repo_name: brenhinkeller/SecularCooling

max_forks_repo_name: brenhinkeller/SecularCooling


text:
using Documenter, PETSc

makedocs(modules = [PETSc])

deploydocs(
  deps = Deps.pip("mkdocs", "python-markdown-math"),
  repo = "github.com/JuliaParallel/PETSc.jl.git",
  julia = "0.4",
)

ext: jl


text:
# "smart" summary and REPL printing

function Base.summary(sys::System)
    return "System (base power $(get_base_power(sys)))"
end

function Base.show(io::IO, ::MIME"text/plain", sys::System)
    show_system_table(io, sys, backend = :auto)

    if IS.get_num_components(sys.data.components) > 0
        show_components_table(io, sys, backend = :auto)
    end

    println(io)
    IS.show_time_series_data(io, sys.data, backend = :auto)
end

function Base.show(io::IO, ::MIME"text/html", sys::System)
    show_system_table(io, sys, backend = :html)

    if IS.get_num_components(sys.data.components) > 0
        show_components_table(io, sys, backend = :html)
    end

    println(io)
    IS.show_time_series_data(io, sys.data, backend = :html)
end

function show_system_table(io::IO, sys::System; kwargs...)
    header = ["Property", "Value"]
    num_components = IS.get_num_components(sys.data.components)
    table = [
        "System Units Base" string(get_units_base(sys))
        "Base Power" string(get_base_power(sys))
        "Base Frequency" string(get_frequency(sys))
        "Num Components" string(num_components)
    ]
    PrettyTables.pretty_table(
        io,
        table,
        header;
        title = "System",
        alignment = :l,
        kwargs...,
    )
end

function show_components_table(io::IO, sys::System; kwargs...)
    static_header = ["Type", "Count", "Has Static Time Series", "Has Forecasts"]
    dynamic_header = ["Type", "Count"]
    components = sys.data.components

    static_types = Vector{DataType}()
    dynamic_types = Vector{DataType}()
    for component_type in keys(components.data)
        if component_type <: DynamicInjection
            push!(dynamic_types, component_type)
        else
            push!(static_types, component_type)
        end
    end
    static_data = Array{Any, 2}(undef, length(static_types), length(static_header))
    dynamic_data = Array{Any, 2}(undef, length(dynamic_types), length(dynamic_header))

    static_type_names = [(IS.strip_module_name(string(x)), x) for x in static_types]
    sort!(static_type_names, by = x -> x[1])
    for (i, (type_name, type)) in enumerate(static_type_names)
        vals = components.data[type]
        has_sts = false
        has_forecasts = false
        for val in values(vals)
            if has_time_series(val, StaticTimeSeries)
                has_sts = true
            end
            if has_time_series(val, Forecast)
                has_forecasts = true
            end
            if has_sts && has_forecasts
                break
            end
        end
        static_data[i, 1] = type_name
        static_data[i, 2] = length(vals)
        static_data[i, 3] = has_sts
        static_data[i, 4] = has_forecasts
    end

    if !isempty(static_types)
        println(io)
        PrettyTables.pretty_table(
            io,
            static_data,
            static_header,
            title = "Static Components",
            alignment = :l;
            kwargs...,
        )
    end

    dynamic_type_names = [(IS.strip_module_name(string(x)), x) for x in dynamic_types]
    sort!(dynamic_type_names, by = x -> x[1])
    for (i, (type_name, type)) in enumerate(dynamic_type_names)
        vals = components.data[type]
        dynamic_data[i, 1] = type_name
        dynamic_data[i, 2] = length(vals)
    end

    if !isempty(dynamic_types)
        println(io)
        PrettyTables.pretty_table(
            io,
            dynamic_data,
            dynamic_header,
            title = "Dynamic Components",
            alignment = :l;
            kwargs...,
        )
    end
end

function Base.summary(tech::DeviceParameter)
    return "$(typeof(tech))"
end

function Base.show(io::IO, ::MIME"text/plain", data::PowerSystemTableData)
    println(io, "$(typeof(data)):")
    println(io, "  directory:  $(data.directory)")
    if !isnothing(data.timeseries_metadata_file)
        println(io, "  timeseries_metadata_file:  $(data.timeseries_metadata_file)")
    end
    println(io, "  base_power:  $(data.base_power)")
    for (field, df) in data.category_to_df
        print(io, "  $field:  ")
        println(io, "$(summary(df))")
    end
end

function Base.show(io::IO, ist::Component)
    print(io, string(nameof(typeof(ist))), "(")
    is_first = true
    for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
        getter_name = Symbol("get_$name")
        if field_type <: IS.TimeSeriesContainer ||
           field_type <: InfrastructureSystemsInternal
            continue
        elseif hasproperty(PowerSystems, getter_name)
            getter_func = getproperty(PowerSystems, getter_name)
            val = getter_func(ist)
        else
            val = getproperty(ist, name)
        end
        if is_first
            is_first = false
        else
            print(io, ", ")
        end
        print(io, val)
    end
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", ist::Component)
    default_units = false
    if !has_units_setting(ist)
        print(io, "\n")
        @warn(
            "SystemUnitSetting not defined, using NATURAL_UNITS for displaying device specification."
        )
        set_units_setting!(
            ist,
            SystemUnitsSettings(100.0, UNIT_SYSTEM_MAPPING["NATURAL_UNITS"]),
        )
        default_units = true
    end
    try
        print(io, summary(ist), ":")
        for (name, field_type) in zip(fieldnames(typeof(ist)), fieldtypes(typeof(ist)))
            obj = getproperty(ist, name)
            getter_name = Symbol("get_$name")
            if (obj isa InfrastructureSystemsInternal) && !default_units
                print(io, "\n   ")
                show(io, MIME"text/plain"(), obj.units_info)
                continue
            elseif obj isa IS.TimeSeriesContainer ||
                   obj isa InfrastructureSystemsType ||
                   obj isa Vector{<:InfrastructureSystemsComponent}
                val = summary(getproperty(ist, name))
            elseif hasproperty(PowerSystems, getter_name)
                getter_func = getproperty(PowerSystems, getter_name)
                val = getter_func(ist)
            else
                val = getproperty(ist, name)
            end
            # Not allowed to print `nothing`
            if isnothing(val)
                val = "nothing"
            end
            print(io, "\n   ", name, ": ", val)
        end
    finally
        if default_units
            set_units_setting!(ist, nothing)
        end
    end
end

"""
Show all components of the given type in a table.

# Arguments
- `sys::System`: System containing the components.
- `component_type::Type{<:Component}`: Type to display. Must be a concrete type.
- `additional_columns::Union{Dict, Vector}`: Additional columns to display.
  The Dict option is a mapping of column name to function. The function must accept
  a component.
  The Vector option is an array of field names for the `component_type`.

Extra keyword arguments are forwarded to PrettyTables.pretty_table.

# Examples
```julia
show_components(sys, ThermalStandard)
show_components(sys, ThermalStandard, Dict("has_time_series" => x -> has_time_series(x)))
show_components(sys, ThermalStandard, [:active_power, :reactive_power])
```
"""
function show_components(
    sys::System,
    component_type::Type{<:Component},
    additional_columns::Union{Dict, Vector} = Dict();
    kwargs...,
)
    show_components(stdout, sys, component_type, additional_columns; kwargs...)
end

function show_components(
    io::IO,
    sys::System,
    component_type::Type{<:Component},
    additional_columns::Union{Dict, Vector} = Dict();
    kwargs...,
)
    IS.show_components(
        io,
        sys.data.components,
        component_type,
        additional_columns;
        kwargs...,
    )
end


text:
const EMPTY_SPOT = ""
# This just prepare the metnet to hold more elements by making
# all the numerical fields larger
function expanded_model(metnet::MetNet{T}, newM::Int, newN::Int) where T
    M, N = size(metnet)
    @assert all((newM, newN) .>= (M, N))
    check_dims(metnet)

    function _similar_copy(col, fill, newdim)
        L = length(col)
        newcol = similar(col, newdim)
        newcol[L:end] .= fill
        newcol[1:L] .= col[1:L]
        return newcol
    end

    net = Dict()

    net[:S] = similar(metnet.S, newM, newN)
    net[:S][M:end, :] .= zero(T)
    net[:S][:, N:end] .= zero(T)
    net[:S][1:M, 1:N] .= metnet.S

    net[:b] = _similar_copy(metnet.b, 0.0, newM)
    net[:c] = _similar_copy(metnet.c, 0.0, newN)
    net[:lb] = _similar_copy(metnet.lb, 0.0, newN)
    net[:ub] = _similar_copy(metnet.ub, 0.0, newN)
    net[:rxns] = _similar_copy(metnet.rxns, EMPTY_SPOT, newN)
    net[:mets] = _similar_copy(metnet.mets, EMPTY_SPOT, newM)

    return MetNet(metnet; reshape = false, net...)
end

function findempty(metnet::MetNet, field::Symbol; check = false)
    idx = findfirst(isequal(EMPTY_SPOT), getfield(metnet, field))
    check && isnothing(idx) && error("no empty spot found in $(field) ")
    idx
end

is_compacted(metnet::MetNet) = !isnothing(findempty(metnet, :mets)) && !isnothing(findempty(metnet, :rxns))

function compacted_model(metnet::MetNet)

    empty_mets = findall(metnet.mets .== EMPTY_SPOT)
    met_idxs = trues(size(metnet, 1))
    met_idxs[empty_mets] .= false

    empty_rxns = findall(metnet.rxns .== EMPTY_SPOT)
    rxn_idxs = trues(size(metnet, 2))
    rxn_idxs[empty_rxns] .= false

    net = Dict()
    net[:S] = metnet.S[met_idxs, rxn_idxs]
    net[:b] = metnet.b[met_idxs]
    net[:c] = metnet.c[rxn_idxs]
    net[:lb] = metnet.lb[rxn_idxs]
    net[:ub] = metnet.ub[rxn_idxs]
    net[:rxns] = metnet.rxns[rxn_idxs]
    net[:mets] = metnet.mets[met_idxs]

    return MetNet(metnet; reshape = false, net...)

end


text:
typealias SpellRecords Dict{Int, IndividualData}

"""
Read spells from `hv_pn`. New indexes are kept with probability
`keep_probability`. Records are inserted into `records` (empty by
default), which is returned.
"""
function read_person_spells!(io::IO, keep_probability::Real;
                             maxiter = Nullable{Int}(),
                             records = SpellRecords())
    readlines(io, maxiter) do line
        (id_field, other_fields) = split(line, ';'; limit=2)
        record = random_push!(records, parse(Int, id_field), keep_probability,
                              () -> IndividualData())
        if !isnull(record)
            spell = tryparse(PersonSpellParser, other_fields)
            if isnull(spell)
                warn("Could not parse line:\n$line\n")
            else
                push!(get(record).person_spells, get(spell))
            end
        end
    end
    records
end

"""
Read AMP spells (eg from `hv_epi_uni_roh`). When the person ID is in
`records`, insert the spell there, otherwise don't parse.
"""
function read_AMP_spells!(io::IO, records::SpellRecords;
                          maxiter = Nullable{Int}())
    readlines(io, maxiter) do line
        (id_field, other_fields) = split(line, ';'; limit = 2)
        person_id = parse(Int, id_field)
        if haskey(records, person_id)
            spell = tryparse(AMPSpellParser, other_fields)
            if !isnull(spell)
                push!(records[person_id].AMP_spells, get(spell))
            end
        end
    end
    records
end

"""
Read VMZ spells (eg from `vmz_v2_int`). When the person ID is in
`records`, insert the spell there, otherwise don't parse.  Will ignore
lines with missing person ID.
"""
function read_VMZ_spells!(io::IO, records::SpellRecords;
                         maxiter = Nullable{Int}())
    readlines(io, maxiter) do line
        (_, id_field, other_fields) = split(line, ';'; limit = 3)
        person_id = tryparse(Int, id_field)
        if !isnull(person_id) && haskey(records, get(person_id))
            # println("parsing line: $(line)")
            # println("other_fields: $(other_fields)")
            # println("id: $(get(person_id))")
            spell = tryparse(VMZSpellParser, other_fields)
            if !isnull(spell)
                push!(records[get(person_id)].VMZ_spells, get(spell))
            end
        end
    end
    records
end


text:
function makejk(n)

m = Int64((n-1)*n/2);
indj = zeros(Int64,m);
indk = zeros(Int64,m);
k1 = 1;
for i = 1:(n-1)
    ind = k1:(k1+(n-i)-1);
    indj[ind] .= i;
    indk[ind] = (i+1):n;
    k1 = k1 + (n-i);
end

return indj, indk

end

function neightripletlist(pairlist, pairnumsum, ilist)

pairnum = Int64.(pairnumsum[2:end]-pairnumsum[1:end-1]);
tripletnum  = Int64.((pairnum.-1).*pairnum/2);

inum = length(ilist);
tripletlistj = [];
tripletlistk = [];

for ii = 1:inum
    n1 = pairnumsum[ii];
    n2 = pairnumsum[ii+1];
    g = pairlist[(n1+1):n2]; # a list of neighbors around atom i

    n = pairnum[ii]
    indj, indk = makejk(n);
    tripletlistj = [tripletlistj; g[indj]];
    tripletlistk = [tripletlistk; g[indk]];
end

tripletlist = [tripletlistj[:] tripletlistk[:]];
tripletnum = [0; cumsum(tripletnum)];

tripletlist = Int64.(tripletlist)
tripletnum  = Int64.(tripletnum)

return tripletlist, tripletnum

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Pkg.Artifacts
using ghr_jll

# See this issue: https://github.com/JuliaGizmos/WebIO.jl/issues/422
interact_base_tag = "v0.10.6"
username = "arthursoprana"

function build_js_bundle(dir)
	deps = Dict(
	    "all.js" => "https://use.fontawesome.com/releases/v5.0.7/js/all.js",
	    "prism.css" => "https://cdn.jsdelivr.net/gh/piever/InteractResources@0.1.0/highlight/prism.css",
	    "prism.js" => "https://cdn.jsdelivr.net/gh/piever/InteractResources@0.1.0/highlight/prism.js",
	    "nouislider.min.js" => "https://cdnjs.cloudflare.com/ajax/libs/noUiSlider/11.1.0/nouislider.min.js",
	    "nouislider.min.css" => "https://cdnjs.cloudflare.com/ajax/libs/noUiSlider/11.1.0/nouislider.min.css",
	    "katex.min.js" => "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.9.0/katex.min.js",
	    "katex.min.css" => "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.9.0/katex.min.css",
	)

    for (dep_name, dep_url) in deps
		download(
            dep_url,
            joinpath(dir, dep_name),
        )
    end
end

hash = create_artifact() do dir
    # Build all JS files in `dir`
	build_js_bundle(dir)
end

# Now that the artifact is built on our machine and we know its content hash,
# let's package it up into a tarball and host it on GitHub releases:
function archive_and_upload()
	local tarball_hash
	mktempdir() do dir
		tarball_hash = archive_artifact(hash, joinpath(dir, "jsbundle.tar.gz"))
		# Use `ghr` from `ghr_jll` to upload this file to the GitHub releases of WebIO.jl
		# Note that I think you will need to export `GITHUB_TOKEN=<token>` before running this.
		# Alternatively, you can skip this step and manually upload the tarball using the browser.
		ghr() do ghr
			run(`$(ghr) -recreate -u $(username) -r InteractBase.jl $(interact_base_tag) $(dir)`)
		end
	end
	return tarball_hash
end
tarball_hash = archive_and_upload()

# Finally, let's bind this out to our Artifact file.
# This relative pathing is assuming that this script lives within `deps/`.
artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")
url = "https://github.com/$(username)/InteractBase.jl/releases/download/$(interact_base_tag)/jsbundle.tar.gz"
bind_artifact!(
	# Where the TOML file is that we're writing this binding into
	artifacts_toml,
	# The name of the artifact, used by the `@artifact_str()` macro later
	"jsbundle",
	# The content-hash of the artifact
	hash;
	# Information about how to download this artifact
	download_info = [
		(url, tarball_hash),
	],
)


text:
using Distributions
using StatsBase

# Let defines new name space
let state = 0
    global _get_next_id() = (state += 1)
 end

mutable struct ZeroIntelligentInvestor <: AbstractInvestor
    name::String
    positions::Dict{AbstractAsset, Int64}
    reserved::Dict{AbstractAsset, Int64}
    std_price::Float64

    side_weights::Dict{Type, Union{Int64, Float64}}

    active_orders:: Dict{AbstractMarket, LimitOrder}
    function ZeroIntelligentInvestor(positions::Dict{AbstractAsset, Int64}; deviation::Union{Int64, Float64}=10, bid_weight::Real=1, ask_weight::Real=1, na_weight::Real=1) # , assets::Dict{}
        new(
            string(_get_next_id()),
            positions,
            Dict(),
            Float64(deviation),
            Dict{Type, Union{Int64, Float64}}(BidLimitOrder=>bid_weight, AskLimitOrder=>ask_weight, Nothing=>na_weight),
            Dict()
            )
    end
    function ZeroIntelligentInvestor(cash, stock; deviation::Union{Int64, Float64}) # , assets::Dict{}
        new(
            _get_next_id(),
            Dict(
                    generic_currency=>cash,
                    generic_stock=>stock),
                    Dict(),
                    Float64(deviation),
                    Dict{Type, Union{Int64, Float64}}(),
                    Dict()
                )
    end
end

"Decide the side of the order the trader does"
function get_side(trader::ZeroIntelligentInvestor, market::AbstractMarket)
    options = [AskLimitOrder, BidLimitOrder, Nothing]
    weights = Weights([get(trader.side_weights, option, 1) for option in options]) #[Dict("a" => 5, "b" => 10, "c" => 1)[elem] for elem in ["a", "b", "c"]]
    return sample(options, weights)
end

"Not strictly the order quantity but the quantity of
the asset the investor is trading away"
function get_order_quantity(trader::ZeroIntelligentInvestor, market::AbstractMarket, asset::AbstractAsset)

    max_tradeable = get_unreserved(trader, asset, exclude=market)
    # Note: max_tradeable is automatically in the asset which
    # is traded away (Sell => stock, buy => currency)

    if max_tradeable <= 0
        return NaN
        #throw(DomainError(max_tradeable, "Max quantity cannot be <0: $trader"))
    end
    distr = 1:Int(max_tradeable)
    #println(max_tradeable)
    return rand(distr)
end

function get_order_price(trader::ZeroIntelligentInvestor, market::AbstractMarket, asset::AbstractAsset, min_price::Union{Int64,Float64}, max_price::Union{Int64,Float64})

    # Define the order price
    if isnan(market.last_price)
        # We pick something
        #mean_price = mean([min_price, get_unreserved(trader, asset, exclude=market)])
        mean_price = get_unreserved(trader, market.currency, exclude=market) / get_unreserved(trader, market.traded_asset, exclude=market)
    else
        mean_price = market.last_price
    end
    std_price = trader.std_price #mean([min_price, max_price]) / 100

    distr = truncated(Normal(mean_price, std_price), min_price, max_price)

    return rand(distr)
end

"Decide the order to do to the market"
function get_order(trader::ZeroIntelligentInvestor, market::AbstractMarket)

    # Define side
    side = get_side(trader, market)
    if side == Nothing
        return
    end

    # Define price limits and assets
    min_price = 1
    if side == AskLimitOrder
        from_asset = market.traded_asset
        to_asset = market.currency
        max_price = Inf

    elseif side == BidLimitOrder
        from_asset = market.currency
        to_asset = market.traded_asset
        max_price = get_unreserved(trader, from_asset, exclude=market) # Absolute maximum is such that the investor can buy one unit

    else
        throw(DomainError(side, "Not implmented placement"))
    end

    if max_price <= min_price
        #println("Seems investor $(trader.name) is out of asset $from_asset: $(get_unreserved(trader, from_asset, exclude=market)) (min: $min_price max: $max_price)")
        return
    end

    # Max quantity the trader can trade the asset it
    # is transacting away; either the market currency
    # or market asset

    # First we define the price from normal distribution
    # then the quantity from uniform distribution
    order_price = get_order_price(trader, market, from_asset, min_price, max_price)
    order_quantity = get_order_quantity(trader, market, from_asset)

    if isnan(order_price) | isnan(order_quantity)
        # Cannot make order
        return nothing
    end

    price = Int64(round(order_price))
    if side == BidLimitOrder
        # Convert the allocating currency
        # to the quantity of the order
        # (amount of stock to buy)

        order_quantity = order_quantity / price

    end

    quantity = Int64(floor(order_quantity))

    #println("Creating order $side $price * $quantity (has $(trader.positions))")
    return side(trader, price=price, quantity=quantity, market=market)
end


text:
module TestAsType

using Base.Test
using ArbFloats

af0 = zero(ArbFloat)
af1 = one(ArbFloat)
af2 = af1 + af1
afhalf = af1 / 2
afNaN = ArbFloat(NaN)
afInf = ArbFloat(Inf)

@test isnan(afNaN)
@test isinf(afInf)
@test isfinite(af1)
@test !isfinite(afInf)

@test af1 === af1
@test af1 == af1
@test af1 != af2

@test typeof(ArbFloat{53}(af2)) == ArbFloat{53}
@test typeof(af2) == ArbFloat{precision(ArbFloat)}

@test string(ArbFloat{24}(pi)) == "3.141593"
@test string(ArbFloat{53}(catalan)) == "0.915965594177219"

@test promote(ArbFloat(4), 8//2) == (ArbFloat(4), ArbFloat(4))

end # module


text:
function ComposeMatrices{T<:Signed, F<:AbstractFloat}(N::T, Y::Array{F}; useCachedMatrices::Bool = false)

    readtype = typeof(Y[1])
    if readtype == Float32
        N = Int32(N)
    else
        N = Int64(N)
    end

    if useCachedMatrices
        check = isfile(string(N, "mu.csv")) && isfile(string(N, "sigma.csv")) && isfile(string(N, "U.csv")) && isfile(string(N, "L.csv")) && isfile(string(N, "X.csv")) &&  isfile(string(N, "GMat.csv")) && isfile(string(N, "sigma2.csv"))

        if check
            println("Reading pre-calculated matrices.")

            mu1         = readcsv(string(N, "mu.csv"), readtype)
            mu1         = mu1[1]
            sigma1      = readcsv(string(N, "sigma.csv"), readtype)
            sigma1      = sigma1[1]

            U           = readcsv(string(N, "U.csv"), readtype)
            L           = readcsv(string(N, "L.csv"), readtype)
            X           = readcsv(string(N, "X.csv"), readtype)
            GMat        = readcsv(string(N, "GMat.csv"), readtype)
            arraySigma2 = readcsv(string(N, "sigma2.csv"), readtype)
        else
            println("Pre-Calculating Matrices.. This might take a while.")

            ## Calculate Matrices
            U, L, mu1, sigma1, X, GMat, arraySigma2 = calculate_matrices(N)
            gc()

            writecsv(string(N, "U.csv"), U)
            writecsv(string(N, "L.csv"), L)
            writecsv(string(N, "mu.csv"), mu1)
            writecsv(string(N, "sigma.csv"), sigma1)
            writecsv(string(N, "X.csv"), X)
            writecsv(string(N, "GMat.csv"), GMat)
            writecsv(string(N, "sigma2.csv"), arraySigma2)
        end
    else
        println("Pre-Calculating Matrices.. This might take a while.")
        U, L, mu1, sigma1, X, GMat, arraySigma2 = calculate_matrices(N)
        gc()
    end

    println("Matrices finished. Ready to run LASSO path.");

    ## Compute maximum lambda
    arrayInnerY = smart_innerY(Y, N, U, L, sigma1, mu1)
    # TODO - verify lambda_max
    lambda_max, lambda_max_var = findmax(arraySigma2.*abs(arrayInnerY))
    lambda_max = lambda_max / arraySigma2[lambda_max_var] / N

    return X, GMat, arrayInnerY, arraySigma2, lambda_max
end


text:
using Callsite

@callsite typed function spirv(body, args)
    checked_arg_expr = check_and_convert_for_gpu(args)
    spirv = juliatyped2spirv(body)
    kernel = spirv2binary(spirv)
    quote
        # <- insert some setup here, maybe dependant on the args
        launch($kernel, $(checked_arg_expr.args...))
        # <- insert some teardown here, maybe dependant on the args
    end
end

#similarly
@callsite llvm function ptx(func_llvm::String)
    kernel = ptx2binary(llvmir2ptx(func_llvm))
    quote
        # <- insert some setup here, maybe dependant on the args
        launch($kernel, $(checked_arg_expr.args...))
        # <- insert some teardown here, maybe dependant on the args
    end
end

function map(f, a::CUDAArray, b::CUDAArray)
    _call(ptx, f, a, b)
end

function map(f, a::VulkanArray, b::VulkanArray)
    _call(spirv, f, a, b)
end


text:
import Base: close, unsafe_read, unsafe_write

mutable struct RecordWithSubrecords{C} <: Record
   io        :: IO      # underlying I/O stream
   maxsrlen  :: Int32   # maximum subrecord length
   subreclen :: Int32   # length of current subrecord
   subleft   :: Int32   # bytes left in current subrecord
   more      :: Bool    # whether more subrecords will follow
   writable  :: Bool    # whether this record is writable
   isfirst   :: Bool    # if writable: whether this is the first subrecord
   totleft   :: Int64   # if writable: bytes left in this record
   convert   :: C       # convert method
end

function rdmarker(io, fn)
   reclen = read(io, Int32) |> fn
   more = (reclen < 0)
   return (abs(reclen), more)
end
rdmarker(rec::Record) = rdmarker(rec.io, rec.convert.onread)

function mkmarker(towrite, msl)
   subreclen = min(towrite, msl)
   more = (subreclen < towrite)
   return (Int32(subreclen), more)
end

function wrmarker(io, fn, subreclen, sign)
   marker = subreclen * (sign ? -1 : 1)
   write(io, convert(Int32, marker) |> fn)
end
wrmarker(rec::Record, subreclen, sign) = wrmarker(rec.io, rec.convert.onwrite, subreclen, sign)

function Record( f::FortranFile{SequentialAccess{WithSubrecords}} )
## constructor for readable records
   conv = f.convert
   msl = f.acctyp.recmrktyp.max_subrecord_length
   subreclen, more = rdmarker(f.io, conv.onread) # read leading record marker
   RecordWithSubrecords(f.io, msl, subreclen, subreclen, more, false, false, 0, conv)
end

function Record( f::FortranFile{SequentialAccess{WithSubrecords}}, towrite::Integer )
## constructor for writable records
   conv = f.convert
   msl = f.acctyp.recmrktyp.max_subrecord_length
   subreclen, more = mkmarker(towrite, msl)
   wrmarker(f.io, conv.onwrite, subreclen, more) # write leading record marker
   RecordWithSubrecords(f.io, msl, subreclen, subreclen, more, true, true, towrite, conv)
end

function advance!( rec::RecordWithSubrecords )
   subreclen, sign = rdmarker(rec) # read trailing record marker
   if subreclen != rec.subreclen
      fthrow("trailing subrecord marker doesn't match")
   end
   if rec.more
      subreclen, more = rdmarker(rec) # read next leading subrecord marker
      rec.subleft = rec.subreclen = subreclen
      rec.more = more
   end
   nothing
end

function unsafe_read( rec::RecordWithSubrecords, p::Ptr{UInt8}, n::UInt )
   while (n>0)
      if rec.subleft==0
         fthrow("attempting to read beyond record end")
      end
      toread = min(n, rec.subleft)
      unsafe_read(rec.io, p, toread)
      p += toread
      n -= toread
      rec.subleft -= toread
      if rec.subleft==0; advance!(rec); end
   end
   nothing
end

function unsafe_write( rec::RecordWithSubrecords, p::Ptr{UInt8}, n::UInt )
   nwritten = 0
   while (n>0)
      if rec.totleft==0; fthrow("attempting to write beyond record end"); end
      towrite = min(n, rec.subleft)
      nwritten += unsafe_write(rec.io, p, towrite)
      p += towrite
      n -= towrite
      rec.subleft -= towrite
      rec.totleft -= towrite
      if rec.subleft==0
         # the current subrecord is full
         # if there are more subrecords, write the record markers
         if rec.more
            # write trailing subrecord marker
            wrmarker(rec, rec.subreclen, !rec.isfirst)
            # write next leading subrecord marker
            subreclen, more = mkmarker(rec.totleft, rec.maxsrlen)
            wrmarker(rec, subreclen, more)
            # update our Record
            rec.subleft = rec.subreclen = subreclen
            rec.more = more
            rec.isfirst = false
         end
      end
   end
   return nwritten
end

function close( rec::RecordWithSubrecords )
   if rec.writable
      @assert rec.totleft == 0
      @assert rec.subleft == 0
      @assert !rec.more
      wrmarker(rec, rec.subreclen, !rec.isfirst)
   else
      while rec.subleft > 0
         skip(rec.io, rec.subleft)
         rec.subleft = 0
         advance!(rec)
      end
      @assert !rec.more
   end
   nothing
end


text:
using FRGRealTime,Plots

u0 = [-8.0, 116585.3]

sol=FRGRealTime.ZeroLPASolve(
    145.0,
    4.0,
    1.0,
    800.0,
    u0 = [-8.0, 116585.3],
)

-sol[end][2]

writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_lpa_m2k.dat",-hcat(sol.u...)[2,:])
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_lpa_lamdak.dat",2*hcat(sol.u...)[1,:])
writedlm("/home/tyy/Documents/CTP-fRG-Test/Real-Time-data/realtime_lpa_k.dat",sol.t)


text:
"""
    solve!()

"""
function solve!(model::CPModel, strategy::T=DFSearch(); variableHeuristic=MinDomainVariableSelection(), valueSelection=BasicHeuristic(), out_solver::Bool=false) where T <: SearchStrategy
    return search!(model, strategy, variableHeuristic, valueSelection, out_solver=out_solver)
end


text:
module AbstractMCMC

import BangBang
import ConsoleProgressMonitor
import LoggingExtras
import ProgressLogging
import StatsBase
import TerminalLoggers
import Transducers
using Setfield

import Distributed
import Logging
import Random

# Reexport sample
using StatsBase: sample
export sample

# Parallel sampling types
export MCMCThreads, MCMCDistributed, MCMCSerial

"""
    AbstractChains

`AbstractChains` is an abstract type for an object that stores
parameter samples generated through a MCMC process.
"""
abstract type AbstractChains end

"""
    AbstractSampler

The `AbstractSampler` type is intended to be inherited from when
implementing a custom sampler. Any persistent state information should be
saved in a subtype of `AbstractSampler`.

When defining a new sampler, you should also overload the function
`transition_type`, which tells the `sample` function what type of parameter
it should expect to receive.
"""
abstract type AbstractSampler end

"""
    AbstractModel

An `AbstractModel` represents a generic model type that can be used to perform inference.
"""
abstract type AbstractModel end

"""
    AbstractMCMCEnsemble

An `AbstractMCMCEnsemble` algorithm represents a specific algorithm for sampling MCMC chains
in parallel.
"""
abstract type AbstractMCMCEnsemble end

"""
    MCMCThreads

The `MCMCThreads` algorithm allows users to sample MCMC chains in parallel using multiple
threads.
"""
struct MCMCThreads <: AbstractMCMCEnsemble end

"""
    MCMCDistributed

The `MCMCDistributed` algorithm allows users to sample MCMC chains in parallel using multiple
processes.
"""
struct MCMCDistributed <: AbstractMCMCEnsemble end

"""
    MCMCSerial

The `MCMCSerial` algorithm allows users to sample serially, with no thread or process parallelism.
"""
struct MCMCSerial <: AbstractMCMCEnsemble end

include("samplingstats.jl")
include("logging.jl")
include("interface.jl")
include("sample.jl")
include("stepper.jl")
include("transducer.jl")
include("deprecations.jl")

end # module AbstractMCMC


text:
# Error Codes
const CL_SUCCESS                                      =  0
const CL_DEVICE_NOT_FOUND                             = -1
const CL_DEVICE_NOT_AVAILABLE                         = -2
const CL_COMPILER_NOT_AVAILABLE                       = -3
const CL_MEM_OBJECT_ALLOCATION_FAILURE                = -4
const CL_OUT_OF_RESOURCES                             = -5
const CL_OUT_OF_HOST_MEMORY                           = -6
const CL_PROFILING_INFO_NOT_AVAILABLE                 = -7
const CL_MEM_COPY_OVERLAP                             = -8
const CL_IMAGE_FORMAT_MISMATCH                        = -9
const CL_IMAGE_FORMAT_NOT_SUPPORTED                   = -10
const CL_BUILD_PROGRAM_FAILURE                        = -11
const CL_MAP_FAILURE                                  = -12
const CL_MISALIGNED_SUB_BUFFER_OFFSET                 = -13
const CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST    = -14
const CL_COMPILE_PROGRAM_FAILURE                      = -15
const CL_LINKER_NOT_AVAILABLE                         = -15
const CL_LINK_PROGRAM_FAILURE                         = -17
const CL_DEVICE_PARTITION_FAILED                      = -18
const CL_KERNEL_ARG_INFO_NOT_AVAILABLE                = -19

const CL_INVALID_VALUE                                = -30
const CL_INVALID_DEVICE_TYPE                          = -31
const CL_INVALID_PLATFORM                             = -32
const CL_INVALID_DEVICE                               = -33
const CL_INVALID_CONTEXT                              = -34
const CL_INVALID_QUEUE_PROPERTIES                     = -35
const CL_INVALID_COMMAND_QUEUE                        = -36
const CL_INVALID_HOST_PTR                             = -37
const CL_INVALID_MEM_OBJECT                           = -38
const CL_INVALID_IMAGE_FORMAT_DESCRIPTOR              = -39
const CL_INVALID_IMAGE_SIZE                           = -40
const CL_INVALID_SAMPLER                              = -41
const CL_INVALID_BINARY                               = -42
const CL_INVALID_BUILD_OPTIONS                        = -43
const CL_INVALID_PROGRAM                              = -44
const CL_INVALID_PROGRAM_EXECUTABLE                   = -45
const CL_INVALID_KERNEL_NAME                          = -46
const CL_INVALID_KERNEL_DEFINITION                    = -47
const CL_INVALID_KERNEL                               = -48
const CL_INVALID_ARG_INDEX                            = -49
const CL_INVALID_ARG_VALUE                            = -50
const CL_INVALID_ARG_SIZE                             = -51
const CL_INVALID_KERNEL_ARGS                          = -52
const CL_INVALID_WORK_DIMENSION                       = -53
const CL_INVALID_WORK_GROUP_SIZE                      = -54
const CL_INVALID_WORK_ITEM_SIZE                       = -55
const CL_INVALID_GLOBAL_OFFSET                        = -56
const CL_INVALID_EVENT_WAIT_LIST                      = -57
const CL_INVALID_EVENT                                = -58
const CL_INVALID_OPERATION                            = -59
const CL_INVALID_GL_OBJECT                            = -60
const CL_INVALID_BUFFER_SIZE                          = -61
const CL_INVALID_MIP_LEVEL                            = -62
const CL_INVALID_GLOBAL_WORK_SIZE                     = -63
const CL_INVALID_PROPERTY                             = -64
const CL_INVALID_IMAGE_DESCRIPTOR                     = -65
const CL_INVALID_COMPILER_OPTIONS                     = -66
const CL_INVALID_LINKER_OPTIONS                       = -67
const CL_INVALID_DEVICE_PARTITION_COUNT               = -68
const CL_INVALID_PIPE_SIZE                            = -69
const CL_INVALID_DEVICE_QUEUE                         = -70

# OpenCL Version
const CL_VERSION_1_0                                  = cl_bool(1)
const CL_VERSION_1_1                                  = cl_bool(1)
const CL_VERSION_1_2                                  = cl_bool(1)
const CL_VERSION_2_0                                  = cl_bool(1)

# cl_bool
const CL_FALSE                                        = cl_bool(0)
const CL_TRUE                                         = cl_bool(1)
const CL_BLOCKING                                     = CL_TRUE
const CL_NON_BLOCKING                                 = CL_FALSE

# cl_platform_info
const CL_PLATFORM_PROFILE                             = cl_uint(0x0900)
const CL_PLATFORM_VERSION                             = cl_uint(0x0901)
const CL_PLATFORM_NAME                                = cl_uint(0x0902)
const CL_PLATFORM_VENDOR                              = cl_uint(0x0903)
const CL_PLATFORM_EXTENSIONS                          = cl_uint(0x0904)

# cl_device_type - bitfield
const CL_DEVICE_TYPE_DEFAULT                          = cl_bitfield(1 << 0)
const CL_DEVICE_TYPE_CPU                              = cl_bitfield(1 << 1)
const CL_DEVICE_TYPE_GPU                              = cl_bitfield(1 << 2)
const CL_DEVICE_TYPE_ACCELERATOR                      = cl_bitfield(1 << 3)
const CL_DEVICE_TYPE_CUSTOM                           = cl_bitfield(1 << 4)
const CL_DEVICE_TYPE_ALL                              = cl_bitfield(0xFFFFFFFF)

# cl_device_info
const CL_DEVICE_TYPE                                  = cl_uint(0x1000)
const CL_DEVICE_VENDOR_ID                             = cl_uint(0x1001)
const CL_DEVICE_MAX_COMPUTE_UNITS                     = cl_uint(0x1002)
const CL_DEVICE_MAX_WORK_ITEM_DIMENSIONS              = cl_uint(0x1003)
const CL_DEVICE_MAX_WORK_GROUP_SIZE                   = cl_uint(0x1004)
const CL_DEVICE_MAX_WORK_ITEM_SIZES                   = cl_uint(0x1005)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_CHAR           = cl_uint(0x1006)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_SHORT          = cl_uint(0x1007)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_INT            = cl_uint(0x1008)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_LONG           = cl_uint(0x1009)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_FLOAT          = cl_uint(0x100A)
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_DOUBLE         = cl_uint(0x100B)
const CL_DEVICE_MAX_CLOCK_FREQUENCY                   = cl_uint(0x100C)
const CL_DEVICE_ADDRESS_BITS                          = cl_uint(0x100D)
const CL_DEVICE_MAX_READ_IMAGE_ARGS                   = cl_uint(0x100E)
const CL_DEVICE_MAX_WRITE_IMAGE_ARGS                  = cl_uint(0x100F)
const CL_DEVICE_MAX_MEM_ALLOC_SIZE                    = cl_uint(0x1010)
const CL_DEVICE_IMAGE2D_MAX_WIDTH                     = cl_uint(0x1011)
const CL_DEVICE_IMAGE2D_MAX_HEIGHT                    = cl_uint(0x1012)
const CL_DEVICE_IMAGE3D_MAX_WIDTH                     = cl_uint(0x1013)
const CL_DEVICE_IMAGE3D_MAX_HEIGHT                    = cl_uint(0x1014)
const CL_DEVICE_IMAGE3D_MAX_DEPTH                     = cl_uint(0x1015)
const CL_DEVICE_IMAGE_SUPPORT                         = cl_uint(0x1016)
const CL_DEVICE_MAX_PARAMETER_SIZE                    = cl_uint(0x1017)
const CL_DEVICE_MAX_SAMPLERS                          = cl_uint(0x1018)
const CL_DEVICE_MEM_BASE_ADDR_ALIGN                   = cl_uint(0x1019)
const CL_DEVICE_MIN_DATA_TYPE_ALIGN_SIZE              = cl_uint(0x101A)
const CL_DEVICE_SINGLE_FP_CONFIG                      = cl_uint(0x101B)
const CL_DEVICE_GLOBAL_MEM_CACHE_TYPE                 = cl_uint(0x101C)
const CL_DEVICE_GLOBAL_MEM_CACHELINE_SIZE             = cl_uint(0x101D)
const CL_DEVICE_GLOBAL_MEM_CACHE_SIZE                 = cl_uint(0x101E)
const CL_DEVICE_GLOBAL_MEM_SIZE                       = cl_uint(0x101F)
const CL_DEVICE_MAX_CONSTANT_BUFFER_SIZE              = cl_uint(0x1020)
const CL_DEVICE_MAX_CONSTANT_ARGS                     = cl_uint(0x1021)
const CL_DEVICE_LOCAL_MEM_TYPE                        = cl_uint(0x1022)
const CL_DEVICE_LOCAL_MEM_SIZE                        = cl_uint(0x1023)
const CL_DEVICE_ERROR_CORRECTION_SUPPORT              = cl_uint(0x1024)
const CL_DEVICE_PROFILING_TIMER_RESOLUTION            = cl_uint(0x1025)
const CL_DEVICE_ENDIAN_LITTLE                         = cl_uint(0x1026)
const CL_DEVICE_AVAILABLE                             = cl_uint(0x1027)
const CL_DEVICE_COMPILER_AVAILABLE                    = cl_uint(0x1028)
const CL_DEVICE_EXECUTION_CAPABILITIES                = cl_uint(0x1029)
const CL_DEVICE_QUEUE_PROPERTIES                      = cl_uint(0x102A) # deprecated
const CL_DEVICE_QUEUE_ON_HOST_PROPERTIES              = cl_uint(0x102A)
const CL_DEVICE_NAME                                  = cl_uint(0x102B)
const CL_DEVICE_VENDOR                                = cl_uint(0x102C)
const CL_DRIVER_VERSION                               = cl_uint(0x102D)
const CL_DEVICE_PROFILE                               = cl_uint(0x102E)
const CL_DEVICE_VERSION                               = cl_uint(0x102F)
const CL_DEVICE_EXTENSIONS                            = cl_uint(0x1030)
const CL_DEVICE_PLATFORM                              = cl_uint(0x1031)
const CL_DEVICE_DOUBLE_FP_CONFIG                      = cl_uint(0x1032)

# 0x1033 reserved for CL_DEVICE_HALF_FP_CONFIG
const CL_DEVICE_PREFERRED_VECTOR_WIDTH_HALF           = cl_uint(0x1034)
const CL_DEVICE_HOST_UNIFIED_MEMORY                   = cl_uint(0x1035) # deprecated
const CL_DEVICE_NATIVE_VECTOR_WIDTH_CHAR              = cl_uint(0x1036)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_SHORT             = cl_uint(0x1037)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_INT               = cl_uint(0x1038)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_LONG              = cl_uint(0x1039)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_FLOAT             = cl_uint(0x103A)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_DOUBLE            = cl_uint(0x103B)
const CL_DEVICE_NATIVE_VECTOR_WIDTH_HALF              = cl_uint(0x103C)
const CL_DEVICE_OPENCL_C_VERSION                      = cl_uint(0x103D)
const CL_DEVICE_LINKER_AVAILABLE                      = cl_uint(0x103E)
const CL_DEVICE_BUILT_IN_KERNELS                      = cl_uint(0x103F)
const CL_DEVICE_IMAGE_MAX_BUFFER_SIZE                 = cl_uint(0x1040)
const CL_DEVICE_IMAGE_MAX_ARRAY_SIZE                  = cl_uint(0x1041)
const CL_DEVICE_PARENT_DEVICE                         = cl_uint(0x1042)
const CL_DEVICE_PARTITION_MAX_SUB_DEVICES             = cl_uint(0x1043)
const CL_DEVICE_PARTITION_PROPERTIES                  = cl_uint(0x1044)
const CL_DEVICE_PARTITION_AFFINITY_DOMAIN             = cl_uint(0x1045)
const CL_DEVICE_PARTITION_TYPE                        = cl_uint(0x1046)
const CL_DEVICE_REFERENCE_COUNT                       = cl_uint(0x1047)
const CL_DEVICE_PREFERRED_INTEROP_USER_SYNC           = cl_uint(0x1048)
const CL_DEVICE_PRINTF_BUFFER_SIZE                    = cl_uint(0x1049)
const CL_DEVICE_IMAGE_PITCH_ALIGNMENT                 = cl_uint(0x104A)
const CL_DEVICE_IMAGE_BASE_ADDRESS_ALIGNMENT          = cl_uint(0x104B)
const CL_DEVICE_MAX_READ_WRITE_IMAGE_ARGS             = cl_uint(0x104C)
const CL_DEVICE_MAX_GLOBAL_VARIABLE_SIZE              = cl_uint(0x104D)
const CL_DEVICE_QUEUE_ON_DEVICE_PROPERTIES            = cl_uint(0x104E)
const CL_DEVICE_QUEUE_ON_DEVICE_PREFERRED_SIZE        = cl_uint(0x104F)
const CL_DEVICE_QUEUE_ON_DEVICE_MAX_SIZE              = cl_uint(0x1050)
const CL_DEVICE_MAX_ON_DEVICE_QUEUES                  = cl_uint(0x1051)
const CL_DEVICE_MAX_ON_DEVICE_EVENTS                  = cl_uint(0x1052)
const CL_DEVICE_SVM_CAPABILITIES                      = cl_uint(0x1053)
const CL_DEVICE_GLOBAL_VARIABLE_PREFERRED_TOTAL_SIZE  = cl_uint(0x1054)
const CL_DEVICE_MAX_PIPE_ARGS                         = cl_uint(0x1055)
const CL_DEVICE_PIPE_MAX_ACTIVE_RESERVATIONS          = cl_uint(0x1056)
const CL_DEVICE_PIPE_MAX_PACKET_SIZE                  = cl_uint(0x1057)
const CL_DEVICE_PREFERRED_PLATFORM_ATOMIC_ALIGNMENT   = cl_uint(0x1058)
const CL_DEVICE_PREFERRED_GLOBAL_ATOMIC_ALIGNMENT     = cl_uint(0x1059)
const CL_DEVICE_PREFERRED_LOCAL_ATOMIC_ALIGNMENT      = cl_uint(0x105A)

# cl_device_fp_config - bitfield
const CL_FP_DENORM                                    = cl_bitfield(1 << 0)
const CL_FP_INF_NAN                                   = cl_bitfield(1 << 1)
const CL_FP_ROUND_TO_NEAREST                          = cl_bitfield(1 << 2)
const CL_FP_ROUND_TO_ZERO                             = cl_bitfield(1 << 3)
const CL_FP_ROUND_TO_INF                              = cl_bitfield(1 << 4)
const CL_FP_FMA                                       = cl_bitfield(1 << 5)
const CL_FP_SOFT_FLOAT                                = cl_bitfield(1 << 6)
const CL_FP_CORRECTLY_ROUNDED_DIVIDE_SQRT             = cl_bitfield(1 << 7)

# cl_device_mem_cache_type
const CL_NONE                                         = cl_uint(0x0)
const CL_READ_ONLY_CACHE                              = cl_uint(0x1)
const CL_READ_WRITE_CACHE                             = cl_uint(0x2)

# cl_device_local_mem_type
const CL_LOCAL                                        = cl_uint(0x1)
const CL_GLOBAL                                       = cl_uint(0x2)

# cl_device_exec_capabilities - bitfield
const CL_EXEC_KERNEL                                  = cl_bitfield(1 << 0)
const CL_EXEC_NATIVE_KERNEL                           = cl_bitfield(1 << 1)

# cl_command_queue_properties - bitfield
const CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE          = cl_bitfield(1 << 0)
const CL_QUEUE_PROFILING_ENABLE                       = cl_bitfield(1 << 1)
const CL_QUEUE_ON_DEVICE                              = cl_bitfield(1 << 2)
const CL_QUEUE_ON_DEVICE_DEFAULT                      = cl_bitfield(1 << 3)

# cl_context_info
const CL_CONTEXT_REFERENCE_COUNT                      = cl_uint(0x1080)
const CL_CONTEXT_DEVICES                              = cl_uint(0x1081)
const CL_CONTEXT_PROPERTIES                           = cl_uint(0x1082)
const CL_CONTEXT_NUM_DEVICES                          = cl_uint(0x1083)

# cl_context_properties
const CL_CONTEXT_PLATFORM                             = cl_uint(0x1084)
const CL_CONTEXT_INTEROP_USER_SYNC                    = cl_uint(0x1085)

# cl_device_partition_property
const CL_DEVICE_PARTITION_EQUALLY                     = cl_uint(0x1086)
const CL_DEVICE_PARTITION_BY_COUNTS                   = cl_uint(0x1087)
const CL_DEVICE_PARTITION_BY_COUNTS_LIST_END          = cl_uint(0x0)
const CL_DEVICE_PARTITION_BY_AFFINITY_DOMAIN          = cl_uint(0x1088)

# cl_device_affinity_domain
const CL_DEVICE_AFFINITY_DOMAIN_NUMA                  = cl_bitfield(1 << 0)
const CL_DEVICE_AFFINITY_DOMAIN_L4_CACHE              = cl_bitfield(1 << 1)
const CL_DEVICE_AFFINITY_DOMAIN_L3_CACHE              = cl_bitfield(1 << 2)
const CL_DEVICE_AFFINITY_DOMAIN_L2_CACHE              = cl_bitfield(1 << 3)
const CL_DEVICE_AFFINITY_DOMAIN_L1_CACHE              = cl_bitfield(1 << 4)
const CL_DEVICE_AFFINITY_DOMAIN_NEXT_PARTITIONABLE    = cl_bitfield(1 << 5)

# cl_device_svm_capabilities
const CL_DEVICE_SVM_COARSE_GRAIN_BUFFER               = cl_bitfield(1 << 0)
const CL_DEVICE_SVM_FINE_GRAIN_BUFFER                 = cl_bitfield(1 << 1)
const CL_DEVICE_SVM_FINE_GRAIN_SYSTEM                 = cl_bitfield(1 << 2)
const CL_DEVICE_SVM_ATOMICS                           = cl_bitfield(1 << 3)

# cl_command_queue_info
const CL_QUEUE_CONTEXT                                = cl_uint(0x1090)
const CL_QUEUE_DEVICE                                 = cl_uint(0x1091)
const CL_QUEUE_REFERENCE_COUNT                        = cl_uint(0x1092)
const CL_QUEUE_PROPERTIES                             = cl_uint(0x1093)
const CL_QUEUE_SIZE                                   = cl_uint(0x1094)

# cl_mem_flags and cl_svm_mem_flags - bitfield
const CL_MEM_READ_WRITE                               = cl_bitfield(1 << 0)
const CL_MEM_WRITE_ONLY                               = cl_bitfield(1 << 1)
const CL_MEM_READ_ONLY                                = cl_bitfield(1 << 2)
const CL_MEM_USE_HOST_PTR                             = cl_bitfield(1 << 3)
const CL_MEM_ALLOC_HOST_PTR                           = cl_bitfield(1 << 4)
const CL_MEM_COPY_HOST_PTR                            = cl_bitfield(1 << 5)
# //reserved                                            = (1 << 6)
const CL_MEM_HOST_WRITE_ONLY                          = cl_bitfield(1 << 7)
const CL_MEM_HOST_READ_ONLY                           = cl_bitfield(1 << 8)
const CL_MEM_HOST_NO_ACCESS                           = cl_bitfield(1 << 9)
const CL_MEM_SVM_FINE_GRAIN_BUFFER                    = cl_bitfield(1 << 10) # used by cl_svm_mem_flags only
const CL_MEM_SVM_ATOMICS                              = cl_bitfield(1 << 11) # used by cl_svm_mem_flags only
const CL_MEM_KERNEL_READ_AND_WRITE                    = cl_bitfield(1 << 12)

# cl_mem_migration_flags - bitfield
const CL_MIGRATE_MEM_OBJECT_HOST                      = cl_bitfield(1 << 0)
const CL_MIGRATE_MEM_OBJECT_CONTENT_UNDEFINED         = cl_bitfield(1 << 1)

# cl_channel_order
const CL_R                                            = cl_uint(0x10B0)
const CL_A                                            = cl_uint(0x10B1)
const CL_RG                                           = cl_uint(0x10B2)
const CL_RA                                           = cl_uint(0x10B3)
const CL_RGB                                          = cl_uint(0x10B4)
const CL_RGBA                                         = cl_uint(0x10B5)
const CL_BGRA                                         = cl_uint(0x10B6)
const CL_ARGB                                         = cl_uint(0x10B7)
const CL_INTENSITY                                    = cl_uint(0x10B8)
const CL_LUMINANCE                                    = cl_uint(0x10B9)
const CL_Rx                                           = cl_uint(0x10BA)
const CL_RGx                                          = cl_uint(0x10BB)
const CL_RGBx                                         = cl_uint(0x10BC)
const CL_DEPTH                                        = cl_uint(0x10BD)
const CL_DEPTH_STENCIL                                = cl_uint(0x10BE)
const CL_sRGB                                         = cl_uint(0x10BF)
const CL_sRGBx                                        = cl_uint(0x10C0)
const CL_sRGBA                                        = cl_uint(0x10C1)
const CL_sBGRA                                        = cl_uint(0x10C2)
const CL_ABGR                                         = cl_uint(0x10C3)

# cl_channel_type
const CL_SNORM_INT8                                   = cl_uint(0x10D0)
const CL_SNORM_INT16                                  = cl_uint(0x10D1)
const CL_UNORM_INT8                                   = cl_uint(0x10D2)
const CL_UNORM_INT16                                  = cl_uint(0x10D3)
const CL_UNORM_SHORT_565                              = cl_uint(0x10D4)
const CL_UNORM_SHORT_555                              = cl_uint(0x10D5)
const CL_UNORM_INT_101010                             = cl_uint(0x10D6)
const CL_SIGNED_INT8                                  = cl_uint(0x10D7)
const CL_SIGNED_INT16                                 = cl_uint(0x10D8)
const CL_SIGNED_INT32                                 = cl_uint(0x10D9)
const CL_UNSIGNED_INT8                                = cl_uint(0x10DA)
const CL_UNSIGNED_INT16                               = cl_uint(0x10DB)
const CL_UNSIGNED_INT32                               = cl_uint(0x10DC)
const CL_HALF_FLOAT                                   = cl_uint(0x10DD)
const CL_FLOAT                                        = cl_uint(0x10DE)
const CL_UNORM_INT24                                  = cl_uint(0x10DF)

# cl_mem_object_type
const CL_MEM_OBJECT_BUFFER                            = cl_uint(0x10F0)
const CL_MEM_OBJECT_IMAGE2D                           = cl_uint(0x10F1)
const CL_MEM_OBJECT_IMAGE3D                           = cl_uint(0x10F2)
const CL_MEM_OBJECT_IMAGE2D_ARRAY                     = cl_uint(0x10F3)
const CL_MEM_OBJECT_IMAGE1D                           = cl_uint(0x10F4)
const CL_MEM_OBJECT_IMAGE1D_ARRAY                     = cl_uint(0x10F5)
const CL_MEM_OBJECT_IMAGE1D_BUFFER                    = cl_uint(0x10F6)
const CL_MEM_OBJECT_PIPE                              = cl_uint(0x10F7)

# cl_mem_info
const CL_MEM_TYPE                                     = cl_uint(0x1100)
const CL_MEM_FLAGS                                    = cl_uint(0x1101)
const CL_MEM_SIZE                                     = cl_uint(0x1102)
const CL_MEM_HOST_PTR                                 = cl_uint(0x1103)
const CL_MEM_MAP_COUNT                                = cl_uint(0x1104)
const CL_MEM_REFERENCE_COUNT                          = cl_uint(0x1105)
const CL_MEM_CONTEXT                                  = cl_uint(0x1106)
const CL_MEM_ASSOCIATED_MEMOBJECT                     = cl_uint(0x1107)
const CL_MEM_OFFSET                                   = cl_uint(0x1108)
const CL_MEM_USES_SVM_POINTER                         = cl_uint(0x1109)

# cl_image_info
const CL_IMAGE_FORMAT                                 = cl_uint(0x1110)
const CL_IMAGE_ELEMENT_SIZE                           = cl_uint(0x1111)
const CL_IMAGE_ROW_PITCH                              = cl_uint(0x1112)
const CL_IMAGE_SLICE_PITCH                            = cl_uint(0x1113)
const CL_IMAGE_WIDTH                                  = cl_uint(0x1114)
const CL_IMAGE_HEIGHT                                 = cl_uint(0x1115)
const CL_IMAGE_DEPTH                                  = cl_uint(0x1116)
const CL_IMAGE_ARRAY_SIZE                             = cl_uint(0x1117)
const CL_IMAGE_BUFFER                                 = cl_uint(0x1118)
const CL_IMAGE_NUM_MIP_LEVELS                         = cl_uint(0x1119)
const CL_IMAGE_NUM_SAMPLES                            = cl_uint(0x111A)

# cl_pipe_info
const CL_PIPE_PACKET_SIZE                             = cl_uint(0x1120)
const CL_PIPE_MAX_PACKETS                             = cl_uint(0x1121)

# cl_addressing_mode
const CL_ADDRESS_NONE                                 = cl_uint(0x1130)
const CL_ADDRESS_CLAMP_TO_EDGE                        = cl_uint(0x1131)
const CL_ADDRESS_CLAMP                                = cl_uint(0x1132)
const CL_ADDRESS_REPEAT                               = cl_uint(0x1133)
const CL_ADDRESS_MIRRORED_REPEAT                      = cl_uint(0x1134)

# cl_filter_mode
const CL_FILTER_NEAREST                               = cl_uint(0x1140)
const CL_FILTER_LINEAR                                = cl_uint(0x1141)

# cl_sampler_info
const CL_SAMPLER_REFERENCE_COUNT                      = cl_uint(0x1150)
const CL_SAMPLER_CONTEXT                              = cl_uint(0x1151)
const CL_SAMPLER_NORMALIZED_COORDS                    = cl_uint(0x1152)
const CL_SAMPLER_ADDRESSING_MODE                      = cl_uint(0x1153)
const CL_SAMPLER_FILTER_MODE                          = cl_uint(0x1154)
const CL_SAMPLER_MIP_FILTER_MODE                      = cl_uint(0x1155)
const CL_SAMPLER_LOD_MIN                              = cl_uint(0x1156)
const CL_SAMPLER_LOD_MAX                              = cl_uint(0x1157)

# cl_map_flags - bitfield
const CL_MAP_READ                                     = cl_bitfield(1 << 0)
const CL_MAP_WRITE                                    = cl_bitfield(1 << 1)
const CL_MAP_WRITE_INVALIDATE_REGION                  = cl_bitfield(1 << 2)

# cl_program_info
const CL_PROGRAM_REFERENCE_COUNT                      = cl_uint(0x1160)
const CL_PROGRAM_CONTEXT                              = cl_uint(0x1161)
const CL_PROGRAM_NUM_DEVICES                          = cl_uint(0x1162)
const CL_PROGRAM_DEVICES                              = cl_uint(0x1163)
const CL_PROGRAM_SOURCE                               = cl_uint(0x1164)
const CL_PROGRAM_BINARY_SIZES                         = cl_uint(0x1165)
const CL_PROGRAM_BINARIES                             = cl_uint(0x1166)
const CL_PROGRAM_NUM_KERNELS                          = cl_uint(0x1167)
const CL_PROGRAM_KERNEL_NAMES                         = cl_uint(0x1168)

# cl_program_build_info
const CL_PROGRAM_BUILD_STATUS                         = cl_uint(0x1181)
const CL_PROGRAM_BUILD_OPTIONS                        = cl_uint(0x1182)
const CL_PROGRAM_BUILD_LOG                            = cl_uint(0x1183)
const CL_PROGRAM_BINARY_TYPE                          = cl_uint(0x1184)
const CL_PROGRAM_BUILD_GLOBAL_VARIABLE_TOTAL_SIZE     = cl_uint(0x1185)

# cl_program_binary_type
const CL_PROGRAM_BINARY_TYPE_NONE                     = cl_uint(0x0)
const CL_PROGRAM_BINARY_TYPE_COMPILED_OBJECT          = cl_uint(0x1)
const CL_PROGRAM_BINARY_TYPE_LIBRARY                  = cl_uint(0x2)
const CL_PROGRAM_BINARY_TYPE_EXECUTABLE               = cl_uint(0x4)

# cl_build_status
const CL_BUILD_SUCCESS                                =  0
const CL_BUILD_NONE                                   = -1
const CL_BUILD_ERROR                                  = -2
const CL_BUILD_IN_PROGRESS                            = -3

# cl_kernel_info
const CL_KERNEL_FUNCTION_NAME                         = cl_uint(0x1190)
const CL_KERNEL_NUM_ARGS                              = cl_uint(0x1191)
const CL_KERNEL_REFERENCE_COUNT                       = cl_uint(0x1192)
const CL_KERNEL_CONTEXT                               = cl_uint(0x1193)
const CL_KERNEL_PROGRAM                               = cl_uint(0x1194)
const CL_KERNEL_ATTRIBUTES                            = cl_uint(0x1195)

# cl_kernel_arg_info
const CL_KERNEL_ARG_ADDRESS_QUALIFIER                 = cl_uint(0x1196)
const CL_KERNEL_ARG_ACCESS_QUALIFIER                  = cl_uint(0x1197)
const CL_KERNEL_ARG_TYPE_NAME                         = cl_uint(0x1198)
const CL_KERNEL_ARG_TYPE_QUALIFIER                    = cl_uint(0x1199)
const CL_KERNEL_ARG_NAME                              = cl_uint(0x119A)

# cl_kernel_arg_address_qualifier
const CL_KERNEL_ARG_ADDRESS_GLOBAL                    = cl_uint(0x119B)
const CL_KERNEL_ARG_ADDRESS_LOCAL                     = cl_uint(0x119C)
const CL_KERNEL_ARG_ADDRESS_CONSTANT                  = cl_uint(0x119D)
const CL_KERNEL_ARG_ADDRESS_PRIVATE                   = cl_uint(0x119E)

# cl_kernel_arg_access_qualifier
const CL_KERNEL_ARG_ACCESS_READ_ONLY                  = cl_uint(0x11A0)
const CL_KERNEL_ARG_ACCESS_WRITE_ONLY                 = cl_uint(0x11A1)
const CL_KERNEL_ARG_ACCESS_READ_WRITE                 = cl_uint(0x11A2)
const CL_KERNEL_ARG_ACCESS_NONE                       = cl_uint(0x11A3)

# cl_kernel_arg_type_qualifer
const CL_KERNEL_ARG_TYPE_NONE                         = cl_bitfield(0)
const CL_KERNEL_ARG_TYPE_CONST                        = cl_bitfield(1 << 0)
const CL_KERNEL_ARG_TYPE_RESTRICT                     = cl_bitfield(1 << 1)
const CL_KERNEL_ARG_TYPE_VOLATILE                     = cl_bitfield(1 << 2)

# cl_kernel_work_group_info
const CL_KERNEL_WORK_GROUP_SIZE                       = cl_uint(0x11B0)
const CL_KERNEL_COMPILE_WORK_GROUP_SIZE               = cl_uint(0x11B1)
const CL_KERNEL_LOCAL_MEM_SIZE                        = cl_uint(0x11B2)
const CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE    = cl_uint(0x11B3)
const CL_KERNEL_PRIVATE_MEM_SIZE                      = cl_uint(0x11B4)
const CL_KERNEL_GLOBAL_WORK_SIZE                      = cl_uint(0x11B5)

# cl_kernel_exec_info
const CL_KERNEL_EXEC_INFO_SVM_PTRS                    = cl_uint(0x11B6)
const CL_KERNEL_EXEC_INFO_SVM_FINE_GRAIN_SYSTEM       = cl_uint(0x11B7)

# cl_event_info
const CL_EVENT_COMMAND_QUEUE                          = cl_uint(0x11D0)
const CL_EVENT_COMMAND_TYPE                           = cl_uint(0x11D1)
const CL_EVENT_REFERENCE_COUNT                        = cl_uint(0x11D2)
const CL_EVENT_COMMAND_EXECUTION_STATUS               = cl_uint(0x11D3)
const CL_EVENT_CONTEXT                                = cl_uint(0x11D4)

# cl_command_type
const CL_COMMAND_NDRANGE_KERNEL                       = cl_uint(0x11F0)
const CL_COMMAND_TASK                                 = cl_uint(0x11F1)
const CL_COMMAND_NATIVE_KERNEL                        = cl_uint(0x11F2)
const CL_COMMAND_READ_BUFFER                          = cl_uint(0x11F3)
const CL_COMMAND_WRITE_BUFFER                         = cl_uint(0x11F4)
const CL_COMMAND_COPY_BUFFER                          = cl_uint(0x11F5)
const CL_COMMAND_READ_IMAGE                           = cl_uint(0x11F6)
const CL_COMMAND_WRITE_IMAGE                          = cl_uint(0x11F7)
const CL_COMMAND_COPY_IMAGE                           = cl_uint(0x11F8)
const CL_COMMAND_COPY_IMAGE_TO_BUFFER                 = cl_uint(0x11F9)
const CL_COMMAND_COPY_BUFFER_TO_IMAGE                 = cl_uint(0x11FA)
const CL_COMMAND_MAP_BUFFER                           = cl_uint(0x11FB)
const CL_COMMAND_MAP_IMAGE                            = cl_uint(0x11FC)
const CL_COMMAND_UNMAP_MEM_OBJECT                     = cl_uint(0x11FD)
const CL_COMMAND_MARKER                               = cl_uint(0x11FE)
const CL_COMMAND_ACQUIRE_GL_OBJECTS                   = cl_uint(0x11FF)
const CL_COMMAND_RELEASE_GL_OBJECTS                   = cl_uint(0x1200)
const CL_COMMAND_READ_BUFFER_RECT                     = cl_uint(0x1201)
const CL_COMMAND_WRITE_BUFFER_RECT                    = cl_uint(0x1202)
const CL_COMMAND_COPY_BUFFER_RECT                     = cl_uint(0x1203)
const CL_COMMAND_USER                                 = cl_uint(0x1204)
const CL_COMMAND_BARRIER                              = cl_uint(0x1205)
const CL_COMMAND_MIGRATE_MEM_OBJECTS                  = cl_uint(0x1206)
const CL_COMMAND_FILL_BUFFER                          = cl_uint(0x1207)
const CL_COMMAND_FILL_IMAGE                           = cl_uint(0x1208)
const CL_COMMAND_SVM_FREE                             = cl_uint(0x1209)
const CL_COMMAND_SVM_MEMCPY                           = cl_uint(0x120A)
const CL_COMMAND_SVM_MEMFILL                          = cl_uint(0x120B)
const CL_COMMAND_SVM_MAP                              = cl_uint(0x120C)
const CL_COMMAND_SVM_UNMAP                            = cl_uint(0x120D)

# command execution status
const CL_COMPLETE                                     = cl_uint(0x0)
const CL_RUNNING                                      = cl_uint(0x1)
const CL_SUBMITTED                                    = cl_uint(0x2)
const CL_QUEUED                                       = cl_uint(0x3)

# cl_buffer_create_type
const CL_BUFFER_CREATE_TYPE_REGION                    = cl_uint(0x1220)

# cl_profiling_info
const CL_PROFILING_COMMAND_QUEUED                     = cl_uint(0x1280)
const CL_PROFILING_COMMAND_SUBMIT                     = cl_uint(0x1281)
const CL_PROFILING_COMMAND_START                      = cl_uint(0x1282)
const CL_PROFILING_COMMAND_END                        = cl_uint(0x1283)
const CL_PROFILING_COMMAND_COMPLETE                   = cl_uint(0x1284)

# OpenCL OpenGL Constants

# cl_gl_context_info
const CL_CURRENT_DEVICE_FOR_GL_CONTEXT_KHR            = cl_uint(0x2006)
const CL_DEVICES_FOR_GL_CONTEXT_KHR                   = cl_uint(0x2007)

# Additional cl_context_properties
const CL_GL_CONTEXT_KHR                               = cl_uint(0x2008)
const CL_EGL_DISPLAY_KHR                              = cl_uint(0x2009)
const CL_GLX_DISPLAY_KHR                              = cl_uint(0x200A)
const CL_WGL_HDC_KHR                                  = cl_uint(0x200B)
const CL_CGL_SHAREGROUP_KHR                           = cl_uint(0x200C)

if Sys.isapple()
const CL_CONTEXT_PROPERTY_USE_CGL_SHAREGROUP_APPLE = cl_uint(0x10000000)
end


text:
# Linear array sonar constraints

# These types should be consolidated in some form -- still exploring for good factorization
"""
$(TYPEDEF)
"""
mutable struct reuseLBRA
  wTb::SE3
  M::Array{Float64,2}
  E::Euler
  inp::Vector{Float64}
  outp::Vector{Float64}
  rbe::Vector{Float64}
  reuseLBRA()=new()
  reuseLBRA(::Int)=new(SE3(0),Matrix{Float64}(LinearAlgebra.I, 4,4),Euler(0),ones(4),ones(4), zeros(3))
  reuseLBRA(a,b,c,d,e,f)=new(a,b,c,d,e,f)
end

"""
$(TYPEDEF)
"""
mutable struct LinearRangeBearingElevation <: IIF.AbstractRelativeMinimize
  range::Normal
  bearing::Normal
  elev::Uniform
  reuse::Vector{reuseLBRA}
  LinearRangeBearingElevation() = new()
  LinearRangeBearingElevation( r::Tuple{Float64,Float64}, b::Tuple{Float64,Float64}; elev=Uniform(-0.25133,0.25133)) = new(Normal(r...),Normal(b...),elev, reuseLBRA[reuseLBRA(0) for i in 1:Threads.nthreads()] )
end
function (cfo::CalcFactor{<:LinearRangeBearingElevation})(meas, pose, landm)
  return residualLRBE!(meas, pose, landm, cfo.factor.reuse[Threads.threadid()])
end

function getSample!(y::Array{Float64,2}, las::LinearRangeBearingElevation, idx::Int )
  y[1,idx] = rand(las.range)
  y[2,idx] = rand(las.bearing)
  y[3,idx] = rand(las.elev)
  nothing
end
function getSample( cfo::CalcFactor{<:LinearRangeBearingElevation}, N::Int=1 )
  y = zeros(3,N)
  for i in 1:N
    getSample!(y, cfo.factor, i)
  end
  return (y,)
end

# returns [Range Bearing Elevation] manifold difference between pose X ominus landmark L
function ominus(::Type{LinearRangeBearingElevation},
                X::AbstractVector{<:Real},
                L::AbstractVector{<:Real})
  # rangeBearing3(X, L)
  wTb = SE3(X[1:3], Euler(X[4:6]...))
  bTl = matrix(wTb)\[L[1:3];1.0]
  b = atan(bTl[2],	bTl[1])
  el = -atan(bTl[3], bTl[1])
  return [norm(bTl[1:3]); b; el]
end

function ominus!( reuse::reuseLBRA,
                  X::AbstractVector{<:Real},
                  L::AbstractVector{<:Real} )
  copyto!(reuse.wTb.t, X[1:3])
  reuse.E.R, reuse.E.P, reuse.E.Y = X[4], X[5], X[6]
  convert!(reuse.wTb.R, reuse.E)  # costly
  matrix!(reuse.M, reuse.wTb)
  reuse.inp[1:3] = L[1:3]
  reuse.outp[1:4] = reuse.M\reuse.inp  # bTl  # costly
  reuse.rbe[1] = norm(reuse.outp[1:3])
  reuse.rbe[2] = atan(reuse.outp[2],  reuse.outp[1])
  reuse.rbe[3] = -atan(reuse.outp[3], reuse.outp[1])
  nothing
end

# residual should equal zero when system is in balance
# measurement z is measurement vector with [range; bearing; elevation]
# variables are tuple (pose X [dim6], landmark L [dim3])
# function handle follows required parameter list
function residualLRBE!( z::AbstractVector{<:Real},
                        X::AbstractVector{<:Real},
                        L::AbstractVector{<:Real},
                        reuse::reuseLBRA )
  #
  # TODO upgrade so the - sign here is used on a manifold too, ominus(z,  ominus(tt, variables...)  )
  # TODO just switch directly to parameterized function

  ominus!(reuse, X, L)
  return  z .- reuse.rbe
end

#


text:
context("example simplices") do
    context("2d simplex in 2d") do
        s = Simplex((Vec(1,0.), Vec(0,1.), Vec(0,0.)))
        @fact (@inferred min_euclidean(Vec(0., 0.), s)) --> roughly(0)
        @fact min_euclidean(Vec(0.5, 0.),s) --> roughly(0)
        @fact min_euclidean(Vec(-1, -1.),s) --> roughly(√(2))
        @fact min_euclidean(Vec(-7, 0.5),s) --> roughly(7)
        @fact min_euclidean(Vec(1., 1.), s) --> roughly(√(0.5))
        @fact (@inferred volume(s)) --> roughly(1/2)
        # containment
        @fact contains(s, Vec(0.1, 0.)) --> true
        for v in vertices(s)
            @fact contains(s, v) --> true
        end
        @fact contains(s, Vec(1,0.1)) --> false
        @fact contains(s, Vec(1,0.1), atol=0.1) --> true
    end

    context("counterexample") do
        # There is a common false believe that in the proj_sqdist algorithm
        # one may throw away all negative weights. See for example
        # https://www.researchgate.net/publication/267801141_-_An_Algorithm_to_Compute_the_Distance_from_a_Point_to_a_Simplex
        # Here is a counterexample

        s = Simplex((Vec(-1,0.), Vec(0,0.), Vec(1,1.)))
        pt = Vec(2.,-1)
        @fact min_euclidean(pt, s) --> roughly(√(4.5))

        pt_proj, sqd = GeometryTypes.proj_sqdist(pt, s)
        @fact pt_proj ≈ Vec(0.5, 0.5)--> true
        @fact !isapprox(pt_proj, Vec(0,0.), atol=1e-1) --> true
        @fact sqd --> roughly(4.5)

    end

    context("3d simplex in 3d") do
        s = Simplex((Vec(1,0.,0), Vec(0,1.,0), Vec(0,0,1.), Vec(0,0.,0)))
        @fact (@inferred min_euclidean(Vec(0., 0.,0), s)) --> roughly(0)
        @fact min_euclidean(Vec(0.5, 0.,0),s) --> roughly(0)
        @fact min_euclidean(Vec(-1, -1.,0),s) --> roughly(√(2))
        @fact min_euclidean(Vec(-7, 0.5,0),s) --> roughly(7)
        @fact min_euclidean(Vec(1., 1.,0), s) --> roughly(√(0.5))
        @fact min_euclidean(Vec(1., 1.,1), s) --> roughly(√(3)*(2/3))
        @fact (@inferred volume(s)) --> roughly(1/6)
        @fact contains(s, Vec(1,0,0.1)) --> false
        @fact contains(s, Vec(0.1,0,0.1)) --> true
    end

    context("1d simplex in 2d") do
        s = Simplex((Vec(-1, 1.), Vec(1,1.)))
        proj(pt) = GeometryTypes.proj_sqdist(pt, s)[1]
        @fact proj(Vec( 0.,  2)) ≈ Vec( 0, 1.) --> true
        @fact proj(Vec( 0., -2)) ≈ Vec( 0, 1.) --> true
        @fact proj(Vec( 56., 2)) ≈ Vec( 1, 1.) --> true
        @fact proj(Vec(-56., 2)) ≈ Vec(-1, 1.) --> true
    end

    context("1d simplex in 3d") do
        s = Simplex((Vec(0,0.,0), Vec(1,1.,1)))
        @fact (@inferred min_euclidean(Vec(0., 0.,0), s)) --> roughly(0)
        @fact min_euclidean(Vec(.5, .5, .5), s) --> roughly(0)
        @fact min_euclidean(Vec(-1,0,0.), s) --> roughly(1)
        @fact min_euclidean(Vec(1,0,0.), s) --> roughly(√(2/3))
        @fact (@inferred volume(s)) --> roughly(√(3))

        @fact contains(s, Vec(0.2,0.2,0.2)) --> true
        @fact contains(s, Vec(0.1,0.0999,0.1)) --> false
    end
end


text:
# #===== PANEL library ====#
#
# function update_panels()
#     ccall((:update_panels, npanel), Void, () )
# end
#
# function doupdate()
#     ccall((:doupdate, npanel ), Void, () )
# end
#
# function new_panel( win::Ptr{Void} )
#     ccall((:new_panel, npanel ), Ptr{Void}, ( Ptr{Void}, ), win )
# end
#
# function top_panel( pan::Ptr{Void} )
#     ccall((:top_panel, npanel ), Int, ( Ptr{Void}, ), pan )
# end
#
# function bottom_panel( pan::Ptr{Void} )
#     ccall((:bottom_panel, npanel ), Int, ( Ptr{Void}, ), pan )
# end
#
# function move_panel( pan::Ptr{Void}, starty::Int, startx::Int )
#     ccall((:move_panel, npanel ), Ptr{Void}, ( Ptr{Void}, Int, Int ), pan, starty, startx )
# end
#
# function del_panel( panel::Ptr{Void} )
#     ccall((:del_panel, npanel ), Void, (Ptr{Void}, ), panel )
# end
#
# function show_panel( panel::Ptr{Void} )
#     ccall((:show_panel, npanel ), Int, (Ptr{Void}, ), panel )
# end
#
# function hide_panel( panel::Ptr{Void} )
#     ccall((:hide_panel ), Int, (Ptr{Void}, ), panel )
# end
#
# function panel_hidden( panel::Ptr{Void } )
#     ccall((:panel_hidden, npanel ), Int, (Ptr{Void}, ), panel ) != 0
# end
#
# function replace_panel( panel::Ptr{Void}, window::Ptr{Void} )
#     ccall((:replace_panel ), Void, (Ptr{Void}, Ptr{Void} ), panel, window )
# end
#
# function set_panel_userptr( p1::Ptr{Void}, p2::Ptr{Void} )
#     ccall((:set_panel_userptr), Void, (Ptr{Void}, Ptr{Void}), p1, p2 )
# end
#
# function wborder( win::Ptr{Void},ls::chtype,rs::chtype,ts::chtype,bs::chtype,tl::chtype,tr::chtype,bl::chtype,br::chtype)
#     ccall((:wborder, npanel), Int, (Ptr{Void}, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8), ls, rs,ts, bs, tl, tr, bl, br)
# end
#
# export update_panels, doupdate, new_panel, top_panel, bottom_panel, hide_panel, move_panel, del_panel, show_panel, panel_hidden, replace_panel, set_panel_userptr, wborder
#
#
# # # typedef struct panel
# # # {
# # #   WINDOW *win;#===== PANEL library ====#
# #
# # function update_panels()
# #     ccall((:update_panels, npanel), Void, () )
# # end
# #
# # function doupdate()
# #     ccall((:doupdate, npanel ), Void, () )
# # end
# #
# # function new_panel( win::Ptr{Void} )
# #     ccall((:new_panel, npanel ), Ptr{Void}, ( Ptr{Void}, ), win )
# # end
# #
# # function top_panel( pan::Ptr{Void} )
# #     ccall((:top_panel, npanel ), Int, ( Ptr{Void}, ), pan )
# # end
# #
# # function bottom_panel( pan::Ptr{Void} )
# #     ccall((:bottom_panel, npanel ), Int, ( Ptr{Void}, ), pan )
# # end
# #
# # function move_panel( pan::Ptr{Void}, starty::Int, startx::Int )
# #     ccall((:move_panel, npanel ), Ptr{Void}, ( Ptr{Void}, Int, Int ), pan, starty, startx )
# # end
# #
# # function del_panel( panel::Ptr{Void} )
# #     ccall((:del_panel, npanel ), Void, (Ptr{Void}, ), panel )
# # end
# #
# # function show_panel( panel::Ptr{Void} )
# #     ccall((:show_panel, npanel ), Int, (Ptr{Void}, ), panel )
# # end
# #
# # function hide_panel( panel::Ptr{Void} )
# #     ccall((:hide_panel ), Int, (Ptr{Void}, ), panel )
# # end
# #
# # function panel_hidden( panel::Ptr{Void } )
# #     ccall((:panel_hidden, npanel ), Int, (Ptr{Void}, ), panel ) != 0
# # end
# #
# # function replace_panel( panel::Ptr{Void}, window::Ptr{Void} )
# #     ccall((:replace_panel ), Void, (Ptr{Void}, Ptr{Void} ), panel, window )
# # end
# #
# # function set_panel_userptr( p1::Ptr{Void}, p2::Ptr{Void} )
# #     ccall((:set_panel_userptr), Void, (Ptr{Void}, Ptr{Void}), p1, p2 )
# # end
# #
# # function wborder( win::Ptr{Void},ls::chtype,rs::chtype,ts::chtype,bs::chtype,tl::chtype,tr::chtype,bl::chtype,br::chtype)
# #     ccall((:wborder, npanel), Int, (Ptr{Void}, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8), ls, rs,ts, bs, tl, tr, bl, br)
# # end
# #
# # export update_panels, doupdate, new_panel, top_panel, bottom_panel, hide_panel, move_panel, del_panel, show_panel, panel_hidden, replace_panel, set_panel_userptr, wborder
# # #   struct panel *below;
# # #   struct panel *above;
# # #   NCURSES_CONST void *user;
# # # } PANEL;
# # #
# # # #if	defined(__cplusplus)
# # # extern "C" {
# # # #endif
# # #
# # # extern NCURSES_EXPORT(WINDOW*) panel_window (const PANEL *);
# # # extern NCURSES_EXPORT(void)    update_panels (void);
# # # extern NCURSES_EXPORT(int)     hide_panel (PANEL *);
# # # extern NCURSES_EXPORT(int)     show_panel (PANEL *);
# # # extern NCURSES_EXPORT(int)     del_panel (PANEL *);
# # # extern NCURSES_EXPORT(int)     top_panel (PANEL *);
# # # extern NCURSES_EXPORT(int)     bottom_panel (PANEL *);
# # # extern NCURSES_EXPORT(PANEL*)  new_panel (WINDOW *);
# # # extern NCURSES_EXPORT(PANEL*)  panel_above (const PANEL *);
# # # extern NCURSES_EXPORT(PANEL*)  panel_below (const PANEL *);
# # # extern NCURSES_EXPORT(int)     set_panel_userptr (PANEL *, NCURSES_CONST void *);
# # # extern NCURSES_EXPORT(NCURSES_CONST void*) panel_userptr (const PANEL *);
# # # extern NCURSES_EXPORT(int)     move_panel (PANEL *, int, int);
# # # extern NCURSES_EXPORT(int)     replace_panel (PANEL *,WINDOW *);
# # # extern NCURSES_EXPORT(int)     panel_hidden (const PANEL *);
# # #
# # # #if NCURSES_SP_FUNCS
# # # extern NCURSES_EXPORT(PANEL *) ground_panel(SCREEN *);
# # # extern NCURSES_EXPORT(PANEL *) ceiling_panel(SCREEN *);
# # #
# # # extern NCURSES_EXPORT(void)    NCURSES_SP_NAME(update_panels) (SCREEN*);
# # # #endif
# # #
# # # #if	defined(__cplusplus)
# # # }
# # # #endif


text:
#by Lazaro Alonso
using CairoMakie
CairoMakie.activate!()
let
    x(t; v₀ = 270,θ = 60) = v₀*t*cos(θ*π/180)
    y(t; v₀ = 270,θ = 60, g = 9.8) = v₀*t*sin(θ*π/180) - g*t^2/2
    t = 0:0.5:50
    xcoord = x.(t)
    ycoord = y.(t)

    xnode = Node(xcoord[1:2])
    ynode = Node(ycoord[1:2])

    xnodePt = Node([xcoord[2]])
    ynodePt = Node([ycoord[2]])
    # first frame, initial plot
    fig, ax = lines(xnode, ynode, color = :black, linewidth = 2,
        figure = (resolution=(700,450),), axis = (xlabel = "x", ylabel = "y",))
    scatter!(xnodePt, ynodePt, color = :red, markersize = 15)
    # the animation is done by updating the nodes values
    record(fig,joinpath(@__DIR__, "output", "animLinePoint.mp4")) do io
        for i in 3:length(t)
            push!(xnode[], xcoord[i])
            push!(ynode[], ycoord[i])
            xnodePt[] = [xcoord[i]]
            ynodePt[] = [ycoord[i]]
            xnode[] = xnode[] # trigger all updates for the new frame
            autolimits!(ax)
            ax.title = "t = $(t[i])"
            ylims!(ax, -1000, 3000)
            recordframe!(io)  # record a new frame
        end
    end
end

using Pkg # HIDE
Pkg.status("CairoMakie") # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Test, LinearAlgebra

using XPORTA: write_ieq, IEQ

using BellScenario

@testset "N-(N-2)-N scenario" begin

@time for N in 4:7
    d = N-2

    filepath = "./data/quick_adjacency_decomposition/"
    filename = "$N-$d-$N"

    println("$filename scenario")

    PM = LocalSignaling(N,N,d)

    BG_seed = BellGame(cat(
        cat(ones(Int64,(1,1)), zeros(Int64, (1,N-1)), dims=2),
        cat(Matrix{Int64}(I,(N-1,N-1)), zeros(Int64,(N-1,1)), dims=2),
        dims =  1
    ), d)

    vertices = LocalPolytope.vertices(PM)
    canonical_facets = LocalPolytope.adjacency_decomposition(vertices, BG_seed, PM, max_vertices = 100)
    games = collect(keys(canonical_facets))

    pretty_print_txt(games, filepath*filename)
    write_ieq(filepath*"ieq/"*filename, convert(IEQ, games))
end

end


text:
"""

Given the index of the atom in the vector of coordinates of the solute or the solvent,
returns the type of the atom, that is, the index of this atom within the molecule
(goes from 1 to natomspermol)

"""
function itype(iatom::Int, natomspermol::Int)
  itype = iatom%natomspermol
  if itype == 0
    itype = natomspermol
  end
  return itype
end

# Calling using the structures of Solute and Solvent, to clear up the code above
itype(iatom::Int, s::Selection) = itype(iatom,s.natomspermol)


text:
function exp_map(psi::Vector, v::Vector)
    v_norm = l2_norm(v)
    expgam = cos(v_norm) .* psi + sin(v_norm) .* v / v_norm

    return expgam
end

function inv_exp_map(Psi::Vector, psi::Vector)
    dot = inner_product(Psi,psi)
    if dot > 1
        dot = 1
    end
    if dot < -1
        dot = -1
    end

    theta = acos(dot)

    if (theta < 1e-10)
        exp_inv = zeros(length(psi))
    else
        exp_inv = theta / sin(theta) .* (psi .- cos(theta) .* Psi)
    end

    return exp_inv
end

function l2_norm(psi::Vector; timet=LinRange(0,1,length(psi)))
    timet = collect(timet)
    l2norm = sqrt(trapz(timet, psi.*psi))
    return l2norm
end

function inner_product(psi1, psi2; timet=LinRange(0,1,length(psi1)))
    timet = collect(timet)
    ip = trapz(timet, psi1.*psi2)
    return ip
end

function f_exp1(g::func)
    area = f_L2norm(g)
    y = cos(area) .+ sin(area)/area*g.y
    if (area==0)
        y = ones(length(g.x))
    end
    out = func(g.x,y)
    return out
end

function f_L2norm(f::func)
    y = f.y
    x = f.x
    ind = sortperm(x)
    area = trapz(x[ind], (y[ind]).^2)
    return sqrt(area)
end

function f_psimean(x, y; e1=0.001, e2=0.3)
    rmy = mean(y,dims=2)
    tmp = func(x, rmy[:])
    result = rmy[:] / f_L2norm(tmp)
    out = func(x, result)
    return out
end

function f_phiinv(psi::func)
    f_domain = psi.x

    result = zeros(length(f_domain))
    for i in 2:length(result)
        tmp = func(f_domain[1:i], psi.y[1:i])
        result[i] = f_L2norm(tmp)^2
    end
    out = func(f_domain,result)
    return out
end


text:
using .TerminalExtensions

function putatend(idisplay, display::iTerm2.InlineDisplay)
  iterm2display = splice!(Base.Multimedia.displays, idisplay)
  push!(Base.Multimedia.displays, iterm2display)
  true
end
putatend(idisplay, display) = false

for (idisplay, display) in enumerate(Base.Multimedia.displays)
  putatend(idisplay, display) && break
end

function display(d::iTerm2.InlineDisplay, p::Union{Plot,Compose.Context})
    draw(PNG(), p)
    print("\r")
end


text:
#Fienga et al. (A&A, 2006), Eqs. (2) and (4)
function μ_star_fun(μ, q, i)
    N = (length(q)-13) ÷ 6 # total number of bodies (Sun+planets+Moon+Pluto+asts)
    sum1 = zero(q[1])
    sum2 = zero(q[1])
    for j in 1:N
        i ==j && continue
        #@show j
        ith_body_ind = nbodyind(N, i)
        jth_body_ind = nbodyind(N, j)
        rvec_i = q[ith_body_ind[1:3]]
        rvec_j = q[jth_body_ind[1:3]]
        rvec_ij = rvec_i - rvec_j
        r_ij = sqrt(rvec_ij[1]^2 + rvec_ij[2]^2 + rvec_ij[3]^2)
        sum1 += μ[j]/r_ij
        vvec_i = q[ith_body_ind[4:6]]
        vvec_j = q[jth_body_ind[4:6]]
        vvec_ipj = vvec_i + vvec_j
        sum2 -= μ[j]*( dot(rvec_ij, vvec_ipj) )/(r_ij^3)
    end
    #@show sum1, sum2
    #rvec_i = q[3i-2:3i]
    vvec_i = q[3(N+i)-2:3(N+i)]
    v_i2 = dot(vvec_i, vvec_i)
    mustar = μ[i]*( 1+v_i2/(2c_au_per_day^2)-sum1/(2c_au_per_day^2) )
    mustardot = μ[i]*sum2/(2c_au_per_day^2)
    #@show μ[i]
    return mustar, mustardot
end

#Fienga et al. (A&A, 2006), LHS of Eq. (5)
function ssb_posvel_pN(μ, q)
    N = (length(q)-13) ÷ 6 # total number of bodies (Sun+planets+Moon+Pluto+asts)
    rvec_ssb = zeros(3)
    vvec_ssb = zeros(3)
    μ_star_SSB = zero(q[1])
    for i in 1:N
        μ_star_i, μ_star_dot_i = μ_star_fun(μ, q, i)
        ith_body_ind = nbodyind(N, i)
        rvec_i = q[ith_body_ind[1:3]]
        vvec_i = q[ith_body_ind[4:6]]
        rvec_ssb += μ_star_i*rvec_i
        vvec_ssb += μ_star_i*vvec_i + μ_star_dot_i*rvec_i
        μ_star_SSB += μ_star_i
    end
    return rvec_ssb, vvec_ssb, μ_star_SSB
end

#Fienga et al. (A&A, 2006), Eq. (5), solved for Sun position and velocity
function sun_posvel_pN(μ, q)
    N = (length(q)-13) ÷ 6 # total number of bodies (Sun+planets+Moon+Pluto+asts)
    rvec_ss_wout_sun = zeros(3)
    vvec_ss_wout_sun = zeros(3)
    for i in 2:N
        μ_star_i, μ_star_dot_i = μ_star_fun(μ, q, i)
        ith_body_ind = nbodyind(N, i)
        rvec_i = q[ith_body_ind[1:3]]
        vvec_i = q[ith_body_ind[4:6]]
        rvec_ss_wout_sun += μ_star_i*rvec_i
        vvec_ss_wout_sun += μ_star_i*vvec_i + μ_star_dot_i*rvec_i
    end
    μ_star_1, μ_star_dot_1 = μ_star_fun(μ, q, 1)
    rvec_1 = q[1:3]
    vvec_ss_wout_sun += μ_star_dot_1*rvec_1
    return -rvec_ss_wout_sun/μ_star_1, -vvec_ss_wout_sun/μ_star_1
end


text:
using Nabla
using Nabla: unbox, pos, tape, oneslike, zeroslike
using BenchmarkTools
using ChainRulesCore
using Distributions
using ForwardDiff
using LinearAlgebra
using Random
using SpecialFunctions
using Statistics
using Test

# Helper function for comparing `Ref`s, since they don't compare equal under `==`
ref_equal(a::Ref{T}, b::Ref{T}) where {T} = a[] == b[]
ref_equal(a::Ref, b::Ref) = false

# for comparing against scalar rules
derivative_via_frule(f, x) = last(frule((NoTangent(), 1.0), f, x))
# Sensiblity checkes that his is defined right
@test derivative_via_frule(cos, 0) == 0
@test derivative_via_frule(sin, 0) == 1
@test derivative_via_frule(sin, 1.2) == derivative_via_frule(sin, 2π + 1.2)

# These are the core scalar sensitives Nabla expects to have defined
# we test against them both for sensitives/scalar.jl and in sensitivities/functional.jl

const UNARY_SCALAR_SENSITIVITIES = [
    # Base:
    +, -, abs, abs2, acos, acosd, acosh, acot, acotd, acoth, acsc, acscd, acsch, asec,
    asecd, asech, asin, asind, asinh, atand, atanh, cbrt, cos, cosd, cosh, cospi, cot,cotd,
    coth, csc, cscd, csch, deg2rad, exp, exp10, exp2, expm1, inv, log, log10, log2,
    rad2deg, sec, secd, sech, sin, sind, sinh, sinpi, sqrt, tan, tand, tanh, transpose,
    # SpecialFunctions.jl:
    airyai, airyaiprime, airybi, airybiprime, besselj0, besselj1, bessely0, bessely1,
    dawson, digamma, erf, erfc, erfcinv, erfcx, erfi, erfinv, gamma, invdigamma, trigamma,
]

const DIFF_IN_FIRST_ANND_SECOND_ARG_SENSITIVITIES = [ #Base
    # Base:
    *, +, -, /, \, ^, hypot, max, min,
    # SpecialFunctions.jl:
    beta,
]
const ONLY_DIFF_IN_SECOND_ARG_SENSITIVITIES = [
    # SpecialFunctions.jl:
    besseli, besselj, besselk, bessely, polygamma,
]
const BINARY_SCALAR_SENSITIVITIES = [
    DIFF_IN_FIRST_ANND_SECOND_ARG_SENSITIVITIES;
    ONLY_DIFF_IN_SECOND_ARG_SENSITIVITIES;
]

@testset "Nabla.jl" begin

@testset "Core" begin
    include("core.jl")
    include("code_transformation/util.jl")
    include("code_transformation/differentiable.jl")
    include("sensitivity.jl")
end

@testset "Sensitivities" begin
    include("finite_differencing.jl")

    include("sensitivities/ambiguity_resolution.jl")
    # Test sensitivities for the basics.
    include("sensitivities/indexing.jl")
    include("sensitivities/scalar.jl")
    include("sensitivities/array.jl")

    # Test sensitivities for functionals.
    @testset "Functional" begin
        include("sensitivities/functional/functional.jl")
        include("sensitivities/functional/reduce.jl")
        include("sensitivities/functional/reducedim.jl")
    end

    # Test sensitivities for linear algebra optimisations.
    @testset "Linear algebra" begin
        include("sensitivities/linalg/generic.jl")
        include("sensitivities/linalg/symmetric.jl")
        include("sensitivities/linalg/uniformscaling.jl")
        include("sensitivities/linalg/diagonal.jl")
        include("sensitivities/linalg/triangular.jl")
        include("sensitivities/linalg/blas.jl")

        @testset "Factorisations" begin
            include("sensitivities/linalg/factorization/cholesky.jl")
            include("sensitivities/linalg/factorization/svd.jl")
        end
    end

    include("checkpointing.jl")
end

end