text:
"""
Computes node betweenness centrality from a given source.

Input:
-----
* G::SparseMatrixCSC
* root::Int

Output:
------
* bc_scores::Vector{Float32}

"""
function bc(G::SparseMatrixCSC, root::Int)
	rows = G.colptr - 1
    cols = G.rowval - 1
    nodes = length(rows) - 1
    edges = length(cols)
    rows = map(Int32, rows)
    cols = map(Int32, cols)
    bc_scores = Array(Cfloat, nodes)
	ccall((:bc, gunrock_lib), Void, (Ptr{Cfloat}, Int, Int, Ptr{Cint}, Ptr{Cint}, Int), bc_scores, nodes, edges, rows, cols, root - 1)
	bc_scores
end


text:
# Get status of the registration
function status(r::HTTP.Request)
    r.method == "GET" || return json(405; error="Method not allowed")

    # Only allow authenticated users
    state = getcookie(r, "state")
    if isempty(state) || !haskey(USERS, state)
        return json(400; error="Missing or invalid state cookie")
    end

    id = getquery(r, "id")
    haskey(REGISTRATIONS, id) || return json(404; state="unknown")
    regstate = REGISTRATIONS[id]
    regstate.state === :pending || delete!(REGISTRATIONS, id)
    return json(; state=regstate.state, message=regstate.msg)
end

ext: jl


text:
# adapted from https://discourse.julialang.org/t/how-to-copy-all-fields-without-changing-the-referece/945/5
@inline @generated function particleCopy!(dest, src)
  fieldNames = fieldnames(dest)
  fieldTypes = dest.types
  numFields = length(fieldNames)
  expressions = Array{Expr}(undef, numFields)

  for i = 1:numFields
    fieldName = fieldNames[i]
    fieldType = fieldTypes[i]
    @assert !fieldType.mutable || hasmethod(copyto!, (fieldType,
      fieldType)) "$fieldName::$fieldType : copyto! must exist for mutable Particle fields"
    if fieldType.mutable
      @inbounds expressions[i] = :(copyto!(dest.$fieldName, src.$fieldName))
    else
      @inbounds expressions[i] = :(dest.$fieldName = src.$fieldName)
    end
  end
  body = Expr(:block, expressions...)

  quote
    $body
    return
  end
end

## better not to inline ; enough work should be done in the loop
function _copyParticles!(out::Vector{Particle}, in::Vector{Particle}) where Particle
  for i in eachindex(out)
    @inbounds particleCopy!(out[i], in[i])
  end
end

## better not to inline ; enough work should be done in the loop
function _copyParticles!(zetaAncs::Vector{Particle}, zetas::Vector{Particle},
  as::Vector{Int64}) where Particle
  for i in eachindex(zetaAncs)
    @inbounds particleCopy!(zetaAncs[i], zetas[as[i]])
  end
end

## better not to inline ; enough work should be done in the loop
function _mutateParticles!(zetas::Vector{Particle}, rng::RNG, p::Int64, M!::F,
  zetaAncs::Vector{Particle}, pScratch::ParticleScratch) where
  {Particle, F<:Function, ParticleScratch}
  for j in eachindex(zetas)
    @inbounds M!(zetas[j], rng, p, zetaAncs[j], pScratch)
  end
end

## better not to inline ; enough work should be done in the loop
function _mutateParticles!(zetas::Vector{Particle}, rng::RNG, p::Int64, M!::F,
  zetaAncs::Vector{Particle}, pScratch::ParticleScratch, xref::Particle) where {Particle,
  F<:Function, ParticleScratch}
  particleCopy!(zetas[1], xref)
  for j in 2:length(zetas)
    @inbounds M!(zetas[j], rng, p, zetaAncs[j], pScratch)
  end
end

## better not to inline ; enough work should be done in the loop
function _logWeightParticles!(lws::Vector{Float64}, p::Int64, lG::F,
  zetas::Vector{Particle}, pScratch::ParticleScratch) where
  {F<:Function, Particle, ParticleScratch}
  for j in eachindex(zetas)
    @inbounds lws[j] = lG(p, zetas[j], pScratch)
  end
end

@inline function _intermediateOutput!(smcio::SMCIO, p::Int64)
  _copyParticles!(smcio.allZetas[p], smcio.zetas)
  @inbounds smcio.allWs[p] .= smcio.ws
  @inbounds p != smcio.n && (smcio.allAs[p] .= smcio.internal.as)
end

@inline function pickParticle!(path::Vector{Particle}, smcio::SMCIO{Particle}) where Particle
  @assert length(path) == smcio.n
  n::Int64 = smcio.n
  allAs::Vector{Vector{Int64}} = smcio.allAs
  allZetas::Vector{Vector{Particle}} = smcio.allZetas
  k::Int64 = sampleCategorical(smcio.ws, getRNG())

  @inbounds particleCopy!(path[n], allZetas[n][k])
  for p = n-1:-1:1
    @inbounds k = allAs[p][k]
    @inbounds particleCopy!(path[p], allZetas[p][k])
  end
  return
end

@inline function pickParticleBS!(path::Vector{Particle},
  smcio::SMCIO{Particle, ParticleScratch}, lM::F) where {Particle,
  ParticleScratch, F<:Function}
  @assert length(path) == smcio.n
  n::Int64 = smcio.n
  N::Int64 = smcio.N
  allAs::Vector{Vector{Int64}} = smcio.allAs
  allZetas::Vector{Vector{Particle}} = smcio.allZetas
  allWs::Vector{Vector{Float64}} = smcio.allWs
  bws::Vector{Float64} = smcio.internal.scratch1
  pScratch::ParticleScratch = smcio.internal.particleScratch

  rng = getRNG()
  k::Int64 = sampleCategorical(smcio.ws, rng)

  @inbounds particleCopy!(path[n], allZetas[n][k])
  for p = n-1:-1:1
    @inbounds bws .= log.(allWs[p])
    for j in 1:N
      @inbounds bws[j] += lM(p+1, allZetas[p][j], path[p+1], pScratch)
    end
    m::Float64 = maximum(bws)
    bws .= exp.(bws .- m)
    k = sampleCategorical(bws, rng)
    @inbounds particleCopy!(path[p], allZetas[p][k])
  end
end


text:
using CVChannel
using Test
"""
This script verifies the Siddhu channel when
tensored with the depolarizing channel is multiplicative.
"""

print("\nHere, we look at the communication value of running the")
print(" Siddhu and the d-dimensional depolarizing channel in parallel for d = {2,3}.")
@testset "Verifying Multiplicativity of Siddhu Channel with d∈[2,3] Depolarizing Channel" begin
    for d in [2:3;]
        println("~~~~Now scanning for d = ",d)
        s_scan_range =[0:0.1:0.5;]
        q_scan_range = [0:0.1:1;]
        max_val = 0
        is_mult = true
        for s in s_scan_range
            println("---Now scanning for s = ", s, "-----")
            for q in q_scan_range
                target_val = 2*(d*(1-q)+q)

                depolqChan(X) = depolarizingChannel(X,q)
                sidchan(X) = siddhuChannel(X,s)

                depolq_chan = Choi(depolqChan,d,d)
                sid_chan= Choi(sidchan,3,3)

                par_cv, = pptCV(parChoi(sid_chan, depolq_chan), :dual)

                non_mult = par_cv - target_val
                !isapprox(non_mult,0,atol=3e-6) ? is_mult = false : nothing
                max_val < abs(non_mult) ? max_val = abs(non_mult) : nothing
            end
        end
        @test is_mult
    end
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
function fraction_to_repeated(time_fraction::Real, n_repeats::Int, yoyo::Bool, prewait, postwait)
    # the 1 values should be reached, so no simple modulo
    # that means there are no zeros after the first, but 1 at the end is more important

    if time_fraction == 0
        return 0
    elseif time_fraction == 1
        # if there is an even yoyo value it needs to come down to zero at the end
        if yoyo && n_repeats % 2 == 0
            return 0
        else
            return 1
        end
    elseif time_fraction > 1 || time_fraction < 0
        error("Time fraction is $time_fraction but should be from 0 to 1")
    end

    if !(0 <= prewait < 1)
        error("Pre-wait fraction $prewait is invalid")
    end
    if !(0 <= postwait < 1)
        error("Post-wait fraction $postwait is invalid")
    end
    if postwait + prewait >= 1
        error("Pre-wait is $prewait and post-wait is $postwait, together larger or equal to 1.")
    end

    multiplied = time_fraction * n_repeats
    multiples, rest = divrem(multiplied, 1)
    nth_repeat = Int(multiples) + 1 # because 1st repeat is 0 for divrem

    # this causes the value to be 0 for the prewait interval and 1 for the postwait interval
    repeat_fraction = clamp((rest - prewait) / ((1 - postwait) - prewait), 0, 1)

    if yoyo
        # 1: 0 to 1
        # 2: 1 to 0
        # 3: 0 to 1
        # etc
        if isodd(nth_repeat) # first up
            return repeat_fraction
        else
            return 1 - repeat_fraction # then down
        end
    else
        # this is different than standard modulo
        return rest == 0 ? 1 : repeat_fraction
    end
end

function interpolate(easing::Easing, t::Real, k1::Keyframe{T}, k2::Keyframe{T}) where T

    # the fraction of the keyframe interval we're at
    time_fraction_unclamped = (t - k1.t) / (k2.t - k1.t)

    # this is so the even yoyos work, the time value needs to go through the easing
    # but at least it's clamped
    time_fraction = clamp(time_fraction_unclamped, 0, 1)

    # if time_fraction <= 0
    #     return k1.value
    # elseif time_fraction >= 1
    #     return k2.value
    # end

    # handle repetitions
    # on the previous interval of 0 to 1 there are now n 0 to 1 intervals
    # this is the interval from 0 to 1 on which the easing function will be applied
    interp_fraction = fraction_to_repeated(time_fraction, easing.n, easing.yoyo, easing.prewait, easing.postwait)

    # this is the actual ratio of the two keyframe values that comes out of the esas
    interp_ratio = interpolation_ratio(easing.easing, interp_fraction)

    # these checks enable to return early if values are 0 or 1, which is why
    # NoEasing EasingType can be used for non-interpolateable values like strings
    if interp_ratio == 0
        return k1.value
    elseif interp_ratio == 1
        return k2.value
    end

    interp_value = linear_interpolate(interp_ratio, k1.value, k2.value)
end

# this should be overloaded for weird types
function linear_interpolate(fraction::Real, value1::T, value2::T) where T
    (value2 - value1) * fraction + value1
end

# array version with broadcasting
function linear_interpolate(fraction::Real, array1::T, array2::T) where T <: AbstractArray
    linear_interpolate.(fraction, array1, array2)
end

function linear_interpolate(fraction::Real, c1::T, c2::T) where T <: Colors.Colorant
    Colors.weighted_color_mean(1 - fraction, c1, c2)
end

function interpolation_ratio(easing::FuncEasing, fraction)
    easing.f(fraction, easing.args...)
end

function interpolation_ratio(easing::MixedEasing, fraction)
    e1 = easing.e1
    e2 = easing.e2
    interp_fraction1 = fraction_to_repeated(fraction, e1.n, e1.yoyo, e1.prewait, e1.postwait)
    interp_fraction2 = fraction_to_repeated(fraction, e2.n, e2.yoyo, e2.prewait, e2.postwait)
    interpolation_ratio(easing.e1.easing, interp_fraction1) * easing.mix + interpolation_ratio(easing.e2.easing, interp_fraction2) * (1 - easing.mix)
end

function interpolation_ratio(easing::MultipliedEasing, fraction)
    e1 = easing.e1
    e2 = easing.e2
    interp_fraction1 = fraction_to_repeated(fraction, e1.n, e1.yoyo, e1.prewait, e1.postwait)
    interp_fraction2 = fraction_to_repeated(fraction, e2.n, e2.yoyo, e2.prewait, e2.postwait)
    interpolation_ratio(easing.e1.easing, interp_fraction1) * interpolation_ratio(easing.e2.easing, interp_fraction2)
end

function interpolation_ratio(easing::EasedEasing, fraction)
    e1 = easing.e1
    e2 = easing.e2
    ease = easing.easing
    interp_fraction1 = fraction_to_repeated(fraction, e1.n, e1.yoyo, e1.prewait, e1.postwait)
    interp_fraction2 = fraction_to_repeated(fraction, e2.n, e2.yoyo, e2.prewait, e2.postwait)
    interp_fraction_ease = fraction_to_repeated(fraction, ease.n, ease.yoyo, ease.prewait, ease.postwait)
    final = interpolation_ratio(easing.e1.easing, interp_fraction1) * (1 - interp_fraction_ease) +
            interpolation_ratio(easing.e2.easing, interp_fraction2) * interp_fraction_ease
end


text:
module Arblib

include(joinpath(@__DIR__, "..", "deps", "deps.jl"))

function __init__()
    check_deps()
end

export Arf, Arb, Acb

import Base: isfinite, isinf, isinteger, isnan, isone, isreal, iszero

