text:
using JuMP; using NCDatasets; using AxisArrays; using CalliopeJuMP.Util;

function load_conversion_plus_constraints(model_dict)
    sets = model_dict["sets"]
    constraint_dict = Dict()

    constraint_dict["balance_conversion_plus_primary_constraint"] = (
        build_constraint(
            ["loc_techs_balance_conversion_plus_primary_constraint",
             "timesteps"], "balance_conversion_plus_primary", model_dict
        )
    )

    constraint_dict["carrier_production_max_conversion_plus_constraint"] = (
        build_constraint(
            ["loc_techs_carrier_production_max_conversion_plus_constraint",
             "timesteps"], "carrier_production_max_conversion_plus", model_dict
        )
    )

    if haskey(sets, "loc_techs_carrier_production_min_conversion_plus_constraint")
        constraint_dict["carrier_production_min_conversion_plus_constraint"] = (
            build_constraint(
                ["loc_techs_carrier_production_min_conversion_plus_constraint",
                "timesteps"], "carrier_production_min_conversion_plus", model_dict
            )
        )
    end

    if haskey(sets, "loc_techs_cost_var_conversion_plus_constraint")
        constraint_dict["cost_var_conversion_plus_constraint"] = (
            build_constraint(
                ["costs", "loc_techs_cost_var_conversion_plus_constraint", "timesteps"],
                "cost_var_conversion_plus", model_dict
            )
        )
    end

    if haskey(sets, "loc_techs_balance_conversion_plus_in_2_constraint")
        constraint_dict["balance_conversion_plus_in_2_constraint"] = (
            build_constraint(
                ["loc_techs_balance_conversion_plus_in_2_constraint", "timesteps"],
                "balance_conversion_plus_tiers", model_dict, tier="in_2"
            )
        )
    end

    if haskey(sets, "loc_techs_balance_conversion_plus_in_3_constraint")
        constraint_dict["balance_conversion_plus_in_3_constraint"] = (
            build_constraint(
                ["loc_techs_balance_conversion_plus_in_3_constraint", "timesteps"],
                "balance_conversion_plus_tiers", model_dict, tier="in_3"
            )
        )
    end

    if haskey(sets, "loc_techs_balance_conversion_plus_out_2_constraint")
        constraint_dict["balance_conversion_plus_out_2_constraint"] = (
            build_constraint(
                ["loc_techs_balance_conversion_plus_out_2_constraint", "timesteps"],
                "balance_conversion_plus_tiers", model_dict, tier="out_2"
            )
        )
    end

    if haskey(sets, "loc_techs_balance_conversion_plus_out_3_constraint")
        constraint_dict["balance_conversion_plus_out_3_constraint"] = (
            build_constraint(
                ["loc_techs_balance_conversion_plus_out_3_constraint", "timesteps"],
                "balance_conversion_plus_tiers", model_dict, tier="out_3"
            )
        )
    end

    return constraint_dict
end

function balance_conversion_plus_primary_constraint_rule(backend_model, set_indices, model_dict)
    """
    Balance energy carrier consumption and production for carrier_in and carrier_out
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    parameters = model_dict["parameters"]

    loc_tech_carriers_out = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, "out"]
    )
    loc_tech_carriers_in = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, "in"]
    )

    energy_eff = get_param(model_dict, "energy_eff", [loc_tech], timestep)

    carrier_prod = sum(
        variables["carrier_prod"][loc_tech_carrier, timestep] /
        get_param(model_dict, "carrier_ratios", [loc_tech_carrier, "out"])
        for loc_tech_carrier in loc_tech_carriers_out
    )
    carrier_con = sum(
        variables["carrier_con"][loc_tech_carrier, timestep] *
        get_param(model_dict, "carrier_ratios", [loc_tech_carrier, "in"])
        for loc_tech_carrier in loc_tech_carriers_in
    )

    return @constraint(backend_model, carrier_prod == -1 * carrier_con * energy_eff)
end

function carrier_production_max_conversion_plus_constraint_rule(backend_model, set_indices, model_dict)
    """
    Set maximum conversion_plus carrier production.
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    parameters = model_dict["parameters"]

    timestep_resolution = model_dict["parameters"]["timestep_resolution"][timestep]
    loc_tech_carriers_out = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, "out"]
    )

    carrier_prod = sum(variables["carrier_prod"][loc_tech_carrier, timestep]
                 for loc_tech_carrier in loc_tech_carriers_out)

    return @constraint(backend_model, carrier_prod <= timestep_resolution * variables["energy_cap"][loc_tech])
end

function carrier_production_min_conversion_plus_constraint_rule(backend_model, set_indices, model_dict)
    """
    Set minimum conversion_plus carrier production.
    """
    loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    parameters = model_dict["parameters"]

    timestep_resolution = model_dict["parameters"]["timestep_resolution"][timestep]
    min_use = get_param(model_dict, "energy_cap_min_use", [loc_tech], timestep)

    loc_tech_carriers_out = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, "out"]
    )

    carrier_prod = sum(variables["carrier_prod"][loc_tech_carrier, timestep]
                       for loc_tech_carrier in loc_tech_carriers_out)

    return @constraint(backend_model, carrier_prod >=
        timestep_resolution * variables["energy_cap"][loc_tech] * min_use
    )
end

function cost_var_conversion_plus_constraint_rule(backend_model, set_indices, model_dict)
    """
    Add time-varying conversion_plus technology costs
    """
    cost, loc_tech, timestep = set_indices
    variables = model_dict["variables"]
    parameters = model_dict["parameters"]
    sets = model_dict["sets"]
    expressions = model_dict["expressions"]

    weight = parameters["timestep_weights"][timestep]

    loc_tech_carrier = parameters["lookup_primary_loc_tech_carriers"][loc_tech]

    var_cost = 0

    if loc_tech_carrier in sets["loc_tech_carriers_prod"]
        cost_om_prod = get_param(model_dict, "cost_om_prod", [loc_tech, cost], timestep)
        if cost_om_prod > 0
            var_cost += (
                cost_om_prod * weight *
                variables["carrier_prod"][loc_tech_carrier, timestep]
            )
        end
    end

    if loc_tech_carrier in sets["loc_tech_carriers_con"]
        cost_om_con = get_param(model_dict, "cost_om_con", [loc_tech, cost], timestep)
        if cost_om_con > 0
            var_cost += (
                cost_om_con * weight *
                variables["carrier_con"][loc_tech_carrier, timestep]
            )
        end
    end

    expressions["cost_var_rhs"][cost, loc_tech, timestep] = var_cost

    return @constraint(backend_model,
        variables["cost_var"][cost, loc_tech, timestep]
        == expressions["cost_var_rhs"][cost, loc_tech, timestep]
    )
end

function balance_conversion_plus_tiers_constraint_rule(backend_model, set_indices, model_dict, tier)
    """
    Force all carrier_in_2/carrier_in_3 and carrier_out_2/carrier_out_3 to follow
    carrier_in and carrier_out (respectively).
    """
    loc_tech, timestep = set_indices
    parameters = model_dict["parameters"]
    primary_tier, decision_variable = get_conversion_plus_io(model_dict, tier)

    loc_tech_carriers_1 = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, primary_tier]
    )
    loc_tech_carriers_2 = split_comma_list(
        parameters["lookup_loc_techs_conversion_plus"][loc_tech, tier]
    )

    c_1 = sum(decision_variable[loc_tech_carrier, timestep]
        / get_param(model_dict, "carrier_ratios", [loc_tech_carrier, primary_tier])
        for loc_tech_carrier in loc_tech_carriers_1)
    c_2 = sum(decision_variable[loc_tech_carrier, timestep]
        / get_param(model_dict, "carrier_ratios", [loc_tech_carrier, tier])
        for loc_tech_carrier in loc_tech_carriers_2)
    c_min = parameters["carrier_ratios_min"][loc_tech, tier]

    return @constraint(backend_model, c_1 * c_min == c_2)
end

ext: jl

max_stars_repo_name: calliope-project/CalliopeJuMP

max_issues_repo_name: calliope-project/CalliopeJuMP

max_forks_repo_name: calliope-project/CalliopeJuMP


text:
#=
A palindromic number reads the same both ways. The largest palindrome made from the product of two 2-digit numbers is 9009 = 91 × 99.

Find the largest palindrome made from the product of two 3-digit numbers.
=#

