text:
#=
disk-phantom.jl
generate random disk phantoms
2019-07-03, Jeff Fessler, University of Michigan
=#

export disk_phantom_params

# using MIRT: jim
using Random: seed!

"""
`params = disk_phantom_params( ; ...)`

Generate ellipse phantom parameters for a head-sized disk
plus many disks within it,
designed so that the disks have some minimum separation `minsep`
to avoid overlap and to simplify patch-based model fitting.
# Arguments
- `fov::Real = 240` image field of view in mm
- `rhead::Real = 100` background radius for "head" [mm]
- `muhead::Real = 1000` "mu" (intensity) value for background head disk
- `rmin::Real = 10` min radius for random disks
- `rmax::Real = 20` max radius for random disks
- `mumin::Real = 100` range of "mu" values for disks
- `mumax::Real = 300`
- `ndisk::Integer = 10` # of random disks
- `minsep::Real = 8` minimum disk separation in mm
- `maxtry::Integer = 500` give up on adding more disks if this is reached
- `warn::Bool = false` warn if maxtry reached?
- `seed::Integer = 0`if nonzero then use this seed
"""
function disk_phantom_params( ;
		fov::Real = 240,
		rhead::Real = 100, # background radius for "head" [mm]
		muhead::Real = 1000, # "mu" value for background head
		rmin::Real = 10, # min radius for disks
		rmax::Real = 20, # max radius for disks
		mumin::Real = 100, # range of "mu" values for disks
		mumax::Real = 300,
		ndisk::Integer = 10,
		minsep::Real = 8, # minimum separation in mm
		maxtry::Integer = 500, # give up on adding more disks if this is reached
		warn::Bool = false, # warn if maxtry reached?
		seed::Integer = 0,
	)

	params = zeros(Float32, ndisk+1, 6)
	params[end,:] = [0, 0, rhead, rhead, 0, muhead]

	randu = (a, b; f = x->x) -> a + f(rand()) * (b - a)
	seper = (a,b,r1,x,y,r2) -> sqrt.((x .- a).^2 + (y .- b).^2) .- r1 .- r2

	(seed != 0) && seed!(seed)

	for id=1:ndisk
		trial = 0
		for ii=1:maxtry
			rc = randu(0, rhead - rmin - minsep, f = x->x^0.5)
			phi = rand() * 2*pi
			rad = randu(rmin, rmax)
			val = randu(mumin, mumax)
			trial = [rc*cos(phi), rc*sin(phi), rad, rad, 0, val]

			id == 1 && break # first one is always fine

			# see if trial is too close to another one
			c = 1:(id-1) # check these
			sep = seper(params[c,1], params[c,2], params[c,3],
				trial[1], trial[2], trial[3])

			(minimum(sep) > minsep) && (rhead - (rc + rad) > minsep) && break
			ii == maxtry && warn && @warn("need more tries")
			ii == maxtry && return params[[1:id; end],:]
		end

#=
		if id > 1
			c = 1:(id-1) # check these
			sep = seper(params[c,1], params[c,2], params[c,3],
				trial[1], trial[2], trial[3])
			@show minimum(sep)
		end
		@show id, trial
=#

		params[id,:] = trial
	end

	return params
end

"""
`disk_phantom_params(:test)`
self test

`disk_phantom_params(:show)`
show example
"""
function disk_phantom_params(test::Symbol)
	ig = image_geom(nx = 128, fov = 240)
	params = disk_phantom_params( ; minsep = ig.dx*8)

	if test == :show
		tmp = ellipse_im(ig, params, oversample=3)
		return jim(ig.x, ig.y, tmp)
	end

	test != :test && throw(ArgumentError("test $test"))
	disk_phantom_params(:show)
	true
end

# disk_phantom_params(:test)


text:
using Documenter, Example

makedocs(sitename="My Documentation")

ext: jl


text:
using InteractionNetworkModels

move = PathPermutationMove(ν=2)

V = 1:10
S_curr = [rand(V, i) for i in 1:5]
S_prop = deepcopy(S_curr)
pointers = [Int[] for i in 1:100]

prop_sample!(S_curr, S_prop, move, pointers, V)

S_curr
S_prop


text:
using Gallium
using Gallium.X86_64: basic_regs, BasicRegs, dwarf_numbering
function rebuild_lib(file, triple, oname)
    contents = readstring(file)
    for idx in basic_regs
        pat = string("UC_MCONTEXT_GREGS_",uppercase(string(dwarf_numbering[idx])))
        @show pat
        contents = replace(contents,
            pat, string(idx*sizeof(Ptr{Void})))
    end
    xsave_area_size = sizeof(fieldtype(Gallium.X86_64.ExtendedRegs,:xsave_state))
    basic_size = length(basic_regs)*sizeof(Ptr{Void})
    contents = replace(contents, "UC_MCONTEXT_SIZE", basic_size)
    contents = replace(contents, "UC_MCONTEXT_TOTAL_SIZE", xsave_area_size + basic_size)
    print(contents)
    cmd = `$(joinpath(JULIA_HOME,"../tools/llvm-mc")) -triple=$triple -filetype=obj - -o $oname`
    (stdout,stdin,process) = readandwrite(cmd)
    write(stdin, contents)
    close(stdin)
    wait(process)
    println(readstring(stdout))
end

rebuild_lib("getcontext-x86_64-elf.s","x86_64-pc-linux-gnu","elfhook.o")
rebuild_lib("jumpto-x86_64-elf.s","x86_64-pc-linux-gnu","elfjump.o")
run(`clang -target x86_64-pc-linux-gnu -fPIC -c -o elf-callback.o callback.c`)
@linux_only run(`clang -target x86_64-pc-linux-gnu -shared -o hooking.so elfjump.o elfhook.o elf-callback.o`)

rebuild_lib("getcontext-x86_64-macho.s","x86_64-apple-darwin15.3.0","machohook.o")
rebuild_lib("jumpto-x86_64-macho.s","x86_64-apple-darwin15.3.0","machojump.o")
run(`clang -target x86_64-apple-darwin15.3.0 -fPIC -c -o macho-callback.o callback.c`)
@osx_only run(`clang -target x86_64-apple-darwin15.3.0 -shared -o hooking.dylib machojump.o machohook.o macho-callback.o`)

windows_triple = "x86_64-pc-windows-gnu"
rebuild_lib("getcontext-x86_64-coff.s",windows_triple,"coffhook.o")
rebuild_lib("jumpto-x86_64-coff.s",windows_triple,"coffjump.o")
run(`clang -target $windows_triple -fPIC -c -o coff-callback.o callback.c`)
@windows_only run(`clang -target $windows_triple -shared -o hooking.dll coffjump.o coffhook.o coff-callback.o`)


text:
#!/usr/bin/env julia

using DataFrames

import Celeste.AccuracyBenchmark
import Celeste.ArgumentParse

parser = ArgumentParse.ArgumentParser()
ArgumentParse.add_argument(
    parser,
    "--write-prediction-csv",
    help="Write matched predictions to the given CSV file",
)
ArgumentParse.add_argument(
    parser,
    "ground_truth",
    help="Ground truth CSV catalog",
)
ArgumentParse.add_argument(
    parser,
    "first_predictions",
    help="CSV catalog of predictions",
)
ArgumentParse.add_argument(
    parser,
    "second_predictions",
    help="Second CSV catalog of predictions, to compare to the first",
    required=false,
)
parsed_args = ArgumentParse.parse_args(parser, ARGS)

truth = AccuracyBenchmark.read_catalog(parsed_args["ground_truth"])
prediction_dfs = [AccuracyBenchmark.read_catalog(parsed_args["first_predictions"])]
if haskey(parsed_args, "second_predictions")
    push!(prediction_dfs, AccuracyBenchmark.read_catalog(parsed_args["second_predictions"]))
end
scores = AccuracyBenchmark.score_predictions(truth, prediction_dfs)
println(repr(scores))

if haskey(parsed_args, "write-prediction-csv")
    matched_truth, matched_prediction_dfs = AccuracyBenchmark.match_catalogs(truth, prediction_dfs)
    for prediction_df in matched_prediction_dfs
        for name in names(prediction_df)
            if !in(name, names(truth))
                delete!(prediction_df, name)
            end
        end
    end
    matched_truth[:index] = collect(1:nrow(matched_truth))
    matched_truth[:source] = fill("truth", size(matched_truth, 1))
    for (i, prediction_df) in enumerate(matched_prediction_dfs)
        prediction_df[:index] = collect(1:nrow(prediction_df))
        prediction_df[:source] = fill("prediction $i", size(prediction_df, 1))
    end
    all_predictions = vcat(matched_truth, matched_prediction_dfs...)
    long_df = melt(all_predictions, [:index, :source])
    long_df[:index_var] = [
        join([idx, variable], " ")
        for (idx, variable) in zip(long_df[:index], long_df[:variable])
    ]
    delete!(long_df, :index)
    delete!(long_df, :variable)
    final_df = unstack(long_df, :index_var, :source, :value)
    writetable(parsed_args["write-prediction-csv"], final_df)
end


text:
# set model type
mutable struct IRG
    alpha::Array{Float64,1}
    beta::Array{Float64,1}
    gamma::Array{Float64,1}
    sigma::Array{Float64, 1}
    psi::Float64
    delta::Float64
end

# functions in this model
# vector valued functions
function u(model::IRG, s::Array{Float64, 2}, x::Array{Float64, 2}, p::Int64)
    return ((s[:, p] - x[:, p]).^(1 - model.alpha[p]))./(1 - model.alpha[p])
end

function h(model::IRG, x::Array{Float64, 2}, epsilon::Float64, p::Int64)
    return model.gamma[p] * x[:, p] + epsilon * (x[:, p].^(model.beta[p]))
end

function s(model::IRG, x::Array{Float64, 2}, epsilon::Array{Float64, 1})
    box = zeros((size(x)[1], 2))
    box[:, 1] = (1 - model.psi)*h(model, x, epsilon[1], 1) + model.psi*h(model, x, epsilon[2], 2)
    box[:, 2] = (1 - model.psi)*h(model, x, epsilon[2], 2) + model.psi*h(model, x, epsilon[1], 1)
    return box
end

# derivatives of  payoff, transition functions
# utility function
function ux(model::IRG, s::Array{Float64, 2}, x::Array{Float64, 2}, p::Int64)
    return -(s[:, p] - x[:, p]).^(-model.alpha[p])
end

function uxx(model::IRG, s::Array{Float64, 2}, x::Array{Float64, 2}, p::Int64)
    return -model.alpha[p]*(s[:, p] - x[:, p]).^(-model.alpha[p]-1)
end

# transition function
function sx(model::IRG, x::Array{Float64, 2}, epsilon::Array{Float64, 1},  p::Int64)
    return (1 - model.psi)*(model.gamma[p] + epsilon[p]*model.beta[p]*x[:, p].^(model.beta[p] - 1))
end

function sxx(model::IRG, x::Array{Float64, 2}, epsilon::Array{Float64, 1},  p::Int64)
    return (1 - model.psi)*epsilon[p]*model.beta[p]*(model.beta[p] - 1)*x[:, p].^(model.beta[p] - 2)
end

# Gaussian Quadrature
function gq(model::IRG, num_nodes::Int64)
    return qnwlogn([num_nodes, num_nodes], [0,0], diagm([sigma^2, sigma^2]))
end

