text:
## Group Members : Ashish Kumar, Pallavi Wats and Firat Melih Yilmaz

## Problem Set 2

using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables
using Distributions
using LaTeXStrings

## INTODUCTION TO NON-LINEAR OPTIMISATION IN JULIA

### Q1 ###
#optim is a function minimizer -- to find max of f(x) we need to minimize -f(x)
#optim packages gives us optimize() -- it requires three inputs : objective function, a string value and an optimizn algorithm

f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
negf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)   # random starting value
result = optimize(negf, startval, LBFGS())

## optimal value of x is -7.38 and the maximum is[-(-9.643*10^2) = 964.3]
m = -9.643134e+02

### Q2 ###
# use optim to compute OLS estimates of a simple linear regression using actual data
nlsw = CSV.read("/Users/prachi/Desktop/Econometrics6343/PS1/nlsw88.csv")
df = convert(DataFrame, nlsw)
save("/Users/prachi/Desktop/Econometrics6343/PS1/nlsw88.jld", "df", df)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)
#coefficients in order of intercept, age, white, collgrad

#alternatively,
using GLM
estimate = inv(X'*X)*X'*y
df.white = df.race.==1
estimate_lm = lm(@formula(married ~ age + white + collgrad), df)

#### Q3 ###
#use optim to estimate logit likelihood -- for this, I will need to pass Optim the negative of the likelihood function

function logit(alpha, X, d)
    a = log.(1 .+ exp.(X*alpha))
    b = (d.*(X*alpha))
    loglike =  sum(a .- b)
    return loglike
end
beta_hat_logit = optimize(b-> logit(b,X,y),rand(size(X,2)),LBFGS(),Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_logit.minimizer)

## minimum : 1.416951e+03 (1416.951)

### Q4 ###
# use glm(formula, data, family-Binomial, link-LogitLink) to check answer
α = glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
-loglikelihood(α) # answer = 1416.951

## Q5 ###
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)

#we also need to aggregate some of the occupation categories (category 10,11,12,13 -- 9)
df[df.occupation.==10,:occupation] .= 9
df[df.occupation.==11,:occupation] .= 9
df[df.occupation.==12,:occupation] .= 9
df[df.occupation.==13,:occupation] .= 9
freqtable(df, :occupation) # problem solved

#we changed the number of rows of df, we will need to re-define our X and y objects
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

#function mlogit(alpha, X, d)

    # your turn

    return loglike
#end


text:
# NOTE: PlanarAzimuthalSymmetry() does not include all possible wavenumbers
function wavenumbers_path(ω::T, medium::PhysicalMedium{Dim}, species::Species{T,Dim};
        symmetry::AbstractSymmetry{Dim} = PlanarAzimuthalSymmetry{Dim}(),
        tol::T = 1e-5,
        mesh_points::Int = 3,
        mesh_size::Number = one(T) / T(2),
        num_wavenumbers::Int = 3,
        verbose::Bool = false,
        inner_optimizer = NelderMead(; parameters = NelderMeadparameters()),
        optimoptions::Optim.Options{T} = Optim.Options(g_tol = tol^T(3),x_tol=tol^T(2)),
        k_effs::Vector{Complex{T}} = Complex{T}[],
        numberofparticles::T = Inf,
        kws...) where {T,Dim}

    # check parameters
        mesh_points = iseven(mesh_points) ? mesh_points + 1 : mesh_points

    # find at least one root to use as a scale for dk_x and dk_y
        eff_medium = effective_medium(medium, species; numberofparticles = numberofparticles)
        k0 = ω / eff_medium.c
        if isnan(k0) k0 = ω + T(0)*im end
        # abs(k0) can be used to non-dimensionlise k_vec
        # kscale = abs(k0)
        kscale = one(T)

        kφ = wavenumber_low_volumefraction(ω, medium, species;
            verbose = false, numberofparticles = numberofparticles
        )

        # guess initial mesh for lowest attenuating wavenumbers
        x_step = T(2) * max(abs(real(kφ)), real(k0), sqrt(eps(T)))
        ys = [imag(kφ), sqrt(eps(T))]

        k_dim_vecs = [[x,y] for x in LinRange(-x_step,x_step,2 * mesh_points + 3) for y in ys]
        # k_vecs is non-dimensional
        k_vecs = k_dim_vecs ./ kscale

        low_tol = min(1e-4, sqrt(tol))

    # Generate some asymptotic wavenumbers
        k_asyms = asymptotic_monopole_wavenumbers(ω, medium, species; num_wavenumbers = num_wavenumbers)
        k_asyms = k_asyms ./ kscale

    # Use asymptotic results to estimate maximum possible imaginary part.
        max_Imk = maximum(imag.(k_asyms))

    # The dispersion equation is given by: `dispersion([k1,k2]) = 0` where k_eff = k1 + im*k2.
        dispersion_dim = dispersion_equation(ω, medium, species, symmetry;
            tol = low_tol, numberofparticles = numberofparticles,
            kws...)
        dispersion(vec::Vector{T}) = dispersion_dim((vec[1] + vec[2]*im) .* kscale)

        k_vecs = [
            optimize(dispersion, kvec, #inner_optimizer,
            Optim.Options(x_tol=low_tol, g_tol = low_tol^3)
            ).minimizer
        for kvec in k_vecs]
        k_vecs = reduce_kvecs(k_vecs, low_tol/10)

        # add any specified keffs
        k_vecs = [
            [[real(kp),imag(kp)] for kp in k_asyms];
            k_vecs;
            [[real(kp),imag(kp)] ./ kscale for kp in k_effs]
        ]

        k_vecs = map(k_vecs) do k_vec
           res = optimize(dispersion, k_vec, inner_optimizer, optimoptions)
           if res.minimum < T(100)*tol || (Optim.converged(res) && res.minimum < T(10)*low_tol)
               res.minimizer
           else
               [zero(T),-one(T)]
           end
        end
        deleteat!(k_vecs, findall(v-> v == [zero(T),-one(T)], k_vecs) )
        k_vecs = reduce_kvecs(k_vecs, T(10)*tol)

        # Delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
        # deleteat!(k_vecs, findall([k_vec[2] < -sqrt(tol) for k_vec in k_vecs]))
        # delete wave travelling in wrong direction with small attenuation
        deleteat!(k_vecs, findall([-low_tol < abs(k_vec[2])/k_vec[1] < zero(T) for k_vec in k_vecs]))
        # deleteat!(k_vecs, find(k_vec[2] < tol && k_vec[1] < tol for k_vec in k_vecs))

    sort!(k_vecs, by= kv -> kv[2])

    # The best estimation of the mesh size is from asymptotic results

    dk = k_asyms[end] - k_asyms[end-1]

    dk_x = abs(real(dk)) * mesh_size
    dk_y = abs(imag(dk)) * mesh_size
    dk_xs = LinRange(-dk_x, dk_x, mesh_points)

    # Find two more roots, one with larger and smaller imaginary parts than the primary root kin
        two_roots = false # in case the roots are not on both branches of the tree, which is tested below
        kxs = dk_xs .+ k_vecs[1][1]
        ky =  k_vecs[1][2]
        invert = true

        # Find the first two roots that lead to the two branches of the root tree
        while !two_roots && (length(k_vecs) < num_wavenumbers) && ky <= max_Imk
            hits = [
                optimize(dispersion, [kx, ky], inner_optimizer, Optim.Options(x_tol = low_tol, g_tol = low_tol^3)).minimizer
            for kx in kxs]
            hits = reduce_kvecs(hits, low_tol)

            # Here we refine the hits
            hits = map(hits) do k_vec
                # res = optimize(dispersion, k_vec; g_tol = tol^2.0, f_tol = tol^4.0, x_tol=tol)
                res = optimize(dispersion, k_vec, inner_optimizer, optimoptions)
                if res.minimum < T(100)*tol || (Optim.converged(res) && res.minimum < T(10)*low_tol)
                    res.minimizer
                else
                    [zero(T),-one(T)]
                end
            end
            # deleteat!(hits, findall([-low_tol < abs(k_vec[2])/k_vec[1] < zero(T) for k_vec in hits]))
            deleteat!(hits, findall(v-> v == [zero(T),-one(T)], hits) )
            k_vecs = reduce_kvecs([hits;k_vecs], T(10)*tol)
            sort!(k_vecs, by= kv -> kv[2])

            if length(k_vecs) >= 2
                min1 = minimum(kv[1] for kv in k_vecs[2:end])
                max1 = maximum(kv[1] for kv in k_vecs[2:end])
                if min1*max1 < 0 || (max1-k_vecs[1][1])*(min1-k_vecs[1][1]) < 0
                    two_roots=true
                end
                dk_x += dk_x
                kxs = -k_vecs[2][1] .+ LinRange(-dk_x, dk_x, mesh_points)
                if invert
                    dk_y = k_vecs[2][2]/10
                    ky = k_vecs[2][2] - dk_y
                    invert = false
                elseif ky > 4*maximum(ks -> ks[2],k_vecs)
                    two_roots=true
                end
            elseif length(k_vecs) == 1
                dk_y += dk_y
            end
            ky = ky + dk_y
        end

    # Find more roots starting from the roots with the smaller imaginary part
        sort!(k_vecs, by = kv->kv[2])
        targets = k_vecs[2:end]
        dys = LinRange(0.5,1.0 + mesh_size, mesh_points)

        counter = 0
        while (num_wavenumbers > 2) && (counter < num_wavenumbers) && length(targets) > 0

            # procedure will be run at most once the desired number of wavenumbers
            counter = counter + 1

            # find roots following on from the smallest attenuating root
                sort!(targets, by = kv->kv[2])
                if verbose println("\n New target: $(targets[1])") end

            # find closest neighbour with smaller imaginary part
                i = findmin(
                    [(targets[1][2] != kvec[2]) ? norm(targets[1] - kvec) : Inf for kvec in k_vecs]
                )[2]
                # find at least one neighbour below, as the tree should grow
                j = (targets[1][2] > k_vecs[i][2]) ?
                    findmin(
                        [(targets[1][2] != kvec[2] && kvec[2] != k_vecs[i][2]) ? norm(targets[1] - kvec) : Inf for kvec in k_vecs]
                    )[2] :
                    findmin(
                        [(targets[1][2] > kvec[2] && kvec[2] != k_vecs[i][2]) ? norm(targets[1] - kvec) : Inf for kvec in k_vecs]
                    )[2]
                if j == Inf j = i end

                if verbose println("Closest neighbours: $(k_vecs[i]) and $(k_vecs[j])") end

            # create a mesh from this
                meshes = map([i,j]) do ind
                    tks = [k_vecs[ind] + dy.*(targets[1] - k_vecs[ind]) for dy in dys]
                    v = reverse(targets[1] - k_vecs[ind])
                    v[2] = -v[2]
                    hcat([
                        [tk + x.*v for x in LinRange(-mesh_size, mesh_size, mesh_points)]
                    for tk in tks]...)[:]
                end
                mesh = vcat(meshes...)
                deleteat!(mesh, findall([vec[2] < 0 || vec[2] > max_Imk for vec in mesh]))

            # search for roots from this mesh
                new_targets = map(mesh) do kin
                   optimize(dispersion, kin, inner_optimizer, Optim.Options(x_tol = low_tol, g_tol = low_tol^3)).minimizer
                end
                new_targets = reduce_kvecs(new_targets, low_tol/10)
                deleteat!(new_targets, findall(dispersion.(new_targets) .> low_tol))
                deleteat!(new_targets, findall([vec[2] < 0 || vec[2] > max_Imk for vec in new_targets]))

                # only keep targets which are not already in k_vecs
                new_targets = Vector{T}[
                        (findmin([norm(h - kvec) for kvec in k_vecs])[1] > low_tol) ? h : [zero(T),-one(T)]
                for h in new_targets]
                deleteat!(new_targets, findall(v-> v == [zero(T),-one(T)], new_targets))
                new_targets = reduce_kvecs(new_targets, low_tol)

                # Here we refine the new roots
                new_targets = map(new_targets) do k_vec
                    # res = optimize(dispersion, k_vec; g_tol = tol^2.0, f_tol = tol^4.0, x_tol=tol)
                    res = optimize(dispersion, k_vec, inner_optimizer, optimoptions)
                    if res.minimum < T(100)*tol || (Optim.converged(res) && res.minimum < T(10)*low_tol)
                        res.minimizer
                    else
                        [zero(T),-one(T)]
                    end
                end

                new_targets = [
                        (findmin([norm(h - kvec) for kvec in k_vecs])[1] > 10*tol) ? h : [zero(T),-one(T)]
                for h in new_targets]
                deleteat!(new_targets, findall(v-> v == [zero(T),-one(T)], new_targets))
                new_targets = reduce_kvecs(new_targets, T(10)*tol)

                if verbose println("New roots: $(new_targets)") end

                deleteat!(targets, 1)
                targets = Vector{T}[targets; new_targets]
                # group together wavenumbers which are closer than tol
                k_vecs = reduce_kvecs([new_targets; k_vecs], T(10)*tol)
                targets = filter(t -> t[2] <= max_Imk, targets)
        end

    # Finally delete unphysical waves, including waves travelling backwards with almost no attenuation. This only is important in the limit of very low frequency or very weak scatterers.
    # deleteat!(k_vecs, findall([k_vec[2] < -sqrt(tol) for k_vec in k_vecs]))
    # deleteat!(k_vecs, find(k_vec[2] < tol && k_vec[1] < tol for k_vec in k_vecs))
    deleteat!(k_vecs, findall([-low_tol < abs(k_vec[2])/k_vec[1] < zero(T) for k_vec in k_vecs]))

    # k_effs is dimensional
    k_effs = kscale .* [ k_vec[1] + k_vec[2]*im for k_vec in k_vecs]
    k_effs = sort(k_effs, by=imag)

    return k_effs
