text:
# Copyright (c) 2019 Simon Castano
# Distributed under the MIT software license, see the accompanying
# file COPYING or http://www.opensource.org/licenses/mit-license.php.

@testset "Transaction" begin
    @testset "TxIn" begin
        txin_hex = "7b1eabe0209b1fe794124575ef807057c77ada2138ae4fa8d6c4de0398a14f3f00000000494830450221008949f0cb400094ad2b5eb399d59d01c14d73d8fe6e96df1a7150deb388ab8935022079656090d7f6bac4c9a94e0aad311a4268e082a725f8aeae0573fb12ff866a5f01ffffffff"
        raw = hex2bytes(txin_hex)
        txin = TxIn(IOBuffer(raw))
        @test typeof(txin.prevout) == Outpoint
        @test bytes2hex(txin.prevout.txid) == "3f4fa19803dec4d6a84fae3821da7ac7577080ef75451294e71f9b20e0ab1e7b"
        @test txin.prevout.index == 0
        want = hex2bytes("30450221008949f0cb400094ad2b5eb399d59d01c14d73d8fe6e96df1a7150deb388ab8935022079656090d7f6bac4c9a94e0aad311a4268e082a725f8aeae0573fb12ff866a5f01")
        @test txin.scriptsig.data == [want]
        @test txin.sequence == 0xffffffff
    end
    @testset "TxOut" begin
        txout_hex = "f0ca052a010000001976a914cbc20a7664f2f69e5355aa427045bc15e7c6c77288ac"
        raw = hex2bytes(txout_hex)
        txout = TxOut(IOBuffer(raw))
        @test txout.value == 4999990000
        @test txout.scriptpubkey.data == [[0x76], [0xa9], hex2bytes("cbc20a7664f2f69e5355aa427045bc15e7c6c772"), [0x88], [0xac]]
    end

    raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
    raw_bip141 = hex2bytes("0100000000010115e180dc28a2327e687facc33f10f2a20da717e5548406f7ae8b4c811072f85601000000000fffffff0100b4f505000000001976a9141d7cd6c75c2e86f4cbf98eaed221b30bd9a0b92888ac02483045022100df7b7e5cda14ddf91290e02ea10786e03eb11ee36ec02dd862fe9a326bbcb7fd02203f5b4496b667e6e281cc654a2da9e4f08660c620a1051337fa8965f727eb19190121038262a6c6cec93c2d3ecd6c6072efea86d02ff8e3328bbd0242b20af3425990ac00000000")
    @testset "Parsing" begin
        @testset "Legacy" begin
            tx = Tx(IOBuffer(raw_tx))
            @test tx.version == 1
            @test tx.locktime == 0x00064319
            @testset "Inputs" begin
                @test length(tx.inputs) == 1
                want = hex2bytes("d1c789a9c60383bf715f3f6ad9d14b91fe55f3deb369fe5d9280cb1a01793f81")
                @test tx.inputs[1].prevout.txid == want
                @test tx.inputs[1].prevout.index == 0
                want = hex2bytes("6b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278a")
                @test serialize(tx.inputs[1].scriptsig) == want
                @test tx.inputs[1].sequence == 0xfffffffe
            end
            @testset "Outputs" begin
                @test length(tx.outputs) == 2
                want = 32454049
                @test tx.outputs[1].value == want
                want = hex2bytes("1976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac")
                @test serialize(tx.outputs[1].scriptpubkey) == want
                want = 10011545
                @test tx.outputs[2].value == want
                want = hex2bytes("1976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac")
                @test serialize(tx.outputs[2].scriptpubkey) == want
            end
        end
        @testset "Segwit" begin
            tx = Tx(IOBuffer(raw_bip141))
            @test tx.version == 0x00000001
            @test tx.flag == 0x01
            hex1 = "3045022100df7b7e5cda14ddf91290e02ea10786e03eb11ee36ec02dd862fe9a326bbcb7fd02203f5b4496b667e6e281cc654a2da9e4f08660c620a1051337fa8965f727eb191901"
            hex2 = "038262a6c6cec93c2d3ecd6c6072efea86d02ff8e3328bbd0242b20af3425990ac"
            @test tx.witnesses[1].data[1] == hex2bytes(hex1)
            @test tx.witnesses[1].data[2] == hex2bytes(hex2)
        end
    end
    @testset "Serialize" begin
        tx = Tx(IOBuffer(raw_tx))
        @test serialize(tx) == raw_tx

        tx = Tx(IOBuffer(raw_bip141))
        @test serialize(tx) == raw_bip141
    end
    @testset "Is CoinbaseTx" begin
        raw_tx = hex2bytes("01000000010000000000000000000000000000000000000000000000000000000000000000ffffffff5e03d71b07254d696e656420627920416e74506f6f6c20626a31312f4542312f4144362f43205914293101fabe6d6d678e2c8c34afc36896e7d9402824ed38e856676ee94bfdb0c6c4bcd8b2e5666a0400000000000000c7270000a5e00e00ffffffff01faf20b58000000001976a914338c84849423992471bffb1a54a8d9b1d69dc28a88ac00000000")
        stream = IOBuffer(raw_tx)
        tx = Tx(stream)
        @test iscoinbase(tx)
    end
    @testset "coinbase_height" begin
        raw_tx = hex2bytes("01000000010000000000000000000000000000000000000000000000000000000000000000ffffffff5e03d71b07254d696e656420627920416e74506f6f6c20626a31312f4542312f4144362f43205914293101fabe6d6d678e2c8c34afc36896e7d9402824ed38e856676ee94bfdb0c6c4bcd8b2e5666a0400000000000000c7270000a5e00e00ffffffff01faf20b58000000001976a914338c84849423992471bffb1a54a8d9b1d69dc28a88ac00000000")
        stream = IOBuffer(raw_tx)
        tx = Tx(stream)
        @test coinbase_height(tx) == 465879
        raw_tx = hex2bytes("0100000001813f79011acb80925dfe69b3def355fe914bd1d96a3f5f71bf8303c6a989c7d1000000006b483045022100ed81ff192e75a3fd2304004dcadb746fa5e24c5031ccfcf21320b0277457c98f02207a986d955c6e0cb35d446a89d3f56100f4d7f67801c31967743a9c8e10615bed01210349fc4e631e3624a545de3f89f5d8684c7b8138bd94bdd531d2e213bf016b278afeffffff02a135ef01000000001976a914bc3b654dca7e56b04dca18f2566cdaf02e8d9ada88ac99c39800000000001976a9141c4bc762dd5423e332166702cb75f40df79fea1288ac19430600")
        stream = IOBuffer(raw_tx)
        tx = Tx(stream)
        @test_throws AssertionError coinbase_height(tx)
    end
end


text:
using Distributions
using LinearAlgebra
using BenchmarkTools

# ------------------------------------------------------------------------------

function find_Q(A, l)
    #=
    Given an m × n matrix A, and an integer l, compute an m × l orthonormal
    matrix Q whose range approximates the range of A.
    =#
    m, n = size(A)
    Ω = rand(Normal(), n, l)
    Y = A * Ω
    Q, R = qr(Y)
    return Q
end

# ------------------------------------------------------------------------------

function randomized_SVD(A, k)
    #=
    Given an m × n matrix A, a target number k of singular vectors, and an
    exponent q (say q = 1 or q = 2), this procedure computes an approximate
    rank-2k factorization UΣVt, where U and V are orthonormal and Σ is
    nonnegative and diagonal.
    =#
    Q = find_Q(A, 2*k)
    B = Q' * A
    S, Σ, Vt = svd(B)
    U = Q * S
    return U, Σ, Vt
end

# ------------------------------------------------------------------------------

m = 2000
n = 20000
k = 10
# Construct low-rank matrix
A = rand(m, k) * rand(k, n)
println("Rank of A: ", rank(A))
println("Size of A: ", size(A))

println("Throwaway test:")
@time svd(A)
@time randomized_SVD(A, k)
println("Actual test:")
@benchmark svd(A)
@benchmark randomized_SVD(A, k)

println("Completed")

# https://stackoverflow.com/questions/53162206/why-does-my-randomized-svd-implementation-use-so-much-memory

ext: jl

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using DeepForest, Flux, ProgressMeter, Plots, StatsPlots
using StatsBase: mode

x = rand(Float32, 2, 1000);
x[1, :] .*= 2*π;
x[1, :] .-= π;
x[2, :] .*= 3;
x[2, :] .-= 1.5;
y = Int64.(x[2, :] .< sin.(x[1, :])) .+ 1;

features = Dict(
    7 => [1, 2],
    6 => [1, 2],
    5 => [1, 2],
    4 => [1, 2],
    3 => [1, 2],
    2 => [1, 2],
    1 => [1, 2]
)

model = Node(features, 10, 2, 1)

tree_train!(2500, model, x, y)

ŷ = predict(model, x)
scatter(x[1, :], x[2, :], color=ŷ, legend=false)
plot!((x) -> sin(x), xlims=(-π, π), color=:green, linewidth=5)

imp = importance(model, x)
ρ = sum(values(imp))
for (k, v) ∈ imp
    imp[k] = v / ρ
end
plot_importance(imp, x)


text:
# Implementations of LR and Mom policies (some depend on coffee-break.jl)

export LRPolicy, get_learning_rate, MomPolicy, get_momentum

############################################################
# Learning rate policy
############################################################

module LRPolicy
using ..Mocha
using Compat
type Fixed <: LearningRatePolicy
  base_lr :: AbstractFloat
end

# base_lr * gamma ^ (floor(iter / stepsize))
type Step <: LearningRatePolicy
  base_lr  :: AbstractFloat
  gamma    :: AbstractFloat
  stepsize :: Int
end

# base_lr * gamma ^ iter
type Exp <: LearningRatePolicy
  base_lr :: AbstractFloat
  gamma   :: AbstractFloat
end

type Inv <: LearningRatePolicy
  base_lr :: AbstractFloat
  gamma   :: AbstractFloat
  power   :: AbstractFloat
end

# curr_lr *= gamma whenever performance
# drops on the validation set
function decay_on_validation_listener(policy, key::AbstractString, coffee_lounge::CoffeeLounge, net::Net, state::SolverState)
  stats = get_statistics(coffee_lounge, key)
  index = sort(collect(keys(stats)))
  if length(index) > 1
    if (policy.higher_better && stats[index[end]] < stats[index[end-1]]) ||
      (!policy.higher_better && stats[index[end]] > stats[index[end-1]])
      # performance drop
      Mocha.info(@sprintf("lr decay %e -> %e", policy.curr_lr, policy.curr_lr*policy.gamma))
      policy.curr_lr *= policy.gamma

      # revert to a previously saved "good" snapshot
      if isa(policy.solver, Solver)
        Mocha.info("reverting to previous saved snapshot")
        solver_state = load_snapshot(net, policy.solver.params[:load_from], state)
        Mocha.info("snapshot at iteration $(solver_state.iter) loaded")
        copy_solver_state!(state, solver_state)
      end
    end
  end
end

type DecayOnValidation <: LearningRatePolicy
  gamma       :: AbstractFloat

  key           :: AbstractString
  curr_lr       :: AbstractFloat
  min_lr        :: AbstractFloat
  listener      :: Function
  solver        :: Any
  initialized   :: Bool
  higher_better :: Bool # set to false if performance score is the lower the better

  DecayOnValidation(base_lr, key, gamma=0.5, min_lr=1e-8; higher_better=true) = begin
    policy = new(gamma, key, base_lr, min_lr)
    policy.solver = nothing
    policy.listener = (coffee_lounge,net,state) -> begin
      if policy.curr_lr < policy.min_lr
        # do nothing if we already fall below the minimal learning rate
        return
      end
      decay_on_validation_listener(policy, key, coffee_lounge, net, state)
    end
    policy.initialized = false
    policy.higher_better = higher_better

    policy
  end
end

using Compat
type Staged <: LearningRatePolicy
  stages     :: Vector{@compat(Tuple{Int, LearningRatePolicy})}
  curr_stage :: Int

  Staged(stages...) = begin
    accum_stages = Array{@compat(Tuple{Int, LearningRatePolicy})}(length(stages))
    accum_iter = 0
    for i = 1:length(stages)
      (n, lrp) = stages[i]
      accum_iter += n
      accum_stages[i] = (accum_iter, convert(LearningRatePolicy, lrp))
    end

    new(accum_stages, 1)
  end
end

end # module LRPolicy

get_learning_rate(policy::LearningRatePolicy) = policy.base_lr # Need an initial rate to create the state
get_learning_rate(policy::LRPolicy.Staged) = get_learning_rate(policy.stages[policy.curr_stage][2])
get_learning_rate(policy::LRPolicy.DecayOnValidation) = policy.curr_lr

get_learning_rate(policy::LRPolicy.Fixed, state::SolverState) = policy.base_lr
get_learning_rate(policy::LRPolicy.Step, state::SolverState) =
    policy.base_lr * policy.gamma ^ (floor(state.iter / policy.stepsize))
get_learning_rate(policy::LRPolicy.Exp, state::SolverState) =
    policy.base_lr * policy.gamma ^ state.iter
