text:
using Stheno, Plots, Random
using Stheno: @model

###########################  Define our model  ###########################

# Define a distribution over f₁, f₂, and f₃, where f₃(x) = f₁(x) + f₂(x).
@model function model()
    f₁ = GP(ConstantMean(randn()), EQ(), )
    f₂ = GP(EQ(), )
    f₃ = f₁ + f₂
    return f₁, f₂, f₃
end

# Randomly sample `N₁` locations at which to measure `f` using `y1`, and `N2` locations
# at which to measure `f` using `y2`.
rng, N₁, N₃ = MersenneTwister(123546), 10, 11;
X₁, X₃ = sort(rand(rng, N₁) * 10), sort(rand(rng, N₃) * 10);
f₁, f₂, f₃ = model();

# Generate some toy observations of `f₁` and `f₃`.
ŷ₁, ŷ₃ = rand(rng, [f₁(X₁), f₃(X₃)]);

# Compute the posterior processes.
(f₁′, f₂′, f₃′) = (f₁, f₂, f₃) | (f₁(X₁)←ŷ₁, f₃(X₃)←ŷ₃);

# Define some plotting stuff.
Np, S = 500, 25;
Xp = range(-2.5, stop=12.5, length=Np);

# Sample jointly from the posterior over each process.
f₁′Xp, f₂′Xp, f₃′Xp = rand(rng, [f₁′(Xp), f₂′(Xp), f₃′(Xp)], S);

# Compute posterior marginals.
μf₁′, σf₁′ = marginals(f₁′(Xp));
μf₂′, σf₂′ = marginals(f₂′(Xp));
μf₃′, σf₃′ = marginals(f₃′(Xp));

###########################  Plot results  ###########################

plotly();
posterior_plot = plot();

# Plot posterior marginal variances
plot!(posterior_plot, Xp, [μf₁′ μf₁′];
    linewidth=0.0,
    fillrange=[μf₁′ .- 3 .* σf₁′, μf₁′ .+ 3 * σf₁′],
    fillalpha=0.3,
    fillcolor=:red,
    label="");
plot!(posterior_plot, Xp, [μf₂′ μf₂′];
    linewidth=0.0,
    fillrange=[μf₂′ .- 3 .* σf₂′, μf₂′ .+ 3 * σf₂′],
    fillalpha=0.3,
    fillcolor=:green,
    label="");
plot!(posterior_plot, Xp, [μf₃′ μf₃′];
    linewidth=0.0,
    fillrange=[μf₃′ .- 3 .* σf₃′, μf₃′ .+ 3 * σf₃′],
    fillalpha=0.3,
    fillcolor=:blue,
    label="");

# Plot joint posterior samples
plot!(posterior_plot, Xp, f₁′Xp,
    linecolor=:red,
    linealpha=0.2,
    label="");
plot!(posterior_plot, Xp, f₂′Xp,
    linecolor=:green,
    linealpha=0.2,
    label="");
plot!(posterior_plot, Xp, f₃′Xp,
    linecolor=:blue,
    linealpha=0.2,
    label="");

# Plot posterior means
plot!(posterior_plot, Xp, μf₁′;
    linecolor=:red,
    linewidth=2.0,
    label="f1");
plot!(posterior_plot, Xp, μf₂′;
    linecolor=:green,
    linewidth=2.0,
    label="f2");
plot!(posterior_plot, Xp, μf₃′;
    linecolor=:blue,
    linewidth=2.0,
    label="f3");