end

ext: jl

max_stars_repo_name: arturgower/EffectiveWaves


text:
# This file was generated, do not modify it. # hide
mutable struct KNNRidgeBlend <: DeterministicComposite
    knn_model::KNNRegressor
    ridge_model::RidgeRegressor
    knn_weight::Float64
end

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#pag 61(81)
'x'
Int('x')
Char(120)
Char(0x110000)
isvalid(Char, 0x110000)
'\u0'
'\u78'
'\u2200'
'\U10ffff'
str = "Hello, world.\n"
str[begin]
str[1]
str[end]

"""Contains "quote" characters"""
firstindex(str)
lastindex(str)
str[end-1]
str[end÷2]

str[begin-1]
str[end+1]
#ambas dão erro
str[4:9]
str[6]
str[6:6]

substr = SubString(str, 1, 5)
#retorna uma visualização da sttring pai
typeof(substr)

#a[i,j,...] é convertido pelo comp em getindex(a, i, j, ...)
#Unicode
s = "\u2200 x \u2203 y"
s[1]
s[2]
s[3]
s[4]
nextind(s, 1)
#proximo indice valido
nextind(s, nextind(s, 4))
s[end-1]
s[end-2]
#erro pois o anterior é multibyte
s[prevind(s, end, 2)]
prevind(s,lastindex(s), 2)

s[1:1]
s[1:2]
s[1:4]
length(s)
#7 caracteres e não 7 bytes
length(s) <= lastindex(s)
lastindex(s)

#ineficiente mas funciona
#onde não der erro tem caracter
for i = firstindex(s):lastindex(s)
    try
        println(s[i])
    catch
        #ignore
    end
end

#menos pior
for c in s  #string como obj iterável
    println(c)
end
#tambem pode usar nextind e prevind
#para obter índices válidos, ou...
collect(eachindex(s))
ncodeunits(s)
for i = 1:ncodeunits(s)
    println(codeunit(s, i))
end

s = "\xc0\xa0\xe2\x88\xe2|"
foreach(display, s)
#1 ou mais caracteres invalidos nao invalidam o todo
isvalid.(collect(s))
#só o ultimo pode ser impresso
s2 = "\xf7\xbf\xbf\xbf"
foreach(display, s2)
#'\U1fffff': Unicode U+1FFFFF (category In: Invalid, too high)

#Julia uses the UTF-8 encoding by default, and support for new encodings can
#be added by packages. For example,
#the LegacyStrings.jl package implements UTF16String and UTF32String types.

greet^3

#Concatenação
greet = "Hello"
whom = "world"
string(greet, ", ", whom, ".\n")
greet * ", " * whom * ".\n"
#na matemática, + significa uma operação comutativa
#onde a ordem dos operandos não importa.
# A + B == B + A
#Porém * é tipicamente uma operação ñ comutativa
#onde a ordem dos operandos importa
#Ex.: multiplicação de matrizes, A * B != B * A
#como aqui a ordem importa então se usa *,
#assim ficando de acordo com o uso na matematica
# * ou ⋅ forma um free monoid (S, *).
a, b = "\xe2\x88", "\x80"
c = a*b
collect.([a, b, c])
length.([a, b, c])

#Interpolação
# como em perl usa $ para interpolar
"$greet, $whom.\n"
#o sistema irá reescrever para string(greet, ", ", whom, ".\n")
"1 + 2 = $(1+2)"
v = [1, 2, 3]
#sem , é Array. Com , é vector
"v: $v"
print("I have \$100 in my account.\n")
str = """
Hello,
 world.
"""

"""\n This
is
  a test"""

"""hello"""
"""
Hello,
world."""

"a CRLF line ending\r\n"

#Operações comuns
"abracadabra" < "xylophone"
"abracadabra" == "xylophone"
"abracadabra" != "xylophone"
"1 + 2 = 3" == "1 + 2 = $(1+2)"

findfirst(isequal('o'), "xylophone")
findlast(isequal('o'), "xylophone")
findfirst(isequal('z'), "xylophone") #nada

findnext(isequal('o'), "xylophone", 1)
findnext(isequal('o'), "xylophone", 5)
findprev(isequal('o'), "xylophone", 5)
findnext(isequal('o'), "xylophone", 8)

occursin("world", "Hello, world.")
occursin("o", "Xylophon")
occursin('a', "Xylophon")

repeat(".:Z:.", 10)
join(["apples", "bananas", "pineapples"], ", ", " and ")

firstindex(str)
lastindex(str)
length(str)
length(str, 1, 5) # nº de char validos de i a j
ncodeunits(str) # nº de code units na string
codeunit(str,1)
thisind(str, 5)
nextind(str, 5, 1)
prevind(str, 5, 1)

#Regula Expressions
#Julia has Perl-compatible regular expressions (regexes),
#as provided by the PCRE library
r"^\s*(?:#|$)"
occursin(r"^\s*(?:#|$)", "not a comment")
occursin(r"^\s*(?:#|$)", "# a comment")
match(r"^\s*(?:#|$)", "not a comment")
match(r"^\s*(?:#|$)", "# a comment")
#similar a occursin, mas mostra o que combina nas strings
m = match(r"^\s*(?:#|$)", "# a comment")
if m === nothing
    println("not a comment")
else
    println("blank or comment")
end
#necessário para mostrar algo na tela.

m = match(r"^\s*(?:#\s*(.*?)\s*$|$)", "# a comment ")
m = match(r"[0-9]","aaaa1aaaa2aaaa3",1)
# 1 é o indice por onde começará a busca, pega a 1° após o indice
m.match
m.captures #substrings capturadas
m.offset #onde começa o que combinou
m.offsets #offset de substrings capturadas como vetor

m= match(r"(a|b)(c)?(d)", "acd")
m.match
m.captures
m.offset
m.offsets
first, second, third = m.captures; second
println(m)
#pag 76(96)

max_stars_repo_name: Mutanne/hiworld

max_issues_repo_name: Mutanne/hiworld

max_forks_repo_name: Mutanne/hiworld


text:
"""
Extended to deal with empty Dict.
"""
function _namedtuple(x::Dict)
    if x == Dict()
        return NamedTuple()  # empty NamedTuple
    else
        return namedtuple(x)
    end
end

"""
# Notes
- Currently, only iip (isinplace) method is supported.
- Default solver has been deprecated.
"""
function sim(state0, dyn, p=nothing;
        solver=nothing,  # DifferentialEquations.jl will find a default solver
        t0=0.0, tf=1.0,
        callback::DECallback=CallbackSet(),
        log_off=false,
        saveat=nothing,
        savestep=nothing,
        kwargs...,
    )
    if saveat == nothing && savestep == nothing
        saveat = []  # default
    elseif saveat != nothing && savestep == nothing
        # nothing
    elseif saveat == nothing && savestep != nothing
        saveat = t0:savestep:tf
    elseif saveat != nothing && savestep != nothing
        error("Assign values of either `saveat` or `savestep`")
    end
    tspan = (t0, tf)
    iip = true
    __dyn = (dx, x, p, t) -> dyn(dx, x, p, t)
    prob = ODEProblem{iip}(__dyn, state0, tspan, p)  # true: isinplace
    saved_values = SavedValues(Float64, Dict)
    cb_save = nothing
    if log_off == false
        # logging function
        if isinplace(prob)
            __log_indicator__ = __LOG_INDICATOR__()  # just an indicator for logging
            if hasmethod(dyn, Tuple{Any, Any, Any, Any, __LOG_INDICATOR__})
                log_func = function (x, t, integrator::DEIntegrator; kwargs...)
                    x = copy(x)  # `x` merely denotes a "view"
                    dyn(zero.(x), x, integrator.p, t, __log_indicator__; kwargs...)
                end
                cb_save = SavingCallback(log_func, saved_values;
                                         saveat=saveat, tdir=Int(sign(tspan[2]-tspan[1])))
            end
        else
            error("Not tested")
        end
        callback = CallbackSet(callback, cb_save)  # save values "after all other callbacks"
    end
    sol = solve(prob, solver;
                callback=callback,
                kwargs...)
    prob, sol
    if log_off == true
        return prob, sol
    else
        # recursive NamedTuple conversion from Dict; https://discourse.julialang.org/t/how-to-make-a-named-tuple-from-a-dictionary/10899/34?u=ihany
        recursive_namedtuple(x::Any) = x
        recursive_namedtuple(d::Dict) = _namedtuple(
                                                    Dict(k => recursive_namedtuple(v) for (k, v) in d)
                                                   )
        df = DataFrame(
                       time = saved_values.t,
                       sol = saved_values.saveval |> Map(recursive_namedtuple) |> collect,
                      )
        return prob, df
    end
end

"""
# Notes
- The basic concept is borrowed from [an MRAC example](https://jonniedie.github.io/ComponentArrays.jl/stable/examples/adaptive_control/).
- It is modified to be compatible with [SimulationLogger.jl](https://github.com/JinraeKim/SimulationLogger.jl).
# Limitations
- Conditional method definition is troublesome; see [#16](https://github.com/JinraeKim/FlightSims.jl/issues/16).
Instead of it, I decided to merely define a new method with argument `__log_indicator__::__LOG_INDICATOR__`,
which will provide "empty Dict" in the case of no logging.
"""
maybe_apply(f::Function, x, p, t) = f(x, p, t)
maybe_apply(f, x, p, t) = f
function apply_inputs(func; kwargs...)
    function simfunc(dx, x, p, t)
        func(dx, x, p, t; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...)
    end
    function simfunc(dx, x, p, t, __log_indicator__::__LOG_INDICATOR__)
        if hasmethod(func, Tuple{Any, Any, Any, Any, __LOG_INDICATOR__})
            return func(dx, x, p, t, __log_indicator__; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...)
        else
            return Dict()  # see the above method `NamedTupleTools.namedtuple`
        end
    end
    simfunc