# collocation function
function vmax(model::IRG, colnodes::Array{Float64, 2}, action::Array{Float64, 2}, coef::Array{Float64, 2}, epss::Array{Float64, 2}, weights::Array{Float64, 1})
    xnew = action
    v = zeros((size(colnodes)[1], 2))
    for p in 1:2
        xl, xu = 0.0, colnodes[:, p]
        order1 = [0 0]
        order1[1, p] = 1
        order2 = [0 0]
        order2[1, p] = 2
        for it in 1:maxit
            util, util_der1, util_der2 = u(model, colnodes, action, p), ux(model, colnodes, action, p), uxx(model, colnodes, action, p)
            # println(util[1])
            # println(util_der1[1])
            # println(util_der2[1])
            println(action)
            Ev, Evx, Evxx = 0.0, 0.0, 0.0
            for k in 1:num_nodes^2
                println(k)
                eps, weight= epss[k, :], weights[k]
                transition, transition_der1, transition_der2 = s(model, action, eps), sx(model, action, eps, p), sxx(model, action, eps, p)
                #println(transition)
                #println(transition_der1)
                #println(transition_der2)
                # coefによって結果が変わらないんだけど、係数はここにしか出てないからこいつがやっぱり怪しいよな
                # 係数は別に良いっぽいダメなのはfuneval？
                vn = funeval(coef[:, p], basis, transition)
                vnder1 =  funeval(coef[:, p], basis, transition, order1)
                vnder2 = funeval(coef[:, p], basis, transition, order2)
                # vnが係数と一致してるんだけどなにこれ
                #println(vn)
                #println(vnder1)
                #println(vnder2)
                Ev += weight * vn
                Evx += weight* vnder1.* transition_der1
                Evxx += weight * (vnder1.*transition_der2 + vnder2 .* (transition_der1.^2))
            end
            # println(Ev[1])
            # println(Evx[1])
            # println(Evxx[1])
            v[:, p] = util + Ev
            delx = -(util_der1 + model.delta * Evx) ./ (util_der2 + model.delta*Evxx)
            # delxがでかすぎる
            #println(delx)
            # println(xl-action[:, p])
            # println(xu-action[:, p])
            delx = min.(max.(delx, xl-action[:, p]), xu-action[:, p])
            action[:, p] = action[:, p] + delx
            if norm(delx) < tol
                break
            end
        end
        xnew[:, p] = action[:, p]
    end
    return v, xnew
end

max_stars_repo_name: keiikegami/IRG

max_issues_repo_name: keiikegami/IRG

max_forks_repo_name: keiikegami/IRG


text:
#= ================ Code copied from BinaryBuilder.jl ============== =#
#=
Since this code is not part of the public BinaryBuilder interface, it is safer to copy it than to use it directly by adding a dependency to BinaryBuider
=#

# This is a global github authentication token that is set the first time
# we authenticate and then reused
const _github_auth = Ref{GitHub.Authorization}()

function github_auth(;allow_anonymous::Bool=true)
    if !isassigned(_github_auth) || !allow_anonymous && isa(_github_auth[], GitHub.AnonymousAuth)
        # If the user is feeding us a GITHUB_TOKEN token, use it!
        if length(get(ENV, "GITHUB_TOKEN", "")) == 40
            _github_auth[] = GitHub.authenticate(ENV["GITHUB_TOKEN"])
        else
            if allow_anonymous
                _github_auth[] = GitHub.AnonymousAuth()
            else
                # If we're not allowed to return an anonymous authorization,
                # then let's create a token.
                _github_auth[] = GitHub.authenticate(obtain_token())
            end
        end
    end

    return _github_auth[]
end

function obtain_token(; outs=stdout, github_api=GitHub.DEFAULT_API)
    println(outs)
    printstyled(outs, "Authenticating with GitHub\n", bold=true)

    @label retry

    headers = Dict{String, String}("User-Agent"=>"BinaryBuilder-jl",
        "Accept"=>"application/json")

    # Request device authentication flow for BinaryBuilder OAauth APP
    resp = HTTP.post("https://github.com/login/device/code", headers,
        "client_id=8c496f428c48a1015b9e&scope=public_repo")
    if resp.status != 200
        GitHub.handle_response_error(resp)
    end
    reply = JSON.parse(HTTP.payload(resp, String))

    println(outs, """
    To continue, we need to authenticate you with GitHub. Please navigate to
    the following page in your browser and enter the code below:

         $(HTTP.URIs.unescapeuri(reply["verification_uri"]))

         #############
         # $(reply["user_code"]) #
         #############
    """)

    interval = reply["interval"]
    device_code = reply["device_code"]
    while true
        # Poll for completion
        sleep(interval)
        resp = HTTP.post("https://github.com/login/oauth/access_token", headers,
            "client_id=8c496f428c48a1015b9e&grant_type=urn:ietf:params:oauth:grant-type:device_code&device_code=$device_code")

        if resp.status != 200
            GitHub.handle_response_error(resp)
        end

        token_reply = JSON.parse(HTTP.payload(resp, String))
        if haskey(token_reply, "error")
            error_kind = token_reply["error"]
            if error_kind == "authorization_pending"
                continue
            elseif error_kind == "slow_down"
                @warn "GitHub Auth rate limit exceeded. Waiting 10s. (This shouldn't happen)"
                sleep(10)
            elseif error_kind == "expired_token"
                @error "Token request expired. Starting over!"
                @goto retry
            elseif error_kind == "access_denied"
                @error "Authentication request canceled by user. Starting over!"
                @goto retry
            elseif error_kind in ("unsupported_grant_type",
                "incorrect_client_credentials", "incorrect_device_code")
                error("Received error kind $(error_kind). Please file an issue.")
            else
                error("Unexpected GitHub login error $(error_kind)")
            end
        end

        token = token_reply["access_token"]

        print(outs, """
        Successfully obtained GitHub authorization token!
        This token will be used for the rest of this BB session.
        You will have to re-authenticate for any future session.
        However, if you wish to bypass this step, you may create a
        personal access token at """)
        printstyled("https://github.com/settings/tokens"; bold=true)
        println("\n and add the token to the")
        printstyled(outs, "~/.julia/config/startup.jl"; bold=true)
        println(outs, " file as:")
        println(outs)

        printstyled(outs, "    ENV[\"GITHUB_TOKEN\"] = <token>"; bold=true)
        println(outs)

        println(outs, "This token is sensitive, so only do this in a computing environment you trust.")
        println(outs)

        return token
    end
end


text:
using Base.Test

import CompScienceMeshes; CM = CompScienceMeshes
import BEAST; BE = BEAST

ω = 1.0
ui, wi = legendre(4, 0.0, 1.0); ui = transpose(ui)
uo, wo = legendre(3, 0.0, 1.0); uo = transpose(uo)
universe = Universe(ω, ui, wi, uo, wo)

fn = Pkg.dir("BEAST", "test", "sphere.in")
mesh = CM.meshfromfile(fn)

rwg = BE.raviartthomas(m);
bc = BE.buffachristiansen(m);

T = BE.MWSingleLayer3D(universe);

trr = BE.assemble(rwg, rwg, t);
tbb = BE.assemble(bc, bc, t);


text:
function plot_net(df::DataFrame)

    # drop the helper columns
    select!(df, 1:findfirst(contains("Column"), names(df))-1)

    # figure
    plot()
    totals = ["GBP", "USD", "EUR", "total (GBP)"]
    individuals = filter!(c -> !(c in totals || c == "date"), names(df))
    for ind in individuals
        plot!(df.date, df[!, ind]*0.001; label=ind, ls=:dash, lw=2.0)
    end
    for tot in totals
        plot!(df.date, df[!, tot]*0.001; label=tot, lw=3)
    end

    # compute days to FI
    Δ_wealth = df[end, "total (GBP)"] - df[1, "total (GBP)"]
    Δ_days = df[end, "date"] - df[1, "date"]
    rate_of_change = Δ_wealth / Δ_days.value
    missing_wealth = 500000 - df[end, "total (GBP)"]
    days2fi = Day(round(missing_wealth / rate_of_change))
    fi_date = today() + days2fi
    bday = Date(1993, 6, 30)
    age_fi = "$(year(fi_date) - year(bday))y, $(month(fi_date) - month(bday))m"

    # cosmetics
    ylabel!("Net worth (k)")
    ylims!((0, ylims()[2]))
    title!("$days2fi to FI (£500k on $fi_date), aged $age_fi")
    plot!(legend=:topleft)

    # save
    savefig("figures/networth.pdf")
end

function plot_categories(expenses::DataFrame)

    # sort the categories by total spent
    categories = Exchequer.aggregate(expenses, "category")
    sort!(categories, :amount_sum, rev=true)

    nmonths = months(expenses)
    currency = expenses.currency[1]

    # figure
    bar(categories.category, categories.amount_sum / nmonths,
        xrotation=20, fillcolor=COLORS, legend=nothing, width=0.5)

    # cosmetics
    ylabel!("Monthly spending ($currency)")
    ylims!((0, ylims()[2]))

    # save
    savefig("figures/categories.pdf")

end

function plot_yearly(expenses::DataFrame, incomes::DataFrame)

    # sort the categories by total spent
    categories = Exchequer.aggregate(expenses, "category")
    sort!(categories, :amount_sum, rev=true)

    # group expenses by year and month and category
    exp = Exchequer.aggregate(expenses, "year", "category")

    # and fill in zeros to categories
    fill_zeros!(exp)

    # group incomes by year and month
    inc = Exchequer.aggregate(incomes, "year")

    # compute number of months in each year (first and last != 12)
    months_in_year = Dict(y=>length(unique(expenses[expenses.year.==y, "month"]))
                          for y in unique(expenses.year))
end

function plot_monthly(expenses::DataFrame, incomes::DataFrame)
    # sort the categories by total spent
    categories = Exchequer.aggregate(expenses, "category")
    sort!(categories, :amount_sum, rev=true)

    # group expenses by year and month and category
    exp_groups = Exchequer.aggregate(expenses, "year", "month", "category")
    fill_zeros!(exp_groups)

    # group incomes by year and month
    inc_groups = Exchequer.aggregate(incomes, "year", "month")
    fill_zeros!(inc_groups)

    # make sure they have the same starting and ending point in time
    fill_zeros!(exp_groups, inc_groups)

    # make figures
    currency = expenses.currency[1]
    plot_unstacked(categories, exp_groups, currency)
    plot_stacked(categories, exp_groups, inc_groups, currency)

    return inc_groups
end

function plot_unstacked(categories::DataFrame, groups::DataFrame, currency::String)
    # clear the figure
    plot()

    # and show every category separately
    spending = groups[groups.category .== "sport", :]
    xpos = [i for i in 1:1:nrow(spending)]
    for (cat, color) in zip(categories.category, COLORS)
        # select the spending and construct dates
        spending = groups[groups.category .== cat, :]
        plot!(xpos, spending.amount_sum, label=cat, color=color, linewidth=2)
    end

    # xticks cosmetics
    yearmonths = zip(spending.year, spending.month)
    labels = [m == 1 ? "$(MONTHS[m])\n$y" : "$(MONTHS[m])" for (y,m) in yearmonths]
    xticks!(xpos, labels)

    # other
    ylabel!("Monthly spending ($currency)")
    ylims!(0, ylims()[2])
    xlims!(1, xlims()[2])
    plot!(legend=:topright)

    # save
    savefig("figures/monthly_unstacked.pdf")
end

function plot_stacked(categories::DataFrame,
                      exp_groups::DataFrame,
                      inc_groups::DataFrame,
                      currency::String)
    # clear the figure
    plot()

    # get the upper and lower limits for each category
    ncat = nrow(categories)
    nmonth = nrow(exp_groups) ÷ ncat
    ys = zeros(ncat+1, nmonth)
    for (i, c) in enumerate(categories.category)
        ys[i+1, :] = exp_groups[exp_groups.category .== c, "amount_sum"]
    end
    cs = cumsum(ys, dims=1)
    lowers = cs[1:end-1, :]
    uppers = cs[2:end, :]

    total_expense = uppers[end,:]
    total_income = inc_groups.amount_sum

    # plot mains
    xpos = [i for i in 1:1:nmonth]
    for (i, c) in enumerate(categories.category)
        plot!(xpos, lowers[i, :], fillrange=uppers[i, :],
              label=c, color=COLORS[i], alpha=0.8, width=0.01)
    end

    # plot income
    plot!(xpos, inc_groups.amount_sum,
          linewidth=2, linecolor="black", linestyle=:dot, label="income")

    # plot net numbers
    net = [trunc(Int, n) for n in total_income - total_expense]
    ypos = [max(a, b)+50 for (a, b) in zip(total_income, total_expense)]
    texts = [(x, y, text(n, 7, :left, n > 0 ? :green : :red; rotation=60,))
             for (x, y, n) in zip(xpos, ypos, net)]
    annotate!(texts)

    # xticks cosmetics
    spending = exp_groups[exp_groups.category .== "sport", :]
    yearmonths = zip(spending.year, spending.month)
    labels = [m == 1 ? "$(MONTHS[m])\n$y" : "$(MONTHS[m])" for (y,m) in yearmonths]
    xticks!(xpos, labels)

    # other cosmetics
    ylabel!("Monthly spending ($currency)")
    ylims!(0, ylims()[2]*1.05)
    xlims!(1, xlims()[2])
    plot!(legend=:topright)

    savefig("figures/monthly_stacked.pdf")