# Plot observations
scatter!(posterior_plot, X₁, ŷ₁;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");
scatter!(posterior_plot, X₃, ŷ₃;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="");

display(posterior_plot);


text:
# # Stuff I Like About The Julia Programming Language

# `git_club 2020-06-19` - Hannes

# [View this notebook online](https://nbviewer.jupyter.org/github/Hasnep/stuff-i-like-about-julia/blob/master/stuff-i-like-about-julia.ipynb)

# > Julia's unofficial tagline is "Looks like Python, feels like Lisp, runs like Fortran."
# > I've been learning Julia for slightly less than a year now, and I'd like to share some of the things I've enjoyed about it.
# > I'll give an overview of the language's main features with code examples and discuss whether it really is the future scientific computing.

# Technically:
# > Julia is a high level, JIT compiled, dynamic language designed with multiple dispatch, automatic differentiation and metaprogramming.

# In practice:
# > Julia finds a balance between being fast and easy to use with lots of features you'll miss when you use another language.

# ## History

# - Julia was started in 2009 by Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and Alan Edelman
# - Announced in 2012
# - Version 1.0 was released in 2018
# - Now on version 1.4, heading for 1.5 soon

# ## Julia is easy to use

# Let's write a function from scratch to sum an array of numbers, first in Python:
# ```python
# def my_sum(array):
#     """
#     Sum a list.
#     """
#     total = 0
#     for x in array:
#         total += x
#     return total
# ```

# Then in Julia:

"""
Sum an array.
"""
function my_sum(array)
    total = 0
    for x in array
        total += x
    end
    return total
end

# This syntax will be familliar if you've used Python.

my_sum(1:10)

# Benchmarks for this function in different languages:
# - C: ~10ms
# - Python: ~500ms
# - Julia: ~10ms
# > Source: [An Introduction to Julia (Beginner Level) | SciPy 2018 Tutorial | Jane Herriman, Sacha Verweij](https://www.youtube.com/watch?v=b5xvVyzUnXI)

# ## It's fast

# Julia was designed to solve the "two language problem", where researchers use a slower, high-level language for research and then have to use a slower, low-level language once they hit a bottleneck or for production.
# Julia is both languages at the same time!

# Some benchmarks (from the Julia website) of different languages relative to C.
# ![Benchmarks of different languages relative to C](https://julialang.org/assets/benchmarks/benchmarks.svg)
# > Source: https://julialang.org/benchmarks/

# These benchmarks try to compare implementations of the same algorithm.
# Julia is comparable to compiled languages like Rust, Go, Fortran, etc. and is sometimes faster than C.
# Python is sometimes 100x slower than C.
# R is a bit slower than that.

# How is it so fast?

# Generally languages fall into two categories:

# - Compiled languages like C or Rust compile all the code before you run
# - Interpreted languages like Python or R don't compile

# Julia uses a _Just-In-Time_ (JIT) compiler which compiles a function the first time it is called.

# In my opinion, speed doesn't matter all that much

# - JIT compilation takes a short while
# - My time is more valuable than the computer's time

# Now, back to some code!

# ## Dynamically typed

# Let's write a function to calculate the $n$th Fibonacci number.
# We want the input `n` to only allow integers and we want to specify the output should only be integers.

"""
Calculate the nth Fibonacci number.
"""
function fib(n::Integer)::Integer
    if n < 2
        return n
    else
        return fib(n - 1) + fib(n - 2)
    end
end

# If we try to run the function with a non-integer input

# ```julia
# julia> fib(0.5)
# ```

# Julia will throw an error:

# ```
# ERROR: MethodError: no method matching fib(::Float64)
# Closest candidates are:
#   fib(::Integer)
# ```

fib(20)

# There is also a one-line function notation that is useful for short functions

short_fib(n::Integer)::Integer = n < 2 ? n : short_fib(n - 1) + short_fib(n - 2)

short_fib(20)

# ## It's free!

# Julia is MIT licenced.

# ## Broadcasting

# In Python, most functions accept one element.
# Running this code:
# ```python
# >>> import math
# >>> math.sin([1, 2, 3])
# ```
# will give this error:
# ```
# TypeError: must be real number, not list
# ```
# Pythonistas would probably use a list comprehension:
# ```python
# >>> [math.sin(x) for x in [1, 2, 3]]
# ```
# or a map:
# ```python
# >>> map(math.sin, [1, 2, 3])
# ```
# ```
# [0.8414709848078965, 0.9092974268256817, 0.1411200080598672]
# ```

# In R, most functions are vectorised
# ```R
# > sin(c(1, 2, 3))
# ```
# ```
# [1] 0.8414710 0.9092974 0.1411200
# ```

# In Julia, functions are not vectorised, for example, this line:
# ```julia
# julia> sin([1, 2, 3])
# ```
# gives this error:
# ```
# ERROR: MethodError: no method matching sin(::Array{Int64,1})
# ```

# Using the broadcast operator, the funciton is applied elementwise!

sin.([1, 2, 3])

# Broadcasting even works for user functions

fib.(1:10)

# And for operators:

## What are the first 10 square numbers?
(1:10).^2

# This is powerful, but sometimes tricky syntax.
# For example, adding a dot makes the length function broadcast over the array:

length(split("How many words are in this sentence?"))
# +
length.(split("How many characters are each of these words?"))

# ## Useful Unicode characters

# Julia lets you type lots of symbols using LaTeX-y abbreviations, e.g. type `\alpha` and press tab for α or `\sqrt` and tab for √.

α = 37

# The square root symbol is an abbreviation for the `sqrt()` function:

√α

# Some constants like π for pi and ℯ for Euler's constant are predefined:

## Use the approximate equals sign ≈ because this calculation is not exact
ℯ^(im * π) ≈ -1

# You can define your own operators using either built in functions:

const ⊂ = issubset
[2, 5] ⊂ [1, 2, 3, 4, 5]

# or user functions:

const ∑ = my_sum
∑(1:10)

# Because Julia supports almost any symbol you can type as a variable name, that includes emojis!

🔥 = 10
🐶 = 20
🌭 = 30
🔥 + 🐶 == 🌭

# ## Julia is (mostly) written in Julia

# If you can read Julia code, you can also read Julia's source code to understand what it does.
# I looked at the most recent PR as an example:

tensor(A::AbstractArray, B::AbstractArray) = [a * b for a in A, b in B]
const ⊗ = tensor
#-
[1, 2] ⊗ [3, 4, 5]

# Python's numpy is fast because it's mostly written in C/C++ (51.4%), but if you want to do something that numpy can't do, you need to either use C++ or write slower Python code.

# This bridges the gap between a Julia user and a Julia developer.
# For every user of Julia, there's another possible contributer.

# ## Multiple dispatch

# As an example, let's build a small Julia package based on [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl/).

"""
A number with some error.
"""
struct Uncertain <: Real
    val::Real
    err::Real
end

# Define the standard gravity on earth.

g = Uncertain(9.8, 0.1)

# Wouldn't it be nicer to show an uncertain number with the plus/minus symbol?
# Let's write a show function that dispatches on the Uncertain type:

Base.show(io::IO, x::Uncertain) = print(io, "$(x.val) ∓ $(x.err)")

g

# Let's define the plus/minus operator to create `Uncertain` numbers:

∓(a::Real, b::Real) = Uncertain(a, b)
my_height = 190 ∓ 1

# How do you add two uncertain measurements?

my_brothers_height = 175 ∓ 2

# ```julia
# julia> my_height + my_brothers_height
# ```
# ```
# ERROR: + not defined for Uncertain
# ```

# $$
# Q = a + b \\
# {\delta Q} = \sqrt{(\delta a)^2 + (\delta b)^2}
# $$

Base.:+(a::Uncertain, b::Uncertain) = (a.val + b.val) ∓ sqrt(a.err^2 + b.err^2)
my_height + my_brothers_height

# Similar for subtraction.

Base.:-(a::Uncertain, b::Uncertain) = (a.val - b.val) ∓ sqrt(a.err^2 + b.err^2)
my_height - my_brothers_height

# Slightly more complicated for multiplcation and division.

function Base.:*(a::Uncertain, b::Uncertain)
    total_value = a.val * b.val
    total_error = total_value * sqrt((a.err / a.val)^2 + (b.err / b.val)^2)
    return Uncertain(total_value, total_error)
end

function Base.:/(a::Uncertain, b::Uncertain)
    total_value = a.val / b.val
    total_error = total_value * sqrt((a.err / a.val)^2 + (b.err / b.val)^2)
    return Uncertain(total_value, total_error)
end
#-
my_height * my_brothers_height
#-
my_height / my_brothers_height

# Finally powers, again the exact formula is not important.

Base.:^(a::Uncertain, b::Real) = (a.val^b) ∓ (abs(b) * a.val^(b - 1) * a.err)
my_height^2.0

# Finally, we need to tell Julia what to do if we give it an `Uncertain` number and some other number in the same operation.

# ```julia
# julia> 1 + g
# ```
# ```
# ERROR: promotion of types Int64 and Uncertain failed to change any arguments
# ```

# We want to convert both numbers to our `Uncertain` type:

Base.promote_rule(::Type{Uncertain}, ::Type{T}) where T <: Real = Uncertain

# Coverting a real number to our Uncertain type just means we add an error of 0:

Base.convert(::Type{Uncertain}, x::Real) = Uncertain(x, 0)

# When Julia wants to convert an Uncertain number it doesn't need to do anything :

Base.convert(::Type{Uncertain}, x::Uncertain) = x
#-
1 + g

# Now we can solve a simple problem, how much time would it take if I dropped a ball from my height and my brother caught it at his height?

# Solve for t:
# $$
# t = \frac{\sqrt{2 a s + u^2} - u}{a} = \frac{\sqrt{2 g (h_1 - h_2)}}{g}
# $$

t = ((2 * g * (my_height - my_brothers_height))^0.5) / g

# Let's use the actual Measurements.jl package

using Measurements: Measurement, ±

# Let's solve the same problem as before and make sure we get the same result:

g = 9.8 ± 0.1
my_height = 190 ± 1
my_brothers_height = 175 ± 2

t = (2 * g * (my_height - my_brothers_height))^0.5 / g

# This next part is inspired by a JuliaCon talk called [_The Unreasonable Effectiveness of Multiple Dispatch_](https://www.youtube.com/watch?v=kc9HwsxE1OY) by one of the Julia co-founders, Stefan Karpinski.

# All the "methods" are outside the class definition, that can include being in a completely different package.

# Let's solve the same problem using differential equations!

# Calculate the position of the ball between time $t = 0$ and $t = 3$
# $$
# v = \frac{ds}{dt}
# $$
# $$
# f(t) = -g t = - (9.8 \pm 0.1) t
# $$
# With initial conditions $s_0 = 190 \pm 1$

using DifferentialEquations
#-
f(s, p, t) = -g * t
s₀ = my_height
time_span = (0.0, 3.0)
problem = ODEProblem(f, s₀, time_span)
#-
solution = solve(problem, Tsit5(),  saveat=0.1)
#-
using Plots
#-
plot(
    solution.t,
    solution.u,
    title="Solution to the ODE",
    xaxis="Time (t) in seconds",
    yaxis="Displacement s(t) in metres"
)

# ## Interoperability

# We can import any python module using the PyCall package:

using PyCall
#-
math = pyimport("math")
#-
math.sqrt(100)

# We can mix Python functions and variables with Julia code:

math.sin.([π,  math.pi, 2π,    2 * math.pi])

# We can define Python functions that call Julia funcitons

py"""
def pyfib(n, fun):
    if n < 2:
        return n
    else:
        return fun(n - 1, pyfib) + fun(n - 2, pyfib)
"""
#-
function jlfib(n, fun)
    if n < 2
        return n
    else
        return fun(n - 1, jlfib) + fun(n - 2, jlfib)
    end
end
#-
jlfib(20, py"pyfib")

# If you have some code written in another language (like Python), there can be a smooth transition to using Julia.

# # Summary

# - Julia looks a bit like Python
# - It's fast (with some caveats)
# - It's free!
# - It's written in Julia
# - It has powerful features like:
#     - Broadcasting syntax
#     - Unicode variables
#     - Multiple dispatch
#     - Good interoperability

# Other stuff I haven't mentioned:
# - Automatic differentiation (Swift also has this now)
# - Metaprogramming
# - Data science stuff (The first two letters of Jupyter Notebooks are named after Julia!)

ext: jl

max_stars_repo_name: Hasnep/stuff-i-like-about-julia

max_issues_repo_name: Hasnep/stuff-i-like-about-julia

max_forks_repo_name: Hasnep/stuff-i-like-about-julia


text:
"Special Additive Runge Kutta integrator."
immutable IntegratorSARK{T} <: Integrator{T}
    equation::Equation
    tableau::TableauSARK

end

"Integrate DAE with Special Additive Runge Kutta integrator."
function integrate!(int::IntegratorSARK, s::SolutionDAE)
    # TODO
end


text:
# All necessary variables (ie, setup vars) should be defined here, to be
# shared between any disparate processes
# push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
# push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
# push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

module consonance

export pars, nframes, nframes_per_proc, nchan, vels, lams, vmax, moviedir
export vmax_d, vmax_g, vmax_v

using DiskJockey.constants
using DiskJockey.model

# For odyssey use
# homedir = "/n/home07/iczekala/DiskJockey/"
# scratchdir = "/scratch/"
# outdir = "/n/home07/iczekala/DiskJockey/output/movie/" #Where the images.out and images.png are stored

# For local use
moviedir = "/home/ian/Grad/Research/Disks/DiskJockey/movie/"

# For cfa CF use
# homedir = "/pool/scout0/DiskJockey/"
# scratchdir = "/pool/cf/iczekala/scratch/"
# outdir = "/pool/scout0/DiskJockey/output/movie/"

# How many frames per process?
nframes_per_proc = 50

global const nchan = 5
global const vels = linspace(-1.5, 1.5, nchan) # [km/s]
# CO 2-1 rest frame
lam0 = cc/230.538e9 * 1e4 # [microns]
# convert velocities to wavelengths
lams = lam0 * (vels/c_kms + 1)

M_star = 1.75 # [M_sun] stellar mass
r_c = 45. # [AU] characteristic radius
r_in = 1.0 # [AU]
r_cav = 2.0 # [AU]
delta = 1.0 # no change
T_10 = 115. # [K] temperature at 10 AU
q = 0.63 # temperature gradient exponent
gamma = 1.0 # surface density gradient
logM_CO = 0.2 # [M_earth] disk mass of CO
ksi = 0.14 # [km/s] microturbulence
dpc = 73. # [pc] distance
incl = 45. # [degrees] inclination
PA = 0. # [degrees] position angle
vel = 0.0 # [km/s]
mu_RA = 0.0 # [arcsec] centroid location
mu_DEC = 0.0 # [arcsec]

# Given a starting parameter, vary it by steps dp to the low bound, then high
# bound, then back to where we started

# In the new edition, keep it at the top and bottom values for some period
function smooth_vary(start, low, high, dp, npause)
    ndown1 = round(Integer, (start - low)/dp) + 1
    nup = round(Integer, (high - low)/dp) + 1
    ndown2 = round(Integer, (high - start)/dp) + 1

    return [linspace(start, low, ndown1)' low * ones(npause)' linspace(low, high, nup)' high * ones(npause)' linspace(high, start, ndown2)' start * ones(npause)']
end

# Create a master parameter list
# First, adjust in radius
# radiuses = smooth_vary(r_c, 25., 65., 0.5)
# nradiuses = length(radiuses)

# First, adjust in inclination
incls = smooth_vary(incl, 0., 90., 1., 15)
nincls = length(incls)

# then, adjust in mass
masses = smooth_vary(M_star, 1.0, 2.5, 0.02, 15)
nmasses = length(masses)

# nframes = nincls + nmasses
#
# # Now create a giant array of Parameters objects
# pars = Array(Parameters, nframes)
#
# for i=1:nframes
#     if i <= nincls
#         pars[i] = Parameters(M_star, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, incls[i], PA, vel, mu_RA, mu_DEC)
#     else
#         pars[i] = Parameters(masses[i - nincls], r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC)
#     end
# end
#
# println("There are ", nframes, " frames to be generated.")
# println("Inclinations ", incls)
# println("Masses ", masses)
# println("Radiuses ", radiuses)

# Return a normalized instance that is symmetric about 0
# function scale(data)
#     s = maximum(abs(data))
#     return norm = plt.Normalize(vmin=-s, vmax=s, clip=false)
# end

# used as arg to imshow

# Poster plot material
nframes = 5
pars = [Parameters(1.0, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC),
Parameters(1.75, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC),
Parameters(2.50, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC),
Parameters(1.75, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, 0.0, PA, vel, mu_RA, mu_DEC),
Parameters(1.75, r_c, r_in, r_cav, delta, T_10, q, gamma, 10^logM_CO, ksi, dpc, 90.0, PA, vel, mu_RA, mu_DEC)
]

vmin_d = 0.0
vmax_d = 6.415994694015499e9

vmin_g = 0.0
vmax_g = 2.9207242096782996e11

vmin_v = 0.0
vmax_v = 45.0193854891226

end # module

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


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

using Markdown
using InteractiveUtils

# ╔═╡ 7c5ac541-689a-43af-8002-1f75f2658ac1
begin
	using CSV
	using DataFrames
	using Plots
	using StatsBase
	using StatsPlots
	using Pkg
	using Plots.Measures
end

# ╔═╡ a5c685cf-083c-447b-8940-4af7e8aaf069
begin
ENV["GKS_ENCODING"]="utf8"
pyplot()
end

# ╔═╡ c78ad002-b7db-11eb-3ed2-adfed9423e36
cd("/home/fcharih/Projects/PhD/peptide_screening/")

# ╔═╡ 45cdb33c-0610-4b80-87fd-b885d09a2bd1
targets_df = DataFrame(
	CSV.File("data/peptide_subset.csv")
)[!, [:compound_name, :targets]]

# ╔═╡ e7ea380f-a7b7-41a8-ae0c-9278fdf3ec19
targets = Dict()

# ╔═╡ 7b571d93-d6bc-4f91-b8da-67b462de3353
for row in eachrow(targets_df)
	targets[row["compound_name"]] = split(row["targets"], ";")
end

# ╔═╡ e136e461-d657-4f4d-9997-aa292f52ad45
targets

# ╔═╡ fdff33aa-d8ea-4a30-83ee-b0be15d6a0cd
pop!(targets, "albiglutide")

# ╔═╡ 9f280a6a-dafc-4cc9-80ea-c5d13494c25d
begin
pipr_minus = DataFrame(CSV.File("data/predictions/predictions_minus_160621_pipr.tsv", delim=" ", header=["peptide", "protein", "score"]))
pipr_plus = DataFrame(CSV.File("data/predictions/predictions_plus_160621_pipr.tsv", delim=" ", header=["peptide", "protein", "score"]))
sprint_minus = DataFrame(CSV.File("data/predictions/predictions_minus_old_sprint.tsv", delim=" ", header=["peptide", "protein","score"]))
sprint_plus = DataFrame(CSV.File("data/predictions/predictions_plus_160621_sprint.tsv", delim=" ", header=["peptide", "protein", "score"]))
dscript = DataFrame(CSV.File("data/predictions/dscript_base_predictions.tsv", delim=" ", header=["peptide", "protein", "score"]))
end

# ╔═╡ 23782d12-f751-4bc8-a386-dac9529b44a6
scores = vcat(
	select(pipr_minus, :, :score => ByRow(x -> "pipr_minus")),
	select(pipr_plus, :, :score => ByRow(x -> "pipr_plus")),
	select(sprint_minus, :, :score => ByRow(x -> "sprint_minus")),
	select(sprint_plus, :, :score => ByRow(x -> "sprint_plus")),
	select(dscript, :, :score => ByRow(x -> "dscript")),
)

# ╔═╡ df4d9a56-38cf-4dbe-bc9c-7a1d886c95f1
 filter!(x -> x.peptide != "albiglutide", scores)

# ╔═╡ a87fc545-8e72-4dd1-a860-5914832a9b28
begin
	methods = ["dscript", "pipr_minus", "pipr_plus", "sprint_minus", "sprint_plus"]
	data = []
	for (peptide, peptide_targets) in targets
		for method in methods
			# Assign ranks
			sub = filter(x -> x.peptide == peptide && x.score_function == method, scores)
			sub["rank"] = competerank(sub[!, :score], rev = true)
			for target in peptide_targets
				selection = sub[sub.protein .== target, :]
				named_tup = Tables.rowtable(selection)
				push!(data, named_tup)
			end
		end
	end
end

# ╔═╡ f8bb23c8-fe74-47fb-8bba-670cf669b0cd
println(targets)

# ╔═╡ 750d3d9d-f86e-45e5-96c5-bd60549098b7
begin
	df2 = DataFrame(peptide=[], protein=[], score=[], score_function=[], rank=[])
	for ele in data
		element = ele[1]
		push!(df2, [element.peptide, element.protein, element.score, element.score_function, element.rank])

	end
end

# ╔═╡ e8c6f943-0226-460a-b4e0-5d0e67ff76b4
df2

# ╔═╡ bf4dc451-679e-451b-9202-e19dd7eedbbb
function plot_curve(dataframe, peptide, method, targets)
	peptide_df = filter(x -> x.peptide == peptide && x.score_function == method, dataframe)
	sort!(peptide_df, :score, rev=true)
	peptide_df["rank"] = 1:nrow(peptide_df)

	targets_df = filter(x -> x.protein ∈ targets[peptide], peptide_df)

	o2a = scatter(peptide_df.rank, peptide_df.score, framestyle=:box, legend=false, xlabel="Rank", ylabel="Score")
	scatter!(targets_df.rank, targets_df.score, color=:red, marker=8, formatter=:plain)

	o2a
end

# ╔═╡ b8d118cb-2fb8-466e-91ef-d4ea47457d6c
plot_curve(scores, "albiglutide", "sprint_minus", targets)

# ╔═╡ 57b61759-6056-493c-89f8-2eaad2347c81
ranking_matrix = vcat([
	df2[df2.score_function .== method, :]["rank"]
	for method in methods
])

# ╔═╡ 956958b7-69e3-4087-a231-b28119e8fb58
begin
	labels2 = ["PIPR (-)" "PIPR (+)" "SPRINT (-)" "SPRINT (+)" "D-SCRIPT\n(pre-trained)"]
	box_plot2 = boxplot(labels2, ranking_matrix, leg = false, framestyle=:box, outliers=false, formatter=:plain, xlabel="Method", ylabel="Rank")
	dotplot!(labels2, ranking_matrix)
end

# ╔═╡ 14dd7609-8784-4322-9200-f0fc160769b4
scores_receptors = vcat(
	select(pipr_minus, :, :score => ByRow(x -> "pipr_minus")),
	select(pipr_plus, :, :score => ByRow(x -> "pipr_plus")),
	select(sprint_minus, :, :score => ByRow(x -> "sprint_minus")),
	select(sprint_plus, :, :score => ByRow(x -> "sprint_plus")),
	select(dscript, :, :score => ByRow(x -> "dscript")),
)

# ╔═╡ fbd43006-fd8f-4eff-b6b0-ca64aa5106bd
begin
	surface_proteins = Set(readlines("data/surfaceome.txt"))
	scores_surface = filter(row -> (row.protein ∈ surface_proteins) || (row.protein ∈ targets[row.peptide]), scores)
	filter!(x -> x.peptide != "albiglutide", scores_surface)
end

# ╔═╡ c8551873-daff-42d6-9647-9ddf3b99a068

# ╔═╡ 87d71a66-9709-4f83-a68f-e34e1c0d337e
sort(filter(x -> x.peptide == "carperitide" && x.score_function == "sprint_plus", scores_surface), :score, rev=true)

# ╔═╡ b71647db-b07d-4875-9477-93ac07539cd8
length(surface_proteins)

# ╔═╡ 0d608224-268c-4386-bac4-4ab41fb47a49
filter(x -> x.peptide == "bivalirudin" && x.protein == "P00734", scores)

# ╔═╡ 3b072371-c283-45e6-a817-c5715a0583a4
begin
	methods2 = ["dscript", "pipr_minus", "pipr_plus", "sprint_minus", "sprint_plus"]
	data2 = []
	for (peptide, peptide_targets) in targets
		for method in methods2
			# Assign ranks
			sub = filter(x -> x.peptide == peptide && x.score_function == method, scores_surface)
			sub["rank"] = competerank(sub[!, :score], rev = true)
			println(peptide)
			for target in peptide_targets

				selection = sub[sub.protein .== target, :]
				println(selection)
				named_tup = Tables.rowtable(selection)
				push!(data2, named_tup)
			end
		end
	end
end

# ╔═╡ 37fe6fdf-ff7b-4986-84d5-820b2fc2da3e
data2

# ╔═╡ 6a0e91da-4d03-439b-b64f-5da983d2ae0c
begin
	df3 = DataFrame(peptide=[], protein=[], score=[], score_function=[], rank=[])
	for ele in data2

		if length(ele) == 0
			println("SKIP")
			continue
		end

		element = ele[1]

		push!(df3, [element.peptide, element.protein, element.score, element.score_function, element.rank])

	end
end

# ╔═╡ 73f0d7c2-f833-450d-a49d-6b0b63a05b54
df3

# ╔═╡ 0bafaf4b-3ea3-4cb3-88df-5331c80507e2
ranking_matrix2 = vcat([
	df3[df3.score_function .== method, :]["rank"]
	for method in methods
])

# ╔═╡ 7d9821d4-f1a3-4afc-8496-667aa67cfc9c
begin
	labels3 = ["D-SCRIPT\n(pre-trained)" "PIPR (-)" "PIPR (+)" "SPRINT (-)" "SPRINT (+)" ]
	box_plot3 = boxplot(labels3, ranking_matrix2, leg = false, framestyle=:box, outliers=false, formatter=:plain, xlabel="Method", ylabel="Rank")
	dotplot!(labels3, ranking_matrix2)
end

# ╔═╡ 45a45178-55c8-4f43-a2f6-a54309c267a0
savefig(box_plot3, "figures/surface_proteins_boxplot.pdf")

# ╔═╡ 82086054-e306-4564-87c7-5db01ee6752b
df3[df3.score_function .== "sprint_plus", :]

# ╔═╡ 592bc548-2f5b-4271-a556-a14f1f18f281
begin
	num_in_top20 = Dict()
	for method in methods
		num_in_top20[method] = nrow(filter(x -> x.rank <= 20 && x.score_function == method, df3))
	end
end

# ╔═╡ 13637bc2-ba9d-4d82-ade3-b2adae750a71
num_in_top20

# ╔═╡ ad50b95b-b876-4f97-8f41-7b84af38445c
begin
bar(["PIPR (-)", "PIPR (+)", "SPRINT (-)", "SPRINT (+)", "D-SCRIPT\n(pre-trained)"], [0, 2, 5, 15, 0], legend=false, framestyle=:box, xlabel="Method", ylabel="Count")
	savefig("figures/top20.pdf")
end

# ╔═╡ ad5e3ba8-3158-494c-bd6b-6cec1487828b

# ╔═╡ 5b7673f3-cbd1-4df9-ac22-e20c358924e0
begin
plot_curve(scores_surface, "bivalirudin", "dscript", targets)
savefig("figures/bivalirudin.pdf")
end

# ╔═╡ e7b92ae5-02e8-4e19-8594-5f929c592d19
targets

# ╔═╡ 8daa52dd-d29d-4e16-83dd-b5b7c84b4e5b

# ╔═╡ 220d4861-c3c9-409d-8487-0f60d507d674
begin
	pp = "teduglutide"
	plot_curve(scores_surface, "$pp", "sprint_plus", targets)
	savefig("figures/$(pp)_sprint_plus_surface.pdf")
end

# ╔═╡ 945781f8-c827-4cd4-8a85-cf90ae91ad12
plot_curve(scores_surface, "$pp", "dscript", targets)

# ╔═╡ 409b4844-021c-4bb2-be5a-4a6d1982cf1d
plot_curve(scores_surface, "$pp", "sprint_plus", targets)

# ╔═╡ d3c32902-d420-43f9-a179-3a91d4bd82cd
begin
	pp2 = "carperitide"
	plot_curve(scores_surface, "$(pp2)", "sprint_minus", targets)
	savefig("figures/$(pp2)_sprint_minus_surface.pdf")
	plot_curve(scores_surface, "$(pp2)", "sprint_plus", targets)
	savefig("figures/$(pp2)_sprint_plus_surface.pdf")
	plot_curve(scores_surface, "$(pp2)", "pipr_minus", targets)
	savefig("figures/$(pp2)_pipr_minus_surface.pdf")
	plot_curve(scores_surface, "$(pp2)", "pipr_plus", targets)
	savefig("figures/$(pp2)_pipr_plus_surface.pdf")
	plot_curve(scores_surface, "$(pp2)", "dscript", targets)
	savefig("figures/$(pp2)_dscript_surface.pdf")
end

# ╔═╡ b6c4aace-610b-4917-9918-7ff2730aa303
filter(x -> x.peptide == "lixisenatide", scores_surface)

# ╔═╡ 3c6371af-2f54-4a31-9746-a1f70facbd50
function get_target_ranks_for_peptide(dataframe, peptide, targets)
	df = filter(x -> x.peptide == peptide, dataframe)
	sort!(df, [:score], rev = true)
	df[!, :rank] = competerank(df[!, :score], rev = true)
	filter(row -> row.protein ∈ targets[peptide], df)[!, :rank]
end

# ╔═╡ 870c1de8-596d-454b-a661-837fb17c16a2
function get_ranks_for_method(dataframe, targets)
	ranks = []
	for peptide in keys(targets)
		peptide_ranks = get_target_ranks_for_peptide(dataframe, peptide, targets)
		append!(ranks, peptide_ranks)
	end
	return ranks
end

# ╔═╡ 54d7153a-3a86-4f5f-92f7-ac9b48c9eca1
get_ranks_for_method(dscript, targets)

# ╔═╡ f6ba3f69-fee5-42b7-bbc9-3113da34bd9a
ranks = DataFrame([
	get_ranks_for_method(pipr_minus, targets),
	get_ranks_for_method(pipr_plus, targets),
	get_ranks_for_method(sprint_minus, targets),
	get_ranks_for_method(sprint_plus, targets),
	get_ranks_for_method(dscript, targets)
])

# ╔═╡ 35dc5c45-6d0e-4b73-915d-56f2c8e06795
rename!(ranks, Symbol.(["pipr_minus", "pipr_plus", "sprint_minus", "sprint_plus", "dscript"]))

# ╔═╡ 8f7bddcd-71a5-4c3b-ac27-36c40c429455
begin
	ranks_matrix = Matrix{Int64}(ranks)
	labels = ["PIPR (-)" "PIPR (+)" "SPRINT (-)" "SPRINT (+)" "D-SCRIPT\n(pre-trained)"]
	box_plot = boxplot(labels, ranks_matrix, leg = false, framestyle=:box, outliers=false, formatter=:plain, xlabel="Method", ylabel="Rank")
	dotplot!(labels, ranks_matrix)
end

# ╔═╡ bd2e27af-efdd-4512-a262-c705d925666c
savefig(box_plot, "figures/boxplot.pdf")

# ╔═╡ 0124ff9f-cc00-4dd4-bd34-c9fe7df1c2fd
function get_hist_series(ranks)
	x = range(1, 20400, step=1)
	y = []
	for i in x
		push!(y, length(filter(x -> x < i, ranks)))
	end
	return (x, y)
end

# ╔═╡ 2f6e1830-1b02-458e-9ae1-38908222917a
begin
rising_curve = plot(get_hist_series(ranks.pipr_minus), label="PIPR (-)", framestyle=:box, xlabel="Rank cutoff", ylabel="Targets below cutoff", formatter=:plain)
plot!(get_hist_series(ranks.pipr_plus), label="PIPR (+)")
plot!(get_hist_series(ranks.sprint_minus), label="SPRINT (-)")
plot!(get_hist_series(ranks.sprint_plus), label="SPRINT (+)")
plot!(get_hist_series(ranks.dscript), label="D-SCRIPT (pre-trained)")

plot!(
    [
	get_hist_series(ranks.pipr_minus),
	get_hist_series(ranks.pipr_plus),
	get_hist_series(ranks.sprint_minus),
	get_hist_series(ranks.sprint_plus),
	get_hist_series(ranks.dscript)
	],
    inset = bbox(0.3, 0.2, 0.3, 0.3, :center),
    subplot = 2,
	framestyle=:box,
	formatter=:plain,
	legend=nothing,
	xrange=(0,50),
	yrange=(0,20)
)
end

# ╔═╡ 6cb631f5-bcc1-45f6-b974-e7015b888384
savefig(rising_curve, "figures/targets_below_rank_plot.pdf")

# ╔═╡ 60b2faaa-71d2-4f77-850a-4910c72e04d3
histogram(sprint_plus.score, bins=100, yrange=(0,1000))

# ╔═╡ Cell order:
# ╠═7c5ac541-689a-43af-8002-1f75f2658ac1
# ╠═a5c685cf-083c-447b-8940-4af7e8aaf069
# ╠═c78ad002-b7db-11eb-3ed2-adfed9423e36
# ╠═45cdb33c-0610-4b80-87fd-b885d09a2bd1
# ╠═e7ea380f-a7b7-41a8-ae0c-9278fdf3ec19
# ╠═7b571d93-d6bc-4f91-b8da-67b462de3353
# ╠═e136e461-d657-4f4d-9997-aa292f52ad45
# ╠═fdff33aa-d8ea-4a30-83ee-b0be15d6a0cd
# ╠═9f280a6a-dafc-4cc9-80ea-c5d13494c25d
# ╠═23782d12-f751-4bc8-a386-dac9529b44a6
# ╠═df4d9a56-38cf-4dbe-bc9c-7a1d886c95f1
# ╠═a87fc545-8e72-4dd1-a860-5914832a9b28
# ╠═f8bb23c8-fe74-47fb-8bba-670cf669b0cd
# ╠═750d3d9d-f86e-45e5-96c5-bd60549098b7
# ╠═e8c6f943-0226-460a-b4e0-5d0e67ff76b4
# ╠═bf4dc451-679e-451b-9202-e19dd7eedbbb
# ╠═b8d118cb-2fb8-466e-91ef-d4ea47457d6c
# ╠═57b61759-6056-493c-89f8-2eaad2347c81
# ╠═956958b7-69e3-4087-a231-b28119e8fb58
# ╠═14dd7609-8784-4322-9200-f0fc160769b4
# ╠═fbd43006-fd8f-4eff-b6b0-ca64aa5106bd
# ╠═c8551873-daff-42d6-9647-9ddf3b99a068
# ╠═87d71a66-9709-4f83-a68f-e34e1c0d337e
# ╠═945781f8-c827-4cd4-8a85-cf90ae91ad12
# ╠═b71647db-b07d-4875-9477-93ac07539cd8
# ╠═0d608224-268c-4386-bac4-4ab41fb47a49
# ╠═3b072371-c283-45e6-a817-c5715a0583a4
# ╠═37fe6fdf-ff7b-4986-84d5-820b2fc2da3e
# ╠═6a0e91da-4d03-439b-b64f-5da983d2ae0c
# ╠═73f0d7c2-f833-450d-a49d-6b0b63a05b54
# ╠═0bafaf4b-3ea3-4cb3-88df-5331c80507e2
# ╠═7d9821d4-f1a3-4afc-8496-667aa67cfc9c
# ╠═45a45178-55c8-4f43-a2f6-a54309c267a0
# ╠═82086054-e306-4564-87c7-5db01ee6752b
# ╠═592bc548-2f5b-4271-a556-a14f1f18f281
# ╠═13637bc2-ba9d-4d82-ade3-b2adae750a71
# ╠═ad50b95b-b876-4f97-8f41-7b84af38445c
# ╠═ad5e3ba8-3158-494c-bd6b-6cec1487828b
# ╠═5b7673f3-cbd1-4df9-ac22-e20c358924e0
# ╠═e7b92ae5-02e8-4e19-8594-5f929c592d19
# ╠═8daa52dd-d29d-4e16-83dd-b5b7c84b4e5b
# ╠═220d4861-c3c9-409d-8487-0f60d507d674
# ╠═409b4844-021c-4bb2-be5a-4a6d1982cf1d
# ╠═d3c32902-d420-43f9-a179-3a91d4bd82cd
# ╠═b6c4aace-610b-4917-9918-7ff2730aa303
# ╠═3c6371af-2f54-4a31-9746-a1f70facbd50
# ╠═870c1de8-596d-454b-a661-837fb17c16a2
# ╠═54d7153a-3a86-4f5f-92f7-ac9b48c9eca1
# ╠═f6ba3f69-fee5-42b7-bbc9-3113da34bd9a
# ╠═35dc5c45-6d0e-4b73-915d-56f2c8e06795
# ╠═8f7bddcd-71a5-4c3b-ac27-36c40c429455
# ╠═bd2e27af-efdd-4512-a262-c705d925666c
# ╠═0124ff9f-cc00-4dd4-bd34-c9fe7df1c2fd
# ╠═2f6e1830-1b02-458e-9ae1-38908222917a
# ╠═6cb631f5-bcc1-45f6-b974-e7015b888384
# ╠═60b2faaa-71d2-4f77-850a-4910c72e04d3

max_stars_repo_name: GreenCUBIC/PeptideScreening

max_issues_repo_name: GreenCUBIC/PeptideScreening

max_forks_repo_name: GreenCUBIC/PeptideScreening


text:
using Test
using DataValues

@testset "DataValueArray: TypeDefs" begin

    x = DataValueArray([1, 2, 3],
    [false, false, true])

    y = DataValueArray([
        1 2;
        3 4;
    ],
    [
        false false;
        true false;
    ],
)

    @test isa(x, DataValueArray{Int,1})
    @test isa(x, DataValueVector{Int})

    @test isa(y, DataValueArray{Int,2})
    @test isa(y, DataValueMatrix{Int})

end


text:
@testset "airspace.jl" begin
    p1 = Point(0.0, 0.0)
    p2 = Point(0.0, 10.0)
    p3 = Point(10.0, 10.0)
    p4 = Point(10.0, 0.0)

    ap = Airspace("test", [p1, p2, p3, p4, p1])

    mp = midpoint(p3, p4)
    @test ap.name == "test"
    @test ap.polygon[1] == p1
    @test ap.polygon[end] == p1
    @test ap.bounding_box[3].ϕ ≈ mp.ϕ atol = 0.0001

    @test isinside(Point(5.0, 5.0), ap.bounding_box)
    @test isinside(Point(15.0, 5.0), ap.bounding_box) == false
    @test isinside(Point(5.0, 15.0), ap.bounding_box) == false

    @test intersection_point(Point(15.0, 5.0), Point(-90.0, 0.0),
    ap.polygon).ϕ ≈ 10.03742304591071 atol = 0.0001
    @test intersection_point(Point(5.0, 5.0), Point(-90.0, 0.0),
    ap.polygon).ϕ ≈ 0.0 atol = 0.0001

    @test intersection_point(Point(15.0, 5.0), Point(-90.0, 0.0), ap).ϕ ≈
    10.03742304591071 atol = 0.0001
    @test intersection_point(Point(5.0, 5.0), Point(-90.0, 0.0), ap).ϕ ≈
    0.0 atol = 0.0001
    @test isnan(intersection_point(Point(5.0, 15.0), Point(-90.0, 0.0), ap).ϕ)

    ap_triangle = Airspace("triangle", [p1, p2, p3, p1])
    @test isinside(Point(9.0, 1.0), ap_triangle.bounding_box) == true
    @test isinside(Point(9.0, 1.0), ap_triangle.polygon) == false

    @test isnan(intersection_point(Point(0.0, -1.0), Point(11.0, 10.0), ap_triangle).ϕ)
    @test intersection_point(Point(11.0, -1.0), Point(-1.0, 11.0), ap_triangle).ϕ ≈ 5.0703 atol = 0.001
    @test intersection_point(Point(-1.0, 11.0), Point(11.0, -1.0), ap_triangle).λ ≈ 10.000   atol = 0.001
end


text:
module PopulationDynamics
using Reexport
using utils
@reexport using EPhys

if VERSION >= v"0.7.0"
    using Random
    using Statistics
    using LinearAlgebra
end

using MultivariateStats, StatsBase
using RecursiveArrayTools, Munkres
#TensorDecompositions

export fit!, weights
export track, track_partial, bootstrap, confint
export similarity, stability
export factors

abstract type PopulationModel end

include("utils.jl")
include("pcaica.jl")
include("types.jl")
#include("tensor_decompositions.jl")

include("ephysIO.jl")
#include("coactivation.jl")
#include("coactivation_surrogates.jl")
#include("crosscor.jl")

end


text:
function Combinations(arg0::jint, arg1::jint)
    return Combinations((jint, jint), arg0, arg1)
end

function comparator(obj::Combinations)
    return jcall(obj, "comparator", Comparator, ())
end

function for_each(obj::Iterable, arg0::Consumer)
    return jcall(obj, "forEach", void, (Consumer,), arg0)
end

function get_k(obj::Combinations)
    return jcall(obj, "getK", jint, ())
end

function get_n(obj::Combinations)
    return jcall(obj, "getN", jint, ())
end

function iterator(obj::Combinations)
    return jcall(obj, "iterator", JIterator, ())
end

function spliterator(obj::Iterable)
    return jcall(obj, "spliterator", Spliterator, ())
end


text:
function ref_add_goc!(pm::AbstractPowerModel)
    if _IM.ismultinetwork(pm.data)
        nws_data = pm.data["nw"]
    else
        nws_data = Dict("0" => pm.data)
    end

    for (n, nw_data) in nws_data
        nw_id = parse(Int, n)
        ref = pm.ref[:nw][nw_id]

        ref[:branch_sm_active] = Dict(x for x in ref[:branch] if haskey(x.second, "rate_a"))

        ref[:shunt_const] = Dict(x for x in ref[:shunt] if (!haskey(x.second, "dispatchable") || !x.second["dispatchable"]))
        ref[:shunt_var] = Dict(x for x in ref[:shunt] if (haskey(x.second, "dispatchable") && x.second["dispatchable"]))

        bus_shunts_const = Dict((i, []) for (i,bus) in ref[:bus])
        for (i,shunt) in ref[:shunt_const]
            push!(bus_shunts_const[shunt["shunt_bus"]], i)
        end
        ref[:bus_shunts_const] = bus_shunts_const

        bus_shunts_var = Dict((i, []) for (i,bus) in ref[:bus])
        for (i,shunt) in ref[:shunt_var]
            push!(bus_shunts_var[shunt["shunt_bus"]], i)
        end
        ref[:bus_shunts_var] = bus_shunts_var
    end
end


text:
include("costs.jl")
include("approximators/approximators.jl")


text:
function mm_spec_jaeschkethesis2_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(5, 0, 0, 3, 7, 9, 0)
    # x = [ c_A c_B c_C F_A F_B]
    # p = [ k_1 k_2 (-\Delta H_1) (-\Delta H_2) c_{A,in} c_{B,in} V F_{max}, q_{max} ]
    return dimspec
end

function mm_spec_jaeschkethesis2_defn(x, pr, ps)

    f = ( (-1) * (x[4] + x[5])^2 * (x[3]^2) ) / (x[4] * pr[5])

    ce = Vector{Num}(undef, 3)
    ce[1] = (x[4] * pr[5]) - ((x[4] + x[5]) * x[1]) - (pr[1] * x[1] * x[2] * pr[7])
    ce[2] = (x[5] * pr[6]) - ((x[4] + x[5]) * x[2]) - (pr[1] * x[1] * x[2] * pr[7]) - (2 * pr[2] * x[2]^2 * pr[7])
    ce[3] = -((x[4] + x[5]) * x[3]) + (pr[1] * x[1] * x[2] * pr[7])

    ci = Vector{Num}(undef, 7)
    ci[1] = x[1] + 1e-15
    ci[2] = x[2] + 1e-15
    ci[3] = x[3] + 1e-15
    ci[4] = x[4] + 1e-15
    ci[5] = x[5] + 1e-15
    ci[6] = pr[8] - (x[4] + x[5])
    ci[7] = pr[9] - (pr[1] * x[1] * x[2] * pr[7] * pr[3]) - (2 * pr[2] * x[2]^2 * pr[7] * pr[4])

    F = Matrix{Num}(undef, 0, 0)
    # F[1,1] = x[1]
    # F[2,1] = x[2]

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_jaeschkethesis2_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_jaeschkethesis2_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_jaeschkethesis2_knownsols()

    return [ ]
end

function mm_spec_jaeschkethesis2_parameterisations(t)

    defn1 = MPCCParameterisationDefn(
            Vector{Num}([ t, 0.014, 7e4, 5e4, 2, 1.5, 500, 22, 1e6 ]),
            (0.3, 1.5),
            "Standard"
        )

    return Vector{MPCCParameterisationDefn}( [ defn1 ] )
end


text:
module GEE

using LinearAlgebra, Distributions, StatsBase, GLM, Tables
using StatsBase: CoefTable, StatisticalModel, RegressionModel
using LinearAlgebra: BlasReal
using GLM: Link, LinPred

import StatsBase: fit, coef
import GLM: linkfun

export GeneralizedEstimatingEquationsModel,
       GeneralizedEstimatingEquation,
       CorrelationStructure,
       Independence,
       Exchangeable,
       Unstructered,
       fit,
       vcov,
       coef,
       stderror,
       Normal,
       Binomial,
       Poisson,
       LogitLink,
       IdentityLink,
       LogLink,
       @formula

abstract type GeneralizedEstimatingEquationModel <: StatisticalModel end

include("workingstructures.jl")
include("generalizedestimatingequation.jl")
include("correlation.jl")
include("sandwich.jl")

end


text:
@doc raw"""
    AbstractKernel

Integral kernel `K(x, y)`.

Abstract base type for an integral kernel, i.e. a AbstractFloat binary function
``K(x, y)`` used in a Fredhold integral equation of the first kind:
```math
    u(x) = ∫ K(x, y) v(y) dy
```
where ``x ∈ [x_\mathrm{min}, x_\mathrm{max}]`` and
``y ∈ [y_\mathrm{min}, y_\mathrm{max}]``.  For its SVE to exist,
the kernel must be square-integrable, for its singular values to decay
exponentially, it must be smooth.

In general, the kernel is applied to a scaled spectral function ``ρ'(y)`` as:
```math
    ∫ K(x, y) ρ'(y) dy,
```
where ``ρ'(y) = w(y) ρ(y)``.
"""
abstract type AbstractKernel end

@doc raw"""
    LogisticKernel <: AbstractKernel

Fermionic/bosonic analytical continuation kernel.

In dimensionless variables ``x = 2 τ/β - 1``, ``y = β ω/Λ``,
the integral kernel is a function on ``[-1, 1] × [-1, 1]``:
```math
    K(x, y) = \frac{e^{-Λ y (x + 1) / 2}}{1 + e^{-Λ y}}
```
LogisticKernel is a fermionic analytic continuation kernel.
Nevertheless, one can model the ``τ`` dependence of
a bosonic correlation function as follows:
```math
    ∫ \frac{e^{-Λ y (x + 1) / 2}}{1 - e^{-Λ y}} ρ(y) dy = ∫ K(x, y) ρ'(y) dy,
```
with
```math
    ρ'(y) = w(y) ρ(y),
```
where the weight function is given by
```math
    w(y) = \frac{1}{\tanh(Λ y/2)}.
```
"""
struct LogisticKernel <: AbstractKernel
    Λ::Float64
    function LogisticKernel(Λ)
        Λ ≥ 0 || throw(DomainError(Λ, "Kernel cutoff Λ must be non-negative"))
        return new(Λ)
    end
end

@doc raw"""
    RegularizedBoseKernel <: AbstractKernel

Regularized bosonic analytical continuation kernel.

In dimensionless variables ``x = 2 τ/β - 1``, ``y = β ω/Λ``, the fermionic
integral kernel is a function on ``[-1, 1] × [-1, 1]``:
```math
    K(x, y) = y \frac{e^{-Λ y (x + 1) / 2}}{e^{-Λ y} - 1}
```
Care has to be taken in evaluating this expression around ``y = 0``.
"""
struct RegularizedBoseKernel <: AbstractKernel
    Λ::Float64
    function RegularizedBoseKernel(Λ)
        Λ ≥ 0 || throw(DomainError(Λ, "Kernel cutoff Λ must be non-negative"))
        return new(Λ)
    end
end

"""
    AbstractSVEHints

Discretization hints for singular value expansion of a given kernel.
"""
abstract type AbstractSVEHints end

struct SVEHintsLogistic{T} <: AbstractSVEHints
    kernel::LogisticKernel
    ε::T
end

struct SVEHintsRegularizedBose{T} <: AbstractSVEHints
    kernel::RegularizedBoseKernel
    ε::T
end

struct SVEHintsReduced{T<:AbstractSVEHints} <: AbstractSVEHints
    inner_hints::T
end

abstract type AbstractReducedKernel <: AbstractKernel end

@doc raw"""
    ReducedKernel

Restriction of centrosymmetric kernel to positive interval.

For a kernel ``K`` on ``[-1, 1] × [-1, 1]`` that is centrosymmetric, i.e.
``K(x, y) = K(-x, -y)``, it is straight-forward to show that the left/right
singular vectors can be chosen as either odd or even functions.

Consequentially, they are singular functions of a reduced kernel ``K_\mathrm{red}``
on ``[0, 1] × [0, 1]`` that is given as either:
```math
    K_\mathrm{red}(x, y) = K(x, y) \pm K(x, -y)
```
This kernel is what this type represents.  The full singular functions can
be reconstructed by (anti-)symmetrically continuing them to the negative
axis.
"""
struct ReducedKernel{K<:AbstractKernel} <: AbstractReducedKernel
    inner::K
    sign::Int
end

@doc raw"""
    LogisticKernelOdd <: AbstractReducedKernel

Fermionic analytical continuation kernel, odd.

In dimensionless variables ``x = 2τ/β - 1``, ``y = βω/Λ``, the fermionic
integral kernel is a function on ``[-1, 1] × [-1, 1]``:
```math
    K(x, y) = -\frac{\sinh(Λ x y / 2)}{\cosh(Λ y / 2)}
```
"""
struct LogisticKernelOdd <: AbstractReducedKernel
    inner::LogisticKernel
    sign::Int

    function LogisticKernelOdd(inner::LogisticKernel, sign)
        iscentrosymmetric(inner) || error("inner kernel must be centrosymmetric")
        abs(sign) == 1 || throw(DomainError(sign, "sign must be -1 or 1"))
        return new(inner, sign)
    end
end

@doc raw"""
    RegularizedBoseKernelOdd <: AbstractReducedKernel

Bosonic analytical continuation kernel, odd.

In dimensionless variables ``x = 2 τ / β - 1``, ``y = β ω / Λ``, the fermionic
integral kernel is a function on ``[-1, 1] × [-1, 1]``:
```math
    K(x, y) = -y \frac{\sinh(Λ x y / 2)}{\sinh(Λ y / 2)}
```
"""
struct RegularizedBoseKernelOdd <: AbstractReducedKernel
    inner::RegularizedBoseKernel
    sign::Int

    function RegularizedBoseKernelOdd(inner::RegularizedBoseKernel, sign)
        iscentrosymmetric(inner) || error("inner kernel must be centrosymmetric")
        abs(sign) == 1 || throw(DomainError(sign, "sign must be -1 or 1"))
        return new(inner, sign)
    end
end

@doc raw"""
    xrange(kernel)

Return a tuple ``(x_\mathrm{min}, x_\mathrm{max})`` delimiting the range
of allowed `x` values.
"""
xrange(::AbstractKernel) = (-1, 1)
xrange(kernel::AbstractReducedKernel) = (0, last(xrange(kernel.inner)))

@doc raw"""
    yrange(kernel)

Return a tuple ``(y_\mathrm{min}, y_\mathrm{max})`` delimiting the range
 of allowed `y` values.
"""
yrange(::AbstractKernel) = (-1, 1)
yrange(kernel::AbstractReducedKernel) = (0, last(yrange(kernel.inner)))

function compute(::LogisticKernel, u₊, u₋, v)
    # By introducing u_± = (1 ± x)/2 and v = Λ * y, we can write
    # the kernel in the following two ways:
    #
    #    k = exp(-u₊ * v) / (exp(-v) + 1)
    #      = exp(-u₋ * -v) / (exp(v) + 1)
    #
    # We need to use the upper equation for v ≥ 0 and the lower one for
    # v < 0 to avoid overflowing both numerator and denominator

    enum = exp(-abs(v) * (v ≥ 0 ? u₊ : u₋))
    denom = 1 + exp(-abs(v))
    return enum / denom
end

function compute(kernel::RegularizedBoseKernel, u₊, u₋, v)
    # With "reduced variables" u, v we have:
    #
    #   K = -1/lambda * exp(-u_+ * v) * v / (exp(-v) - 1)
    #     = -1/lambda * exp(-u_- * -v) * (-v) / (exp(v) - 1)
    #
    # where we again need to use the upper equation for v ≥ 0 and the
    # lower one for v < 0 to avoid overflow.
    absv = abs(v)
    enum = exp(-absv * (v ≥ 0 ? u₊ : u₋))
    T = eltype(v)

    # The expression ``v / (exp(v) - 1)`` is tricky to evaluate: firstly,
    # it has a singularity at v=0, which can be cured by treating that case
    # separately.  Secondly, the denominator loses precision around 0 since
    # exp(v) = 1 + v + ..., which can be avoided using expm1(...)
    denom = absv ≥ 1e-200 ? absv / expm1(-absv) : -one(absv)
    return -1 / T(kernel.Λ) * enum * denom
end

"""
    segments_x(kernel)

Segments for piecewise polynomials on the ``x`` axis.

List of segments on the ``x`` axis for the associated piecewise
polynomial. Should reflect the approximate position of roots of a
high-order singular function in ``x``.
"""
function segments_x end

"""
    segments_y(kernel)

Segments for piecewise polynomials on the ``y`` axis.

List of segments on the ``y`` axis for the associated piecewise
polynomial. Should reflect the approximate position of roots of a
high-order singular function in ``y``.
"""
function segments_y end

function segments_x(hints::SVEHintsLogistic)
    nzeros = max(round(Int, 15 * log10(hints.kernel.Λ)), 1)
    diffs = 1 ./ cosh.(0.143 * range(0; length=nzeros))
    zeros = cumsum(diffs)
    zeros ./= last(zeros)
    return [-reverse(zeros); 0; zeros]
end

function segments_y(hints::SVEHintsLogistic)
    nzeros = max(round(Int, 20 * log10(hints.kernel.Λ)), 2)

    # Zeros around -1 and 1 are distributed asymptotically identically
    leading_diffs = [0.01523, 0.03314, 0.04848, 0.05987, 0.06703, 0.07028, 0.07030, 0.06791,
        0.06391, 0.05896, 0.05358, 0.04814, 0.04288, 0.03795, 0.03342, 0.02932, 0.02565,
        0.02239, 0.01951, 0.01699][1:min(nzeros, 20)]

    diffs = [leading_diffs; 0.25 ./ exp.(0.141 * (20:(nzeros - 1)))]
    zeros = cumsum(diffs)
    zeros ./= pop!(zeros)
    zeros .-= 1
    return [-1; zeros; 0; -reverse(zeros); 1]
end

function segments_x(hints::SVEHintsRegularizedBose)
    # Somewhat less accurate...
    nzeros = max(round(Int, 15 * log10(hints.kernel.Λ)), 15)
    diffs = 1 ./ cosh.(0.18 * range(0; length=nzeros))
    zeros = cumsum(diffs)
    zeros ./= last(zeros)
    return [-reverse(zeros); 0; zeros]
end

function segments_y(hints::SVEHintsRegularizedBose)
    nzeros = max(round(Int, 20 * log10(hints.kernel.Λ)), 20)
    i = range(0; length=nzeros)
    diffs = @. 0.12 / exp(0.0337 * i * log(i + 1))
    zeros = cumsum(diffs)
    zeros ./= pop!(zeros)
    zeros .-= 1
    return [-1; zeros; 0; -reverse(zeros); 1]
end

"""
    matrix_from_gauss(kernel, gauss_x, gauss_y)

Compute matrix for kernel from Gauss rules.
"""
function matrix_from_gauss(kernel, gauss_x, gauss_y)
    # (1 ± x) is problematic around x = -1 and x = 1, where the quadrature
    # nodes are clustered most tightly.  Thus we have the need for the
    # matrix method.
    return @inbounds kernel.(
        gauss_x.x, transpose(gauss_y.x), gauss_x.x_forward, gauss_x.x_backward
    )
end

"""
    check_domain(kernel, x, y)

Check that `(x, y)` lies within `kernel`'s domain and return it.
"""
function check_domain(kernel, x, y)
    xmin, xmax = xrange(kernel)
    xmin ≤ x ≤ xmax || throw(DomainError(x, "x value not in range [$xmin, $xmax]"))

    ymin, ymax = yrange(kernel)
    ymin ≤ y ≤ ymax || throw(DomainError(y, "y value not in range [$ymin, $ymax]"))

    return x, y
end

function compute_uv(Λ, x, y, x₊=1 + x, x₋=1 - x)
    u₊ = x₊ / 2
    u₋ = x₋ / 2
    v = Λ * y
    return u₊, u₋, v
end

"""
    get_symmetrized(kernel, sign)

Construct a symmetrized version of `kernel`, i.e. `kernel(x, y) + sign * kernel(x, -y)`.

!!! warning "Beware!"

    By default, this returns a simple wrapper over the current instance
    which naively performs the sum.  You may want to override this
    to avoid cancellation.
"""
get_symmetrized(kernel::AbstractKernel, sign) = ReducedKernel(kernel, sign)

function get_symmetrized(kernel::LogisticKernel, sign)
    sign == -1 && return LogisticKernelOdd(kernel, sign)
    return Base.invoke(get_symmetrized, Tuple{AbstractKernel,typeof(sign)}, kernel, sign)
end

function get_symmetrized(kernel::RegularizedBoseKernel, sign)
    sign == -1 && return RegularizedBoseKernelOdd(kernel, sign)
    return Base.invoke(get_symmetrized, Tuple{AbstractKernel,typeof(sign)}, kernel, sign)
end

get_symmetrized(::AbstractReducedKernel, sign) = error("cannot symmetrize twice")

function callreduced(kernel::AbstractReducedKernel, x, y, x₊, x₋)
    x, y = check_domain(kernel, x, y)

    # The reduced kernel is defined only over the interval [0, 1], which
    # means we must add one to get the x_plus for the inner kernels.  We
    # can compute this as 1 + x, since we are away from -1.
    x₊ = 1 + x₊

    K₊ = kernel.inner(x, +y, x₊, x₋)
    K₋ = kernel.inner(x, -y, x₊, x₋)
    return K₊ + kernel.sign * K₋
end

(kernel::ReducedKernel)(x, y, x₊, x₋) = callreduced(kernel, x, y, x₊, x₋)

"""
    is_centrosymmetric(kernel)

Return `true` if `kernel(x, y) == kernel(-x, -y)` for all values of `x` and `y`
in range. This allows the kernel to be block-diagonalized,
speeding up the singular value expansion by a factor of 4.  Defaults
to `false`.
"""
iscentrosymmetric(::AbstractKernel) = false
iscentrosymmetric(::LogisticKernel) = true
iscentrosymmetric(::RegularizedBoseKernel) = true
iscentrosymmetric(::AbstractReducedKernel) = false

"""
    (kernel::AbstractKernel)(x, y[, x₊, x₋])

Evaluate `kernel` at point `(x, y)`.

The parameters `x₊` and `x₋`, if given, shall contain the
values of `x - xₘᵢₙ` and `xₘₐₓ - x`, respectively.  This is useful
if either difference is to be formed and cancellation expected.
"""
function (kernel::AbstractKernel)(
    x, y, x₊=x - first(xrange(kernel)), x₋=last(xrange(kernel)) - x
)
    @boundscheck x, y = check_domain(kernel, x, y)
    u₊, u₋, v = compute_uv(kernel.Λ, x, y, x₊, x₋)
    return compute(kernel, u₊, u₋, v)
end

function (kernel::LogisticKernelOdd)(
    x, y, x₊=x - first(xrange(kernel)), x₋=last(xrange(kernel)) - x
)
    result = callreduced(kernel, x, y, x₊, x₋)

    # For x * y around 0, antisymmetrization introduces cancellation, which
    # reduces the relative precision. To combat this, we replace the
    # values with the explicit form
    v_half = kernel.inner.Λ / 2 * y
    xy_small = x * v_half < 1
    cosh_finite = v_half < 85
    return (xy_small && cosh_finite) ? -sinh(v_half * x) / cosh(v_half) : result
end

function (kernel::RegularizedBoseKernelOdd)(
    x, y, x₊=x - first(xrange(kernel)), x₋=last(xrange(kernel)) - x
)
    result = callreduced(kernel, x, y, x₊, x₋)

    # For x * y around 0, antisymmetrization introduces cancellation, which
    # reduces the relative precision.  To combat this, we replace the
    # values with the explicit form
    v_half = kernel.inner.Λ / 2 * y
    xv_half = x * v_half
    xy_small = xv_half < 1
    sinh_range = 1e-200 < v_half < 85
    return xy_small && sinh_range ? -y * sinh(xv_half) / sinh(v_half) : result
end

segments_x(hints::SVEHintsReduced) = symm_segments(segments_x(hints.inner_hints))
segments_y(hints::SVEHintsReduced) = symm_segments(segments_y(hints.inner_hints))

function symm_segments(x)
    x ≈ -reverse(x) || error("segments must be symmetric")
    xpos = x[(begin + length(x) ÷ 2):end]
    iszero(first(xpos)) || (xpos = [0; xpos])
    return xpos
end

"""
    sve_hints(kernel, ε)

Provide discretisation hints for the SVE routines.

Advises the SVE routines of discretisation parameters suitable in
tranforming the (infinite) SVE into an (finite) SVD problem.

See also [`AbstractSVEHints`](@ref).
"""
function sve_hints end

sve_hints(kernel::LogisticKernel, ε) = SVEHintsLogistic(kernel, ε)
sve_hints(kernel::RegularizedBoseKernel, ε) = SVEHintsRegularizedBose(kernel, ε)
sve_hints(kernel::AbstractReducedKernel, ε) = SVEHintsReduced(sve_hints(kernel.inner, ε))

"""
    nsvals(hints)

Upper bound for number of singular values.

Upper bound on the number of singular values above the given threshold, i.e. where
`s[l] ≥ ε * first(s)`.
"""
function nsvals end
function nsvals(hints::SVEHintsLogistic)
    log10_Λ = max(1, log10(hints.kernel.Λ))
    return round(Int, (25 + log10_Λ) * log10_Λ)
end
function nsvals(hints::SVEHintsRegularizedBose)
    log10_Λ = max(1, log10(hints.kernel.Λ))
    return round(Int, 28 * log10_Λ)
end
function nsvals(hints::SVEHintsReduced)
    return (nsvals(hints.inner_hints) + 1) ÷ 2
end

"""
    ngauss(hints)

Gauss-Legendre order to use to guarantee accuracy.
"""
ngauss(hints::SVEHintsLogistic) = hints.ε ≥ 1e-8 ? 10 : 16
ngauss(hints::SVEHintsRegularizedBose) = hints.ε ≥ 1e-8 ? 10 : 16
ngauss(hints::SVEHintsReduced) = ngauss(hints.inner_hints)

"""
    ypower(kernel)

Power with which the ``y`` coordinate scales.
"""
ypower(::AbstractKernel) = 0
ypower(::RegularizedBoseKernel) = 1
ypower(kernel::AbstractReducedKernel) = ypower(kernel.inner)

"""
    conv_radius(kernel)

Convergence radius of the Matsubara basis asymptotic model.

For improved relative numerical accuracy, the IR basis functions on the
Matsubara axis `uhat(basis, n)` can be evaluated from an asymptotic
expression for `abs(n) > conv_radius`.  If `isinf(conv_radius)`, then
the asymptotics are unused (the default).
"""
conv_radius(::AbstractKernel) = Inf
conv_radius(kernel::LogisticKernel) = 40 * kernel.Λ
conv_radius(kernel::RegularizedBoseKernel) = 40 * kernel.Λ
conv_radius(kernel::AbstractReducedKernel) = conv_radius(kernel.inner)

"""
    weight_func(kernel, statistics::Statistics)

Return the weight function for the given statistics.
"""
function weight_func(::AbstractKernel, ::Statistics)
    return x -> ones(eltype(x), size(x))
end
function weight_func(kernel::LogisticKernel, statistics::Statistics)
    if statistics == fermion
        return y -> ones(eltype(y), size(y))
    else
        return y -> 1 ./ tanh.(0.5 * kernel.Λ * y)
    end
end
function weight_func(::RegularizedBoseKernel, statistics::Statistics)
    statistics == boson || error("Kernel is designed for bosonic functions")
    return y -> 1 / y
end


text:
## NOTE: It appears that collect always copies. convert does nothing if they type is already correct.
## Goal here is apply collect when it is needed. But not needlessly copy an array
"""
    tomargs(a::Vector)

converts `a` to type `MxprArgs`. `tomargs` is the identity
if `a` is of type `MxprArgs`.
"""
tomargs(a::Vector) = convert(MxprArgs,a)
tomargs(a::AbstractVector) = convert(MxprArgs,collect(MxprArgT,a)) ## convert should be redundant here
tomargs(a::MxprArgs) = a

## Following includes iterables.
## TODO: We could collect directly into type Any rather than copying
## TODO: test if as is iterable or s.t. There was a discussion about this
tomargs(a) = collect(MxprArgT,a)

"""
    tolist(a::AbstractArray)

returns a Syamta list with elements in `a`. If `a` is already
of type `MxprArgs`, a copy is *not* made. If `a` is of type `ListT`,
`tolist` is the identity.
"""
tolist(a::AbstractArray) = MListA(tomargs(a))
tolist(p::Pair) = MListA(MxprArgType[first(p), last(p)])
tolist(a::ListT) = a

### TODO: decide whether this copies or not.
"""
    tolistoflists(a)

convert a Julia object that can be indexed on two levels to a `List` of `Lists`.
`a` might be, for instance, `Array{Array{T, 1}, 1}` or nested `Tuple`s.
"""
function tolistoflists(ain)
    a = collect(ain)
    nargs = newargs(a)
    for (i,v) in enumerate(a)
       nargs[i] = tolist(v)
    end
    MListA(nargs)
end

tolistfixed(x) = setfixed(tolist(x))

## Both levels are fixed. We may want more fine grained function some day.
function tolistoflistsfixed(ain)
    a = collect(ain)
    nargs = newargs(a)
    for (i,v) in enumerate(a)
       nargs[i] = tolistfixed(v)
    end
    setfixed(MListA(nargs))
end

function tolistoflistsfixed(dict::Dict)
    setfixed(MListA(tomargsoflistfixed(dict)))
end

function tomargsoflistfixed(dict::Dict)
    args = newargs(length(dict))
    i = 0
    for (k,v) in dict
        i += 1
        args[i] = mxprcf(:List,k,v)
    end
    return args
end

"""
   maplist(f, ls::ListT)

maps `f` over `ls`, returing a `ListT`
"""
function maplist(f, ls::ListT)
    nargs = newargs(ls)
    for (i,x) in enumerate(ls)
        nargs[i] = f(x)
    end
    MListA(nargs)
end

"""
    mapmargs(f, inargs::AbstractArray)

map `f` over `inargs` returning an `Array{Any}`.
Arguments to `Mxpr`s are stored in this type. Julia v0.6
no longer creates a container of the same type as `inargs`,
but rather a container of the least common supertype. We don't
want that behavior. So we use this function.
"""
function mapmargs(f, inargs::AbstractArray)
    nargs = newargs(length(inargs))
    map!(f, nargs, inargs)
    return nargs
end

function mapmargs(f, inargs::Tuple)
    nargs = newargs()
    for x in inargs
        push!(nargs, f(x))
    end
    return nargs
end

# not used yet
#totuple(mx::Mxpr{:List}) = tuple(margs(mx)...)

"""
    recursive_copy(x)

just calls `deepcopy`.
"""
recursive_copy(x) = deepcopy(x)

sjcopy(s::T) where {T<:Union{AbstractString,Symbol}} = identity(s)
sjcopy(x) = copy(x)

"""
    mxpr_count_heads(mx::Mxpr, head)

Count the number of first-level elements of mx that are of type Mxpr{head}.
"""
mxpr_count_heads(mx::Mxpr, head) = count(x -> isa(x,Mxpr{head}), mx)

### MList

@doc """
    MList(a...)

is equivalent to mxpr(:List,a...)
"""

@doc """
    MListA(a...)

is equivalent to mxpra(:List,a...)
"""

## This is not really utility, should be in mxpr.jl
@doc """
    MPlus(a...)

is equivalent to mxpr(:Plus,a...)
"""

for sym in ( :List, :Plus, :Times, :Power, :Blank )
    f = Symbol(:M,sym)
    fa = Symbol(:M,sym,:A)
    s = QuoteNode(sym)
    @eval ($f)(a) = mxpr($s,a...)  # Is this really what we want ?
    @eval ($f)(a...) = mxpr($s,a...)
    @eval ($fa)(a) = mxpra($s,a)
end

###

# This is applied at toplevel after expression has been constructed.
# It is applied now to binomial expansion. Needs to be applied more
# generally.
function apply_upvalues_to_args!(mx::Mxpr)
    syms = listsyms(mx)
    goodsyms = Array{SJSym}(undef, 0)
    for sym in syms
        if has_upvalues(sym)
            push!(goodsyms,sym)
        end
    end
    args = margs(mx)
    for sym in goodsyms
        mx = deep_apply_upvalues!(mx,sym)
    end
    return mx
end

function deep_apply_upvalues!(mx::Mxpr,sym)
    args = margs(mx)
    for i in 1:length(mx)
        args[i] = deep_apply_upvalues!(args[i],sym)
    end
    return applyupvalues(mx,sym)
end
deep_apply_upvalues!(x,sym) = x

#using Base.Test
import Base.Test
using Symata

# For use in ../test/
macro testex(expr)
    mx = Expr(:macrocall, Symbol("@symfull"), expr)
    result = eval(mx)
    retresult::Bool = true
    if typeof(result) <: Bool
        retresult = result
    else
        retresult = false
    end
    retresult || @symwarn("Test failed: ", mx, " evaluated to ", retresult)
    Expr(:macrocall, Meta.parse("Test.@test"), retresult)
end


text:
function transform_img(; img, flip_lr=false, flip_ud=false, negative=false)
    if flip_lr
        img = reverse(img, dims=2)
    end

    if flip_ud
        img = reverse(img, dims=1)
    end

    if negative
        img = map(x->typeof(x)(1-x.r,1-x.g,1-x.b), img) # works probably for RGB images only
    end

    return (; bigarrs = (;img, ))
end

fns = (; transform_img)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "UBPF"
version = v"0.0.1" # Temporary fake version, until ubpf provides releases

# Collection of sources required to complete build
sources = [
    # Use jpsamaroo's repo until upstream (iovisor/ubpf) becomes active
    GitSource("https://github.com/jpsamaroo/ubpf.git", "173a58790af1fdefc67aea6e22aaec0d6d8b4d67")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/ubpf/vm
if [[ "${target}" == *mingw* ]]; then
sed -i 's/-fPIC//g' Makefile
elif [[ "${target}" == *bsd* ]]; then
sed -i 's#endian.h#sys/endian.h#g' *.c *.h
fi
make -j $(nproc)
make PREFIX=${prefix} install
cd ..
install_license LICENSE-APACHE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(x->!(Sys.iswindows(x)||Sys.isapple(x)), supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libubpf", :libubpf)
]

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

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


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

__precompile__(true)

module SolidStateDetectors

using LinearAlgebra
using Random
using Statistics

using Adapt
using ArraysOfArrays
using FillArrays
using Formatting
using GPUArrays
using Interpolations
using IntervalSets
using JSON
using KernelAbstractions
using LaTeXStrings
using ParallelProcessingTools
using ProgressMeter
using RadiationDetectorSignals
using RecipesBase
using Requires
using Rotations
using StaticArrays
using StatsBase
using Unitful
using UnitfulAtomic
using UnitfulRecipes
using YAML

include("ConstructiveSolidGeometry/ConstructiveSolidGeometry.jl")
using .ConstructiveSolidGeometry
using .ConstructiveSolidGeometry:
            CylindricalPoint, CartesianPoint, AbstractCoordinatePoint, _convert_point,
            CartesianVector, CylindricalVector, AbstractCoordinateVector,
            Cartesian, Cylindrical, AbstractCoordinateSystem, CoordinateSystemType,
            Geometry, AbstractGeometry, AbstractSurfacePrimitive,
            parse_rotation_matrix, parse_translate_vector, parse_CSG_transformation,
            transform, CSG_dict, Transformations, combine_transformations,
            ConfigFileError, _parse_value,
            LengthQuantity, AngleQuantity

import .ConstructiveSolidGeometry: sample, to_internal_units, from_internal_units
export CartesianPoint, CartesianVector, CylindricalPoint

import Clustering
import DataStructures
import Distributions
import IntervalSets
import Tables
import TypedTables

import Base: size, sizeof, length, getindex, setindex!, axes, getproperty, broadcast,
             range, ndims, eachindex, enumerate, iterate, IndexStyle, eltype, in, convert,
             show, print, println, display, +, -, &

export SolidStateDetector
export SSD_examples

export Grid, GridPoint

export ElectricPotential, PointTypes, EffectiveChargeDensity, DielectricDistribution, WeightingPotential, ElectricField
export apply_initial_state!
export calculate_electric_potential!, calculate_weighting_potential!, calculate_electric_field!, calculate_drift_fields!
export ElectricFieldChargeDriftModel, ADLChargeDriftModel
export get_active_volume, is_depleted
export calculate_stored_energy, calculate_mutual_capacitance, calculate_capacitance_matrix
export simulate_waveforms
export Simulation, simulate!
export Event, drift_charges!
export add_baseline_and_extend_tail
export NBodyChargeCloud

const SSDFloat = Union{Float16, Float32, Float64}

include("Units.jl")
include("Axes/DiscreteAxis.jl")
include("World/World.jl")
include("Grids/Grids.jl")

include("Types/Types.jl")

include("MaterialProperties/MaterialProperties.jl")
include("Config/Config.jl")
include("ChargeDensities/ChargeDensities.jl")
include("ImpurityDensities/ImpurityDensities.jl")
include("ChargeDriftModels/ChargeDriftModels.jl")
include("SolidStateDetector/DetectorGeometries.jl")

include("PotentialSimulation/PotentialSimulation.jl")
include("PotentialSimulation/ConvergenceGPU.jl")

include("ElectricField/ElectricField.jl")

include("ChargeCloudModels/ChargeCloudModels.jl")
include("ChargeDrift/ChargeDrift.jl")
include("SignalGeneration/SignalGeneration.jl")

include("ChargeStatistics/ChargeStatistics.jl")
include("ChargeClustering/ChargeClustering.jl")

include("Simulation/Simulation.jl")
include("Event/Event.jl")
include("MCEventsProcessing/MCEventsProcessing.jl")

include("IO/IO.jl")

include("examples.jl")

include("PlotRecipes/PlotRecipes.jl")
export @P_str # protected strings to overwrite plot labels with units

function __init__()
    @require HDF5="f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" begin
        @require LegendHDF5IO="c9265ca6-b027-5446-b1a4-febfa8dd10b0" begin
            include("IO/hdf5_specific.jl")
        end
        include("MCEventsProcessing/MCEventsProcessing_hdf5.jl")
    end
    @require CUDAKernels="72cfdca4-0801-4ab0-bf6a-d52aa10adc57" begin
        get_device(::Type{CUDAKernels.CUDA.CuArray}) = CUDAKernels.CUDADevice()
    end
    @require ROCKernels="7eb9e9f0-4bd3-4c4c-8bef-26bd9629d9b9" begin
        get_device(::Type{ROCKernels.AMDGPU.ROCArray}) = CUDAKernels.ROCDevice()
    end
end

end # module


text:
struct WeekDayEffect1{Tα} <: ContinuousMultivariateDistribution
	α::Tα
end
WeekDayEffect1() = WeekDayEffect1(ones(Int64, 7))

Distributions.rand(rng::AbstractRNG, d::WeekDayEffect1{Tα}) where {Tα} = begin
	weekdayeffect_simplex = rand( Dirichlet(d.α) )
	weekdayeffect_simplex .* length(d.α)
end

Distributions.logpdf(d::WeekDayEffect1{Tα}, x::AbstractVector{T}) where {Tα,T} =
    logpdf( Dirichlet( d.α ), x ./ length(d.α) )

Bijectors.bijector(d::WeekDayEffect1) =
	Bijectors.bijector(Dirichlet(d.α)) ∘ Bijectors.Scale(1/7, dim=Val(1))

#Bijectors.bijector(Dirichlet( d.α ))
@model function test(y, d, ::Type{TV} = Vector{Float64}) where TV
	dwe ~ WeekDayEffect1()
	e ~ Normal(4, 0.1)
	ŷ = TV(undef, length(y))
	for i in eachindex(ŷ)
		ŷ = e + dwe[d[i]]
	end
	ŷ ~ MvNormal(y, 0.1)
	return ŷ
end

m = 10
d = repeat(1:7, outer=m)
n = length(d)
y = repeat(1:7, outer=m) + randn(n) * 0.1

m = test(y, d)
plot!(m())
c = sample(m, NUTS(1000, 0.8), 2000; progress=true)

gq = generated_quantities(m, c[1001:end])|>vec
ŷ = mean(gq)
plot(y)
plot!(ŷ)

plot(group(c, "dwe"))

logpdf(Dirichlet(ones(7)), ones(7) / 7.)
logpdf(WeekDayEffect1(ones(7)), ones(7))

bs = Bijectors.Scale(1/7, dim=Val(1))
bs isa Bijectors.Bijector
x = 1.
y = bs(x)
inv(bs)(y) ≈ x

bd = bijector(Dirichlet(ones(7)))
bd isa Bijectors.Bijector
x = ones(7)/7
y = bd(x)
inv(bd)(y) ≈ x

b = bd ∘ bs
x = ones(7)
y = b(x)
inv(b)(y) ≈ x

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
# === Misc Utility Functions === #
const LOG_2_PI = log(2 * pi)

"""Log-gamma function."""
lgamma(x) = logabsgamma(x)[1]

"""
Sample from a list of log probabilities. Overwrites
vector so that no new memory is allocated.
"""
sample_logprobs!(log_probs) = sample(pweights(softmax!(log_probs)))

"""
Log normalizer of 1D normal distribution in information form.
given potential `m` and precision `v`.

N(x | m, v) = exp{ m⋅x - .5v⋅x² - log Z(m, v) }
"""
gauss_info_logZ(m, v) = 0.5 * (LOG_2_PI - log(v) + m * m / v)


text:
using BinDeps
WORD_SIZE = 32
@BinDeps.setup

group = library_group("libav")

deps = [
    libavcodec = library_dependency("libavcodec", aliases=[["libavcodec-ffmpeg.so.","libavcodec.","libavcodec.so.","libavcodec.ffmpeg.so.","avcodec-"].*["53" "54" "55" "56"]...], group = group)
    libavformat = library_dependency("libavformat", aliases=[["libavformat-ffmpeg.so.","libavformat.","libavformat.so.","libavformat.ffmpeg.so.","avformat-"].*["53" "54" "55" "56"]...], group = group)
    libavutil = library_dependency("libavutil", aliases=[["libavutil-ffmpeg.so.", "libavutil.","libavutil.so.", "libavutil.ffmpeg.so.", "avutil-"].*["51" "52" "54"]...], group = group)
    libswscale = library_dependency("libswscale", aliases=[["libswscale-ffmpeg.so.","libswscale.","libswscale.so.","libswscale.ffmpeg.so.","swscale-"].*["2" "3"]...], group = group)
    libavfilter = library_dependency("libavfilter", aliases=[["libavfilter-ffmpeg.so.","libavfilter.","libavfilter.so.","libavfilter.ffmpeg.so.","avfilter-"].*["2" "3" "4" "5"]...], group = group)
    libavdevice = library_dependency("libavdevice", aliases=[["libavdevice-ffmpeg.so.","libavdevice.","libavdevice.so.","libavdevice.ffmpeg.so.","avdevice-"].*["53" "54" "55" "56"]...], group = group)
]

libav_libs = [libavutil, libavcodec, libavformat, libavfilter, libswscale, libavdevice]

# if (have_avresample = dlopen_e("libavresample") != C_NULL)
#     libavresample = library_dependency("libavresample", aliases=[["libavresample.so.","avresample-"].*["1"]...], group = group)
#     push!(deps, libavresample)
# end

# if (have_swresample = dlopen_e("libswresample") != C_NULL)
#     libswresample = library_dependency("libswresample", aliases=[["libswresample.so.","libswresample.ffmpeg.so.","swresample-"].*["0"]...], group = group)
#     push!(deps, libswresample)
#     push!(libav_libs, libswresample)
# end

if is_windows()
    provides(Binaries, URI("http://www.daniel-thiem.de/ffmpeg/ffmpeg-2.2.3-win64-shared.7z"),
             libav_libs,
             os = :Windows,
             unpacked_dir="ffmpeg-2.2.3-win64-shared/bin")
end

if is_apple()
    using Homebrew
    provides( Homebrew.HB, "ffmpeg", libav_libs, os = :Darwin )
end

# System Package Managers
apt_packages = Dict(
    "libavcodec-extra-53"   => libavcodec,
    "libavcodec53"          => libavcodec,
    "libavcodec-extra-54"   => libavcodec,
    "libavcodec54"          => libavcodec,
    "libavcodec-extra-55"   => libavcodec,
    "libavcodec55"          => libavcodec,
    "libavdevice53"         => libavdevice,
    "libavfilter2"          => libavfilter,
    "libavfilter3"          => libavfilter,
    "libavfilter4"          => libavfilter,
    "libavformat53"         => libavformat,
    "libavformat54"         => libavformat,
    "libavformat55"         => libavformat,
    #"libavresample1"        => libavresample,
    "libavutil51"           => libavutil,
    "libavutil52"           => libavutil,
    "libswscale2"           => libswscale,

    ## Available from https://launchpad.net/~jon-severinsson/+archive/ubuntu/ffmpeg
    "libavcodec55-ffmpeg"   => libavcodec,
    "libavdevice55-ffmpeg"  => libavdevice,
    "libavfilter4-ffmpeg"   => libavfilter,
    "libavformat55-ffmpeg"  => libavformat,
    "libavutil52-ffmpeg"    => libavutil,
    #"libswresample0-ffmpeg" => libswresample,
    "libswscale2-ffmpeg"    => libswscale,
    # ffmpeg is again available in standard packages as of ubuntu 15.04
    "libavcodec-ffmpeg56"   => libavcodec,
    "libavdevice-ffmpeg56"  => libavdevice,
    "libavfilter-ffmpeg5"   => libavfilter,
    "libavformat-ffmpeg56"  => libavformat,
    "libavutil-ffmpeg54"    => libavutil,
    "libswscale-ffmpeg3"    => libswscale,
    #"libavresample-ffmpeg2" => libavresample,
)

#if have_swresample
#    apt_packages["libswresample0-ffmpeg"] = libswresample
#end

#if have_avresample
#    apt_packages["libavresample1"] = libavresample
#end

provides(AptGet,
         apt_packages)

provides(Yum,
         Dict("ffmpeg" => libav_libs))

provides(Pacman,
         Dict("ffmpeg" => libav_libs))

provides(Sources,
         URI("http://www.ffmpeg.org/releases/ffmpeg-2.3.2.tar.gz"),
         libav_libs)

provides(BuildProcess, Autotools(configure_options=["--enable-gpl"]), libav_libs, os = :Unix)

@BinDeps.install Dict(
    :libavcodec => :libavcodec,
    :libavformat => :libavformat,
    :libavutil => :libavutil,
    :libswscale => :libswscale,
    :libavfilter => :libavfilter,
    :libavdevice => :libavdevice
)

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
# Adapted from:
# SOSDEMO6 --- MAX CUT
# Section 3.6 of SOSTOOLS User's Manual

using Test
using SumOfSquares
using DynamicPolynomials

function maxcut_test(optimizer, config::MOIT.Config)
    atol = config.atol
    rtol = config.rtol

    @polyvar x[1:5]

    # Number of cuts
    f = 2.5 - 0.5*x[1]*x[2] - 0.5*x[2]*x[3] - 0.5*x[3]*x[4] - 0.5*x[4]*x[5] - 0.5*x[5]*x[1]

    # Boolean constraints
    bc = vec(x).^2 .- 1

    model = _model(optimizer)

    Z = monomials(x, 0:1)
    @variable(model, p1, SOSPoly(Z))

    Z = monomials(x, 0:2)
    @variable(model, p[1:5], Poly(Z))

    @testset "with γ=$γ it should be $(feasible ? "feasible" : "infeasible")" for (γ, feasible) in [(3.9, false), (4.1, true)]
        cref = @constraint(model, p1*(γ-f) + dot(p, bc) - (γ-f)^2 in SOSCone())

        JuMP.optimize!(model)

        if feasible
            @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
        else
            @test JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE
        end
        delete(model, cref)
    end
end
sd_tests["maxcut"] = maxcut_test

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test
using MozillaCACerts_jll

@testset "MozillaCACerts_jll" begin
    @test isfile(MozillaCACerts_jll.cacert)
end


text:
"""
	degree(s,node)

Return the degree of the given node in the givn stream.

```math
d \\left(v \\right)= \\left| \\mathcal{N} \\left(v\\right) \\right|
```

Where

```math
\\mathcal{N} \\left(v\\right)=\\left\\{ \\left(t,u\\right),\\ \\left(t,uv \\right) \\in E}
```

is the neighborhood of node v.
Note: If the stream has a duration of 0, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function degree(s::AbstractUndirectedStream,node::AbstractString)
	duration(s)!=0 ? duration(neighborhood(s,node))/duration(s) : 0.0
end

"""
	degree(s,node,tc)

Return the degree of the given node in the given stream using a TimeCursor.
Note: Not implemented yet...
"""
function degree(s::AbstractUndirectedStream,node::AbstractString,tc::TimeCursor)
	throw("Not Implemented")
end

"""
	degree(s,node,t)

Return the instantaneous degree of the given node at the given time in the given stream.
```math
d_t \\left(v \\right)= \\left| \\mathcal{N}_t \\left(v\\right) \\right|
```

Where

```math
\\mathcal{N}_t\\left(v\\right)=\\left\\{ u,\\ \\left(t,uv \\right) \\in E}
```

is the instantaneous neighborhood of node v.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function degree(s::AbstractUndirectedStream,node::AbstractString,t::Float64)
	length(neighborhood(s,node,t))
end

"""
	degree(s,t)

Return the degree at t of the stream s.
```math
d\\left(t\\right) = \\sum_v \\frac{d_t \\left(v \\right)}{\\left|V\\right||}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function degree(s::AbstractUndirectedStream,t::Float64)
	length(s.V) != 0 ? 1.0/length(s.V)*sum([degree(s,v,t) for v in s.V]) : 0.0
end

"""
	degree(s)

Return the degree of the stream s.
```math
d\\left(S\\right)= \\sum_v \\frac{1}{\\left|V\\right|}d\\left(v\\right)
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function degree(s::AbstractUndirectedStream)
	length(s.V) != 0 ? 1.0/length(s.V)*sum([degree(s,v) for v in s.V]) : 0.0
end

"""
	degree(s,tc)

Return the degree of stream s using a TimeCursor.
In this case, the degree is computed using:
```math
d\\left(S\\right)= \\int_t \\frac{1}{\\left|T\\right|}d\\left(t\\right) dt
```
Note: If the stream has no node or a duration of 0, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function degree(s::AbstractUndirectedStream,tc::TimeCursor)
	((duration(s)==0) | (length(s.V)==0)) && 0.0
	start!(tc)
	δ=collect(values(degrees(tc.S)))
	d::Float64=sum(δ)*duration(tc.S)
	while tc.S.t1 < s.T.list[1][2]
        next!(tc)
        δ=collect(values(degrees(tc.S)))
        d+=sum(δ)*duration(tc.S)
    end
    1.0/(duration(s)*length(s.V))*d
end

"""
	average_node_degree(s)

Return the average node degree of the given stream graph.
```math
d\\left(V\\right)=\\frac{1}{n} \\sum_{v \\in V} n_v d\\left( v \\right)
```
Note: If card(W)=0, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function average_node_degree(s::StreamGraph)
	duration(s.W)!=0 ? 1.0/duration(s.W)*sum([length(times(s,v))*degree(s,v) for v in s.V]) : 0.0
end

"""
	average_node_degree(ls)

Return the average node degree of the given link stream.
For link streams, the average node degree simplifies to:
```math
d\\left(L\\right) = \\frac{2 m}{n}
```
Note: If the link stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function average_node_degree(ls::LinkStream)
	number_of_nodes(ls)!=0 ? 2.0*number_of_links(ls)/number_of_nodes(ls) : 0.0
end

"""
	average_time_degree(s)

Return the average time degree of the given stream.
```math
d\\left(T\\right)=\\int_t \\frac{\\left|V_t\\right|}{\\left|W\\right|} d\\left(t\\right) dt
```
Note: If the node duration is 0, this function returns 0.

Warning: This function is very inefficient, use `average_time_degree(s,tc)` if possible.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function average_time_degree(s::AbstractStream)
    k = node_duration(s)
    τ = times(s)
    k != 0 ? 1.0/k*sum([node_contribution(s,0.5*(t[2]+t[1]))*degree(s,0.5*(t[2]+t[1]))*(t[2]-t[1]) for t in zip(τ[1:end-1],τ[2:end])]) : 0.0
end

"""
	average_time_degree(s,tc)

Return the average time degree of the given stream using a TimeCursor.
"""
function average_time_degree(s::AbstractStream,tc::TimeCursor)
	length(s.V)==0 && 0.0
	dW=duration(s.W)
	dW==0 && 0.0
	start!(tc)
	δ=collect(values(degrees(tc.S)))
	d::Float64=sum(δ)/length(s.V)*number_of_nodes(tc)*duration(tc.S)
    while tc.S.t1 < s.T.list[1][2]
        next!(tc)
        δ=collect(values(degrees(tc.S)))
        d+=sum(δ)/length(s.V)*number_of_nodes(tc)*duration(tc.S)
    end
    1.0/dW*d
end

"""
	degrees(s)

Return the degree distribution of the graph corresponding to the given state.
"""
function degrees(S::State)
	l=String[]
	for ll in S.links
    	push!(l,ll[1])
    	push!(l,ll[2])
	end
	ul=unique(l)
	Dict([(i,count(x->x==i,l)) for i in ul])
end

"""
	average_degree(s)

Return the average degree of the graph corresponding to the given state.
"""
function average_degree(S::State)
	d=collect(values(degrees(S)))
	length(d)==0 && 0.0
	sum(d)/length(d)
end

"""
	degrees(tc)

Return the degree distribution of the cursor's current state.
"""
degrees(tc::TimeCursor)=degrees(tc.S)

"""
	average_degree(tc)

Return the average degree of the cursor's current state.
"""
average_degree(tc::TimeCursor)=average_degree(tc.S)

"""
	expected_degree(s,node)

Return the expected degree of the given node in the given stream.
Note: Not implemented yet...
"""
function expected_degree(s::AbstractStream,node::AbstractString)
	throw("Not Implemented")
end

"""
	expected_degree(s,node,tc)

Return the expected degree of the given node in the given stream using a TimeCursor.
```math
\\hat{d}\\left(v\\right)=\\int_t \\frac{d_t \\left(v\\right)}{\\left|T_v\\right|} dt
```

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function expected_degree(s::AbstractStream,node::AbstractString,tc::TimeCursor)
	Tv=times(s,node)
	length(Tv)==0 && 0.0
	start!(tc)
	δ=degrees(tc.S)
	d::Float64=0
	if haskey(δ,node)
		d+=δ[node]*duration(tc.S)
	end
    while tc.S.t1 < s.T.list[1][2]
        next!(tc)
        δ=degrees(tc.S)
        if haskey(δ,node)
        	d+=δ[node]*duration(tc.S)
        end
    end
    1.0/length(Tv)*d
end

"""
	expected_degree(s,node,t)

Return the expected degree of the given node at the given time in the given stream.
Note: Not implemented yet...
"""
function expected_degree(s::AbstractStream,node::AbstractString,t::Float64)
	throw("Not Implemented")
end

"""
	expected_degree(s,node,t,tc)

Return the expected degree of the given node at the given time in the given
stream using a TimeCursor.
Note: Not implemented yet...
"""
function expected_degree(s::AbstractStream,node::AbstractString,t::Float64,tc::TimeCursor)
	throw("Not Implemented")
end

"""
	average_expected_degree(s)

Return the average expected degree of the given stream.
```math
\\hat{d}\\left(S\\right)=\\frac{2 m}{n}
```
Note: If the stream has no node, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function average_expected_degree(s::AbstractStream)
	number_of_nodes(s)!=0 ? 2 * number_of_links(s) / number_of_nodes(s) : 0.0
end


text:
abstract type PartiallySpecifiedRange{T} <: AbstractUnitRange{Int} end
abstract type ZeroClampedRange{T} <: PartiallySpecifiedRange{T} end
Base.isempty(::PartiallySpecifiedRange) = false

"""
    SingleValuedRange(n::Int)

The range `n:n`.
"""
struct SingleValuedRange <: AbstractUnitRange{Int}
    n :: Int
end
Base.isempty(::SingleValuedRange) = false
Base.first(x::SingleValuedRange) = x.n
Base.last(x::SingleValuedRange) = x.n
Base.length(x::SingleValuedRange) = 1
Base.show(io::IO, x::SingleValuedRange) = print(io, repr(x.n), ":", repr(x.n))

function SingleValuedRange(r::AbstractUnitRange{<:Integer})
    length(r) == 1 || throw(ArgumentError("range must contain only one value"))
    SingleValuedRange(first(r))
end

"""
    ZeroTo(l::Int)

The range `0:l` for an `l ≥ 0`.
"""
struct ZeroTo{T} <: ZeroClampedRange{T}
    l :: Int

    function ZeroTo{T}(l) where {T}
        ensure_nonnegative(l)
        new{T}(l)
    end
end
ZeroTo(l) = ZeroTo{false}(l)
Base.first(::ZeroTo) = 0
Base.last(r::ZeroTo) = r.l
Base.show(io::IO, r::ZeroTo) = print(io, "0:", repr(r.l))

# conversion to false is always possible
ZeroTo{false}(z::ZeroTo) = ZeroTo{false}(z.l)
function ZeroTo{false}(z::AbstractUnitRange{<:Integer})
    first(z) == 0 || throw(ArgumentError("range does not begin at zero"))
    ZeroTo{false}(last(z))
end

function ZeroTo(r::AbstractUnitRange{<:Integer})
    first(r) == 0 || throw(ArgumentError("range must start at zero"))
    ZeroTo(last(r))
end

"""
    ToZero(l::Int)

The range `-l:0` for an `l ≥ 0`.
"""
struct ToZero{T} <: ZeroClampedRange{T}
    l :: Int

    function ToZero{T}(l) where {T}
        ensure_nonnegative(l)
        new{T}(l)
    end
end
ToZero(l) = ToZero{false}(l)
Base.first(r::ToZero) = -r.l
Base.last(::ToZero) = 0
Base.show(io::IO, r::ToZero) = print(io,"-",repr(r.l),":0")

# conversion to false is always possible
ToZero{false}(z::ToZero) = ToZero{false}(z.l)
function ToZero{false}(z::AbstractUnitRange{<:Integer})
    last(z) == 0 || throw(ArgumentError("range does not end at zero"))
    ToZero{false}(-first(z))
end

function ToZero(r::AbstractUnitRange{<:Integer})
    last(r) == 0 || throw(ArgumentError("range must end at zero"))
    ToZero(-first(r))
end

"""
    FullRange(l::Int)

The range `-l:l` for an `l ≥ 0`.
"""
struct FullRange{T} <: PartiallySpecifiedRange{T}
    l :: Int

    function FullRange{T}(l) where {T}
        ensure_nonnegative(l)
        new{T}(l)
    end
end
FullRange(l) = FullRange{false}(l)
Base.first(r::FullRange) = -r.l
Base.last(r::FullRange) = r.l
Base.show(io::IO, r::FullRange) = print(io, repr(-r.l),":",repr(r.l))

# conversion to false is always possible
FullRange{false}(z::FullRange) = FullRange{false}(z.l)
function FullRange{false}(z::AbstractUnitRange{<:Integer})
    last(z) == -first(z) || throw(ArgumentError("first element of range must be negative the last element"))
    FullRange{false}(last(z))
end

function FullRange(r::AbstractUnitRange{<:Integer})
    first(r) == -last(r) || throw(ArgumentError("first element of range must be negative the last element"))
    FullRange(last(r))
end

Base.intersect(a::FullRange, b::FullRange) = FullRange(min(maximum(a), maximum(b)))
Base.intersect(a::T, b::FullRange) where {T<:ZeroClampedRange} = T(min(a.l, b.l))
Base.intersect(a::FullRange, b::T) where {T<:ZeroClampedRange} = T(min(a.l, b.l))

Base.intersect(a::T, b::T) where {T<:ZeroClampedRange} = T(min(a.l, b.l))
Base.intersect(a::ZeroClampedRange, b::ZeroClampedRange) = ZeroTo(0)

function Base.intersect(a::SingleValuedRange, b::AbstractUnitRange{<:Integer})
    first(a) in b || return nothing
    return a
end
Base.intersect(b::AbstractUnitRange{<:Integer}, a::SingleValuedRange) = intersect(a, b)
function Base.intersect(a::SingleValuedRange, b::SingleValuedRange)
    first(a) == first(b) || return nothing
    return a
end

Base.intersect(a::SingleValuedRange, b::PartiallySpecifiedRange) = intersect(a, UnitRange(b))
Base.intersect(a::PartiallySpecifiedRange, b::SingleValuedRange) = intersect(UnitRange(a), b)

Base.promote_rule(::Type{<:PartiallySpecifiedRange}, ::Type{<:PartiallySpecifiedRange}) = UnitRange{Int}
Base.promote_rule(::Type{<:PartiallySpecifiedRange}, ::Type{SingleValuedRange}) = UnitRange{Int}


text:
using Test, PolynomialBases, FastGaussQuadrature

if !haskey(ENV, "JULIA_PKGEVAL") # sympy is not install on https://github.com/JuliaComputing/NewPkgEval.jl
  import SymPy
  x = SymPy.symbols("x")

  @test 0 == SymPy.simplify( legendre(x, 0) - 1 )
  @test 0 == SymPy.simplify( legendre(x, 1) - x )
  @test 0 == SymPy.simplify( legendre(x, 2) - ( 3x^2 - 1 ) / 2 )
  @test 0 == SymPy.simplify( legendre(x, 3) - ( 5x^3 - 3x) / 2 )
  @test 0 == SymPy.simplify( legendre(x, 4) - ( 35x^4 - 30x^2 + 3 ) / 8 )
  @test 0 == SymPy.simplify( legendre(x, 5) - ( 63x^5 - 70x^3 + 15x ) / 8 )
  @test 0 == SymPy.simplify( legendre(x, 6) - ( 231x^6 - 315x^4 + 105x^2 - 5 ) / 16 )

  @test_skip legendre(x, 4)
end
@inferred legendre(10., 4)

# Gauss Legendre nodes and weights
for p in 0:20
    x1, w1 = gausslegendre(p+1)
    @inferred PolynomialBases.gauss_legendre_nodes_and_weights(p)
    x2, w2 = PolynomialBases.gauss_legendre_nodes_and_weights(p)
    @test x1 ≈ x2 atol=1.e-14
    @test w1 ≈ w2 atol=1.e-14
end

# Lobatto Legendre nodes and weights
for p in 1:20
    x1, w1 = gausslobatto(p+1)
    @inferred PolynomialBases.lobatto_legendre_nodes_and_weights(p)
    x2, w2 = PolynomialBases.lobatto_legendre_nodes_and_weights(p)
    @test x1 ≈ x2 atol=1.e-14
    @test w1 ≈ w2 atol=1.e-14
end

# Derivative and mass matrices for Gauss Legendre nodes and weights
# NOTE: Only some tests for low polynomial degrees since the Vandermonde matrix
#       becomes very ill-conditioned for higher values of p.
for p in 0:8
    basis = GaussLegendre(p)
    V = legendre_vandermonde(basis)

    Dhat = legendre_D(p)
    @test basis.D ≈ V*Dhat/V atol=1.e-13

    Mhat = legendre_M(p)
    @test Diagonal(basis.weights) ≈ V'\Mhat/V atol=1.e-13
end


text:
# binocdf is the cdf for the Binomial distribution and binoccdf is its complement.
# Distributions.jl implements the cdf of its Binomial by calling out to Rmath which
# is a Julia package that wraps the corresponding functions used in R. Since we do not want to
# load the Distributions lib when starting Autotest we want a pure Julia implementation
# of binocdf. We base this code on the pure Ruby code for incomplete beta function available here:
#   http://rubydoc.info/gems/rubystats/0.2.3/Rubystats/SpecialMath#incomplete_beta-instance_method
#

# Cumulative Distribution Function evaluated at x of the Binomial(n, p) distribution:
function binocdf(x, n, p)
  x = ifloor(x)
  if x == n
    return 1.0
  else
    incompletebeta(1-p, n - x, x + 1)
  end
end

# Complementary cdf of Binomial distr func:
binoccdf(x, n, p) = 1.0 - binocdf(x, n, p)

# incompletebeta(x, a, b) evaluates incomplete beta function, here a, b > 0 and 0 <= x <= 1.
function incompletebeta(x, a, b)
  if x == 0
    return 0.0
  elseif x == 1
    return 1.0
  else
    beta_gam = exp( a * log(x) + b * log(1.0 - x) - log_beta(a, b) )
    if x < ((a+1) / (a+b+2))
      return (beta_gam * beta_fraction(x, a, b) / a)
    else
      return (1.0 - (beta_gam * beta_fraction(1.0 - x, b, a) / b))
    end
  end
end

const LOG_GAMMA_X_MAX_VALUE = 2.55e292

function log_beta(a, b)
  if (a <= 0.0) || (b <= 0.0) || (a + b) > LOG_GAMMA_X_MAX_VALUE
    return 0.0
  else
    return (lgamma(a) + lgamma(b) - lgamma(a + b))
  end
end

const XMININ = 2.23e-303
const PRECISION = 8.88e-016

function beta_fraction(x, a, b, MaxIterations = 200)
  c = 1.0
  sum_ab  = a + b
  a_plus  = a + 1.0
  a_minus = a - 1.0
  h = 1.0 - sum_ab * x / a_plus
  if abs(h) < XMININ
    h = XMININ
  end
  h     = 1.0 / h
  frac  = h
  m     = 1
  delta = 0.0

  while (m <= MaxIterations) && (abs(delta - 1.0) > PRECISION)
    m2 = 2 * m
    # even index for d
    d = m * (b - m) * x / ( (a_minus + m2) * (a + m2))
    h = 1.0 + d * h
    if abs(h) < XMININ
      h = XMININ
    end
    h = 1.0 / h
    c = 1.0 + d / c
    if abs(c) < XMININ
      c = XMININ
    end
    frac *= (h * c)
    # odd index for d
    d = -(a + m) * (sum_ab + m) * x / ((a + m2) * (a_plus + m2))
    h = 1.0 + d * h
    if abs(h) < XMININ
      h = XMININ
    end
    h = 1.0 / h
    c = 1.0 + d / c
    if abs(c) < XMININ
      c = XMININ
    end
    delta = h * c
    frac *= delta
    m += 1
  end

  return frac
end

# Compare to Distributions.jl
using Distributions

# call once to compile things
binocdf(1, 10, 0.3)
binocdf(1.0, 10.0, 0.3)
binocdf(2.0, 10, 0.3)

# Sample lower values with higher prob through a log transform.
# Often more realistic since people tend to use lower ranges/values
# frequently.
function logsample(min, max)
  shift = 1.0 - min
  logmin = log(min + shift)
  logmax = log(max + shift)
  exp(logmin + (logmax - logmin) * rand()) - shift
end

NumReps = 1e6
ourtimesum = 0.0
dtimesum = 0.0

for rep in 1:ifloor(NumReps)
  n = ifloor(logsample(2, 1e6))
  x = rand() * n # rand(0:n)
  p = rand()

  tic()
  ourres = binocdf(x, n, p)
  ourtimesum += toq()

  tic()
  dres = cdf(Binomial(n, p), x)
  dtimesum += toq()

  if !isapprox(ourres, dres)
    println("Failed for x = $x, n = $n, p = $p: our = $ourres, Distr.jl = $dres")
  end
end

println("Our average binocdf performance: ", ourtimesum / NumReps)
println("Distributions.jl average binocdf performance: ", dtimesum / NumReps)
println("Performance diff: ", round(ourtimesum / dtimesum, 2), "x")


text:
import Flux: Data.DataLoader, stack

generate_2d_data(T) = generate_data(T; stacked=false)
generate_3d_data(T) = generate_data(T; stacked=true)
generate_3d_ct_data(T) = generate_data(T; stacked=true, ct=true)

function generate_ts_data(T::DataType=Float32)
  N = 48
  tsx = T.(collect(0:N-1))
  tsy = T.(collect(1:N))
  data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
  data_x = [reshape([T(data_x[1][i]),T(data_x[2][i])],2,1) for i in 1:N]# |> f32
  data_y = [reshape([T(y)],1,1) for y in sin.(range(0,stop=6π,length=N))]# |> f32
  dl = DataLoader((zip(data_x,tsx), (data_y,tsy)), batchsize=N)
end

function generate_data(T::DataType=Float32; stacked=true, ct=false)
    in_features = 2
    out_features = 1
    N = 48
    data_x = [sin.(range(0,stop=3π,length=N)), cos.(range(0,stop=3π,length=N))]
    data_x = [reshape([T(data_x[1][i]),T(data_x[2][i])],2,1) for i in 1:N]# |> f32
    data_y = [reshape([T(y)],1,1) for y in sin.(range(0,stop=6π,length=N))]# |> f32
    dl = DataLoader((data_x, data_y), batchsize=N)
    @show length(dl)
    fx, fy = first(dl)
    @show size(fx)
    @show size(fx[1])
    @show size(fy)
    @show size(fy[1])
    fig = plot([x[1,1] for x in fx], label="x1")
    plot!(fig, [x[2,1] for x in fx], label="x2")
    plot!(fig, [y[1,1] for y in fy], label="y1")
    display(fig)

    if stacked == false && ct == false
      return dl
    end

    data_x, data_y = stacked == true ? (stack(data_x,2), stack(data_y,2)) : (data_x, data_y)
    saveat = collect(T(0):size(data_x,2)-1)
    #data_x, data_y = ct == true ? ((data_x, saveat), data_y) : (data_x, data_y)
    batchsize = stacked == true ? 1 : N
    batchsize = 48
    dx = [(data_x[:,s:s+batchsize-1,:],saveat[s:s+batchsize-1]) for s in 1:batchsize:size(data_x,2)-batchsize+1]
    dy = [data_y[:,s:s+batchsize-1,:] for s in 1:batchsize:size(data_y,2)-batchsize+1]

    @show summary(dx)
    @show summary(dy)
    dl = Flux.Data.DataLoader((dx, dy), batchsize=1)

    fxt, fy = first(dl)
    fxt, fy = fxt[1], fy[1]
    @show summary(fxt)
    @show summary(fy)
    fx,tx = fxt
    @show size(fx)
    @show size(tx)
    @show size(fy)
    @show size(fy)

    dl
end

function mycb(p,l,ŷ,y;doplot=true)
  #display(l)
  if doplot
    y = ndims(y) < 3 ? stack(y,2) : y
    ŷ = ndims(ŷ) < 3 ? stack(ŷ,2) : ŷ

    if size(y,1) > 3
      hy = heatmap(1:size(y,2),1:size(y,1),(@view y[:,:,1]),title="data")
      hŷ = heatmap(1:size(ŷ,2),1:size(ŷ,1),(@view ŷ[:,:,1]),title="prediction")
      he = heatmap(1:size(y,2),1:size(y,1),(@view y[:,:,1]).-(@view ŷ[:,:,1]),title="error", c=:curl)
      p = plot((@view y[1,:,1]), label="y[1,:,1]")
      plot!(p, (@view ŷ[1,:,1]), label="ŷ[1,:,1]")
      fig = plot(hy,he,hŷ,p)
      display(fig)
    else
      p = plot((@view y[1,:,1]), label="y1")
      plot!(p, (@view ŷ[1,:,1]), label="ŷ1")
      for i in 2:size(y,1)
        plot!(p, (@view y[i,:,1]), label="y$(i)")
        plot!(p, (@view ŷ[i,:,1]), label="ŷ$(i)")
      end
      display(p)

      if size(y,3) > 1
        p = plot((@view y[1,:,2]), label="y1")
        plot!(p, (@view ŷ[1,:,2]), label="ŷ1")
        for i in 2:size(y,1)
          plot!(p, (@view y[i,:,2]), label="y$(i)")
          plot!(p, (@view ŷ[i,:,2]), label="ŷ$(i)")
        end
        display(p)
      end
    end

    # fig = plot(y[1,:,1], label="y1")
    # size(y,1) > 1 && plot!(fig, y[2,:,1], label="y2")
    # plot!(fig, ŷ[1,:,1], label="ŷ1")
    # size(y,1) > 1 && plot!(fig, ŷ[2,:,1], label="ŷ2")
    # display(fig)
    #
    # fig = plot(y[end-1,:,1], label="ye1")
    # size(y,1) > 1 && plot!(fig, y[end,:,1], label="ye2")
    # plot!(fig, ŷ[end-1,:,1], label="ŷe1")
    # size(y,1) > 1 && plot!(fig, ŷ[end,:,1], label="ŷe2")
    # display(fig)
  end
  return false
end


text:
module NetworkInference

using InformationMeasures
using Distributions
using Distributed
using Pkg
using DelimitedFiles
using SharedArrays

export
    # Common types
    Node,
    Edge,
    InferredNetwork,
    # Network inference algorithms
    AbstractNetworkInference,
    MINetworkInference,
    CLRNetworkInference,
    PUCNetworkInference,
    PIDCNetworkInference,
    # Functions for inferring networks
    get_nodes,
    write_network_file,
    read_network_file,
    get_adjacency_matrix,
    infer_network

include("common.jl")
include("network_inference.jl")
include("infer_network.jl")
include("empirical_bayes_glue.jl")

# Optional exports
if EB_EXISTS
export
    # Empirical Bayes glue functions
    to_index,
    make_priors,
    empirical_bayes
end

end # module