end


text:
for N in [Float64, Rational{Int}, Float32]
    # random empty set
    rand(EmptySet)

    E = EmptySet{N}(2)
    B = BallInf(ones(N, 2), N(1))

    # testing that the empty set is an absorbing element for the cartesian product
    @test B * E isa EmptySet && E * B isa EmptySet
    # testing the mathematical alias ×
    @test B × E isa EmptySet && E × B isa EmptySet

    cpa = CartesianProductArray([B, N(2) * B, N(3) * B])
    @test cpa * E isa EmptySet && E * cpa isa EmptySet
    @test cpa × E isa EmptySet && E × cpa isa EmptySet

    # testing cp of empty set with itself
    @test E * E == E

    # testing that the empty set is an absorbing element for the Minkowski sum
    @test B + E isa EmptySet && E + B isa EmptySet
    # testing the mathematical alias ⊕
    @test B ⊕ E isa EmptySet && E ⊕ B isa EmptySet

    msa = MinkowskiSumArray([B, N(2) * B, N(3) * B])
    @test msa + E isa EmptySet && E + msa isa EmptySet
    @test msa ⊕ E isa EmptySet && E ⊕ msa isa EmptySet

    # testing M-sum of empty set with itself
    @test E + E == E

    # testing that the emptyset is neutral for the convex hull
    @test CH(B, E) == B
    @test CH(E, B) == B

    # test convex hull of empty set with itself
    @test CH(E, E) == E

    # dim
    @test dim(E) == 2

    # support vector
    @test_throws ErrorException σ(N[0], E)

    # boundedness
    @test isbounded(E) && isboundedtype(typeof(E))

    # isuniversal
    res, w = isuniversal(E, true)
    @test !isuniversal(E) && !res && w ∉ E

    # membership
    @test N[0] ∉ E
    @test N[0, 0] ∉ E

    # subset
    @test E ⊆ B && ⊆(E, B, true)[1]
    subset, point = ⊆(B, E, true)
    @test B ⊈ E && !subset && point ∈ B
    @test E ⊆ E && ⊆(E, E, true)[1]

    # emptiness check
    @test isempty(E)

    # an_element/norm/radius/diameter functions
    @test_throws ErrorException an_element(E)
    @test_throws ErrorException norm(E)
    @test_throws ErrorException radius(E)
    @test_throws ErrorException diameter(E)

    # vertices / vertices_list
    @test collect(vertices(E)) == vertices_list(E) == Vector{Vector{N}}()
    @test typeof(collect(vertices(E))) == typeof(vertices_list(E)) ==
          Vector{Vector{N}}

    # linear map of an empty set
    linear_map(ones(N, 2, 2), E) == E

    # translation
    @test translate(E, N[1, 2]) == E
    @test translate!(E, N[1, 2]) == E

    # disjointness
    for X in [B, Singleton(N[0, 0])]
        @test isdisjoint(E, X) && isdisjoint(X, E)
    end

    # projection
    @test project(EmptySet{N}(5), [1, 4, 5]) == EmptySet{N}(3)
end

# tests that only work with Float64 and Float32
for N in [Float64, Float32]
    B = Ball2(N[0, 0], N(1))
    E = EmptySet{N}(2)

    @test isdisjoint(E, B) && isdisjoint(B, E)
end

# default Float64 constructor
@test EmptySet(2) == ∅(2) == EmptySet{Float64}(2)

# intersection
for X in LazySets.subtypes(LazySet, true)
    if X <: RotatedHyperrectangle || isoperationtype(X)  # TODO #2391
        continue
    end
    Y = rand(X)
    E = EmptySet(dim(Y))
    @test intersection(Y, E) == intersection(E, Y) == E
end


text:
# Julia wrapper for header: /Users/dls/jl/openssl-1.0.1e/include/openssl/rand.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Int32 RAND_set_rand_method (Ptr{RAND_METHOD},) libcrypto
@c Ptr{RAND_METHOD} RAND_get_rand_method () libcrypto
@c Int32 RAND_set_rand_engine (Ptr{ENGINE},) libcrypto
@c Ptr{RAND_METHOD} RAND_SSLeay () libcrypto
@c None RAND_cleanup () libcrypto
@c Int32 RAND_bytes (Ptr{Uint8}, Int32) libcrypto
@c Int32 RAND_pseudo_bytes (Ptr{Uint8}, Int32) libcrypto
@c None RAND_seed (Ptr{None}, Int32) libcrypto
@c None RAND_add (Ptr{None}, Int32, Float64) libcrypto
@c Int32 RAND_load_file (Ptr{Uint8}, Int64) libcrypto
@c Int32 RAND_write_file (Ptr{Uint8},) libcrypto
@c Ptr{Uint8} RAND_file_name (Ptr{Uint8}, Csize_t) libcrypto
@c Int32 RAND_status () libcrypto
@c Int32 RAND_query_egd_bytes (Ptr{Uint8}, Ptr{Uint8}, Int32) libcrypto
@c Int32 RAND_egd (Ptr{Uint8},) libcrypto
@c Int32 RAND_egd_bytes (Ptr{Uint8}, Int32) libcrypto
@c Int32 RAND_poll () libcrypto
@c None ERR_load_RAND_strings () libcrypto