end

function plot_weekly(expenses::DataFrame)

    # sort the categories by total spent
    categories = Exchequer.aggregate(expenses, "category")
    sort!(categories, :amount_sum, rev=true)

    # day of week, normalised to per day spend
    day_exp = Exchequer.aggregate(expenses, "category", "dayofweek")
    ndays = Dates.value(maximum(expenses.date) - minimum(expenses.date))
    day_exp.amount_sum /= (ndays/7)

    # prepare data for plotting
    data = zeros(7, nrow(categories))
    for (i, c) in enumerate(categories.category)
        for d in 1:1:7
            cut = (day_exp.category.==c) .& (day_exp.dayofweek.==d)
            if nrow(day_exp[cut, :]) != 0
                data[d, i] = day_exp[cut, "amount_sum"][1]
            end
        end
    end

    # plot histograms
    bins = [i for i in 0.5:1:7.5]
    lowers = bins[1:end-1]
    highers = bins[2:end]
    centres = 0.5*(lowers+highers)

    ncat = nrow(categories)
    width = 0.8 / ncat

    plot()
    for (i, c) in enumerate(categories.category)
        shift = (-(ncat/2)+(i-0.5))*width
        histogram!(centres.+shift, bins=bins.+shift, weights=data[:, i],
                   bar_width=width, label=c, color=COLORS[i], width=0.5)
    end

    # total spend
    plot!(centres, sum(data, dims=2), color="black", label="total")

    # cosmetics
    currency = expenses.currency[1]
    ylabel!("Daily spending ($currency)")
    days = Dict(1=>"Mon", 2=>"Tue", 3=>"Wed", 4=>"Thu", 5=>"Fri", 6=>"Sat", 7=>"Sun")
    xticks!(centres, [days[d] for d in centres])
    savefig("figures/weekly.pdf")

end

max_stars_repo_name: mzgubic/exchequer

max_issues_repo_name: mzgubic/exchequer

max_forks_repo_name: mzgubic/exchequer


text:
#!/usr/bin/env julia

function packagelt(a::String, b::String)
    a == "julia" && b != "julia" && return true
    a != "julia" && b == "julia" && return false
    return lowercase(a) < lowercase(b)
end

≲(v::VersionNumber, t::NTuple{0,Int}) = true
≲(v::VersionNumber, t::NTuple{1,Int}) = v.major ≤ t[1]
≲(v::VersionNumber, t::NTuple{2,Int}) = v.major < t[1] ||
                                        v.major ≤ t[1] && v.minor ≤ t[2]
≲(v::VersionNumber, t::NTuple{3,Int}) = v.major < t[1] ||
                                        v.major ≤ t[1] && v.minor < t[2] ||
                                        v.major ≤ t[1] && v.minor ≤ t[2] && v.patch ≤ t[3]

≲(t::NTuple{0,Int}, v::VersionNumber) = true
≲(t::NTuple{1,Int}, v::VersionNumber) = t[1] ≤ v.major
≲(t::NTuple{2,Int}, v::VersionNumber) = t[1] < v.major ||
                                        t[1] ≤ v.major && t[2] ≤ v.minor
≲(t::NTuple{3,Int}, v::VersionNumber) = t[1] < v.major ||
                                        t[1] ≤ v.major && t[2] < v.minor ||
                                        t[1] ≤ v.major && t[2] ≤ v.minor && t[3] ≤ v.patch

function compress_versions(inc::Vector{VersionNumber}, from::Vector{VersionNumber})
    issorted(inc) || (inc = sort(inc))
    exc = sort!(setdiff(from, inc))
    pairs = []
    if isempty(exc)
        lo, hi = first(inc), last(inc)
        push!(pairs, (lo.major, lo.minor) => (hi.major, hi.minor))
    else
        for v in inc
            t = (v.major, v.minor)
            if any(t ≲ w ≲ t for w in exc)
                t = (v.major, v.minor, v.patch)
            end
            if isempty(pairs) || any(pairs[end][1] ≲ w ≲ t for w in exc)
                push!(pairs, t => t) # need a new interval
            else
                pairs[end] = pairs[end][1] => t # can be merged with last
            end
        end
    end
    @assert all(any(p[1] ≲ v ≲ p[2] for p ∈ pairs) for v ∈ inc)
    @assert all(!any(p[1] ≲ v ≲ p[2] for p ∈ pairs) for v ∈ exc)
    return pairs
end
compress_versions(f::Function, from::Vector{VersionNumber}) =
    compress_versions(filter(f, from), from)
compress_versions(vi::VersionInterval, from::Vector{VersionNumber}) =
    compress_versions(v->v in vi, from)
compress_versions(inc, from) = compress_versions(inc, collect(from))

versions_string(p::Pair) = versions_string(p...)
versions_string(a::Tuple{}, b::Tuple{}) = "*"
versions_string(a::NTuple{m,Int}, b::NTuple{n,Int}) where {m,n} =
    a == b ? join(a, '.') : "$(join(a, '.'))-$(join(b, '.'))"

versions_repr(x) = repr(versions_string(x))
versions_repr(v::Vector) = length(v) == 1 ? repr(versions_string(v[1])) :
    "[" * join(map(repr∘versions_string, v), ", ") * "]"

## compress per-version data ##

function compress_versions_data(
    d::Dict{VersionNumber,Dict{String,String}},
    versions::Vector{VersionNumber},
)
    kvs = Dict{Pair{String,String},Vector{VersionNumber}}()
    for (ver, x) in d, (k, v) in x
        push!(get!(kvs, k => v, VersionNumber[]), ver)
    end
    tx = Tuple{Any,String,String}[]
    for (kv, vs) in kvs
        for v in compress_versions(sort!(vs), versions)
            push!(tx, (v, kv...))
        end
    end
    sort!(tx, by=t->(t[1][1], t[1][2], t[2], t[3]))
    return Tuple{String,String,String}[
        (versions_string(t[1]), t[2], t[3]) for t in tx
    ]
end

## dict utility functions ##

function invert_map(fwd::Dict{K,V}) where {K,V}
    rev = Dict{V,Vector{K}}()
    for (k, v) in fwd
        push!(get!(rev, v, K[]), k)
    end
    return rev
end

function invert_map(fwd::Dict{Vector{K},V}) where {K,V}
    rev = Dict{V,Vector{K}}()
    for (k, v) in fwd
        append!(get!(rev, v, K[]), k)
    end
    return rev
end

flatten_keys(d::Dict{Vector{K},V}) where {K,V} =
    isempty(d) ? Dict{K,V}() : Dict{K,V}(k => v for (ks, v) in d for k in ks)


text:
module Fast_Hertz

using JuLIP
using JuLIP.Potentials: pairs
import JuLIP: energy, forces, cutoff
using NeighbourLists
using StaticArrays

const SVecG{T} = SVector{3, T}
const SMatG{T} = SMatrix{3, 3, T, 9}

export Fast_Hertz_Potential

@pot struct Fast_HP <: AbstractCalculator
    rcut::Float64    # cutoff when calculating neighbourlist
    rtol2::Float64   # if max_displacement > rtol, recalculate the neighbourlist rtol2 = rtol^2
    rsum::Vector{Float64} # rsum = [2*rmin; rmin + rmax; 2*rmax]
    Lhalf::Float64
end

function Fast_Hertz_Potential(L, r0::Float64, r1::Float64; rcutfact = 2.5)
    rcut = rcutfact * r0 # rcut is rcutfact*2*R_{max}, used when searching possible pairs
    rtol = 0.5 * rcut    # rtol is the neighbourlist recomputing tolerance
    rtol2 = rtol * rtol
    rsum = [r1; (r0 + r1)/2; r0]
    Lhalf = L/2.0
    return Fast_HP(rcut, rtol2, rsum, Lhalf)
end

function energy(V::Fast_HP, at::AbstractAtoms)
    nlist_i, nlist_j, nlist_Z = update_neighbourlist!(V, at)
    bins = at.cell'
    inv_bins = inv(bins)
    return energy_inner(nlist_i, nlist_j, nlist_Z, at.X, V, bins, inv_bins)
end
# compute the maximum displacement of all particles
@inline max_dr2(X::Vector{JVecF}, X_Old::Vector{JVecF}) =
                    maximum( [ sum( (X[i] - X_Old[i]) .* (X[i] - X_Old[i]) ) for i = 1:length(X) ] )

@inline PairType(i::Vector{Int}, j::Vector{Int}, Z::Vector{Int}) = [ Z[i[k]] + Z[j[k]] - 1 for k = 1:length(i) ]

function update_neighbourlist!(V::Fast_HP, at::AbstractAtoms)
    if has_data(at, :nlist_X)
        X_Old = get_data(at, :nlist_X)::Vector{JVecF}

        X = at.X::Vector{JVecF}
        max_d = max_dr2(X, X_Old)
        if max_d < V.rtol2
            nlist_i = get_data(at, :nlist_i)::Vector{Int}
            nlist_j = get_data(at, :nlist_j)::Vector{Int}
            nlist_Z = get_data(at, :nlist_Z)::Vector{Int}
            return nlist_i, nlist_j, nlist_Z
        end
    end
    nlist = neighbourlist(at, V.rcut)
    nlist_Z = PairType(nlist.i, nlist.j, at.Z)
    set_data!(at, :nlist_X, copy(at.X))
    set_data!(at, :nlist_i, nlist.i)
    set_data!(at, :nlist_j, nlist.j)
    set_data!(at, :nlist_Z, nlist_Z)
    return nlist.i, nlist.j, nlist_Z
end

@inline wrap_to_unit(inv_bins::SMatG{T}, x::SVecG{T}) where {T <: Real} = mod.(inv_bins * x, 1)
@inline adjust_dx(bins::SMatG{T}, dx::SVecG{T}) where {T <: Real} = bins * (dx - (dx .> 0.5))
@inline dx2dx(dx::SVecG{T}, inv_bins::SMatG{T}, bins::SMatG{T}) where {T <: Real} =
                                    adjust_dx(bins, wrap_to_unit(inv_bins, dx))

function energy_inner(nlist_i, nlist_j, nlist_Z, X, V::Fast_HP, bins, inv_bins)
    E = 0.0
    for k = 1:length(nlist_i)
        dx = X[nlist_j[k]] - X[nlist_i[k]]
        r = norm(dx)
        if r > V.Lhalf
            dx = dx2dx(dx, inv_bins, bins)
            r = norm(dx)
        end
        nlist_s = r - V.rsum[nlist_Z[k]]
        E += 1.0 * nlist_s * nlist_s * (nlist_s < 0)
    end
    return 0.5*E
end

end


text:
__precompile__()
module FMMLIB2D

using Compat

export lfmm2d, rfmm2d, hfmm2d, zfmm2d
export lfmm2dparttarg, lfmm2dpartself
export rfmm2dparttarg, rfmm2dpartself
export hfmm2dparttarg
export zfmm2dparttarg

const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("FMMLIB2D is not properly installed. Please run Pkg.build(\"FMMLIB2D\").")
end

struct FMMOutput
    pot
    grad
    hess
    pottarg
    gradtarg
    hesstarg
end

## KEYWORD INTERFACES