get_learning_rate(policy::LRPolicy.Inv, state::SolverState) =
    policy.base_lr * (1 + policy.gamma * state.iter) ^ (-policy.power)

function setup(policy::LRPolicy.DecayOnValidation, validation::ValidationPerformance, solver::Solver)
  register(validation, policy.listener)
  policy.solver = solver
end

get_learning_rate(policy::LRPolicy.DecayOnValidation, state::SolverState) = begin
    # TODO this is really for a SolverState{SGDSolverState} but that's a forward declaration
  if !policy.initialized
    if state.internal.learning_rate > 0
      # state.learning_rate is initialized to 0, if it is non-zero, then this might
      # be loaded from some saved snapshot, we try to align with that
      @info("Switching to base learning rate $(state.specific.learning_rate)")
      policy.curr_lr = state.internal.learning_rate
    end
    policy.initialized = true
  end

  policy.curr_lr
end

function get_learning_rate(policy::LRPolicy.Staged, state::SolverState)
  if policy.curr_stage == length(policy.stages)
    # already in the last stage, stick there forever
  else
    maxiter = policy.stages[policy.curr_stage][1]
    while state.iter >= maxiter && policy.curr_stage < length(policy.stages)
      policy.curr_stage += 1
      @info("Staged learning rate policy: switching to stage $(policy.curr_stage)")
      maxiter = policy.stages[policy.curr_stage][1]
    end
  end
  return get_learning_rate(policy.stages[policy.curr_stage][2], state)
end

############################################################
# Momentum policy
############################################################

module MomPolicy
using ..Mocha.MomentumPolicy
using Compat
type Fixed <: MomentumPolicy
  base_mom :: AbstractFloat
end

# min(base_mom * gamma ^ (floor(iter / stepsize)), max_mom)
type Step <: MomentumPolicy
  base_mom :: AbstractFloat
  gamma    :: AbstractFloat
  stepsize :: Int
  max_mom  :: AbstractFloat
end

type Linear <: MomentumPolicy
  base_mom :: AbstractFloat
  gamma    :: AbstractFloat
  stepsize :: Int
  max_mom  :: AbstractFloat
end

using Compat
type Staged <: MomentumPolicy
  stages     :: Vector{@compat(Tuple{Int, MomentumPolicy})}
  curr_stage :: Int

  Staged(stages...) = begin
    accum_stages = Array{@compat(Tuple{Int, MomentumPolicy})}(length(stages))
    accum_iter = 0
    for i = 1:length(stages)
      (n, mmp) = stages[i]
      accum_iter += n
      accum_stages[i] = (accum_iter, convert(MomentumPolicy, mmp))
    end

    new(accum_stages, 1)
  end
end

end # module MomPolicy

get_momentum(policy::MomentumPolicy) = policy.base_mom # Need an initial rate to create the state
get_momentum(policy::MomPolicy.Fixed, state::SolverState) = policy.base_mom
get_momentum(policy::MomPolicy.Step, state::SolverState) =
    min(policy.base_mom * policy.gamma ^ (floor(state.iter / policy.stepsize)), policy.max_mom)
get_momentum(policy::MomPolicy.Linear, state::SolverState) =
    min(policy.base_mom + floor(state.iter / policy.stepsize) * policy.gamma, policy.max_mom)

function get_momentum(policy::MomPolicy.Staged, state::SolverState)
  if policy.curr_stage == length(policy.stages)
    # already in the last stage, stick there forever
  else
    maxiter = policy.stages[policy.curr_stage][1]
    while state.iter >= maxiter && policy.curr_stage < length(policy.stages)
      policy.curr_stage += 1
      @info("Staged momentum policy: switching to stage $(policy.curr_stage)")
      maxiter = policy.stages[policy.curr_stage][1]
    end
  end
  return get_momentum(policy.stages[policy.curr_stage][2], state)
end

####


text:
#This example fits a GP regression model to the Mauna Loa CO2 data set. This data set is regularly updated and can be found at co2now.org/images/stories/data/co2-mlo-monthly-noaa-esrl.xls
#
#This example follows from Chapter 5 of Gaussian Processes for Machine Learning, Rasmussen and Williams (2006)
############################################################################################

using Gadfly, GaussianProcesses

data = readcsv("CO2_data.csv")

year = data[:,1]; co2 = data[:,2];
x = year[year.<2004]; y = co2[year.<2004];
xpred = year[year.>=2004]; ypred = co2[year.>=2004];

mConst = MeanConst(mean(y))       #Fit the constant mean function

#Kernel is represented as a sum of kernels
kernel = SE(4.0,4.0) + Periodic(0.0,1.0,0.0)*SE(4.0,0.0) + RQ(0.0,0.0,-1.0) + SE(-2.0,-2.0)

gp = GP(x,y,mConst,kernel,-2.0)   #Fit the GP

plot(gp,clim=(2004.0,2024.0))  #Gadfly can take a while to load


text:
using Manopt, ManifoldsBase, Manifolds, LinearAlgebra, BenchmarkTools, Test

M = Euclidean(2)
x = [4.0, 2.0]
x0 = [5.0, 2.0]
f(M, y) = distance(M, y, x)
function ∂f(M, y)
    if distance(M, x, y) == 0
        return zero_vector(M, y)
    end
    return -2 * log(M, y, x) / max(10 * eps(Float64), distance(M, x, y))
end
x1 = subgradient_method(M, f, ∂f, x0)
@btime subgradient_method($M, $f, $∂f, $x0)

function ∂f!(M, X, y)
    d = distance(M, x, y)
    if d == 0
        return zero_vector!(M, X, y)
    end
    log!(M, X, y, x)
    X .*= -2 / max(10 * eps(Float64), d)
    return X
end
x2 = copy(x0)
subgradient_method!(M, f, ∂f!, x2; evaluation=MutatingEvaluation())
@btime subgradient_method!($M, $f, $∂f!, x3; evaluation=$(MutatingEvaluation())) setup = (
    x3 = deepcopy($x0)
)

@test x1 == x2


text:
using MetaGen
using Gen
using PyPlot
const plt = PyPlot
using Random

Random.seed!(1234);

# Generates a v matrix
params = Video_Params()
v = make_visual_system(params)
print(v)

# There's just one receptive field
receptive_fields = make_receptive_fields()
num_receptive_fields = length(receptive_fields)

# Generates scenes and detections.
num_videos = 2
num_frames = 2

fs = fill(num_frames, num_videos) #number of frames per video
ps = fill(params, num_videos)
vs = fill(v, num_videos)
receptive_fieldses = fill(receptive_fields, num_videos)

Gen.load_generated_functions()
scenes_trace = Gen.simulate(video_map, (fs, ps, vs, receptive_fieldses))

#Gen.get_choices(scenes_trace)

#converts trace to arrays required for `unfold_particle_filter`
objects_observed = Matrix{Array{Array{Detection2D}}}(undef, num_videos, num_frames)
camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)

for v=1:num_videos
    for f=1:num_frames
        temp = scenes_trace[v => :frame_chain => f => 1 => :observations_2D]
	    #temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
	    #objects_observed[v, f] = temp_sorted_into_rfs
        objects_observed[v, f] = temp
        camera_trajectories[v, f] = scenes_trace[v => :frame_chain => f => :camera]
    end
end

# Does inference

num_particles = 1
inference_traces = unfold_particle_filter(num_particles, objects_observed, camera_trajectories, num_receptive_fields)

# Gets best particle
top_index = 0
top_score = -Inf
for (i, trace) in enumerate(inference_traces)
    score = get_score(trace)
    if score > top_score
        global top_score = score
        global top_index = i
    end
end

best_trace = inference_traces[top_index]

# evaluates inferred visual_system.
v_hat = Matrix{Real}(undef, length(params.possible_objects), 2)

#Gen.get_choices(best_trace)
for j=1:length(params.possible_objects)
    v_hat[j, 1] = best_trace[:v_matrix => (:lambda_fa, j)]
    v_hat[j, 2] = best_trace[:v_matrix => (:miss_rate, j)]
end

L1_errors = vec(abs.(v - v_hat))

plt.clf()
plt.title("V matrix L1 errors")
bins = [0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
plt.hist(L1_errors, bins=bins)
plt.xlim(0, 1)
plt.savefig("errors.pdf")

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using XGBoost, Printf

export predictdemand, trainmodel, crossvalidate, defaultvariables, getcvdata, plotcv

const defaultvariables = [:localhour, :weekend01, :temp_monthly, :ranked_month, :temp_top3,
                            :temp1_mean, :temp1_qlow, :temp1_qhigh, :demandpercapita, :gdppercapita]

function predictdemand(; variables=defaultvariables, gisregion="Europe8",
            sspscenario="ssp2-34", sspyear=2050, era_year=2018, numcenters=3, mindist=3.3,
            nrounds=100, max_depth=7, eta=0.05, subsample=0.75, metrics=["mae"], more_xgoptions...)
    df, offsets, pop = buildtrainingdata(; gisregion=gisregion, sspscenario=sspscenario,
            sspyear=sspyear, era_year=era_year, numcenters=numcenters, mindist=mindist)
    regionlist = unique(df[:, :country])
    numhours = 24*daysinyear(era_year)
    demandpercapita = df[1:numhours:end, :demandpercapita]      # MWh/year/capita
    select!(df, variables)
    traindata = Matrix(df)
    model = trainmodel(; nrounds=nrounds, max_depth=max_depth, eta=eta, subsample=subsample, metrics=metrics, more_xgoptions...)
    normdemand = XGBoost.predict(model, traindata)              # mean(normdemand) == 1
    numreg = length(regionlist)
    demand = reshape(normdemand, (numhours, numreg)) .* (demandpercapita/8760)' .* pop'    # MW
    println("\nConverting synthetic demand to UTC...")
    for r = 1:numreg
        demand[:,r] = circshift(demand[:,r], round(Int, -offsets[r]))
    end
    println("\nSaving...")
    JLD.save(in_datafolder("output",
            "SyntheticDemand_$(gisregion)_$sspscenario-$(sspyear)_$era_year.jld"),
            "demand", demand, compress=true)
    nothing
end

function trainmodel(; variables=defaultvariables, nrounds=100, xgoptions...)
    df_train, offsets = loadtrainingdata()
    println("\nTraining model...")
    select!(df_train, variables)
    traindata = Matrix(df_train)
    normdemand = loaddemanddata()[:, :normdemand]

    model = xgboost(traindata, nrounds; label=normdemand, xgoptions...)
end

function crossvalidate(; variables=defaultvariables, nrounds=100, max_depth=7, eta=0.05, subsample=0.75, metrics=["mae"], more_xgoptions...)
    df_train, offsets = loadtrainingdata()
    regionlist = unique(df_train[:, :country])
    select!(df_train, variables)
    traindata = Matrix(df_train)
    normdemand = loaddemanddata()[:, :normdemand]

    numreg = length(regionlist)
    params = Any["max_depth"=>round(Int, max_depth), "eta"=>eta, "subsample"=>subsample, "metrics"=>metrics, more_xgoptions...]

    models = nfold_cv_return(traindata, nrounds, numreg; label=normdemand, metrics=metrics, param=params)   # "rmse" or "mae"
    display(importance(models[1], string.(variables)))

    return models
end

function getcvdata(models; variables=defaultvariables, era_year=2018)
    df_train, offsets = loadtrainingdata()
    regionlist = unique(df_train[:, :country])
    numhours = 24*daysinyear(era_year)
    numreg = length(models)
    demandpercapita = df_train[1:numhours:end, :demandpercapita]      # MWh/year/capita
    select!(df_train, variables)
    traindata = Matrix(df_train)
    normdemand = loaddemanddata()[:, :normdemand]
    normdemand_predicted = similar(normdemand)
    for reg = 1:numreg
        rows = numhours*(reg-1) + 1 : numhours*reg
        normdemand_predicted[rows] = XGBoost.predict(models[reg], traindata[rows,:])          # mean(normdemand) == 1
    end
    return reshape(normdemand, (numhours,numreg)), reshape(normdemand_predicted, (numhours,numreg)), regionlist
end

# Same as XGBoost.nfold_cv() but deterministic: split the training data into nfold equal parts using sequential indices.
# Also returns the models used in the last iteration.
function nfold_cv_return(data, num_boost_round::Integer = 10, nfold::Integer = 3; label = Union{},
                  param=[], metrics=[], obj = Union{}, feval = Union{}, fpreproc = Union{},
                  show_stdv = true, seed::Integer = 0, kwargs...)
    dtrain = XGBoost.makeDMatrix(data, label)
    results = String[]
    cvfolds = mknfold_deterministic(dtrain, nfold, param, metrics, fpreproc=fpreproc, kwargs = kwargs)
    for i in 1:num_boost_round
        for f in cvfolds
            XGBoost.update(f.bst, 1, f.dtrain, obj = obj)
        end
        res = XGBoost.aggcv([XGBoost.eval_set(f.bst, f.watchlist, i, feval = feval) for f in cvfolds],
                    show_stdv = show_stdv)
        push!(results, res)
        @printf(stderr, "%s\n", res)
    end
    models = [f.bst for f in cvfolds]
    return models
end

# Same as XGBoost.mknfold() but deterministic: split the training data into nfold equal parts using sequential indices
function mknfold_deterministic(dall::DMatrix, nfold::Integer, param, evals=[]; fpreproc = Union{},
                 kwargs = [])
    idx = collect(1:XGBoost.XGDMatrixNumRow(dall.handle))
    kstep = size(idx)[1] / nfold
    idset = [idx[round(Int64, (i-1) * kstep) + 1 : min(size(idx)[1],round(Int64, i * kstep))] for i in 1:nfold]
    ret = XGBoost.CVPack[]
    for k in 1:nfold
        selected = Int[]
        for i in 1:nfold
            if k != i
                selected = vcat(selected, idset[i])
            end
        end
        dtrain = XGBoost.slice(dall, selected)
        dtest = XGBoost.slice(dall, idset[k])
        if typeof(fpreproc) == Function
            dtrain, dtest, tparam = XGBoost.fpreproc(dtrain, dtest, deepcopy(param))
        else
            tparam = param
        end
        plst = vcat([itm for itm in param], [("eval_metric", itm) for itm in evals])
        plst = vcat(plst, [(string(itm[1]), string(itm[2])) for itm in kwargs])
        push!(ret, XGBoost.CVPack(dtrain, dtest, plst))
    end
    return ret
end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
@testset "713.subarray-product-less-than-k.jl" begin
    @test num_subarray_productLess_than_k([10, 5, 2, 6], 100) == 8
    @test num_subarray_productLess_than_k([40, 33, 8, 49, 36, 7, 49, 46, 32, 7, 20, 28, 43,
                                           2, 41, 35, 2, 2, 19, 28, 16, 9, 30, 40, 23, 17,
                                           47, 33, 18, 41, 28, 28, 7, 16, 36, 40, 48, 50,
                                           32, 16, 5, 2, 47, 40, 39, 12, 50, 49, 13, 4],
                                          300) == 70
end


text:
# # Ocean convection example
#
# In this example, two-dimensional convection into a stratified fluid
# mixes a phytoplankton-like tracer. This example demonstrates how
#
#   * to set boundary conditions;
#   * to defined and insert a user-defined forcing function into a simulation.
#   * to use the `TimeStepWizard` to manage and adapt the simulation time-step.
#
# To begin, we load Oceananigans, a plotting package, and a few miscellaneous useful packages.

using Oceananigans, PyPlot, Random, Printf

# ## Parameters
#
# We choose a modest two-dimensional resolution of 128² in a 64² m² domain ,
# implying a resolution of 0.5 m. Our fluid is initially stratified with
# a squared buoyancy frequency
#
# $ N² = 10⁻⁵ \rm{s⁻²} $
#
# and a surface buoyancy flux
#
# $ Q_b = 10⁻⁸ \rm{m³ s⁻²} $
#
# Because we use the physics-based convection whereby buoyancy flux by a
# positive vertical velocity implies positive flux, a positive buoyancy flux
# at the top of the domain carries buoyancy out of the fluid and causes convection.
# Finally, we end the simulation after 1 day.

Nz = 128
Lz = 64.0
N² = 1e-5
Qb = 1e-8
end_time = day / 2

# ## Creating boundary conditions
#
# Create boundary conditions. Note that temperature is buoyancy in our problem.
#

buoyancy_bcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Flux, Qb),
                                       bottom = BoundaryCondition(Gradient, N²))