text:
"""
Non-negative CANDECOMP tensor decomposition.
"""
function nncp(tnsr::StridedArray,
              r::Integer;
              tol::Float64=1e-4,
              maxiter::Integer=100,
              compute_error::Bool=false,
              verbose::Bool=true)

    minimum(tnsr) >= 0 || error("Input tensor must be nonnegative.")
    num_modes = _check_tensor(tnsr, r)
    T_norm = vecnorm(tnsr)

    factors = [abs(F) * (T_norm ^ (1/num_modes) / vecnorm(F)) for F::Matrix{Float64} in _random_factors(size(tnsr), r)]
    factors_old = deepcopy(factors)
    factors_exp = deepcopy(factors)
    gram = Matrix{Float64}[F'F for F in factors]

    niters = 0
    converged = false
    obj_old = 0.0
    t = 1.0
    t_old = 1.0

    LB = ones(num_modes)
    LB_old = ones(num_modes)

    U = Array(Float64, num_modes, num_modes)
    M = Array(Float64, r, r)
    while !converged && niters < maxiter
        LB_old = copy(LB)
        for i in 1:num_modes
            idx = [num_modes:-1:i+1; i-1:-1:1]
            U = reduce(.*, gram[idx])
            LB[i] = vecnorm(U)
            M = _row_unfold(tnsr, i) * reduce(khatrirao, factors[idx])
            factors[i] = max(0, factors_exp[i] - (factors_exp[i] * U - M) * (1 / LB[i]))
            At_mul_B!(gram[i], factors[i], factors[i])
        end

        obj = sum(gram[num_modes] .* U) - 4 * sum(factors[num_modes] .* M)

        if obj > obj_old
            factors_exp = factors_old
        else
            t = 1 + sqrt(1 + 4 * t_old^2)
            weights = min((t_old - 2) / t, sqrt( LB_old ./ LB ))
            factors_exp = map((F, F_old, w) -> F + w * (F - F_old), factors, factors_old, weights)
            factors_old = deepcopy(factors)
            converged = obj_old - obj < tol * abs(obj_old)
            t_old = t
            obj_old = obj
        end
        niters += 1

    end

    verbose && _iter_status(converged, niters, maxiter)

    res = CANDECOMP((factors...), ones(r))
    if compute_error
      _set_rel_residue(res, tnsr)
    end
    return res
end


text:
"""
    OpConnectionIdentity

An iterable structure containing all the changes from one state to another, and
an associated value.
This effectively encodes a full row of an operator.
"""
struct OpConnection{A<:AbstractArray,B<:AbstractArray,C<:AbstractArray} <: AbsOpConnection
    mel::A
    changes::Vector{StateChanges{B,C}}

    length::Ref{Int}
end

function OpConnection(op::AbsLinearOperator)
    t = conn_type(op)
    v1, v2, v3 = t.parameters
    return OpConnection(v1(), Vector{StateChanges{v2,v3}}(), Ref(0))
end

function OpConnection{A,B,C}() where {A,B,C}
    return OpConnection(A(), Vector{typeof(StateChanges(B(),C()))}(), Ref(0))
end

function Base.resize!(c::OpConnection, n)
    # diminish size
    if n <= length(c)
        c.length[] = n
    else
        # increase size
        c.length[] = n
        resize!(c.mel,     n)
        resize!(c.changes, n)
    end

    return c
end

Base.eltype(::Type{OpConnection{A,B,C}}) where {A,B,C} = Tuple{eltype(A), StateChanges{B,C}}
Base.length(c::OpConnection) = c.length[]
Base.size(c::OpConnection) = (length(c), )
capacity(c::OpConnection) = length(c.mel)

Base.:(==)(a::OpConnection, b::OpConnection) = (a.mel == b.mel &&
                                              a.changes == b.changes)

@inline Base.push!(c::OpConnection, (m_els, to_change, new_values)) =
    push!(c, (m_els, StateChanges(to_change, new_values)))

function Base.push!(c::OpConnection, (mel, cngs)::Tuple{Number, StateChanges})
    if length(c) < capacity(c)
        resize!(c, length(c)+1)
        @inbounds setindex!(c, length(c), (mel, cngs))
    else
        c.length[] += 1
        push!(c.mel, mel)
        push!(c.changes, cngs)
    end
    return c
end

@inline function Base.setindex!(c::OpConnection, i, (mel, cngs)::Tuple{Number, StateChanges})
    @boundscheck checkbounds(1:length(c), i)
    @inbounds c.mel[i]     = mel
    @inbounds c.changes[i] = cngs
end

function Base.append!(c::OpConnection, (m_els, to_change, new_values)::Tuple{Vector, Vector, Vector})
    cngs = Vector{StateChanges{eltype(to_change), eltype(new_values)}}()
    for (tcng, nwvls) = zip(to_change, new_values)
        push!(cngs, StateChanges(tcng, nwvls))
    end
    return append!(c, (m_els, cngs))
end

function Base.append!(c::OpConnection, (m_els, cngs)::Tuple)
    l = length(c)
    resize!(c, l+length(m_els))

    for i=1:length(m_els)
        @inbounds setindex!(c, l+i, (m_els[i], cngs[i]))
    end
    return c
end

function Base.append!(c::OpConnection, new::OpConnection)
    L = length(new)

    L == 0 && return c

    # If this accumulator is empty, then simply append
    if length(c) == 0
        append!(c, (new.mel, new.changes))
    else
        # We assume that the first element is always for len(statechanges) = 0
        c.mel[1] += new.mel[1]

        L == 1 && return c

        append!(c, (uview(new.mel, 2:L), uview(new.changes, 2:L )))
    end

    return c
end

Base.@propagate_inbounds Base.getindex(c::OpConnection, i) = (c.mel[i], c.changes[i])

function Base.iterate(iter::OpConnection, state=1)
    state > length(iter) && return nothing
    @inbounds conn = iter[state]
    return (conn, state + 1)
end

function add!(c::OpConnection, mel, cngs::StateChanges)
    i = findfirst(isequal(cngs), c.changes)
    isnothing(i) ? push!(c, (mel, cngs)) : c.mel[i] += mel
    return c
end

function add!(c::OpConnection, o::OpConnection)
    for (mel, cngs)=o
        add!(c, mel, cngs)
    end
end

Base.conj!(op::OpConnection) = conj!(op.mel)

function clear_duplicates(c::OpConnection)
    u_cngs = unique(c.changes)
    u_mel  = similar(c.mel, 0)
    for el=u_cngs
        ids = findall(isequal(el), c.changes)
        push!(u_mel, sum(c.mel[ids]))
    end
    return OpConnection(u_mel, u_cngs)
end

#showing
function Base.show(io::IO, c::OpConnection{A,B,C}) where {A,B,C}
    print(io, "$(length(c))-elements - OpConnection{$(eltype(A)),StateChanges{$B,$C}}:\n")
    Base.print_matrix(IOContext(io, :compact=>true), collect(c))
    return io
end


text:
function _search(numbers::Vector, search::Number, start_idx::Int64, end_idx::Int64, lt)
    (start_idx < 1 || start_idx > end_idx) && return start_idx:end_idx
    middle_idx = div(start_idx + end_idx, 2)
    numbers[middle_idx] == search && return middle_idx:middle_idx
    lt(numbers[middle_idx], search) && return _search(numbers, search, middle_idx + 1, end_idx, lt)
    return _search(numbers, search, start_idx, middle_idx - 1, lt)
end

function binarysearch(numbers::Vector, search::Number; rev = false, lt = <, by = identity)
    rev && (lt = >)
    by != identity && (numbers = by.(numbers); search = by(search))
    _search(numbers, search, 1, length(numbers), lt)
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
type UniformMutation{R, T} <: Mutation
  stage::AbstractString
  representation::R
  min_value::T
  max_value::T
  rate::Float

  UniformMutation(stage::AbstractString, rep::R, rate::Float) =
    new(stage, rep, minimum_value(rep), maximum_value(rep), rate)
end

"""
Provides a definition for a uniform mutation operator.
"""
type UniformMutationDefinition <: MutationDefinition
  stage::AbstractString
  rate::Float
end

"""
Composes a uniform mutation operator from its definition.
"""
function compose!(d::UniformMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  UniformMutation(d.stage, d.rate, r)
end

"""
TODO: Describe uniform mutation.

**Parameters:**
* `stage::AbstractString`, the name of the developmental stage that this
  operator should be applied to. Defaults to the genotype if no stage is
  specified.
* `rate::Float`, the probability of a bit flip at any given index.
  Defaults to 0.01 if no rate is provided.
"""
uniform() = uniform(0.01)
uniform(rate::Float) = uniform("", rate)
uniform(stage::AbstractString) = uniform(stage, 0.01)
uniform(stage::AbstractString, rate::Float) =
  UniformMutationDefinition(stage, rate)

"""
Returns the number of inputs to this operation.
"""
num_inputs(o::UniformMutation) = 1

"""
Returns the number of outputs produced by this operation.
"""
num_outputs(o::UniformMutation) = 1

function mutate!{T}(o::UniformMutation, input::Vector{T})
  # If we knew the length of our genome in advance, how much faster
  # would this operation be?
  for i in 1:length(input)
    if rand() <= o.rate
      input[i] = rand(o.min_value:o.max_value)
    end
  end
  return input
end


text:
#
# simpleAcquisition.jl -
#
# This is the simplest one from the available examples.  It shows the user how
# to start the acquisition of multiples images, how to read the images and to
# detect if an overrun or even a timeout happened.
#
# This example also explains how to save acquired images.
#

module SimpleAcquisition

# In this example, we want to use the low-level interface as it is the closest
# to the C API.
import NuvuCameras: NC

function main()
    saveFormat = NC.FITS

    # Opens the acquisition channel using the automatic detection and 4 loop
    # buffers (recommended).
    myCam = NC.open(NC.Cam, NC.AUTO_UNIT, NC.AUTO_CHANNEL, 4)

    # For this example we will not bother checking for the available readout
    # modes on the camera and simply assume that the first one is valid.
    NC.setReadoutMode(myCam, 1)

    # Recover the readout time for use later.
    readoutTime = NC.getReadoutTime(myCam)

    # We set a non-zero exposure so that the first image is not blank.
    NC.setExposureTime(myCam, readoutTime)

    # Recover the exposure time for use later.
    exposureTime = NC.getExposureTime(myCam, true)

    # We set a reasonable waiting time.
    NC.setWaitingTime(myCam, 0.1*exposureTime)

    # Recover the waiting time for use later.
    waitingTime = NC.getWaitingTime(myCam, true)

    # Set a reasonable timeout on reading an image.  The delay between images
    # is the sum of the waiting time, the exposure and the readout time (if the
    # waiting time is non-zero).
    NC.setTimeout(myCam, waitingTime + exposureTime + readoutTime + 1000.0)

    # Opens the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Launches an acquisition on the frame grabber and requests images from the
    # camera (this function does not wait for the acquisition to be complete
    # before returning).
    NC.start(myCam, 0)

    # Reads a single received image: if a timeout occurs an error code will be
    # returned.
    myImage = NC.read(myCam)

    # Uses the "abort" function to tell the frame grabber to stop acquiring
    # frames from the camera, also tells the camera to stop sending images to
    # the frame grabber.
    NC.abort(myCam)

    # Closes the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    # Saves the image acquired.
    NC.saveImage(myCam, myImage, "FirstImage", saveFormat,
                 "This is my first image grabbed", true)

    println("You just acquired one image, by pressing Enter you will acquire 15 new images")
    readline(STDIN)

    # Opens the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Launches 15 acquisitions on the frame grabber and requests 15 images from
    # the camera (this function does not wait for the acquisition to be
    # complete before returning).
    NC.start(myCam, 15)

    # Loop in which acquired images are read.
    for i in 1:15
	# Reads the image received, if a timeout occurs an error code will be
	# returned.
	myImage = NC.read(myCam)

	# Checks if an overrun occured on the last image (implying that the
	# buffer we are reading has been overwritten prior to this "read"
	# call).
	overrun = NC.getOverrun(myCam)

	# If you wish to proceed with some processing on each image read, it
	# should be done here (in this example we will only save the image).

	# Saves each image acquired, at the end of each name the loop index
	# will be added.
	imageName = "Image_$i"
	NC.saveImage(myCam, myImage, imageName, saveFormat,
                     "This is my first series of images grabbed", true)
	println("Saved ", imageName)
    end

    # Closes the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    # Closes the acquisition channel no longer in use.
    NC.close(myCam)

    if is_windows()
	print("\nHit the Return key to close the command window.")
        flush(STDOUT)
        readline(STDIN)
    else
	println("The program ended.")
    end
end

end # module


text:
include("Yfind.jl")
include("is_it_taylor.jl")

H_test(u)=( (u[3]-u[1])^2+(u[2]-u[4])^2 )*( (u[3]+u[1])^2+(u[2]+u[4])^2 )/ ((u[2]^4+2*u[2]^2*(u[1]^2-1)+(1+u[1]^2)^2 )*(u[3]^4+2*u[3]^2*(u[4]^2+1)+(u[4]^2-1)^2 ))
ODE1(z,w)=conj(  im * w.*( 1 ./(z.^2-w.^2)+1 ./(1+w.^2) ))
ODE2(z,w)=conj(  im * z.*( 1 ./(w.^2-z.^2)+1 ./(1+z.^2) ))
function Eq_of_M!(t,u,du)
    du[1]=real(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[2]=imag(ODE1(u[1]+im*u[2],u[3]+im*u[4]))
    du[3]=real(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    du[4]=imag(ODE2(u[1]+im*u[2],u[3]+im*u[4]))
    return nothing
end
function g(t,u,du)
    px_ = constant_term(u[4])
 # if px > 0...
 if px_ > zero(px_)
     # ...return x
     return u[1]
 else
     #otherwise, discard the crossing
     return nothing
 end
end

function PSS_function_taylor(Q,P, Energy,  t_end)
    Q_0=0
    P_0=0

    Q_1=0
    P_1=0

    Q_2=0
    P_2=0

    Q_3=0
    P_3=0

    tol_dist=1e-2
    H=(2*Energy)^2

    Tabs_Q=zeros(0)
    Tabs_P=zeros(0)
    Y=Yfind(Q,P,H)
    if ~isempty(Y) && Y>0
        u0=zeros(4)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        # counter variable
        #are the initial conditions right?
        # prob = ODEProblem(Eq_of_M,u0,(0., t_end))
        tv_i, xv_i, tvS_i, xvS_i, gvS_i=taylorinteg(Eq_of_M!, g, u0, 0., t_end, 12, 1e-12,maxsteps=300000, newtoniter=30 )
        if is_it_taylor(xv_i,H,tol_dist)
            if isempty(xvS_i)
                Q_0=Q
                P_0=P
            else
                uf=zeros(5)
                uf[1]=xvS_i[end,1] #X
                uf[2]=xvS_i[end,2] #P
                uf[3]=xvS_i[end,3] #Q
                uf[4]=xvS_i[end,4] #Y
                dH=abs(H_test(uf)-H)
                if dH<1e-5
                    Tabs_Q=xvS_i[:,3]
                    Tabs_P=xvS_i[:,2]
                    if length(xvS_i[:1])<5
                        if length(Tabs_Q)==1
                            #never hits, exits
                            Q_0=Tabs_Q[end]
                            P_0=Tabs_P[end]
                        elseif length(Tabs_Q)==2
                            #hits once, exits
                            Q_0=Tabs_Q[end]
                            P_0=Tabs_P[end]
                            Q_1=Tabs_Q[end-1]
                            P_1=Tabs_P[end-1]
                        elseif length(Tabs_Q)==3
                            #hits twice, exits
                            Q_0=Tabs_Q[end]
                            P_0=Tabs_P[end]
                            Q_1=Tabs_Q[end-1]
                            P_1=Tabs_P[end-1]
                            Q_2=Tabs_Q[end-2]
                            P_2=Tabs_P[end-2]
                        elseif  length(Tabs_Q)==4
                            Q_0=Tabs_Q[end]
                            P_0=Tabs_P[end]
                            Q_1=Tabs_Q[end-1]
                            P_1=Tabs_P[end-1]
                            Q_2=Tabs_Q[end-2]
                            P_2=Tabs_P[end-2]
                            Q_3=Tabs_Q[end-3]
                            P_3=Tabs_P[end-3]
                        end
                    else
                    Q_0=Tabs_Q[end]
                    P_0=Tabs_P[end]
                    Q_1=Tabs_Q[end-1]
                    P_1=Tabs_P[end-1]
                    Q_2=Tabs_Q[end-2]
                    P_2=Tabs_P[end-2]
                    Q_3=Tabs_Q[end-3]
                    P_3=Tabs_P[end-3]
                    end
                end
            end
        end
    end
return     Q_0,P_0,Q_1,P_1,Q_2,P_2,Q_3,P_3
end


text:
test = """
    abc

    a
    b
    c

    ab
    ac

    a
    a
    a
    a

    b
    """

compute(fn, io) = sum(fn, split(read(io, String), "\n\n"))

function q1(group::AbstractString)::Integer
    return split(replace(group, r"\s+" => ""), "") |> unique |> length
end

function q2(group::AbstractString)::Integer
    people = split(group, r"\n"; keepempty = false)
    return split.(people, "") |> Base.splat(intersect) |> length
end

@assert compute(q1, IOBuffer(test)) == 11
@assert compute(q2, IOBuffer(test)) == 6

println("Q1: ", compute(q1, "day06.in"))
println("Q2: ", compute(q2, "day06.in"))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
"""
    setindex!(x, value, i, j)

Assign the `SummarizedExperiment` `value` to a subset of `SummarizedExperiment` `x` by the rows or columns based on `i` and `j`, respectively.
Types for the arguments to `i` and `j` are similar to those for arrays:

- An integer `Vector` containing indices.
- An `Int` containing a single index.
- A boolean `Vector` of length equal to the relevant dimension,
  indicating whether each entry of that dimension should be retained.
- A `:` operator to retain the entirety of a dimension's extent.

On assignment, the assay values in the specified subset of `x` will be replaced by the corresponding values in `value`.
Rows of the `rowdata(x)` and `coldata(x)` will be replaced by those in `value`, according to `i` and `j` respectively.
Metadata fields in `metadata(value)` will be added to or overwrite those in `metadata(x)`.

It is assumed that `x` and `value` contain the same name, order and type of columns in their `rowdata` and `coldata`.
Similarly, both objects should contain the same name, order and type of arrays in their `assays`.

# Examples
```jldoctest
julia> using SummarizedExperiments

julia> x = exampleobject(20, 10);

julia> x[:,1] = x[:,2];

julia> sn = coldata(x)[!,"name"];

julia> sn[1] == sn[2]
true

julia> y = assay(x);

julia> y[:,1] == y[:,2]
true
```
"""
function Base.setindex!(x::SummarizedExperiment, value::SummarizedExperiment, i, j)
    indices = to_indices(x, (i, j))

    if isa(indices[1], Int)
        index_rows = [indices[1]]
    else
        index_rows = indices[1]
    end

    if isa(indices[2], Int)
        index_cols = [indices[2]]
    else
        index_cols = indices[2]
    end

    # Make sure we don't modify 'x' by reference until the very end.
    rd_copy = copy(x.rowdata)
    rd_copy[index_rows,:] = value.rowdata
    cd_copy = copy(x.coldata)
    cd_copy[index_cols,:] = value.coldata

    if collect(keys(x.assays)) != collect(keys(value.assays))
        throw(ErrorException("'x' and 'value' should have the same assay names"))
    end

    output = deepcopy(x.assays) # no way around this, as any additions in progress
                                # will operate by reference and overwrite the source.

    for (k, v) in output
        subdex = Array{Any,1}(undef, length(size(v)))
        fill!(subdex, :)
        subdex[1] = index_rows
        subdex[2] = index_cols
        setindex!(v, value.assays[k], subdex...)
    end

    # Setting everything to their modified states. Everything at this point
    # is assumed to be a no-throw, so we are guaranteed to return.
    x.rowdata = rd_copy
    x.coldata = cd_copy
    x.assays = output

    for (m, val) in value.metadata
        x.metadata[m] = val
    end

    return x
end


text:
using Test, ArgParse
using SignalingDimension

try print_test_results catch err
    include("../utilities.jl")
    using .ScriptUtilities: print_test_results
end
# parsing command line arguments
arg_table = ArgParseSettings()
@add_arg_table! arg_table begin
    "-Y"
        help = "max number of outputs to consider"
        arg_type = Int64
        default = 24
    "--dir"
        help = "directory to print results .txt to"
        arg_type = String
        default = "./script/facet_verifications/results/"
    "--stdout"
        help = "Print output to STDOUT. If not flagged, output is printed to `.txt` file."
        action = :store_true
end
args = parse_args(ARGS, arg_table)

# verification script: verifies all k-guessing facets with Y ∈ [4, Y_max] and k = 2
function verify_k_guessing_facet_k2(Y_max :: Int64)
    @testset "Testing all k-guessing facets with `k = 2` and `Y ≤ $Y_max`." begin
        @testset "Verifying Y = $Y, d = $(Y - 2), and k = 2" for Y in 4:Y_max
            k = 2
            d = Y - k
            vertices = aff_ind_k_guessing_vertices(Y, d, k)
            facet = k_guessing_facet(Y, d, k)

            @test verify_facet(d, facet, vertices)
        end
    end
end

# outputting result to file or stdout
print_test_results(
    verify_k_guessing_facet_k2,
    params=[args["Y"]],
    stdout=args["stdout"],
    dir=args["dir"]
)


text:
const signal_by_id = Dict{UInt64, WeakRef}()

export Signal

type Signal
    handle :: hsa_signal_t

    consumers :: Set{Agent}

    is_alive :: Bool

    function Signal(h::hsa_signal_t, consumers::Set{Agent} = Set{Agent}())
        assert_runtime_alive()

        if haskey(signal_by_id, h.handle)
            existing = signal_by_id[h.handle].value
            if existing != nothing
                return existing
            end
        end

        s = new(h, consumers, true)

        signal_by_id[h.handle] = WeakRef(s)

        return s
    end
end

# Wait Expectancy Hints
const WaitActive = HSA_WAIT_STATE_ACTIVE
const WaitBlocked = HSA_WAIT_STATE_BLOCKED

function Signal(; value::hsa_signal_value_t = 0, consumers = Set{Agent}())
    h = signal_create(value, consumers)
    c = Set{Agent}(consumers)

    s = Signal(h, c)

    finalizer(s, finalize_signal)

    return s
end

function finalize_signal(s :: Signal)
    if s.is_alive
        hsa_signal_destroy(s)
        s.is_alive = false
    end
end

import Base.convert

convert(::Type{hsa_signal_t}, s::UInt64) = hsa_signal_t(s)
convert(::Type{UInt64}, s::hsa_signal_t) = s.handle

convert(::Type{hsa_signal_t}, s::Signal) = s.handle
convert(::Type{UInt64}, s::Signal) = s.handle.handle

convert(::Type{Signal}, h :: UInt64) = Signal(hsa_signal_t(h))
convert(::Type{Signal}, h :: hsa_signal_t) = Signal(h)

import Base.get

get(s::Signal) = load(s)

function signal_create(initial_value::hsa_signal_value_t, consumers)
    local num_consumers, consumer_arr, consumer_ptr
    if !isempty(consumers)
        num_consumers = size(consumers,1)
        consumer_arr = hsa_agent_t[ a.handle for a in consumers ]
        consumer_ptr = pointer(consumer_arr)
    else
        num_consumers = 0x0
        consumer_ptr = C_NULL
    end

    res = Ref(hsa_signal_t(0))

    err = ccall((:hsa_signal_create, libhsa), hsa_status_t, (hsa_signal_value_t, UInt32, Ptr{hsa_agent_t}, Ptr{hsa_signal_t}),
                initial_value, num_consumers, consumer_ptr, res)

    test_status(err)

    return res.x
end

import Base.wait

function wait(s::Signal, cond::Symbol, compare_value::hsa_signal_value_t;
    timeout_hint = typemax(UInt64),
    wait_state_hint::hsa_wait_state_t = WaitBlocked
    )
    const cond_map = Dict(
        :(==) => HSA_SIGNAL_CONDITION_EQ,
        :(!=) => HSA_SIGNAL_CONDITION_NE,
        :(<) => HSA_SIGNAL_CONDITION_LT,
        :(>=) => HSA_SIGNAL_CONDITION_GTE
        )

    if !haskey(cond_map, cond)
        error("unsupported signal wait condition $cond")
    end

    hsa_cond = cond_map[cond]

    timeout_hint = convert(UInt64, timeout_hint)

    return wait(s, hsa_cond, compare_value, timeout_hint, wait_state_hint)
end


text:
using Bridge
using Random
using DataFrames
using CSV

SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
include(joinpath(SRC_DIR, "auxiliary", "data_simulation_fns.jl"))
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "fitzHughNagumo.jl"))
param = :complexConjug
FILENAME_OUT = joinpath(OUT_DIR,
                     "test_path_part_obs_"*String(param)*".csv")

K = 4

P = FitzhughDiffusion(param, 10.0, -8.0, 15.0, 0.0, 3.0)
# starting point under :regular parametrisation
# translate to conjugate parametrisation
x0 = [regularToConjug(ℝ{2}(-rand(), rand()), P.ϵ, 0.0) for k in 1:K]

dt = 1/50000
T = 10.0
tt = 0.0:dt:T
num_obs = 100
skip_ = div(length(tt), num_obs)

Random.seed!(4)
XX = [simulateSegment(0.0, x0[k], P, tt)[1][1:skip_:length(tt)-1] for k in 1:K]

#df = DataFrame(time=XX[1].tt, x1=[x[1] for x in XX[1].yy],
#            x2=[(i==1 ? x0[2] : NaN) for (i,t) in enumerate(Time)])
#CSV.write(FILENAME_OUT, df)


text:
###### Weight vector #####

immutable WeightVec{W,Vec<:RealVector}
    values::Vec
    sum::W
end

"""
    WeightVec(vs, [wsum])

Construct a `WeightVec` with weight values `vs` and sum of weights `wsum`.
If omitted, `wsum` is computed.
"""
WeightVec{Vec<:RealVector,W<:Real}(vs::Vec,wsum::W) = WeightVec{W,Vec}(vs, wsum)
WeightVec(vs::RealVector) = WeightVec(vs, sum(vs))

"""
    weights(vs)

Construct a `WeightVec` from a given array.
"""
weights(vs::RealVector) = WeightVec(vs)
weights(vs::RealArray) = WeightVec(vec(vs))

eltype(wv::WeightVec) = eltype(wv.values)
length(wv::WeightVec) = length(wv.values)
values(wv::WeightVec) = wv.values
sum(wv::WeightVec) = wv.sum
isempty(wv::WeightVec) = isempty(wv.values)

Base.getindex(wv::WeightVec, i) = getindex(wv.values, i)

##### Weighted sum #####

## weighted sum over vectors

"""
    wsum(v, w::AbstractVector, [dim])

Compute the weighted sum of an array `v` with weights `w`, optionally over the dimension `dim`.
"""
wsum(v::AbstractVector, w::AbstractVector) = dot(v, w)
wsum(v::AbstractArray, w::AbstractVector) = dot(vec(v), w)

# Note: the methods for BitArray and SparseMatrixCSC are to avoid ambiguities
Base.sum(v::BitArray, w::WeightVec) = wsum(v, values(w))
Base.sum(v::SparseMatrixCSC, w::WeightVec) = wsum(v, values(w))
Base.sum(v::AbstractArray, w::WeightVec) = dot(v, values(w))

## wsum along dimension
#
#  Brief explanation of the algorithm:
#  ------------------------------------
#
#  1. _wsum! provides the core implementation, which assumes that
#     the dimensions of all input arguments are consistent, and no
#     dimension checking is performed therein.
#
#     wsum and wsum! perform argument checking and call _wsum!
#     internally.
#
#  2. _wsum! adopt a Cartesian based implementation for general
#     sub types of AbstractArray. Particularly, a faster routine
#     that keeps a local accumulator will be used when dim = 1.
#
#     The internal function that implements this is _wsum_general!
#
#  3. _wsum! is specialized for following cases:
#     (a) A is a vector: we invoke the vector version wsum above.
#         The internal function that implements this is _wsum1!
#
#     (b) A is a dense matrix with eltype <: BlasReal: we call gemv!
#         The internal function that implements this is _wsum2_blas!
#
#     (c) A is a contiguous array with eltype <: BlasReal:
#         dim == 1: treat A like a matrix of size (d1, d2 x ... x dN)
#         dim == N: treat A like a matrix of size (d1 x ... x d(N-1), dN)
#         otherwise: decompose A into multiple pages, and apply _wsum2!
#         for each
#
#     (d) A is a general dense array with eltype <: BlasReal:
#         dim <= 2: delegate to (a) and (b)
#         otherwise, decompose A into multiple pages
#

function _wsum1!(R::AbstractArray, A::AbstractVector, w::AbstractVector, init::Bool)
    r = wsum(A, w)
    if init
        R[1] = r
    else
        R[1] += r
    end
    return R
end

function _wsum2_blas!{T<:BlasReal}(R::StridedVector{T}, A::StridedMatrix{T}, w::StridedVector{T}, dim::Int, init::Bool)
    beta = ifelse(init, zero(T), one(T))
    trans = dim == 1 ? 'T' : 'N'
    BLAS.gemv!(trans, one(T), A, w, beta, R)
    return R
end

function _wsumN!{T<:BlasReal,N}(R::StridedArray{T}, A::StridedArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool)
    if dim == 1
        m = size(A, 1)
        n = div(length(A), m)
        _wsum2_blas!(view(R,:), reshape(A, (m, n)), w, 1, init)
    elseif dim == N
        n = size(A, N)
        m = div(length(A), n)
        _wsum2_blas!(view(R,:), reshape(A, (m, n)), w, 2, init)
    else # 1 < dim < N
        m = 1
        for i = 1:dim-1; m *= size(A, i); end
        n = size(A, dim)
        k = 1
        for i = dim+1:N; k *= size(A, i); end
        Av = reshape(A, (m, n, k))
        Rv = reshape(R, (m, k))
        for i = 1:k
            _wsum2_blas!(view(Rv,:,i), view(Av,:,:,i), w, 2, init)
        end
    end
    return R
end

function _wsumN!{T<:BlasReal,N}(R::StridedArray{T}, A::DenseArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool)
    @assert N >= 3
    if dim <= 2
        m = size(A, 1)
        n = size(A, 2)
        npages = 1
        for i = 3:N
            npages *= size(A, i)
        end
        rlen = ifelse(dim == 1, n, m)
        Rv = reshape(R, (rlen, npages))
        for i = 1:npages
            _wsum2_blas!(view(Rv,:,i), view(A,:,:,i), w, dim, init)
        end
    else
        _wsum_general!(R, @functorize(identity), A, w, dim, init)
    end
    return R
end

# General Cartesian-based weighted sum across dimensions
@ngenerate N typeof(R) function _wsum_general!{T,RT,WT,N}(R::AbstractArray{RT}, f::supertype(typeof(@functorize(abs))),
                                                          A::AbstractArray{T,N}, w::AbstractVector{WT}, dim::Int, init::Bool)
    init && fill!(R, zero(RT))
    wi = zero(WT)
    if dim == 1
        @nextract N sizeR d->size(R,d)
        sizA1 = size(A, 1)
        @nloops N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
            @inbounds r = (@nref N R j)
            for i_1 = 1:sizA1
                @inbounds r += f(@nref N A i) * w[i_1]
            end
            @inbounds (@nref N R j) = r
        end
    else
        @nloops N i A d->(if d == dim
                               wi = w[i_d]
                               j_d = 1
                           else
                               j_d = i_d
                           end) @inbounds (@nref N R j) += f(@nref N A i) * wi
    end
    return R
end

@ngenerate N typeof(R) function _wsum_centralize!{T,RT,WT,N}(R::AbstractArray{RT}, f::supertype(typeof(@functorize(abs))),
                                                             A::AbstractArray{T,N}, w::AbstractVector{WT}, means,
                                                             dim::Int, init::Bool)
    init && fill!(R, zero(RT))
    wi = zero(WT)
    if dim == 1
        @nextract N sizeR d->size(R,d)
        sizA1 = size(A, 1)
        @nloops N i d->(d>1? (1:size(A,d)) : (1:1)) d->(j_d = sizeR_d==1 ? 1 : i_d) begin
            @inbounds r = (@nref N R j)
            @inbounds m = (@nref N means j)
            for i_1 = 1:sizA1
                @inbounds r += f((@nref N A i) - m) * w[i_1]
            end
            @inbounds (@nref N R j) = r
        end
    else
        @nloops N i A d->(if d == dim
                               wi = w[i_d]
                               j_d = 1
                           else
                               j_d = i_d
                           end) @inbounds (@nref N R j) += f((@nref N A i) - (@nref N means j)) * wi
    end
    return R
end

# N = 1
_wsum!{T<:BlasReal}(R::StridedArray{T}, A::DenseArray{T,1}, w::StridedVector{T}, dim::Int, init::Bool) =
    _wsum1!(R, A, w, init)

# N = 2
_wsum!{T<:BlasReal}(R::StridedArray{T}, A::DenseArray{T,2}, w::StridedVector{T}, dim::Int, init::Bool) =
    (_wsum2_blas!(view(R,:), A, w, dim, init); R)

# N >= 3
_wsum!{T<:BlasReal,N}(R::StridedArray{T}, A::DenseArray{T,N}, w::StridedVector{T}, dim::Int, init::Bool) =
    _wsumN!(R, A, w, dim, init)

_wsum!(R::AbstractArray, A::AbstractArray, w::AbstractVector, dim::Int, init::Bool) = _wsum_general!(R, @functorize(identity), A, w, dim, init)

## wsum! and wsum

wsumtype{T,W}(::Type{T}, ::Type{W}) = typeof(zero(T) * zero(W) + zero(T) * zero(W))
wsumtype{T<:BlasReal}(::Type{T}, ::Type{T}) = T

"""
    wsum!(R, A, w, dim; init=true)

Compute the weighted sum of `A` with weights `w` over the dimension `dim` and store
the result in `R`. If `init=false`, the sum is added to `R` rather than starting
from zero.
"""
function wsum!{T,N}(R::AbstractArray, A::AbstractArray{T,N}, w::AbstractVector, dim::Int; init::Bool=true)
    1 <= dim <= N || error("dim should be within [1, $N]")
    ndims(R) <= N || error("ndims(R) should not exceed $N")
    length(w) == size(A,dim) || throw(DimensionMismatch("Inconsistent array dimension."))
    # TODO: more careful examination of R's size
    _wsum!(R, A, w, dim, init)
end

function wsum{T<:Number,W<:Real}(A::AbstractArray{T}, w::AbstractVector{W}, dim::Int)
    length(w) == size(A,dim) || throw(DimensionMismatch("Inconsistent array dimension."))
    _wsum!(Array(wsumtype(T,W), Base.reduced_dims(size(A), dim)), A, w, dim, true)
end

# extended sum! and wsum

Base.sum!{W<:Real}(R::AbstractArray, A::AbstractArray, w::WeightVec{W}, dim::Int; init::Bool=true) =
    wsum!(R, A, values(w), dim; init=init)

Base.sum{T<:Number,W<:Real}(A::AbstractArray{T}, w::WeightVec{W}, dim::Int) = wsum(A, values(w), dim)

###### Weighted means #####

"""
    wmean(v, w::AbstractVector)

Compute the weighted mean of an array `v` with weights `w`.
"""
function wmean{T<:Number}(v::AbstractArray{T}, w::AbstractVector)
    Base.depwarn("wmean is deprecated, use mean(v, weights(w)) instead.", :wmean)
    mean(v, weights(w))
end

Base.mean(v::AbstractArray, w::WeightVec) = sum(v, w) / sum(w)

Base.mean!(R::AbstractArray, A::AbstractArray, w::WeightVec, dim::Int) =
    scale!(Base.sum!(R, A, w, dim), inv(sum(w)))

wmeantype{T,W}(::Type{T}, ::Type{W}) = typeof((zero(T)*zero(W) + zero(T)*zero(W)) / one(W))
wmeantype{T<:BlasReal}(::Type{T}, ::Type{T}) = T

Base.mean{T<:Number,W<:Real}(A::AbstractArray{T}, w::WeightVec{W}, dim::Int) =
    mean!(Array(wmeantype(T, W), Base.reduced_dims(size(A), dim)), A, w, dim)

###### Weighted median #####

function Base.median{W<:Real}(v::RealVector, w::WeightVec{W})
    isempty(v) && error("median of an empty array is undefined")
    if length(v) != length(w)
        error("data and weight vectors must be the same size")
    end
    @inbounds for x in w.values
        isnan(x) && error("weight vector cannot contain NaN entries")
    end
    @inbounds for x in v
        isnan(x) && return x
    end
    mask = w.values .!= 0
    if !any(mask)
        error("all weights are zero")
    end
    if all(w.values .<= 0)
        error("no positive weights found")
    end
    v = v[mask]
    wt = w[mask]
    midpoint = w.sum / 2
    maxval, maxind = findmax(wt)
    if maxval > midpoint
        v[maxind]
    else
        permute = sortperm(v)
        cumulative_weight = zero(eltype(wt))
        i = 0
        for (i, p) in enumerate(permute)
            if cumulative_weight == midpoint
                i += 1
                break
            elseif cumulative_weight > midpoint
                cumulative_weight -= wt[p]
                break
            end
            cumulative_weight += wt[p]
        end
        if cumulative_weight == midpoint
            middle(v[permute[i-2]], v[permute[i-1]])
        else
            middle(v[permute[i-1]])
        end
    end
end

"""
    wmedian(v, w)

Compute the weighted median of an array `v` with weights `w`, given as either a
vector or `WeightVec`.
"""
wmedian(v::RealVector, w::RealVector) = median(v, weights(w))
wmedian{W<:Real}(v::RealVector, w::WeightVec{W}) = median(v, w)

###### Weighted quantile #####

# http://stats.stackexchange.com/questions/13169/defining-quantiles-over-a-weighted-sample
# In the non weighted version, we compute a vector of index h(N, p)
# and take interpolation between floor and ceil of this index
# Here there is a supplementary function from index to weighted index k -> Sk

"""
    quantile(v, w::WeightVec, p)

Compute `p`th quantile(s) of `v` with weights `w`.
"""
function quantile{V, W <: Real}(v::RealVector{V}, w::WeightVec{W}, p::RealVector)

    # checks
    isempty(v) && error("quantile of an empty array is undefined")
    isempty(p) && throw(ArgumentError("empty quantile array"))

    w.sum == 0 && error("weight vector cannot sum to zero")
    length(v) == length(w) || error("data and weight vectors must be the same size, got $(length(v)) and $(length(w))")
    for x in w.values
        isnan(x) && error("weight vector cannot contain NaN entries")
        x < 0 && error("weight vector cannot contain negative entries")
    end

    # full sort
    vw = sort!(collect(zip(v, w.values)))

    wsum = w.sum

    # prepare percentiles
    ppermute = sortperm(p)
    p = p[ppermute]
    p = bound_quantiles(p)

    # prepare out vector
    N = length(vw)
    out = Array(typeof(zero(V)/1), length(p))
    fill!(out, vw[end][1])

    # start looping on quantiles
    cumulative_weight, Sk, Skold =  zero(W), zero(W), zero(W)
    vk, vkold = zero(V), zero(V)
    k = 1
    for i in 1:length(p)
        h = p[i] * (N - 1) * wsum
        if h == 0
            # happens when N or p or wsum equal zero
            out[ppermute[i]] = vw[1][1]
        else
            while Sk <= h
                # happens in particular when k == 1
                vk, wk = vw[k]
                cumulative_weight += wk
                if k >= N
                    # out was initialized with maximum v
                    return out
                end
                k += 1
                Skold, vkold = Sk, vk
                vk, wk = vw[k]
                Sk = (k - 1) * wk + (N - 1) * cumulative_weight
            end
            # in particular, Sk is different from Skold
            g = (h - Skold) / (Sk - Skold)
            out[ppermute[i]] = vkold + g * (vk - vkold)
        end
    end
    return out
end

# similarly to statistics.jl in Base
function bound_quantiles{T <: Real}(qs::AbstractVector{T})
    epsilon = 100 * eps()
    if (any(qs .< -epsilon) || any(qs .> 1+epsilon))
        throw(ArgumentError("quantiles out of [0,1] range"))
    end
    T[min(one(T), max(zero(T), q)) for q = qs]
end

quantile{W <: Real}(v::RealVector, w::WeightVec{W}, p::Number) = quantile(v, w, [p])[1]

"""
    wquantile(v, w, p)

Compute the `p`th quantile(s) of `v` with weights `w`, given as either a vector
or a `WeightVec`.
"""
wquantile{W <: Real}(v::RealVector, w::WeightVec{W}, p::RealVector) = quantile(v, w, p)
wquantile{W <: Real}(v::RealVector, w::WeightVec{W}, p::Number) = quantile(v, w, [p])[1]
wquantile(v::RealVector, w::RealVector, p::RealVector) = quantile(v, weights(w), p)
wquantile(v::RealVector, w::RealVector, p::Number) = quantile(v, weights(w), [p])[1]


text:
#
println("Ca) Test of the Einstein module with ASF from an internally generated multiplet.")
#
setDefaults("print summary: open", "zzz-Einstein.sum")
grid=JAC.Radial.Grid(true)
## grid = Radial.Grid(Radial.Grid(true), rnt = 2.0e-5,h = 1.0e-2, hp = 1.0e-2, NoPoints = 2000)
wa = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(26.),
                        properties=[JAC.EinsteinX()],
                        configs=[Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d"), Configuration("[Ne] 3s^2 3p^5")],
                        einsteinSettings=Einstein.Settings([E1], true, LineSelection(false), 0., 0., 10000. ) )

wb = perform(wa; output=true)
setDefaults("print summary: close", "")


text:
using ICU
using Base.Test

# Tests for not overrunning buffer
str = "\u3b0"
upp = "\u3a5\u308\u301"

@test u_strToUpper(utf16(str^8)) == utf16(upp^8)


text:
# BLUP test
using ProgenyTestingTools
using DelimitedFiles
using SparseArrays

# mu, var_p, h2_poly, h2_qtl, rep
par = PTParameters(50, 100, 0.36, 0, 0, 1.0)
hp = generate_population(par,nm=20,nf=20)
pop = generate_population(par)
bulls1 = migrate_from_hp!(hp,pop,random_sampling(hp,2,male=true))
cows1 = migrate_from_hp!(hp,pop,random_sampling(hp,3,female=true))
mgroup1 = generate_group(pop, sires=bulls1)
fgroup1 = generate_group(pop, dams=cows1)
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=1, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 1]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=2, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 2]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=3, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 3]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=4, repeated=true)
calving!(fgroup1)