include("macros.jl")
include("arb_types.jl")
include("types.jl")
include("rounding.jl")
include("precision.jl")

include("constructors.jl")

include("predicates.jl")
include("show.jl")

include("arbcall.jl")

include("arb.jl")

end # module


text:
function HermiteReduce(f::AbstractAlgebra.ResFieldElem{P}, DE::AlgebraicExtensionDerivation) where
    {T<:FieldElement, P<:PolyElem{T}}
    # "Lazy" Hermite reduction, see Section 2.1 of:
    # Manuel Bronstein. Symbolic integration tutorial. ISSAC’98, 1998.
    # http://www-sop.inria.fr/cafe/Manuel.Bronstein/publications/issac98.pdf
    iscompatible(f, DE) || error("rational function f must be in the domain of derivation D")

    E = parent(f)
    F = base_ring(base_ring(E))
    y = E(gen(base_ring(E)))
    p = modulus(E)
    n = degree(p)
    M_n_n = MatrixSpace(F, n, n)
    M_n = MatrixSpace(F, n, 1)
    g = zero(E)

    if all([iszero(coeff(p, j)) for j=1:n-1]) # simple radical extension
        # Simpler method for simple radical extensions, see Section 2.2. of Bronstein's tutorial.
        @assert isone(leading_coefficient(p))
        a = constant_coefficient(p)
        A = numerator(a)
        D = denominator(a)
        As = Squarefree(A*D^(n-1))
        k = length(As)
        qrs = [divrem(i, n) for i=1:k]
        qs = [q for (q, r) in qrs]
        rs = [r for (q, r) in qrs]
        F = prod([As[i]^qs[i] for i=1:k])
        H = prod([As[i]^rs[i] for i=1:k])
        z = y*D//F
        Hs = Squarefree(H)
        m = length(Hs)
        ws = [z^(i-1)*1//prod([Hs[j]^div((i-1)*j, n) for j=1:m]) for i=1:n]
        Ds = [prod([Hs[j]^div((i-1)*j, n) for j=1:m]) for i=1:n] # note: index shift in i
        dDs_over_Ds = [derivative(D)//D for D in Ds]
        dH_over_H = derivative(H)//H
        Winv = inv(M_n_n([coeff(data(w), j) for j=0:n-1, w in ws]))
        while true
            fs = Winv*M_n([coeff(data(f), i) for i=0:n-1])
            D = lcm(denominator.(fs)[:,1])
            As = [numerator(f*D) for f in fs]
            if degree(gcd(D, derivative(D)))<= 0 #  D squarefree?
                return g, f, ws, As, D
            end

            Ds1 = Squarefree(D) # "Ds1" because "Ds" already in use and needed below
            m = length(Ds1)
            V = Ds1[m]
            U = divexact(D, V^m)
            dV = derivative(V)

            fs = [As[i]//(U*(V*((i-1)//n*dH_over_H - dDs_over_Ds[i]) - (m-1)*dV)) for i=1:n]
            Q = lcm(denominator.(fs)[:,1])
            Ts = [numerator(f*Q) for f in fs]
            _, _, R = gcdx(V, Q)
            B = sum([rem(R*Ts[i], V)*ws[i] for i=1:n])*1//V^(m - 1)
            g += B
            f -= DE(B)
        end
    end

    ws = [y^j for j=0:n-1]
    Winv = M_n_n(1)

    while true
        fs = Winv*M_n([coeff(data(f), i) for i=0:n-1])
        D = lcm(denominator.(fs)[:,1])
        As = [numerator(f*D) for f in fs]
        if degree(gcd(D, derivative(D)))<= 0 #  D squarefree?
            return g, f, ws, As, D
        end

        Ds = Squarefree(D)
        m = length(Ds)
        V = Ds[m]
        U = divexact(D, V^m)

        Ss = [U*V^m*DE(w*1//V^(m-1)) for w in ws]

        SM = M_n_n([coeff(data(S), j) for j=0:n-1, S in Ss])
        d, N = nullspace(SM)
        if d>0 # Theorem 1
            _, N = nullspace(SM)
            w0 = U//V*sum(N[i,1]*ws[i] for i=1:n)
        else
            fs = solve(SM, M_n(As))
            Q = lcm(denominator.(fs)[:,1])
            Ts = [numerator(f*Q) for f in fs]
            w0 = zero(parent(f))
            gcdVQ = gcd(V, Q)
            if degree(gcdVQ)>=1 # Theorem 2
                w0 = U*divexact(V, gcdVQ)//gcdVQ*sum([Ts[i]*ws[i] for i=1:n])
            else
                j = 1
                while j<=n
                    dw = DE(ws[j])
                    dws = Winv*M_n([coeff(data(dw), i) for i=0:n-1])
                    F = lcm(denominator.(dws)[:,1])
                    if degree(gcd(F, derivative(F)))>=1 # denominator F of ws[j]' not squarefree, apply Theorem 3
                        w0 = prod(Squarefree(F))*dw
                        break
                    end
                    j += 1
                end
                if j>n # all denominators of the ws[j]' squarefree, do the reduction
                    _, _, R = gcdx(V, Q)
                    B = sum([rem(R*Ts[i], V)*ws[i] for i=1:n])*1//V^(m - 1)
                    g += B
                    f -= DE(B)
                    ws = [y^j for j=0:n-1]
                    Winv = M_n_n(1)
                    continue
                end
            end
        end

        C0s = solve(M_n_n([coeff(data(w), j) for j=0:n-1, w in ws]), M_n([coeff(data(w0), i) for i=0:n-1]))
        F = lcm(denominator.(C0s)[:,1])
        Cs = [numerator(C0*F) for C0 in C0s]

        C = zero_matrix(parent(F), n+1, n)
        for j=1:n
            C[j,j] = F
        end
        for j=1:n
            C[n+1,j] = Cs[j]
        end
        C = hnf(C)

        ws = [sum([C[i,j]*ws[j] for j=1:n])*1//F for i=1:n]
        Winv = inv(M_n_n([coeff(data(w), j) for j=0:n-1, w in ws]))
    end
end

function IntegralBasis(E::AbstractAlgebra.ResField{P}) where {T<:FieldElement, P<:PolyElem{T}}
    # Trager's algorithm, see Chapter 2 of
    # B.M. Trager. On the integration of algebraic functions. PhD thesis, MIT, Computer Science, 1984.
    # https://dspace.mit.edu/bitstream/handle/1721.1/15391/12487590-MIT.pdf
    Ky = base_ring(E)
    K = base_ring(Ky)
    y = E(gen(Ky))
    f = modulus(E)
    n = degree(f)
    M_n_n = MatrixSpace(K, n, n)
    MP_n_n = MatrixSpace(base_ring(K), n, n)
    D = resultant(f, derivative(f))
    @assert isone(denominator(D))
    D = numerator(D)
    D = 1//leading_coefficient(D) * D
    k = D
    bs = [y^j for j=0:n-1]
    B = M_n_n([coeff(data(bs[j]), i) for i=0:n-1, j=1:n])
    while true
        Ds = SI.Squarefree(D)
        if length(Ds)<2
            return bs
        end
        Q = gcd(prod(SI.Squarefree(D)[2:end]), k)
        if degree(Q)<=0
            return bs
        end

        # Compute J, the Q-trace radical of V
        ZE = zero(E)
        S = [i<=j ? bs[i]*bs[j] : ZE  for i=1:n, j=1:n] # only upper triangle needed
        ZK = zero(K)
        # TODO: optimize computation of TM, which is one of the bottlenecks.
        # Compute only upper triangle of TM, lower triangle by symmetry:
        TM = [i<=j ? numerator(sum([coeff(data(S[i, j]*y^l), l) for l=0:n-1])) : ZK for i=1:n, j=1:n]
        for i=1:n
            for j=1:i-1
                TM[i, j] = TM[j, i]
            end
        end
        H = hnf(vcat(MP_n_n(TM), MP_n_n(Q)))[1:n,:] # Hermite normal form
        J = inv(map(x->x//1, H)) # map(x->x/1, H) is H as element in FractionField

        # Compute the idealizer of J
        JB = J*B
        ms = [E(Ky(reshape(Matrix(JB[:,k]), n))) for k=1:n]
        Minv = inv(M_n_n([coeff(data(m), i) for i=0:n-1, m in ms]))
        MM = 0
        for k=1:n
            mbs = [data(ms[k]*b) for b in bs]
            Mk = Minv*M_n_n([coeff(mbs[j], i) for i=0:n-1, j=1:n])
            MM = k==1 ? Mk : vcat(MM, Mk)
        end
        Mhat = hnf(map(x->numerator(x), MM))[1:n,:] # Hermite normal form
        Mhatinv = inv(map(x->x//1, Mhat))
        k = det(Mhat)

        if degree(k)<=0 # is k a unit ?
            return bs
        end

        #Update bs by applying the change of basis
        B = transpose(Mhatinv)*B
        bs = [E(Ky(reshape(Matrix(B[:,k]), n))) for k=1:n]
        D = divexact(D, k^2)
    end
end


text:
lines = readlines("input.txt")
field = zeros(Int, (1024,1024))
re = r"(\d+)"

for line = lines
    x1,y1,x2,y2 = parse.(Int, m.match for m = eachmatch(re, line))
    if x1 != x2 && y1 != y2
        continue
    end

    idx1,idx2,idy1,idy2 = min(x1,x2)+1, max(x1,x2)+1, min(y1,y2)+1, max(y1,y2)+1
    field[idy1:idy2, idx1:idx2] .+= 1
end

println(sum(@. field > 1))


text:
# In the future we wil validate namelists using JSONSchema

"""
    isvalid(namelist)

Check that `namelist` is valid against Harmonie.schema
"""
isvalid(namelist::Dict)  = JSONSchema.isvalid(namelist,schema)

"""
    diagnose(namelist)

Check that `namelist` is valid against Harmonie.schema. If
valid return `nothing`, and if not, return a diagnostic String containing a
selection of one or more likely causes of failure.
"""
diagnose(namelist::Dict) = JSONSchema.diagnose(namelist,schema)


text:
module OverlapIntegrals

export fact2, binomial_prefactor, overlap1d, dist2, product_center_1d, tho66

function fact2(n)
    if n <= 0
        return 1
    else
        return n * fact2(n - 2)
    end
end

function binomial_prefactor(s::Integer, ia::Integer, ib::Integer, xpa::Number, xpb::Number)
    acc = 0
    for t = 0 : s + 1
        if (s - ia) <= t && t <= ib
            acc += binomial(ia, s - t) * binomial(ib, t) * (xpa ^ (ia - s + t)) * (xpb ^ (ib - t))
        end
    end
    acc
end

function overlap1d(l1::Integer, l2::Integer, pax::Number, pbx::Number, gamma::Number)
    acc = 0.0
    for i = 0 : 1 + floor(0.5 * (l1 + l2))
        acc += binomial_prefactor(Integer(2 * i), l1, l2, pax, pbx) * fact2(2 * i - 1) / ((2 * gamma) ^ i)
    end
    acc
end

function dist2(xa, ya, za, xb, yb, zb)
    (xa - xb) * (xa - xb) + (ya - yb) * (ya - yb) + (za - zb) * (za - zb)
end

function product_center_1d(za, xa, zb, xb)
    ((za * xa) + (zb * xb)) / (za + zb)
end

function tho66(alpha1::Number, alpha2::Number, ra, rb, la, lb)
    gamma = alpha1 + alpha2
    xa, ya, za = ra
    xb, yb, zb = rb
    rab2 = dist2(xa, ya, za, xb, yb, zb)
    l1, m1, n1 = la
    l2, m2, n2 = lb
    xp = product_center_1d(alpha1, xa, alpha2, xb)
    yp = product_center_1d(alpha1, ya, alpha2, yb)
    zp = product_center_1d(alpha1, za, alpha2, zb)

    pre = exp(-alpha1 * alpha2 * rab2 / gamma) * (pi / gamma) ^ 1.5

    wx = overlap1d(l1, l2, xp - xa, xp - xb, gamma)
    wy = overlap1d(m1, m2, yp - ya, yp - yb, gamma)
    wz = overlap1d(n1, n2, zp - za, zp - zb, gamma)

    pre * wx * wy * wz
end

end

max_stars_repo_name: berquist/overlapintegrals

max_issues_repo_name: berquist/overlapintegrals

max_forks_repo_name: berquist/overlapintegrals


text:
function readdocumentrecord(io::IO)
    n = read(io, Int32)
    skip(io, 80 * n)
    return read(io, Int32)
end


text:
import DataFrames
import GLM
import StatsModels

function GLMModel(
        formula::StatsModels.AbstractTerm,
        family::GLM.Distribution,
        link::GLM.Link;
        name::AbstractString = "",
        isclassificationmodel::Bool = false,
        isregressionmodel::Bool = false,
        )
    underlyingglm = FitNotYetRunUnderlyingObject()
    result = GLMModel(
        name,
        isclassificationmodel,
        isregressionmodel,
        formula,
        family,
        link,
        underlyingglm,
        )
    return result
end

"""
"""
function get_underlying(x::GLMModel;
                        saving::Bool = false,
                        loading::Bool = false)
    return x.underlyingglm
end

"""
"""
function fit!(
        estimator::GLMModel,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        )
    labelsandfeatures_df = hcat(labels_df, features_df)
    @info(string("Starting to train GLM model."))
    glm = try
        GLM.glm(
            estimator.formula,
            labelsandfeatures_df,
            estimator.family,
            estimator.link,
            )
    catch e
        @warn(
            string(
                "while training GLM model, ignored error: ",
                e,
                )
            )
        FitFailedUnderlyingObject()
    end
    # glm =
    @info(string("Finished training GLM model."))
    estimator.underlyingglm = glm
    return estimator
end

"""
"""
function predict(
        estimator::GLMModel,
        features_df::DataFrames.AbstractDataFrame,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            features_df,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc
            )
        result = DataFrames.DataFrame()
        label_name = estimator.formula.lhs.sym
        result[label_name] = predictionsvector
        return result
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        if isa(estimator.underlyingglm, AbstractNonExistentUnderlyingObject)
            glmpredictoutput = fill(Float64(0), size(features_df,1))
        else
            glmpredictoutput = GLM.predict(
                estimator.underlyingglm,
                features_df,
                )
        end
        result = DataFrames.DataFrame()
        label_name = estimator.formula.lhs.sym
        result[label_name] = glmpredictoutput
        return result
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function predict(
        estimator::GLMModel,
        features_df::DataFrames.AbstractDataFrame,
        positive_class::Integer,
        threshold::AbstractFloat,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        probabilitiesassoc = predict_proba(
            estimator,
            features_df,
            )
        predictionsvector = single_labelprobabilitiestopredictions(
            probabilitiesassoc,
            positive_class,
            threshold,
            )
        result = DataFrames.DataFrame()
        label_name = estimator.formula.lhs.sym
        result[label_name] = predictionsvector
        return result
    else
        error(
            "Can only use the `threshold` argument with classification models"
            )
    end
end

"""
"""
function predict_proba(
        estimator::GLMModel,
        features_df::DataFrames.AbstractDataFrame,
        )
    if estimator.isclassificationmodel && !estimator.isregressionmodel
        if isa(estimator.underlyingglm, AbstractNonExistentUnderlyingObject)
            glmpredictoutput = fill(Float64(0), size(features_df, 1))
        else
            glmpredictoutput = GLM.predict(
                estimator.underlyingglm,
                features_df,
                )
        end
        result = Dict()
        result[1] = glmpredictoutput
        result[0] = 1 .- glmpredictoutput
        result = fix_type(result)
        return result
    elseif !estimator.isclassificationmodel && estimator.isregressionmodel
        error("predict_proba is not defined for regression models")
    else
        error(
            "Could not figure out if model is classification or regression"
            )
    end
end

"""
"""
function singlelabelbinaryclassdataframelogisticclassifier_GLM(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        intercept::Bool = true,
        interactions::Integer = 1,
        name::AbstractString = "",
        )
    negative_class = single_label_levels[1]
    positive_class = single_label_levels[2]
    formula = generate_formula(
        [single_label_name],
        feature_names;
        intercept = intercept,
        interactions = interactions,
        )
    dftransformer =
        ImmutableDataFrame2GLMSingleLabelBinaryClassTransformer(
            single_label_name,
            positive_class,
            )
    glmestimator = GLMModel(
        formula,
        GLM.Binomial(),
        GLM.LogitLink();
        isclassificationmodel = true,
        isregressionmodel = false,
        )
    predictlabelfixer =
        ImmutablePredictionsSingleLabelInt2StringTransformer(
            0,
            single_label_levels,
            )
    predprobalabelfixer =
        ImmutablePredictProbaSingleLabelInt2StringTransformer(
            0,
            single_label_levels,
            )
    probapackager = ImmutablePackageSingleLabelPredictProbaTransformer(
        single_label_name,
        )
    finalpipeline = dftransformer |>
                    glmestimator |>
                    predictlabelfixer |>
                    predprobalabelfixer |>
                    probapackager
    finalpipeline.name = name
    return finalpipeline
end

"""
"""
function singlelabelbinaryclassdataframelogisticclassifier(
        feature_names::AbstractVector,
        single_label_name::Symbol,
        single_label_levels::AbstractVector;
        package::Symbol = :none,
        intercept::Bool = true,
        interactions::Integer = 1,
        name::AbstractString = "",
        )
    if package == :GLM
        result = singlelabelbinaryclassdataframelogisticclassifier_GLM(
            feature_names,
            single_label_name,
            single_label_levels;
            intercept = intercept,
            interactions = interactions,
            name = name,
            )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end

"""
"""
# function singlelabelbinaryclassdataframeprobitclassifier_GLM(
#         feature_names::AbstractVector,
#         single_label_name::Symbol,
#         single_label_levels::AbstractVector;
#         intercept::Bool = true,
#         interactions::Integer = 1,
#         name::AbstractString = "",
#         )
#     negative_class = single_label_levels[1]
#     positive_class = single_label_levels[2]
#     formula = generate_formula(
#         [single_label_name],
#         feature_names;
#         intercept = intercept,
#         interactions = interactions,
#         )
#     dftransformer =
#         ImmutableDataFrame2GLMSingleLabelBinaryClassTransformer(
#             single_label_name,
#             positive_class,
#             )
#     glmestimator = GLMModel(
#         formula,
#         GLM.Binomial(),
#         GLM.ProbitLink();
#         isclassificationmodel = true,
#         isregressionmodel = false,
#         )
#     predictlabelfixer =
#         ImmutablePredictionsSingleLabelInt2StringTransformer(
#             0,
#             single_label_levels,
#             )
#     predprobalabelfixer =
#         ImmutablePredictProbaSingleLabelInt2StringTransformer(
#             0,
#             single_label_levels,
#             )
#     probapackager = ImmutablePackageSingleLabelPredictProbaTransformer(
#         single_label_name,
#         )
#     finalpipeline = dftransformer |>
#                     glmestimator |>
#                     predictlabelfixer |>
#                     predprobalabelfixer |>
#                     probapackager
#     finalpipeline.name = name
#     return finalpipeline
# end

"""
"""
# function singlelabelbinaryclassdataframeprobitclassifier(
#         feature_names::AbstractVector,
#         single_label_name::Symbol,
#         single_label_levels::AbstractVector;
#         package::Symbol = :none,
#         intercept::Bool = true,
#         interactions::Integer = 1,
#         name::AbstractString = "",
#         )
#     if package == :GLM
#         result = singlelabelbinaryclassdataframeprobitclassifier_GLM(
#             feature_names,
#             single_label_name,
#             single_label_levels;
#             intercept = intercept,
#             interactions = interactions,
#             name = name,
#             )
#         return result
#     else
#         error("$(package) is not a valid value for package")
#     end
# end

"""
"""
function single_labeldataframelinearregression_GLM(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        intercept::Bool = true,
        interactions::Integer = 1,
        name::AbstractString = "",
        )
    formula = generate_formula(
        [single_label_name],
        feature_names;
        intercept = intercept,
        interactions = interactions,
        )
    glmestimator = GLMModel(
        formula,
        GLM.Normal(),
        GLM.IdentityLink();
        isclassificationmodel = false,
        isregressionmodel = true,
        )
    finalpipeline = SimplePipeline(
        AbstractFittable[
            glmestimator,
            ];
        name = name,
        )
    return finalpipeline
end

"""
"""
function single_labeldataframelinearregression(
        feature_names::AbstractVector,
        single_label_name::Symbol;
        package::Symbol = :none,
        intercept::Bool = true,
        interactions::Integer = 1,
        name::AbstractString = "",
        )
    if package == :GLM
        result = single_labeldataframelinearregression_GLM(
            feature_names,
            single_label_name;
            intercept = intercept,
            interactions = interactions,
            name = name,
            )
        return result
    else
        error("$(package) is not a valid value for package")
    end
end


text:
using SelfPropelledModel
using Test
using SPTrackingToolkit
using DelimitedFiles

mean(X)=foldl((x,y)->(x[1]+1,(x[2]*x[1]+y)/(x[1]+1)),X,init=(0,zero(first(X))/1))[2]

function simulation(Np,Nf)
    sv=sample_particles(Np,()->500.0.*(rand(2).-0.5),()->(rand(2).-0.5))
    frames=[]

    for i in 1:Nf
        X=foldl(hcat,[[p.x...] for p in sv])
        push!(frames,X)
        spm_step!(sv,0.4,0.2,0.2,1.5)
    end

    specs = SPT(maxtimegap=5.0,maxdist=5.0,dims=2,rest=0,verbose=true)

    tr=track(specs,frames)

    errors=zeros(length(tr))
    Threads.@threads for j in eachindex(tr)
        trace = tr[j]
        _ , Ntf=size(trace)
        op=nothing
        toterrors=-1
        for i in 1:Ntf
            t=trunc(Int,trace[1,i])
            X=trace[2:3,i]
            p=findfirst(isequal(X),eachcol(frames[t])|>collect)
            p==op || (toterrors+=1)
            op=p
        end
        errors[j]=toterrors/Ntf*100
    end

    (errors .> 0)|>sum,mean(errors)
end

function tohdsim(Np,Nf)
    R=[]
    for i in 1:300
        push!(R,simulation(Np,Nf))
    end
    AR=foldl(hcat,[x...] for x in R)'
    writedlm("results/P$(Np)_F$(Nf).dat",AR)
end
Ns=[20,40,60,80,100,120,140,160]
for p in Ns
    tohdsim(p,170)
end

using DelimitedFiles
using PyPlot
using PyCall
using Statistics
cd(dirname(@__FILE__));
include("mplstyle.jl")
cd("../results")

pygui(true)
matplotlibstyle()
Ns=20:20:160
A=[(N,readdlm("P$(N)_F170.dat")) for N in Ns]

avgA=foldl(hcat,[N,Statistics.mean(x,dims=1)...,std(x,dims=1)...] for (N,x) in A)

py"setfonts(12)"
figure(figsize=py"mplfigsize(1)".*1.8,dpi=300)
errorbar(avgA[1,:],avgA[2,:],yerr=avgA[4,:]/10,fmt="o",lw=2,color="royalblue",
    markersize=7,markeredgewidth=2,marker="o",markerfacecolor="None",
    markeredgecolor="royalblue")
xticks(Ns)
xlabel(L"N^\circ \rm{particles}")
ylabel(L"N^\circ \rm{erroneous\,\,tracks}")
tight_layout()
savefig("n_errors_vs_n_particles.pdf")

figure(figsize=py"mplfigsize(1)".*1.8,dpi=300)
errorbar(avgA[1,:],avgA[3,:],yerr=avgA[5,:]/10,fmt="o",lw=2,color="royalblue",
    markersize=7,markeredgewidth=2,marker="o",markerfacecolor="None",
    markeredgecolor="royalblue")
xticks(Ns)
xlabel(L"N^\circ \rm{particles}")
ylabel(L"\rm{Percentual\,\,errors\,\,}\%")
tight_layout()
savefig("rel_errors_vs_n_particles.pdf")

using LsqFit
using LinearAlgebra

wt = inv.(avgA[4,:]/10).^2
m(n, p) = @. p[1]*p[1]*n*n + p[2]*p[1]*n + p[3]
fit = curve_fit(m, avgA[1,:],avgA[2,:], wt, [0.01,0.01,0.01])
cov = estimate_covar(fit)
P=fit.param

diag(cov).^0.5

py"setfonts(12)"
figure(figsize=py"mplfigsize(1.4)".*1.8,dpi=300)
subplot(2,1,1)
errorbar(avgA[1,:],avgA[2,:],yerr=avgA[4,:]/10,fmt="o",lw=1,color="C0",
    markersize=6,markeredgewidth=1.333,marker="o",markerfacecolor="None",
    markeredgecolor="C0")
plot(0:180,m(0:180,P),lw=1.333)
xticks(Ns)
xlim(10,170)
xlabel(L"N^\circ \rm{particles}")
ylabel(L"N^\circ \rm{erroneous\,\,tracks}")

χ²=sum(abs2,(avgA[2,:]-m(avgA[1,:],P))./(avgA[4,:]/10))/(length(avgA[2,:]) -3 )

text(20, 9, "χ² = $(round(χ²*100)/100)", fontsize=12)

subplot(2,1,2)
scatter(avgA[1,:],(avgA[2,:]-m(avgA[1,:],P))./(avgA[4,:]/10))
xlabel(L"N^\circ \rm{particles}")
ylabel(L"\rm{Std. Residuals}")
xticks(Ns)
xlim(10,170)
ylim(-1,1)
tight_layout(rect=[0,0,0.95,0.95])
savefig("n_errors_vs_n_particles.pdf")

wt = inv.(avgA[5,:]/10).^2
m(n, p) = @. p[1]*p[1]*n*n + p[2]*p[1]*n + p[3]
fit = curve_fit(m, avgA[1,:],avgA[3,:], wt, [0.01,0.01,0.01])
cov = estimate_covar(fit)
P=fit.param

py"setfonts(12)"
figure(figsize=py"mplfigsize(1.4)".*1.8,dpi=300)
subplot(2,1,1)
errorbar(avgA[1,:],avgA[3,:],yerr=avgA[5,:]/10,fmt="o",lw=1,color="C0",
    markersize=6,markeredgewidth=1.333,marker="o",markerfacecolor="None",
    markeredgecolor="C0")
plot(0:180,m(0:180,P),lw=1.333)
xticks(Ns)
xlim(10,170)
xlabel(L"N^\circ \rm{particles}")
ylabel(L"\rm{Percentual\,\,errors\,\,}\%")

χ²=sum(abs2,(avgA[3,:]-m(avgA[1,:],P))./(avgA[5,:]/10))/(length(avgA[2,:]) -3 )

text(20, 0.04, "χ² = $(round(χ²*100)/100)", fontsize=12)

subplot(2,1,2)
scatter(avgA[1,:],(avgA[3,:]-m(avgA[1,:],P))./(avgA[5,:]/10))
xlabel(L"N^\circ \rm{particles}")
ylabel(L"\rm{Std. Residuals}")
xticks(Ns)
xlim(10,170)
ylim(-1,1)
tight_layout(rect=[0,0,0.95,0.95])

tight_layout()
savefig("rel_errors_vs_n_particles.pdf")

py"setfonts(12)"
figure(figsize=py"mplfigsize(1)".*1.4,dpi=300)

A=randn(10000).*2.2 .+ 16.3
B=randn(10000).*0.6 .+ 19.5
AB=[A;B]

samples=(2,2)
dhp=[abs(mean(rand(A,samples[1]))-mean(rand(B,samples[2]))) for i in 1:120000]
dnull=[abs(mean(rand(AB,samples[1]))-mean(rand(AB,samples[2]))) for i in 1:120000]
using Statistics
qnull=quantile(dnull,0.95)
V=[dhp dnull]

hist(V,60,density=true, histtype="step", stacked=false,label=["Truth", "Null"], fill=false,color=["C3","C0"])
vlines(qnull,0,0.29,linestyle="dashed")
text(qnull, 0.3, "95%", fontsize=10, horizontalalignment="center")
text(mean(B)-mean(A), 0.3, "obs", fontsize=10, horizontalalignment="center")
vlines(mean(B)-mean(A),0,0.29)
legend()
xlabel(L"ΔT = |Δτ_{mcf7} - Δτ_{mda}|")
ylabel(L"p.d.f.")
xlim(0,10)
savefig("timecomp.pdf")

sv=sample_particles(160,()->500.0.*(rand(2).-0.5),()->(rand(2).-0.5))
frames=[]

for i in 1:170
    X=foldl(hcat,[[p.x...] for p in sv])
    push!(frames,X)
    spm_step!(sv,0.4,0.2,0.2,1.5)
end

pygui(true)
figure(figsize=(6,6))
for f in frames
    scatter(f[1,:],f[2,:],color="black",s=0.125)
end


text:
module SpectralMeasures
using Base, LinearAlgebra, ApproxFun, RecipesBase, RatFun, BandedMatrices, BlockArrays,
        InfiniteArrays

import Base: +, -, *, /, getindex, inv
import LinearAlgebra: adjoint, transpose, issymmetric, SymTridiagonal, eigvals, eigen

import ApproxFun: Operator, ToeplitzOperator, DiracSpace, IdentityOperator,
            TridiagonalOperator, setdomain, resizedata!, bandwidths, PointSpace,
            BandedMatrix, TimesOperator, SpaceOperator,
            MatrixSpace, ∞, ℓ⁰, domainspace, rangespace, domain, mul_coefficients,
            ldiv_coefficients, InterlaceOperator

import BlockArrays: nblocks

import InfiniteArrays: AbstractInfUnitRange

export spectralmeasure, discreteeigs, principalresolvent, discresolvent, validatedspectrum, spectrum
export connectioncoeffsoperator, applyconversion, SymTriOperator, SymTriPertToeplitz
export connectioncoeffsmatrix
export tripleplot
export freejacobioperator, jacobioperator, ql

include("HessenbergUnitary.jl")
include("PertToeplitz.jl")
include("helper.jl")
include("ql.jl")

function spectralmeasure(a,b)
    TT = promote_type(eltype(a),eltype(b))
  # Chop the a and b down
  a = chop!(a); b = 0.5 .+ chop!(b .- 0.5)
  n = max(2,length(a),length(b)+1)
  a = [a;zeros(TT,n-length(a))]; b = [b; 0.5 .+ zeros(TT,n-length(b))]

  # Finds C such that J*C = C*Toeplitz([0,1/2])
  C = SpectralMeasures.connectioncoeffsoperator(a,b)
  c = Fun(Taylor,C.T.nonnegative)
  f = Fun(C*(C'*[1]),Ultraspherical(1))

  # Check for discrete eigenvalues
  z = sort(real(filter!(z->abs(z)<1 && abs(imag(z)) ≤ 10eps(TT)  && !isapprox(abs(z),1),
                        complexroots(c))))
  if length(z) > 0
     Cmu = SpectralMeasures.connectioncoeffsoperator(a[2:end],b[2:end]) # Technically not Cmu from the paper
     cmu = Fun(Taylor,[0;Cmu.T.nonnegative]/b[1]) # this is cmu from the paper
     cprime = differentiate(c)
     eigs=real(map(joukowsky,z))
     weights = (1 .- 1 ./ z.^2).*(real.(cmu.(z))./(2 .* real.(cprime.(z))))
     p = Fun(DiracSpace(eigs),weights) + Fun(JacobiWeight(0.5,0.5,Ultraspherical(1)),[2/TT(pi)])
     q = Fun(PointSpace(eigs),ones(TT,length(eigs))) + f
     μ = RationalFun(p,q)
  else
    μ = RationalFun(Fun(JacobiWeight(0.5,0.5,Ultraspherical(1)),[2/TT(pi)]),f)
  end
  μ
end

function principalresolvent(a,b)
  # Chop the a and b down
  a = chop!(a); b = .5+chop!(b-.5)
  n = max(2,length(a),length(b)+1)
  a = [a;zeros(n-length(a))]; b = [b;.5+zeros(n-length(b))]

  # Compute the necessary polynomials
  C = connectioncoeffsoperator(a,b)
  Cmu = connectioncoeffsoperator(a[2:end],b[2:end]) # Technically not Cmu from the paper
  f = Fun((C*(C'*[1])),Ultraspherical(1))
  fmu = Fun(Ultraspherical(1),coefficients(Cmu*((C'*[1]).coefficients[2:end])/b[1]))

  # Return the resolvent
  λ->(2*sqrt(complex(λ-1)).*sqrt(complex(λ+1))-2*λ-extrapolate(fmu,λ))./extrapolate(f,λ)
end

function discresolvent(a,b)
  # Chop the a and b down
  a = chop!(a); b = .5+chop!(b-.5)
  n = max(2,length(a),length(b)+1)
  a = [a;zeros(n-length(a))]; b = [b;.5+zeros(n-length(b))]

  # Compute the necessary polynomials
  C = SpectralMeasures.connectioncoeffsoperator(a,b)
  Cmu = SpectralMeasures.connectioncoeffsoperator(a[2:end],b[2:end]) # Technically not Cmu from the paper
  c = Fun(Taylor,C.T.nonnegative)
  cmu = Fun(Taylor,[0;Cmu.T.nonnegative]/b[1]) # this is the cmu from the paper

  # Return the rational function
  z->-cmu(z)./c(z)
end

function discreteeigs(a,b)
  a = chop!(a); b = 0.5 .+ chop!(b .- 0.5)
  n = max(2,length(a),length(b)+1)
  a = [a;zeros(n-length(a))]; b = [b;.5 .+ zeros(n-length(b))]
  # Finds C such that C*J = Toeplitz([0,1/2])*C
  C = connectioncoeffsoperator(a,b)
  Tfun = Fun(Taylor,C.T.nonnegative)
  sort(real(map(joukowsky,filter!(z->abs(z)<1 && isreal(z) && !isapprox(abs(z),1),complexroots(Tfun)))))
end

#Finds C such that C'(U_k(s)) =  (P_k(s)),
# where P_k has Jacobi coeffs a,b and U_k is Chebyshev U
function connectioncoeffsoperator(a,b)
  n = max(2,length(a),length(b)+1)
  N = 2*n #This is sufficient only because we go from Chebyshev U
  a = [a;zeros(N-length(a))]; b = [b;.5 .+ zeros(N-length(b))]

  elType = eltype(a)
  ToeplitzVec = zeros(elType,N)
  K = BandedMatrix(Zeros{elType}(n,N), (0,N+1))
  K[1,1] = 1
  K[1,2] = -a[1]/b[1]
  K[2,2] = .5/b[1]
  # The recurrence for the first n+1 cols depend on a and b
  for j = 3:n+1
    K[1,j] = (-a[j-1]*K[1,j-1] + .5*K[2,j-1] - b[j-2]*K[1,j-2])/b[j-1]
    for i = 2:j-2
      K[i,j] = (.5*K[i-1,j-1] -a[j-1]*K[i,j-1] + .5*K[i+1,j-1] - b[j-2]*K[i,j-2])/b[j-1]
    end
    K[j-1,j] = (.5*K[j-2,j-1] -a[j-1]*K[j-1,j-1] - b[j-2]*K[j-1,j-2])/b[j-1]
    if j<n+1
      K[j,j] = .5*K[j-1,j-1]/b[j-1]
    end
  end
  ToeplitzVec[1] = K[n,n]
  ToeplitzVec[2] = K[n,n+1]
  # The recurrence for rows n+2 to 2n do not depend on a and b
  for j = n+2:N
    K[1,j] = K[2,j-1] - K[1,j-2]
    for i = 2:N-j
      K[i,j] = K[i-1,j-1] + K[i+1,j-1] - K[i,j-2]
    end
    if j < N
      K[N+1-j,j] = K[N-j,j-1] + K[N+2-j,j-1] - K[N+1-j,j-2]
    end
    ToeplitzVec[2*(j-n)-1] = K[N+1-j,j-1]
    ToeplitzVec[2*(j-n)] = K[N+1-j,j]
  end
  T = ToeplitzOperator(elType[],chop!(ToeplitzVec))
  for j = 1:N
    for i = 1:min(j,N+1-j)
      K[i,j]-=T[i,j]
    end
  end
  T+FiniteOperator(K,ℓ⁰,ℓ⁰)
end
end  #Module


text:
@testset "trajectories" begin

    @testset "VectorialTrajectory" begin
        b = VectorialTrajectory(; state = Vector{Int}, reward = Float64)

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == true

        t1 = (state = [1, 2], reward = 0.0)
        push!(b; t1...)

        @test length(b) == 1
        @test size(b) == (1,)
        @test b[1] == b[end] == t1
        @test isempty(b) == false
        @test get_trace(b, :state) == [t1.state]
        @test get_trace(b, :reward) == [t1.reward]

        t2 = (state = [3, 4], reward = -1)
        push!(b; t2...)

        @test length(b) == 2
        @test size(b) == (2,)
        @test b[2] == b[end] == t2
        @test isempty(b) == false
        @test get_trace(b, :state) == [t1.state, t2.state]
        @test get_trace(b, :reward) == [t1.reward, t2.reward]

        empty!(b)

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == true
    end

    @testset "VectorialCompactSARTSATrajectory" begin
        b = VectorialCompactSARTSATrajectory()

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == true

        t1 = (state = 1, action = 2)
        push!(b; t1...)
        t2 = (reward = 1.0, terminal = false, state = 2, action = 3)
        push!(b; t2...)

        @test length(b) == 1
    end

    @testset "EpisodicCompactSARTSATrajectory" begin
        b = EpisodicCompactSARTSATrajectory()

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == true

        t1 = (state = 1, action = 2)
        push!(b; t1...)

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1]
        @test get_trace(b, :action) == [2]
        @test get_trace(b, :reward) == []
        @test get_trace(b, :terminal) == []
        @test get_trace(b, :next_state) == []
        @test get_trace(b, :next_action) == []

        t2 = (reward = 1.0, terminal = false, state = 2, action = 3)
        push!(b; t2...)

        @test length(b) == 1
        @test size(b) == (1,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1]
        @test get_trace(b, :action) == [2]
        @test get_trace(b, :reward) == [1.0]
        @test get_trace(b, :terminal) == [false]
        @test get_trace(b, :next_state) == [2]
        @test get_trace(b, :next_action) == [3]
        @test b[1] ==
              b[end] ==
              (
                  state = 1,
                  action = 2,
                  reward = 1.0f0,
                  terminal = false,
                  next_state = 2,
                  next_action = 3,
              )

        t3 = (reward = 2.0, terminal = true, state = 3, action = 4)
        push!(b; t3...)

        @test length(b) == 2
        @test size(b) == (2,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1, 2]
        @test get_trace(b, :action) == [2, 3]
        @test get_trace(b, :reward) == [1.0, 2.0]
        @test get_trace(b, :terminal) == [false, true]
        @test get_trace(b, :next_state) == [2, 3]
        @test get_trace(b, :next_action) == [3, 4]
        @test b[end] == (
            state = 2,
            action = 3,
            reward = 2.0f0,
            terminal = true,
            next_state = 3,
            next_action = 4,
        )

        pop!(b)

        @test length(b) == 1
        @test size(b) == (1,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1]
        @test get_trace(b, :action) == [2]
        @test get_trace(b, :reward) == [1.0]
        @test get_trace(b, :terminal) == [false]
        @test get_trace(b, :next_state) == [2]
        @test get_trace(b, :next_action) == [3]
    end

    @testset "CircularCompactSARTSATrajectory" begin
        b = CircularCompactSARTSATrajectory(; capacity = 3)

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == true

        t1 = (state = 1, action = 2)
        push!(b; t1...)

        @test length(b) == 0
        @test size(b) == (0,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1]
        @test get_trace(b, :action) == [2]
        @test get_trace(b, :reward) == []
        @test get_trace(b, :terminal) == []
        @test get_trace(b, :next_state) == []
        @test get_trace(b, :next_action) == []

        t2 = (reward = 1.0, terminal = false, state = 2, action = 3)
        push!(b; t2...)

        @test length(b) == 1
        @test size(b) == (1,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1]
        @test get_trace(b, :action) == [2]
        @test get_trace(b, :reward) == [1.0]
        @test get_trace(b, :terminal) == [false]
        @test get_trace(b, :next_state) == [2]
        @test get_trace(b, :next_action) == [3]
        @test b[1] ==
              b[end] ==
              (
                  state = 1,
                  action = 2,
                  reward = 1.0f0,
                  terminal = false,
                  next_state = 2,
                  next_action = 3,
              )

        t3 = (reward = 2.0, terminal = true, state = 3, action = 4)
        push!(b; t3...)

        @test length(b) == 2
        @test size(b) == (2,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1, 2]
        @test get_trace(b, :action) == [2, 3]
        @test get_trace(b, :reward) == [1.0, 2.0]
        @test get_trace(b, :terminal) == [false, true]
        @test get_trace(b, :next_state) == [2, 3]
        @test get_trace(b, :next_action) == [3, 4]
        @test b[end] == (
            state = 2,
            action = 3,
            reward = 2.0f0,
            terminal = true,
            next_state = 3,
            next_action = 4,
        )

        pop!(b, :state, :action)
        push!(b, state = 4, action = 5)

        @test b[end] == (
            state = 2,
            action = 3,
            reward = 2.0f0,
            terminal = true,
            next_state = 4,
            next_action = 5,
        )

        pop!(b)

        @test b[1] ==
              b[end] ==
              (
                  state = 1,
                  action = 2,
                  reward = 1.0f0,
                  terminal = false,
                  next_state = 2,
                  next_action = 3,
              )
        @test length(b) == 1

        t4 = (reward = 3.0, terminal = false, state = 4, action = 5)
        push!(b; t4...)

        @test length(b) == 2
        @test size(b) == (2,)
        @test isempty(b) == false
        @test get_trace(b, :state) == [1, 2]
        @test get_trace(b, :action) == [2, 3]
        @test get_trace(b, :reward) == [1.0, 3.0]
        @test get_trace(b, :terminal) == [false, false]
        @test get_trace(b, :next_state) == [2, 4]
        @test get_trace(b, :next_action) == [3, 5]
        @test b[end] == (
            state = 2,
            action = 3,
            reward = 3.0f0,
            terminal = false,
            next_state = 4,
            next_action = 5,
        )
    end

end


text:
using LinearAlgebra
using SparseArrays
using Plots

"This routine solves the advection equation with an inflow boundary condition using
Forward Euler in time and finite differences in space."

m = 100 # number of points

# define spatial grid
x = LinRange(-1,1,m+2)
xint = x[2:end]
h = x[2]-x[1]

α(t) = sin(pi*t)
dt = h # timestep
T = 4.0 # final time

# initial condition and forcing
u0(x) = 0.0
f(x,t) = 0.0
function F(α,xint,t,h)
    fvec = f.(xint,t)
    fvec[1] += α(t)/h
    return fvec
end

# upwind approximation of first derivative
Q = (1/h) * spdiagm(0=>ones(m+1),-1=>-ones(m))

u = u0.(xint)
Nsteps = ceil(Int,T/dt)
dt = T / Nsteps

interval = 1
plot()
@gif for k = 1:Nsteps
    tk = k*dt
    u .= u + dt * (F(α,xint,tk,h) - Q*u)
    if k % interval==0
        plot(xint,u,linewidth=2,label="Solution")
        println("on timestep $k out of $Nsteps.")
    end
end every interval

# plot(xint,u,linewidth=2,label="Solution")
# plot!(x,u0.(x .- T),linewidth=2,label="Exact solution",ylims=(-1.0,2.0))


text:
using Plots
using ImplicitEquations

f(x,y) = y^2 - x^3 + x

p1 = plot(Eq(f,0))

xs = -1:0.01:4
ys = -4:0.01:4
zs = f.(xs',ys)

@show size(zs)
p1 = plot(Eq(f,0))
p2 = plot(xs,ys,zs,st=:surface,camera=(10,30))

lay = @layout [a b]
plot(p1,p2,layout=lay,xlabel="x",ylabel="y")

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function NamedTupleFromStruct(str)
    nt_struct = NamedTupleTools.ntfromstruct(str)
    nt_type = (type = typeof(str))
    return (nt_type = nt_type, nt_struct = nt_struct)
end

function StructFromNamedTuple(nt::NamedTuple)
    return NamedTupleTools.structfromnt(nt.nt_type, nt.nt_struct)
end

function NamedTuple(::Missing)
    return (object = missing,)
end

function NamedTuple(d::Dict)
    return (dict_json_string = json(d),)
end
Base.convert(T::Type{NamedTuple}, x::Dict) = T(x)

function Dict(nt::NamedTuple)
    JSON.parse(nt.dict_json_string)
end
Base.convert(T::Type{Dict}, x::NamedTuple) = T(x)


text:
import Pkg
Pkg.activate("..")

using Distributions: Beta, Bernoulli, Normal, LogNormal, truncated
using CrimsonSkyline
using Random: seed!
using Logging
using StatsBase: mean

seed!(2021)

function coin_model(t::Trace, data::Vector{Bool})
    bias = sample(t, "bias", Beta(3.0, 3.0))
    coin = Bernoulli(bias)
    for (i, d) in enumerate(data)
        observe(t, "flip $i", coin, d)
    end
end

function plated_coin_model(t::Trace, data::Vector{Bool})
    bias = sample(t, "bias", Beta(3.0, 3.0))
    plate(t, observe, "flips", Bernoulli(bias), data)
end

function normal_model(t::Trace, data::Vector{Float64})
    loc = sample(t, "loc", Normal())
    scale = sample(t, "scale", LogNormal())
    for (i, d) in enumerate(data)
        observe(t, "obs $i", Normal(loc, scale), d)
    end
end

function plated_normal_model(t::Trace, data::Vector{Float64})
    loc = sample(t, "loc", Normal())
    scale = sample(t, "scale", LogNormal())
    plate(t, observe, "obs", Normal(loc, scale), data)
end

loc_proposal(t0, t1, params...) = propose(t1, "loc", Normal(t0["loc"].value, 0.25))
scale_proposal(t0, t1, params...) = propose(t1, "scale", truncated(Normal(t0["scale"].value, 0.25), 0.01, Inf))

function slow()
    @warn "\n~~~Coin model: not using plate~~~"
    data = [true, true, false, false, true, true, true, true]
    _ = mh(coin_model; params = (data,))  # to compile
    @time results = mh(coin_model; params = (data,))
    @info "Mean posterior bias estimate: $(mean(results, "bias"))\n"

    @warn "\n~~~Normal model, not using plate or typed trace"
    seed!(2021)
    data = randn(10) .* 2.0 .+ 2.0
    kernels = [loc_proposal, scale_proposal]
    _ = mh(normal_model, kernels; params = (data,), burn=1000, thin=50, num_iterations=1, inverse_verbosity=Inf)  # to compile
    @time results = mh(normal_model, kernels; params = (data,), burn=1000, thin=50, num_iterations=10000, inverse_verbosity=Inf)
end

function fast()
    @warn "\n~~~Coin model: using plate~~~\n"
    data = [true, true, false, false, true, true, true, true]
    _ = mh(plated_coin_model; params = (data,))  # to compile
    @time results = mh(plated_coin_model; params = (data,))
    @info "Mean posterior bias estimate: $(mean(results, "bias"))\n"

    @warn "\n~~~Normal model, using typed trace~~~"
    seed!(2021)
    data = randn(10) .* 2.0 .+ 2.0
    kernels = [loc_proposal, scale_proposal]
    types = (String, Float64)
    _ = mh(normal_model, kernels, types; params = (data,), burn=1000, thin=50, num_iterations=1, inverse_verbosity=Inf)  # to compile
    @time results = mh(normal_model, kernels, types; params = (data,), burn=1000, thin=50, num_iterations=10000, inverse_verbosity=Inf)

    @warn "\n~~~Normal model, using plate~~~"
    seed!(2021)
    data = randn(10) .* 2.0 .+ 2.0
    kernels = [loc_proposal, scale_proposal]
    _ = mh(plated_normal_model, kernels; params = (data,), burn=1000, thin=50, num_iterations=1, inverse_verbosity=Inf)  # to compile
    @time results = mh(plated_normal_model, kernels; params = (data,), burn=1000, thin=50, num_iterations=10000, inverse_verbosity=Inf)
end

function main()
    slow()
    fast()
end

main()


text:
  #
  # Here is clonning algorithm. It was created as a replacement of embedded
  # deepcopy() function
  # TODO: describe my copy() version
  # TODO: describe meta map structure
  # TODO: describe mem addresses
  # TODO: describe requirement, that meta info and expr tree should be well prepared
  #
  include("..\\src\\ImportFolders.jl")
  import Manager
  import Creature.Block
  import Creature.Func
  #
  # Type, which describes a function in original Expr tree. It contains
  #
  type MetaFunc
    expr::Expr
    blocks::Dict{Int, Expr}
  end
  type MetaCode
    curFn::Int
    funcs::Dict{Int, MetaFunc}
  end
  #
  # Makes a clone of "e" expression and updates meta information
  # in "funcs" Array
  #
  function clone(e::Expr, funcs::Array{Func, 1})
    local f::Func
    local b::Block
    local mf::MetaFunc
    local h::Int = Int(pointer_from_objref(e))
    local mOrg::MetaCode = MetaCode(h,Dict{Int, MetaFunc}(h => MetaFunc(e, Dict{Int, Expr}(Int(pointer_from_objref(e.args[2])) => e))))
    local ex::Expr = _copy(e, mOrg)
    #
    # This is how we update meta information
    #
    for f in funcs
      mf     = mOrg.funcs[Int(pointer_from_objref(f.code))]
      f.code = mf.expr
      for b in f.blocks b.expr = mf.blocks[Int(pointer_from_objref(b.expr))] end
    end
    #
    # References main function and it's block should be set separately
    # because our recursive code doesn't work well with main (first) function
    #
    funcs[1].code = ex                    # Main function
    funcs[1].blocks[1].expr = ex.args[2]  # Main function block

    ex
  end
  function _copy(e::Expr, mOrg::MetaCode)
    local n::Expr = Expr(e.head);
    n.args = _copyExprArgs(e.args, mOrg);
    n.typ = e.typ;
    n
  end
  function _copyExprArgs(arr::Array{Any,1}, mOrg::MetaCode)
    local a::Any
    local arrCopy::Array{Any, 1} = []
    for a in arr push!(arrCopy, _copyExprs(a, mOrg)) end
    arrCopy
  end
  function _copyExprs(a::ANY, mOrg::MetaCode) a end
  #
  # Don't change the order of these two if blocks
  #
  function _copyExprs(e::Expr, mOrg::MetaCode)
    local fn::MetaFunc
    local blocks::Dict{Int, Expr}
    local blockId::Int
    #
    # Don't touch this code. It's optimized for speed
    #
    if e.head === :function
      push!(mOrg.funcs, (mOrg.curFn = Int(pointer_from_objref(e))) => (fn = MetaFunc(e, Dict{Int, Expr}())))
      return fn.expr = _copy(e, mOrg)
    elseif e.head === :block
      blocks  = mOrg.funcs[mOrg.curFn].blocks
      blockId = Int(pointer_from_objref(e))
      push!(blocks, blockId => e)
      return blocks[blockId] = _copy(e, mOrg)
    end

    _copy(e, mOrg)
  end
  #
  # entry point
  #
  cfg = Config.create()
  c   = Creature.create(cfg)
  for i = 1:1000 Mutator.mutate(cfg, c) end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
"""
    abstract type AbstractMesh{N,T}

An abstract mesh structure in dimension `N` with primite data of type `T` (e.g.
`Float64` for double precision representation).

Concrete subtypes of `AbstractMesh` should implement [`ElementIterator`](@ref)
for accessing the mesh elements, and/or [`NodeIterator`](@ref) for accesing the
mesh nodes.

# See also: [`GenericMesh`](@ref), [`UniformCartesianMesh`](@ref)
"""
abstract type AbstractMesh{N,T} end

ambient_dimension(M::AbstractMesh{N}) where {N} = N

# we define the geometric dimension of the parent to be the largest of the geometric
# dimension of its entities.
geometric_dimension(M::AbstractMesh) = maximum(x -> geometric_dimension(x), entities(M))

primitive_type(M::AbstractMesh{N,T}) where {N,T} = T

function Base.show(io::IO,msh::AbstractMesh)
    print(io,"$(typeof(msh)) containing:")
    for E in keys(msh)
        iter = ElementIterator(msh,E)
        print(io,"\n\t $(length(iter)) elements of type ",E,)
    end
    return io
end

"""
    struct ElementIterator{E,M}

Return an iterator for all elements of type `E` on a mesh of type `M`.

Besides the methods listed in the [iterator
iterface](https://docs.julialang.org/en/v1/manual/interfaces/) of `Julia`, some
functions also require the `getindex(iter,i::Int)` method for accessing the
`i`-th element directly.
"""
struct ElementIterator{E,M}
    mesh::M
end

mesh(iter::ElementIterator) = iter.mesh

Base.eltype(::SType{ElementIterator{E}}) where {E} = E

ElementIterator{E}(parent::M) where {E,M <: AbstractMesh} = ElementIterator{E,M}(parent)

ElementIterator(parent,E) = ElementIterator{E}(parent)

# indexing a mesh with an element type return an interator over that element
# (like a dict)
Base.getindex(m::AbstractMesh,E::DataType) = ElementIterator(m, E)

function Base.length(iter::ElementIterator)
    prod(size(iter))
end

"""
    struct NodeIterator{M}

Iterator for all the nodes in a mesh of type `M`.
"""
struct NodeIterator{M}
    mesh::M
end

mesh(iter::NodeIterator) = iter.mesh

# nodes of `AbstractMesh{N,T}` should be points of type `SVector{N,T}`
Base.eltype(::SType{NodeIterator{<:AbstractMesh{N,T}}}) where {N,T} = SVector{N,T}

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# m5.5s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

# ### snippet 5.29

println()
df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
df = filter(row -> !(row[:neocortex_perc] == "NA"), df);
df[!, :neocortex_perc] = parse.(Float64, df[:, :neocortex_perc])
df[!, :lmass] = log.(df[:, :mass])
scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])

stan5_5 = "
data {
    int < lower = 1 > N; // Sample size
    vector[N] K; // Outcome
    vector[N] NC; // Predictor
}

parameters {
    real a; // Intercept
    real bN; // Slope (regression coefficients)
    real < lower = 0 > sigma;    // Error SD
}

model {
    vector[N] mu;               // mu is a vector
    a ~ normal(0, 0.2);           //Priors
    bN ~ normal(0, 0.5);
    sigma ~ exponential(1);
    mu = a + bN * NC;
    K ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m5_5s = SampleModel("m5.5", stan5_5);

# Input data for cmdstan

m5_5_data = Dict("N" => size(df, 1), "NC" => df[!, :neocortex_perc_s],
    "K" => df[!, :kcal_per_g_s]);

# Sample using StanSample

rc5_5s = stan_sample(m5_5s, data=m5_5_data);

if success(rc5_5s)

    # Describe the draws

    part5_5s = read_samples(m5_5s, :particles)
    part5_5s |> display

    retinking = "
            mean   sd  5.5% 94.5%
    a     0.04 0.15 -0.21  0.29
    bN    0.13 0.22 -0.22  0.49
    sigma 1.00 0.16  0.74  1.26
    "

end


text:
# #############################################################################
# #############################################################################
#Helpful functions for computation

#Small nonnegative integer powers implemented in GSL
function gsl_pow_int(x::Float64,n::Int64)
    return ccall((:gsl_pow_int,"libgsl"),Cdouble,(Cdouble,Cint),x,convert(Int32,n))
end

#Squares and cubes are more common, simply implemented directly in Julia with shorter names, and also vectorized
sq{T<:Number}(x::T) = x*x::T
cb{T<:Number}(x::T) = x*x*x::T
@vectorize_1arg Number sq
@vectorize_1arg Number cb

#Custom Base64 representation of a number
const base64digits = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ+-"
function dec2base64(n::Integer,pad::Int64=0)
    return lpad(base64digits[reverse(digits(n,64)+1)],pad,'0')
end

#Generic hashing of composite types based on typename, fields and values
function fieldhashing(obj)
    #Initiate hash with typename
    h::UInt64 = hash(summary(obj))
    for field in fieldnames(obj)
        #Use xor to combine hashes from all fields so that order does not matter
        #Hash the field name and value together as a tuple
        h $= hash((field,obj.(field)))
    end
    return h
end

#Generic equality of composite types based on typename, fields and values
function fieldequality(o1,o2)
    if typeof(o1) != typeof(o2)
        return false
    end
    for field in fieldnames(o1)
        if o1.(field) != o2.(field)
            return false
        end
    end
    return true
end

#Shrinking a matrix by deleting the last few rows or columns as much as possible without losing any nonzero value. Useful for matrices of coefficients, where position is important but zero values can be disregarded.
#This saves time interating through trailing zero coefficients, and also allows equivalent cavity objects with different coefficient matrix sizes to be hashed identically.
#Since this is run only at the creation of each instance of a Boundary or RefractiveIndex object, it does not affect computation time.
function coeffmatrix{T}(m::Array{T,1})
    I::Array{Int64,1},J::Array{Int64,1},V::Array{T,1} = findnz(cat(2,m))
    if isempty(I)
        return Array(T,0)
    else
        return full(sparsevec(I,V))[:]
    end
end
function coeffmatrix{T}(m::Array{T,2})
    I::Array{Int64,1},J::Array{Int64,1},V::Array{T,1} = findnz(m)
    if isempty(I)
        return Array(T,0,0)
    else
        return full(sparse(I,J,V))
    end
end


text:
#=
  Demos for standard glm algos
=#

using glmSolverjl
using Random: seed!

seed!(0);
gammaX, gammaY = simulateData(Float64, GammaDistribution(), LogLink(), 30, 10_000);
gammaBlockX = matrixToBlock(gammaX, 10);
gammaBlockY = matrixToBlock(gammaY, 10);

binomialX, binomialY = simulateData(Float64, BinomialDistribution(), LogLink(), 30, 10_000);
binomialBlockX = matrixToBlock(binomialX, 10);
binomialBlockY = matrixToBlock(binomialY, 10);

poissonX, poissonY = simulateData(Float64, BinomialDistribution(), LogLink(), 30, 10_000);
poissonBlockX = matrixToBlock(poissonX, 10);
poissonBlockY = matrixToBlock(poissonY, 10);

# Fit the models
gamma_distrib_log_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse());
gamma_distrib_inverse_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse());
gamma_distrib_identity_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), IdentityLink(), GESVSolver(),  inverse = GETRIInverse());

gamma_distrib_power_model_1 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), PowerLink(0.0), GESVSolver(), inverse = GETRIInverse());
gamma_distrib_power_model_2 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), PowerLink(1/3), GESVSolver(), inverse = GETRIInverse());
gamma_distrib_negative_binomial_link_1 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), NegativeBinomialLink(1.0), GESVSolver(), inverse = GETRIInverse());
gamma_distrib_negative_binomial_link_2 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), NegativeBinomialLink(2.0), GESVSolver(), inverse = GETRIInverse());

#= Binomial Distribution With Logit Link =#
binomial_logit_link = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse());
#= Binomial Distribution With Probit Link =#
binomial_probit_link = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), ProbitLink(), GESVSolver(), inverse = GETRIInverse());
binomial_cauchit_link = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), CauchitLink(), GESVSolver(), inverse = GETRIInverse());
binomial_oddspower_link_1 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse());

binomial_distrib_odds_power_link_1 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse());
# binomial_distrib_odds_power_link_2 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())