"""
    U = lfmm2d(;source, target, charge, dipstr, dipvec, tol,
                ifpot, ifgrad, ifhess,
                ifpottarg, ifgradtarg, ifhesstarg)

Laplace particle target FMM in R^2 (complex). Keyword interface.

Returns: U::FMMLIB2D.FMMOutput
"""
function lfmm2d(;
    source::Array{Float64} = zeros(2, 0),
    target::Array{Float64} = zeros(2, 0),
    charge::Array{ComplexF64} = zeros(ComplexF64, 0),
    dipstr::Array{ComplexF64} = zeros(ComplexF64, 0),
    dipvec::Array{Float64} = zeros(2, 0),
    tol::Float64 = 1e-15,
    ifpot::Bool = true,
    ifgrad::Bool = false,
    ifhess::Bool = false,
    ifpottarg::Bool = true,
    ifgradtarg::Bool = false,
    ifhesstarg::Bool = false
                )
    # Parse keywords
    iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
    target, ifpottarg, ifgradtarg, ifhesstarg =
        parse_keywords(tol, source, charge, dipstr, dipvec, ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    # Call FMM and pack output
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        lfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot,
                       ifgrad, ifhess, target, ifpottarg, ifgradtarg, ifhesstarg)
    return FMMOutput(pot, grad, hess, pottarg, gradtarg, hesstarg)
end

"""
    U = rfmm2d(;source, target, charge, dipstr, dipvec, tol,
                ifpot, ifgrad, ifhess,
                ifpottarg, ifgradtarg, ifhesstarg)

Laplace particle target FMM in R^2 (real). Keyword interface.

Returns: U::FMMLIB2D.FMMOutput
"""
function rfmm2d(;
    source::Array{Float64} = zeros(2, 0),
    target::Array{Float64} = zeros(2, 0),
    charge::Array{Float64} = zeros(Float64, 0),
    dipstr::Array{Float64} = zeros(Float64, 0),
    dipvec::Array{Float64} = zeros(2, 0),
    tol::Float64 = 1e-15,
    ifpot::Bool = true,
    ifgrad::Bool = false,
    ifhess::Bool = false,
    ifpottarg::Bool = true,
    ifgradtarg::Bool = false,
    ifhesstarg::Bool = false
                )
    # Parse keywords
    iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
    target, ifpottarg, ifgradtarg, ifhesstarg =
        parse_keywords(tol, source, charge, dipstr, dipvec, ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    # Call FMM and pack output
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        rfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot,
                       ifgrad, ifhess, target, ifpottarg, ifgradtarg, ifhesstarg)
    return FMMOutput(pot, grad, hess, pottarg, gradtarg, hesstarg)
end