write_files_for_blup(pop,"phenotype.txt","pedigree.txt", repeated=false)

(lhs,rhs) = ProgenyTestingTools.build_mme(pop,"phenotype.txt","pedigree.txt",false,false)
sol = lhs \ rhs

open("solutions.txt", "w") do io
   writedlm(io, sol)
end

# repeatability model
# mu, var_p, h2_poly, h2_qtl, rep
par = PTParameters(50, 100, 0.36, 0, 0.45, 1.0)
hp = generate_population(par,nm=20,nf=20)
pop = generate_population(par)
bulls1 = migrate_from_hp!(hp,pop,random_sampling(hp,2,male=true))
cows1 = migrate_from_hp!(hp,pop,random_sampling(hp,3,female=true))
mgroup1 = generate_group(pop, sires=bulls1)
fgroup1 = generate_group(pop, dams=cows1)
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=1, repeated=true)
assign_phenotype!(fgroup1, gen=1, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 1]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=2, repeated=true)
assign_phenotype!(fgroup1, gen=2, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 2]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=3, repeated=true)
assign_phenotype!(fgroup1, gen=3, repeated=true)
calving!(fgroup1)

fgroup1.dams = fgroup1.id[fgroup1.generation .== 3]
mating!(mgroup1, fgroup1, "fgroup", method="dairy_standard_ai")
assign_phenotype!(fgroup1, gen=4, repeated=true)
assign_phenotype!(fgroup1, gen=4, repeated=true)
calving!(fgroup1)