# binomial_oddspower_link_2 = glm(RegularData(), binomialX, binomialY, BinomialDistribution, OddsPowerLink{1})

# bernoulli_logcomplementary = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogComplimentLink(), GESVSolver(), inverse = GETRIInverse())
bernoulli_loglog = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogLogLink(), GESVSolver(), inverse = GETRIInverse());
# bernoulli_complementaryloglog = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), ComplementaryLogLogLink(), GESVSolver(), inverse = GETRIInverse())

# LogLink With Gaussian Distribution
log_link_gaussian_distrib = glm(RegularData(), gammaX, gammaY, GaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())

# LogLink With Gamma Distribution
log_link_gamma_distribution = glm(RegularData(), gammaX, gammaY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
log_link_inversegaussian_distribution = glm(RegularData(), gammaX, gammaY, InverseGaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())

log_link_poisson_distribution = glm(RegularData(), poissonX, poissonY, PoissonDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
logit_link_bernoulli_distrib = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())

log_link_negative_bernoulli_distrib = glm(RegularData(), gammaX, gammaY, NegativeBernoulliDistribution(0.5), LogLink(), GESVSolver(), inverse = GETRIInverse())
cauchit_link_binomial_distribution = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), CauchitLink(), GESVSolver(), inverse = GETRIInverse())