"""
    U = hfmm2d(;source, target, charge, dipstr, dipvec,
                tol, zk,
                ifpot, ifgrad, ifhess,
                ifpottarg, ifgradtarg, ifhesstarg)

Helmholtz particle target FMM in R^2. Keyword interface.

Returns: U::FMMLIB2D.FMMOutput
"""
function hfmm2d(;
                source::Array{Float64} = zeros(2, 0),
                target::Array{Float64} = zeros(2, 0),
                charge::Array{ComplexF64} = zeros(ComplexF64, 0),
                dipstr::Array{ComplexF64} = zeros(ComplexF64, 0),
                dipvec::Array{Float64} = zeros(2, 0),
                tol::Float64 = 1e-15,
                zk::ComplexF64 = complex(0.0),
                ifpot::Bool = true,
                ifgrad::Bool = false,
                ifhess::Bool = false,
                ifpottarg::Bool = true,
                ifgradtarg::Bool = false,
                ifhesstarg::Bool = false
                )
    @assert zk != 0
    # Parse keywords
    iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
    target, ifpottarg, ifgradtarg, ifhesstarg =
        parse_keywords(tol, source, charge, dipstr, dipvec, ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    # Call FMM and pack output
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        hfmm2dparttarg(iprec, zk, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot,
                       ifgrad, ifhess, target, ifpottarg, ifgradtarg, ifhesstarg)
    return FMMOutput(pot, grad, hess, pottarg, gradtarg, hesstarg)
end

"""
    U = zfmm2d(;source, target, dipstr, tol,
                ifpot, ifgrad, ifhess,
                ifpottarg, ifgradtarg, ifhesstarg)

Complex particle target FMM in R^2. Keyword interface.

Returns: U::FMMLIB2D.FMMOutput
"""
function zfmm2d(;
                source::Array{ComplexF64} = complex(zeros(0)),
                target::Array{ComplexF64} = complex(zeros(0)),
                dipstr::Array{ComplexF64} = complex(zeros(0)),
                tol::Float64 = 1e-15,
                ifpot::Bool = true,
                ifgrad::Bool = false,
                ifhess::Bool = false,
                ifpottarg::Bool = true,
                ifgradtarg::Bool = false,
                ifhesstarg::Bool = false
                )
    # Parse keywords
    iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
    target, ifpottarg, ifgradtarg, ifhesstarg =
        parse_keywords(tol, source, [], dipstr, [], ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    # Call FMM and pack output
    pot, grad, hess, pottarg, gradtarg, hesstarg = zfmm2dparttarg(iprec, source, dipstr,
                                                                  ifpot, ifgrad, ifhess,
                                                                  target, ifpottarg,
                                                                  ifgradtarg, ifhesstarg)
    return FMMOutput(pot, grad, hess, pottarg, gradtarg, hesstarg)
end

## DIRECT INTERFACES

"""
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        lfmm2dparttarg( iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec,
                        ifpot, ifgrad, ifhess,
                        target, ifpottarg, ifgradtarg, ifhesstarg)

    Laplace particle target FMM in R^2 (complex). Direct library interface.
"""
function lfmm2dparttarg(iprec::Int64,
                        source::Array{Float64},
                        ifcharge::Int64,
                        charge::Array{ComplexF64},
                        ifdipole::Int64,
                        dipstr::Array{ComplexF64},
                        dipvec::Array{Float64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64,
                        target::Array{Float64},
                        ifpottarg::Int64,
                        ifgradtarg::Int64,
                        ifhesstarg::Int64)
    # Size checks
    nsource = size(source, 2)
    @assert size(source, 1) == 2
    if ifcharge != 0
        @assert length(charge) == nsource
    end
    if ifdipole != 0
        @assert length(dipstr) == nsource
        @assert size(dipvec) == (2, nsource)
    end
    ntarget = size(target, 2)
    @assert size(source, 1) == 2
    # Prepare output structures, only allocate those needed
    ier = 0
    pot = zeros(ComplexF64, nsource * (ifpot!=0 ? 1 : 0) )
    grad = zeros(ComplexF64, 2, nsource * (ifgrad!=0 ? 1 : 0) )
    hess = zeros(ComplexF64, 3, nsource * (ifhess!=0 ? 1 : 0) )
    pottarg = Array{ComplexF64}(undef, ntarget * (ifpottarg!=0 ? 1 : 0) )
    gradtarg = Array{ComplexF64}(undef, 2, ntarget * (ifgradtarg!=0 ? 1 : 0) )
    hesstarg = Array{ComplexF64}(undef, 3, ntarget * (ifhesstarg!=0 ? 1 : 0) )
    # Library call
    ccall( (:lfmm2dparttarg_, fmmlib2d), Nothing,
           (Ref{Int64}, # ier
            Ref{Int64}, # iprec
            Ref{Int64}, # nsource
            Ref{Float64}, # source
            Ref{Int64}, # ifcharge
            Ref{ComplexF64}, #charge
            Ref{Int64}, # ifdipole
            Ref{ComplexF64}, #dipstr
            Ref{Float64}, #dipvec
            Ref{Int64}, # ifpot
            Ref{ComplexF64}, # pot
            Ref{Int64}, # ifgrad
            Ref{ComplexF64}, #grad
            Ref{Int64}, # ifhess
            Ref{ComplexF64}, #hess
            Ref{Int64}, # ntarget
            Ref{Float64}, # target
            Ref{Int64}, # ifpottarg
            Ref{ComplexF64}, # pottarg
            Ref{Int64}, # ifgradtarg
            Ref{ComplexF64}, # gradtarg
            Ref{Int64}, # ifhesstarg
            Ref{ComplexF64}, # hesstarg
            ),
           ier,iprec,nsource,source,
           ifcharge,charge,ifdipole,dipstr,dipvec,
           ifpot,pot,ifgrad,grad,ifhess,hess,
           ntarget,target,ifpottarg,pottarg,ifgradtarg,gradtarg,ifhesstarg,hesstarg
           )
    # Check and return
    @assert ier == 0
    return pot, grad, hess, pottarg, gradtarg, hesstarg
end

"""
    pot, grad, hess =
        lfmm2dpartself( iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec,
                        ifpot, ifgrad, ifhess)

Laplace particle FMM in R^2 (complex). Direct library interface.
"""
function lfmm2dpartself(iprec::Int64,
                        source::Array{Float64},
                        ifcharge::Int64,
                        charge::Array{ComplexF64},
                        ifdipole::Int64,
                        dipstr::Array{ComplexF64},
                        dipvec::Array{Float64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64)
    target = zeros(0,0)
    ifpottarg = 0
    ifgradtarg = 0
    ifhesstarg = 0
    pot, grad, hess, _, _, _ =
        lfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    return pot, grad, hess
end

"""
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        rfmm2dparttarg( iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec,
                        ifpot, ifgrad, ifhess,
                        target, ifpottarg, ifgradtarg, ifhesstarg)

Laplace particle target FMM in R^2 (real). Direct library interface.
"""
function rfmm2dparttarg(iprec::Int64,
                        source::Array{Float64},
                        ifcharge::Int64,
                        charge::Array{Float64},
                        ifdipole::Int64,
                        dipstr::Array{Float64},
                        dipvec::Array{Float64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64,
                        target::Array{Float64},
                        ifpottarg::Int64,
                        ifgradtarg::Int64,
                        ifhesstarg::Int64)
    # Size checks
    nsource = size(source, 2)
    @assert size(source, 1) == 2
    if ifcharge != 0
        @assert length(charge) == nsource
    end
    if ifdipole != 0
        @assert length(dipstr) == nsource
        @assert size(dipvec) == (2, nsource)
    end
    ntarget = size(target, 2)
    @assert size(source, 1) == 2
    # Prepare output structures, only allocate those needed
    ier = 0
    pot = zeros(Float64, nsource * (ifpot!=0 ? 1 : 0) )
    grad = zeros(Float64, 2, nsource * (ifgrad!=0 ? 1 : 0) )
    hess = zeros(Float64, 3, nsource * (ifhess!=0 ? 1 : 0) )
    pottarg = Array{Float64}(undef, ntarget * (ifpottarg!=0 ? 1 : 0) )
    gradtarg = Array{Float64}(undef, 2, ntarget * (ifgradtarg!=0 ? 1 : 0) )
    hesstarg = Array{Float64}(undef, 3, ntarget * (ifhesstarg!=0 ? 1 : 0) )
    # Library call
    ccall( (:rfmm2dparttarg_, fmmlib2d), Nothing,
           (Ref{Int64}, # ier
            Ref{Int64}, # iprec
            Ref{Int64}, # nsource
            Ref{Float64}, # source
            Ref{Int64}, # ifcharge
            Ref{Float64}, #charge
            Ref{Int64}, # ifdipole
            Ref{Float64}, #dipstr
            Ref{Float64}, #dipvec
            Ref{Int64}, # ifpot
            Ref{Float64}, # pot
            Ref{Int64}, # ifgrad
            Ref{Float64}, #grad
            Ref{Int64}, # ifhess
            Ref{Float64}, #hess
            Ref{Int64}, # ntarget
            Ref{Float64}, # target
            Ref{Int64}, # ifpottarg
            Ref{Float64}, # pottarg
            Ref{Int64}, # ifgradtarg
            Ref{Float64}, # gradtarg
            Ref{Int64}, # ifhesstarg
            Ref{Float64}, # hesstarg
            ),
           ier,iprec,nsource,source,
           ifcharge,charge,ifdipole,dipstr,dipvec,
           ifpot,pot,ifgrad,grad,ifhess,hess,
           ntarget,target,ifpottarg,pottarg,ifgradtarg,gradtarg,ifhesstarg,hesstarg
           )
    # Check and return
    @assert ier == 0
    return pot, grad, hess, pottarg, gradtarg, hesstarg
end

"""
    pot, grad, hess =
        rfmm2dpartself( iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec,
                        ifpot, ifgrad, ifhess)

Laplace particle FMM in R^2 (real). Direct library interface.
"""
function rfmm2dpartself(iprec::Int64,
                        source::Array{Float64},
                        ifcharge::Int64,
                        charge::Array{Float64},
                        ifdipole::Int64,
                        dipstr::Array{Float64},
                        dipvec::Array{Float64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64)
    target = zeros(0,0)
    ifpottarg = 0
    ifgradtarg = 0
    ifhesstarg = 0
    pot, grad, hess, _, _, _ =
        rfmm2dparttarg(iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad, ifhess,
                       target, ifpottarg, ifgradtarg, ifhesstarg)
    return pot, grad, hess
end

"""
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        hfmm2dparttarg( iprec, zk, source, ifcharge, charge,
                        ifdipole, dipstr, dipvec,
                        ifpot, ifgrad, ifhess,
                        target, ifpottarg, ifgradtarg, ifhesstarg)

Helmholtz particle target FMM in R^2. Direct library interface.
"""
function hfmm2dparttarg(iprec::Int64,
                        zk::ComplexF64,
                        source::Array{Float64},
                        ifcharge::Int64,
                        charge::Array{ComplexF64},
                        ifdipole::Int64,
                        dipstr::Array{ComplexF64},
                        dipvec::Array{Float64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64,
                        target::Array{Float64},
                        ifpottarg::Int64,
                        ifgradtarg::Int64,
                        ifhesstarg::Int64)
    # Size checks
    nsource = size(source, 2)
    @assert size(source, 1) == 2
    if ifcharge != 0
        @assert length(charge) == nsource
    end
    if ifdipole != 0
        @assert length(dipstr) == nsource
        @assert size(dipvec) == (2, nsource)
    end
    ntarget = size(target, 2)
    @assert size(source, 1) == 2
    # Prepare output structures, only allocate those needed
    ier = 0
    pot = zeros(ComplexF64, nsource * (ifpot!=0 ? 1 : 0) )
    grad = zeros(ComplexF64, 2, nsource * (ifgrad!=0 ? 1 : 0) )
    hess = zeros(ComplexF64, 3, nsource * (ifhess!=0 ? 1 : 0) )
    pottarg = Array{ComplexF64}(undef, ntarget * (ifpottarg!=0 ? 1 : 0) )
    gradtarg = Array{ComplexF64}(undef, 2, ntarget * (ifgradtarg!=0 ? 1 : 0) )
    hesstarg = Array{ComplexF64}(undef, 3, ntarget * (ifhesstarg!=0 ? 1 : 0) )
    # Library call
    ccall( (:hfmm2dparttarg_, fmmlib2d), Nothing,
           (Ref{Int64}, # ier
            Ref{Int64}, # iprec
            Ref{ComplexF64}, # zk
            Ref{Int64}, # nsource
            Ref{Float64}, # source
            Ref{Int64}, # ifcharge
            Ref{ComplexF64}, #charge
            Ref{Int64}, # ifdipole
            Ref{ComplexF64}, #dipstr
            Ref{Float64}, #dipvec
            Ref{Int64}, # ifpot
            Ref{ComplexF64}, # pot
            Ref{Int64}, # ifgrad
            Ref{ComplexF64}, #grad
            Ref{Int64}, # ifhess
            Ref{ComplexF64}, #hess
            Ref{Int64}, # ntarget
            Ref{Float64}, # target
            Ref{Int64}, # ifpottarg
            Ref{ComplexF64}, # pottarg
            Ref{Int64}, # ifgradtarg
            Ref{ComplexF64}, # gradtarg
            Ref{Int64}, # ifhesstarg
            Ref{ComplexF64}, # hesstarg
            ),
           ier,iprec,zk,nsource,source,
           ifcharge,charge,ifdipole,dipstr,dipvec,
           ifpot,pot,ifgrad,grad,ifhess,hess,
           ntarget,target,ifpottarg,pottarg,ifgradtarg,gradtarg,ifhesstarg,hesstarg
           )
    # Check and return
    @assert ier == 0
    return pot, grad, hess, pottarg, gradtarg, hesstarg
end

"""
    pot, grad, hess, pottarg, gradtarg, hesstarg =
        zfmm2dparttarg(iprec,
                        source,
                        dipstr,
                        ifpot,
                        ifgrad,
                        ifhess,
                        target,
                        ifpottarg,
                        ifgradtarg,
                        ifhesstarg)

Complex particle target FMM in R^2. Direct library interface.
"""
function zfmm2dparttarg(iprec::Int64,
                        source::Array{ComplexF64},
                        dipstr::Array{ComplexF64},
                        ifpot::Int64,
                        ifgrad::Int64,
                        ifhess::Int64,
                        target::Array{ComplexF64},
                        ifpottarg::Int64,
                        ifgradtarg::Int64,
                        ifhesstarg::Int64)
    # Size checks
    nsource = length(source)
    @assert length(dipstr) == nsource
    ntarget = length(target)
    # Prepare output structures, only allocate those needed
    ier = 0
    pot  = zeros(ComplexF64, nsource * (ifpot !=0 ? 1 : 0) )
    grad = zeros(ComplexF64, nsource * (ifgrad!=0 ? 1 : 0) )
    hess = zeros(ComplexF64, nsource * (ifhess!=0 ? 1 : 0) )
    pottarg  = zeros(ComplexF64, ntarget * (ifpottarg !=0 ? 1 : 0) )
    gradtarg = zeros(ComplexF64, ntarget * (ifgradtarg!=0 ? 1 : 0) )
    hesstarg = zeros(ComplexF64, ntarget * (ifhesstarg!=0 ? 1 : 0) )
    # Library call
    ccall( (:zfmm2dparttarg_, fmmlib2d), Nothing,
           (Ref{Int64}, # ier
            Ref{Int64}, # iprec
            Ref{Int64}, # nsource
            Ref{ComplexF64}, # source
            Ref{ComplexF64}, #dipstr
            Ref{Int64}, # ifpot
            Ref{ComplexF64}, # pot
            Ref{Int64}, # ifgrad
            Ref{ComplexF64}, #grad
            Ref{Int64}, # ifhess
            Ref{ComplexF64}, #hess
            Ref{Int64}, # ntarget
            Ref{ComplexF64}, # target
            Ref{Int64}, # ifpottarg
            Ref{ComplexF64}, # pottarg
            Ref{Int64}, # ifgradtarg
            Ref{ComplexF64}, # gradtarg
            Ref{Int64}, # ifhesstarg
            Ref{ComplexF64}, # hesstarg
            ),
           ier, iprec, nsource, source, dipstr,
           ifpot, pot, ifgrad, grad, ifhess, hess,
           ntarget, target, ifpottarg, pottarg, ifgradtarg, gradtarg, ifhesstarg, hesstarg
           )
    # Check and return
    @assert ier == 0
    return pot, grad, hess, pottarg, gradtarg, hesstarg
end

## HELPER FUNCTIONS

function parse_keywords(tol, source, charge, dipstr, dipvec, ifpot, ifgrad, ifhess,
                        target, ifpottarg, ifgradtarg, ifhesstarg)
    # Set correct integer flags
    if length(target)==0
        ifpottarg = false
        ifgradtarg = false
        ifhesstarg = false
    end
    ifcharge = length(charge)>0 ? 1 : 0
    ifdipole = length(dipstr)>0 ? 1 : 0
    ifpot = ifpot ? 1 : 0
    ifgrad = ifgrad ? 1 : 0
    ifhess = ifhess ? 1 : 0
    ifpottarg = ifpottarg ? 1 : 0
    ifgradtarg = ifgradtarg ? 1 : 0
    ifhesstarg = ifhesstarg ? 1 : 0
    # Select the right precision flag
    iprec = tol2iprec(tol)
    return iprec, source, ifcharge, charge, ifdipole, dipstr, dipvec, ifpot, ifgrad,
    ifhess, target, ifpottarg, ifgradtarg, ifhesstarg
end

function tol2iprec(tol)
    if tol < 0.5e-14
        iprec = 5
    elseif tol < 0.5e-12
        iprec = 4
    elseif tol < 0.5e-9
        iprec = 3
    elseif tol < 0.5e-6
        iprec = 2
    elseif tol < 0.5e-3
        iprec = 1
    elseif tol < 0.5e-2
        iprec =
            0
    elseif tol < 0.5e-1
        iprec = -1
    else
        iprec = -2
    end
end

end # module


text:
"""
This function updates the expected profit of the project from all products.
Returns nothing.
"""
function update_expected_profit!(project::P,
                                 scenario_data::Vector{Scenario},
                                 market_prices::MarketPrices,
                                 price_years::NamedTuple{(:start_year, :end_year),
                                              Tuple{Int64, Int64}},
                                 update_years::NamedTuple{(:start_year, :end_year),
                                              Tuple{Int64, Int64}},
                                 rep_hour_weight::Vector{Float64},
                                 queue_cost::Vector{Float64},
                                 capacity_forward_years::Int64,
                                 iteration_year::Int64,
                                 solver::JuMP.MOI.OptimizerWithAttributes) where P <: Project{<: BuildPhase}

    extend_profit_array(project, get_name.(scenario_data), iteration_year)

    profit_array_length = price_years[:end_year] - price_years[:start_year] + 1

    operating_profit_array = AxisArrays.AxisArray{Float64, 2}[] # Operating profit array for each scenario
    capacity_revenue_array = Vector{Float64}[]  # Capacity revenue array for each scenario
    rec_revenue_array = Vector{Float64}[]  # REC revenue array for each scenario

    expected_energy_production = zeros(profit_array_length)

    capacity_prices = get_capacity_price(market_prices)
    rec_prices = get_rec_price(market_prices)

    sz = price_years[:end_year] - price_years[:start_year] + 1

    total_annual_cost = calculate_annualized_costs(project,
                                                   sz,
                                                   queue_cost)

    for scenario in scenario_data
        operating_profit,
        energy_production = calculate_operating_profit(project,
                                                      get_name(scenario),
                                                      market_prices,
                                                      price_years[:start_year],
                                                      price_years[:end_year],
                                                      rep_hour_weight,
                                                      solver)
        push!(operating_profit_array, operating_profit)

        # Calculate expected energy production
        expected_energy_production += get_probability(scenario) * energy_production

        capacity_revenues = zeros(profit_array_length)
        rec_revenues = zeros(profit_array_length)

        for product in get_products(project)
            if !isnothing(capacity_prices)
                capacity_revenues = update_capacity_revenues!(
                                                              product,
                                                              get_name(scenario),
                                                              price_years,
                                                              capacity_revenues,
                                                              capacity_prices,
                                                              get_maxcap(project)
                                                            )
            end

            if !isnothing(rec_prices)
                rec_revenues = update_rec_revenues!(
                                                    product,
                                                    price_years,
                                                    rec_revenues,
                                                    rec_prices[get_name(scenario)],
                                                    energy_production
                                                    )
            end
        end

        annual_operating_profit =  [sum(operating_profit[:, y])
                               for y in price_years[:start_year]:price_years[:end_year]]

        capacity_going_forward_cost = total_annual_cost - annual_operating_profit - rec_revenues

        for product in get_products(project)
            update_profit!(project,
                           get_name(scenario),
                           product,
                           operating_profit,
                           capacity_revenues,
                           capacity_forward_years,
                           rec_revenues,
                           price_years,
                           update_years,
                           iteration_year)
        end

        push!(capacity_revenue_array, capacity_revenues)
        push!(rec_revenue_array, rec_revenues)
    end

    profit_axis_values = AxisArrays.axisvalues(operating_profit_array[1])

    expected_operating_profit = operating_profit_array[1]

    for x in profit_axis_values[1]
        for y in profit_axis_values[2]
            expected_operating_profit[x, y] = sum(get_probability(scenario_data[i]) * operating_profit_array[i][x, y] for i in 1:length(operating_profit_array))
        end
    end

    # Calculate expected profits
    expected_annual_operating_profit = [sum(expected_operating_profit[x, y] for x in profit_axis_values[1]) for y in profit_axis_values[2]]

    expected_capacity_revenues = [sum(get_probability(scenario_data[i]) * capacity_revenue_array[i][y] for i in 1:length(capacity_revenue_array))
                                 for y in 1:profit_array_length]

    expected_rec_revenues = [sum(get_probability(scenario_data[i]) * rec_revenue_array[i][y] for i in 1:length(rec_revenue_array))
                                 for y in 1:profit_array_length]

    #Calculate capacity market going forward cost.
    expected_capacity_going_forward_cost = total_annual_cost - expected_annual_operating_profit - expected_rec_revenues
    capacity_market_bid = 0.0

    capacity_revenue_start_year = max(1, capacity_forward_years - price_years[:start_year] + 1)

    if capacity_revenue_start_year <= length(expected_capacity_going_forward_cost)
        capacity_market_bid = round(max(0.0, expected_capacity_going_forward_cost[capacity_revenue_start_year] / get_maxcap(project)), digits = 4)
    end

    #Calculate REC market bid.
    rec_market_bid = max(0.0, (total_annual_cost[1] - expected_annual_operating_profit[1] - expected_capacity_revenues[1]) / max(1e-3, expected_energy_production[1]))

    for product in get_products(project)
        update_bid!(product,
                    capacity_market_bid,
                    rec_market_bid,
                    expected_energy_production[1]
                    )
    end

    # Sets the project profit values beyond the horizon equal to the last value.
    if get_end_life_year(project) > update_years[:end_year]
        for scenario_name in get_name.(scenario_data)
            for product_name in get_name.(get_products(project))
                for y in (update_years[:end_year] + 1):get_end_life_year(project)
                    profit = get_scenario_profit(get_finance_data(project))[scenario_name][iteration_year][Symbol(product_name),
                                                                                update_years[:end_year]]
                    set_scenario_profit!(get_finance_data(project), scenario_name, product_name, iteration_year, y, profit)
                end
            end
        end
    end

    return
end

"""
This function does nothing if the product is neither
an operating market or capacity or REC product.
The update_profit functions for each product is defined in files pertaning to those products.
"""
function update_profit!(project::P,
                        scenario_name::String,
                        product::T,
                        operating_profit::AxisArrays.AxisArray{Float64, 2},
                        capacity_revenues::Vector{Float64},
                        capacity_forward_years::Int64,
                        rec_revenues::Vector{Float64},
                        price_years::NamedTuple{(:start_year, :end_year),
                                    Tuple{Int64, Int64}},
                        update_years::NamedTuple{(:start_year, :end_year),
                                     Tuple{Int64, Int64}},
                        iteration_year::Int64) where {P <: Project{<: BuildPhase}, T <: Product}
    return
end

"""
This function does nothing if project is not an option type.
Returns nothing
"""
function extend_profit_array(project::P, scenario_names::Vector{String}, iteration_year::Int64) where P <: Project{<: BuildPhase}
    return
end

"""
This function extends the profit array of option projects
if their end life year exceeds existing length of the profit array.
"""
function extend_profit_array(project::P, scenario_names::Vector{String}, iteration_year::Int64) where P <: Project{Option}

    if get_end_life_year(project) > size(get_scenario_profit(get_finance_data(project))[scenario_names[1]][iteration_year], 2)
        for scenario_name in scenario_names
            profit = AxisArrays.AxisArray(zeros(length(get_products(project)), get_end_life_year(project)),
                                AxisArrays.Axis{:prod}(get_name.(get_products(project))),
                                AxisArrays.Axis{:year}(1:1:get_end_life_year(project)))
            set_scenario_profit!(get_finance_data(project), scenario_name, iteration_year, profit)
        end
    end
    return
end


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

using MLJ, PrettyPrinting
MLJ.color_off() # hide
@load KNNRegressor
# input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
# target
height = [178, 194, 165, 173, 168];

scitype(X.age)

pipe = @pipeline MyPipe(X -> coerce(X, :age=>Continuous),
                       hot = OneHotEncoder(),
                       knn = KNNRegressor(K=3),
                       target = UnivariateStandardizer());

pipe.knn.K = 2
pipe.hot.drop_last = true;

evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# functions only defined for windows

 include("ft_getdriverversion.jl")
 include("ft_getlibraryversion.jl")
 include("ft_getcomportnumber.jl")
 include("ft_resetport.jl")
 include("ft_cycleport.jl")
 include("ft_rescan.jl")
 include("ft_reload.jl")
 include("ft_setresetpiperetrycount.jl")
 include("ft_eeprom_read__ft_eeprom_program.jl")


text:
### chessboard

tile_size = 50
board_size = 4 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]

chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
img = chessboard

img = gen_chess();

###

f(x :: Array{<:Real, 2}) = println(x)

Array{Float64, 2} <: Array{Real, 2}

Float64 <: Real

function g()
    ar::Array{Real, 2} = ones(1,1)
end

g()

f(ones(1, 1))

###

using TestImages
using PyPlot

img = testimage("mandril_gray");
img = reverse(img, dims=1)
img = chessboard

point_count = 100;
spacing = 30;

coverage = 1;

(pi * spacing^2 * point_count)/length(img)

r = sqrt(( * ))

points, mag, mag_points = LAP_julia.gradient_points.find_keypoints_from_gradients(img, sigma = 5, number = point_count, spacing = spacing)

imgshow(mag_points)
imgshow(mag)

pos_x = [point.pos[1] for point in points]
pos_y = [point.pos[2] for point in points]

imgshow(img)
PyPlot.scatter(pos_y, pos_x, marker = :x)
gcf()

### probability
using PyPlot

img = testimage("mandril_gray");
img = chessboard

# all(vec(reshape(img, (1, :))) == img[:])

points, mag, mag_points = LAP_julia.gradient_points.find_keypoints_from_gradients_p_field(img, sigma = 1, number = 100, spacing = 30)

wsample(Array(1:2), [-0, 1], 12)

imgshow(mag_points)
imgshow(mag)

pos_x = [point.pos[1] for point in points]
pos_y = [point.pos[2] for point in points]

imgshow(img)
PyPlot.scatter(pos_y, pos_x, marker = :x)
gcf()

###
using LAP_julia

img = gen_chess()

# genrate flow
flow = gen_tiled_flow(size(img), 20, 1000);

showflow(flow)

# generate warpped image
imgw = LAP_julia.interpolation.warp_img(img, real(flow), imag(flow));

# params:
#lap
image_1 = img;
image_2 = imgw;
filter_num = 3;
filter_half_size = 32;
window_size = [65, 65];
#points
point_count = 100;
spacing = 10;

# get points
points, mag, mag_points = find_keypoints_from_gradients(img, sigma = 1, number = point_count, spacing = spacing);

using PyPlot
# see points
pos_x = [point.pos[1] for point in points]
pos_y = [point.pos[2] for point in points]

pos = transpose([pos_x pos_y])
pos[:, 1]

imgshow(img)
PyPlot.scatter(pos_y, pos_x, marker = :x)
gcf()

u_est, coeffs = single_lap(image_1, image_2, filter_half_size, window_size, 3, pos);

showflow(u_est)

### try imfilter with indices:
img = chessboard;

filter_size = 21
filter_half_size = (filter_size - 1) / 2
# Calculate separable filters from basis:
sigma = (filter_half_size + 2) / 4
gaus = ImageFiltering.KernelFactors.gaussian(sigma, filter_size)

using PyPlot
using Plots
# 1
kernel = kernelfactors((gaus, gaus))
kernp1 = broadcast(*, kernel...)
surface([-filter_half_size:filter_half_size], -filter_half_size:filter_half_size, kernp1)

@time imgfilt1 = ImageFiltering.imfilter(chessboard, kernel, "symmetric")
imgshow(imgfilt1)
surface(imgfilt1)

# @time imgfilt1 = ImageFiltering.imfilter(chessboard, kernel, "symmetric", inds=[1, 1])

imgfilt = similar(img)

inds = CartesianIndex(50, 50)

imfilter!(CPU1(Algorithm.FIRTiled()), imgfilt, img, kernel, "symmetric", inds)

imfilter!(CPU1(Algorithm.FIR()), imgfilt, img, (), NoPad(), inds)

imgshow(imgfilt)
surface(imgfilt1)

### exapmle
using ImageFiltering, Images, ComputationalResources

# Create a sample image:
tile_size = 50
board_size = 4 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]
chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
img = chessboard

# Create a filter:
sigma = 3
filter_size = 21
gaus = ImageFiltering.KernelFactors.gaussian(sigma, filter_size)
kernel = kernelfactors((gaus, gaus))

typeof(kernel) <: ImageFiltering.ProcessedKernel
typeof(kernel) <: Union{ImageFiltering.ArrayLike, ImageFiltering.Laplacian}

# params of imfilter function
imgfilt = copy(img)
# inds = CartesianIndex(50, 50)

# The filtering that errors:
# imfilter!(CPU1(Algorithm.FIRTiled()), imgfilt, img, kernel, "symmetric", (30:40, 30:40))
imfilter!(CPU1(Algorithm.FIRTiled()), imgfilt, img, kernel, Pad(:Symmetric), (50:50, 50:50))

imgshow(imgfilt)

### alternative to imfilter inds
function f(img_pix)
    sum(img_pix .* [0 1/10 0;
                    1/10 6/10 1/10;
                    0 1/10 0])
end

mapped_img = mapwindow(f, img, (3, 3), "symmetric", (2:8, 2:8))

imgshow(mapped_img)

###

flow = gen_tiled_flow(size(img), 30, 1000)

one = [1 + 0im, 2 - im]
two = [0 + 6im, 3 + im]

one = [1 + 0im]
two = [0 + im]

import LAP_julia: helpers.angle_rmse

angle_rmse(one, two)

rad2deg.(angle.(one).-angle.(two))

### mse

function mse2(x, y)
    @assert (size(x) == size(y)) "sizes of $x and $y don't match"
    return 1/length(x) * sum((x .- y).^2)
end

function mse1(x, y)
    @assert (size(x) == size(y)) "sizes of $x and $y don't match"
    return mean((x .- y).^2)
end

function mean(x)
    sum(x)/length(x)
end

ran = rand(123)
dan = rand(123)

mse1(ran, dan)
mse2(ran, dan)
mse1(ran, dan) ≈ mse2(ran, dan)

####


text:
# Volume properties of an allipsoid
# =================================
# densite defines the mass density;
# ellipsoid is defined by the centre x, radii a, b, c, and orientation R;
# the subroutine outpuss mass, volume, Euler tensor E, and inertia tensor J;
# __________________________________________________________________________
function ellipsoid_properties(a, b, c, R, density)

  volume = (4.0/3.0) * pi * a * b * c;
  mass = density * volume;

  J = [mass*(b*b + c*c)/5.0 0.0 0.0;
       0.0 mass*(a*a + c*c)/5.0 0.0;
       0.0 0.0 mass*(b*b + a*a)/5.0];

  # note that Inertia = Trace(Euler)*Identity - Euler,
  # hence Euler = 0.5 * Trace (Inertia)*Identity - Inertia,
  # as Trace(Inertia) = 3*Trace(Euler) - Trace(Euler)

  Id = [1. 0. 0.; 0. 1. 0.; 0. 0. 1.]
  E0 = 0.5*(J[1,1]+J[2,2]+J[3,3])*Id - J;

  E = R*E0*R';

  return mass, volume, E, J
end

max_stars_repo_name: tkoziara/edyna

max_issues_repo_name: tkoziara/edyna

max_forks_repo_name: tkoziara/edyna


text:
using StrideArrays
using Documenter

makedocs(;
    modules=[StrideArrays],
    authors="Chris Elrod",
    repo="https://github.com/chriselrod/StrideArrays.jl/blob/{commit}{path}#L{line}",
    sitename="StrideArrays.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://chriselrod.github.io/StrideArrays.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Architecture Benchmarks" => [
            "arches/cascadelake.md",
            "arches/tigerlake.md",
            "arches/haswell.md"
        ],
        "Random Number Generation" => "rng.md",
        "Broadcasting" => "broadcasting.md",
        "Stack Allocattion" => "stack_allocation.md"
    ],
    strict=false,
)

deploydocs(;
    repo="github.com/chriselrod/StrideArrays.jl",
)


text:
module STLCuttersTests

using Test

@testset "STLs" begin include("STLsTests.jl") end
@testset "SimplexFaces" begin include("SimplexFacesTests.jl") end
@testset "Polyhedra" begin include("PolyhedraTests.jl") end
@testset "SubTriangulations" begin include("SubTriangulationsTests.jl") end
@testset "Rotations" begin include("Rotations.jl") end
@testset "SampleGeometries" begin include("SampleGeometriesTests.jl") end

@testset "Integration" begin include("Integration.jl") end
@testset "Poisson" begin include("Poisson.jl") end
@testset "PoissonAgFEM" begin include("PoissonAgFEM.jl") end

include(joinpath(@__DIR__,"..","examples","runexamples.jl"))

end # module


text:
"""
    gen_initial_state(problem,parameters,information,options)

Generate an initial state, i.e., compute uniformly distributed random vectors in bounds,
after that are evaluated in objective function. This method require that `parameters.N`
is valid attribute.
"""
gen_initial_state(problem,parameters,information,options,status::State{Any}) = gen_initial_state(problem,parameters,information,options)

function gen_initial_state(problem,parameters,information,options, status)
    if parameters.N != length(status.population)
        error("Population size in provided State differs from that in parameters")
    end

    if size(problem.bounds,2) != length(get_position(status.best_sol))
        error("Invalid population (dimension does not match with bounds)")
    end

    return State(status.best_sol, status.population)

end

function gen_initial_state(problem,parameters,information,options)
    # population array
    population = generate_population(parameters.N, problem,ε=options.h_tol)

    # best solution
    best_solution = get_best(population)

    return State(best_solution, population; f_calls = length(population), iteration=1)
end

function Base.show(io::IO, parameters::AbstractParameters)
    s = typeof(parameters)

    vals = string.(map(f -> getfield(parameters, f), fieldnames(s)))
    str = string(s) * "(" * join(string.(fieldnames(s)) .* "=" .* vals, ", ") * ")"

    print(io, str)
end


text:
module DiscreteChoice

using StatsBase
using ForwardDiff, Optim

export llk_sample, ccp, opt_homo, opt_hetero,
    estDCModel!, initial_guess,
    DCModel, updateCoef!

# Write your package code here.

include("struct.jl")
include("util.jl")
include("main.jl")

end


text:
module Grisu

const SHORTEST = 1
const FIXED = 2
const PRECISION = 3

const BUFFER = Array(Uint8,309+17)

include("float.jl")
include("fastdtoa.jl")
include("fastfixedtoa.jl")
include("bignum-dtoa.jl")
#TODO
 #add Float16 support/tests
 #add grisu folder, rework grisu.jl in base
 #remove double-conversion from makefiles
 #what to do with tests? everything but mega files?

function grisu(v::FloatingPoint,mode,requested_digits,buffer=BUFFER)
    if signbit(v)
        neg = true
        v = -v
    else
        neg = false
    end
    if mode == PRECISION && requested_digits == 0
        buffer[1] = 0x00
        len = 0
        return 0, 0, neg, buffer
    end
    if v == 0.0
        buffer[1] = 0x30
        buffer[2] = 0x00
        len = point = 1
        return len, point, neg, buffer
    end
    if mode == SHORTEST
        status,len,point,buffer = fastdtoa(v,SHORTEST,0,buffer)
    elseif mode == FIXED
        status,len,point,buffer = fastfixedtoa(v,0,requested_digits,buffer)
    elseif mode == PRECISION
        status,len,point,buffer = fastdtoa(v,PRECISION,requested_digits,buffer)
    end
    status && return len, point, neg, buffer
    status, len, point, buffer = bignumdtoa(v,mode,requested_digits,buffer)
    return len, point, neg, buffer
end

end # module


text:
#!/usr/local/bin/julia

pushfirst!(LOAD_PATH, "/home/matt/Documents/3dp/3dprinting/")

using Mesh

n = Mesh.Net()

function slice!(n::Net, r1, r2, z1, z2, rstp, steps)
	r = r1
	for z in z1:z2
		vo = vertex!(n, 0.0, 0.0, z)
		for t = rstp:rstp:2pi
			vertex!(n, r * sin(t), r * cos(t), z)
		end

		face!(n, vo, vo+steps, vo+1)
		for v in vo+2:vo+steps
			face!(n, vo, v-1, v)
		end

		if z > z1
			face!(n, vo+1, vo+steps, vo-1)
			face!(n, vo+1, vo-steps, vo-1)
			for s in 1:steps-1
				face!(n, vo+s+1, vo+s, vo+s-steps-1)
				face!(n, vo+s-steps-1, vo+s-steps, vo+s+1 )
			end
		end
		r = r2
	end
end

function solid(n::Net)
	r = 30
	steps = 5
	rstp = 2pi / steps
	zstp = 21

	for z in 0.0:zstp:360.0
		r = 30 + 15 * sin(pi * 2z/360)
		vo = vertex!(n, 0.0, 0.0, z)
		ve = vo
		for t = rstp:rstp:2pi
			ve = vertex!(n, r * sin(t), r * cos(t), z)
		end
		if z < zstp || z > 300 - zstp
			face!(n, vo, vo+steps, vo+1)
			for v in vo+2:vo+steps
				face!(n, vo, v-1, v)
			end
		end

		if z > 0
			face!(n, vo+1, vo+steps, vo-1)
			face!(n, vo+1, vo-steps, vo-1)
			for s in 1:steps-1
				face!(n, vo+s+1, vo+s, vo+s-steps-1)
				face!(n, vo+s-steps-1, vo+s-steps, vo+s+1 )
			end
		end
	end

	STL_ASCII(n, "circ.stl")
end

function hollow(n::Net)
	r = 30
	steps = 14
	rstp = 2pi / steps
	zstp = 21
	thk = 5
	vso = veo = vsi = vei = 0

	slice!(n, 30, 30, -thk, 0, rstp, steps)

	z = 0

	for z in 0:zstp:360.0
		r = 30 + 15 * sin(pi * 2z/360)
		vso = length(n.Vertices)+1
		veo = vso
		t = 0
		while t < 2pi
			veo = vertex!(n, r * sin(t), r * cos(t), z)
			t += rstp
		end
		vsi= veo+1
		vei = vsi
		t = 0
		while t < 2pi
			vei = vertex!(n, (r-thk) * sin(t), (r-thk) * cos(t), z)
			t += rstp
		end

		if z > 0
			# outside
			for v in vso:veo-1
				face!(n, v-2steps, v, v-2steps+1)
				face!(n, v, v+1, v-2steps+1)
			end
			face!(n, veo, vso, vso- 2steps)
			face!(n, veo-2steps, veo, vso- 2steps)

			# inside
			for v in vsi:vei-1
				face!(n, v-2steps, v, v-2steps+1)
				face!(n, v, v+1, v-2steps+1)
			end
			face!(n, vei, vsi, vsi- 2steps)
			face!(n, vei-2steps, vei, vsi- 2steps)
		end
	end

	for k in 0:steps-2
		face!(n, vso+k, vsi+k, vsi+k+1)
		face!(n, vsi+k+1, vso+k+1, vso+k)
	end
	face!(n, vso+steps-1, vsi+steps-1, vso)
	face!(n, vso, vsi+steps-1, vsi)

	STL_ASCII(n, "circH.stl")

end

solid(n)


text:
const Evt{T} = NamedTuple{(:pos, :val),Tuple{Int,T}}
Evt(pos::Int, val::T) where T = (pos = pos, val = val)
1
"""
`xmax, xmin = movmaxmin(x, w)`

Moving maximum-minimum filter of window `w`
"""
movmaxmin(x::AbstractArray{T}, w) where T = movmaxmin!(similar(x), similar(x), x, w)
# the fastest - with circular buffers inlined, but less obvious
function _movmaxmin!(xmax::AbstractArray{T}, xmin::AbstractArray{T}, x::AbstractArray{T}, w::Int) where T
    U_buf 	   = Vector{Evt{T}}(undef, w)
    U_first    = Int(1)
    U_last     = Int(w)
    U_capacity = Int(0)

    L_buf 	   = Vector{Evt{T}}(undef, w)
    L_first    = Int(1)
    L_last     = Int(w)
    L_capacity = Int(0)

    Len = length(x)
    xi = x[1]
    @inbounds for i=2:Len
        xi_1 = xi
        xi = x[i]
        xmax[i-1] = ifelse(U_capacity > 0, U_buf[U_first].val, xi_1) ### x[get_front[U]]
        xmin[i-1] = ifelse(L_capacity > 0, L_buf[L_first].val, xi_1) ### x[get_front[L]]

        if xi > xi_1
            ### L = push_back[L, i-1];
            L_last = ifelse(L_last < w, L_last+1, 1)
            L_buf[L_last] = Evt(i-1, xi_1)
            L_capacity += 1
            ###
            if i == w + L_buf[L_first].pos
                ### L = pop_front[L];
                L_first = ifelse(L_first < w, L_first+1, 1)
                L_capacity -= 1
                ###
            end
            while U_capacity > 0
                if xi <= U_buf[U_last].val
                    if i == w + U_buf[U_first].pos
                        ### U = pop_front[U];
                        U_first = ifelse(U_first < w, U_first+1, 1)
                        U_capacity -= 1
                        ###
                    end
                    break
                end
                ### U = pop_back[U];
                U_last = ifelse(U_last > 1, U_last-1, w)
                U_capacity -= 1
                ###
            end
        else
            ### U = push_back[U, i-1];
            U_last = ifelse(U_last < w, U_last+1, 1)
            U_buf[U_last] = Evt(i-1, xi_1)
            U_capacity += 1
            ###
            if i == w + U_buf[U_first].pos
                ### U = pop_front[U];
                U_first = ifelse(U_first < w, U_first+1, 1)
                U_capacity -= 1
                ###
            end
            while L_capacity > 0
                if xi >= L_buf[L_last].val
                    if i == w + L_buf[L_first].pos
                        ### L = pop_front[L];
                        L_first = ifelse(L_first < w, L_first+1, 1)
                        L_capacity -= 1
                        ###
                    end
                    break
                end
                ### L = pop_back[L];
                L_last = ifelse(L_last > 1, L_last-1, w)
                L_capacity -= 1
                ###
            end
        end
    end
    xmax[Len] = ifelse(U_capacity > 0, U_buf[U_first].val, x[Len]) ### x[get_front[U]]
    xmin[Len] = ifelse(L_capacity > 0, L_buf[L_first].val, x[Len]) ### x[get_front[L]]

    xmin, xmax
end
# clearer code, but a bit slower
function movmaxmin!(xmax::AbstractArray{T}, xmin::AbstractArray{T}, x::AbstractArray{T}, w::Int) where T
    U_buf = RingBuffer{Evt{T}}(w)
    L_buf = RingBuffer{Evt{T}}(w)
    Len = length(x)
    xi = x[1]
    @inbounds for i=2:Len
        xi_1 = xi
        xi = x[i]
        xmax[i-1] = ifelse(length(U_buf) > 0, first(U_buf).val, xi_1) ### x[get_front[U]]
        xmin[i-1] = ifelse(length(L_buf) > 0, first(L_buf).val, xi_1) ### x[get_front[L]]

        if xi > xi_1
            push!(L_buf, Evt(i-1, xi_1))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, xi_1))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
    end

    xmax[Len] = ifelse(length(U_buf) > 0, first(U_buf).val, x[Len])
    xmin[Len] = ifelse(length(L_buf) > 0, first(L_buf).val, x[Len])

    xmax, xmin
end

"""
`xrange = movrange(x, w)`

Moving max-min range filter of window `w`.
"""
movrange(x::AbstractArray{T}, w) where T = movrange!(similar(x), x, w)
function movrange!(xrange::AbstractArray{T}, x::AbstractArray{T}, w::Int) where T
    U_buf = RingBuffer{Evt{T}}(w)
    L_buf = RingBuffer{Evt{T}}(w)
    Len = length(x)
    xi = x[1]
    @inbounds for i=2:Len
        xi_1 = xi
        xi = x[i]
        mx = ifelse(length(U_buf) > 0, first(U_buf).val, xi_1) ### x[get_front[U]]
        mn = ifelse(length(L_buf) > 0, first(L_buf).val, xi_1) ### x[get_front[L]]
        xrange[i-1] = mx - mn

        if xi > xi_1
            push!(L_buf, Evt(i-1, xi_1))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, xi_1))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
    end

    mx = ifelse(length(U_buf) > 0, first(U_buf).val, x[Len])
    mn = ifelse(length(L_buf) > 0, first(L_buf).val, x[Len])
    xrange[Len] = mx - mn

    xrange