write_files_for_blup(pop,"phenotype2.txt","pedigree2.txt", repeated=true)

(lhs,rhs) = ProgenyTestingTools.build_mme(pop,"phenotype2.txt","pedigree2.txt",true,false, verbose=true)
sol = lhs \ rhs

open("solutions2.txt", "w") do io
   writedlm(io, sol)
end


text:
@reexport module BaselineAlgos

using TSML.TSMLTypes
using TSML.Utils

import TSML.TSMLTypes.fit!
import TSML.TSMLTypes.transform!

export Baseline,Identity
       fit!,
       transform!

using StatsBase: mode

"""
    Baseline(
       default_args = Dict(
          :output => :class,
          :strat => mode
       )
    )

Baseline model that returns the mode during classification.
"""
mutable struct Baseline <: TSLearner
    model
    args

    function Baseline(args=Dict())
        default_args = Dict(
            :output => :class,
            :strat => mode
        )
        new(nothing,mergedict(default_args,args))
    end
end

"""
    fit!(bsl::Baseline,x::Matrix,y::Vector)

Get the mode of the training data.
"""
function fit!(bsl::Baseline,x::Matrix,y::Vector)
    bsl.model = bsl.args[:strat](y)
end

"""
    transform!(bsl::Baseline,x::Matrix)

Return the mode in classification.
"""
function transform!(bsl::Baseline,x::Matrix)
    fill(bsl.model,size(x,1))