# ## Define a forcing function
#
# Our forcing function roughly corresponds to the growth of phytoplankton in light
# (with a penetration depth of 16 meters here), and death due to natural mortality
# at a rate of 1 phytoplankton unit per second.

growth_and_decay = SimpleForcing((x, y, z, t) -> exp(z/16) - 1)

## Instantiate the model
model = Model(
                   grid = RegularCartesianGrid(size = (Nz, 1, Nz), length = (Lz, Lz, Lz)),
                closure = ConstantIsotropicDiffusivity(ν=1e-4, κ=1e-4),
               coriolis = FPlane(f=1e-4),
                tracers = (:b, :plankton),
               buoyancy = BuoyancyTracer(),
                forcing = ModelForcing(plankton=growth_and_decay),
    boundary_conditions = BoundaryConditions(b=buoyancy_bcs)
)

# Set makeplot = true to live-update a plot of vertical velocity, temperature, and salinity
# as the simulation runs.

makeplot = false

## Set initial condition. Initial velocity and salinity fluctuations needed for AMD.
Ξ(z) = randn() * z / Lz * (1 + z / Lz) # noise
b₀(x, y, z) = N² * z + N² * Lz * 1e-6 * Ξ(z)
set!(model, b=b₀)

## A wizard for managing the simulation time-step.
wizard = TimeStepWizard(cfl=0.1, Δt=1.0, max_change=1.1, max_Δt=90.0)

## Create a plot
fig, axs = subplots(ncols=2, figsize=(10, 6))

"""
    makeplot!(axs, model)

Make side-by-side x-z slices of vertical velocity and plankton associated with `model` in `axs`.
"""
function makeplot!(axs, model)
    xC = repeat(model.grid.xC, 1, model.grid.Nz)
    zF = repeat(reshape(model.grid.zF[1:end-1], 1, model.grid.Nz), model.grid.Nx, 1)
    zC = repeat(reshape(model.grid.zC, 1, model.grid.Nz), model.grid.Nx, 1)

    sca(axs[1]); cla()
    ## Calling the Array() constructor here allows plots to be made for GPU models:
    pcolormesh(xC, zF, Array(interior(model.velocities.w))[:, 1, :])
    title("Vertical velocity")
    xlabel("\$ x \$ (m)")
    ylabel("\$ z \$ (m)")

    sca(axs[2]); cla()
    ## Calling the Array() constructor here allows plots to be made for GPU models:
    pcolormesh(xC, zC, Array(interior(model.tracers.plankton))[:, 1, :])
    title("Phytoplankton concentration")
    xlabel("\$ x \$ (m)")
    axs[2].tick_params(left=false, labelleft=false)

    suptitle(@sprintf("\$ t = %.2f\$ hours", model.clock.time / hour))
    [ax.set_aspect(1) for ax in axs]
    pause(0.01)
    gcf()
    return nothing
end

# Run the model:

while model.clock.time < end_time
    update_Δt!(wizard, model)
    walltime = @elapsed time_step!(model, 100, wizard.Δt)

    ## Print a progress message
    @printf("progress: %.1f %%, i: %04d, t: %s, Δt: %s, wall time: %s\n",
            model.clock.time / end_time * 100, model.clock.iteration,
            prettytime(model.clock.time), prettytime(wizard.Δt), prettytime(walltime))

    makeplot && makeplot!(axs, model)
end

# Plot the result at the end
makeplot!(axs, model)
gcf()


text:
using PyPlot;
using MAT;
using jInv.Mesh
using jInvVis
file = matopen("hand.mat")
B = read(file, "B");
close(file);
plotModel(convert(Array{Float64,3},B));

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
module ShaleDrillingLikelihood_Royalty_Test

# using Revise

# using Juno
using ShaleDrillingLikelihood
using Test
using StatsFuns
using Calculus
using Optim
using Random
using InteractiveUtils
using BenchmarkTools
using Distributions: _tnmom1

using ShaleDrillingLikelihood: RoyaltyModelNoHet,
    idx_royalty_ρ, idx_royalty_ψ, idx_royalty_β, idx_royalty_κ,
    theta_royalty_ρ, theta_royalty_ψ, theta_royalty_β, theta_royalty_κ,
    η12,
    lik_royalty, lik_loglik_royalty,
    theta_royalty_check,
    simloglik_royalty!,
    grad_simloglik_royalty!,
    # _LLm,
    logsumexp_and_softmax!,
    ObservationRoyalty, DataRoyalty, _y, _x, _xbeta, num_choices, _num_x,
    SimulationDraws,
    update_ψ1!, update_dψ1dθρ!, _psi1, _u, _v, _dψ1dθρ,
    update_xbeta!, _am, _bm, _cm, _qm,
    llthreads!, _i, _nparm,
    logsumexp!,
    theta_royalty_level_to_cumsum,
    theta_royalty_cumsum_to_level,
    kappa_level_to_cumsum,
    kappa_cumsum_to_level

@testset "RoyaltyModelNoHet" begin
    k = 3
    nobs = 1_000
    L = 3
    RM = RoyaltyModelNoHet()
    Random.seed!(1234)

    X       = randn(k,nobs)
    epsilon = randn(nobs)

    # I forgot a negative sign originally
    @testset "Check tnmom" begin
        a,b = (-2,1)
        dlogdiffcdf = (normpdf(b) - normpdf(a)) / (normcdf(b) - normcdf(a))
        @test dlogdiffcdf ≈ -_tnmom1(a,b)
    end

    theta0 = [-2.0, 2.0, 2.0, -0.5, 0.5]  # β, κ
    theta = vcat(theta0[1:end-1], sqrt(2*(theta0[end] - theta0[end-1])))
    @test kappa_level_to_cumsum(theta0[end-1:end]) ≈ theta[end-1:end]
    @test kappa_cumsum_to_level(theta[end-1:end]) ≈ theta0[end-1:end]
    @test length(theta) == k + L - 1

    rstar = X'*theta[1:k] .+ epsilon
    l = map((r) ->  searchsortedfirst(theta0[k+1:end], r), rstar)
    data = DataRoyalty(RM, l, X, 1:L)
    update_xbeta!(data, theta[1:k])

    @test theta ≈ theta_royalty_level_to_cumsum(data, theta0)
    @test theta0 ≈ theta_royalty_cumsum_to_level(data, theta)

    @test length(theta) == _nparm(data) == _nparm(first(first(data)))

    @test theta_royalty_ρ(    data, theta) == Float64[]
    @test theta_royalty_ψ(    data, theta) == Float64[]
    @test all(theta_royalty_β(data, theta) .== theta[1:3])
    @test theta_royalty_κ(    data, theta, 1) == theta[4]
    @test theta_royalty_κ(    data, theta, 2) == theta[5]

    @test idx_royalty_κ(data, 1) == 4
    @test idx_royalty_κ(data, 2) == 5
    @test_throws DomainError idx_royalty_κ(data, L+1)
    @test_throws DomainError idx_royalty_κ(data, 0)

    M = 10
    RT = SimulationDraws(M,nobs,0)
    am = _am(RT)
    bm = _bm(RT)
    cm = _cm(RT)

    # test η12s
    η12s = [η12(first(grp), theta, r) for (grp,r) in zip(data, rstar)]
    @test all(x[1] < x[2] for x in η12s )

    # form functions for stuff
    grad = similar(theta)
    # @code_warntype llthreads!(grad, theta, RM, data, true)

    f(parm) = - llthreads!(grad, parm, data, false)

    function fg!(g, parm)
        LL = llthreads!(g, parm, data, true)
        g .*= -1
        return -LL
    end

    # check gradient
    println("Print me to prevent error on v1.2")
    fg!(grad, theta)
    @test grad ≈ Calculus.gradient(f, theta)

    # check that it solves
    res = optimize(OnceDifferentiable(f, fg!, fg!, theta), theta*0.1, BFGS(), Optim.Options(time_limit = 1.0))
    @show res
    @show res.minimizer
    @test maximum(abs.(res.minimizer .- theta)) < 0.25
end

@testset "RoyaltyModel" begin

    nobs = 500  # observations
    k = 3       # number of x
    L = 3       # ordered choices
    M = 5_000     # simulations

    RM = RoyaltyModel()

    # make data
    X      = randn(k,nobs)
    eps    = randn(nobs)
    theta0  = [0.1, 1.0,    -2.0, 2.0, 2.0,    -0.6, 0.6]  # dψdρ, ψ, β, κ
    theta = vcat(theta0[1:end-2], kappa_level_to_cumsum(theta0[end-1:end]))

    data = DataRoyalty(500, theta, L)

    # check indexing is correct
    @test theta_royalty_ρ(data, theta) == theta[1]
    @test theta_royalty_ψ(data, theta) == theta[2]
    @test all(theta_royalty_β(data, theta) .== theta[2 .+ (1:k)])
    @test all(theta_royalty_κ(data, theta) .== theta[end-L+2:end])

    # simulations
    uv = SimulationDraws(M,nobs)

    # let obs = first(first(data)), simi = view(uv,1)
    #     @code_warntype simloglik_royalty!(obs, theta, simi, false)
    # end

    let obs = first(first(data)), simi = view(uv,1)
        print("")
        # @show @benchmark simloglik_royalty!($obs, $theta, $simi, true)
        print("")
    end

    function fg!(grad::AbstractVector, θ::AbstractVector, dograd::Bool=true)
        LL = zero(eltype(θ))
        update_ψ1!(uv, theta_royalty_ρ(data,θ))
        update_dψ1dθρ!(uv, theta_royalty_ρ(data,θ))

        qm = _qm(uv)

        update_xbeta!(data, theta_royalty_β(data,θ))

        for grp in data
            for obs in grp
                fill!(qm, 0)                    # b/c might do other stuff to LLm
                simi = view(uv,_i(grp))
                simloglik_royalty!(obs, θ, simi, dograd)    # update LLm
                LL += logsumexp!(qm) - log(M)
                if dograd
                    grad_simloglik_royalty!(grad, obs, θ, simi)
                end
            end
        end
        return LL
    end

    gradtmp = zeros(size(theta))
    fg!(gradtmp, theta, false)
    fg!(gradtmp, theta, true)

    fd = Calculus.gradient(x -> fg!(gradtmp, x, false), theta)
    fill!(gradtmp, 0)
    fg!(gradtmp, theta, true)
    @test gradtmp ≈ fd