#=======================================================================================#
igaussianLogModel = glm(RegularData(), gammaX, gammaY, InverseGaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
igaussianInverseModel = glm(RegularData(), gammaX, gammaY, InverseGaussianDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse())

bernoulliLogit = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
bernoulliProbit = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), ProbitLink(), GESVSolver(), inverse = GETRIInverse())

bernoulliOddsPower0 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse())

# bernoulliLogComplementary = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogComplimentLink(), GESVSolver(), inverse = GETRIInverse())
bernoulliLogLog = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogLogLink(), GESVSolver(), inverse = GETRIInverse())

#
# bernoulliComplementaryLogLog = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), ComplementaryLogLogLink(), GESVSolver(), inverse = GETRIInverse())

# Count Data
poissonLog = glm(RegularData(), poissonX, poissonY, PoissonDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
# Here the parameter is the inverse of what it is in R
# negativeBernoulliLog = glm(RegularData(), poissonX, poissonY, NegativeBernoulliDistribution(1/2), LogLink(), GESVSolver(), inverse = GETRIInverse())

# Education for full binomial data
binomialLogitModel = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())

# Testing the Cauchit Link function
binomialCauchit = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), CauchitLink(), GESVSolver(), inverse = GETRIInverse())

#=
TODO:

1. Do test/examples for offset and weights
=#