end

"""
`xmax = movmax(x, w)`

Moving maximum filter of window `w`.
"""
movmax(x::AbstractArray{T}, w) where T = movmax!(similar(x), x, w)
function movmax!(xmax::AbstractArray{T}, x::AbstractArray{T}, w::Int) where T
    U_buf = RingBuffer{Evt{T}}(w)
    Len = length(x)
    xi = x[1]
    @inbounds for i=2:Len
        xi_1 = xi
        xi = x[i]
        xmax[i-1] = ifelse(length(U_buf) > 0, first(U_buf).val, xi_1) ### x[get_front[U]]
        ### ...
        if xi > xi_1
            ### ...
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, xi_1))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            ### ...
        end
    end
    xmax[Len] = ifelse(length(U_buf) > 0, first(U_buf).val, x[Len])

    xmax
end

"""
`xmin = movmin(x, w)`

Moving minimum filter of window `w`.
"""
movmin(x::AbstractArray{T}, w) where T = movmin!(similar(x), x, w)
function movmin!(xmin::AbstractArray{T}, x::AbstractArray{T}, w::Int) where T
    L_buf = RingBuffer{Evt{T}}(w)
    Len = length(x)
    xi = x[1]
    @inbounds for i=2:Len
        xi_1 = xi
        xi = x[i]
        # ...
        xmin[i-1] = ifelse(length(L_buf) > 0, first(L_buf).val, xi_1) ### x[get_front[L]]

        if xi > xi_1
            push!(L_buf, Evt(i-1, xi_1))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            # ...
        else
            # ...
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
    end
    xmin[Len] = ifelse(length(L_buf) > 0, first(L_buf).val, x[Len])

    xmin