end

flush(stdout)

end # module


text:
# helper functions for defining adjoints
# why be sensible and readable when you can use metaprogramming?

"""symbol or sub-expr replacement in expressions"""
function replace_expr!(e::Expr, pairs...)

  pairs = Dict(pairs)

  if e ∈ keys(pairs)
    return pairs[e]
  end

  for (i, e₀) in enumerate(e.args)
      if e₀ ∈ keys(pairs)
          e.args[i] = pairs[e₀]
      elseif e₀ isa Expr
          replace_expr!(e₀, pairs...)
      end
  end

  return e
end

"""expression case"""
function replace_expr!(e::Symbol, pairs...)

  pairs = Dict(pairs)

  return e ∈ keys(pairs) ? pairs[e] : e
end

"""non-(symbol/expr) cases"""
replace_expr!(e, pairs...) = e


text:
function read_xyz(ifile::String)
    header = Array{String, 1}()
    atom_labels = Array{Array{String, 1}, 1}()
    geoms = Array{Array{Float64, 2}, 1}()
    open(ifile) do io
        for line in eachline(io)
            if isa(tryparse(Int, line), Int)
                # allocate the geometry for this frame
                N = parse(Int, line)
                labels = String[]
                geom = zeros((3, N))
                # store the header for this frame
                head = string(line, readline(io))
                push!(header, head)
                # loop through the geometry storing the vectors and atom labels as you go
                for i = 1:N
                    coords = split(readline(io))
                    push!(labels, coords[1])
                    geom[:,i] = parse.(Float64, coords[2:end])
                end
                push!(atom_labels, labels)
                push!(geoms, geom)
            end
        end
    end
return header, atom_labels, geoms
end

max_stars_repo_name: heindelj/pyMD

max_issues_repo_name: heindelj/pyMD

max_forks_repo_name: heindelj/pyMD


text:
using Blosc2:Lib
@testset "storage" begin
    s = Storage()
    @test s.contiguous == false
    @test isnothing(s.urlpath)
    @test isnothing(s.cparams)
    @test isnothing(s.dparams)

    Blosc2.with_preserved_blosc2_storage(s) do bs
        @test bs.contiguous == false
        @test bs.urlpath == C_NULL
        @test bs.cparams == C_NULL
        @test bs.dparams == C_NULL
        @test bs.io == C_NULL
    end

    s = Storage(contiguous = true,
        urlpath = "test",
        cparams = CompressionParams(Int32),
        dparams = DecompressionParams(nthreads = 2)
    )

    @test s.contiguous == true
    @test s.urlpath == "test"
    @test s.cparams.typesize == 4
    @test s.dparams.nthreads == 2

    Blosc2.with_preserved_blosc2_storage(s) do bs
        @test bs.contiguous == true
        @test unsafe_string(bs.urlpath) == "test"
        @test bs.cparams != C_NULL
        @test bs.dparams != C_NULL
        @test unsafe_load(bs.cparams.typesize) == 4
        @test unsafe_load(bs.dparams.nthreads) == 2
        @test bs.io == C_NULL
    end
end

@testset "superchunk create" begin
    schunk = SChunk(Storage())
    @test clevel(schunk) == 5
    @test typesize(schunk) == 8

    schunk = SChunk(Storage(cparams = CompressionParams(Int16, level = 9)))
    @test clevel(schunk) == 9
    @test typesize(schunk) == 2
end

@testset "schunk append buffer" begin
    rm("test.bl2", force = true, recursive = true)
    schunk = SChunk(Storage(urlpath = "test.bl2"))
    n = 100000
    data = rand(1:1000, n)
    r = GC.@preserve data begin
        unsafe_append_buffer!(schunk, pointer(data), sizeof(data))
    end
    @test r == 1

    data2 = rand(1:1000, floor(Int64, n/2))
    r = GC.@preserve data begin
        unsafe_append_buffer!(schunk, pointer(data2), sizeof(data2))
    end
    @test r == 2

    @test nchunks(schunk) == 2

    res = Vector{Int64}(undef, n)
    r = GC.@preserve res begin
        unsafe_decompress_chunk(schunk, 2, pointer(res), sizeof(res))
    end
    @test r == length(data2)
    @test res[1:r] == data2

    r = GC.@preserve res begin
        unsafe_decompress_chunk(schunk, 1, pointer(res), sizeof(res))
    end
    @test r == n
    @test res == data

    @test_throws BoundsError unsafe_decompress_chunk(schunk, 3, pointer(res), sizeof(res))
    rm("test.bl2", force = true, recursive = true)
end

@testset "schunk unsafe chunk manipulations" begin
    schunk = SChunk(Storage())
    n = 100000
    data = rand(1:1000, n)
    r = GC.@preserve data begin
        unsafe_append_buffer!(schunk, pointer(data), sizeof(data))
    end
    println(chunksize(schunk))
    @test r == 1

    data2 = rand(1:1000, n)
    chunk = compress(data2)
    r = GC.@preserve chunk begin
        unsafe_append_chunk!(schunk, pointer(chunk))
    end
    @test r == 2

    data3 = rand(1:1000, n)
    chunk = compress(data3)
    r = GC.@preserve chunk begin
        unsafe_insert_chunk!(schunk, 2, pointer(chunk))
    end
    @test r == 3

    res = Vector{Int64}(undef, n)
    r = GC.@preserve res begin
        unsafe_decompress_chunk(schunk, 2, pointer(res), sizeof(res))
    end
    @test res == data3

    data4 = rand(1:1000, n)
    chunk = compress(data4)
    r = GC.@preserve chunk begin
        unsafe_update_chunk!(schunk, 2, pointer(chunk))
    end
    @test r == 3
    r = GC.@preserve res begin
        unsafe_decompress_chunk(schunk, 2, pointer(res), sizeof(res))
    end
    @test res == data4

    @test delete_chunk!(schunk, 2) == 2

    r = GC.@preserve res begin
        unsafe_decompress_chunk(schunk, 2, pointer(res), sizeof(res))
    end
    @test res == data2

    chunk2 = unsafe_get_chunk(schunk, 1)
    @test decompress(Int64, chunk2) == data

end


text:
export NeutralBranch

# ExpandingBranch
@compat abstract type AbstractBranch{D<:Domain,R<:Domain}; end

@compat abstract type ExpandingBranch{D<:Domain,R<:Domain} <: AbstractBranch{D,R}; end

Base.summary(b::ExpandingBranch) =  string(typeof(b).name.name)*":"*string(domain(b))*"↦"*string(rangedomain(b)) #branches??
ApproxFun.prectype(b::ExpandingBranch) = prectype(rangedomain(b))
eltype(b::ExpandingBranch) = eltype(rangedomain(b))
# Base.show(io::IO,b::ExpandingBranch) = print(io,typeof(b)) #temporary

@compat struct FwdExpandingBranch{ff,gg,D<:Domain,R<:Domain} <: ExpandingBranch{D,R}
  f::ff
  dfdx::gg
  domain::D
  rangedomain::R
  #  sgn::T
  # function FwdExpandingBranch(fc,dfdxc,dom,ran)
    # @assert all([in(fc(p),∂(ran)) for p in ∂(dom)])
  #   new(fc,dfdxc,dom,ran)
  # end
end
function FwdExpandingBranch(f,dfdx,dom,ran)
  domd = convert(Domain,dom); randm = convert(Domain,ran)
  FwdExpandingBranch{typeof(f),typeof(dfdx),typeof(domd),typeof(randm)}(f,dfdx,domd,randm)
end

unsafe_call(b::FwdExpandingBranch,x) = b.f(x)
unsafe_mapD(b::FwdExpandingBranch,x) = b.dfdx(x)
unsafe_mapP(b::FwdExpandingBranch,x) = (unsafe_call(b,x),unsafe_mapD(b,x))
unsafe_mapinv(b::FwdExpandingBranch,x) = domain_newton(b.f,b.dfdx,x,b.domain,domain_guess(x,b.domain,b.rangedomain))
unsafe_mapinvD(b::FwdExpandingBranch,x) = inv(b.dfdx(unsafe_mapinv(b,x)))
function unsafe_mapinvP(b::FwdExpandingBranch,x)
  vx = unsafe_mapinv(b,x)
  (vx,inv(unsafe_mapD(b,vx)))
end

# RevExpandingBranch

@compat struct RevExpandingBranch{ff,gg,D<:Domain,R<:Domain} <: ExpandingBranch{D,R}
  v::ff
  dvdx::gg
  domain::D
  rangedomain::R
  #   sgn::T
  # function RevExpandingBranch(vc,dvdxc,dom,ran)
  #   new(vc,dvdxc,dom,ran)
  # end
end
function RevExpandingBranch(v,dvdx,dom,ran)
    domd = convert(Domain,dom); randm = convert(Domain,ran)
    RevExpandingBranch{typeof(v),typeof(dvdx),typeof(domd),typeof(randm)}(v,dvdx,domd,randm)
end

unsafe_call(b::RevExpandingBranch,x) = domain_newton(b.v,b.dvdx,x,b.rangedomain,domain_guess(x,b.rangedomain,b.domain))
unsafe_mapD(b::RevExpandingBranch,x) =  inv(b.dvdx(unsafe_call(b,x)))
function unsafe_mapP(b::RevExpandingBranch,x)
  fx = unsafe_call(b,x)
  (fx,inv(unsafe_mapinvD(b,fx)))
end
unsafe_mapinv(b::RevExpandingBranch,x) = b.v(x)
unsafe_mapinvD(b::RevExpandingBranch,x) = b.dvdx(x)
unsafe_mapinvP(b::RevExpandingBranch,x) = (unsafe_mapinv(b,x),unsafe_mapinvD(b,x))

for B in (FwdExpandingBranch,RevExpandingBranch)
  for (M,UNS_M) in ((:mapinv,:unsafe_mapinv),(:mapinvD,:unsafe_mapinvD),(:mapinvP,:unsafe_mapinvP))
    @eval $M(b::$B,x) = begin @assert approx_in(x,rangedomain(b)); $UNS_M(b,x); end # TODO: sort out approx_in stuff
  end

  @eval (b::$B)(x::T) where T<:AbstractInterval = #begin println(leftendpoint(x)); println(rightendpoint(x));
    T(extrema((b(leftendpoint(x)),b(rightendpoint(x))))...) # end
  #@eval (b::$B)(x::Interval) = Interval(sort(b.(extrema(x)))...)
end

mapinv(b::ExpandingBranch,x::T) where T<:AbstractInterval = T(extrema((mapinv(b,leftendpoint(x)),mapinv(b,rightendpoint(x))))...)
# mapinv(b::ExpandingBranch,x::Interval) = Interval(sort(mapinv.(b,extrema(x)))...)

# UNDER CONSTRUCTION: NeutralBranch
@compat struct NeutralBranch
end

# branch constructors

autodiff(f,d) = autodiff_dual(f,ApproxFun.checkpoints(convert(Domain,d)))
autodiff(f::Fun,d) = f'

function autodiff_dual(f,bi)
  fd = FunctionDerivative(f)
  try
    for b in bi
      fd(b)
    end
  catch e
    isa(e,MethodError) && error("To use automatic differentiation, your function must accept DualNumbers")
    throw(e)
  end
  fd
end

function branch(f,dom,ran,diff=autodiff(f,(dir==Forward ? dom : ran));dir=Forward,
                      ftype=typeof(f),difftype=typeof(diff))
  domd = convert(Domain,dom); randm  = convert(Domain,ran);
  dir==Forward ? FwdExpandingBranch{ftype,difftype,typeof(domd),typeof(randm)}(f,diff,domd,randm) :
        RevExpandingBranch{ftype,difftype,typeof(domd),typeof(randm)}(f,diff,domd,randm)
end
branch(f,dom,ran,diff::Nothing; dir =Forward) = branch(f,dom,ran;dir=dir)

@deprecate branch(f,dfdx,dom::Domain,ran::Domain;dir=Forward) branch(f,dom,ran,diff;dir=dir)

for TYP in (:FwdExpandingBranch,:RevExpandingBranch,:NeutralBranch)
  @eval (b::$TYP)(x) = x ∈ b.domain ? unsafe_call(b,x) : error("DomainError: $x ∉ $(b.domain)")
  @eval mapD(b::$TYP,x) = x ∈ b.domain ? unsafe_mapD(b,x) : error("DomainError: $x ∉ $(b.domain)")
  @eval mapP(b::$TYP,x) = x ∈ b.domain ? unsafe_mapP(b,x) : error("DomainError: $x ∉ $(b.domain)")
  @eval domain(b::$TYP) = b.domain
  @eval rangedomain(b::$TYP) = b.rangedomain
end

# transferbranch

function transferbranch_int_edges(x,y,b::ExpandingBranch)
  x ∈ rangedomain(b) && y ∈ rangedomain(b) && (return x,y)
  iv = Interval(x,y)∩rangedomain(b)
  leftendpoint(iv), rightendpoint(iv)