# Block Demos
block_gamma_distrib_log_link = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
block_gamma_distrib_inverse_link = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse())
block_gamma_distrib_identity_link = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), IdentityLink(), GESVSolver(), inverse = GETRIInverse())
block_gamma_distrib_power_model_1 = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), PowerLink(0.0), GESVSolver(), inverse = GETRIInverse())
block_gamma_distrib_negative_binomial_link_2 = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), NegativeBinomialLink(2.0), GESVSolver(), inverse = GETRIInverse())
block_log_link_gaussian_distrib = glm(Block1D(), gammaBlockX, gammaBlockY, GaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
block_igaussianLogModel = glm(Block1D(), gammaBlockX, gammaBlockY, InverseGaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
block_igaussianInverseModel = glm(Block1D(), gammaBlockX, gammaBlockY, InverseGaussianDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse())
# block_binomial_oddspower_link_1 = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())
# block_binomial_distrib_odds_power_link_2 = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())
block_cauchit_link_binomial_distribution = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), CauchitLink(), GESVSolver(), inverse = GETRIInverse())
block_poissonLog = glm(Block1D(), poissonBlockX, poissonBlockY, PoissonDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())

block_logit_link_bernoulli_distrib = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
block_binomial_distrib_odds_power_link_1 = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse())
block_bernoulliLogit = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
block_bernoulliOddsPower0 = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse())