end

"""
`xenvelope = movenvelope(x, w)`

Moving max-min envelope filter.
Envelope is defined as consecutive `movrange` of window `w1` (1-d dilation)
and `movmin` of window `w2` (1-d erosion). By default `w2 = w1`.
Filter delay = `(w1 + w2) ÷ 2`
"""
movenvelope(x::AbstractArray{T}, w1, w2 = w1 - 1) where T = movenvelope!(similar(x), x, w1, w2)
function movenvelope!(xenvelope::AbstractArray{T}, x::AbstractArray{T}, w1::Int, w2::Int = w1 - 1) where T
    movrange!(xenvelope, x, w1)
    if (w2 > 0)
        movmin!(xenvelope, xenvelope, w2)
    end
    xenvelope
end


text:
using DataDeps

register(DataDep(
  "thoracic-surgery",
  "http://archive.ics.uci.edu/ml/datasets/Thoracic+Surgery+Data",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00277/ThoraricSurgery.arff",
  "a92b446b681138619ac3e210d5641744fe77b039c2b7219456a21533e36f3015",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=17,
        feature_indices=1:16,
        categoric_indices=[1; 4:15],
        skipto=22,
        delim=",",
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

struct FwddiffFwd{F<:Base.Callable,i} <: Function
    f::F
end
FwddiffFwd(f::F, ::Val{i}) where {F<:Base.Callable,i} = FwddiffFwd{F,i}(f)

(fwd::FwddiffFwd{F,i})(xs...) where {F,i} = forwarddiff_fwd(fwd.f, xs, Val(i))

struct FwddiffBack{TX<:Any} <: Function end

(bck::FwddiffBack{TX})(ΔΩ, y_dual) where TX = forwarddiff_back(TX, ΔΩ, y_dual)

function forwarddiff_bc_fwd_back(f::Base.Callable, Xs::Tuple, ::Val{i}, ΔΩA::Any) where i
    # @info "RUN forwarddiff_bc_fwd_back(f, Xs, Val($i), ΔΩA)"
    fwd = FwddiffFwd(f, Val(i))
    TX = eltype(Xs[i])
    bck = FwddiffBack{TX}()
    bck.(ΔΩA, fwd.(Xs...))
end


text:
"""
```julia
LogisticLikelihood()
```

Bernoulli likelihood with a logistic link for the Bernoulli likelihood
```math
    p(y|f) = \\sigma(yf) = \\frac{1}{1+\\exp(-yf)},
```
(for more info see : [wiki page](https://en.wikipedia.org/wiki/Logistic_function))

---

For the analytic version the likelihood, it is augmented via:
```math
    p(y|f,ω) = exp(0.5(yf - (yf)^2 ω))
```
where ``ω ~ PG(ω | 1, 0)``, and `PG` is the Polya-Gamma distribution
See paper : [Efficient Gaussian Process Classification Using Polya-Gamma Data Augmentation](https://arxiv.org/abs/1802.06383)
"""
struct LogisticLikelihood{T<:Real} <: ClassificationLikelihood{T}
    c::AbstractVector{T}
    θ::AbstractVector{T}
    function LogisticLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function LogisticLikelihood{T}(c::AbstractVector{<:Real},θ::AbstractVector{<:Real}) where {T<:Real}
        new{T}(c,θ)
    end
end

function LogisticLikelihood()
    LogisticLikelihood{Float64}()
end

function init_likelihood(likelihood::LogisticLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        LogisticLikelihood{T}(abs.(rand(T,nSamplesUsed)),zeros(T,nSamplesUsed))
    else
        LogisticLikelihood{T}()
    end
end

function pdf(l::LogisticLikelihood,y::Real,f::Real)
    logistic(y*f)
end

function logpdf(l::LogisticLikelihood,y::T,f::T) where {T<:Real}
    -log(one(T)+exp(-y*f))
end

function Base.show(io::IO,model::LogisticLikelihood{T}) where T
    print(io,"Bernoulli Likelihood with Logistic Link")
end

function compute_proba(l::LogisticLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    pred = zeros(T,N)
    σ²_pred = zeros(T,N)
    for i in 1:N
        x = pred_nodes.*sqrt(max(σ²[i],zero(T))).+μ[i]
        pred[i] = dot(pred_weights,logistic.(x))
        σ²_pred[i] = max(dot(pred_weights,logistic.(x).^2)-pred[i]^2,zero(T))
    end
    return pred, σ²_pred
end

### Local Updates Section ###

function local_updates!(l::LogisticLikelihood{T},y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    l.c .= sqrt.(diag_cov+abs2.(μ))
    l.θ .= 0.5*tanh.(0.5.*l.c)./l.c
end

function sample_local!(l::LogisticLikelihood,y::AbstractVector,f::AbstractVector) where {T}
    pg = PolyaGammaDist()
    set_ω!(l,draw.([pg],[1.0],f))
end

### Natural Gradient Section ###

@inline ∇E_μ(l::LogisticLikelihood,::AOptimizer,y::AbstractVector) where {T} = (0.5*y,)
@inline ∇E_Σ(l::LogisticLikelihood,::AOptimizer,y::AbstractVector) where {T} = (0.5*l.θ,)

### ELBO Section ###

function expec_log_likelihood(l::LogisticLikelihood{T},i::AnalyticVI,y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    tot = -(0.5*length(y)*logtwo)
    tot += 0.5.*(dot(μ,y)-dot(l.θ,diag_cov)-dot(l.θ,μ))
    return tot
end

AugmentedKL(l::LogisticLikelihood,::AbstractVector) = PolyaGammaKL(l)

function PolyaGammaKL(l::LogisticLikelihood{T}) where {T}
    sum(broadcast(PolyaGammaKL,ones(T,length(l.c)),l.c,l.θ))
end

### Gradient Section ###

@inline grad_log_pdf(::LogisticLikelihood{T},y::Real,f::Real) where {T} = y*logistic(-y*f)

@inline hessian_log_pdf(::LogisticLikelihood{T},y::Real,f::Real) where {T<:Real} = -exp(y*f)/logistic(-y*f)^2


text:
using CVRPLIB, TSPLIB
using HGSCVRP_jll

const LIBHGSCVRP = HGSCVRP_jll.get_libhgscvrp_path()

include("c_api.jl")
include("cvrp.jl")
include("tsp.jl")

const C_DBL_MAX = floatmax(Cdouble)
const C_INT_MAX = typemax(Cint)

function reporting(visited_customers)
    routes = deepcopy(visited_customers)
    for r in routes
        r .-= 1
    end
    return routes
end


text:
"""
    ClimaComms

Abstracts the communications interface for the various CliMA components
in order to:
- enable the use of different backends as transports (MPI, SharedArrays,
  etc.), and
- transparently support single or double buffering for GPUs, depending
  on whether the transport has the ability to access GPU memory.

Use one of the `ClimaComms` backends, currently:
- `ClimaCommsMPI` uses MPI-2 asynchronous primitives through `MPI.jl`
  and supports a single buffer if the MPI implementation is CUDA-aware.
- `ClimaCommsSA` uses Julia's `Distributed` capability and
  `SharedArrays.jl` for within-node communication.
"""
module ClimaComms

"""
    AbstractCommsContext

The base type for a communications context. Each backend defines a
concrete subtype of this.
"""
abstract type AbstractCommsContext end

"""
    (pid, nprocs) = init(ctx::AbstractCommsContext)

Perform any necessary initialization for the specified backend. Return a
tuple of the processor ID and the number of participating processors.
"""
function init end

"""
    mypid(ctx::AbstractCommsContext)

Return the processor ID.
"""
function mypid end

"""
    iamroot(ctx::AbstractCommsContext)

Return `true` if the calling processor is the root processor.
"""
function iamroot end

"""
    nprocs(ctx::AbstractCommsContext)

Return the number of participating processors.
"""
function nprocs end

"""
    barrier(ctx::CC) where {CC <: AbstractCommsContext}

Perform a global synchronization across all participating processors.
"""
function barrier end
barrier(::Nothing) = nothing

"""
    reduce(ctx::CC, val, op) where {CC <: AbstractCommsContext}

Perform a reduction across all participating processors, using `op` as
the reduction operator and `val` as this rank's reduction value. Return
the result to the first processor only.
"""
function reduce end
reduce(::Nothing, val, op) = val

"""
    gather(ctx::AbstractCommsContext, array)

Gather an array of values from all processors into a single array,
concattenating along the last dimension.
"""
gather(::Nothing, array) = array

"""
    abort(ctx::CC, status::Int) where {CC <: AbstractCommsContext}

Terminate the caller and all participating processors with the specified
`status`.
"""
function abort end
abort(::Nothing, status::Int) = exit(status)

"""
    AbstractGraphContext

A context for communicating between processes in a graph.
"""
abstract type AbstractGraphContext end

"""
    graph_context(context::AbstractCommsContext, sendarray, sendpids, sendlengths, recvarray, recvpids, recvlengths)

Construct a communication context for exchanging neighbor data via a graph.

Arguments:
- `context`: the communication context on which to construct the graph context.
- `sendarray`: array containing data to send
- `sendpids`: list of processor IDs to send
- `sendlengths`: list of lengths of data to send to each process ID
- `recvarray`: array to receive data into
- `recvpids`: list of processor IDs to receive from
- `recvlengths`: list of lengths of data to receive from each process ID

This should return an `AbstractGraphContext` object.
"""
function graph_context end

"""
    start(ctx::AbstractGraphContext)

Initiate graph data exchange.
"""
function start end

"""
    progress(ctx::AbstractGraphContext)

Drive communication. Call after `start` to ensure that communication
proceeds asynchronously.
"""
function progress end

"""
    finish(ctx::AbstractGraphContext)

Complete the communications step begun by `start()`. After this returns,
data received from all neighbors will be available in the stage areas of
each neighbor's receive buffer.
"""
function finish end

include("singleton.jl")

end # module


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "mixtures sharing structure tests" begin
    spc = fully_factorized_circuit(SharedProbCircuit, 10, 3)
    pc = fully_factorized_circuit(ProbCircuit, 10)
    @test spc isa SharedProbCircuit
    @test spc isa ProbCircuit
    @test issum(spc)
    @test ismul(children(spc)[1])
    @test num_components(spc) == 3
    @test num_parameters_node(spc) == num_parameters_node(pc) * 3 == 3
    @test num_parameters(spc) == num_parameters(pc) * 3 == 63
    @test num_nodes(spc) == num_nodes(pc) == 32
    @test num_edges(spc) == num_edges(pc) == 31
end