end

function transferbranch(x,b::ExpandingBranch,f)
  x ∉ rangedomain(b) && return zero(typeof(x))
  (v,dvdx) = mapinvP(b,x)
  abs(det(dvdx))*f(v)
end
function transferbranch_int(x,y,b::ExpandingBranch,f)
  x ∉ rangedomain(b) && y ∉ rangedomain(b) && (return zero(typeof(x)))
  x,y = transferbranch_int_edges(x,y,b)
  csf = cumsum(f)
  vy = unsafe_mapinv(b,y); vx = unsafe_mapinv(b,x)
  sgn = sign((vy-vx)/(y-x))
  sgn*(csf(vy)-csf(vx))
end

function transferbranch(x,b::ExpandingBranch,sk::BasisFun)
  x ∉ rangedomain(b) && return zero(typeof(x))
  (v,dvdx) = unsafe_mapinvP(b,x)
  abs(det(dvdx))*getbasisfun(v,sk)
end
function transferbranch_int(x,y,b::ExpandingBranch,sk::BasisFun)
  x ∉ rangedomain(b) && y ∉ rangedomain(b) && (return zero(typeof(x)))
  x,y = transferbranch_int_edges(x,y,b)
  vy = unsafe_mapinv(b,y); unsafe_vx = unsafe_mapinv(b,x)
  sgn = sign((vy-vx)/(y-x))
  sgn*(getbasisfun_int(vy,sk)-getbasisfun_int(vx,sk))
end


text:
include(joinpath(@__DIR__, "../../src/FLS.jl"))

"""
Infinite sequence of a's and b's
"""

inflistschema = @acset LabeledGraph begin
    V = 3
    E = 4
    vlabel = [s1, :d, :l]
    elabel = [:a,:b,:head,:tail]
    src = [1,1,3,3]
    tgt = [2,2,2,3]
end

add_cone!(inflist, @acset LabeledGraph begin
    V = 1
    vlabel = [s1]
end)

lconed =  @acset LabeledGraph begin
    V = 2
    vlabel = [:d, :l]
end

inflist = FLS(:inflist, inflistschema, [
    Cone(c2d, s1,[]),
    Cone(lconed, :l, [1=>:head,2=>:tail])], [])


text:
"""
    edm_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}

Returns a Euclidean distance matrix completion of `A` if it is EDM-completable
using Algorithm 11.1 in [VA15].

# Reference
[VA15] [Chordal Graphs and Semidefinite Optimization](https://www.seas.ucla.edu/~vandenbe/publications/chordalsdp.pdf)
by Lieven Vandenberghe and Martin Andersen
"""
function edm_completion(A::SparseMatrixCSC{Tv, Ti}) where {Tv <: AbstractFloat, Ti <: Integer}
    !issymmetric(A) && error(ArgumentError("A must be symmetric"))
    n = size(A, 1)

    sp = sparsity_pattern(A)
    !is_chordal(A; peo=1:n) && error(ArgumentError("A must have a chordal sparsity pattern"))
    perm, iperm, L = get_chordal_extension(sp; perm=nothing, verbose=false)
    ct = CliqueTree(L)
    order_snds!(ct)
    iperm = invperm(ct.perm)

    !is_edm_completable(A[ct.perm, ct.perm], ct) && error("A is not EDM completable")

    W = Matrix(A[ct.perm, ct.perm])

    for j in (length(ct.snds)-1):-1:1
        vrep_ind = ct.postordering[j]
        vrep = ct.vreps[vrep_ind]

        ν = ct.snds[vrep_ind]
        α = ct.seps[vrep_ind]
        η = filter(x->(!(x in ν) && !(x in α)), vrep+1:n)

        aj = first(α)

        # See eq (11.4)
        # (I - 1*eₖ') * A * (I - eₖ*1') = A - A[k,k]*1*1' - 1*ãₖ' - aₖ*1'
        Yνν = -0.5*W[ν, ν] .- 0.5*W[aj, aj]
        Yνν .+= 0.5 .* W[ν, aj]
        Yνν .+= 0.5 .* W[aj, ν]'

        Yαα = -0.5*W[α, α] .- 0.5*W[aj, aj]
        Yαα .+= 0.5 .* W[α, aj]
        Yαα .+= 0.5 .* W[aj, α]'

        Yηη = -0.5*W[η, η] .- 0.5*W[aj, aj]
        Yηη .+= 0.5 .* W[η, aj]
        Yηη .+= 0.5 .* W[aj, η]'

        Yηα = -0.5*W[η, α] .- 0.5*W[aj, aj]
        @. Yηα += 0.5 * W[η, aj]
        @. Yηα += 0.5 * W[aj, α]'

        Yαν = -0.5*W[α, ν] .- 0.5*W[aj, aj]
        @. Yαν += 0.5 * W[α, aj]
        @. Yαν += 0.5 * W[aj, ν]'

        d, V = LAPACK.syev!('V', 'L', Yαα)
        λmax = maximum(d)
        Yαα_pinv = V * Diagonal(map(x->(x > λmax*eps() ? 1.0/x : 0.0), d)) * V'
        cache = Yαα_pinv * Yαν
        @views mul!(W[η, ν], Yηα, cache)
        W[η, ν] .*= -2
        W[η, ν] .+= diag(Yηη)
        W[η, ν] .+= diag(Yνν)'

        # W[η, ν] = -2*Yηα*Yαα_pinv*Yαν + ones(length(η))*diag(Yνν)' + diag(Yηη)*ones(length(ν))'
        # BLAS.ger!(1.0, ones(length(η)), diag(Yνν), W[η, ν])
        # BLAS.ger!(1.0, diag(Yηη), ones(length(ν)), W[η, ν])

        # μ = sparsevec([aj], [1.0], n)
        # tmp  = -0.5*(I - ones(n)*μ')*W*(I - μ*ones(n)')
        # rtol = sqrt(eps(real(float(one(eltype(tmp[α, α]))))))
        # W[η, ν] .= -2*tmp[η, α]*pinv(tmp[α, α]; rtol=rtol)*tmp[α, ν] +
        #           ones(length(η))*diag(tmp[ν,ν])' + diag(tmp[η,η])*ones(length(ν))'

        W[ν, η] .= W[η, ν]'
        # display(Matrix(W[iperm, iperm] .* .!(sp .> 0)))
        # break
    end

    return W[iperm, iperm]

end

"""
    is_edm_completable(A::SparseMatrixCSC, ct::CliqueTree)

Checks if `A` is EDM completable. Requires the associated CliqueTree `ct`.
"""
function is_edm_completable(A::SparseMatrixCSC, ct::CliqueTree)
    cliques = get_cliques(ct)
    for c in cliques
        !is_edm(A[c,c]) && return false
    end
    return true
end