block_gamma_distrib_log_link.coefficients |> println
block_gamma_distrib_inverse_link.coefficients |> println
block_gamma_distrib_identity_link.coefficients |> println
block_gamma_distrib_power_model_1.coefficients |> println
block_gamma_distrib_negative_binomial_link_2.coefficients |> println
block_log_link_gaussian_distrib.coefficients |> println
block_igaussianLogModel.coefficients |> println
block_igaussianInverseModel.coefficients |> println
block_cauchit_link_binomial_distribution.coefficients |> println
block_poissonLog.coefficients |> println

block_logit_link_bernoulli_distrib.coefficients |> println
block_binomial_distrib_odds_power_link_1.coefficients |> println
block_bernoulliLogit.coefficients |> println
block_bernoulliOddsPower0.coefficients |> println

gamma_distrib_log_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
gamma_distrib_inverse_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse())
gamma_distrib_identity_link = glm(RegularData(), gammaX, gammaY, GammaDistribution(), IdentityLink(), GESVSolver(), inverse = GETRIInverse())
gamma_distrib_power_model_1 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), PowerLink(0.0), GESVSolver(), inverse = GETRIInverse())
gamma_distrib_negative_binomial_link_2 = glm(RegularData(), gammaX, gammaY, GammaDistribution(), NegativeBinomialLink(2.0), GESVSolver(), inverse = GETRIInverse())
log_link_gaussian_distrib = glm(RegularData(), gammaX, gammaY, GaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
igaussianLogModel = glm(RegularData(), gammaX, gammaY, InverseGaussianDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
igaussianInverseModel = glm(RegularData(), gammaX, gammaY, InverseGaussianDistribution(), InverseLink(), GESVSolver(), inverse = GETRIInverse())
# binomial_oddspower_link_1 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())
# binomial_distrib_odds_power_link_2 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())
cauchit_link_binomial_distribution = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), CauchitLink(), GESVSolver(), inverse = GETRIInverse())
poissonLog = glm(RegularData(), poissonX, poissonY, PoissonDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())