function calc()
  min = 100; max = 999
  maxPalindrome = 0
  for i in min:max, j in min:max
      current = i * j
      if string(current) == reverse(string(current))
        maxPalindrome = current > maxPalindrome ? current : maxPalindrome
      end
  end
  maxPalindrome
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Print function of the html backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Low-level function to print the table using the text backend.
function _pt_html(io, pinfo;
                  tf::HTMLTableFormat = html_default,
                  highlighters::Union{HTMLHighlighter,Tuple} = (),
                  linebreaks::Bool = false,
                  noheader::Bool = false,
                  nosubheader::Bool = false,
                  show_row_number::Bool = false,
                  standalone::Bool = true,
                   # Deprecated
                   formatter = nothing)

    @unpack_PrintInfo pinfo
    @unpack_HTMLTableFormat tf

    # Let's create a `IOBuffer` to write everything and then transfer to `io`.
    buf_io = IOBuffer()
    buf    = IOContext(buf_io)

    !noheader && num_cols != header_num_cols &&
    error("The header length must be equal to the number of columns.")

    # Additional processing necessary if the user wants to print the header.
    if !noheader
        # If the user do not want to print the sub-header but wants to print the
        # header, then just force the number of rows in header to be 1.
        if nosubheader
            # Now, `header` will be a view of the first line of the matrix that
            # has the header.
            header = @view header[1:header_num_rows:end]
            header_num_rows = 1
        end
    end

    # Make sure that `highlighters` is always a tuple.
    !(highlighters isa Tuple) && (highlighters = (highlighters,))

    # Get the string which is printed when `print` is called in each element of
    # the matrix. Notice that we must create only the matrix with the printed
    # rows and columns.
    header_str = Matrix{String}(undef, header_num_rows, num_printed_cols)
    data_str   = Matrix{AbstractString}(undef, num_printed_rows, num_printed_cols)

    @inbounds for i = 1:num_printed_cols
        # Index of the i-th printed column in `data`.
        ic = id_cols[i]

        if !noheader
            for j = 1:header_num_rows
                header_str[j,i] = _str_escaped(sprint(print, header[(ic-1)*header_num_rows + j]))
            end
        end

        for j = 1:num_printed_rows
            # Index of the j-th printed row in `data`.
            jr = id_rows[j]

            # Apply the formatters.
            data_ij = data[jr,ic]

            for f in formatters
                data_ij = f(data_ij, jr, ic)
            end

            # Handle `nothing` and `missing`.
            if ismissing(data_ij)
                data_str_ij = "missing"
            elseif data_ij == nothing
                data_str_ij = "nothing"
            else
                data_str_ij = sprint(print, data_ij)
            end

            # If `linebreaks` is true, then replace `\n` to `<BR>`.
            if linebreaks
                data_str_ij = replace(data_str_ij, "\n" => "
")
            end

            data_str_ij_esc = _str_escaped(data_str_ij)
            data_str[j,i]   = data_str_ij_esc
        end
    end

    # Print HTML header
    # ==========================================================================

    if standalone
        println(buf, """
                <!DOCTYPE html>

                <meta charset=\"UTF-8\">

                """)
    end

    println(buf, "
")

    # Data header
    # ==========================================================================

    # Header and sub-header texts
    # --------------------------------------------------------------------------

    if !noheader
        @inbounds @views for i = 1:header_num_rows
            if (i == 1) && (header_num_rows == 1)
                println(buf, "
")
            elseif i == 1
                println(buf, "
")
            elseif i == header_num_rows
                println(buf, "
")
            else
                println(buf, "
")
            end

            # The text "Row" must appear only on the first line.
            if show_row_number
                if i == 1
                    println(buf, "Row	")
                else
                    println(buf, "	")
                end
            end

            for j = 1:num_printed_cols
                # Index of the j-th printed column in `data`.
                jc = id_cols[j]

                # Alignment of this cell.
                style = Dict{String,String}("text-align" => _html_alignment[alignment[j]])

                println(buf, _styled_html("th", header_str[i,j], style))
            end

            println(buf, "
")
        end
    end

    # Data
    # ==========================================================================

    @inbounds @views for i = 1:num_printed_rows
        ir = id_rows[i]

        println(buf, "
")

        if show_row_number
            println(buf, "" * string(ir) * "	")
        end

        for j = 1:num_printed_cols
            jc = id_cols[j]

            # Check the alignment of this cell.
            alignment_ij = alignment[jc]

            for f in cell_alignment
                aux = f(data, ir, jc)

                if aux ∈ [:l, :c, :r, :L, :C, :R]
                    alignment_ij = aux
                    break
                end
            end

            # Alignment of this cell.
            style = Dict{String,String}("text-align" => _html_alignment[alignment_ij])

            # If we have highlighters defined, then we need to verify if this
            # data should be highlight.
            for h in highlighters
                if h.f(data, ir, jc)
                    merge!(style, Dict(h.fd(h,data,i,j)))
                    break
                end
            end

            println(buf, _styled_html("td", data_str[i,j], style))
        end

        println(buf, "
")
    end

    # Print HTML footer
    # ==========================================================================

    println(buf, "
")
    if standalone
        println(buf, """

                """)
    end

    # Print the buffer into the io.
    # ==========================================================================

    if io == stdout
        display("text/html", HTML(String(take!(buf_io))))
    else
        print(io, String(take!(buf_io)))
    end

    return nothing
end


text:
export @E_cmd, Experiment

import Parsers

macro E_cmd(s)
    Experiment(s)
end

function try_parse(s, TS=(Bool, Int, Float32, Float64))
    if s == "nothing"
        nothing
    else
        for T in TS
            res = Parsers.tryparse(T, s)
            if !isnothing(res)
                return res
            end
        end
        s
    end
end

function try_parse_kw(s)
    kw = []
    # !!! obviously, it's not correct when a value is string and contains ","
    for part in split(s, ",")
        kv = split(part, "=")
        @assert length(kv) == 2
        k, v = kv
        push!(kw, Symbol(strip(k)) => try_parse(strip(v)))
    end
    NamedTuple(kw)
end

struct Experiment
    policy::Any
    env::Any
    stop_condition::Any
    hook::Any
end

function Experiment(s::String)
    m = match(r"(?<source>\w+)_(?<method>\w+)_(?<env>\w+)(\((?<game>.*)\))?", s)
    isnothing(m) && throw(
        ArgumentError(
            "invalid format, got $s, expected format is JuliaRL_DQN_Atari(game=\"pong\")`",
        ),
    )
    source = m[:source]
    method = m[:method]
    env = m[:env]
    kw_args = isnothing(m[:game]) ? (;) : try_parse_kw(m[:game])
    Experiment(Val(Symbol(source)), Val(Symbol(method)), Val(Symbol(env)); kw_args...)
end

Base.run(ex::Experiment) = run(ex.policy, ex.env, ex.stop_condition, ex.hook)

function Base.run(
    policy::AbstractPolicy,
    env::AbstractEnv,
    stop_condition=StopAfterEpisode(1),
    hook=EmptyHook(),
)
    policy, env = check(policy, env)
    _run(policy, env, stop_condition, hook)
end

"Inject some customized checkings here by overwriting this function"
check(policy, env) = policy, env

function _run(policy::AbstractPolicy, env::AbstractEnv, stop_condition, hook)

    hook(PreExperimentStage(), policy, env)
    policy(PreExperimentStage(), env)
    is_stop = false
    while !is_stop
        reset!(env)
        policy(PreEpisodeStage(), env)
        hook(PreEpisodeStage(), policy, env)

        while !is_terminated(env) # one episode
            policy(PreActStage(), env)
            hook(PreActStage(), policy, env)

            env |> policy |> env
            optimise!(policy)

            policy(PostActStage(), env)
            hook(PostActStage(), policy, env)

            if stop_condition(policy, env)
                is_stop = true
                policy(PreActStage(), env)
                hook(PreActStage(), policy, env)
                policy(env)  # let the policy see the last observation
                break
            end
        end # end of an episode

        if is_terminated(env)
            policy(PostEpisodeStage(), env)  # let the policy see the last observation
            hook(PostEpisodeStage(), policy, env)
        end
    end
    policy(PostExperimentStage(), env)
    hook(PostExperimentStage(), policy, env)
    hook
end


text:
using POMDPGym, POMDPs, Test
mdp = ContinuumWorldMDP()

render(mdp)

for i=1:1000
    @test abs(reward(mdp, rand(initialstate(mdp)))) < 1
end
s = rand(initialstate(mdp))
svec = convert_s(AbstractArray, s, mdp)
@test all(s.==svec)
@test size(svec) == (4,)

s, r = gen(mdp, s, (0.01, 0.01))
c = render(mdp, POMDPGym.Vec4(0.9, 0.9, 0.05f0, 0.05f0), POMDPGym.Vec2(0.05, 0.05))

mdp_cost = ContinuumWorldMDP(rewards = Dict(Circle(Vec2(-0.5f0,0.5f0), 0.2f0)=>1f0), costs=Dict(Circle(Vec2(-0.2f0,-0.4f0), 0.1)=>1f0))

@test POMDPGym.cost(mdp_cost, POMDPGym.Vec4(0.9, 0.9, 0.05f0, 0.05f0)) == 0
@test POMDPGym.cost(mdp_cost, POMDPGym.Vec4(-0.2, -0.4, 0.05f0, 0.05f0)) == 1

render(mdp_cost, POMDPGym.Vec4(0.9, 0.9, 0.05f0, 0.05f0))
render(mdp_cost, POMDPGym.Vec4(-0.2, -0.4, 0.05f0, 0.05f0))

max_forks_repo_name: ancorso/POMDPGym


text:
using TypeDBClient: delete_database, create_database, contains_database, define, commit, transaction, CoreSession, Proto

function delete_db_create_db(client::TypeDBClient.AbstractCoreClient, database_name::String)
    if contains_database(client, database_name)
        delete_database(client, database_name)
    end
    create_database(client, database_name)

    sess = CoreSession(client, database_name , Proto.Session_Type.SCHEMA, request_timeout=Inf)
    trans = transaction(sess, Proto.Transaction_Type.WRITE)

    define_string = read(open("typedb_defines/disease-netowrk-schema.gql"), String)

    define(trans, define_string)

    commit(trans)
    close(sess)
end


text:
function readtile(f)
    M = 8
    title = readline(f)
    id = parse(Int, split(title, " ")[2][1:4])
    mat = Array{Char,2}(undef, M+2, M+2)
    for i = 1:M+2
        line = readline(f)
        mat[i,:] = collect(line)
    end
    readline(f)
    (id, mat)
end

function edges(mat)
    edgelist = Array{Array{Char,1},1}()
    push!(edgelist, mat[1,:])
    push!(edgelist, mat[end,:])
    push!(edgelist, mat[:,1])
    push!(edgelist, mat[:,end])
    edgelist
end

function rotate(mat)
    m = copy(mat)
    m = permutedims(m,(2,1))
    for i=1:size(m,1)
        m[i,:] = reverse(m[i,:])
    end
    m
end

function flip(mat)
    m = copy(mat)
    for i=1:size(m,1)
        m[i,:] = reverse(m[i,:])
    end
    m
end

# prints solution of both parts
function x()
	f = open("20.txt")
    N = 12
    M = 8
    list = []
    cnt = Dict{Int32, Int64}()
    tiles = Dict{Int32, Array{Char,2}}()
	for i=1:N*N
		(id, tile) = readtile(f)
        tiles[id] = tile
        cnt[id] = 0
        edgelist = edges(tile)
        for edge in edgelist
            push!(list, (id, edge))
        end
	end
    for i=1:length(list), j=1:length(list)
        if i != j && list[i][1] != list[j][1] && (list[i][2] == list[j][2] || reverse(list[i][2]) == list[j][2])
            cnt[list[i][1]] += 1
        end
    end
    y=1
    starttile = 0
    for (k, v) in cnt
        if v == 2
            y *= k
            starttile = k
        end
    end
    println(y)

    img = Array{Char,2}(undef, M*N, M*N)
    lastpadding = undef
    paddings = Array{Array{Char,1},1}(undef, N)
    fill!(img, 0)
    for i=1:N
        for j=1:N
            id = undef
            tile = undef
            if i==1 && j==1
                id = starttile
                tile = copy(tiles[id])
                while true
                    e = edges(tile)
                    a = count(f->f==e[2], map(a->a[2], list)) + count(f->reverse(f)==e[2], map(a->a[2], list))
                    b = count(f->f==e[4], map(a->a[2], list)) + count(f->reverse(f)==e[4], map(a->a[2], list))
                    if a > 1 && b > 1
                        break
                    end
                    tile = rotate(tile)
                    if tile == tiles[id]
                        tile = flip(tile)
                    end
                end
                tile=permutedims(tile,(2,1))
            elseif i == 1
                id = undef
                tile = undef
                for (iid,t) in tiles
                    f = false
                    while true
                        e = edges(t)
                        a = count(f->f==e[1], map(a->a[2], list)) + count(f->reverse(f)==e[1], map(a->a[2], list))
                        if a == 1 && e[3][1:2+M] == lastpadding
                            f = true
                            tile = copy(t)
                            break
                        end
                        t = rotate(t)
                        if t == tiles[iid]
                            t = flip(t)
                        elseif t == flip(tiles[iid])
                            break
                        end
                    end
                    if f
                        id = iid
                        break
                    end
                end
            elseif j == 1
                id = undef
                tile = undef
                for (iid,t) in tiles
                    f = false
                    while true
                        e = edges(t)
                        a = count(f->f==e[3], map(a->a[2], list)) + count(f->reverse(f)==e[3], map(a->a[2], list))
                        if a == 1 && e[1][1:2+M] == paddings[j]
                            f = true
                            tile = copy(t)
                            break
                        end
                        t = rotate(t)
                        if t == tiles[iid]
                            t = flip(t)
                        elseif t == flip(tiles[iid])
                            break
                        end
                    end
                    if f
                        id = iid
                        break
                    end
                end
            else
                id = undef
                tile = undef
                for (iid,t) in tiles
                    f = false
                    while true
                        e = edges(t)
                        if e[1][1:2+M] == paddings[j] && e[3][1:2+M] == lastpadding
                            f = true
                            tile = copy(t)
                            break
                        end
                        t = rotate(t)
                        if t == tiles[iid]
                            t = flip(t)
                        elseif t == flip(tiles[iid])
                            break
                        end
                    end
                    if f
                        id = iid
                        break
                    end
                end
            end
            ind1 = (i-1)*M+1
            ind2 = (j-1)*M+1
            img[ind1:ind1+M-1, ind2:ind2+M-1] = tile[2:1+M, 2:1+M]
            lastpadding = tile[:,end]
            paddings[j] = tile[end,:]
            delete!(tiles, id)
        end
    end

    work = copy(img)

    pattern =
"""                  #
#    ##    ##    ###
 #  #  #  #  #  #   """
    dragon = permutedims(hcat(collect.(split(pattern, "\n"))...), (2,1))
    num = 0
    while true
        for i = 1:N*M-2, j = 1:N*M-19
            match = true
            for a = 1:3, b = 1:20
                if dragon[a,b] == '#' && work[i+a-1,j+b-1] != '#'
                    match = false
                    break
                end
            end
            if match
                num -= count(c->c=='#', dragon)
            end
        end
        work = rotate(work)
        if work == img
            work = flip(work)
        elseif work == flip(img)
            break
        end
    end
    num + count(c->c=='#', img)
end


text:
###################################################################
# Positive Poisson GLM
###################################################################

###################################################################
# Distributions.jl extensions
###################################################################
"""
    PositivePoisson(λ)

A *PositivePoisson distribution* (aka zero-truncated Poisson, ZTP) descibes the number of
independent events occurring within a unit time interval, given the average rate of occurrence `λ`
and, importantly, given that the number is not zero.

``P(X = k) = \\frac{λ^k}{k!(1-e^{-λ})} e^{-λ}, \\quad \\text{ for } k = 1,2,\\ldots.``

```julia
PositivePoisson()        # PositivePoisson distribution with rate parameter 1
PositivePoisson(lambda)       # PositivePoisson distribution with rate parameter lambda

params(d)        # Get the parameters, i.e. (λ,)
mean(d)          # Get the mean arrival rate, i.e. λ
```

External links:

* [PositivePoisson distribution on Wikipedia](https://en.wikipedia.org/wiki/Zero-truncated_Poisson_distribution)

"""
struct PositivePoisson <: DiscreteUnivariateDistribution
    λ::Float64

    PositivePoisson(λ::Real) = new(λ)
    PositivePoisson() = new(1.0)
end

import Base.minimum, Base.maximum
Distributions.@distr_support PositivePoisson 1 (d.λ == 0.0 ? 0 : Inf)

### Parameters
Distributions.params(d::PositivePoisson) = (d.λ,)
Distributions.rate(d::PositivePoisson) = d.λ
### Statistics
Distributions.mean(d::PositivePoisson) = d.λ / (1.0-exp(-d.λ))

function logfactorialapprox(n::Integer)
  x=n+1.0::Float64
  (x-0.5)*log(x) - x + 0.5*log(2π) + 1.0/(12.0*x)
end

function logfactorial(n::Integer)
  if n<2
    return 0.0
  end
  x=0.0::Float64
  for i=2:n
    x+=log(i)
  end
  x
end

"""
Evaluates log(exp(x)-1.0) even when exp(x) blows up and the answer is simply x
(solution is related to more general logsumexp trick)
"""
function logexpm1(x::T) where {T}
  x + log(one(T)-exp(-x))
end

# Note that the factorial will blow up with moderately large x
Distributions.pdf(d::PositivePoisson, x::Int) = d.λ^x / ((exp(d.λ)-1.0) * factorial(x))

# calculating the log factorial is too expensive, we use an approximation
# given here: http://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
Distributions.logpdf(d::PositivePoisson, x::Int) = x*log(d.λ) - logexpm1(d.λ) - logfactorialapprox(x)

# for comparison we define the exact logpdf for the positive poisson
logpdf_exact(d::PositivePoisson, x::Integer) = x*log(d.λ) - log(exp(d.λ)-1.0) - logfactorial(x)

# similar calculation for the regular poisson
logpdf_approx(d::Poisson, x::Integer) = x*log(d.λ) - d.λ - logfactorialapprox(x)

###################################################################
# GLM.jl extensions
###################################################################
struct LogProductLogLink <: Link end

# denoting by λ=e^η the intensity of the untruncated poisson
LAMBERTW_USE_NAN=true
λfn(μ) = μ + lambertw(-exp(-μ)*μ)

GLM.linkfun(::LogProductLogLink, μ) = log(λfn(μ))

# exp(η) / (1.0 - exp(-exp(η))) = 1.0 / (exp(-η) - exp(-η)exp(-exp(η)))
#                               = 1.0 / (exp(-η) - exp(-exp(η)-η))
function GLM.linkinv(::LogProductLogLink, η::T) where {T}
  λ = exp(η)
  λ / (one(T) - exp(-λ))
end

function GLM.mueta(::LogProductLogLink, η::T) where {T}
  λ = exp(η)
  expmλ = exp(-λ)
  μ = λ / (one(T) - expmλ)
  μ * (one(T) - μ*expmλ)
end

function GLM.inverselink(::LogProductLogLink, η::T) where {T}
    λ = exp(η)
    expmλ = exp(-λ)
    μ = λ / (one(T) - expmλ)
    dμdη = μ * (one(T) - μ*expmλ)
    μ, dμdη, dμdη # third arg is the glmvar which in this case equals the derivative
end

GLM.canonicallink(::PositivePoisson) = LogProductLogLink()
# GLM.glmvar(::PositivePoisson, ::LogProductLogLink, μ, η) = μ * (1.0 - μ*exp(-exp(η)))
# GLM.glmvar(::PositivePoisson, μ) = μ * (1.0 - μ*exp(-exp(log(λfn(μ)))))
GLM.glmvar(::PositivePoisson, μ::T) where {T} = μ * (one(T) - μ*exp(-λfn(μ)))
GLM.mustart(::PositivePoisson, y, wt) = y + oftype(y, 0.1)

function GLM.devresid(::PositivePoisson, y::T, μ::T) where {T}
  if y>one(T)
    if μ>one(T)
      λμ = λfn(μ)
      λy = λfn(y)
      2 * (y*(log(λy)-log(λμ)) - logexpm1(λy) + logexpm1(λμ))
    else
      typemax(T) # +∞
    end
  else
    if μ>one(T)
      2 * (λfn(μ) - log(μ))
    else
      zero(T)
    end
  end
end

GLM.dispersion_parameter(::PositivePoisson) = false

function GLM.loglik_obs(::PositivePoisson, y, μ, wt, ϕ)
  wt*logpdf(PositivePoisson(λfn(μ)), y)
end

# function GLM.dispersion(m::GLM.AbstractGLM, sqr::Bool=false)
#     wrkwts = m.rr.wrkwts
#     wrkresid = m.rr.wrkresid
#
#     if isa(m.rr.d, @compat Union{Binomial, Poisson, PositivePoisson})
#         return one(eltype(wrkwts))
#     end
#
#     s = zero(eltype(wrkwts))
#     @inbounds @simd for i = eachindex(wrkwts,wrkresid)
#         s += wrkwts[i]*abs2(wrkresid[i])
#     end
#     s /= df_residual(m)
#     sqr ? s : sqrt(s)
# end


text:
println("test matrix <-> vec conversions")
X = rand(3, 3)
Y = vecs(X) |> mat
@test X == Y
Y[1] = -1.0
@test X == Y


text:
using Revise
using FastGroupBy, BenchmarkTools
import FastGroupBy: fcollect, load_bits, roughhash
const N = Int(1e7);const K =100;
srand(1);
byvec = rand(["id"*dec(k,10) for k = 1:N÷K],N);
valvec = rand(N);
# @time FastGroupBy.radixsort_ntoh!(byvec)
@btime radixsort!($byvec)
@btime radixsort8!($byvec)

@time x = ntoh.(unsafe_load.(Ptr{UInt}.(pointer.(byvec))))

def(byvec) =  ntoh.(unsafe_load.(Ptr{UInt}.(pointer.(byvec))))
@time def(byvec)
@time byvec2 = grouptwo!(byvec, valvec); #42

@time FastGroupBy._contiguousby(sum, byvec2[1], byvec2[2]);

#rl is random length
byvec_rl = rand(["id"*dec(k,rand(1:10)) for k = 1:N÷K],N);

@time hh = hash.(byvec);
@time grouptwo!(hh, byvec)

@time load_bits.(byvec);
@time roughhash.(byvec);

@time x = unsafe_load(Ptr{UInt}(pointer("abc")))
@time ntoh(x)

@time hash.(byvec_rl);
@time load_bits.(byvec_rl);
@time load_bits.(byvec_rl, 8);
@time load_bits2.(byvec_rl);
@time roughhash.(byvec_rl);

# roughhash(byvec[1])
# load_bits(byvec[1])

@time unsafe_load.(Ptr{UInt}.(pointer.(byvec)))
@time load_bits.(byvec)

@time byvec2 = grouptwo!(byvec, valvec); #42
@time FastGroupBy._contiguousby(sum, byvec2[1], byvec2[2]);

@time radixsort!(byvec)

# using StatsBase
# srand(1);
# byvec = rand(["id"*dec(k,10) for k = 1:N÷K],N);
# valvec = rand(N);
# @time cm = countmap(byvec, weights(valvec)); # 41.7

# @code_warntype  grouptwo!(byvec, valvec);
\


text:
using BinaryBuilder

# Collection of sources required to build LLVM
llvm_ver = "6.0.0"
sources = [
    "http://releases.llvm.org/$(llvm_ver)/llvm-$(llvm_ver).src.tar.xz" =>
    "1ff53c915b4e761ef400b803f07261ade637b0c269d99569f18040f3dcee4408",
    "http://releases.llvm.org/$(llvm_ver)/cfe-$(llvm_ver).src.tar.xz" =>
    "e07d6dd8d9ef196cfc8e8bb131cbd6a2ed0b1caf1715f9d05b0f0eeaddb6df32",
    "http://releases.llvm.org/$(llvm_ver)/compiler-rt-$(llvm_ver).src.tar.xz" =>
    "d0cc1342cf57e9a8d52f5498da47a3b28d24ac0d39cbc92308781b3ee0cea79a",
    #"http://releases.llvm.org/$(llvm_ver)/lldb-$(llvm_ver).src.tar.xz" =>
    #"46f54c1d7adcd047d87c0179f7b6fa751614f339f4f87e60abceaa45f414d454",
    #"http://releases.llvm.org/$(llvm_ver)/libcxx-$(llvm_ver).src.tar.xz" =>
    #"70931a87bde9d358af6cb7869e7535ec6b015f7e6df64def6d2ecdd954040dd9",
    #"http://releases.llvm.org/$(llvm_ver)/libcxxabi-$(llvm_ver).src.tar.xz" =>
    #"91c6d9c5426306ce28d0627d6a4448e7d164d6a3f64b01cb1d196003b16d641b",
    #"http://releases.llvm.org/$(llvm_ver)/polly-$(llvm_ver).src.tar.xz" =>
    #"47e493a799dca35bc68ca2ceaeed27c5ca09b12241f87f7220b5f5882194f59c",
    #"http://releases.llvm.org/$(llvm_ver)/libunwind-$(llvm_ver).src.tar.xz" =>
    #"256c4ed971191bde42208386c8d39e5143fa4afd098e03bd2c140c878c63f1d6",

    # Include our LLVM patches
    "patches",
]

# Since we kind of do this LLVM setup twice, this is the shared setup start:
script_setup = raw"""
cd $WORKSPACE/srcdir/

# First, move our other projects into llvm/projects
for f in *.src; do
    # Don't symlink llvm itself into llvm/projects...
    if [[ ${f} == llvm-*.src ]]; then
        continue
    fi

    # clang lives in tools/clang and not projects/cfe
    if [[ ${f} == cfe-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/clang
    elif [[ ${f} == polly-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/polly
    else
        mv $(pwd)/${f} $(echo llvm-*.src)/projects/${f%-*}
    fi
done

# Next, boogie on down to llvm town
cd llvm-*.src

# Update config.guess/config.sub stuff
update_configure_scripts

# Apply all our patches
for f in $WORKSPACE/srcdir/llvm_patches/*.patch; do
    patch -p1 < ${f}
done
"""

# The very first thing we need to do is to build llvm-tblgen for x86_64-linux-gnu
# This is because LLVM's cross-compile setup is kind of borked, so we just
# build the tools natively ourselves, directly.  :/
script = script_setup * raw"""
# Build llvm-tblgen, clang-tblgen, and llvm-config
mkdir build && cd build
CMAKE_FLAGS="-DLLVM_TARGETS_TO_BUILD:STRING=host"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_CXX_FLAGS=-std=c++0x"
cmake .. ${CMAKE_FLAGS}
make -j${nproc} llvm-tblgen clang-tblgen llvm-config

# Copy the tblgens and llvm-config into our destination `bin` folder:
mkdir -p $prefix/bin
mv bin/llvm-tblgen $prefix/bin/
mv bin/clang-tblgen $prefix/bin/
mv bin/llvm-config $prefix/bin/
"""

# We'll do this build for x86_64-linux-gnu only, as that's the arch we're building on
platforms = [
    Linux(:x86_64),
]

# We only care about llvm-tblgen and clang-tblgen
products(prefix) = [
    ExecutableProduct(prefix, "llvm-tblgen", :llvm_tblgen)
    ExecutableProduct(prefix, "clang-tblgen", :clang_tblgen)
    ExecutableProduct(prefix, "llvm-config", :llvm_config)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarball, overriding ARGS so that the user doesn't shoot themselves in the foot,
# but only do this if we don't already have a Tblgen tarball available:
tblgen_tarball = joinpath("products", "tblgen.x86_64-linux-gnu.tar.gz")
if !isfile(tblgen_tarball)
    tblgen_ARGS = ["x86_64-linux-gnu"]
    if "--verbose" in ARGS
        push!(tblgen_ARGS, "--verbose")
    end
    product_hashes = build_tarballs(tblgen_ARGS, "tblgen", sources, script, platforms, products, dependencies)

    # Extract path information to the built tblgen tarball and its hash
    tblgen_tarball, tblgen_hash = product_hashes["x86_64-linux-gnu"]
    tblgen_tarball = joinpath("products", tblgen_tarball)
else
    info("Using pre-built tblgen tarball at $(tblgen_tarball)")
    using SHA: sha256
    tblgen_hash = open(tblgen_tarball) do f
        bytes2hex(sha256(f))
    end
end

# Take that tarball, feed it into our next build as another "source".
push!(sources, tblgen_tarball => tblgen_hash)

# Next, we will Bash recipe for building across all platforms
script = script_setup * raw"""
# This value is really useful later
LLVM_DIR=$(pwd)

# Let's do the actual build within the `build` subdirectory
mkdir build && cd build

# Accumulate these flags outside CMAKE_FLAGS,
# they will be added at the end.
CMAKE_CPP_FLAGS=""
CMAKE_CXX_FLAGS=""
CMAKE_C_FLAGS=""

# Start the cmake flags off with building for both our host arch, NVidia and AMD
CMAKE_FLAGS="-DLLVM_TARGETS_TO_BUILD:STRING=\"host\;NVPTX\;AMDGPU\""

# Also target Wasm because Javascript is the Platform Of The Future (TM)
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_EXPERIMENTAL_TARGETS_TO_BUILD:STRING=\"WebAssembly\""

# We want a relase build with no bindings
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_BUILD_TYPE=Release"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_BINDINGS_LIST=\"\" "

# Disable useless things like docs, terminfo, etc....
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_INCLUDE_DOCS=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_ENABLE_TERMINFO=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DHAVE_HISTEDIT_H=Off"
CMAKE_FLAGS="${CMAKE_FLAGS} -DHAVE_LIBEDIT=Off"

# We want a shared library
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_BUILD_LLVM_DYLIB:BOOL=ON"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_LINK_LLVM_DYLIB:BOOL=ON"

# Install things into $prefix, and make sure it knows we're cross-compiling
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_INSTALL_PREFIX=${prefix}"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_CROSSCOMPILING=True"

# Tell LLVM where our pre-built tblgen tools are
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TABLEGEN=${WORKSPACE}/srcdir/bin/llvm-tblgen"
CMAKE_FLAGS="${CMAKE_FLAGS} -DCLANG_TABLEGEN=${WORKSPACE}/srcdir/bin/clang-tblgen"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_CONFIG_PATH=${WORKSPACE}/srcdir/bin/llvm-config"

# Explicitly use our cmake toolchain file
CMAKE_FLAGS="${CMAKE_FLAGS} -DCMAKE_TOOLCHAIN_FILE=/opt/${target}/${target}.toolchain"

# Manually set the host triplet, as otherwise on some platforms it tries to guess using
# `ld -v`, which is hilariously wrong.
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_HOST_TRIPLE=${target}"

# Maximize homogeneity across platforms as this is necessary for windows and ppc64le
CMAKE_FLAGS="${CMAKE_FLAGS} -DLIBCXX_ENABLE_THREADS=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLIBCXX_ENABLE_MONOTONIC_CLOCK=OFF"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLIBCXXABI_ENABLE_THREADS=OFF"

if [[ "${target}" == *apple* ]]; then
    # On OSX, we need to override LLVM's looking around for our SDK
    CMAKE_FLAGS="${CMAKE_FLAGS} -DDARWIN_macosx_CACHED_SYSROOT:STRING=/opt/${target}/MacOSX10.10.sdk"

    # LLVM actually won't build against 10.8, so we bump ourselves up slightly to 10.9
    export MACOSX_DEPLOYMENT_TARGET=10.9
    export LDFLAGS=-mmacosx-version-min=10.9
fi

if [[ "${target}" == *apple* ]] || [[ "${target}" == *freebsd* ]]; then
    # On clang-based platforms we need to override the check for ffs because it doesn't work with `clang`.
    export ac_cv_have_decl___builtin_ffs=yes
fi

if [[ "${target}" == *mingw* ]]; then
    CMAKE_CPP_FLAGS="${CMAKE_CPP_FLAGS} -remap -D__USING_SJLJ_EXCEPTIONS__ -D__CRT__NO_INLINE"
    # Windows is case-insensitive and some dependencies take full advantage of that
    echo "BaseTsd.h basetsd.h" >> /opt/${target}/${target}/include/header.gcc

    # We don't build libc++, libc++abi or Polly on windows
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBCXXABI_BUILD=OFF"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBCXX_BUILD=OFF"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_TOOL_LIBUNWIND_BUILD=OFF"
    CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVM_POLLY_BUILD=OFF"
fi

# Build!
cmake .. ${CMAKE_FLAGS} -DCMAKE_C_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_C_FLAGS}" -DCMAKE_CXX_FLAGS="${CMAKE_CPP_FLAGS} ${CMAKE_CXX_FLAGS}"
cmake -LA
make -j${nproc} VERBOSE=1

# Install!
make install -j${nproc} VERBOSE=1
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
    BinaryProvider.Linux(:x86_64, :glibc),
    BinaryProvider.Linux(:aarch64, :glibc),
    BinaryProvider.Linux(:armv7l, :glibc),
    BinaryProvider.Linux(:powerpc64le, :glibc),
    BinaryProvider.MacOS(),
    BinaryProvider.Windows(:i686),
    BinaryProvider.Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    # libraries
    LibraryProduct(prefix, "libLLVM",  :libLLVM)
    LibraryProduct(prefix, "libLTO",   :libLTO)
    LibraryProduct(prefix, "libclang", :libclang)
    # tools
    ExecutableProduct(prefix, "llvm-config", :llvm_config)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "LLVM", sources, script, platforms, products, dependencies)

# Remove tblgen tarball as it's no longer useful, and we don't want to upload them.
rm(tblgen_tarball; force=true)

max_stars_repo_name: DokFaust/LLVMBuilder

max_issues_repo_name: DokFaust/LLVMBuilder

max_forks_repo_name: DokFaust/LLVMBuilder


text:
#!/usr/bin/env julia
using Covid19Inc, PGFPlotsX
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{fontspec}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\setmainfont{Asana-Math}")
p0531 = Covid19Inc.covid_pop(752, :i40i10o90, "2020-05-31") |> Covid19Inc.covidpl
p0601 = Covid19Inc.covid_pop(752, :i50i10i70i5o90, "2020-06-01") |> Covid19Inc.covidpl
pgfsave("../../images/752_covid19_2020-05-31.svg", p0531)
pgfsave("../../images/752_covid19_2020-06-01.svg", p0601)

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
#!/usr/bin/env julia

using SphericalFourierBesselDecompositions
SFB = SphericalFourierBesselDecompositions

@testset "Theory functions" begin
    rmin = 500.0
    rmax = 1000.0
    amodes = SFB.AnlmModes(3, 5, rmin, rmax)
    cmodes = SFB.ClnnModes(amodes, Δnmax=1)

    # realistic
    pk(k) = 1e4 * (k/1e-2)^(-3.1)
    Clnn = SFB.gen_Clnn_theory(pk, cmodes)
    CnlmNLM = SFB.Clnn2CnlmNLM(Clnn, cmodes)
    Clnn1 = SFB.sum_m_lmeqLM(CnlmNLM, cmodes)
    CnlmNLM1 = SFB.Clnn2CnlmNLM(Clnn1, cmodes)
    Clnn2 = SFB.sum_m_lmeqLM(CnlmNLM1, cmodes)
    @test Clnn ≈ Clnn1
    @test Clnn ≈ Clnn2

    # random pk
    for i=1:10
        Clnn = rand(SFB.getlnnsize(cmodes))
        CnlmNLM = SFB.Clnn2CnlmNLM(Clnn, cmodes)
        Clnn1 = SFB.sum_m_lmeqLM(CnlmNLM, cmodes)
        CnlmNLM1 = SFB.Clnn2CnlmNLM(Clnn1, cmodes)
        Clnn2 = SFB.sum_m_lmeqLM(CnlmNLM1, cmodes)
        @test Clnn ≈ Clnn1
        @test Clnn ≈ Clnn2
    end

    # random pk
    for i=1:10
        CnlmNLM = rand(Complex{Float64}, SFB.getnlmsize(amodes), SFB.getnlmsize(amodes))
        Clnn1 = SFB.sum_m_lmeqLM(CnlmNLM, cmodes)
        CnlmNLM1 = SFB.Clnn2CnlmNLM(Clnn1, cmodes)
        Clnn2 = SFB.sum_m_lmeqLM(CnlmNLM1, cmodes)
        CnlmNLM2 = SFB.Clnn2CnlmNLM(Clnn2, cmodes)
        # Note: cannot compare with CnlmNLM, because information is lost in first conversion.
        @test Clnn1 ≈ Clnn2
        @test CnlmNLM1 ≈ CnlmNLM2
    end
end

# vim: set sw=4 et sts=4 :


text:
"""
This function updates the expected utility of a project for Risk Neutral investors.
"""
function update_project_utility!(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               risk_preference::RiskNeutral,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes) where P <: Project{<: BuildPhase}

    update_project_npv!(project,
                        scenario_data,
                        market_prices,
                        rep_hour_weight,
                        iteration_year,
                        yearly_horizon,
                        queue_cost,
                        capacity_forward_years,
                        solver)

    finance_data = get_finance_data(project)

    for scenario in scenario_data
        scenario_name = get_name(scenario)
        set_scenario_utility!(finance_data, scenario_name, iteration_year, get_scenario_npv(finance_data)[scenario_name][iteration_year])
    end

    set_expected_utility!(finance_data, iteration_year, get_expected_npv(finance_data)[iteration_year])

    return

end

"""
This function updates the expected utility of a project for Risk Averse investors.
"""
function update_project_utility!(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               risk_preference::RiskAverse,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes) where P <: Project{<: BuildPhase}

    update_project_npv!(project,
                        scenario_data,
                        market_prices,
                        rep_hour_weight,
                        iteration_year,
                        yearly_horizon,
                        queue_cost,
                        capacity_forward_years,
                        solver)

    expected_utility = 0.0
    finance_data = get_finance_data(project)

    for scenario in scenario_data
        scenario_name = get_name(scenario)
        scenario_npv = get_scenario_npv(finance_data)[scenario_name][iteration_year]
        scenario_utility = get_constant(risk_preference) -
                  (get_multiplier(risk_preference) * exp(-get_risk_coefficient(risk_preference) * scenario_npv))
        set_scenario_utility!(finance_data, scenario_name, iteration_year, scenario_utility)

        expected_utility += scenario_utility * get_probability(scenario)

    end

    set_expected_utility!(finance_data, iteration_year, expected_utility)

    return

end

"""
This function does nothing if a project is retired.
"""
function update_project_utility!(project::P,
                               scenario_data::Vector{Scenario},
                               market_prices::MarketPrices,
                               risk_preference::R,
                               rep_hour_weight::Vector{Float64},
                               iteration_year::Int64,
                               yearly_horizon::Int64,
                               queue_cost::Vector{Float64},
                               capacity_forward_years::Int64,
                               solver::JuMP.MOI.OptimizerWithAttributes) where {P <: Project{Retired}, R <: RiskPreference}

end


text:
# Thin wrapper type, mainly for dispatch
struct FullBinner{T, A <: AbstractVector{T}} <: AbstractVector{T}
    x::A
end

@forward FullBinner.x (Base.length, Base.size, Base.lastindex, Base.ndims, Base.iterate,
                        Base.getindex, Base.setindex!, Base.view, Base.axes,
                        Base.push!, Base.append!, Base.resize!, Base.sizehint!, Base.empty!,
                        Base.isempty)

FullBinner(::Type{T} = Float64) where T = FullBinner(Vector{T}(undef, 0))

function _print_header(io::IO, B::FullBinner{T,A}) where {T, A}
    print(io, "FullBinner{$(T),$(A)}")
    nothing
end

function _println_body(io::IO, B::FullBinner{T,A}) where {T, A}
    n = length(B)
    println(io)
    print(io, "| Count: ", n)
    if n > 0 && T <: Number
        print(io, "\n| Mean: ", round.(mean(B), digits=5))
        # print(io, "\n| StdError: ", round.(std_error(B), digits=5))
    end
    nothing
end

# short version (shows up in arrays etc.)
Base.show(io::IO, B::FullBinner{T,A}) where {T, A} = print(io, "FullBinner{$(T),$(A)}()")
# verbose version (shows up in the REPL)
Base.show(io::IO, m::MIME"text/plain", B::FullBinner) = (_print_header(io, B); _println_body(io, B))

#####
# Calculation of error coefficient (function) R. (Ch. 3.4 in QMC book)
#####
"""
Groups datapoints in bins of fixed binsize and returns error coefficient R. (Eq. 3.20)
"""
function R_value(X::AbstractVector{T}, binsize::Int) where T<:Real
    N = length(X)
    n_bins = div(N,binsize)
    lastbs = rem(N,binsize)

    @views blockmeans = vec(mean(reshape(X[1:n_bins*binsize], (binsize,n_bins)), dims=1))
    # if lastbs != 0
    #     vcat(blockmeans, mean(X[n_bins*binsize+1:end]))
    #     n_bins += 1
    # end

    blocksigma2 = 1/(n_bins-1)*sum((blockmeans .- mean(X)).^2)
    return blocksigma2 != 0 ? binsize * blocksigma2 / var(X) : 0.0
end

_R2std_error(X, Rvalue) = sqrt(Rvalue*var(X)/length(X))

"""
Groups datapoints in bins of varying size `bs`.
Returns the used binsizes `bss`, the error coefficient function values `R(bss)` (Eq. 3.20), and
the cumulative error coefficients `<R>(bss)`. The function should feature a plateau,
i.e. `R(bs_p) ~ R(bs)` for `bs >= bs_p`. (Fig. 3.3)

Optional keyword `min_nbins`. Only bin sizes used that lead to at least `min_nbins` bins.
"""
function R_function(X::AbstractVector{T}; min_nbins=32) where T<:Real
    max_binsize = floor(Int, length(X)/min_nbins)
    binsizes = 1:max_binsize

    R = Vector{Float64}(undef, length(binsizes))
    @inbounds for bs in binsizes
        R[bs] = R_value(X, bs)
    end

    means = @views mean.([R[1:i] for i in 1:max_binsize])
    return binsizes, R, means
end


text:
include("./kernel.jl")
include("./transform.jl")
include("./incremental.jl")

#####
##### Structs
#####

@doc raw"""
"""
mutable struct SystemData
    U::Matrix{Float64} # p x p subset of locations with unisolvent F
    Fu::LU{Float64, Matrix{Float64}} # unisolvent covariates at U
    X::DataMatrix{Float64} # Observation locations
    W::DataMatrix{Float64} # Fu \ Fx
end

@doc raw"""
"""
struct KernelData{K<:AbstractCorrelation}
    k::K # Kernel Correlation function
    Ku::Matrix{Float64} # Correlation at locations U
    Kux::DataMatrix{Float64} # Cross correlation of observation locations
    rKx::IncrementalCholesky{Float64} # Reduced kernel system at locations X
end

@doc raw"""
"""
struct KernelSystem{K<:AbstractCorrelation}
    sd::SystemData
    kd::KernelData{K}
end

#####
##### Construction
#####

function sdata_structs(capacity, d, p)
    U = Matrix{Float64}(undef, d, p)
    Fu = Matrix{Float64}(undef, p, p)
    X = data_array(Matrix{Float64}(undef, d, capacity), 0)
    Fx = data_array(Matrix{Float64}(undef, p, capacity), 0)
    return U, Fu, X, Fx
end
nothing
function kdata_structs(capacity, d, p)
    Ku = Matrix{Float64}(undef, p, p)
    Kux = data_array(Matrix{Float64}(undef, p, capacity), 0)
    rKx = incremental_cholesky!(Matrix{Float64}(undef, capacity, capacity), 0)
    return Ku, Kux, rKx
end

function kernel_system!(Ku, Kux, rKx, Fu, Fx, k, U, X)
    sd = system_data!(Fu, Fx, U, X)
    return KernelSystem(sd, kernel_data!(Ku, Kux, rKx, k, sd))
end

extend(ks::KernelSystem, k) = extend(ks.sd, k)
function update!(ks::KernelSystem, k, X, Fx)
    ldiv!(ks.sd.Fu, Fx)
    update!(ks.kd, ks.sd, X, Fx)
    update!(ks.sd, k)
    return nothing
end

function system_data!(Fu, Fx, U, X)
    # TODO error checking for input dimensions
    Fu_fact = lu!(Fu)
    ldiv!(Fu_fact, array_view(Fx))
    return SystemData(U, Fu_fact, X, Fx)
end

function valid_update(sd::SystemData, k)
    return valid_update(sd.X, k) && valid_update(sd.W, k)
end

function extend(sd::SystemData, k)
    valid_update(sd, k) || throw(ErrorException("Size cannot exceed capacity"))
    return (X = extend(sd.X, k), Fx = extend(sd.W, k))
end
function update!(sd::SystemData, k)
    valid_update(sd, k) || throw(ErrorException("Size cannot exceed capacity"))
    update!(sd.X, k)
    update!(sd.W, k)
    return nothing
end
function remove!(sd::SystemData, k)
    valid_update(sd, -k) || throw(ErrorException("Size cannot be negative"))
    remove!(sd.X, k)
    remove!(sd.W, k)
    return nothing
end

function system_schur!(Kxx′, Kx′, Ku, Kux, Kux′, W, W′)
    mul!(Kxx′, W', Kux′, -1, 1)
    mul!(Kxx′, Kux', W′, -1, 1)

    mul!(Kx′, W′', Kux′, -1, 1)
    mul!(Kx′, Kux′', W′, -1, 1)

    tmp = Ku * W′
    mul!(Kxx′, W', tmp, 1, 1)
    mul!(Kx′, W′', tmp, 1, 1)
    return nothing
end

function kernel_data!(Ku, Kux, rKx, k, sd)
    # TODO error checking for input dimensions
    correlation!(Ku, k, sd.U)
    cross_correlation!(Kux, k, sd.U, sd.X)

    p, m = size(sd.W)
    Kxx′, Kx′ = extend(rKx, m)
    correlation!(Kx′, k, sd.X)
    system_schur!(Kxx′, Kx′, Ku, ones(p, 0), Kux, ones(p, 0), sd.W)
    update!(rKx, m)

    return KernelData(k, Ku, Kux, rKx)
end

function valid_update(kd::KernelData, k)
    return valid_update(kd.Kux, k) && valid_update(kd.rKx, k)
end

function extend(kd::KernelData, k)
    valid_update(kd, k) || throw(ErrorException("Size cannot exceed capacity"))
    Kxx′, Kx′ = extend(kd.rKx, k)
    return (Kux′ = extend(kd.Kux, k), Kxx′ = Kxx′, Kx′ = Kx′)
end
function update!(kd::KernelData, sd, X′, W′)
    p, k = size(W′)
    Kux′, Kxx′, Kx′ = extend(kd, k)

    cross_correlation!(Kux′, kd.k, sd.U, X′)
    cross_correlation!(Kxx′, kd.k, sd.X, X′)
    correlation!(Kx′, kd.k, X′)

    system_schur!(Kxx′, Kx′, kd.Ku, kd.Kux, Kux′, sd.W, W′)
    update!(kd.Kux, k)
    update!(kd.rKx, k)
    return nothing
end
function remove!(kd::KernelData, k)
    valid_update(kd, -k) || throw(ErrorException("Size cannot be negative"))
    remove!(kd.Kux, k)
    remove!(kd.rKx, k)
    return nothing
end

#####
##### Functions
#####

function reduced_solve!(cx, ks, yu, yx)
    cx .= yx
    mul!(cx, ks.sd.W', yu, -1, 1)
    ldiv!(ks.kd.rKx, cx)
    return nothing
end
function solve!(c, ks, y)
    p = size(ks.sd.W, 1)
    @views cu, cx = c[1:p], c[p+1:end]
    @views yu, yx = y[1:p], y[p+1:end]
    reduced_solve!(cx, ks, yu, yx)
    mul!(cu, ks.sd.W, cx)
    return cu, cx
end
function fit_tail!(d, ks, yu, cu, cx)
    d .= yu
    mul!(d, ks.kd.Ku, cu, -1, 1)
    mul!(d, ks.kd.Kux, cx, -1, 1)
    ldiv!(ks.sd.Fu', d)
    return nothing
end

#####
##### Prediction & Statistics
#####

function det(ks::KernelSystem)
    d, s = logabsdet(ks)
    return s * exp(d)
end
function logdet(ks::KernelSystem)
    d, s = logabsdet(ks)
    s < 0 || throw(DomainError(s))
    return d
end
function logabsdet(ks::KernelSystem)
    s = (-1) ^ (size(ks.sd.Fu, 1) % 2)
    return 2 * logabsdet(ks.sd.Fu)[1] + logdet(ks.kd.rKx), s
end

function predict!(ks::KernelSystem, y, X′, Fx′)
    cu, cx = solve!(similar(y), ks, y)
    d = similar(cu)
    fit_tail!(d, ks, @views y[1:length(d)], cu, cx)
end

function predict_point(kd::KernelData, y, X′, Fx′)
    c1, c2, β = solve_system(kd, y1, y2)
    k1 = Array{Float64}(undef, size(kd.K1, 1))
    k2 = Array{Float64}(undef, kd.n)
    point_correlation!(k1, k2, kd, x)
    m = dot(k1, c1) + dot(k2, c2) + dot(fx, β)
    return m
end

function qmC!(kd::KernelData, y1, y2, x, fx)
    c1, c2, β = solve_system(kd, y1, y2)
    k1 = Array{Float64}(undef, size(kd.K1, 1))
    k2 = Array{Float64}(undef, kd.n)
    point_correlation!(k1, k2, kd, x)
    m = dot(k1, c1) + dot(k2, c2) + dot(fx, β)

    q = dot(y1, c1) + dot(y2, c2)

    add_point!(kd, x, fx)
    C = kd.K2_tilde.R[kd.n, kd.n] ^ 2
    remove_point!(kd)

    return q, m, C
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Verboseprint_mpi

    struct Verbose_print
        level::Int8
        fp::Union{Nothing,IOStream}
        myid::Int64

        function Verbose_print(level;myid = 0,filename=nothing)
            if filename == nothing
                fp = nothing
            else
                if myid == 0
                    fp = open(filename,"w")
                else
                    fp = nothing
                end
            end
            return new(level,fp,myid)
        end
    end

    function println_verbose_level1(v::Verbose_print,val...)
        if v.myid == 0 && v.level >= 1
            println(val...)
            if v.fp != nothing
                println(v.fp,val...)
            end
        end
    end

    function println_verbose_level2(v::Verbose_print,val...)
        if v.myid == 0 && v.level >= 2
            println(val...)
            if v.fp != nothing
                println(v.fp,val...)
            end
        end
    end

    function println_verbose_level3(v::Verbose_print,val...)
        if v.myid == 0 && v.level >= 3
            println(val...)
            if v.fp != nothing
                println(v.fp,val...)
            end
        end
    end

end


text:
## UTILITIES FOR ENCODING AND DECODING MODEL METADATA
# (for serializing/deserializing into TOML format)

# fallback encoding:
function encode_dic(s)
    prestring = string("`", s, "`")
    # hack for objects with gensyms in their string representation:
    str = replace(prestring, '#'=>'_')
    return str
end

encode_dic(s::AbstractString) = string(s)
encode_dic(s::Symbol) = string(":", s)
encode_dic(s::Nothing) = "`nothing`"
encode_dic(v::AbstractVector) = encode_dic.(v)
function encode_dic(d::AbstractDict)
    ret = LittleDict{}()
    for (k, v) in d
        ret[encode_dic(k)] = encode_dic(v)
    end
    return ret
end

function decode_dic(s::String)
    if !isempty(s)
        if  s[1] == ':'
            return Symbol(s[2:end])
        elseif s[1] == '`'
            return eval(Meta.parse(s[2:end-1]))
        else
            return s
        end
    else
        return ""
    end
end

decode_dic(v::Vector) = decode_dic.(v)
function decode_dic(d::AbstractDict)
    ret = LittleDict()
    for (k, v) in d
        ret[decode_dic(k)] = decode_dic(v)
    end
    return ret
end

# the inverse of a multivalued dictionary is a multivalued
# dictionary:
function inverse(d::LittleDict{S,Set{T}}) where {S,T}
    dinv = LittleDict{T,Set{S}}()
    for key in keys(d)
        for val in d[key]
            if val in keys(dinv)
                push!(dinv[val], key)
            else
                dinv[val] = Set([key,])
            end
        end
    end
    return dinv
end
function inverse(d::Dict{S,Set{T}}) where {S,T}
    dinv = Dict{T,Set{S}}()
    for key in keys(d)
        for val in d[key]
            if val in keys(dinv)
                push!(dinv[val], key)
            else
                dinv[val] = Set([key,])
            end
        end
    end
    return dinv
end

function Base.isless(h1::Handle, h2::Handle)
    if isless(h1.name, h2.name)
        return true
    elseif h1.name == h2.name
        return isless(h1.pkg, h2.pkg)
    else
        return false
    end
end

## FUNCTIONS TO BUILD GLOBAL METADATA CONSTANTS IN MLJMODELS
## INITIALIZATION

# to define INFO_GIVEN_HANDLE
function info_given_handle(metadata_file)
    metadata = LittleDict(TOML.parsefile(metadata_file))
    metadata_given_pkg = decode_dic(metadata)

    # build info_given_handle dictionary:
    ret = Dict{Handle}{Any}()
    packages = keys(metadata_given_pkg)
    for pkg in packages
        info_given_name = metadata_given_pkg[pkg]
        for name in keys(info_given_name)
            handle = Handle(name, pkg)
            ret[handle] = info_given_name[name]
        end
    end
    return ret

end

# for use in __init__ to define AMBIGUOUS_NAMES
function ambiguous_names(info_given_handle)
    names_with_duplicates = map(keys(info_given_handle) |> collect) do handle
        handle.name
    end
    frequency_given_name = countmap(names_with_duplicates)
    return filter(keys(frequency_given_name) |> collect) do name
        frequency_given_name[name] > 1
    end
end

# for use in __init__ to define PKGS_GIVEN_NAME
function pkgs_given_name(info_given_handle)
    handles = keys(info_given_handle) |> collect
    ret = Dict{String,Vector{String}}()
    for handle in handles
        if haskey(ret, handle.name)
           push!(ret[handle.name], handle.pkg)
        else
            ret[handle.name] =[handle.pkg, ]
        end
    end
    return ret
end

function model_names(info_given_handle)
    names_allowing_duplicates = map(keys(info_given_handle) |> collect) do handle
        handle.name
    end
    return unique(names_allowing_duplicates)
end

function (::Type{Handle})(name::String)
    if name in AMBIGUOUS_NAMES
        return Handle(name, missing)
    else
        return Handle(name, first(PKGS_GIVEN_NAME[name]))
    end
end

function model_traits_in_registry(info_given_handle)
    first_entry = info_given_handle[Handle("ConstantRegressor")]
    return keys(first_entry) |> collect
end


text:
#TODO: Refactor as summary

function csv(data::ModelSelection.ModelSelectionData, filename::String; resultnum::Int64=1)
    csv(data, filename=filename, resultnum=resultnum)
end

function csv(data::ModelSelection.ModelSelectionData; filename::Union{Nothing, String}=nothing, resultnum::Int64=1)
    addextras(data, :csv, filename, nothing)
    if size(data.results, 1) > 0
        return csv(data, data.results[resultnum], filename=filename)
    end
    return ""
end

function csv(data::ModelSelection.ModelSelectionData, result::ModelSelection.AllSubsetRegression.AllSubsetRegressionResult, filename::String)
    return csv(data, result, filename=filename)
end

function csv(data::ModelSelection.ModelSelectionData, result::ModelSelection.AllSubsetRegression.AllSubsetRegressionResult; filename::Union{Nothing, String}=nothing)
    header = []
    for dataname in result.datanames
        push!(header, String(dataname))
    end

    rows = vcat(permutedims(header), result.data)

    if filename != nothing
        file = open(filename, "w")
        writedlm(file, rows, ',')
        close(file)
    end

    res = ""
    for row in eachrow(rows)
        res *= @sprintf("%s\n", join(row, ','))
    end
    return res
end

function csv(data::ModelSelection.ModelSelectionData, result::ModelSelection.CrossValidation.CrossValidationResult, filename::String)
    return csv(data, result, filename=filename)
end

function csv(data::ModelSelection.ModelSelectionData, result::ModelSelection.CrossValidation.CrossValidationResult; filename::Union{Nothing, String}=nothing)
    header = []
    for dataname in result.datanames
        push!(header, String(dataname))
    end

    rows = vcat(permutedims(header), result.data)

    if filename != nothing
        file = open(filename, "w")
        writedlm(file, rows, ',')
        close(file)
    end

    res = ""
    for row in eachrow(rows)
        res *= @sprintf("%s\n", join(row, ','))
    end
    return res
end


text:
module NetworkProtocols

using Printf
using Sockets
using UnsafeArrays

export MACAddress
export isunicast, ismulticast, isuniversal, islocal

export EthernetHeader, EthernetPacket
export ETHERTYPE_IPV4, ETHERTYPE_ARP
export decode_ethernet, ethertype_string

export IPv4Header, IPv4Packet
export IPPROTOCOL_IGMP, IPPROTOCOL_TCP, IPPROTOCOL_UDP
export decode_ipv4

export UDPHeader, UDPPacket
export decode_udp

export dispatch_ethernet

include("mac_address.jl")
include("ethernet.jl")
include("ip.jl")
include("udp.jl")
include("dispatch.jl")

end


text:
inputs = SNN.Poisson(;N = 1000)
inputs.param = SNN.PoissonParameter(;rate = 15Hz)

neurons = SNN.IF2(;N = 1)
neurons.param = SNN.IF2Parameter(;τm = 10ms, τe = 5ms, El = -74mV, Ee = 0mV, Vt = -54mV, Vr = -60mV)

S = SNN.SpikingSynapse(inputs, neurons, :ge; σ = 0.01, p = 1.0)
S.param = SNN.SpikingSynapseParameter(;Wmax = 0.01)

P = [inputs, neurons]; C = [S]

# histogram(S.W / S.param.Wmax; nbins = 20) |> gui
# SNN.monitor(S, [(:W, [1, 2])])
@time SNN.train!(P, C; duration = 100second)

#scatter(S.W / S.param.Wmax) |> gui
#histogram(S.W / S.param.Wmax; nbins = 20) |> gui
# plot(hcat(SNN.getrecord(S, :W)...)' / S.param.Wmax) |> gui
# heatmap(full(sparse(S.I, S.J, S.W / S.param.Wmax))) |> gui


text:
export iauLdn
"""
For a star, apply light deflection by multiple solar-system bodies,
as part of transforming coordinate direction into natural direction.

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

Status:  support function.

Given:
   n    int           number of bodies (note 1)
   b    iauLDBODY[n]  data for each of the n bodies (Notes 1,2):
   bm   double         mass of the body (solar masses, Note 3)
   dl   double         deflection limiter (Note 4)
   pv   [2][3]         barycentric PV of the body (au, au/day)
   ob   double[3]     barycentric position of the observer (au)
   sc   double[3]     observer to star coord direction (unit vector)

Returned:
   sn    double[3]      observer to deflected star (unit vector)

   1. The array b contains n entries, one for each body to be
   considered.  If n = 0, no gravitational light deflection will be
   applied, not even for the Sun.

   2. The array b should include an entry for the Sun as well as for
   any planet or other body to be taken into account.  The entries
   should be in the order in which the light passes the body.

   3. In the entry in the b array for body i, the mass parameter
   b[i].bm can, as required, be adjusted in order to allow for such
   effects as quadrupole field.

   4. The deflection limiter parameter b[i].dl is phi^2/2, where phi is
   the angular separation (in radians) between star and body at
   which limiting is applied.  As phi shrinks below the chosen
   threshold, the deflection is artificially reduced, reaching zero
   for phi = 0.   Example values suitable for a terrestrial
   observer, together with masses, are as follows:

      body i     b[i].bm        b[i].dl

      Sun        1.0            6e-6
      Jupiter    0.00095435     3e-9
      Saturn     0.00028574     3e-10

   5. For cases where the starlight passes the body before reaching the
   observer, the body is placed back along its barycentric track by
   the light time from that point to the observer.  For cases where
   the body is "behind" the observer no such shift is applied.  If
   a different treatment is preferred, the user has the option of
   instead using the iauLd function.  Similarly, iauLd can be used
   for cases where the source is nearby, not a star.

   6. The returned vector sn is not normalized, but the consequential
   departure from unit magnitude is always negligible.

   7. The arguments sc and sn can be the same array.

   8. For efficiency, validation is omitted.  The supplied masses must
   be greater than zero, the position and velocity vectors must be
   right, and the deflection limiter greater than zero.

Reference:

   Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
   the Astronomical Almanac, 3rd ed., University Science Books
   (2013), Section 7.2.4.

Called:
   iauCp        copy p-vector
   iauPdp       scalar product of two p-vectors
   iauPmp       p-vector minus p-vector
   iauPpsp      p-vector plus scaled p-vector
   iauPn        decompose p-vector into modulus and direction
   iauLd        light deflection by a solar-system body

This revision:   2017 March 16

SOFA release 2018-01-30

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

# void iauLdn(int n, iauLDBODY b[], double ob[3], double sc[3],
#             double sn[3])

function iauLdn(n::Real, b::Array{iauLDBODY, 1},
               ob::AbstractVector{<:Real}, sc::AbstractVector{<:Real})
   ref_b = Ref{iauLDBODY}(b[1])
   sn    = zeros(Float64, 3)

   ccall((:iauLdn, libsofa_c), Cvoid,
         (Cint, Ref{iauLDBODY}, Ptr{Cdouble}, Ptr{Cdouble},
          Ptr{Cdouble}),
          convert(Int32, n),
          ref_b,
          convert(Array{Float64, 1}, ob),  convert(Array{Float64, 1}, sc),
          sn)

   return SVector{3}(sn)
end


text:
Δx = 0.01
xlim = (-2,2)
ylim = (-2,2)
gr = PhysicalGrid(xlim,ylim,Δx)
w = Nodes(Dual,size(gr));
q = Edges(Primal,w);
p = Nodes(Primal,w);

L = plan_laplacian(w,with_inverse=true);

n = 300
θ = range(0,2π,length=n+1)
rad = 1.0
x, y = rad*cos.(θ[1:n]), rad*sin.(θ[1:n])
X = VectorData(x,y);
f = ScalarData(X);
τ = VectorData(X);

reg = Regularize(X,cellsize(gr),I0=origin(gr))

Δxg = 0.001
xg = -2.0:Δxg:2.0;

function evaluate_ddf(ddf)
  #for xi in xg, yi in xg
      ddf(xg,xg)
  #end
end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
using StemCellModels
using DataFrames
using Pkg
Pkg.build("RCall")
using RCall
using ProgressMeter
using StatsBase
using Statistics
using Random
using StatsBase
using ArgParse
R"""
library(ggplot2)
library(cowplot)
library(tidyverse)
library(viridis)
library(jcolors)
""";

# This script contains some functions to do a Maximum Likelihood fit of our model
include("optim.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--examplefitsout"
        help = "File to save example fits"
    "--exampledifferentbins"
        help = "File to save example fits"
    "--powerout"
        help = "File to save power analysis"
    "--nsamples"
        help = "Number of samples to take"
        arg_type = Int
        default = 10
end

parsed_args = parse_args(ARGS, s)

Random.seed!(12345)

#function to calculate interval dN/dS from simulations
function intervaldnds_data(Msel, Mneut, mmin, Nmax; mgap = 5.0, mup = 1.0, mud = 1.0)
    #calculate cumulative dn/ds across cohort, normalizing for mutation rates
    dndscumsum = Float64[]
    dnvec = Float64[]
    dsvec = Float64[]
    driverVAFt = filter(x -> x >= mmin, Msel)
    passengerVAFt = filter(x -> x >= mmin, Mneut)
    freqs = collect(mmin:mgap:Nmax)
    for f in freqs
        dn = length(filter(x -> x< f, driverVAFt))/mud
        ds = length(filter(x -> x< f, passengerVAFt))/mup
        push!(dndscumsum, dn/ds)
        push!(dnvec, dn)
        push!(dsvec, ds)
    end
    return dndscumsum[2:end], freqs[2:end], dnvec[2:end], dsvec[2:end]
end

#function to simulate population of cells with different paramters
function simulatepopulation(;Δ = 0.0, ρ = 100.0, Amin = 1 ./ ρ, gap = 0.02, tend = 20.0,
        Ncells = 1000, r = 0.5, λ = 1.0,
        sample = false, nsample = 10, Amax = 0.0)

    rangeN = 1:10^4
    Nits = 100
    Cntd = zeros(length(rangeN), Nits)
    Cntp = zeros(length(rangeN), Nits)
    Msel = Int64[]
    Mneut = Int64[]
    SMtog = SkinStemCellModel(Ncells, Δmut = Δ, μp = 0.001, μd = 0.001, tend = tend, r = r, λ = λ)
    for i in 1:Nits
        scst = runsimulation(SMtog, progress = false, finish = "time", onedriver = true, restart = true)
        append!(Msel, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[2]), rangeN)))
        append!(Mneut, filter!(x -> x > 0, counts(sort(StemCellModels.cellsconvert(scst.stemcells)[1]), rangeN)))
    end

    Asel = Msel ./ ρ
    Aneut = Mneut ./ ρ
    Asel = filter(x -> x >= Amin, Asel)
    Aneut = filter(x -> x >= Amin, Aneut)

    if sample
        prevlength = length(Asel)
        Asel = Asel[StatsBase.sample(1:length(Asel), nsample, replace = false)]
        Aneut = Aneut[StatsBase.sample(1:length(Aneut),
                Int64(round(nsample * (length(Aneut) / prevlength))), replace = false)]
    end

    if Amax == 0.0
        Amax = maximum(Asel)
    end

    x1 = intervaldnds_data(Asel, Aneut, Amin, Amax; mgap = gap, mup = SMtog.μp, mud = SMtog.μd)

    DFinterval = DataFrame(dnds = x1[1], A = x1[2], dn = x1[3] * SMtog.μd / Nits,
    ds = x1[4] * SMtog.μd / Nits)

    return DFinterval, SMtog, [DFinterval[:dn][end] * Nits, DFinterval[:ds][end] * Nits]
end

println("Generating some example data with fits")

#create empty dataframe to store data
myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda, :rsq, :deltatrue, :lambdartrue], 0)

for Δ in [0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5]
    println(Δ)
    for rlam in [0.25, 0.5]
        println(rlam)
        rlambda = rlam
        DF1, SM = simulatepopulation(;Δ = Δ, tend = 30.0, Amin = 0.05, ρ = 100.0, λ = rlambda, Amax = 25.0)
        x = LLoptimizationresults(DF1[:dnds], DF1[:A]; t = 30.0, Amin = 0.05, ρ = 100.0)
        x.DF[:deltatrue] = Δ
        x.DF[:lambdartrue] = rlambda * SM.r
        append!(myDF, x.DF)
    end
end

println("Writing data to file")
@rput myDF
R"""
write_csv(myDF, $(parsed_args["examplefitsout"]))
""";

#####################################################
# Different bin sizes
#####################################################

println("Generating some example data with fits using different bin sizes")

#create empty dataframe to store data
myDF = DataFrame([Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64],
[:dnds, :A, :dndsfit, :dndsfitlq, :dndsfituq, :deltafit, :lambdarfit,
    :deltafitlq, :lambdarfitlq, :deltafituq, :lambdarfituq, :sedelta, :selambda, :rsq, :deltatrue, :lambdartrue, :stepsize], 0)

for Δ in [0.1, 0.2, 0.3, 0.4, 0.5]
    println(Δ)
    for rlam in [0.5]
        for mygap in [0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2]
            Random.seed!(111)
            println(rlam)
            rlambda = rlam
            DF1, SM = simulatepopulation(;Δ = Δ, tend = 30.0, Amin = 0.05, ρ = 100.0, λ = rlambda, Amax = 25.0, gap = mygap)
            x = LLoptimizationresults(DF1[:dnds], DF1[:A]; t = 30.0, Amin = 0.05, ρ = 100.0)
            x.DF[:deltatrue] = Δ
            x.DF[:lambdartrue] = rlambda * SM.r
            x.DF[:stepsize] = mygap
            append!(myDF, x.DF)
        end
    end
end

println("Writing data to file")
@rput myDF
R"""
write_csv(myDF, $(parsed_args["exampledifferentbins"]))
""";


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Transform" begin
    A = [
        1 2 3 4;
        4 5 6 7;
        8 9 10 11;
        12 13 14 15
    ]
    x,y,z,w = Vec4.([A[:,i] for i in 1:4])
    @test Geometry.matrix(Geometry.Transform(x,y,z,w)) == A

    B = [
        1 2 3 4;
        4 5 6 7;
        8 9 10 11;
        0 0 0 1
    ]
    x,y,z,w = Vec3.([B[1:3,i] for i in 1:4])
    @test B == Geometry.matrix(Geometry.Transform(x,y,z,w))

end # testset Allocations


text:
module BITEModel
using VAWTools
import Interpolations
const Interp = Interpolations
import Proj4
using Suppressor
using DocStringExtensions
import DataStructures: OrderedDict
import Roots

const SRangeL{F} = StepRangeLen{F,Base.TwicePrecision{F},Base.TwicePrecision{F}}

# a few short-hands
const F = Float64 # to switch between Float32 and Float64

# Base types
abstract type APara end #{Name} # Name is a symbol to encode the used model (both forward and inverse)

###############
# Misc helpers
###############

include("misc.jl")

################
# Pre-processing
################

# Various data structures for model
include("data-structs.jl")

# All related to data loading
# (probably should be its own package at some stage)
include("data-loading.jl")
include("data-loading-misc.jl")
include("loaders/data-loaders.jl")

###############
# Actual model
###############
# (F&H 2009)

## Forward model implementation, a variation of the Farinotti et al. (2009) model.
include("TwoD.jl")
include("OneD.jl")

# Forward model post-proc:
# plotting, volume,

###############
# Inverse model
###############
# Do some other time...
#
# - Bayesian Fitting
# - Huss ad-hoc fitting

include("mcmc/mcmc.jl")
include("data-loading-inverse.jl")

# """
# TODO

# The inverse model, i.e. the model fitting function.  It needs to
# specify fitting parameters and measurements.  All 2D input/output
# needs to be on the same gird (really?).

# Inputs:
# -------

# Model input with no errors:
# - dem: glacier DEM with ice mask

# Model input with errors (random fields):

# - bdot: mass balance rate
# - dhdt: elevation change rate

# Fitting parameters:

# - fq: function with fitting parameter
# - A: function with fitting parameter

# Measurements of forward model output:

# - us_obs: surface flow speed (usually 2D)
# - h_obs: ice thickness measurements (usually along radar tracks)

# Outputs:
# -------

# - h: ice thickness
# - h_err: error estimate
# - us + us_err: surface flow speed
# - fq + fq_err
# - A + A_err

# """
# function invm end

######
# Post-processing tools
#####
include("postproc.jl")

######
# Plotting
#####
include("plotting.jl")

end # module


text:
using ExcelReaders

function getdice2016excelparameters(filename)
    p = Dict{Symbol,Any}()

    T = 100

    #Open DICE_2016 Excel File to Read Parameters
    f = openxl(filename)

	p[:a0]			= getparams(f, "B108:B108", :single, "Parameters",1)#Initial level of total factor productivity
    p[:a1]          = getparams(f, "B25:B25", :single, "Base", 1)       #Damage coefficient on temperature
    p[:a2]          = getparams(f, "B26:B26", :single, "Base", 1)       #Damage quadratic term
    p[:a3]          = getparams(f, "B27:B27", :single, "Base", 1)       #Damage exponent
    p[:al]          = getparams(f, "B21:CW21", :all, "Base", T)         #Level of total factor productivity
    p[:b12]         = getparams(f, "B67:B67", :single, "Base", 1)       #Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b23]         = getparams(f, "B70:B70", :single, "Base", 1)       #Carbon cycle transition matrix shallow to deep ocean
    p[:c1]          = getparams(f, "B82:B82", :single, "Base", 1)       #Speed of adjustment parameter for atmospheric temperature (per 5 years)
    p[:c3]          = getparams(f, "B83:B83", :single, "Base", 1)       #Coefficient of heat loss from atmosphere to oceans
    p[:c4]          = getparams(f, "B84:B84", :single, "Base", 1)       #Coefficient of heat gain by deep oceans
    p[:cca0]        = getparams(f, "B92:B92", :single, "Base", T)       #Initial cumulative industrial emissions
    p[:cost1]       = getparams(f, "B32:CW32", :all, "Base", T)         #Abatement cost function coefficient
    p[:cumetree0]   = 100							         			#Initial cumulative emissions from deforestation (see GAMS code)
    p[:damadj]      = getparams(f, "B65:B65", :single, "Parameters", 1) #Adjustment exponent in damage function
	p[:dela]		= getparams(f, "B110:B110", :single, "Parameters",1)#Decline rate of TFP per 5 years
	p[:deland]		= getparams(f, "D64:D64", :single, "Parameters", 1) #Decline rate of land emissions (per period)
    p[:dk]          = getparams(f, "B6:B6", :single, "Base", 1)         #Depreciation rate on capital (per year)
	p[:dsig]		= getparams(f, "B66:B66", :single, "Parameters", 1) #Decline rate of decarbonization (per period)
	p[:eland0]		= getparams(f, "D63:D63", :single, "Parameters", 1)	#Carbon emissions from land 2015 (GtCO2 per year)
	p[:e0]			= getparams(f, "B113:B113", :single, "Base", 1)		#Industrial emissions 2015 (GtCO2 per year)
    p[:elasmu]      = getparams(f, "B19:B19", :single, "Base", 1)       #Elasticity of MU of consumption
    p[:eqmat]       = getparams(f, "B82:B82", :single, "Parameters", 1) #Equilibirum concentration of CO2 in atmosphere (GTC)
    p[:expcost2]    = getparams(f, "B39:B39", :single, "Base", 1)       #Exponent of control cost function
    p[:fco22x]      = getparams(f, "B80:B80", :single, "Base", 1)       #Forcings of equilibrium CO2 doubling (Wm-2)
	p[:fex0]		= getparams(f, "B87:B87", :single, "Parameters", 1) #2015 forcings of non-CO2 GHG (Wm-2)
	p[:fex1]		= getparams(f, "B88:B88", :single, "Parameters", 1) #2100 forcings of non-CO2 GHG (Wm-2)
    p[:fosslim]     = getparams(f, "B57:B57", :single, "Base", 1)       #Maximum carbon resources (Gtc)
	p[:ga0]			= getparams(f, "B109:B109", :single, "Parameters",1)#Initial growth rate for TFP per 5 years
    p[:gama]        = getparams(f, "B5:B5", :single, "Base", 1)         #Capital Share
	p[:gback]		= getparams(f, "B26:B26", :single, "Parameters", 1) #Initial cost decline backstop cost per period
	p[:gsigma1]		= getparams(f, "B15:B15", :single, "Parameters", 1)	#Initial growth of sigma (per year)
    p[:k0]          = getparams(f, "B12:B12", :single, "Base", 1)       #Initial capital
    p[:l]           = getparams(f, "B53:CW53", :all, "Base", T)         #Level of population and labor (millions)
    p[:mat0]        = getparams(f, "B61:B61", :single, "Base", 1)       #Initial Concentration in atmosphere in 2015 (GtC)
	p[:mateq]		= getparams(f, "B82, B82", :single, "Parameters", 1)#Equilibrium concentration atmosphere  (GtC)
    p[:MIU]         = getparams(f, "B135:CW135", :all, "Base", T)       #Optimized emission control rate results from DICE2016R (base case)
    p[:ml0]         = getparams(f, "B63:B63", :single, "Base", 1)       #Initial Concentration in deep oceans 2010 (GtC)
	p[:mleq]		= getparams(f, "B84, B84", :single, "Parameters", 1)#Equilibrium concentration in lower strata (GtC)
    p[:mu0]         = getparams(f, "B62:B62", :single, "Base", 1)       #Initial Concentration in biosphere/shallow oceans 2010 (GtC)
	p[:mueq]		= getparams(f, "B83:B83", :single, "Parameters", 1) #Equilibrium concentration in upper strata (GtC)
    p[:pback]	    = getparams(f, "B10:B10", :single, "Parameters", 1) #Cost of backstop 2010$ per tCO2 2015
    p[:rr]          = getparams(f, "B18:CW18", :all, "Base", T)         #Social Time Preference Factor
    p[:S]           = getparams(f, "B131:CW131", :all, "Base", T)       #Optimized savings rate (fraction of gross output) results from DICE2016 (base case)
    p[:scale1]      = getparams(f, "B49:B49", :single, "Base", 1)       #Multiplicative scaling coefficient
    p[:scale2]      = getparams(f, "B50:B50", :single, "Base", 1)       #Additive scaling coefficient
    p[:t2xco2]      = getparams(f, "B79:B79", :single, "Base", 1)       #Equilibrium temp impact (oC per doubling CO2)
    p[:tatm0]       = getparams(f, "B76:B76", :single, "Base", 1)       #Initial atmospheric temp change 2015 (C from 1940-60)
    p[:tocean0]     = getparams(f, "B77:B77", :single, "Base", 1)       #Initial temperature of deep oceans (deg C above 1940-60)

    return p
end

function getdice2016gamsparameters(filename)
    p = Dict{Symbol,Any}()

    T = 100

    #Open DICE_2016 Excel File to Read Parameters
    f = openxl(filename)
    sheet = "DICE2016_Base"

	p[:a0]			= 5.115											   #Initial level of total factor productivity
    p[:a1]          = getparams(f, "B41:B41", :single, sheet, 1)       #Damage coefficient on temperature
    p[:a2]          = getparams(f, "B42:B42", :single, sheet, 1)       #Damage quadratic term
    p[:a3]          = getparams(f, "B43:B43", :single, sheet, 1)       #Damage exponent
    p[:al]          = getparams(f, "B5:CWI5", :all, sheet, T)          #Level of total factor productivity
    p[:b12]         = getparams(f, "B20:B20", :single, sheet, 1)       #Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b23]         = getparams(f, "B21:B21", :single, sheet, 1)       #Carbon cycle transition matrix shallow to deep ocean
    p[:c1]          = getparams(f, "B36:B36", :single, sheet, 1)       #Speed of adjustment parameter for atmospheric temperature (per 5 years)
    p[:c3]          = getparams(f, "B37:B37", :single, sheet, 1)       #Coefficient of heat loss from atmosphere to oceans
    p[:c4]          = getparams(f, "B38:B38", :single, sheet, 1)       #Coefficient of heat gain by deep oceans
    p[:cca0]        = getparams(f, "B14:B14", :single, sheet, 1)       #Initial cumulative industrial emissions
    p[:cost1]       = getparams(f, "B46:CW46", :all, sheet, T)         #Abatement cost function coefficient
    p[:cumetree0]   = 100											   #Initial cumulative emissions from deforestation (see GAMS code)
    p[:damadj]      = getparams(f, 10.)                                #Adjustment exponent in damage function
	p[:dela]		= 0.0050										   #Decline rate of TFP per 5 years
	p[:deland]		= 0.115											   #Decline rate of land emissions (per period)
    p[:dk]          = getparams(f, "B8:B8", :single, sheet, 1)         #Depreciation rate on capital (per year)
	p[:dsig]		= -0.001										   #Decline rate of decarbonization (per period)
	p[:eland0]		= 2.6										   	   #Carbon emissions from land 2015 (GtCO2 per year)
	p[:e0]			= 35.74471365402390								   #Industrial emissions 2015 (GtCO2 per year)
    p[:elasmu]      = getparams(f, "B54:B54", :single, sheet, 1)       #Elasticity of MU of consumption
    p[:eqmat]       = getparams(f, 588.)                               #Equilibirum concentration of CO2 in atmosphere (GTC)
    p[:expcost2]    = getparams(f, "B48:B48", :single, sheet, 1)       #Exponent of control cost function
    p[:fco22x]      = getparams(f, "B30:B30", :single, sheet, 1)       #Forcings of equilibrium CO2 doubling (Wm-2)
	p[:fex0]		= 0.5											   #2015 forcings of non-CO2 GHG (Wm-2)
	p[:fex1]		= 1.0											   #2100 forcings of non-CO2 GHG (Wm-2)
	p[:ga0]			= 0.0760										   #Initial growth rate for TFP per 5 years
    p[:gama]        = getparams(f, "B7:B7", :single, sheet, 1)         #Capital Share
	p[:gback]		= 0.025										 	   #Initial cost decline backstop cost per period
	p[:gsigma1]		= -0.0152										   #Initial growth of sigma (per year)
    p[:k0]          = getparams(f, "B9:B9", :single, sheet, 1)         #Initial capital
    p[:l]           = getparams(f, "B6:CW6", :all, sheet, T)           #Level of population and labor (millions)
    p[:mat0]        = getparams(f, "B17:B17", :single, sheet, 1)       #Initial Concentration in atmosphere in 2015 (GtC)
	p[:mateq]		= getparams(f, "B82, B82", :single, "Parameters",1)#Equilibrium concentration atmosphere  (GtC)
    p[:MIU]         = getparams(f, "B47:CW47", :all, sheet, T)         #Optimized emission control rate results from DICE2016R (base case)
    p[:ml0]         = getparams(f, "B19:B19", :single, sheet, 1)       #Initial Concentration in deep oceans 2015 (GtC)
	p[:mleq]		= getparams(f, "B84, B84", :single, "Parameters",1)#Equilibrium concentration in lower strata (GtC)
    p[:mu0]         = getparams(f, "B18:B18", :single, sheet, 1)       #Initial Concentration in biosphere/shallow oceans 2015 (GtC)
	p[:mueq]		= getparams(f, "B83:B83", :single, "Parameters", 1)#Equilibrium concentration in upper strata (GtC)
    p[:partfract]   = getparams(f, "B49:CW49", :all, sheet, T)         #Fraction of emissions in control regime
    p[:pback]	    = 550											   #Cost of backstop 2010$ per tCO2 2015
    p[:rr]          = getparams(f, "B55:CW55", :all, sheet, T)         #Social Time Preference Factor
    p[:S]           = getparams(f, "B51:CW51", :all, sheet, T)         #Optimized savings rate (fraction of gross output) results from DICE2016R (base case)
    p[:scale1]      = getparams(f, "B56:B56", :single, sheet, 1)       #Multiplicative scaling coefficient
    p[:scale2]      = getparams(f, "B57:B57", :single, sheet, 1)       #Additive scaling coefficient
    p[:t2xco2]      = getparams(f, "B33:B33", :single, sheet, 1)       #Equilibrium temp impact (oC per doubling CO2)
    p[:tatm0]       = getparams(f, "B34:B34", :single, sheet, 1)       #Initial atmospheric temp change (C from 1900)
    p[:tocean0]     = getparams(f, "B35:B35", :single, sheet, 1)       #Initial temperature of deep oceans (deg C above 1900)

    return p
end


text:
using GeometryTypes, StaticArrays, Colors, GLAbstraction

"""
Hack to quickly make things more consistent inside Makie, without
changing GLVisualize too much! So we need to rewrite the attributes, the names and the
values a bit!
"""
function expand_for_glvisualize(kw_args)
    result = Dict{Symbol, Any}()
    for (k, v) in kw_args
        k in (:marker, :positions, always_skip...) && continue

        if k == :rotations
            k = :rotation
            v = Vec4f0(0, 0, 0, 1)
            result[:billboard] = true
        end
        if k == :markersize
            k = :scale
        end
        if k == :glowwidth
            k = :glow_width
        end
        if k == :glowcolor
            k = :glow_color
        end
        if k == :strokewidth
            k = :stroke_width
        end
        if k == :strokecolor
            k = :stroke_color
        end
        if k == :positions
            k = :position
        end
        result[k] = to_signal(v)
    end
    result[:fxaa] = false
    result
end

# function scatter()
#     marker
#     strokecolor
#     strokewidth
#     glowcolor
#     glowwidth
#     markersize
#     rotations
#     if haskey(scene, :positions)
#         positions = to_positions(positions)
#     elseif haskey(scene, :z)
#         xyz = getindex.(scene, (:x, :y, :z))
#         positions = to_positions(xyz)
#     elseif haskey(scene, :y)
#         xy = getindex.(scene, (:x, :y, :z))
#         positions = to_positions(xy)
#     end
#     begin
#         color = to_color(color)
#     end
#     begin
#         colormap = to_colormap(colormap)
#         intensity = to_intensity(intensity)
#         colornorm = to_colornorm(colornorm, intensity)
#     end
# end
#
function _scatter(scene, kw_args)
    attributes = scatter_defaults(scene, kw_args)
    gl_data = expand_for_glvisualize(attributes)
    shape = to_signal(attributes[:marker])
    main = (shape, to_signal(attributes[:positions]))
    viz = visualize(main, Style(:default), gl_data)
    insert_scene!(scene, :scatter, viz, attributes)
end

function mesh2glvisualize(kw_args)
    result = Dict{Symbol, Any}()
    for (k, v) in kw_args
        k in (:marker, :positions, always_skip...) && continue
        if k == :rotations
            k = :rotation
        end
        if k == :markersize
            k = :scale
        end
        if k == :positions
            k = :position
        end
        result[k] = to_signal(v)
    end
    result[:fxaa] = true
    result
end

function _meshscatter(scene, kw_args)
    attributes = meshscatter_defaults(scene, kw_args)
    gl_data = mesh2glvisualize(attributes)
    shape = to_signal(attributes[:marker])
    main = (shape, to_signal(attributes[:positions]))
    viz = GLVisualize.meshparticle(main, Style(:default), gl_data)
    viz = GLVisualize.assemble_shader(viz).children[]
    insert_scene!(scene, :meshscatter, viz, attributes)
end

for arg in ((:x, :y), (:x, :y, :z), (:positions,))
    insert_expr = map(arg) do elem
        :(attributes[$(QuoteNode(elem))] = $elem)
    end
    @eval begin
        function scatter(scene::makie, $(arg...), attributes::Dict)
            $(insert_expr...)
            _scatter(scene, attributes)
        end
        @eval begin
            function meshscatter(scene::makie, $(arg...), attributes::Dict)
                $(insert_expr...)
                _meshscatter(scene, attributes)
            end
        end
    end
end


text:
using CSDP

#=
# Example copied from `example/example.c`
   An example showing how to call the easy_sdp() interface to CSDP.  In
   this example, we solve the problem

      max tr(C*X)
          tr(A1*X)=1
          tr(A2*X)=2
          X >= 0       (X is PSD)

   where

    C=[2 1
       1 2
           3 0 1
           0 2 0
           1 0 3
                 0
                   0]

   A1=[3 1
       1 3
           0 0 0
           0 0 0
           0 0 0
                 1
                   0]

   A2=[0 0
       0 0
           3 0 1
           0 4 0
           1 0 5
                 0
                   1]

  Notice that all of the matrices have block diagonal structure.  The first
  block is of size 2x2.  The second block is of size 3x3.  The third block
  is a diagonal block of size 2.

=#

C = BlockMatrix(
   [2 1
    1 2],
        [3 0 1
         0 2 0
         1 0 3],
      Diagonal([0,
                  0]))

b = [1.0, 2.0]

A1 = ConstraintMatrix(1,
   [3 1
    1 3],
          [0 0 0
           0 0 0
           0 0 0],
       Diagonal([1,
                   0]))

A2 = ConstraintMatrix(2,
      [0 0
       0 0],
          [3 0 1
           0 4 0
           1 0 5],
       Diagonal([0,
                   1]))

constraints = [A.csdp for A in [A1, A2]]

CSDP.write_prob("prob.dat-s", C, b, constraints)

X, y, Z = initsoln(C, b, constraints)

pobj, dobj = easy_sdp(C, b, constraints, X, y, Z)

# Fails on Windows
#CSDP.write_sol("prob.sol", X, y, Z)