"""
    is_edm(M::AbstractMatrix)

Checks if `M` is a Euclidean Distance Matrix. Note: a matrix is a EDM iff it is
negative semidefinite on the subspace on the subspace orthogonal to the all ones
vector.
"""
function is_edm(M::AbstractMatrix)
    (!iszero(diag(M)) || !issymmetric(M)) && return false
    n = size(M, 1)
    Q = spdiagm(n, n-1, 0 => ones(n-1), -1 => -ones(n-1))

    return isposdef(-Symmetric(Q'*M*Q) + 1e-10*I)
end


text:
module RandomODE
    using Distributions
    using LinearAlgebra
    using Random
    using SparseArrays
    using FFTW

    include("structures.jl")
    export RandomODEProblem, DirichletBC

    include("constructors.jl")
    export RandomODEProblem

    include("sampler.jl")
    export sample_field

    include("assembly.jl")
    export assemble_system!

    include("solve.jl")
    export solve!, solve_ensemble!
end


text:
using OutlierDetectionInterface
using Test

@testset "OutlierDetectionInterface.jl" begin
    # Make sure the import worked
end


text:
using Dates
using Setfield

using TeaFiles.Header: create_item_namedtuple, create_metadata, create_metadata_julia_time

@testset "create" begin
    @testset "create_metadata" begin
        @testset "non bits type" begin
            struct Invalid
                a::Array
                b::UInt8
            end
            @test !isbitstype(Invalid)
            @test_throws ArgumentError create_metadata(Invalid)
        end

        @testset "invalid duration" begin
            struct Moo2 a::Int end
            # should not throw
            create_metadata(Moo2)
            create_metadata(Moo2; tick_duration=Day(1))
            @test_throws ArgumentError create_metadata(Moo2; tick_duration=Day(2))
        end

        @testset "example" begin
            metadata = create_metadata(
                Tick;
                time_fields=[:time],
                content_description="ACME prices",
                name_values=NameValue[NameValue("decimals", Int32(2))]
            )
            reference_metadata = _get_example_metadata()
            @test metadata == reference_metadata
        end

        @testset "example namedtuple" begin
            namedtuple_tick = NamedTuple{
                (:time, :price, :volume),
                Tuple{Int64, Float64, Int64}
            }
            metadata = create_metadata(
                namedtuple_tick;
                time_fields=[:time],
                content_description="ACME prices",
                name_values=NameValue[NameValue("decimals", Int32(2))]
            )
            reference_metadata = _get_example_metadata()
            item_section = get_section(reference_metadata, ItemSection)
            item_section = @set item_section.item_name = "NamedTuple"
            _replace_section!(reference_metadata, item_section)
            @test metadata == reference_metadata
        end
    end

    @testset "create_metadata_julia_time" begin
        @testset "example" begin
            metadata = create_metadata_julia_time(
                DateTimeTick;
                content_description="ACME prices",
                name_values=NameValue[NameValue("decimals", Int32(2))]
            )
            reference_metadata = _get_example_datetime_metadata()
            @test metadata == reference_metadata
        end
    end

    @testset "create_item_namedtuple" begin
        metadata = _get_example_metadata()
        namedtuple = create_item_namedtuple(metadata)

        namedtuple_tick = NamedTuple{
            (:time, :price, :volume),
            Tuple{Int64, Float64, Int64}
        }
        @test namedtuple == namedtuple_tick

        metadata = _get_example_datetime_metadata()
        namedtuple = create_item_namedtuple(metadata)

        namedtuple_datetime_tick = NamedTuple{
            (:time, :price, :volume),
            Tuple{DateTime, Float64, Int64}
        }
        @test namedtuple == namedtuple_datetime_tick
    end
end


text:
using DataFrames
using CSV
using CategoricalArrays

datafile(filename) = joinpath(Pkg.dir("SoftComputingFinal"), "data", filename)

# Glass identification:
# 1. Id number: 1 to 214
# 2. RI: refractive index
# 3. Na: Sodium (unit measurement: weight percent in corresponding oxide, as are attributes 4-10)
# 4. Mg: Magnesium
# 5. Al: Aluminum
# 6. Si: Silicon
# 7. K: Potassium
# 8. Ca: Calcium
# 9. Ba: Barium
# 10. Fe: Iron
# 11. Type of glass: (class attribute)
# -- 1 building_windows_float_processed
# -- 2 building_windows_non_float_processed
# -- 3 vehicle_windows_float_processed
# -- 4 vehicle_windows_non_float_processed (none in this database)
# -- 5 containers
# -- 6 tableware
# -- 7 headlamps

function load_glass()
    glass = CSV.read(datafile("glass-identification.data"),
                     header = [:id, :ri, :na, :mg, :al, :si, :k, :ca, :ba, :fe, :type],
                     types = Dict(11 => String))

    categorical!(glass, :type)
    conversions = map(=>, string.(1:7),
                      ["building_windows_float_processed",
                       "building_windows_non_float_processed",
                       "vehicle_windows_float_processed",
                       "vehicle_windows_non_float_processed",
                       "containers",
                       "tableware",
                       "headlamps"])
    recode!(glass[:type], conversions...)

    delete!(glass, :id)

    return glass
end

# Ionosphere
# First 34 are continuous values
# The 35th attribute is either "good" or "bad" evidence

function load_ionosphere()
    ionosphere = CSV.read(datafile("ionosphere.data"),
                          header = [:id, [Symbol("x", i) for i = 1:34]..., :evidence],
                          types = Dict(2 => Float64, 3 => Float64))
    recode!(ionosphere[:evidence], "g" => "good", "b" => "bad")
    delete!(ionosphere, :id)

    return ionosphere
end

# Image segmentation
# 19 continuous attributes:
#     1.  region-centroid-col:  the column of the center pixel of the region.
#     2.  region-centroid-row:  the row of the center pixel of the region.
#     3.  region-pixel-count:  the number of pixels in a region = 9.
#     4.  short-line-density-5:  the results of a line extractoin algorithm that
#          counts how many lines of length 5 (any orientation) with
#          low contrast, less than or equal to 5, go through the region.
#     5.  short-line-density-2:  same as short-line-density-5 but counts lines
#          of high contrast, greater than 5.
#     6.  vedge-mean:  measure the contrast of horizontally
#          adjacent pixels in the region.  There are 6, the mean and
#          standard deviation are given.  This attribute is used as
#         a vertical edge detector.
#     7.  vegde-sd:  (see 6)
#     8.  hedge-mean:  measures the contrast of vertically adjacent
#           pixels. Used for horizontal line detection.
#     9.  hedge-sd: (see 8).
#     10. intensity-mean:  the average over the region of (R + G + B)/3
#     11. rawred-mean: the average over the region of the R value.
#     12. rawblue-mean: the average over the region of the B value.
#     13. rawgreen-mean: the average over the region of the G value.
#     14. exred-mean: measure the excess red:  (2R - (G + B))
#     15. exblue-mean: measure the excess blue:  (2B - (G + R))
#     16. exgreen-mean: measure the excess green:  (2G - (R + B))
#     17. value-mean:  3-d nonlinear transformation
#          of RGB. (Algorithm can be found in Foley and VanDam, Fundamentalsa
#          of Interactive Computer Graphics)
#     18. saturatoin-mean:  (see 17)
#     19. hue-mean:  (see 17)
# Class Distribution:
#    Classes:  brickface, sky, foliage, cement, window, path, grass.

function load_segmentation()
    segmentation = CSV.read(datafile("image-segmentation.data"),
                            header = [:id,
                                      :region_centroid_col,
                                      :region_centroid_row,
                                      :region_pixel_count,
                                      :short_line_density_5,
                                      :short_line_density_2,
                                      :vedge_mean,
                                      :vedge_sd,
                                      :hedge_mean,
                                      :hedge_sd,
                                      :intensity_mean,
                                      :rawred_mean,
                                      :rawblue_mean,
                                      :rawgreen_mean,
                                      :exred_mean,
                                      :exblue_mean,
                                      :exgreen_mean,
                                      :value_mean,
                                      :saturation_mean,
                                      :hue_mean,
                                      :image],
                            types = Dict(4 => Float64),
                            transforms = Dict(21 => lowercase ∘ string))
    categorical!(segmentation, :image)
    delete!(segmentation, :id)

    return segmentation
end

# Linearly separable 2D data set from https://github.com/cuekoo/Binary-classification-dataset

function load_testdata()
    testdata = CSV.read(datafile("testdata.data"),
                        header = [:x, :y, :class])
    categorical!(testdata, :class)
end

function check_loading(df)
    display(names(df))
    display(eltypes(df))
    display(levels(df[end]))
end

max_stars_repo_name: phipsgabler/soft-computing-final

max_issues_repo_name: phipsgabler/soft-computing-final

max_forks_repo_name: phipsgabler/soft-computing-final


text:
module uwuify
const CHAR_REPLACEMENTS = [
    "l"=>"w","L"=>"W",
    "r"=>"w","R"=>"W"
]
export uwu
function uwu(s::AbstractString)
    reduce(replace, CHAR_REPLACEMENTS,init=s)
end

export @uwu_str
macro uwu_str(s)
    uwu(s)
end
end


text:
function closeness_centrality(g::AbstractGraph,
    distmx::AbstractMatrix=weights(g);
    normalize=true)::Vector{Float64}

    n_v = Int(nv(g))
    closeness = SharedVector{Float64}(n_v)

    @sync @distributed for u in vertices(g)
        if degree(g, u) == 0     # no need to do Dijkstra here
            closeness[u] = 0.0
        else
            d = LightGraphs.dijkstra_shortest_paths(g, u, distmx).dists
            δ = filter(x -> x != typemax(x), d)
            σ = sum(δ)
            l = length(δ) - 1
            if σ > 0
                closeness[u] = l / σ
                if normalize
                    n = l * 1.0 / (n_v - 1)
                    closeness[u] *= n
                end
            end
        end
    end
    return sdata(closeness)
end


text:
module SpringCollab2020MusicAreaChangeTrigger

using ..Ahorn, Maple

@mapdef Trigger "SpringCollab2020/MusicAreaChangeTrigger" MusicAreaChangeTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
	musicParam1::String="", musicParam2::String="", musicParam3::String="", enterValue::Number=1.0, exitValue::Number=0.0)

const placements = Ahorn.PlacementDict(
    "Music Area Change Trigger (Spring Collab 2020)" => Ahorn.EntityPlacement(
        MusicAreaChangeTrigger,
        "rectangle"
    )
)

end


text:
export AS

module AS
    const Generic               = 0 # No Generic address space?
    const Device                = 1 # Checked
    const Constant              = 2 # Checked
    const ThreadGroup           = 3 # Checked
    const Thread                = 4 # Ends up same as Device?
    const ThreadGroup_ImgBlock  = 5 # Like ThreadGroup but only accessible from
    const Ray                   = 6
end


text:
"""
    VariableLinearAdvectionNonperiodicSemidiscretization(D, Di, a, split_form,
                                                         left_bc, right_bc)

A semidiscretization of the linear advection equation
    \$\\partial_t u(t,x) + \\partial_x ( a(x) u(t,x) ) = 0\$
with boundary conditions `left_bc(t)`, `right_bc(t)`.

`D` is an SBP derivative operator, `Di` an associated dissipation operator
or `nothing`, `a(x)` the variable coefficient, and `split_form::Union{Val(false), Val(true)}`
determines whether the canonical split form or the conservative form should be
used.
"""
struct VariableLinearAdvectionNonperiodicSemidiscretization{T,Derivative<:AbstractDerivativeOperator{T},
                                                            Dissipation,
                                                            SplitForm<:Union{Val{false}, Val{true}},
                                                            LeftBC, RightBC} <: AbstractSemidiscretization
    derivative::Derivative
    dissipation::Dissipation
    a::Vector{T}
    tmp1::Vector{T}
    tmp2::Vector{T}
    split_form::SplitForm
    left_bc::LeftBC
    right_bc::RightBC

    function VariableLinearAdvectionNonperiodicSemidiscretization(derivative::Derivative, dissipation::Dissipation, afunc, split_form::SplitForm, left_bc::LeftBC, right_bc::RightBC) where {T, Derivative<:AbstractDerivativeOperator{T}, Dissipation, SplitForm<:Union{Val{false}, Val{true}}, LeftBC, RightBC}
        if dissipation !== nothing
            @argcheck size(derivative) == size(dissipation) DimensionMismatch
            @argcheck grid(derivative) == grid(dissipation) ArgumentError
        end
        N = size(derivative, 2)
        a = compute_coefficients(afunc, derivative)
        tmp1 = Array{T}(undef, N)
        tmp2 = Array{T}(undef, N)
        new{T,Derivative,Dissipation,SplitForm,LeftBC,RightBC}(derivative, dissipation, a, tmp1, tmp2, split_form, left_bc, right_bc)
    end
end

function Base.show(io::IO, semi::VariableLinearAdvectionNonperiodicSemidiscretization)
    if get(io, :compact, false)
        print(io, "Semidiscretization of the linear advection equation (periodic)")
    else
        print(io, "Semidiscretization of the linear advection equation\n")
        print(io, "  \$ \\partial_t u(t,x) + \\partial_x ( a(x) u(t,x) ) = 0 \$ \n")
        print(io, "with nonperiodic boundaries using")
        if semi.split_form == Val{true}()
            print(io, " a split form and: \n")
        else
            print(io, " no split form and: \n")
        end
        print(io, semi.derivative)
        print(io, semi.dissipation)
    end
end

function godunov_flux_variablelinearadvection(uₗ::T, uᵣ::T, a::T) where {T<:Real}
    ifelse(a > 0, a*uₗ, a*uᵣ)
end

function (disc::VariableLinearAdvectionNonperiodicSemidiscretization)(du, u, p, t)
    @unpack a, tmp1, tmp2, derivative, dissipation, split_form, left_bc, right_bc = disc
    @boundscheck begin
        @argcheck length(u) == length(a)
        @argcheck length(u) == length(tmp1)
        @argcheck length(u) == length(tmp2)
        @argcheck length(u) == length(du)
    end

    # volume terms
    if typeof(split_form) <: Val{true}
        m1_2 = -one(eltype(u)) / 2

        ## a * D * u
        mul!(tmp1, derivative, u)
        @. du = m1_2 * a * tmp1
        ## u * D * a
        mul!(tmp1, derivative, a)
        @. du += m1_2 * u * tmp1
        ## D * a*u
        @. tmp2 = a*u
        mul!(tmp1, derivative, tmp2)
        @. du += m1_2 * tmp1
    else
        @. tmp2 = a * u
        mul!(tmp1, derivative, tmp2)
        @. du = -tmp1
    end

    # dissipation
    if dissipation !== nothing
        mul!(tmp1, dissipation, u)
        @. du += tmp1
    end

    # boundary conditions via Godunov's flux
    @inbounds fnum_left = godunov_flux_variablelinearadvection(left_bc(t), u[1], a[1])
    @inbounds du[1] += (fnum_left - a[1]*u[1]) / left_boundary_weight(derivative)
    @inbounds fnum_right = godunov_flux_variablelinearadvection(u[end], right_bc(t), a[end])
    @inbounds du[end] -= (fnum_right - a[end]*u[end]) / right_boundary_weight(derivative)

    nothing
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -3.1544947018849654 <= q <= 4.066926790342683)

                     add_NL_constraint(m, :(log(1 + exp(-0.7544717811096784 + -0.6289249660430758*log(1 + exp(-0.401848603561894 + 0.9879850116297502*log(1 + exp(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5]))) + -0.6833595346541679*log(1 + exp(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))))) + -0.656030629474813*log(1 + exp(-0.9665555719654946 + 0.014645656858022438*log(1 + exp(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5]))) + -0.5583055948877655*log(1 + exp(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))))))) + log(1 + exp(-0.3137145052115069 + -0.6627216812930383*log(1 + exp(-0.401848603561894 + 0.9879850116297502*log(1 + exp(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5]))) + -0.6833595346541679*log(1 + exp(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))))) + -0.39516873731093893*log(1 + exp(-0.9665555719654946 + 0.014645656858022438*log(1 + exp(0.7471741672319268 + 0.1562161389212755*$(x[1]) + -0.06104000939937393*$(x[2]) + -0.24396039580062512*$(x[3]) + -0.35002310972719863*$(x[4]) + -0.47890036372503086*$(x[5]))) + -0.5583055948877655*log(1 + exp(0.07659783311726898 + 0.1400377608602934*$(x[1]) + -0.5920144548543678*$(x[2]) + -0.3715997206073691*$(x[3]) + -0.08423829417739581*$(x[4]) + -0.48744433011047184*$(x[5]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# AXIS
#

mutable struct Axis
    prefix::String # x, y, x2, y2, z
    #
    ticks    ::Ticks          # ticks of the axis
    textstyle::TextStyle      # parent textstyle of axis
    title    ::Option{Title}  # title of the axis
    base     ::Option{F64}    # scale font and ticks
    lwidth   ::Option{F64}    # width of the axis spine
    min      ::Option{F64}    # minimum span of the axis
    max      ::Option{F64}    # maximum span of the axis
    off      ::Bool           # if true, axis is not shown
    log      ::Bool           # log scale
end
Axis(p) = default(Axis, p)

#
# AXES (set of 2 or more axis)
#

abstract type Axes end

mutable struct Axes2D <: Axes
    parent  ::String # id of the parent figure
    xaxis   ::Axis
    x2axis  ::Axis
    yaxis   ::Axis
    y2axis  ::Axis
    drawings::Vector{Drawing2D}
    objects ::Vector{Object2D}
    #
    title ::Option{Title}
    size  ::Option{T2F}     # (width cm, height cm)
    legend::Option{Legend}
    origin::Option{T2F}     # related to layout
    math  ::Bool            # axis crossing (0, 0)
    off   ::Bool
    scale ::String
end
Axes2D(p=""; parent=p) = default(Axes2D,
    parent, Axis("x"), Axis("x2"), Axis("y"), Axis("y2"),
    Drawing2D[], Object2D[]
)

mutable struct Axes3D <: Axes
end

#
# @with_kw mutable struct Axes3D <: Axes
#     parent::String
#     # --
#     xaxis::Axis = Axis("x"; min=0, max=1) # NOTE color = ticks, not spine if box
#     yaxis::Axis = Axis("y"; min=0, max=1) # if nobox, then spine
#     zaxis::Axis = Axis("z"; min=0, max=1)
#     # ---
#     drawings::Vector{Drawing3D} = Vector{Drawing3D}()
#     objects ::Vector{Object3D}  = Vector{Object3D}()
#     # ---
#     title::Option{Title} = ∅
#     size ::T2F           = (20.,10.) # box size
#     # cube
#     nocube   ::Bool      = false # XXX if true, then xaxis have an expressed linestyle
#     cubedims ::T3F       = (20.,20.,10.) # cube sides x,y,z
#     linestyle::LineStyle = LineStyle() # XXX only lstyle, color, see set
#     origin   ::Option{T2F} = ∅
#     rotate   ::Option{T2F} = ∅
#     off      ::Bool = false # do not show
# end
# Axes3D(p) = default(Axes3D,
#     p, Axis("x")
# )


text:
export relatedness, pairwise_relatedness, kinship

#=TODO

    Solve dyadic optimization issues
        combined backtracking failure - need to isolate and find cause
            #N111 N83
            #N63 N64
        ~1/3 - 1/5 not with optimal solution in cats dataset
    Multithread

=#

#= Would be good to include

    Implement alternative relatedness metrics

    Warning if a not implemented (or typo) of method included

    Streamline output

=#

"""
    pr_l_s(x::Tuple, y::Tuple, alleles::Dict)
Calculate the probability of observing the particular allele state given each of
the 9 Jacquard Identity States for a single locus to create Table 1 from
Milligan 2003.
"""
function pr_l_s(x::Tuple, y::Tuple, alleles::Dict)
    #= Example
    cats = nancycats()
    cat1=get_genotype(cats, sample = "N100", locus = "fca23")
    cat2=get_genotype(cats, sample = "N111", locus = "fca23")
    allele = allele_freq(cats.loci.fca23)
    pr_l_s(cat1, cat2, allele)
    =#
    #=
    Calculate Pr(Li | Sj)
    If the allele identity falls into this class (L1-L9), generate the
    probabilities of it belonging to each of the different classes and
    return that array of 9 distinct probabilities
    =#

    ## class L1 -  AᵢAᵢ AᵢAᵢ ##
    if x[1] == x[2] == y[1] == y[2]
        x = string.(x)
        y = string.(y)

        p = alleles[x[1]]
        [p, p^2, p^2, p^3, p^2, p^3, p^2, p^3, p^4]

    ## class L2 - AᵢAᵢ AⱼAⱼ ##
    elseif (x[1] == x[2]) & (y[1] == y[2]) & (x[1] != y[1])
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[y[1]])
        [0, prod(p), 0, prod(p) * p[2], 0, prod(p) * p[1], 0, 0, prod(p) * prod(p)]

    ## class L3a - AᵢAᵢ AᵢAⱼ ## - has issues because of allele order
    elseif ((x[1] == x[2] == y[1]) & (x[1] != y[2]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[y[2]])
        [0, 0, prod(p), 2 * prod(p) * p[1], 0, 0, 0, prod(p) * p[1], 2 * prod(p) * p[1]^2]

    ## class L3b - AᵢAᵢ AⱼAᵢ ## - has issues because of allele order
    elseif ((x[1] == x[2] == y[2]) & (x[1] != y[1]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[y[1]])
        [0, 0, prod(p), 2 * prod(p) * p[1], 0, 0, 0, prod(p) * p[1], 2 * prod(p) * p[1]^2]

    ## class L4 - AᵢAᵢ AⱼAₖ ##
    elseif (x[1] == x[2]) & (y[1] != y[2]) & (x[1] != y[1]) & (x[1] != y[2])
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[y[1]], alleles[y[2]])
        [0, 0, 0, 2 * prod(p), 0, 0, 0, 0, 2 * prod(p) * p[1]]

    ## L5a - AiAj AiAi ## - has issues because of allele order
    elseif ((x[1] == y[1] == y[2]) & (x[1] != x[2]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, prod(p), 2 * prod(p) * p[1], 0, prod(p) *p[1], 2 * prod(p) * p[1]^2]

    ## L5b - AjAi AiAi ## - has issues because of allele order
    elseif (x[2] == y[1] == y[2] & (x[1] != x[2]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[2]], alleles[x[1]])
        [0, 0, 0, 0, prod(p), 2 * prod(p) * p[1], 0, prod(p) *p[1], 2 * prod(p) * p[1]^2]

    ## L6 - AjAk AiAi ##
    elseif (x[1] != x[2]) & (y[1] == y[2]) & (x[1] != y[1]) & (x[2] != y[1])
        x = string.(x)
        y = string.(y)

        p = (alleles[y[1]], alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, 0, 2 * prod(p), 0, 0, 2 * prod(p) * p[1]]

    ## L7 - AiAj AiAj ##
    elseif (x[1] == y[1]) & (x[2] == y[2]) & (x[1] != x[2])
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[x[2]])
        [0, 0, 0, 0, 0, 0, 2 * prod(p), prod(p) * sum(p), 4 * prod(p) * prod(p)]

    ## L8a - AiAj AiAk ##  - has issues because of allele order
    elseif ((x[1] == y[1]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[2] != y[2]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[x[2]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8b - AjAi AkAi ##  - has issues because of allele order
    elseif ((x[2] == y[2]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[1]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[2]], alleles[x[1]], alleles[y[1]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8c - AjAi AiAk ##  - has issues because of allele order
    elseif ((x[2] == y[1]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[2]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[2]], alleles[x[1]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L8d - AiAj AkAi ##  - has issues because of allele order
    elseif ((x[1] == y[2]) & (x[1] != x[2]) & (y[1] != y[2]) & (x[1] != y[1]))
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[x[2]], alleles[y[1]])
        [0, 0, 0, 0, 0, 0, 0, prod(p), 4 * prod(p) * p[1]]

    ## L9 - AiAj AkAl ##
    elseif (x[1] != x[2]) & (x[1] != y[1]) & (x[1] != y[2]) & (x[2] != y[1]) & (x[2] != y[2]) & (y[1] != x[2])
        x = string.(x)
        y = string.(y)

        p = (alleles[x[1]], alleles[x[2]], alleles[y[1]], alleles[y[2]])
        [0, 0, 0, 0, 0, 0, 0, 0, 4 * prod(p)]
    else
        [-9, -9, -9, -9, -9, -9, -9, -9, -9]
    end
end

"""
    all_loci_Pr_L_S(data::PopObj, ind1::String, ind2::String, alleles::Dict)
Calculate the probability of observing the particular allele state given each of
the 9 Jacquard Identity States for all loci Creates Table 1 from Milligan 2002
"""
function all_loci_Pr_L_S(data::PopData, ind1::String, ind2::String, alleles::Dict)
    #Need to skip loci where one or both individuals have missing data
    Pr_L_S = []
    for locus in String.(names(data.loci))
        #Extract the pair of interest's genotypes
        gen1 = get_genotype(data, sample = ind1, locus = locus)
        gen2 = get_genotype(data, sample = ind2, locus = locus)

        if gen1 !== missing && gen2 !== missing
            tmp = pr_l_s(gen1, gen2, alleles[locus])
            return tmp
            push!(Pr_L_S, tmp)
        end
    end
    return transpose(hcat(Pr_L_S...))
end

#Pr_L_S_inbreeding = all_loci_Pr_L_S(ind1, ind2, data, allele_frequencies)

"""
    no_inbreeding(Pr_L_S::LinearAlgebra.Transpose{Float64,Array{Float64,2}})
Remove Jacquard States which can only be the result of inbreeding
"""
function no_inbreeding(Pr_L_S::Transpose{Float64,Array{Float64,2}})
    for i in 1:6
        Pr_L_S[:,i] = 0 .* Pr_L_S[:,i]
    end
    return Pr_L_S
end

#Pr_L_S_noinbreeding = dyadic_ML(data, allele_frequencies) |> no_inbreeding

## Calculate Δ coefficients
#Need to either maximize this sum or use it as the likelihood in a bayesian model and sample from the posterior.
#currently only maximum likelihood optimization

"""
    Δ_optim(Pr_L_S::Transpose{Float64,Array{Float64,2}}, verbose::Bool)
Takes the probability of the allelic state given the identity by descent from
all available loci (either allowing for inbreeding or not) and calculated the
maximum likelihood Δ coefficients
"""
function Δ_optim(Pr_L_S::Transpose{Float64,Array{Float64,2}}, verbose::Bool = true)
    #Δ is what needs to be optimized
    #consist of 9 values between 0 and 1 which must also sum to 1
    #is then used to calculate relatedness

    Δ = Variable(8)
    #problem = maximize(sum(log(Pr_L_S * vcat(1 - sum(Δ), Δ))))
    problem = maximize(sum(Pr_L_S * vcat(1 - sum(Δ), Δ)))
    problem.constraints += 0 <= 1 - sum(Δ)
    problem.constraints += 1 - sum(Δ) <= 1
    problem.constraints += 0 <= Δ[1:8]
    problem.constraints += Δ[1:8] <= 1

    #shifted from actual relatedness calculations because the 1 - sum(Δ) goes at beginning
    # problem.constraints += 2 * ((1 - sum(Δ)) + 0.5 * (Δ[2] + Δ[4] + Δ[6]) + 0.25 * Δ[7]) <= 1
    # problem.constraints += 0 <= 2 * ((1 - sum(Δ)) + 0.5 * (Δ[2] + Δ[4] + Δ[6]) + 0.25 * Δ[7])

    Convex.solve!(problem, ECOSSolver(verbose = verbose, maxit = 100), verbose = verbose) #maxit=100,
    #Convex.solve!(problem, ECOSSolver(verbose = verbose, maxit=100, feastol=5e-6, reltol = 1e-3, reltol_inacc = 5e-2), verbose = verbose)
    #Convex.solve!(problem, SCSSolver(verbose = verbose, max_iters = 100), verbose = verbose)

    vcat(1-sum(Δ.value), Δ.value), problem.status
    # Should probably include some output that confirms that it did in fact
    # converge and/or use multiple random starts to confirm not a local maxima
end

# JuMP variant
function Δ_optim(Pr_L_S::Transpose{Float64,Array{Float64,2}}, verbose::Bool = true)
    #Δ is what needs to be optimized
    #consist of 9 values between 0 and 1 which must also sum to 1
    #is then used to calculate relatedness

    jacquard = Model(ECOS.Optimizer)
    @variable(jacquard, Δ[i=1:9])
    @objective(jacquard, Max, sum(log(Pr_L_S * Δ)))
    @constraint(jacquard, con, sum(Δ) == 1)
    @constraint(jacquard, con_bounds, 0 .<= Δ .<= 1)

end
## Calculate theta and r
"""
    relatedness_dyadicML(Δ::Array{Float64,2})
Takes the Δ coefficents (with or without inbreeding allowed) and calculates the coefficient of relatedness
"""
function relatedness_from_Δ(Δ::Array{Float64,2})
    θ = Δ[1] + 0.5 * (Δ[3] + Δ[5] + Δ[7]) + 0.25 * Δ[8]
    2 * θ
end

#relatedness_dyadicML(Δ_inbreeding)
#relatedness_dyadicML(Δ_noinbreeding)
#Relatedness R package appears to have a bug. When allow.inbreeding = TRUE the relatedness value is the same as when I assume no inbreeding
#when you set allow.inbreeding = FALSE then the relatedness calculated is the same as when I assume there is inbreeding

"""
    dyadicML_relatedness(data::PopObj, ind1::String, ind2::String; alleles::Dict, inbreeding::Bool = true, verbose::Bool = true)
Calculates the dyadic maximum likelihood relatedness using all available loci following
Milligan 2002 dyadic relatedness - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1462494/pdf/12663552.pdf

Inbreeding can either be assumed not to occur (inbreeding = false) or be
included in the calculation of the relatedness coefficient Verbose controls the
verbosity of the optimization process to find the Maximum likelihood Δ coefficents
"""
function dyadicML_relatedness(data::PopObj, ind1::String, ind2::String; alleles::Dict, inbreeding::Bool = true, verbose::Bool = true)

    Pr_L_S = all_loci_Pr_L_S(data, ind1, ind2, alleles)

    if !inbreeding
        no_inbreeding(Pr_L_S)
    end

    Δ,convergence = Δ_optim(Pr_L_S, verbose)

    return relatedness_from_Δ(Δ), Δ, convergence
end

"""
    qg_relatedness(data::PopObj, ind1::String, ind2::String; alleles::Dict)
Calculates the moments based estimator of pairwise relatedness developed by Queller & Goodnight (1989).
- Bases allele frequencies on entire population
- Inbreeding can only be assumed not to exist.
See equation 3 in: https://www.nature.com/articles/hdy201752 for variant of estimator used
"""
function qg_relatedness(data::PopData, ind1::String, ind2::String; alleles::T) where T <: NamedTuple

    #NEED TO CHECK TO CONFIRM EQUATIONS

    n1 = 0.0
    n2 = 0.0
    d1 = 0.0
    d2 = 0.0

    for loc in loci(data)
        #Extract the pair of interest's genotypes
        gen1 = get_genotype(data, sample = ind1, locus = loc)
        gen2 = get_genotype(data, sample = ind2, locus = loc)

        #Skip missing
        if gen1 !== missing && gen2 !== missing
            a,b = gen1
            c,d = gen2
            sym_loc = Symbol(loc)

            n1 += ((a == c) + (a == d) + (b == c) + (b == d)) - 2 * (alleles[sym_loc][a] + alleles[sym_loc][b])
            n2 += ((a == c) + (a == d) + (b == c) + (b == d)) - 2 * (alleles[sym_loc][c] + alleles[sym_loc][d])

            d1 += (2 * (1 + (a==b) - alleles[sym_loc][a] - alleles[sym_loc][b]))
            d2 += (2 * (1 + (c==d) - alleles[sym_loc][c] - alleles[sym_loc][d]))
        end
    end
    return (n1/d1 + n2/d2)/2.0
end

"""
    ritland_relatedness(data::PopObj, ind1::String, ind2::String; alleles::Dict)
Calculates the moments based estimator of pairwise relatedness proposed by Li and Horvitz (1953) and implemented/made popular by Ritland (1996).
- Bases allele frequencies on entire population
- Inbreeding can only be assumed not to exist.
See equation 7 in: https://www.nature.com/articles/hdy201752 for variant of estimator used
Ritland original citation: https://www.cambridge.org/core/journals/genetics-research/article/estimators-for-pairwise-relatedness-and-individual-inbreeding-coefficients/9AE218BF6BF09CCCE18121AA63561CF7
"""
function ritland_relatedness(data::PopData, ind1::String, ind2::String; alleles::T) where T <: NamedTuple

    #NEED TO CHECK TO CONFIRM EQUATIONS

    n = 0.0
    d = 0.0
    for loc in loci(data)
        #Extract the pair of interest's genotypes
        gen1 = get_genotype(data, sample = ind1, locus = loc)
        gen2 = get_genotype(data, sample = ind2, locus = loc)

        #Skip missing
        if gen1 !== missing && gen2 !== missing
            a,b = gen1
            c,d = gen2
            sym_loc = Symbol(loc)

            A = ((alleles[sym_loc] |> length) - 1)

            R = 0.0
            for i in keys(alleles[sym_loc])
                R += ((((a == i) + (b == i)) * ((c == i) + (d == i))) / (4 * alleles[sym_loc][i])) #Individual locus relatedness value (eq 7 in paper)
            end
            R = (2 / A) * (R - 1)

            n += (R * A) #numerator for weighted combination of loci
            d += A #denominator for weighted combination of loci
        end
    end
    return (n / d)
end

"""
    lr_relatedness(data::PopObj, ind1::String, ind2::String; alleles::Dict)
Calculates the moments based estimator of pairwise relatedness by Ritland (1996).
- Bases allele frequencies on entire population
- Inbreeding can only be assumed not to exist.
See equation 10 in: https://www.nature.com/articles/hdy201752 for variant of estimator used
Ritland original citation: https://www.cambridge.org/core/journals/genetics-research/article/estimators-for-pairwise-relatedness-and-individual-inbreeding-coefficients/9AE218BF6BF09CCCE18121AA63561CF7
"""
function lr_relatedness(data::PopData, ind1::String, ind2::String; alleles::T) where T <: NamedTuple
    #NEED TO CHECK TO CONFIRM EQUATIONS

    n = 0.0
    d = 0.0
    for loc in loci(data)
        #Extract the pair of interest's genotypes
        gen1 = get_genotype(data, sample = ind1, locus = loc)
        gen2 = get_genotype(data, sample = ind2, locus = loc)

        #Skip missings
        if gen1 !== missing && gen2 !== missing
            a,b = gen1
            c,d = gen2
            sym_loc = Symbol(loc)

            n1 = alleles[sym_loc][a] * ((b == c) + (b == d)) + alleles[sym_loc][b] * ((a == c) + (a == d)) - 4 * alleles[sym_loc][a] * alleles[sym_loc][b]
            n2 = alleles[sym_loc][c] * ((d == a) + (d == b)) + alleles[sym_loc][d] * ((c == a) + (c == b)) - 4 * alleles[sym_loc][c] * alleles[sym_loc][d]

            d1 = 2 * (1 + (a == b)) * (alleles[sym_loc][a] + alleles[sym_loc][b]) - 8 * alleles[sym_loc][a] * alleles[sym_loc][b]
            d2 = 2 * (1 + (c == d)) * (alleles[sym_loc][c] + alleles[sym_loc][d]) - 8 * alleles[sym_loc][c] * alleles[sym_loc][d]

            RL = (n1 / d1) + (n2 / d2)

            n += RL #JDS - CHECK THIS IS CORRECT
            d += ((alleles[sym_loc] |> length) - 1)
        end
    end
    return (n / d)
end

"""
    pairwise_relatedness(data::PopObj; method::String, inbreeding::Bool = true, verbose::Bool = true)
Calculates various pairwise relatedness measures between all pairs of individuals based on the entire sample population
allele frequency

If verbose is set to false then there is a progress bar. If set to true then there is
estimator specific feedback and statements when an individual has been compared to all other pairs

If the method is able to account for inbreeding in it's calculation then that option may be used

Available methods:
- `"qg"` : Queller & Goodknight 1989  (diploid only)
"""

function pairwise_relatedness(data::PopData; method::String = "qg", inbreeding::Bool = true, verbose::Bool = true)
    # check that dataset is entirely diploid
    all(data.meta.ploidy .== 2) == false && error("Relatedness analyses currently only support diploid samples")

    allele_frequencies = NamedTuple{Tuple(Symbol.(loci(data)))}(
                            Tuple(allele_freq.(locus.(Ref(data), loci(data))))
                        )

    #=
    if !verbose
        n = size(data.samples)[1]
        n = n*(n-1) ÷ 2
        prog = Progress(n, 1)
    end
    #Add switch to slightly change output depending on relatednes metric (e.g. convergence doesn't make sense for Moments Estimators)
    if method == "dyadml"
        output = DataFrame(ind1 = [], ind2 = [], relatedness = [], convergence = [])
    end
    =#
    sample_names = samples(data)
    sample_pairs = [tuple(sample_names[i], sample_names[j]) for i in 1:length(sample_names)-1 for j in i+1:length(sample_names)]
    relate_vec = zeros(length(sample_pairs))
    idx = 0
    if method == "qg"

        @inbounds for (sample_n, ind1) in enumerate(sample_names[1:end-1])
            #Base.Threads.@threads
            @inbounds Base.Threads.@threads for ind2 in sample_names[sample_n+1:end]
                idx += 1
                #=
                if method == "dyadml"
                    dyad_out = dyadicML_relatedness(data, ind1, data.samples.name[ind2], alleles = allele_frequencies, inbreeding = inbreeding, verbose = verbose)
                    append!(output,DataFrame(ind1 = ind1, ind2 = data.samples.name[ind2], relatedness = dyad_out[1], convergence = dyad_out[3]))
                end
                =#
                relate_vec[idx] += qg_relatedness(data, ind1, ind2, alleles = allele_frequencies)
                #=
                if !verbose
                    next!(prog)
                end
                =#
            end
        end

        #=
        if verbose
            println("All pairs with ", ind1, " finished")
        end
        =#
    end
    method_colname = Symbol("relatedness_" * method)
    return DataFrame(:sample_1 => getindex.(sample_pairs, 1), :sample_2 => getindex.(sample_pairs, 2), method_colname => relate_vec)
end

#Multithreaded version - requires "Combinatorics" Package
function pairwise_relatedness(data::PopObj; method::String, inbreeding::Bool = true, verbose::Bool = true)
    # check that dataset is entirely diploid
    all(data.samples.ploidy .== 2) == false && error("Relatedness analyses currently only support diploid samples")

    allele_frequencies = Dict()
    for locus in names(data.loci)
        allele_frequencies[String(locus)] = allele_freq(data.loci[:, locus])
    end

    n = size(data.samples)[1]
    n = n*(n-1) ÷ 2

    if !verbose
        prog = Progress(n, 1)
    end

    #Add switch to slightly change output depending on relatednes metric (e.g. convergence doesn't make sense for Moments Estimators)
    if method == "dyadml"
        output = DataFrame([String, String, Float64, Float64, Float64, Symbol], Symbol.(["Ind1", "Ind2", "I", "thread", "relatedness", "convergence"]), n)
    end
    if method == "qg"
        #output = DataFrame(ind1 = [], ind2 = [], relatedness = [])
        output = DataFrame([String, String, Float64, Float64, Float64], Symbol.(["Ind1", "Ind2", "I", "thread", "relatedness"]), n)
    end

    pairs = combinations(data.samples.name, 2) |> collect #This line requires Combinatorics pacakge

    Threads.@threads for i in 1:n
        ind1 = pairs[i][1]
        ind2 = pairs[i][2]

        output[i, :I] = i
        output[i, :Ind1] = ind1
        output[i, :Ind2] = ind2
        output[i, :thread] = Threads.threadid()
        if method == "qg"
            output[i, :relatedness] = qg_relatedness(data, ind1, ind2, alleles = allele_frequencies)
        end
        if method == "dyadml"
            dyad_out = dyadicML_relatedness(data, ind1, ind2, alleles = allele_frequencies, inbreeding = inbreeding, verbose = verbose)
            output[i, :relatedness] = dyad_out[1]
            output[i, :convergence] = dyad_out[3]
        end

        next!(prog)
    end
    return output
end

const relatedness = pairwise_relatedness
const kinship = pairwise_relatedness

# - `"dyadml"` : Milligan 2002 Dyadic Maximum Likelihood relatedness estimator

#=
cat_rel_noInbreeding = pairwise_relatedness(nancycats(), method = "dyadml", inbreeding = false, verbose = false)
cat_rel_Inbreeding = pairwise_relatedness(nancycats(), method = "dyadml", inbreeding = true, verbose = false)

cat_rel_qg = pairwise_relatedness(nancycats(), method = "qg", verbose = false)
=#

#=
Testing area
cat_rel_noInbreeding = pairwise_relatedness(nancycats(), method = "dyadml", inbreeding = false, verbose = false)

cat_rel_Inbreeding = pairwise_relatedness(nancycats(), method = "dyadml", inbreeding = true, verbose = false)

count(cat_rel_noInbreeding[i,4] == :Optimal for i in 1:27966)
count(cat_rel_Inbreeding[i,4] == :Optimal for i in 1:27966)

"N100"
"N106"
=#


text:
# spring-mass type model.
# to be created


text:
module Hello

function f(x)
    return x^2 -2x + x
end

end

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
module LabRig

using CImGui, DataStructures
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using DataStructures
import ThreadPools.@tspawnat
const ZEISS_TID = 3
const PIPELINE_THREADS = 4:8

include("gamepad.jl")
include("manipulator.jl")
include("pressure.jl")
include("zeiss.jl")
include("daq.jl")
#include("camera.jl")
include("gui.jl")

export Gamepad, Manipulator, Pressure, GUI, CImGui, CSyntax, CStatic, Zeiss, DAQ # Camera

function init_subsystems()
    Pressure.initialize()
    @tspawnat ZEISS_TID Zeiss.server_init()
end

end # module


text:
const pyjlsettype = pynew()

pyjlset_add(x::AbstractSet, v::Py) = (push!(x, pyconvertarg(eltype(x), v, "value")); Py(nothing))

function pyjlset_discard(x::AbstractSet, v_::Py)
    v = @pyconvert(eltype(x), v_, return Py(nothing))
    delete!(x, v)
    Py(nothing)
end

pyjlset_clear(x::AbstractSet) = (empty!(x); Py(nothing))

function pyjlset_pop(x::AbstractSet)
    if isempty(x)
        errset(pybuiltins.KeyError, "pop from an empty set")
        pynew()
    else
        Py(pop!(x))
    end
end

function pyjlset_remove(x::AbstractSet, v_::Py)
    v = @pyconvert eltype(x) v_ begin
        errset(pybuiltins.KeyError, v_)
        return pynew()
    end
    if v in x
        delete!(x, v)
        return Py(nothing)
    else
        errset(pybuiltins.KeyError, v_)
        return pynew()
    end
end

function pyjlset_update(x::AbstractSet, vs_::Py)
    for v_ in vs_
        v = pyconvert(eltype(x), v_)
        push!(x, v)
        pydel!(v_)
    end
    Py(nothing)
end

function pyjlset_difference_update(x::AbstractSet, vs_::Py)
    for v_ in vs_
        v = @pyconvert(eltype(x), v_, continue)
        delete!(x, v)
        pydel!(v_)
    end
    Py(nothing)
end

function pyjlset_intersection_update(x::AbstractSet, vs_::Py)
    vs = Set{eltype(x)}()
    for v_ in vs_
        v = @pyconvert(eltype(x), v_, continue)
        push!(vs, v)
        pydel!(v_)
    end
    intersect!(x, vs)
    Py(nothing)
end

function pyjlset_symmetric_difference_update(x::AbstractSet, vs_::Py)
    vs = Set{eltype(x)}()
    for v_ in vs_
        v = pyconvert(eltype(x), v_)
        push!(vs, v)
        pydel!(v_)
    end
    symdiff!(x, vs)
    Py(nothing)
end

function init_jlwrap_set()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class SetValue(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def add(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_add)), value)
        def discard(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_discard)), value)
        def clear(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_clear)))
        def copy(self):
            return self._jl_callmethod($(pyjl_methodnum(Py ∘ copy)))
        def pop(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_pop)))
        def remove(self, value):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_remove)), value)
        def difference(self, other):
            return set(self).difference(other)
        def intersection(self, other):
            return set(self).intersection(other)
        def symmetric_difference(self, other):
            return set(self).symmetric_difference(other)
        def union(self, other):
            return set(self).union(other)
        def isdisjoint(self, other):
            return set(self).isdisjoint(other)
        def issubset(self, other):
            return set(self).issubset(other)
        def issuperset(self, other):
            return set(self).issuperset(other)
        def difference_update(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_difference_update)), other)
        def intersection_update(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_intersection_update)), other)
        def symmetric_difference_update(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_symmetric_difference_update)), other)
        def update(self, other):
            return self._jl_callmethod($(pyjl_methodnum(pyjlset_update)), other)
    import collections.abc
    collections.abc.MutableSet.register(SetValue)
    del collections
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlsettype, jl.SetValue)
end

pyjl(v::AbstractSet) = pyjl(pyjlsettype, v)