logit_link_bernoulli_distrib = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
binomial_distrib_odds_power_link_1 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse())
bernoulliLogit = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
bernoulliOddsPower0 = glm(RegularData(), binomialX, binomialY, BinomialDistribution(), OddsPowerLink(0.0), GESVSolver(), inverse = GETRIInverse())

gamma_distrib_log_link.coefficients |> println
gamma_distrib_inverse_link.coefficients |> println
gamma_distrib_identity_link.coefficients |> println
gamma_distrib_power_model_1.coefficients |> println
gamma_distrib_negative_binomial_link_2.coefficients |> println
log_link_gaussian_distrib.coefficients |> println
igaussianLogModel.coefficients |> println
igaussianInverseModel.coefficients |> println
binomial_oddspower_link_1.coefficients |> println
cauchit_link_binomial_distribution.coefficients |> println
poissonLog.coefficients |> println

logit_link_bernoulli_distrib.coefficients |> println
binomial_distrib_odds_power_link_1.coefficients |> println
bernoulliLogit.coefficients |> println
bernoulliOddsPower0.coefficients |> println

# Parallel Block Algorithm
block_gamma_distrib_log_link = glm(Block1D(), gammaBlockX, gammaBlockY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())
block_parallel_gamma_distrib_log_link = glm(Block1DParallel(), gammaBlockX, gammaBlockY, GammaDistribution(), LogLink(), GESVSolver(), inverse = GETRIInverse())

# block_binomial_oddspower_link_1 = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())
# block_parallel_binomial_oddspower_link_1 = glm(Block1DParallel(), binomialBlockX, binomialBlockY, BinomialDistribution(), OddsPowerLink(2.0), GESVSolver(), inverse = GETRIInverse())

block_logit_link_bernoulli_distrib = glm(Block1D(), binomialBlockX, binomialBlockY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())
block_parallel_logit_link_bernoulli_distrib = glm(Block1DParallel(), binomialBlockX, binomialBlockY, BinomialDistribution(), LogitLink(), GESVSolver(), inverse = GETRIInverse())

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
   # This is a test file

    #+MyTag1,MyTag2
              #
  # My documentedItem
   #
   function    myFunction (;x::Float64=sin(1))
     return 0
end

    # final comment (no return)


text:
# source:
# https://github.com/JuliaGeometry/GeometryTypes.jl/issues/71#issuecomment-302175458

"""
A `Cone3` is a 3D cylinder defined by its origin point, its extremity
and a radius. `origin`, `extremity` and `r`, must be specified.
"""
struct Cone3{T} <: GeometryPrimitive{3,T}
  origin::Point{3,T}
  extremity::Point{3,T}
  r::T
end

origin(c::Cone3) = c.origin
extremity(c::Cone3) = c.extremity
radius(c::Cone3) = c.r
height(c::Cone3) = norm(c.extremity-c.origin)
direction(c::Cone3) = (c.extremity .- c.origin)./height(c)
function rotation(c::Cone3)
  u = [direction(c)...];
  v = abs(u[1])>0 || abs(u[2])>0 ? [u[2],-u[1],0] : [0,-u[3],u[2]]
  v ./= norm(v)
  w = [u[2]*v[3]-u[3]*v[2],-u[1]*v[3]+u[3]*v[1],u[1]*v[2]-u[2]*v[1]]
  return hcat(v,w,u)
end

isdecomposable(::Type{T}, ::Type{C}) where {T<:Point, C<:Cone3} = true
isdecomposable(::Type{T}, ::Type{C}) where {T<:Face, C<:Cone3} = true

function decompose1(PT::Type{Point{3,T}},c::Cone3{T},facets=12) where {T}
  isodd(facets) ? facets = 2*div(facets,2) : nothing
  facets<8 ? facets = 8 : nothing
  nbv = Int(facets/2)
  M = rotation(c)
  h = height(c)
  position = 1
  vertices = Array{PT}(undef, 2*nbv)
  for j = 1:nbv
    phi                = T((2*pi*(j-1))/nbv)
    vertices[position] = PT(M*[c.r*cos(phi);c.r*sin(phi);0])+PT(c.origin)
    vertices[position+1] = PT([0;0;h])+PT(c.origin)
    position += 2
  end
  vertices
end

function decompose1(::Type{FT},c::Cone3{T},facets=12) where FT<:Face where T
  isodd(facets) ? facets = 2*div(facets,2) : nothing
  facets<8 ? facets = 8 : nothing
  nbv = Int(facets/2)
  indexes = Array{Face{3,Int}}(undef, facets)
  index = 1
  for j = 1:(nbv-1)
    indexes[index] = (index,index+1,index+2)
    indexes[index+1] = (index+2,index+1,index+3)
    index += 2
  end
  indexes[index] = (index,index+1,1)
  indexes[index+1] = (1,index+1,2)
  indexes
end

function ConeMesh(orig, extr, rad, nfacets)
  obj = Cone3(orig, extr, rad)
  vert = decompose1(Point3f0, obj, nfacets)
  fac = decompose1(GLTriangle, obj, nfacets)
  fovMesh = GLNormalMesh(vertices=vert, faces=fac)
end


text:
function save_x_tp(df, t, desired_pressures_gpa, p_sample_gpa, outfile_name)
    qha_basicio_out.save_x_tp(df, t, desired_pressures_gpa, p_sample_gpa, outfile_name)
end

function save_x_tv(df, t, desired_pressures_gpa, p_sample_gpa, outfile_name)
    qha_basicio_out.save_x_tv(df, t, desired_pressures_gpa, p_sample_gpa, outfile_name)
end

save_to_output(fn_output, text) = qha_basicio_out.save_to_output(fn_output, text)

make_starting_string() = qha_basicio_out.make_starting_string()


text:
#=
Given an array of integers, find the first missing positive integer in linear time and constant space. In other words, find the lowest positive integer that does not exist in the array. The array can contain duplicates and negative numbers as well.

For example, the input [3, 4, -1, 1] should give 2. The input [1, 2, 0] should give 3.
=#

using Test
include("Solutions/problem4_lowestmissing.jl")
@test lowestmissing([3, 4, -1, 1]) == 2
@test lowestmissing([1, 2, 0]) == 3

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
module TrussComplianceDemo3D1

using TopOpt, LinearAlgebra, StatsFuns
# using Makie
# using TopOpt.TrussTopOptProblems.TrussVisualization: visualize

# 3D
ndim = 3
node_points, elements, mats, crosssecs, fixities, load_cases = load_truss_json(
    joinpath(@__DIR__, "tim_$(ndim)d.json")
)
ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
loads = load_cases["0"]
problem = TrussProblem(
    Val{:Linear}, node_points, elements, loads, fixities, mats, crosssecs
)

xmin = 0.0001 # minimum density
x0 = fill(1.0, ncells) # initial design
p = 4.0 # penalty
V = 0.5 # maximum volume fraction

solver = FEASolver(Direct, problem; xmin=xmin)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize compliance
    return comp(x)
end
function constr(x)
    # volume fraction constraint
    return sum(x) / length(x) - V
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr)

options = MMAOptions(; maxiter=1000, tol=Tolerance(; kkt=1e-4, f=1e-4))
TopOpt.setpenalty!(solver, p)
@time r = Nonconvex.optimize(
    m, MMA87(; dualoptimizer=ConjugateGradient()), x0; options=options
)

@show r.minimizer
@show obj(r.minimizer)
@show constr(r.minimizer)
# fig = visualize(
#     problem, solver.u; topology = r.minimizer,
#     default_exagg_scale=0.0, exagg_range=100.0,
# )
# Makie.display(fig)

end


text:
function db_create(context, path)
  db = ccall((:grn_db_create, libgroonga),
             Ptr{Void}, (Ptr{Void}, Ptr{UInt8},), context, path)
  if db == C_NULL
    error("Database open failed: $path")
  end
  db
end

function db_open(context, path)
  db = ccall((:grn_db_open, libgroonga),
             Ptr{Void}, (Ptr{Void}, Ptr{UInt8},), context, path)
  if db == C_NULL
    error("Database open failed: $path")
  end
  db
end


text:
# Compatibility shims for old BP versions, when we make a change

for T in [:Linux, :MacOS, :Windows, :FreeBSD]
    @eval $T(arch, libc) = $T(arch; libc=libc)
    @eval $T(arch, libc, call_abi) = $T(arch; libc=libc, call_abi=call_abi)
end

# Compatibility shim to deal with old build.jl files that don't yet use download_info()
function platform_key(machine::AbstractString = Sys.MACHINE)
    Base.depwarn("platform_key() is deprecated, use platform_key_abi() from now on", :binaryprovider_platform_key)
    platkey = platform_key_abi(machine)
    return typeof(platkey)(arch(platkey), libc(platkey), call_abi(platkey))
end

# TODO: fill in better upper bound here when #27674 is backported
# ref: https://github.com/JuliaLang/julia/pull/27674
if v"0.7.0" <= VERSION < v"0.7.1"
    Base.peek(io::Base.AbstractPipe) = Base.peek(Base.pipe_reader(io))
end


text:
function check_min_norm(A, b, x; λ=0.0)
  (nrow, ncol) = size(A);
  if λ > 0.0
    AI = [A sqrt(λ)*I];
    xI = [x ; (b-A*x)/sqrt(λ)];
  else
    AI = A;
    xI = x;
  end
  xmin = AI' * ((AI * AI') \ b);
  xmin_norm = norm(xmin);
  @printf("‖x - xmin‖ = %7.1e\n", norm(xI - xmin));
  return (xI, xmin, xmin_norm)
end


text:
# This file was generated, do not modify it. # hide
hp = X.Horsepower
Xhp = DataFrame(hp1=hp, hp2=hp.^2, hp3=hp.^3);


text:
# Functions for the minimum-eigenvalue method

"""
    splitting(t1, t2, window_start=starttime(t1), window_end=endtime(t1); nphi=$SPLIT_NPHI, ndt=$SPLIT_NDT, dt_max=$SPLIT_DT_MAX, xcorr=true) -> results

Perform a search over a pair of Seis traces, `t1` and `t2`, for the smallest value of the
minimum eigenvalue of the covariance matrix between the traces, for a set of `nphi`×`ndt`
shear wave splitting operators, up to `dt_max` s.

## Output

`results` is a `SeisSplit.Result` containing:

- `phi`: The set of fast shear wave orientations in °
- `dt`: The set of delays times in s
- `lam1`: The larger eigenvalues at each [phi,dt] point
- `lam2`: The smaller eigenvalues at each point
- `phi_best` and `dphi`: The best ϕ and its 1σ uncertainty.  ϕ is measured
  clockwise from local north (towards east) in °.
- `dt_best` and `ddt`: The best δt and its 1σ uncertainty, in s
- `spol` and `dspol`: The source polarisation and an estimate of its uncertainty for the
  best-fitting ϕ and δt.  `spol` is given in ° clockwise of local north.
- `trace1` and `trace2`, the original input traces, where `trace2` is clockwise of `trace1`
- `window_start`, `window_end`, the analysis time window end points.
- `ndf`, the number of degrees of freedom in the signal.

If `xcorr` is `true`, then the rotation correlation map for the pair
of traces is also computed and the following additional fields are
present in `results`:

- `xcorr_phi`: Angles over which rotation correlation was calculated (°)
- `xcorr_dt`: Delay times over which rotation correlation was calculated (s)
- `xcorr_map`: Cross correlation at each [phi,dt] point
- `xcorr_phi_best`: Fast orientation of maximum cross correlation (°)
- `xcorr_dt_best`: Delay time of maximum cross correlation (s)
"""
function splitting(t1::Seis.Trace{T,V}, t2::Seis.Trace{T,V},
                   window_start=Seis.starttime(t1), window_end=Seis.endtime(t1);
                   xcorr=true,
                   nphi=SPLIT_NPHI, ndt=SPLIT_NDT, dt_max=SPLIT_DT_MAX) where {T,V}
    t1, t2 = check_trace_order(t1, t2) # t2 is now clockwise of t1
    n, e = Seis.rotate_through(t1, t2, -t1.sta.azi)
    Seis.cut!.((n, e), window_start, window_end)
    phi = range(-90., stop=90., length=nphi)
    dt = range(0.0, stop=dt_max, length=ndt)
    lam1, lam2 = zeros(T, nphi, ndt), zeros(T, nphi, ndt)
    N, E = deepcopy(n), deepcopy(e)
    @inbounds for j in eachindex(dt)
        # Only need to compute zero-δt once
        if j == firstindex(dt)
            lam1₀, lam2₀ = compute_eigvals(n, e, zero(first(phi)), dt[j], N, E)
            lam1[:,j] .= lam1₀
            lam2[:,j] .= lam2₀
            continue
        end
        for i in eachindex(phi)
            lam1[i,j], lam2[i,j] = compute_eigvals(n, e, phi[i], dt[j], N, E)
        end
    end
    # Best parameters
    ip, idt = minloc2(lam2)
    phi_best = phi[ip]
    dt_best = dt[idt]
    # Find source polarisation and error therein
    spol, dspol = sourcepol(n, e, phi[ip], dt[idt])
    # Errors in splitting parameters
    dphi, ddt, ν = errors_scale_lam2!(lam2, N, E, n, e, phi_best, dt_best, spol, phi, dt)

    # Rotation correlation if desired
    xc_time, xc_phi, xc_map, xc_phi_best, xc_dt_best = if xcorr
        phi_xcorr = range(-90, 90, length=nphi÷2+1)[1:end-1]
        _rotation_correlation!(n, e, phi_xcorr, dt_max)
    else
        nothing, nothing, nothing, nothing, nothing
    end

    # Return result
    Result(phi, dt, lam1, lam2, phi_best, dphi, dt_best, ddt, spol, dspol, t1, t2,
           window_start, window_end, ν,
           # rotation correlation parameters
           xc_phi, xc_time, xc_map, xc_phi_best, xc_dt_best)