end

"""
    Identity(args=Dict())

Returns the input as output.
"""
mutable struct Identity <: Transformer
    model
    args

    function Identity(args=Dict())
        default_args = Dict{Symbol,Any}()
        new(nothing,mergedict(default_args,args))
    end
end

"""
    fit!(idy::Identity,x::Matrix,y::Vector)

Does nothing.
"""
function fit!(idy::Identity,x::Matrix,y::Vector)
    nothing
end

"""
    transform!(idy::Identity,x::Matrix)

Return the input as output.
"""
function transform!(idy::Identity,x::Matrix)
    return x
end

end


text:
# Raw memory management

export Mem, attribute, attribute!, memory_type, is_managed

module Mem

using ..CUDA
using ..CUDA: @enum_without_prefix, CUstream, CUdevice, CuDim3

using Base: @deprecate_binding

using IntervalTrees

#
# untyped buffers
#

abstract type Buffer end

# expected interface:
# - similar()
# - ptr, bytesize and ctx fields
# - convert() to Ptr and CuPtr

Base.pointer(buf::Buffer) = buf.ptr

Base.sizeof(buf::Buffer) = buf.bytesize

CUDA.device(buf::Buffer) = device(buf.ctx)

# ccall integration
#
# taking the pointer of a buffer means returning the underlying pointer,
# and not the pointer of the buffer object itself.
Base.unsafe_convert(T::Type{<:Union{Ptr,CuPtr}}, buf::Buffer) = convert(T, buf)

## device buffer

"""
    Mem.DeviceBuffer
    Mem.Device

A buffer of device memory residing on the GPU.
"""
struct DeviceBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::DeviceBuffer, ptr::CuPtr{Cvoid}=pointer(buf),
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    DeviceBuffer(ptr, bytesize, ctx)

Base.convert(::Type{<:Ptr}, buf::DeviceBuffer) =
    throw(ArgumentError("cannot take the CPU address of a GPU buffer"))

Base.convert(::Type{CuPtr{T}}, buf::DeviceBuffer) where {T} =
    convert(CuPtr{T}, pointer(buf))

"""
    Mem.alloc(DeviceBuffer, bytesize::Integer)

Allocate `bytesize` bytes of memory on the device. This memory is only accessible on the
GPU, and requires explicit calls to `unsafe_copyto!`, which wraps `cuMemcpy`,
for access on the CPU.
"""
function alloc(::Type{DeviceBuffer}, bytesize::Integer)
    bytesize == 0 && return DeviceBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CUDA.CUdeviceptr}()
    CUDA.cuMemAlloc(ptr_ref, bytesize)

    return DeviceBuffer(reinterpret(CuPtr{Cvoid}, ptr_ref[]), bytesize, CuCurrentContext())
end