end

"""
    compute_eigvals(t1, t2, phi, dt, T1, T2) -> λ₁, λ₂

Compute the two eigenvalues for traces for a specific fast orientation, `phi`°,
and delay time, `dt` s.  `t1` and `t2` are assumed to be orthogonal, and `T1`
and `T2` are a pair of temporary storage traces matching `t1` and `t2`, respectively.
"""
function compute_eigvals(t1::Seis.Trace, t2::Seis.Trace, phi::Number, dt::Number,
                         T1::Seis.Trace, T2::Seis.Trace)
    # Note that we swap the delay time, because we are applying the inverse
    # splitting operator.
    copy_trace!(T1, t1)
    copy_trace!(T2, t2)
    apply_split!(T1, T2, phi, -dt)
    c = covar(T1, T2)
    eval = eigvals(c)
    maximum(eval), minimum(eval)
end

"""
    sourcepol(t1, t2, phi, dt) -> spol, dspol

Compute the source polarisation `spol` and error `dspol` for the pair of traces `t1`
and `t2` for the best-fitting splitting parameters `phi`° and `dt` s.
"""
function sourcepol(t1, t2, phi::Number, dt::Number)
    T1, T2 = deepcopy(t1), deepcopy(t2)
    apply_split!(T1, T2, phi, -dt)
    c = covar(T1, T2)
    eval, evec = eigen(c)
    i1 = argmax(eval)
    i2 = 3 - i1
    lam1 = eval[i1]
    lam2 = eval[i2]
    spol = mod(rad2deg(atan(evec[2,i1], evec[1,i1])), 180)
    dspol = rad2deg(atan(lam2/lam1))
    spol, dspol
end

"""
    covar(t1, t2) -> c

Return the covariance matrix `c` for two traces.
"""
function covar(t1::Seis.Trace{T,V,S}, t2::Seis.Trace{T,V,S}) where {T,V,S}
    c11 = c22 = c12 = zero(eltype(V))
    @inbounds for i in 1:Seis.nsamples(t1)
        c11 += t1.t[i]^2
        c22 += t2.t[i]^2
        c12 += t1.t[i]*t2.t[i]
    end
    c = SArray{Tuple{2,2},T}(c11, c12, c12, c22)
    c
end

"""
    errors_scale_lam2!(lam2, N, E, n, e, phi_best, dt_best, spol, phi, dt) -> σϕ, σδt, ν

Return the 1σ errors in ϕ and δt, `σϕ` and `σδt` plus the number of degrees of freedom
`ν`, for the best-fitting splitting parameters
`phi_best` and `dt_best`, and scale the λ₂ surface, `lam2`, such that the 95% confidence contour
has a value of 1.  `spol` is the retrieved source polarisation.
`n` and `e` are the north and east input traces, and `N` and `E` are
temporary traces which will be rotated to the fast-slow direction to compute
the errors.  `phi` and `dt` are the values of ϕ and δt along the axes of `lam2`.
"""
function errors_scale_lam2!(lam2, N, E, n, e, phi_best, dt_best, spol, phi, dt)
    copy_trace!(N, n)
    copy_trace!(E, e)
    # Rotate to fast-slow and remove splitting
    rotate_traces!(N, E, -phi_best)
    shift_trace!(N, dt_best)
    # Rotate to polarisation direction.  E should be ~ 0 and thus represent noise
    rotate_traces!(N, E, phi_best - spol)
    ν = degrees_of_freedom(E.t)
    if ν <= 2
        @warn("Degrees of freedom ≤ 2; resetting to 3.  σ may be inaccurate.")
        ν = 3
    end
    # Critical F-statistic at 95% F(k, ν-2) where k = 2
    Fν = Distributions.quantile(Distributions.FDist(2, ν - 2), 0.95)
    # Scale λ₂
    λ₂_min = minimum(lam2)
    if λ₂_min < eps(λ₂_min)
        @warn("Minimum λ₂ is zero (can occur in synthetics).  " *
              "Normalising by non-zero value.  σ may be inaccurate.")
        λ₂_min = √eps(λ₂_min)
    end
    # Normalise 95% contour to 1
    lam2 .= lam2./(1 + 2Fν/(ν - 2))./λ₂_min
    # Get errors by scanning the error surface
    σ_ϕ, σ_δt = errors_from_normalised_lam2(lam2, phi, dt)
    σ_ϕ, σ_δt, ν
end

"""
    degrees_of_freedom(f) -> ndf

Return an estimation of the degrees of freedom in the noise vector `f`.
"""
function degrees_of_freedom(f)
    # Compute amplitude spectrum
    F = abs.(FFTW.rfft(f))
    E₂ = sum(x->x^2, F)
    E₂ -= (F[1]^2 + F[end]^2)/2
    E₄ = 1/3*(F[1]^4 + F[end]^4)
    for i in 2:(length(F)-1)
        E₄ += 4/3*F[i]^4
    end
    ndf = round(Int, 4E₂^2/E₄ - 2)
end

"""
    errors_from_normalised_lam2(lam2, phi, dt) -> σ_ϕ, σ_δt

Compute the 1σ uncertainties in ϕ and δt for the normalised λ₂ surface `lam2`,
given the values of ϕ `phi` and δt `dt` along the axes of `lam2`.
"""
function errors_from_normalised_lam2(lam2, phi, dt)
    I = findall(x->x<=1, lam2)
    idt1 = minimum(x->x[2], I)
    idt2 = maximum(x->x[2], I)
    σ_δt = max((idt2 - idt1)*step(dt)/4, step(dt)/4)
    nphi = length(unique(getindex.(I, 1)))
    σ_ϕ = max(nphi*step(phi)/4, step(phi)/4)
    σ_ϕ, σ_δt
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
using Ripemd
using Nettle
using BenchmarkTools

reload("Ripemd")
d = b"";
@benchmark Ripemd.ripemd160($d)
@benchmark Nettle.digest("ripemd160", $d)

d = b"a";
@benchmark Ripemd.ripemd160($d)
@benchmark Nettle.digest("ripemd160", $d)

d = [b"a"[1] for j in 1:10000];
@benchmark Ripemd.ripemd160($d)
@benchmark Nettle.digest("ripemd160", $d)

d = *(["a" for j in 1:10000]...);
@benchmark Ripemd.ripemd160($d)
@benchmark Nettle.digest("ripemd160", $d)


text:
# non equivalent common item design

"""
    FreqTab(tab, raw, interval stats)
The basic struct for all equating.

- `table` A data frame of test scales and thier frequaencies.
- `raw` A raw vector. The order of examinees should correspond to another `FreqTab.raw`
- `interval` The interval of the score scale.
- `stats` Basic stats of the frequency table.
"""
mutable struct FreqTab <: EG
    table
    raw
    interval
    stats
end

function basicstats(X)
    N = length(X)
    Nm = sum(X .=== missing)
    mins, maxs = minimum(X), maximum(X)
    (N = N, Missing = Nm, min = mins, max = maxs, μ = mean(X), σ = std(X), kurtosis = kurtosis(X), skewness = skewness(X))
end

# Frequaency table for equivalent group design
"""
    freqtab(X; interval = 1.0, scale = minimum(X):interval:maximum(X))

Create `FreqTab`, which is used for SG design methods in Equate package.

# Arguments

- `X` Vector of raw test score that dose not contain missing value. The order of examinees should correspond to another vector to be equated.
- `interval` The interval size of scale (must be Float64). Default is 1.0. If `scale` was specified by user, ignored.
- `scale` Vector or StepRange represents a scale of test score X.
"""
function freqtab(X; interval = 1.0, scale = minimum(X):interval:maximum(X))
    stats = basicstats(X)
    freq = map(j -> count(i -> i == j, X), scale)
    cumfreq = cumsum(freq)
    cumprob = cumsum(freq) ./ sum(freq)
    res = DataFrame(scale = scale, freq = freq, cumfreq = cumfreq, prob = freq ./ sum(freq), cumprob = cumprob)
    return FreqTab(res, X[X .!== missing], interval, stats)
end

# equivalent group design
mutable struct NEATFreqTab <: NEAT
    tableX
    tableV
    rawX # independent form
    rawV # common form
    intervalX
    intervalV
    marginal # conditional freqency
    statsX
    statsV
end
# Frequency table for nonequivalent gtoup design.
"""
    freqtab(X, V;intervalX = 1.0, intervalV = 1.0, scaleX = minimum(X):intervalX:maximum(X), scaleV = minimum(V):intervalV:maximum(V))

Create `NEATFreqTab` for NEAT design.

# Arguments

- `X` Vector of raw score, except missing values, of target test.
- `V` Vector of raw score, except missing values, of anchor test.

"""
function freqtab(X, V; intervalX = 1.0, intervalV = 1.0, scaleX = minimum(X):intervalX:maximum(X), scaleV = minimum(V):intervalV:maximum(V))
    if length(X) != length(V)
        println("X and V must be same length(test scores of which the same group).")
    end
    # missing data hangling
    if any(ismissing.(X)) || any(ismissing.(V))
        println("There are more than one missing value in X or V. Use the listwise deletion.")
        x_key = X .!== missing
        v_key = V .!== missing
        X = X[x_key .& x_key]
        V = V[v_key .& v_key]
    end
    statsX = basicstats(X)
    statsV = basicstats(V)
    freqx = map(j -> count(i -> i == j, X), scaleX)
    cumprobx = cumsum(freqx)./ sum(freqx)
    freqv = map(j -> count(i -> i == j, V), scaleV)
    cumprobv = cumsum(freqv)./ sum(freqv)
    tableX = DataFrame(scale = scaleX, freq = freqx, cumfreq = cumsum(freqx), cumprob = cumprobx, prob = freqx ./ sum(freqx))
    tableV = DataFrame(scale = scaleV, freq = freqv, cumfreq = cumsum(freqv), cumprob = cumprobv, prob = freqv ./ sum(freqv))
    marginaltable = zeros(Int64, length(scaleX), length(scaleV))
    for (xi, xv) in enumerate(scaleX), (vi, vv) in enumerate(scaleV)
        marginaltable[xi,vi] = count(i -> i == vv, V[X .== xv])
    end
    return NEATFreqTab(tableX, tableV, X, V, intervalX, intervalV, marginaltable, statsX, statsV)
end
# Frequency table from smoothed frequency table.
"""
    freqtab(X::EG, V::EG)

Create `NEATFreqTab` which has been smoothed (Log Linear or Kernel method).

# Arguments

- `X` Smoothed `FreqTab` of target test.
- `V` Smoothed `FreqTab` of anchor test.
"""
function freqtab(X::EG, V::EG)
    marginaltable = zeros(Int64, length(X.table.scale), length(V.table.scale))
    for (xi, xv) in enumerate(X.table.scale), (vi, vv) in enumerate(V.table.scale)
        marginaltable[xi,vi] = count(i -> i == vv, V.raw[X.raw .== xv])
    end
    return NEATFreqTab(X.table, V.table, X.raw, V.raw, X.interval, V.interval, marginaltable)
end


text:
include("lib_for_robot.jl")
function task7(r::Robot)
    x=0
    y=0
    x,y=tostart(r)
    while isborder(r,Ost)==0
        movetonord(r,x,y)
        move!(r,Ost)
        x+=1
        movetosud(r,x,y)
    end
    if ((x+y)%2==0)
        putmarker!(r)
    end
end