function free(buf::DeviceBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFree(buf)
    end
end

## host buffer

"""
    Mem.HostBuffer
    Mem.Host

A buffer of pinned memory on the CPU, possible accessible on the GPU.
"""
struct HostBuffer <: Buffer
    ptr::Ptr{Cvoid}
    bytesize::Int
    ctx::CuContext

    mapped::Bool
end

Base.similar(buf::HostBuffer, ptr::Ptr{Cvoid}=pointer(buf), bytesize::Int=sizeof(buf),
             ctx::CuContext=buf.ctx, mapped::Bool=buf.mapped) =
    HostBuffer(ptr, bytesize, ctx, mapped)

Base.convert(::Type{Ptr{T}}, buf::HostBuffer) where {T} =
    convert(Ptr{T}, pointer(buf))

function Base.convert(::Type{CuPtr{T}}, buf::HostBuffer) where {T}
    if buf.mapped
        pointer(buf) == C_NULL && return convert(CuPtr{T}, CU_NULL)
        ptr_ref = Ref{CuPtr{Cvoid}}()
        CUDA.cuMemHostGetDevicePointer(ptr_ref, pointer(buf), #=flags=# 0)
        convert(CuPtr{T}, ptr_ref[])
    else
        throw(ArgumentError("cannot take the GPU address of a pinned but not mapped CPU buffer"))
    end
end

@deprecate_binding HOSTALLOC_DEFAULT 0 false
const HOSTALLOC_PORTABLE = CUDA.CU_MEMHOSTALLOC_PORTABLE
const HOSTALLOC_DEVICEMAP = CUDA.CU_MEMHOSTALLOC_DEVICEMAP
const HOSTALLOC_WRITECOMBINED = CUDA.CU_MEMHOSTALLOC_WRITECOMBINED

"""
    Mem.alloc(HostBuffer, bytesize::Integer, [flags])

Allocate `bytesize` bytes of page-locked memory on the host. This memory is accessible from
the CPU, and makes it possible to perform faster memory copies to the GPU. Furthermore, if
`flags` is set to `HOSTALLOC_DEVICEMAP` the memory is also accessible from the GPU.
These accesses are direct, and go through the PCI bus.
If `flags` is set to `HOSTALLOC_PORTABLE`, the memory is considered mapped by all CUDA contexts,
not just the one that created the memory, which is useful if the memory needs to be accessed from
multiple devices. Multiple `flags` can be set at one time using a bytewise `OR`:

    flags = HOSTALLOC_PORTABLE | HOSTALLOC_DEVICEMAP

"""
function alloc(::Type{HostBuffer}, bytesize::Integer, flags=0)
    bytesize == 0 && return HostBuffer(C_NULL, 0, CuContext(C_NULL), false)

    ptr_ref = Ref{Ptr{Cvoid}}()
    CUDA.cuMemHostAlloc(ptr_ref, bytesize, flags)

    mapped = (flags & HOSTALLOC_DEVICEMAP) != 0
    return HostBuffer(ptr_ref[], bytesize, CuCurrentContext(), mapped)
end

const HOSTREGISTER_PORTABLE = CUDA.CU_MEMHOSTREGISTER_PORTABLE
const HOSTREGISTER_DEVICEMAP = CUDA.CU_MEMHOSTREGISTER_DEVICEMAP
const HOSTREGISTER_IOMEMORY = CUDA.CU_MEMHOSTREGISTER_IOMEMORY

"""
    Mem.register(HostBuffer, ptr::Ptr, bytesize::Integer, [flags])

Page-lock the host memory pointed to by `ptr`. Subsequent transfers to and from devices will
be faster, and can be executed asynchronously. If the `HOSTREGISTER_DEVICEMAP` flag is
specified, the buffer will also be accessible directly from the GPU.
These accesses are direct, and go through the PCI bus.
If the `HOSTREGISTER_PORTABLE` flag is specified, any CUDA context can access the memory.
"""
function register(::Type{HostBuffer}, ptr::Ptr, bytesize::Integer, flags=0)
    bytesize == 0 && throw(ArgumentError("Cannot register an empty range of memory."))

    CUDA.cuMemHostRegister(ptr, bytesize, flags)

    mapped = (flags & HOSTREGISTER_DEVICEMAP) != 0
    return HostBuffer(ptr, bytesize, CuCurrentContext(), mapped)
end

"""
    Mem.unregister(HostBuffer)

Unregisters a memory range that was registered with [`Mem.register`](@ref).
"""
function unregister(buf::HostBuffer)
    CUDA.cuMemHostUnregister(buf)
end

function free(buf::HostBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFreeHost(buf)
    end
end

## unified buffer

"""
    Mem.UnifiedBuffer
    Mem.Unified

A managed buffer that is accessible on both the CPU and GPU.
"""
struct UnifiedBuffer <: Buffer
    ptr::CuPtr{Cvoid}
    bytesize::Int
    ctx::CuContext
end

Base.similar(buf::UnifiedBuffer, ptr::CuPtr{Cvoid}=pointer(buf),
             bytesize::Int=sizeof(buf), ctx::CuContext=buf.ctx) =
    UnifiedBuffer(ptr, bytesize, ctx)

Base.convert(::Type{Ptr{T}}, buf::UnifiedBuffer) where {T} =
    convert(Ptr{T}, reinterpret(Ptr{Cvoid}, pointer(buf)))

Base.convert(::Type{CuPtr{T}}, buf::UnifiedBuffer) where {T} =
    convert(CuPtr{T}, pointer(buf))

@enum_without_prefix CUDA.CUmemAttach_flags CU_MEM_

"""
    Mem.alloc(UnifiedBuffer, bytesize::Integer, [flags::CUmemAttach_flags])

Allocate `bytesize` bytes of unified memory. This memory is accessible from both the CPU and
GPU, with the CUDA driver automatically copying upon first access.
"""
function alloc(::Type{UnifiedBuffer}, bytesize::Integer,
              flags::CUDA.CUmemAttach_flags=ATTACH_GLOBAL)
    bytesize == 0 && return UnifiedBuffer(CU_NULL, 0, CuContext(C_NULL))

    ptr_ref = Ref{CuPtr{Cvoid}}()
    CUDA.cuMemAllocManaged(ptr_ref, bytesize, flags)

    return UnifiedBuffer(ptr_ref[], bytesize, CuCurrentContext())
end

function free(buf::UnifiedBuffer)
    if pointer(buf) != CU_NULL
        CUDA.cuMemFree(buf)
    end
end

"""
    prefetch(::UnifiedBuffer, [bytes::Integer]; [device::CuDevice], [stream::CuStream])

Prefetches memory to the specified destination device.
"""
function prefetch(buf::UnifiedBuffer, bytes::Integer=sizeof(buf);
                  device::CuDevice=device(buf), stream::CuStream=CuDefaultStream())
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDA.cuMemPrefetchAsync(buf, bytes, device, stream)
end

@enum_without_prefix CUDA.CUmem_advise CU_MEM_

"""
    advise(::UnifiedBuffer, advice::CUDA.CUmem_advise, [bytes::Integer]; [device::CuDevice])

Advise about the usage of a given memory range.
"""
function advise(buf::UnifiedBuffer, advice::CUDA.CUmem_advise, bytes::Integer=sizeof(buf);
                device::CuDevice=device(buf))
    bytes > sizeof(buf) && throw(BoundsError(buf, bytes))
    CUDA.cuMemAdvise(buf, bytes, advice, device)
end

## convenience aliases

const Device  = DeviceBuffer
const Host    = HostBuffer
const Unified = UnifiedBuffer

#
# typed pointers
#

## initialization

"""
    Mem.set!(buf::CuPtr, value::Union{UInt8,UInt16,UInt32}, len::Integer;
             async::Bool=false, stream::CuStream)

Initialize device memory by copying `val` for `len` times. Executed asynchronously if
`async` is true, in which case a valid `stream` is required.
"""
set!

for T in [UInt8, UInt16, UInt32]
    bits = 8*sizeof(T)
    fn_sync = Symbol("cuMemsetD$(bits)")
    fn_async = Symbol("cuMemsetD$(bits)Async")
    @eval function set!(ptr::CuPtr{$T}, value::$T, len::Integer;
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing)
        if async
          stream===nothing &&
              throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDA, fn_async))(ptr, value, len, stream)
        else
          stream===nothing ||
              throw(ArgumentError("Synchronous memory operations cannot be issues on a stream."))
            $(getproperty(CUDA, fn_sync))(ptr, value, len)
        end
    end
end

## copy operations

for (f, srcPtrTy, dstPtrTy) in (("cuMemcpyDtoH", CuPtr, Ptr),
                                ("cuMemcpyHtoD", Ptr,   CuPtr),
                                ("cuMemcpyDtoD", CuPtr, CuPtr),
                               )
    @eval function Base.unsafe_copyto!(dst::$dstPtrTy{T}, src::$srcPtrTy{T}, N::Integer;
                                       stream::Union{Nothing,CuStream}=nothing,
                                       async::Bool=false) where T
        if async
            stream===nothing &&
                throw(ArgumentError("Asynchronous memory operations require a stream."))
            $(getproperty(CUDA, Symbol(f * "Async")))(dst, src, N*sizeof(T), stream)
        else
            stream===nothing ||
                throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
            $(getproperty(CUDA, Symbol(f)))(dst, src, N*sizeof(T))
        end
        return dst
    end
end

"""
    unsafe_copy3d!(dst, dstTyp, src, srcTyp, width, height=1, depth=1;
                   dstPos=(1,1,1), dstPitch=0, dstHeight=0,
                   srcPos=(1,1,1), srcPitch=0, srcHeight=0,
                   async=false, stream=nothing)

Perform a 3D memory copy between pointers `src` and `dst`, at respectively position `srcPos`
and `dstPos` (1-indexed). Both pitch and destination can be specified for both the source
and destination; consult the CUDA documentation for more details. This call is executed
asynchronously if `async` is set, in which case `stream` needs to be a valid CuStream.
"""
function unsafe_copy3d!(dst::Union{Ptr{T},CuPtr{T}}, dstTyp::Type{<:Buffer},
                        src::Union{Ptr{T},CuPtr{T}}, srcTyp::Type{<:Buffer},
                        width::Integer, height::Integer=1, depth::Integer=1;
                        dstPos::CuDim=(1,1,1), srcPos::CuDim=(1,1,1),
                        dstPitch::Integer=0, dstHeight::Integer=0,
                        srcPitch::Integer=0, srcHeight::Integer=0,
                        async::Bool=false, stream::Union{Nothing,CuStream}=nothing) where T
    srcPos = CUDA.CuDim3(srcPos)
    dstPos = CUDA.CuDim3(dstPos)

    srcMemoryType, srcHost, srcDevice = if srcTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        src::Ptr,
        CU_NULL
    elseif srcTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        C_NULL,
        src::CuPtr
    elseif srcTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        C_NULL,
        reinterpret(CuPtr{Cvoid}, src)
    end
    srcArray = C_NULL

    dstMemoryType, dstHost, dstDevice = if dstTyp == Host
        CUDA.CU_MEMORYTYPE_HOST,
        dst::Ptr,
        CU_NULL
    elseif dstTyp == Mem.Device
        CUDA.CU_MEMORYTYPE_DEVICE,
        C_NULL,
        dst::CuPtr
    elseif dstTyp == Mem.Unified
        CUDA.CU_MEMORYTYPE_UNIFIED,
        C_NULL,
        reinterpret(CuPtr{Cvoid}, dst)
    end
    dstArray = C_NULL

    params_ref = Ref(CUDA.CUDA_MEMCPY3D(
        # source
        (srcPos.x-1)*sizeof(T), srcPos.y-1, srcPos.z-1,
        0, # LOD
        srcMemoryType, srcHost, srcDevice, srcArray,
        C_NULL, # reserved
        srcPitch, srcHeight,
        # destination
        (dstPos.x-1)*sizeof(T), dstPos.y-1, dstPos.z-1,
        0, # LOD
        dstMemoryType, dstHost, dstDevice, dstArray,
        C_NULL, # reserved
        dstPitch, dstHeight,
        # extent
        width*sizeof(T), height, depth
    ))
    if async
        stream===nothing &&
            throw(ArgumentError("Asynchronous memory operations require a stream."))
        CUDA.cuMemcpy3DAsync_v2(params_ref, stream)
    else
        stream===nothing ||
            throw(ArgumentError("Synchronous memory operations cannot be issued on a stream."))
        CUDA.cuMemcpy3D_v2(params_ref)
    end
end

#
# auxiliary functionality
#

## memory pinning

const pinned_memory = IntervalTree{Int, Interval{Int}}()

function pin(ptr::Ptr, bytesize::Integer, flags=0)
    low = Int(ptr)
    high = Int(ptr+bytesize)

    # register gaps between already-registered slices
    intersections = [Interval{Int}(0, low),
                     intersect(pinned_memory, Interval{Int}(low, high))...,
                     Interval{Int}(high, 0)]
    for (lower, upper) in zip(intersections[1:end-1], intersections[2:end])
        gap = Interval{Int}(last(lower), first(upper))
        if first(gap) < last(gap)
            gap_ptr = Ptr{Nothing}(first(gap))
            gap_bytesize = last(gap) - first(gap)
            Mem.register(Mem.Host, gap_ptr, gap_bytesize, flags)
        end
    end

    # delete actual intersections
    low = min(low, first(intersections[2]))
    high = max(high, last(intersections[end-1]))
    for range in intersections[2:end-1]
        delete!(pinned_memory, range)
    end

    push!(pinned_memory, Interval{Int}(low, high))
    return
end

pin(a::Union{Array,SubArray{<:Any,<:Any,<:Array},Base.ReshapedArray{<:Any,<:Any,<:Array}}, flags=0) =
    pin(pointer(a), sizeof(a), flags)

## memory info

function info()
    free_ref = Ref{Csize_t}()
    total_ref = Ref{Csize_t}()
    CUDA.cuMemGetInfo(free_ref, total_ref)
    return convert(Int, free_ref[]), convert(Int, total_ref[])
end

end # module Mem

"""
    available_memory()

Returns the available_memory amount of memory (in bytes), available for allocation by the CUDA context.
"""
available_memory() = Mem.info()[1]

"""
    total_memory()

Returns the total amount of memory (in bytes), available for allocation by the CUDA context.
"""
total_memory() = Mem.info()[2]

## pointer attributes

"""
    attribute(X, ptr::Union{Ptr,CuPtr}, attr)

Returns attribute `attr` about pointer `ptr`. The type of the returned value depends on the
attribute, and as such must be passed as the `X` parameter.
"""
function attribute(X::Type, ptr::Union{Ptr{T},CuPtr{T}}, attr::CUpointer_attribute) where {T}
    ptr = reinterpret(CuPtr{T}, ptr)
    data_ref = Ref{X}()
    cuPointerGetAttribute(data_ref, attr, ptr)
    return data_ref[]
end

"""
    attribute!(ptr::Union{Ptr,CuPtr}, attr, val)

Sets attribute` attr` on a pointer `ptr` to `val`.
"""
function attribute!(ptr::Union{Ptr{T},CuPtr{T}}, attr::CUpointer_attribute, val) where {T}
    ptr = reinterpret(CuPtr{T}, ptr)
    cuPointerSetAttribute(Ref(val), attr, ptr)
    return
end

@enum_without_prefix CUpointer_attribute CU_

# some common attributes

@enum_without_prefix CUmemorytype CU_
memory_type(x) = CUmemorytype(attribute(Cuint, x, POINTER_ATTRIBUTE_MEMORY_TYPE))

is_managed(x) = convert(Bool, attribute(Cuint, x, POINTER_ATTRIBUTE_IS_MANAGED))