text:
function plotdata(obs::T, fit=nothing) where T<:R2Data
    plt = plot(framestyle=:box,
        xguide="Magnetic field strength / T",
        yguide="Relaxation rate / s⁻¹",
        legend=:outertopright)

    scatter!(plt, obs.B0, obs.ZQ, label="ZQ")
    isnothing(fit) || plot!(plt, fit.B0, fit.ZQ, primary=false)
    scatter!(plt, obs.B0, obs.DQ, label="DQ")
    isnothing(fit) || plot!(plt, fit.B0, fit.DQ, primary=false)
    scatter!(plt, obs.B0, obs.QQ⁻, label="QQ⁻")
    isnothing(fit) || plot!(plt, fit.B0, fit.QQ⁻, primary=false)
    scatter!(plt, obs.B0, obs.QQ⁺, label="QQ⁺")
    isnothing(fit) || plot!(plt, fit.B0, fit.QQ⁺, primary=false)

    return plt
end

function plotdata(obs::T, fit=nothing) where T<:CPMGData
    plt = plot(framestyle=:box,
        xguide="νCPMG / Hz",
        yguide="Relaxation rate / s⁻¹",
        legend=nothing)

    scatter!(plt, obs.νCPMG, obs.R2eff)
    isnothing(fit) || plot!(plt, fit.νCPMG, fit.R2eff, primary=false)

    return plt
end

function plotdata(obs::T, fits=nothing) where T<:Vector{CPMGData}
    plt = plot(framestyle=:box,
        xguide="νCPMG / Hz",
        yguide="Relaxation rate / s⁻¹",
        legend=:outertopright)

    if isnothing(fits)
        fits = fill(nothing,length(obs))
    end
    for (ob, fit) in zip(obs, fits)
        if isa(ob, CPMGData_MQ)
            type = "MQ"
        elseif isa(ob, CPMGData_SQ1H)
            type = "1H SQ"
        else
            type = "unknown CPMG type"
        end
        label = "$(round(ob.B0,digits=1)) T, $type"
        scatter!(plt, ob.νCPMG, ob.R2eff, label=label)
        isnothing(fit) || plot!(plt, fit.νCPMG, fit.R2eff, primary=false)
    end

    return plt
end

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module MeshIO

using GeometryTypes
using ColorTypes
using Compat
import FileIO

import FileIO: DataFormat, @format_str, Stream, File, filename, stream, skipmagic
@compat import Base.show

include("io/off.jl")
include("io/ply.jl")
include("io/stl.jl")
include("io/obj.jl")
include("io/2dm.jl")

load{format}(fn::File{format}, MeshType=GLNormalMesh) = open(fn) do s
	skipmagic(s)
	load(s, MeshType)
end
save{format}(fn::File{format}, msh::AbstractMesh) = open(fn, "w") do s
	save(s, msh)
end

end # module

ext: jl


text:
using Printf

mutable struct WilliamsBeer
    Iₘᵢₙ::Float64
    Π::Float64
end
Base.zero(::Type{WilliamsBeer}) = WilliamsBeer(0.0, 0.0)

function Base.isequal(p::WilliamsBeer, q::WilliamsBeer)
    isequal(p.Iₘᵢₙ, q.Iₘᵢₙ) && isequal(p.Π, q.Π)
end

Base.:(==)(p::WilliamsBeer, q::WilliamsBeer) = p.Iₘᵢₙ == q.Iₘᵢₙ && p.Π == q.Π

Base.:≈(p::WilliamsBeer, q::WilliamsBeer) = p.Iₘᵢₙ ≈ q.Iₘᵢₙ && p.Π ≈ q.Π

Base.iszero(p::WilliamsBeer) = p.Π == zero(p.Π)

function Base.show(io::IO, ::MIME"text/dot", p::WilliamsBeer)
    @printf io "Π = %0.3f\\nIₘᵢₙ = %0.3f" p.Π p.Iₘᵢₙ
end

function pid!(h::Hasse{<:AbstractVertex{WilliamsBeer}},
             stimulus::AbstractVector{Int},
             responses::AbstractMatrix{Int};
             zero=true)
    if zero
        zero!(h)
    end

    bs = maximum(stimulus)

    L = size(responses, 1)
    ss = subsets(L)
    si = Vector{Vector{Float64}}(undef, length(ss))
    for i in eachindex(ss)
        si[i] = specificinfo(stimulus, responses, ss[i])
    end

    sdist = observe!(Entropy(maximum(stimulus)), stimulus)

    for i in eachindex(h)
        α = h[i]
        for s in 1:bs
            x = si[α[1]][s]
            for k in 2:length(α)
                x = min(x, si[α[k]][s])
            end
            payload(α).Iₘᵢₙ += sdist[s] * x;
        end
        payload(α).Iₘᵢₙ /= sdist.N;
    end

    for i in eachindex(h)
        α = h[i]
        for s in 1:bs
            u = -Inf
            for β in below(α)
                x = si[β[1]][s]
                for k in 2:length(β)
                    x = min(x, si[β[k]][s])
                end
                u = max(u, x)
            end
            if isinf(u)
                u = 0.0
            end
            payload(α).Π += sdist[s] * u
        end
        payload(α).Π = payload(α).Iₘᵢₙ - payload(α).Π / sdist.N
    end

    h
end


text:
# setup for Plots.jl, without plotting does not work
ENV["PLOTS_DEFAULT_BACKEND"]="GR"
ENV["GKSwstype"]="png"

max_stars_repo_name: kipfstuhl/doftiles

max_issues_repo_name: kipfstuhl/doftiles

max_forks_repo_name: kipfstuhl/doftiles


text:
# Juno Activate Enviornment in folder of the project
# using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end
using SemidefiniteProgramming
using Test

@testset "sparsesymmetricmatrix.jl" begin
    SSM = SparseSymmetricMatrix(Float64)
    @test typeof(SSM) == SparseSymmetricMatrix{Float64}
    for i in 1:10
        for j in 1:10
            setindex!(SSM, rand(1)[1], i, j)
        end
    end
    # show(SSM.entries)
    @test size(SSM) == 10
    SSM10 = *(10.0,SSM)
    SSM_copy = copy(SSM)

    @test ==(SSM_copy, SSM) == true
    @test ==(SSM10, *(1.0,SSM)) == false
    @test isapprox(SSM10, *(10.0,SSM)) == true

    delete!(SSM,1,1)
    @test SSM[1,1] == 0
end

@testset "sparsesymmetricmatrix.jl" begin
    SSBM = SparseSymmetricBlockMatrix(Float64)
    @test typeof(SSBM) == SparseSymmetricBlockMatrix{Float64}
    @test  typeof(blocks(SSBM)) == Dict{Any,SparseSymmetricMatrix{Float64}}
    for b in 1:2
        for i in 1:3
            for j in i:3
                setindex!(SSBM, rand(1)[1], b,i, j)
            end
        end
    end
    @test getindex(SSBM,1,1,1) == SSBM.blocks[1][1,1]
    SSBM_copy = copy(SSBM) # What is the point of this copy.
    SSBM_alt = SSBM
    SSBM10 = *(10.0, SSBM)
    # show(SSBM10)
    @test ==(SSBM_copy, SSBM_alt) == true
    @test isapprox(SSBM_copy, *(0.0000001,SSBM_alt)) == false
end

@testset "sparsesdp.jl" begin
    SSDP = SparseSDP(Float64)
    @test ismaximizationproblem(SSDP)
    @test setmaximizationproblem!(SSDP, true)
    @test !normalized(SSDP)

    SSBM = SparseSymmetricBlockMatrix(Float64)
    for b in 1:2
        for i in 1:3
            for j in i:3
                setindex!(SSBM, rand(1)[1], b,i, j)
            end
        end
    end
    setobj!(SSDP,SSBM)
    obj(SSDP)
    obj(SSDP,1)
    obj(SSDP,2)
    freeobj(SSDP)
    setcon!(SSDP,1, 1, 1, 1, 2.3)
    cons(SSDP)

    getcon(SSDP,1)
    setrhs!(SSDP, 1, rand(1)[1])
    rhs(SSDP,1)
    rhsdense(SSDP)
    ncons(SSDP)
end
    ## Errors
    # setfreecon!(SSDP, 1, 1, 1.1)
    # freecons(SSDP)
    # blocksizes(SSDP)
    # SSDP_copy = copy(SSDP)

@testset "normalization.jl" begin
    sanitycheck()
    # IndexMap()
    # newindex!,
    # normalize
end

#
#
#
# @testset "SemidefiniteProgramming.jl" begin
#     # Write your tests here.
#     sdp = SparseSDP(maximize=true)
#
#     setobj!(sdp, 1, 1, 1, 1.0)
#     setobj!(sdp, 2, 1, 1, 2.0)
#     setobj!(sdp, 3, 1, 1, 3.0)
#     setobj!(sdp, 3, 2, 2, 4.0)
#
#     setrhs!(sdp, 1, 10.0)
#     setcon!(sdp, 1, 1, 1, 1, 1.0)
#     setcon!(sdp, 1, 2, 1, 1, 1.0)
#
#     setrhs!(sdp, 2, 20.0)
#     setcon!(sdp, 2, 2, 1, 1, 1.0)
#     setcon!(sdp, 2, 3, 1, 1, 5.0)
#     setcon!(sdp, 2, 3, 1, 2, 2.0)
#     setcon!(sdp, 2, 3, 2, 1, 2.0)
#     setcon!(sdp, 2, 3, 2, 2, 6.0)
#
#     println(obj(solve(sdp, CSDP())))
# end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
# TODO: rename it to something since "Interval Tree" means a specific data structure.
# NestedIntervals? CallTree?
struct IntervalTree
    eventid::EventID
    branches::Union{Vector{IntervalTree},Nothing}
    height::Int  # TODO: maybe `level` is a better name (as its 0-origin)
end

IntervalTree(eventid::EventID) = IntervalTree(eventid, nothing, 0)

function IntervalTree(eventid::EventID, branches::Vector{IntervalTree})
    if isempty(branches)
        return IntervalTree(eventid)
    else
        return IntervalTree(eventid, branches, maximum(i -> i.height, branches) + 1)
    end
end

Base.:(==)(x::IntervalTree, y::IntervalTree) =
    isequal(x.eventid, y.eventid) &&
    isequal(x.height, y.height) &&
    isequal(x.branches, y.branches)

struct IntervalStacks
    trees::Vector{IntervalTree}
    points::Vector{EventID}
    events::Dict{EventID,Event}
    locations::Dict{UInt64,Location}
end

function interval_tree(etbl::AbstractVector{Event}, events::Dict{EventID,Event})
    intervals_data = StructVector(
        (
            approx_taskid = e.approx_taskid,
            eventid = e.eventid,
            start = e.start,
            stop = e.stop,
        ) for e in etbl if e.stop isa TimeNS
    )
    intervals_df = DataFrame(intervals_data; copycols = false)
    trees = IntervalTree[]
    for task in groupby(intervals_df, :approx_taskid)
        intervals = Tables.rowtable(task[:, [:eventid, :start, :stop]])
        sort!(intervals; by = i -> i.start)
        append!(trees, _interval_trees(intervals))
    end
    sort!(trees, by = t -> events[t.eventid].start)
    return trees
end

# Optimization ideas: It's straightforward to parallelize this.  Also, maybe
# galloping search is better than binary search.
function _interval_trees(intervals)
    length(intervals) == 1 && return [IntervalTree(intervals[1].eventid)]
    k1 = firstindex(intervals)
    trees = IntervalTree[]
    while true
        root = intervals[k1]
        offset = searchsortedlast(
            @view(intervals[k1:end]),
            (start = root.stop,),
            by = i -> i.start,
        )
        @assert offset > 0   # no time travel `start <= stop`
        k2 = k1 + offset - 1
        if offset == 1
            t = IntervalTree(root.eventid)
        else
            t = IntervalTree(root.eventid, _interval_trees(view(intervals, k1+1:k2)))
        end
        push!(trees, t)
        k2 == lastindex(intervals) && break
        k1 = k2 + 1
    end
    return trees
end

# [...](@ref EventTrackerBase.@recordinterval) didn't work somehow
"""
    EventTracker.stacks()

Analyze nestings of the recorded intervals to recover the "call trees" of the
functions instrumented with [`@recordinterval`](@ref
EventTracker.@recordinterval) macro.

Currently, the returned object only supports `Plots.plot` API.
"""
EventTracker.stacks
EventTracker.stacks() =
    EventTracker.stacks(EventTracker.event_table(), EventTracker.location_table())

function EventTracker.stacks(etbl::AbstractVector{Event}, ltbl::AbstractVector{Location})
    events = Dict{EventID,Event}(e.eventid => e for e in etbl)
    return IntervalStacks(
        interval_tree(etbl, events),
        EventID[e.eventid for e in etbl if e.eventkind === EVENT_POINT],
        events,
        Dict{UInt64,Location}(l.locationid => l for l in ltbl),
    )
end

function append_owntime!(table, stks::IntervalStacks)
    for tree in stks.trees
        append_owntime!(table, tree, stks.events)
    end
    return table
end

function duration_of(tree::IntervalTree, events::AbstractDict{EventID,Event})
    e = events[tree.eventid]
    return e.stop - e.start
end

function owntime_of(tree::IntervalTree, events::AbstractDict{EventID,Event})
    cumtime = duration_of(tree, events)
    tree.branches === nothing && return cumtime
    subtime = sum(tree.branches) do sub
        duration_of(sub, events)
    end
    return cumtime - subtime
end

function append_owntime!(table, tree::IntervalTree, events::AbstractDict{EventID,Event})
    push!(table, (eventid = tree.eventid, owntime = owntime_of(tree, events)))
    if tree.branches !== nothing
        for sub in tree.branches
            append_owntime!(table, sub, events)
        end
    end
    return table
end

const OwntimeAccumulatorDict = Dict{EventID,typeof((owntime = TimeNS(0), count = 0))}

struct OwntimeAccumulator
    dict::OwntimeAccumulatorDict
end
OwntimeAccumulator() = OwntimeAccumulator(OwntimeAccumulatorDict())

function Base.push!(acc::OwntimeAccumulator, row)
    dict, _ = modify!!(acc.dict, row.eventid) do entry
        if entry === nothing
            Some((owntime = row.owntime, count = 1))
        else
            v = something(entry)
            Some((owntime = row.owntime + v.owntime, count = v.count + 1))
        end
    end
    @assert dict === acc.dict
    return acc
end

function owntime_dict(stks::IntervalStacks)
    acc = append_owntime!(OwntimeAccumulator(), stks)
    return acc.dict
end


text:
#RSDeltaSigma: time-domain signal generators
#-------------------------------------------------------------------------------

#==Input signal generators
===============================================================================#
"""`genTestTone_sin()`

Limiation:
If N is small relative to ftest (you don't generate a whole period), iftest gets rounded to nothing.
Typical ampdB: -3dB
"""
function genTestTone_sin(ampdB::Float64, ftest::Float64; M::Int=1, phi0::Float64=0.0, N::Int=2^12)
	amp = undbv(ampdB) #Test tone amplitude, relative to full-scale.
	iftest = round(Int, ftest*N)
	u = amp*M*sin.((2π/N)*iftest * (0:N-1) .+ phi0)
	return (u, iftest)
end

function genTestTone_quad(ampdB::Float64, ftest::Float64; M::Int=1, phi0::Float64=0.0, N::Int=2^12)
	amp = undbv(ampdB) #Test tone amplitude, relative to full-scale.
	iftest = round(Int, ftest*N)
	u = amp*M*exp.((2π*j/N)*iftest * (0:N-1) .+ phi0)
	return (u, iftest)
end

genTestTone(dsm::RealDSM, ampdB::Float64, ftest::Float64; phi0::Float64=0.0, N::Int=2^12) =
	genTestTone_sin(ampdB, ftest, M=dsm.M, phi0=phi0, N=N)
genTestTone(dsm::QuadratureDSM, ampdB::Float64, ftest::Float64; phi0::Float64=0.0, N::Int=2^12) =
	genTestTone_quad(ampdB, ftest, M=dsm.M, phi0=phi0, N=N)

#==impL1
===============================================================================#
"""`y=impL1(ntf, n=10)`

Compute the impulse response from the comparator
output to the comparator input for the given NTF.
n is the (optional) number of points (10).

This function is useful when verifying the realization
of a NTF with a specified topology.
"""
function impL1(ntf::ZPKData, n::Int=10)
	(z, p) = _zpkdata(ntf)

	lf_den = padr(poly(z), length(p)+1)
	lf_num = lf_den .- poly(p)
	if any(imag.(vcat(lf_num, lf_den)) .!= 0)
		#Complex loop filter
		lfr_den = real.( conv(lf_den, conj.(lf_den)) )
		lfr_num = conv(lf_num, conj.(lf_den))
		lf_i = _tf(real.(lfr_num), lfr_den, 1)
		lf_q = _tf(imag.(lfr_num), lfr_den, 1)
		(yi,) = _impulse(lf_i, n)
		(yq,) = _impulse(lf_q, n)
		y = yi .+ j*yq
	else
		lf_num = real.(lf_num); lf_den = real(lf_den)
		(y,) = _impulse(_tf(lf_num, lf_den, 1), n)
	end
	return y
end

#==pulse
===============================================================================#
"""`y = pulse(S, tp=[0 1], dt=1, tfinal=10, nosum=true)`

Calculate the sampled pulse response of a ct system. `tp` may be an array of
pulse timings, one for each input.

# Outputs
 - `y` The pulse response.

# Inputs
 - `S`: An LTI object specifying the system.
 - `tp`: An nx2 array of pulse timings.
 - `dt`: The time increment.
 - `tfinal`: The time of the last desired sample.
 - `nosum`: A flag indicating that the responses are not to be summed.
"""
function pulse(S, tp=[0 1], dt::Float64=1.0, tfinal::Float64=10.0, nosum::Bool=false)
	#Check the arguments
	if isdiscrete(S)
		error("S must be a cts-time system.")
	end

	#Compute the time increment
	dd = 1
	for i in 1:prod(size(tp))
		(x, di) = rat(tp[i],1e-3)
		dd = lcm(di,dd);
	end
	(x, ddt) = rat(dt,1e-3)
	(x, df) = rat(tfinal,1e-3)
	delta_t = 1 / lcm( dd, lcm(ddt,df) )
	delta_t = max(1e-3, delta_t) #Put a lower limit on delta_t
	(y1, ty1, xy1) = step(S, 0:delta_t:tfinal)

#= #DEBUG
	p=cons(:plot, title="")
	for i in 1:size(y1, 3)
	push!(p, cons(:wfrm, waveform(collect(ty1[:]), y1[:,1,i])))
	end
	displaygui(p)
=#

	nd = round(Int, dt/delta_t)
	nf = round(Int, tfinal/delta_t)
	ndac = size(tp,1)
	ni = size(S.B,2)
	if mod(ni,ndac) != 0
		error("The number of inputs must be divisible by the number of dac timings.")
		#This requirement comes from the complex case, where the number of inputs
		#is 2 times the number of dac timings. I think this could be tidied up.
	end
	nis = div(ni,ndac) #Number of inputs grouped together with a common DAC timing
	# (2 for the complex case)
	if !nosum #Sum the responses due to each input set
		y = zeros(array_round(tfinal/dt+1),size(S.C,1),nis)
	else
		y = zeros(array_round(tfinal/dt+1),size(S.C,1),ni)
	end

	for i = 1:ndac
		n1 = round(Int, tp[i,1]/delta_t)
		n2 = round(Int, tp[i,2]/delta_t)
		z1 = (n1, size(y1,2), nis)
		z2 = (n2, size(y1,2), nis)
		yy = [zeros(z1); y1[1:nf-n1+1,:,(i-1)*nis+1:i*nis]] .-
		     [zeros(z2); y1[1:nf-n2+1,:,(i-1)*nis+1:i*nis]]
		yy = yy[1:nd:end,:,:]
		if !nosum #Sum the responses due to each input set
			y = y .+ yy
		else
			y[:,:,i] = yy
		end
	end

	return y
end

#==pulse
===============================================================================#
"""`sv = ds_therm(v,M)`

`sv` is an `M` by `length(v)` matrix of `+/-1`s
where `sum(sv)=v` and the `-1`s are located in the bottom rows
"""
function ds_therm(v, M::Int)
	sv = -ones(M, length(v))
	for i in 1:length(v)
		iend = array_round((M+v[i])/2)
		sv[1:iend, i] .= 1
	end
	return sv
end

#Last line


text:
using LightGraphs

export
    max_cut,
    tsp_qubo,
    hamiltonian_cycle_qubo

#max cut
function max_cut(g::SimpleGraph, vars::Variables)
    @assert vars.indexcapacities == [nv(g), nv(g)]
    hobo = HOBO(vars)
    for edge = edges(g)
        v = edge.src
        w = edge.dst
        hobo -= (vars[v]-vars[w])^2
    end
    hobo
end

function max_cut(g::SimpleGraph; name::String="Max-Cut")
        max_cut(g, Variables(name, nv(g), "vertices"))
end

# hamiltonian cycle
function hamiltonian_cycle_qubo(n::Int, vars::Variables)
    @assert vars.indexcapacities == [n, n]
    hobo = HOBO(vars)

    for t=1:n
        hobo_tmp = HOBO(vars)
        for j=1:n
            hobo_tmp += vars[t, j]
        end
        hobo += (1. - hobo_tmp)^2
    end

    for j=1:n
        hobo_tmp = HOBO(vars)
        for t=1:n
            hobo_tmp += vars[t, j]
        end
        hobo += (1. - hobo_tmp)^2
    end

    #remove staying
    #for v=1:n, t=1:n
    #    hobo += vars[v, t]*vars[v, (t % n) + 1]
    #end
    hobo
end

function hamiltonian_cycle_qubo(n::Int, name::String="Hamiltonian Cycle Comlete Graph")
    hamiltonian_cycle_qubo(n, Variables(name, [n, n], ["timestep", "vertex"]))
end

function hamiltonian_cycle_qubo(g::SimpleGraph, vars::Variables)
    n = nv(g)
    hobo = hamiltonian_cycle_qubo(n, vars)

    #remove according to graph
    for edge = edges(LightGraphs.complement(g))
        u = edge.src
        v = edge.dst
        hobo += vars[t, u]*vars[t % n + 1, v]
        hobo += vars[t, v]*vars[t % n + 1, u]
    end

    hobo
end

function hamiltonian_cycle_qubo(g::SimpleGraph; name::String="Hamiltonian Cycle")
        hamiltonian_cycle_qubo(g, Variables(name, [nv(g),nv(g)], ["timestep", "vertex"]))
end

# tsp
function tsp_qubo(w::Matrix{T}, vars::Variables; B::Real=1., A::Real=2*B*maximum(w)) where T<:Real
    n = size(w, 1)
    @assert size(w, 1) == size(w, 2)
    @assert vars.indexcapacities == [n, n]
    @assert A >= 0 && B >= 0

    hobo_hamiltonian = hamiltonian_cycle_qubo(n, vars)
    hobo_tsp = HOBO(vars)
    #remove according to graph
    for u = 1:n, v = 1:n, t = 1:n
        if u != v
            hobo_tsp += w[u, v] * vars[t, u] * vars[t % n + 1, v]
        end
    end
    A*hobo_hamiltonian + B*hobo_tsp
end

function tsp_qubo(w::Matrix; name::String="TSP", B::Real=1., A::Real=2*B*maximum(w))
    tsp_qubo(w, Variables(name, collect(size(w)), ["timestep", "city"]), B=B, A=A)
end


text:
import Flux

include("../../src/TensorOps.jl")
include("../../src/Embedding.jl")
include("../../src/PositionalEncoder.jl")
include("../../src/Attention.jl")
include("../../src/SequenceTools.jl")
include("../../src/MultiHead.jl")
include("../../src/FeedForward.jl")

struct SequenceFeed
    T::Integer
    V::Integer
    Emb
end

SequenceFeed(Dims::Integer, VocabSize::Integer, TimeSteps::Integer) = SequenceFeed(TimeSteps, VocabSize, Embedding(VocabSize, Dims))
Flux.@treelike SequenceFeed

(m::SequenceFeed)(x::AbstractArray{T, 1} where T) = SequencePad(m.Emb(x), m.T)

struct Encoding
    S::Integer
    IDrop
    HDrop
    PosEnc
    MHA
    FFN
    Attn
end

Encoding(Dims::Integer, TimeSteps::Integer, MHALayers::Integer, MHAStacks::Integer) = Encoding(
    MHAStacks,
    Dropout(0.1), Dropout(0.1),
    PositionEncoder(TimeSteps, Dims), MultiHeadAttention(Dims, MHALayers),
    PointwiseFeedForward(Dims),
    Attention(Dims))
Flux.@treelike Encoding

function (m::Encoding)(x::AbstractArray{T, 3}, SeqLens::AbstractArray{Int32, 1}) where T
    Mask = SequenceMask(x, SeqLens)

    x = m.PosEnc(x, Mask)
    x = m.IDrop(x)

    # multi-head attention
    for _ in 1:m.S
        x = m.MHA(x, x, x)
    end
    x = m.HDrop(x) .+ x
    x = BatchLayerNorm(x)

    # feed forward network
    x = m.FFN(x)
    x = m.HDrop(x) .+ x
    x = BatchLayerNorm(x)

    context = m.Attn(x .* Mask)
    return x, Mask, context
end

struct Decoding
    S::Integer
    IDrop
    HDrop
    PosEnc
    MHA
    FFN
end

Decoding(Dims::Integer, TimeSteps::Integer, MHALayers::Integer, MHAStacks::Integer) = Decoding(
    MHAStacks,
    Dropout(0.1), Dropout(0.1),
    PositionEncoder(TimeSteps, Dims), MultiHeadAttention(Dims, MHALayers),
    PointwiseFeedForward(Dims))
Flux.@treelike Decoding

function (m::Decoding)(x::AbstractArray{T, 3}, y::AbstractArray{T, 3}, SeqLens::AbstractArray{Int32, 1}, Context::AbstractArray{T, 2}, EncMask::AbstractArray{T, 3}, Attn::Attention) where T
    Mask = SequenceMask(x, SeqLens)

    y = m.PosEnc(y, Mask)
    y = m.IDrop(y)

    # multi-head attention
    for _ in 1:m.S
        y = m.MHA(y, y, y, futuremask=true)
    end
    y = m.HDrop(y) .+ y
    y = BatchLayerNorm(y)

    CrossContext = Attn(x .* EncMask, y .* Mask)
    y = Projection(y, CrossContext)
    y = m.HDrop(y) .+ y
    y = BatchLayerNorm(y)
    y = m.FFN(y)
    y = m.HDrop(y) .+ y
    y = BatchLayerNorm(y)

    y = Projection(y, Context)
    return m.HDrop(y) .+ y
end

struct DenseProjection{T}
    Bias::AbstractArray{T, 2}
end

DenseProjection(VocabSize::Integer) = DenseProjection(Flux.param(Flux.gpu(zeros(Float32, 1, VocabSize))))
Flux.@treelike DenseProjection

function (m::DenseProjection)(x::AbstractArray{T, 3}, W::AbstractArray{T, 2}) where T
    x_out = TensorDot(x, permutedims(W, [2, 1]))
    return x_out .+ m.Bias
end


text:
module GPMLj

using PyCall

export
    gpflow,
    instantiate!,
    minimize!,
    predict_f,
    predict_f_samples,
    GPFlowObject

abstract type GPFlowObject end

function instantiate!(o) end
function instantiate!(o::Nothing) return nothing end
function minimize!(opt, m) end
function predict_f(m, Xnew) end
function predict_f_samples(m, Xnew, num_samples) end

include("gpflow.jl")

end # module


text:
export to_device, to_host

using CUDAapi
using Flux
using Knet
using Adapt
import Adapt: adapt, adapt_storage

to_host(x) = to_device(Val(:cpu), x)
to_device(::Val{:cpu}, x) = x  # cpu(x) is not very efficient! So by default we do nothing here.

if has_cuda()
    using CuArrays
    to_device(::Val{:Zygote_gpu}, x) = Flux.fmap(a -> adapt(CuArray, a), x)

    to_device(::Val{:Zygote_gpu}, x::SubArray) = CuArray(x)  # !!! Do not use `adapt` here! For `SubArray`, doing this will send the parent to gpu and this is not what we want for most cases regarding our implementation of buffers.
    to_device(::Val{:Zygote_gpu}, x::SubArray{T, N, P}) where {T, N, P<:CuArray} = x

    # TODO: use adapt
    to_device(::Val{:Zygote_gpu}, x::Base.ReshapedArray) = CuArray(x)
    to_device(::Val{:Zygote_gpu}, x::Base.ReshapedArray{T, N, P}) where {T, N, P<:CuArray} = x

    to_device(::Val{:cpu}, x::Union{KnetArray, CuArray}) = Array(x)
end

to_device(::Val{:Knet_gpu}, x) = Flux.fmap(a -> adapt(KnetArray{Float32}, a), x)
to_device(::Val{:Knet_gpu}, x::SubArray) = KnetArray{Float32}(x)

adapt_storage(T::Type{<:KnetArray}, x::AbstractArray{<:Real}) = T(x)
adapt_storage(T::Type{<:KnetArray}, x::Knet.Param{<:AbstractArray}) = Knet.Param(T(value(x)))


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: personalize_events
using AWS.Compat
using AWS.UUIDs

"""
    put_events(event_list, session_id, tracking_id)
    put_events(event_list, session_id, tracking_id, params::Dict{String,<:Any})

Records user interaction event data. For more information see Recording Events.

# Arguments
- `event_list`: A list of event data from the session.
- `session_id`: The session ID associated with the user's visit. Your application generates
  the sessionId when a user first visits your website or uses your application. Amazon
  Personalize uses the sessionId to associate events with the user before they log in. For
  more information, see Recording Events.
- `tracking_id`: The tracking ID for the event. The ID is generated by a call to the
  CreateEventTracker API.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"userId"`: The user associated with the event.
"""
function put_events(
    eventList, sessionId, trackingId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return personalize_events(
        "POST",
        "/events",
        Dict{String,Any}(
            "eventList" => eventList, "sessionId" => sessionId, "trackingId" => trackingId
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_events(
    eventList,
    sessionId,
    trackingId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return personalize_events(
        "POST",
        "/events",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "eventList" => eventList,
                    "sessionId" => sessionId,
                    "trackingId" => trackingId,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_items(dataset_arn, items)
    put_items(dataset_arn, items, params::Dict{String,<:Any})

Adds one or more items to an Items dataset. For more information see Importing Items
Incrementally.

# Arguments
- `dataset_arn`: The Amazon Resource Name (ARN) of the Items dataset you are adding the
  item or items to.
- `items`: A list of item data.

"""
function put_items(datasetArn, items; aws_config::AbstractAWSConfig=global_aws_config())
    return personalize_events(
        "POST",
        "/items",
        Dict{String,Any}("datasetArn" => datasetArn, "items" => items);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_items(
    datasetArn,
    items,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return personalize_events(
        "POST",
        "/items",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("datasetArn" => datasetArn, "items" => items),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_users(dataset_arn, users)
    put_users(dataset_arn, users, params::Dict{String,<:Any})

Adds one or more users to a Users dataset. For more information see Importing Users
Incrementally.

# Arguments
- `dataset_arn`: The Amazon Resource Name (ARN) of the Users dataset you are adding the
  user or users to.
- `users`: A list of user data.

"""
function put_users(datasetArn, users; aws_config::AbstractAWSConfig=global_aws_config())
    return personalize_events(
        "POST",
        "/users",
        Dict{String,Any}("datasetArn" => datasetArn, "users" => users);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_users(
    datasetArn,
    users,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return personalize_events(
        "POST",
        "/users",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("datasetArn" => datasetArn, "users" => users),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
"""
    pinthread(cpuid::Integer; warn::Bool = true)

Pin the calling Julia thread to the CPU with id `cpuid`.
"""
function pinthread(cpuid::Integer; warn::Bool=true)
    if warn
        (minimum(cpuids_all()) ≤ cpuid ≤ maximum(cpuids_all())) || throw(
            ArgumentError(
                "cpuid is out of bounds ($(minimum(cpuids_all())) ≤ cpuid ≤ $(maximum(cpuids_all()))).",
            ),
        )
        _check_environment()
    end
    return uv_thread_setaffinity(cpuid)
end

"""
    pinthread(threadid::Integer, cpuid::Integer; kwargs...)

Pin the given Julia thread (`threadid`) to the CPU with ID `cpuid`.
"""
function pinthread(threadid::Integer, cpuid::Integer; kwargs...)
    fetch(@tspawnat threadid pinthread(cpuid; kwargs...))
    return nothing
end

"""
    pinthreads(cpuids::AbstractVector{<:Integer}[; warn])
Pins the first `1:length(cpuids)` Julia threads to the CPUs with ids `cpuids`.
Note that `length(cpuids)` may not be larger than `Threads.nthreads()`.

For more information see `pinthread`.
"""
function pinthreads(cpuids::AbstractVector{<:Integer}; warn::Bool=true)
    warn && _check_environment()
    ncpuids = length(cpuids)
    ncpuids ≤ nthreads() ||
        throw(ArgumentError("length(cpuids) must be ≤ Threads.nthreads()"))
    (minimum(cpuids) ≥ minimum(cpuids_all()) && maximum(cpuids) ≤ maximum(cpuids_all())) ||
        throw(
            ArgumentError(
                "All cpuids must be ≤ $(maximum(cpuids_all)) and ≥ $(minimum(cpuids_all())).",
            ),
        )
    @threads :static for tid in 1:ncpuids
        pinthread(cpuids[tid]; warn=false)
    end
    return nothing
end

"""
    pinthreads(strategy::Symbol[; nthreads, warn, kwargs...])
Pin the first `1:nthreads` Julia threads according to the given pinning `strategy`.
Per default, `nthreads == Threads.nthreads()`

Allowed strategies:
* `:compact`: pins to the first `nthreads` cpu threads while trying to avoid using hyperthreads (i.e. moving to next socket before using hyperthreads). If `hyperthreads=true`, hyperthreads will be used before moving to the next socket, if necessary.
* `:scatter` or `:spread` or `sockets`: pins to all available sockets in an alternating / round robin fashion.
* `:numa`: pins to all available NUMA nodes in an alternating / round robin fashion.
* `:random` or `:rand`: pins threads to random cpu threads (ensures that no cpu thread is double occupied). By default (`hyperthreads=false`), hyperthreads will be ignored.
* `:firstn`: pins to the cpuids `0:nthreads-1`
"""
function pinthreads(
    strategy::Symbol; nthreads=Threads.nthreads(), warn::Bool=true, kwargs...
)
    maybe_gather_sysinfo()
    warn && _check_environment()
    if strategy == :compact
        return _pin_compact(nthreads; kwargs...)
    elseif strategy in (:scatter, :spread, :sockets)
        return _pin_scatter(nthreads; kwargs...)
    elseif strategy == :numa
        return _pin_numa(nthreads; kwargs...)
    elseif strategy in (:rand, :random)
        return _pin_random(nthreads; kwargs...)
    elseif strategy == :firstn
        return _pin_firstn(nthreads)
    else
        throw(ArgumentError("Unknown pinning strategy."))
    end
end

_pin_firstn(nthreads) = pinthreads(0:(nthreads - 1); warn=false)
function _pin_random(nthreads; hyperthreads=false)
    if !hyperthreads
        cpuids = shuffle!(filter(!ishyperthread, cpuids_all()))
    else
        cpuids = shuffle(cpuids_all())
    end
    return pinthreads(@view(cpuids[1:nthreads]); warn=false)
end
function _pin_compact(nthreads; hyperthreads=false)
    if !hyperthreads
        cpuids_noht = filter(!ishyperthread, cpuids_all())
        cpuids_ht = filter(ishyperthread, cpuids_all())
        cpuids = @views vcat(cpuids_noht, cpuids_ht)[1:nthreads]
    else
        cpuids = @views cpuids_all()[1:nthreads]
    end
    return pinthreads(cpuids; warn=false)
end
function _pin_scatter(nthreads)
    cpuids = interweave(cpuids_per_socket()...)
    pinthreads(@view cpuids[1:nthreads]; warn=false)
    return nothing
end
function _pin_numa(nthreads)
    cpuids = interweave(cpuids_per_numa()...)
    pinthreads(@view cpuids[1:nthreads]; warn=false)
    return nothing
end

# Potentially throw warnings if the environment is such that thread pinning might not work.
function _check_environment()
    if Threads.nthreads() > 1 && mkl_is_loaded() && mkl_get_dynamic() == 1
        @warn(
            "Found MKL_DYNAMIC == true. Be aware that calling an MKL function can spoil the pinning of Julia threads! Use `ThreadPinning.mkl_set_dynamic(0)` to be safe. See https://discourse.julialang.org/t/julia-thread-affinity-not-persistent-when-calling-mkl-function/74560/3."
        )
    end
    return nothing
end


text:
# Model Zoo: Convolution Neural Network
# =====================

# This notebook is a quick overview of how to implement a
# [convolutional neural network](https://en.wikipedia.org/wiki/Convolutional_neural_network)
# using [Flux](https://github.com/FluxML/Flux.jl) to classify handwritten characters from the
# [MNIST](http://yann.lecun.com/exdb/mnist/) dataset

# MNIST is a collection of 70,000 images of handwritten numbers (0-9) at a resolution of 28x28.
# Those images are split into 60,000 training images and 10,000 testing images.

# ![A sample from the MNIST dataset](mnistsample.png "Some MNIST Entries")A sample from the MNIST dataset

# In this demo we will be using the training part of the dataset to train a convolutional neural
# network to identify which hand-written number is being shown. The network will then be evaluated
# against the testing images to give a score for how well it is able to classify the images.

# This demo has four steps:
# 0. Dependencies
# 1. Preparing the training data
# 2. Defining our convolutional neural network
# 3. Training and evaluating the model

# 0) Dependencies
# =====================

# This demo uses the standard `Flux` package, but also pulls in `Flux.Data.MNIST` to be able to use
# the MNIST dataset in a convenient form. The dataset could be pulled in from disk manually, but this is a nice and simple
# way to be able to get access to the data quickly.

# We also pull in `onehotbatch` and `onecold` for encoding the label data of the MNIST dataset. Lastly
# from Flux we pull in `crossentropy` and `throttle`. The `crossentropy` function is used in determining the loss at each iteration
# and `throttle` is used to control the rate of output to improve the user experience.

# The `Statisitics` package is used here for access to the 'mean' function and the final import of `partition` from Base is to aid
# in the preparation of our test data.

using Statistics
using Flux, Flux.Data.MNIST
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: partition

# Optional extra:
# Using CuArrays will allow for us to place training data and the model onto the GPU and allow for much
# faster training and evaluation. CuArrays currently only supports NVidia GPUs.
using CuArrays

# 1) Preparing the data
# ======================

# To begin, we start by importing the MNIST dataset through the functions made available by `Flux.Data.MNIST`.
# The images are taken as-is, but for the labels we use [`one-hot encoding`](https://github.com/FluxML/Flux.jl/blob/master/docs/src/data/onehot.md).
# One-hot encoding is a common technique in Machine Learning to simplfy complex data into a classification index for each possible result.
# This means changing the label from a single value for each label, such as `2` into a vector where only the second index is set. In our case we have 10 possible
# labels (0,1,2, ... 8,9), so each label is replaced with a vector with 10 entries and the corresponding index for the value is set to `1` and all other values to
# zero. So, an label value of `2` becomes `[0 1 0 0 0 0 0 0 0]`. One-hot encoding allows our classifier to predict how much an image looks like each of the possible
# outputs rather than simply having to predict a single value as an output.

imgs = MNIST.images()
labels = onehotbatch(MNIST.labels(), 0:9)

# Next, we can partition the training data into batches of 1,000 entries.
# This will allow us to do training with a subset of the data.
train = [(cat(float.(imgs[i])..., dims = 4), labels[:,i])
         for i in partition(1:60_000, 1000)]

# When using Flux with CuArrays we want our training data to be placed on the GPU. To do that we have the `gpu(x)` function
# which takes the data that is passed to it and places it on the GPU. This can later be reversed with a call to the `cpu(x)` function.
# Calling this function without CuArrays will do nothing.
train = gpu.(train)

# Now that the training data has been created, the test data needs to be created in the same way.
# In this case we are preparing 1,000 images to evaluate the model. The test images and labels are created
# in the same format as the training data, which means also one-hot encoding the test labels.

# In this code we use the `|>` operator. This operator takes everything before it and wraps it in the function after it.
# So, `10 |> sqrt` becomes sqrt(10). In our code we are using the `gpu` function to place all the test data on the GPU
# in a non-obtrusive way.

# If you have not included the `using CuArrays` call because you are not running on a GPU then the
# calls to `|> gpu` will pass through without attempting to place the code on the GPU, so do not worry about it being there
# if you are running on the CPU.
tX = cat(float.(MNIST.images(:test)[1:1000])..., dims = 4) |> gpu
tY = onehotbatch(MNIST.labels(:test)[1:1000], 0:9) |> gpu

# 2) Defining the convolutional network
# =====================================

# In Flux multiple layers are connected through the `Chain` function
# which allows for the layers to be applied in sequence on an input.
# To create a network with this function we simply pass each layer as
# an input.

# A convolutional layer in Flux is constructed with the `Conv` function.
# The `Conv` function takes a minimum of two parameters.

# Firstly, the size of the convolutional
# mask as a tuple. In this demo we are using `(2,2)` which represents a 2x2 grid of weights. The second
# input is the mapping of input to output channels.

# The standard convolutional network setup involves applying the convolution to the input, then applying
# an `activation` to normalise the output of that convolution (in this demo we are using ReLU) and finally
# the highest value in each convolutional grid is selected and passed through with the `maxpooling` layer.

# These three steps are shown in the first two inputs to the chain function below (the activation is added as the last input to the `Conv` call)

# As with the data above, we append the call to `Chain` with `|> gpu` which will compile these functions to
# the gpu and allow them to be used with the gpu data we have already setup.

m = Chain(
  Conv((2,2), 1=>16, relu),
  x -> maxpool(x, (2,2)),
  Conv((2,2), 16=>8, relu),
  x -> maxpool(x, (2,2)),
  x -> reshape(x, :, size(x, 4)),
  Dense(288, 10), softmax) |> gpu

# In our demo we have two convolutional layers with max pooling and then a dense connection of the results
# to apply the weights to our 10 `one-hot` output.

# The next bit of code applies our training data to the model
# as well as constructing a simple loss function using `crossentropy`
# and an accuracy function we can use to track the progress of the training
# which compares the generated output from the current model for the inputs
# compared to the correct answers.

m(train[1][1])
loss(x, y) = crossentropy(m(x), y)
accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))
opt = ADAM(params(m))

# 3) Training and evaluating the model
# =====================================

# With the model and data in place we can begin training.
# To train one epoch of the model we call the Flux `train!` function.
# this function takes the loss function, training data, our optimisation function
# and then an optional call back function.

# In this demo we use the callback function to evaluate the current model and output
# the performance as it is being trained.
evalcb = throttle(() -> @show(accuracy(tX, tY)), 10)
Flux.train!(loss, train, opt, cb = evalcb)

max_stars_repo_name: NTimmons/model-zoo

max_issues_repo_name: NTimmons/model-zoo

max_forks_repo_name: NTimmons/model-zoo


text:
using LowRank, Test
# using Plots
# using JLD
# using Random, Distributions, LinearAlgebra

N = 30
d = 2
n = 10
k = 1
err = 0.0
indices = 1:n

A = zeros(d, N)
U = rand(d, k)
V = rand(k, n)
A[:, 1:n] = U*V
A[:, n+1:N] = rand(d, N-n)
data_set = collect(eachcol(A))

# file = load(joinpath(@__DIR__, "..","datasets", "N_30_d_2.jld"))
# data_set, (err, indices) = file["data"], file["n_10_k_1"]
# n = 10
# k = 1

# @time @test find_low_rank_subset_checkall(data_set,n,k) == (err, indices)
@test find_low_rank_subset_sample_rep(data_set,n,k)[1] <= 4*err + 0.0001
@test FastLowRank(data_set, n, k)[1] <= err + 0.0001

# PLotting
# p = scatter([Tuple(x) for x in V], label = "all points")
# scatter!(p, [Tuple(V[i]) for i in indices], label = "optimal points")
# plot!(p, [(0., 0.), -1 .* Tuple(truncated_svd(hcat(V[indices]...), 1).U[:,1])], label = "optimal line" )


text:
export write_newick, write_fasta, write_newick!, write_branchlength

"""
"""
function write_newick(file::String, tree::Tree)
	write_newick(file, tree.root)
end

"""
"""
function write_newick(file::String, root::TreeNode)
	out = write_newick!("", root)
	out *= ';'
	f = open(file, "w")
	write(f, out)
	close(f)
end

"""
"""
function write_newick!(s::String, root::TreeNode)
	if !isempty(root.child)
		s *= '('
		temp = sort(root.child, by=x->length(node_leavesclade(x)))
		for c in temp
			s = write_newick!(s, c)
			s *= ','
		end
		s = s[1:end-1] # Removing trailing ','
		s *= ')'
	end
	s *= root.label
	if !ismissing(root.data.tau)
		s *= ':'
		s *= string(root.data.tau)
	end
	return s
end

"""
"""
function write_newick(root::TreeNode)
	return write_newick!("", root)*";"
end

"""
	 write_fasta(file::String, tree::Tree ; internal = false)
"""
function write_fasta(file::String, tree::Tree ; internal = false)
	write_fasta(file, tree.root, internal = internal)
end

"""
	write_fasta(file::String, root::TreeNode ; internal = false)
"""
function write_fasta(file::String, root::TreeNode ; internal = false)
	out = write_fasta!("", root, internal)
	f = open(file, "w")
	write(f, out)
	close(f)
end

"""
"""
function write_fasta!(s::String, root::TreeNode{EvoData}, internal::Bool)
	if internal || root.isleaf
		# s = s * ">$(root.label)\n$(num2seq(root.data.sequence))\n"
		s = s * ">$(root.label)\n$(prod(root.data.sequence))\n"
	end
	for c in root.child
		s = write_fasta!(s, c, internal)
	end
	return s
end

"""
"""
function num2seq(numseq::Array{Int64,1})
	mapping = "ACGT-NWSMKRYBDHV"
	seq = ""
	for a in numseq
		seq *= mapping[a]
	end
	return seq
end

"""
	write_branchlength(tree::Tree, msaname::String, treename::String)

Create a JSON file storing branch lengths of all branches in `tree`. This is intended for reading with auspice.
"""
function write_branchlength(jsonfile::String, tree::Tree, jsontemplate::String)

end

"""
	write_branchlength(tree::Tree, msaname::String, treename::String)

Create a JSON file storing branch lengths of all branches in `tree`. This is intended for reading with auspice.
"""
function write_branchlength(jsonfile::String, tree::Tree, msaname::String, treename::String ; order = keys(tree.lnodes))
	out = Dict()
	node = Dict();
	for k in order
		n = tree.lnodes[k]
		node[n.label] = Dict("branch_length"=>n.data.tau)
	end
	out["nodes"] = node
	out["input_tree"] = treename
	out["alignment"] = msaname
	open(jsonfile, "w") do f
		JSON.print(f, out)
	end
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
using Pkg
Pkg.add(url="https://github.com/carstenbauer/WorkshopWizard.jl/")
using WorkshopWizard

if Sys.iswindows()
    path = joinpath(homedir(), "Desktop")
else
    path = homedir()
end
success = WorkshopWizard.install(;repo="https://github.com/mfherbst/2022-rwth-julia-workshop",
                                  path,
                                  global_IJulia=true,
                                  auto_overwrite=false)

if success
    @info "That's it. Start the notebook server with 'using IJulia; notebook(dir=\"$path/2022-rwth-julia-workshop\")' ..."
end


text:
function serial_number_4d(n)
    @sprintf "_%.4d" n
end

function pathnorepeat(filepath; suffix_fun = serial_number_4d)
    i = 1
    pathv0 = [splitext(filepath)...]
    while isfile(filepath) || isdir(filepath)
        pathv = copy(pathv0)
        insert!(pathv, 2, suffix_fun(i))
        i = i +1
        filepath = join(pathv, "")
    end
    return filepath
end


text:
module Waveguides_FrFT

include("Ankara.jl")
include("JxC.jl")
include("Kravchuk.jl")

end # module


text:
mutable struct outputcom
    T
    H_min
    H_max
    HR
    K
    B
    P
    λ
    type
    δ
    idx
    Θ
    IR
end

function smoothlocalprojection(packagedinput)
    # Unpackage data
    y        = packagedinput.y;
    x        = packagedinput.x;
    w        = packagedinput.w;
    H_min    = packagedinput.H_min;
    H_max    = packagedinput.H_max;
    type     = packagedinput.type;
    r        = packagedinput.r;
    λ        = packagedinput.λ;

    # Assign δ
    if isempty(w)
        δ = std(x);
    else
        δ = std(x - w * inv(w' * w) * w' * x);
    end

    # Is it the regular regression-type?
    isreg = type == "reg";

    T  = length(y);
    HR = H_max + 1 - H_min;

    # construct the B-spline basis functions
    κ = 3;
    if ~isreg
        B = bspline((H_min:H_max)', H_min, (H_max + 1), (H_max + 1 - H_min), κ);
        K = size(B, 2);
    else
        K = HR;
        B = fill(0.0, (H_max, H_max + κ));
        λ = 0;
    end

    # building up the regression representation of the local projection
    idx = nan((H_max+1)*T, 2);
    Y   = nan((H_max+1)*T, 1);
    Xb  = zeros(Float16, (H_max+1)*T, K);
    Xc  = zeros(Float16, (H_max+1)*T, HR, size(w,2)+1 );

    w = hcat(ones(T,1), w);

    for t ∈ 1:(T - H_min)

        idx_beg = (t - 1)*HR + 1;
        idx_end = t * HR;

        idx[idx_beg:idx_end, 1] .= t;
        idx[idx_beg:idx_end, 2] .= H_min:H_max;

        # y
        y_range = (t + H_min) : min((t + H_max), T)';
        Y[idx_beg:idx_end] = [ y[y_range] ; nan(HR-length(y_range),1) ];

        # x
        if isreg
            Xb[ idx_beg:idx_end , : ] = eye(HR) .* x[t];
        else
            Xb[ idx_beg:idx_end , : ] = B * x[t];
        end

        # w
        for i ∈ 1:size(w,2)
            Xc[idx_beg:idx_end , : , i ] = eye(HR)*w[t,i];
        end
    end

    X = Xb;
    for i ∈ 1:size(w,2)
        X = hcat(X, Xc[:, :, i]);
    end

    select = isfinite.(Y)[:];
    idx = idx[select, :];
    Y   = Y[select];
    X   = X[select,:];

    IR = zeros(Float16, H_max+1, 1);

    if isreg
        Θ = (X' * X)\X'*Y
        IR[(H_min+1):end] = Θ[1:K] .* δ;

        # Parameter not used but still defined,
        B = 0;
        P = zeros(Float16, size(X,2), size(X,2) );
        λ = 0;
    else
        P = zeros(Float16, size(X,2), size(X,2));

        D = eye(K);
        for k ∈ 1:r
            D = diff(D, dims = 1);
        end

        P[1:K,1:K] = D' * D;

        Θ = ( (X' * X) + (λ .* P) )\( X'*Y );

        IR[(1+H_min):end] = B * Θ[1:K] .* δ;
    end

    output = outputcom(T, H_min, H_max, HR, K, B,P,λ,type,δ,idx,Θ,IR) ;

    return(output);
    # debug stuff
    # Bs is for display / debugging pourposes only
    # Bs = bspline( (H_min:0.1:H)' , H_min , H+1 , H+1-H_min , 3 );
    # obj.Bs = Bs;
end

function bspline(x, xl, xr, ndx, bdeg)

    dx = (xr - xl) / ndx;

    t = xl .+ dx .* vec(-bdeg:(ndx - 1))'
    T = ones(Float16, length(x), 1) * t;
    X = x' * ones(Float16, 1, length(t));
    P = (X .- T) ./ dx;

    B = (T .<= X) .& (X .< (T .+ dx));
    r = vcat(2:length(t), 1);

    for k = 1:bdeg
        B = (P .* B + (k + 1 .- P) .* B[:, r]) ./ k;
    end
    return B;
end

function matlag(x, lag)
    t, k = size(x);
    res = zeros(Float16, (t-lag), lag * k)

    for i = 1:lag
        res[:, ((i-1)*k + 1):(i*k)] = x[(lag-i + 1):end-i ,:]
    end

    res = nan(t, lag * k);
    for i = 1:lag
        lag_array = x[1:(end-i),:];
        res[(i+1):end, ((i-1)*k + 1):(i*k)] = lag_array;
    end

    return res
end

eye(n) = Matrix{Float16}(I, n, n)

function nan(x, y)
    temp = zeros(Float16, x, y);
    return replace(temp, 0 => NaN)
end


text:
module BenchDataType

using Benchmarkplotting
using Statistics
using Gadfly
using MLStyle
using DataFrames
import Match

import Base.getindex
getindex(asoc_lst :: Vector{Pair{Symbol, T}}, key ::Symbol) where T =
    for (search_key, value) in asoc_lst
        if search_key === key
            return value
        end
    end

@data NormalData begin
    Normal1(a :: Int, b :: Int, c :: Any)
    Normal2(a :: Int)
end

@data GenericData{A} begin
    Generic1(a :: A, b :: Int)
    Generic2(A)
end

"""
ok cases:
1.   Normal1(c :: Number)
2.   Normal1(Normal2(1))
3.   Generic1(a = 3, b = 3)
4.   Generic2((1, 2))
5.   GenericData{String}
"""
data = [
    Symbol("ok1") => Normal1(2, 3, 1.0),

    Symbol("fail1") => Normal1(2, 3, ()),

    Symbol("ok2") => Normal1(2, 3, Normal2(1)),

    Symbol("fail2") => Normal1(2, 3, Normal2(3)),

    Symbol("ok3") => Generic1(3, 3),

    Symbol("fail3") => Generic1(3, -3),

    Symbol("ok4") => Generic2((1, 2)),

    Symbol("fail4") => Generic2(()),

    Symbol("ok5") => Generic2("5"),

    Symbol("fail5") => Generic2(Int)
]

implementations = [
    :MLStyle => (@λ begin
        Normal1(c = ::T where T <: Number) -> 1
        Normal1(c = Normal2(1)) -> 2
        Generic1(a = 3, b = 3) -> 3
        Generic2((1, 2)) -> 4
        ::GenericData{String} -> 5
        _ -> 0
    end),
    Symbol("Match.jl") => function (x)
        Match.@match x begin
            Normal1(_, _, _ :: Number) => 1
            Normal1(_, _, Normal2(1)) => 2
            Generic1(3, 3) => 3
            Generic2((1, 2)) => 4
            _ :: GenericData{String} => 5
            _ => 0
        end
    end
]

criterion(x) = (meantime = mean(x.times), allocs = 1 + x.allocs)
df = bcompare(criterion, data, implementations, repeat=1)

@info df

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size=5px
)

report_meantime, df_time = report(:meantime, df, Scale.y_log10, theme)
report_allocs, df_allocs = report(:allocs, df, theme)

open("stats/vs-match(datatype).txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/vs-match(datatype)-on-time.svg", 10inch, 4inch), report_meantime);
draw(SVG("stats/vs-match(datatype)-on-allocs.svg", 10inch, 4inch), report_allocs);

end


text:
using NamedArrays

function compute_betas{T <: AbstractFloat}(factors::AbstractArray{T,2}, returns::AbstractArray{T,2})
    # make sure the dimension labels match
    if typeof(factors) <: NamedArrays.NamedArray && typeof(returns) <: NamedArrays.NamedArray
        @assert NamedArrays.names(factors,1) == NamedArrays.names(returns,1)
        @assert factors.dimnames[1] == returns.dimnames[1]
    end

    println("DEBUG: factors $factors")
    println("DEBUG: returns $returns")
    ff = factors' * factors
    fr = factors' * returns
    println("DEBUG: ff $ff")
    println("DEBUG: fr $fr")
    # return ff \ fr
    #TODO revert me after fixed!
    result = ff \ fr
    println("DEBUG: result $result")
    return result
end

nTimes = 1000
nSecurities = 100
nFactors = 5

r = rand(nTimes,nSecurities)
returns = NamedArray(r, (map(x->"Time$x",1:nTimes), map(x->"Sec$x",1:nSecurities)), ("Time","Security"))

f = rand(nTimes,nFactors)
factors = NamedArray(f, (names(returns,1), map(x->"Factor$x",1:nFactors)), (returns.dimnames[1],"Factor") )

betas = compute_betas(factors,returns)

# The labels get messed up during the calculations.  The labels should be:
# dim1: Factor -> Factor1, Factor2, ... Factor5
# dim2: Security -> Sec1, Sec2, ... Sec100
println("DEBUG: betas $betas")


text:
@cfgdef struct HGFBertConfig <: HGFConfig
  vocab_size::Int = 30522
  hidden_size::Int = 768
  num_hidden_layers::Int = 12
  num_attention_heads::Int = 12
  intermediate_size::Int = 3072
  hidden_act::String = "gelu"
  hidden_dropout_prob::Float64 = 0.1
  attention_probs_dropout_prob::Float64 = 0.1
  max_position_embeddings::Int = 512
  type_vocab_size::Int = 2
  initializer_range::Float32 = 0.02
  layer_norm_eps::Float32 = 1e-12
  pad_token_id::Int = 0
end

config_type(::Val{:bert}) = HGFBertConfig


text:
## Utilities

to_tuple(v::Vector) = tuple(v...)
to_tuple(t::Union{Tuple, NTuple}) = t

# Apply a function to each dimension of an array, conditional upon a boolean
function mapdim(func::Function, A::Array, flag)
    if ndims(A) != size(flag,1)
        error("Boolean flag must have as many rows as dimensions of A")
    end
    Am = A
    for j = 1:size(flag, 2)
        for i = 1:size(flag, 1)
            if flag[i,j]
                Am = func(Am, i)
            end
        end
    end
    return Am
end

function filledges!(A::Array, val)
    n_dims = ndims(A)
    ind = Array(Any, n_dims)
    for idim = 1:n_dims
        ind[idim] = 1:size(A, idim)
    end
    for idim = 1:n_dims
        ind[idim] = [1, size(A, idim)]
        A[to_tuple(ind)...] = val
        ind[idim] = 1:size(A, idim)
    end
end

# Devectorized linear combinations
lincomb(coef1, A1, coef2, A2) = coef1*A1 + coef2*A2
lincomb(coef1, A1, coef2, A2, coef3, A3) = coef1*A1 + coef2*A2 + coef3*A3
lincomb(coef1, A1, coef2, A2, coef3, A3, coef4, A4) = coef1*A1 + coef2*A2 + coef3*A3 + coef4*A4
function lincomb(coef1::Number, A1::AbstractArray, coef2::Number, A2::AbstractArray)
    out = similar(A1)
    for i = 1:length(A1)
        out[i] = coef1*A1[i] + coef2*A2[i]
    end
    out
end

function lincomb(coef1::Number, A1::AbstractArray, coef2::Number, A2::AbstractArray, coef3::Number, A3::AbstractArray)
    out = similar(A1)
    for i = 1:length(A1)
        out[i] = coef1*A1[i] + coef2*A2[i] + coef3*A3[i]
    end
    out
end

function lincomb(coef1::Number, A1::AbstractArray, coef2::Number, A2::AbstractArray, coef3::Number, A3::AbstractArray, coef4::Number, A4::AbstractArray)
    out = similar(A1)
    for i = 1:length(A1)
        out[i] = coef1*A1[i] + coef2*A2[i] + coef3*A3[i] + coef4*A4[i]
    end
    out
end

# Converting coordinates to linear indices
function coords2lin(c, strides)
    ind = 1
    for i = 1:length(c)
        ind += (c[i]-1)*strides[i]
    end
    ind
end


text:
using FiniteHorizonPOMDPs
using Test
using POMDPLinter
using POMDPModelTools
using DiscreteValueIteration
using POMDPModels: SimpleGridWorld, BabyPOMDP

@testset "interface" begin
    @test HorizonLength(SimpleGridWorld()) == InfiniteHorizon()
    @test HorizonLength(BabyPOMDP()) == InfiniteHorizon()
end

# include("instances/1DGridWorld.jl")
include("instances/1DFiniteHorizonGridWorld.jl")

@testset "1DGridWorld" begin
    include("1dtest.jl")
end


text:
############################################################################
#
#  Software.jl
#
#  Purpose:    Software file to control a quadrotor
#
#  Author:     Edward Daughtery
#
#  Date:       05 August 2016
#
#
############################################################################
include("Guidance.jl")

#######################################################################
# Intialization
function Initialize_SW(time::Float64, DT::Float64, nav_state::NavStates)
    zz = 0.0;
    zi = UInt(0);
    zvec = [0.0;0.0;0.0];
    dt = DT;
    count = UInt(0);
    mode = UInt(0);

    ## Controls Initialization
    throttle = 0.0;
    alpha = 0.0;
    alphadot = 0.0;
    phi = 0.0;
    phidot = 0.0;
    controls = Controls(throttle,alpha,alphadot,phi,phidot);

    ## Guidance Initialization
#    num_waypoints = 3;
#    GuidInit(wyptlist, num_waypoints)

    ## Software Message
    tgo = 999999.0;
    missvec = [9999999.0;9999999.0;9999999.0];
    miss = norm(missvec);
    swMsg = SoftwareMsg( time, dt, mode, controls, nav_state, zi, zz, zz, zz, zz, zz, tgo, miss, missvec, zz );

    return swMsg;
end

#######################################################################
# Flight Control
function FlightControl( softwareMsg::SoftwareMsg, apcmd::Float64, phicmd::Float64 )

    ## Software Telemetry Variable
    swOut = softwareMsg
    time = swOut.time;
    dt =  swOut.dt;

    ## Control Variables
    alpha = swOut.controls.alpha;
    alphadot = swOut.controls.alphadot;
    phi = swOut.controls.phi;
    phidot = swOut.controls.phidot;
    throttle = swOut.controls.throttle;

    ## IMU Measurements
    axmeas = swOut.navstates.accel_meas[1];
    aymeas = swOut.navstates.accel_meas[2];
    azmeas = swOut.navstates.accel_meas[3];

    # Alpha Gains
    Kpa = 0.20;
    Kpad = 0.3;

    # Roll Gain
    Kr = 0.30;

    # Z Acceleration Error
    alphacmd = -Kpa*(apcmd  - azmeas);

    # Alpha Response
    alphadot = Kpad*(alphacmd - alpha);
    tmpalpha = alpha + dt*alphadot;

    # Limits on the Angle of Attack
    if (tmpalpha*57.3 < -6.0 )
        alpha = -6.0/57.3;
    elseif (tmpalpha*57.3 > 9.0 )
        alpha = 9.0/57.3;
    else
        alpha = tmpalpha;
    end

    # Phi Response
    phidot = Kr*(phicmd - phi);
    tmpphi = phi + dt*phidot;

    # Limits on the Roll Angle
    if ( tmpphi*57.3 > 80.0)
        phi = 80.0/57.3;
    elseif ( tmpphi*57.3 < -80.0)
        phi = -80.0/57.3;
    else
        phi = tmpphi;
    end

    swOut.controls.alpha = alpha;
    swOut.controls.alphadot = alphadot;
    swOut.controls.alphadot = phi;
    swOut.controls.alphadot = phidot;

    return swOut
end

#######################################################################
function Steering(nav_state::NavStates, acc_cmd_ned::Array{Float64,1})

    # Computing NED to Body DCM
    Tr_BN = nav_state.euler.Tr_bn;
    roll = nav_state.euler.roll;
    pitch = nav_state.euler.pitch;
    heading = nav_state.euler.heading;

    # Computing Non-Rolled NED To Body DCM
    Tr_BN_NR = TR_BN(0.0, pitch, heading);

    # Acceleration Commands Body
    acc_cmd_body = Tr_BN*acc_cmd_ned;
    acc_cmd_bodyNR = Tr_BN_NR*acc_cmd_ned;

    # Acceleration Commands Non Rolled
    ay_cmd_ned_nr = acc_cmd_bodyNR[2];
    az_cmd_ned_nr = acc_cmd_bodyNR[3];

    # Zero Denominator Protection
    if ( abs(az_cmd_ned_nr) < 0.00000001 )
        phicmd = 0.5*eta_turn_cmd;
    else
        phicmd = atan(ay_cmd_ned_nr/abs(az_cmd_ned_nr + 0.001) );
    end

    # Phi Limiting
    philim = 60.0*d2r;
    phicmd = Bound(phicmd, -philim, philim);

    # Pitch Acceleration Command
    apcmd = acc_cmd_body[3];

    return  [apcmd; phicmd];
end

#######################################################################
function SpeedControl(nav_state::NavStates, vcmd_coop::Float64, softwareMsg::SoftwareMsg )

    ## Parameter Initialization
    vm_ned = nav_state.velocity.vm_ned;
    gamma = nav_state.velocity.gamma;
    alt = nav_state.wgs_pos.altitude;
    Vmag = nav_state.velocity.vmag;
    vel_err_int = softwareMsg.vel_err_int;
    dt50 = softwareMsg.dt;
    Throttle = softwareMsg.controls.throttle;

    ## Atmosphere Model
    atmdata = Atmosphere(alt, Vmag);
    Press = atmdata[1];
    Temp = atmdata[2];
    rho = atmdata[3];
    aa = atmdata[4];
    Mach = atmdata[5];
    qbar = atmdata[6];

    ## Mach Command converted to speed command
    machcmd = 0.7;
    vcmd = machcmd*aa + vcmd_coop;

    ## Velocity Error and Integrator
    verr = (vcmd - Vmag);
    tmp2 = vel_err_int + dt50*verr;

    ## PI Control to compute Throttle Command
    Kpv = 1.5;
    Kiv = 0.66;
    tmp = Kpv*verr + Kiv*tmp2;

    if ( tmp > 100.0 )
        tmp = 100.0;
    elseif ( tmp < 0.0)
        tmp = 0.0;
    else
        Throttle = tmp;
        vel_err_int = tmp2;
    end

    ## Bounding Throttle to be between 0 and 100
    Throttle = Bound(Throttle, 0.0, 100.0);

    return Throttle;
end

#######################################################################
# Software Update
function Update_SW(swMsg::SoftwareMsg)

    swOut = swMsg;
    alpha  = swOut.controls.alpha;
    phi  = swOut.controls.phi;
    gamma  = swOut.navstates.velocity.gamma;
    azimuth  = swOut.navstates.velocity.azimuth;

    ## Update
    swOut.time = swOut.time + swOut.dt;

    ## Acceleration Commands
    acc_cmd_ned = zeros(3);

    ## Guidance Section
    alt_cmd = 4000.0;
    eta_turn_cmd = 0.0;
    eta_dive_cmd = AltitudeControl(swOut.navstates, alt_cmd);
    println(eta_dive_cmd)

    ## Velocity Commands in Steering Frame (Velocity Frame)
    acc_cmd_vel = zeros(3);
    acc_cmd_vel[1] = 0.0;
    acc_cmd_vel[2] = eta_turn_cmd;
    acc_cmd_vel[3] = eta_dive_cmd;

    Tr_vn = TR_VN(gamma, azimuth);
    Tr_nv  = Tr_vn';
    acc_cmd_ned = Tr_nv*acc_cmd_vel;
    acc_cmd_ned[3] = acc_cmd_ned[3]  - 9.80665;

    ## Steering Function
    steerout = Steering(swOut.navstates, acc_cmd_ned);
    apcmd = steerout[1];
    phicmd = steerout[2];

    ## Flight Controls
    throttle = SpeedControl(swOut.navstates, 0.0, swOut);
    FlightControl( swOut, apcmd, phicmd );

    swOut.controls.throttle = throttle;
    return swOut;
end


text:
using Turing, Distributions, Test
using Turing: link, invlink, logpdf_with_trans

# Standard tests for all distributions involving a single-sample.
function single_sample_tests(dist)
  x = rand(dist)
  @test invlink(dist, link(dist, copy(x))) ≈ x atol=1e-9

  y = link(dist, x)
  @test link(dist, invlink(dist, copy(y))) ≈ y atol=1e-9
  logpdf_with_trans(dist, x, true)
  logpdf_with_trans(dist, x, false)

  @test typeof(x) == typeof(y)
end

# Standard tests for all distributions involving multiple samples. xs should be whatever
# the appropriate repeated version of x is for the distribution in question. ie. for
# univariate distributions, just a vector of identical values. For vector-valued
# distributions, a matrix whose columns are identical.
function multi_sample_tests(dist, x, xs, N)
  ys = link(dist, copy(xs))
  @test invlink(dist, link(dist, copy(xs))) ≈ xs atol=1e-9
  @test link(dist, invlink(dist, copy(ys))) ≈ ys atol=1e-9
  @test logpdf_with_trans(dist, xs, true) == fill(logpdf_with_trans(dist, x, true), N)
  @test logpdf_with_trans(dist, xs, false) == fill(logpdf_with_trans(dist, x, false), N)

  @test typeof(xs) == typeof(ys)
end

# Tests with scalar-valued distributions.
uni_dists = [
  Arcsine(2, 4),
  Beta(2,2),
  BetaPrime(),
  Biweight(),
  Cauchy(),
  Chi(3),
  Chisq(2),
  Cosine(),
  Epanechnikov(),
  Erlang(),
  Exponential(),
  FDist(1, 1),
  Frechet(),
  Gamma(),
  InverseGamma(),
  InverseGaussian(),
  Kolmogorov(),
  Laplace(),
  Levy(),
  Logistic(),
  LogNormal(1.0, 2.5),
  Normal(0.1, 2.5),
  Pareto(),
  Rayleigh(1.0),
  TDist(2),
  TruncatedNormal(0, 1, -Inf, 2),
]
for dist in uni_dists

  single_sample_tests(dist)

  # specialised multi-sample tests.
  N = 10
  x = rand(dist)
  xs = fill(x, N)
  multi_sample_tests(dist, x, xs, N)
end

# Tests with vector-valued distributions.
vector_dists = [
  Dirichlet(2, 3),
  MvNormal(randn(10), exp.(randn(10))),
  MvLogNormal(MvNormal(randn(10), exp.(randn(10)))),
]
for dist in vector_dists

  single_sample_tests(dist)

  # Multi-sample tests. Columns are observations due to Distributions.jl conventions.
  N = 10
  x = rand(dist)
  xs = repeat(x, 1, N)
  multi_sample_tests(dist, x, xs, N)
end

# Tests with matrix-valued distributions.
matrix_dists = [
  Wishart(7, [1 0.5; 0.5 1]),
  InverseWishart(2, [1 0.5; 0.5 1])
]
for dist in matrix_dists

  single_sample_tests(dist)

  # Multi-sample tests comprising vectors of matrices.
  N = 10
  x = rand(dist)
  xs = [x for _ in 1:N]
  multi_sample_tests(dist, x, xs, N)
end

################################## Miscelaneous old tests ##################################

# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1000., -1000., -1000.]), true)
# NaN
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1000., -1000., -1000.], true, true)
# -1999.30685281944
#
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), exp.([-1., -2., -3.]), true)
# -3.006450206744678
# julia> logpdf_with_trans(Dirichlet([1., 1., 1.]), [-1., -2., -3.], true, true)
# -3.006450206744678
d  = Dirichlet([1., 1., 1.])
r  = [-1000., -1000., -1000.]
r2 = [-1., -2., -3.]

# test link
#link(d, r)

# test invlink
@test invlink(d, r) ≈ [0., 0., 1.] atol=1e-9

# test logpdf_with_trans
#@test logpdf_with_trans(d, invlink(d, r), true) -1999.30685281944 1e-9 ≈ # atol=NaN
@test logpdf_with_trans(d, invlink(d, r2), true) ≈ -3.760398892580863 atol=1e-9


text:
function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function free_slip_bc_nodal_displacement(
    mesh,
    basis,
    cellquads,
    facequads,
    interfacequads,
    stiffness,
    theta0,
    penalty;
    eta = 1,
)

    L, W = CutCellDG.widths(mesh)
    lambda1, mu1 = CutCellDG.lame_coefficients(stiffness, +1)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    sysmatrix = CutCellDG.SystemMatrix()
    sysrhs = CutCellDG.SystemRHS()

    CutCellDG.assemble_displacement_bilinear_forms!(
        sysmatrix,
        basis,
        cellquads,
        stiffness,
        mesh,
    )
    CutCellDG.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        mesh,
    )
    CutCellDG.assemble_interelement_condition!(
        sysmatrix,
        basis,
        facequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_interelement_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
    )

    CutCellDG.assemble_incoherent_interface_condition!(
        sysmatrix,
        basis,
        interfacequads,
        stiffness,
        mesh,
        penalty,
        eta,
    )
    CutCellDG.assemble_incoherent_interface_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        mesh,
    )

    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x->0.0,
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onleftboundary(x, L, W),
        [1.0,0.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
        x -> onleftboundary(x, L, W),
        [1.0,0.0]
    )

    CutCellDG.assemble_penalty_displacement_component_bc!(
        sysmatrix,
        sysrhs,
        x->0.0,
        basis,
        facequads,
        stiffness,
        mesh,
        x -> onbottomboundary(x, L, W),
        [0.0,1.0],
        penalty,
    )
    CutCellDG.assemble_penalty_displacement_component_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        mesh,
        x -> onbottomboundary(x, L, W),
        [0.0,1.0]
    )

    matrix = CutCellDG.sparse_displacement_operator(sysmatrix, mesh)
    rhs = CutCellDG.displacement_rhs_vector(sysrhs, mesh)

    solution = matrix \ rhs

    return solution
end

function construct_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    corner,
    numqp;
    tinyratio = 0.2,
)
    cgmesh = CutCellDG.CGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = CutCellDG.LevelSet(
        x -> -corner_distance_function(x, corner),
        cgmesh,
        basis,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset)
    cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
    interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
    facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

    mergedwithcell, hasmergedcells = CutCellDG.merge_tiny_cells_in_mesh!(
        cutmesh,
        cellquads,
        facequads,
        interfacequads,
        tinyratio = tinyratio,
    )
    mergedmesh = CutCellDG.MergedMesh!(cutmesh, mergedwithcell)

    return mergedmesh, cellquads, facequads, interfacequads
end

function construct_unmerged_mesh_and_quadratures(
    meshwidth,
    nelmts,
    basis,
    interfacecenter,
    interfaceradius,
    numqp,
)
    mesh = CutCellDG.DGMesh([0.0, 0.0], meshwidth, [nelmts, nelmts], basis)

    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCellDG.levelset_coefficients(
        x -> -circle_distance_function(x, interfacecenter, interfaceradius),
        mesh,
    )

    cutmesh = CutCellDG.CutMesh(mesh, levelset, levelsetcoeffs)
    cellquads =
        CutCellDG.CellQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    interfacequads =
        CutCellDG.InterfaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)
    facequads =
        CutCellDG.FaceQuadratures(cutmesh, levelset, levelsetcoeffs, numqp)

    return cutmesh, cellquads, facequads, interfacequads
end

function symmetric_displacement_gradient(
    celldisp,
    basis,
    point,
    jac,
    vectosymmconverter,
)
    dim = length(vectosymmconverter)
    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp
    return symmdispgrad
end

function parent_strain_at_reference_points(
    nodaldisplacement,
    basis,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert length(referencecellids) == numpts

    parentstrain = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    for i = 1:numpts
        cellid = referencecellids[i]
        nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[:, i]

        parentstrain[1:3, i] .= symmetric_displacement_gradient(
            celldisp,
            basis,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return parentstrain
end

function parent_stress(symmdispgrad, stiffness)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, -1)

    inplanestress = stiffness[-1] * symmdispgrad[1:3, :]
    s33 = lambda * (symmdispgrad[1, :] + symmdispgrad[2, :])

    return vcat(inplanestress, s33')
end

function product_strain_at_reference_points(
    nodaldisplacement,
    basis,
    theta0,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert length(referencecellids) == numpts

    productstrain = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    for i = 1:numpts
        cellid = referencecellids[i]
        nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[:, i]

        productstrain[1:3, i] .= symmetric_displacement_gradient(
            celldisp,
            basis,
            point,
            jac,
            vectosymmconverter,
        )
    end
    productstrain[[1, 2, 4], :] .-= theta0 / 3
    return productstrain
end

function product_stress(elasticstrain, stiffness, theta0)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, +1)

    inplanestress = stiffness[+1] * elasticstrain[1:3, :]
    inplanecorrection = lambda * elasticstrain[4, :]

    inplanestress[1, :] .+= inplanecorrection
    inplanestress[2, :] .+= inplanecorrection

    s33 =
        (lambda + 2mu) * elasticstrain[4, :] +
        lambda * (elasticstrain[1, :] + elasticstrain[2, :])

    return vcat(inplanestress,s33')
end

function product_stress(
    celldisp,
    basis,
    stiffness,
    transfstress,
    theta0,
    point,
    jac,
    vectosymmconverter,
)

    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, +1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = (stiffness[+1] * symmdispgrad) - transfstress
    s33 =
        lambda * (symmdispgrad[1] + symmdispgrad[2]) -
        (lambda + 2mu / 3) * theta0

    stress = vcat(inplanestress, s33)

    return stress
end

function product_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    transfstress,
    theta0,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert length(referencecellids) == numpts

    productstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    for i = 1:numpts
        cellid = referencecellids[i]
        nodeids = CutCellDG.nodal_connectivity(mesh, +1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[:, i]

        productstress[:, i] .= product_stress(
            celldisp,
            basis,
            stiffness,
            transfstress,
            theta0,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return productstress
end

function parent_stress(
    celldisp,
    basis,
    stiffness,
    point,
    jac,
    vectosymmconverter,
)
    dim = length(vectosymmconverter)
    lambda, mu = CutCellDG.lame_coefficients(stiffness, -1)

    grad = CutCellDG.transform_gradient(gradient(basis, point), jac)
    NK = sum([
        CutCellDG.make_row_matrix(vectosymmconverter[k], grad[:, k]) for
        k = 1:dim
    ])
    symmdispgrad = NK * celldisp

    inplanestress = stiffness[-1] * symmdispgrad
    s33 = lambda * (symmdispgrad[1] + symmdispgrad[2])

    stress = vcat(inplanestress, s33)

    return stress
end

function parent_stress_at_reference_points(
    nodaldisplacement,
    basis,
    stiffness,
    referencepoints,
    referencecellids,
    mesh,
)

    dim = CutCellDG.dimension(mesh)
    dim2, numpts = size(referencepoints)
    @assert dim == dim2
    @assert length(referencecellids) == numpts

    parentstress = zeros(4, numpts)
    vectosymmconverter = CutCellDG.vector_to_symmetric_matrix_converter()
    jac = CutCellDG.jacobian(mesh)

    for i = 1:numpts
        cellid = referencecellids[i]
        nodeids = CutCellDG.nodal_connectivity(mesh, -1, cellid)
        celldofs = CutCellDG.element_dofs(nodeids, dim)
        celldisp = nodaldisplacement[celldofs]
        point = referencepoints[:, i]

        parentstress[:, i] .= parent_stress(
            celldisp,
            basis,
            stiffness,
            point,
            jac,
            vectosymmconverter,
        )
    end
    return parentstress
end

function pressure(stress)
    return -1.0 / 3.0 * (stress[1, :] + stress[2, :] + stress[4, :])
end

function deviatoric_stress(stress, press)
    devstress = copy(stress)
    devstress[1, :] .+= press
    devstress[2, :] .+= press
    devstress[4, :] .+= press
    return devstress
end

function strain_energy(stress,strain)
    product = stress .* strain
    return 0.5*vec(sum(product,dims=1))
end


text:
using CartesianGrids
using Test

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Fields"
  include("fields.jl")
end

if GROUP == "All" || GROUP == "DDF"
  include("ddf.jl")
end

if GROUP == "All" || GROUP == "Points"
  include("points.jl")
end

if GROUP == "All" || GROUP == "GeneratedFields"
  include("generatedfields.jl")
end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
#######################################################################################################################################################################################################
# Changes to these constants
# General
#     2020-May-30: put the 2017 mat file in an artifact
#     2020-Aug-30: add the updated 2021 mat file along with that of 2017 into a new artifact
#
#######################################################################################################################################################################################################const FILE_SUN  = artifact"land_model_spectrum_V1" * "/sun.mat";
const OPTI_2017 = artifact"land_model_spectrum_V1" * "/Optipar2017_ProspectD.mat";
const OPTI_2021 = artifact"land_model_spectrum_V1" * "/Optipar2021_ProspectPRO_CX.mat";
WAVELENGTHS = [collect(400:10:650.1); collect(655:5:770.1); collect(780:25:2400.1)];

#######################################################################################################################################################################################################
# Changes to this structure
# General
#     2021-Aug-10: refactor the structure with renamed fields
#     2021-Aug-10: add a constructor within the structure to avoid external initialization
#     2021-Oct-19: sort variable to prognostic and dignostic catergories
#
#######################################################################################################################################################################################################
"""
$(TYPEDEF)

Structure that stores wave length information.

$(TYPEDFIELDS)

---
# Examples
```julia
wls = WaveLengthSet{FT}();
wls = WaveLengthSet{FT}(collect(FT,400:5:2500));
wls = WaveLengthSet{FT}(collect(FT,400:5:2500); opti=Emerald.OPTI_2017);
```
"""
mutable struct WaveLengthSet{FT<:AbstractFloat}
    # parameters that do not change with time
    "Indicies of Λ_NIR in Λ"
    IΛ_NIR::Vector{Int}
    "Indicies of Λ_PAR in Λ"
    IΛ_PAR::Vector{Int}
    "Indicies of Λ_SIF in Λ"
    IΛ_SIF::Vector{Int}
    "Indicies of Λ_SIFE in Λ"
    IΛ_SIFE::Vector{Int}
    "Number of wavelength bins"
    NΛ::Int
    "Number of wavelength bins for PAR"
    NΛ_PAR::Int
    "Number of wavelength bins for SIF"
    NΛ_SIF::Int
    "Number of wavelength bins for SIF excitation"
    NΛ_SIFE::Int
    "Wavelength limits for NIR `[nm]`"
    WL_NIR::Vector{FT}
    "Wavelength limits for PAR `[nm]`"
    WL_PAR::Vector{FT}
    "Wavelength limits for SIF emission `[nm]`"
    WL_SIF::Vector{FT}
    "Wavelength limits for SIF excitation `[nm]`"
    WL_SIFE::Vector{FT}
    "Standard wavelength (boundaries) `[nm]`"
    SΛ::Vector{FT}
    "Differential wavelength `[nm]`"
    ΔΛ::Vector{FT}
    "Differential wavelength for PAR `[nm]`"
    ΔΛ_PAR::Vector{FT}
    "Differential wavelength for SIF excitation `[nm]`"
    ΔΛ_SIFE::Vector{FT}
    "Wavelength (bins) `[nm]`"
    Λ::Vector{FT}
    "Wavelength bins for PAR `[nm]`"
    Λ_PAR::Vector{FT}
    "Wavelength bins for SIF `[nm]`"
    Λ_SIF::Vector{FT}
    "Wavelength bins for SIF excitation `[nm]`"
    Λ_SIFE::Vector{FT}

    # prognostic variables that change with time

    # dignostic variables that change with time

    # caches to speed up calculations

    # constructors
    WaveLengthSet{FT}(swl::Vector{FT}=FT.(WAVELENGTHS); opti::String=OPTI_2021) where {FT<:AbstractFloat} = (
        _dwl    = diff(swl);
        _λ      = zeros(FT, length(swl)-1);
        _opti   = matread(opti)["optipar"];
        _λ_opti = _opti["wl"];
        @inbounds for _i in 1:length(swl)-1
            _wo    = findall( swl[_i] .<= _λ_opti .< swl[_i+1] );
            _λ[_i] = nanmean(_λ_opti[_wo]);
        end;
        _iλ_nir  = findall( 700 .<= _λ .<= 2500 );
        _iλ_par  = findall( 400 .<= _λ .<= 700  );
        _iλ_sif  = findall( 640 .<= _λ .<= 850  );
        _iλ_sife = findall( 400 .<= _λ .<= 750  );

        return new(_iλ_nir, _iλ_par, _iλ_sif, _iλ_sife, length(_λ), length(_iλ_par), length(_iλ_sif), length(_iλ_sife), FT[700,2500], FT[400,700], FT[640,850], FT[400,750], swl, _dwl, _dwl[_iλ_par],
                   _dwl[_iλ_sife], _λ, _λ[_iλ_par], _λ[_iλ_sif], _λ[_iλ_sife])
    )
end


text:
mutable struct BaseTween
    target::Ranger.AbstractNode

    auto_start::Bool
    auto_remove::Bool

    # General
    step::Int64
    repeat_cnt::Int64
    is_iterationStep::Bool
    is_yoyo::Bool

    # Timings
    # the delay of the tween or timeline. Nothing will happen before
    delay::Int64

    # the duration of a single iteration.
    duration::Int64
    repeat_delay::Int64
    current_time::Int64 #elapsed time, relative to current iteration
    delta_time::Int64
    started::Bool # true when the object is started
    initialized::Bool # true after the delay
    finished::Bool # true when all repetitions are done
    killed::Bool # true if kill() was called
    paused::Bool # true if pause() was called

    state_callback::Function
    callback_triggers::UInt64 # Bit flags

    function BaseTween()
        o = new()
        reset!(o)
        o
    end
end

function start!(base::BaseTween)
    base.current_time = 0
    base.started = true
end

function reset!(base::BaseTween)
    base.step = -2
    base.repeat_cnt = 0
    base.is_iterationStep = false
    base.is_yoyo = false

    base.delay = 0

    base.duration = 0
    base.repeat_delay = 0
    base.current_time = 0
    base.delta_time = 0
    base.started = false
    base.initialized = false
    base.finished = false
    base.killed = false
    base.paused = false

    # base.state_callback = nothing
    base.callback_triggers = TWEEN_CALLBACK_COMPLETE

    base.auto_start = true
    base.auto_remove = true
end

# Kills the tween or timeline. If you are using a [TweenManager],
# this object will be removed automatically.
function kill(base::BaseTween)
    base.killed = true
end

# Returns the complete duration, including initial delay and repetitions.
function full_duration(base::BaseTween)
    if base.repeat_cnt < 0
        return -1
    else
        return base.delay + base.duration + (base.repeat_delay + base.duration) * base.repeat_cnt
    end
end

function force_to_start!(base::BaseTween, tween::AbstractTween)
    base.current_time = -base.delay
    base.step = -1
    base.is_iterationStep = false

    if is_reverse(base, 0)
        force_end_values!(tween)
    else
        force_start_values!(tween)
    end
end

function force_to_end!(base::BaseTween, tween::AbstractTween)
    base.current_time = tween.base.duration - full_duration(base)
    base.step = base.repeat_cnt * 2 + 1
    base.is_iterationStep = false

    if is_reverse(base, base.repeat_cnt * 2)
        force_start_values!(tween)
    else
        force_end_values!(tween)
    end
end

function is_reverse(base::BaseTween, step::Int64)
    base.is_yoyo && abs(base.step % 4) == 2
end

function is_valid(base::BaseTween, step::Int64)
    (step >= 0 && step <= base.repeat_cnt * 2) || base.repeat_cnt < 0
end

function kill_target()

end

# Updates the tween or timeline state. **You may want to use a
# TweenManager to update objects for you.**
#
# Slow motion, fast motion and backward play can be easily achieved by
# tweaking this delta time. Multiply it by -1 to play the animation
# backward, or by 0.5 to play it twice slower than its normal speed.
#
# [delta] The time between now and the last call, in seconds.
function update!(base::BaseTween, delta::Float64)
    println("BaseTween updating")

    if (!base.started || base.paused || base.killed)
        return
    end

    base.delta_time = delta

    if !base.initialized
        initialize!(base)
    end
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
"""
    stabilize_denom(d, [eps = 1f-9])

Replace zero terms of a matrix `d` with `eps`.
"""
function stabilize_denom(d::T, eps=T(1.0f-9)) where {T}
    iszero(d) && (return T(eps))
    return d + sign(d) * T(eps)
end
stabilize_denom(D::AbstractArray{T}, eps=T(1.0f-9)) where {T} = stabilize_denom.(D, eps)

"""
    safedivide(a, b, [eps = 1f-6])

Elementwise division of two matrices avoiding near zero terms
in the denominator by replacing them with `± eps`.
"""
function safedivide(a::AbstractArray{T}, b::AbstractArray{T}, eps=T(1.0f-9)) where {T}
    return a ./ stabilize_denom(b, T(eps))
end

"""
    batch_dim_view(A)

Return a view onto the array `A` that contains an extra singleton batch dimension at the end.
This avoids allocating a new array.

## Example
```juliarepl
julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> batch_dim_view(A)
2×2×1 view(::Array{Int64, 3}, 1:2, 1:2, :) with eltype Int64:
[:, :, 1] =
 1  2
 3  4
```
"""
batch_dim_view(A::AbstractArray{T,N}) where {T,N} = view(A, ntuple(Returns(:), N + 1)...)

"""
    drop_batch_index(I)

Drop batch dimension index (last value) from CartesianIndex.

## Example
julia> drop_batch_index(CartesianIndex(5,3,2))
CartesianIndex(5, 3)
"""
drop_batch_index(C::CartesianIndex) = CartesianIndex(C.I[1:(end - 1)])

# Utils for printing model check summary using PrettyTable.jl
_print_name(layer) = "$layer"
_print_name(layer::Parallel) = "Parallel(...)"
_print_activation(layer) = hasproperty(layer, :σ) ? "$(layer.σ)" : "—"
_print_activation(layer::Parallel) = "—"

function _show_check_summary(
    c::Chain, layer_names, layer_checks, activation_names, activation_checks
)
    hl_pass = Highlighter((data, i, j) -> j in (3, 5) && data[i, j]; foreground=:green)
    hl_fail = Highlighter((data, i, j) -> j in (3, 5) && !data[i, j]; foreground=:red)
    data = hcat(
        collect(1:length(c)),
        layer_names,
        collect(layer_checks),
        activation_names,
        collect(activation_checks),
    )
    pretty_table(
        data;
        header=["", "Layer", "Layer supported", "Activation", "Act. supported"],
        alignment=[:r, :l, :r, :c, :r],
        highlighters=(hl_pass, hl_fail),
    )
    return nothing
end


text:
using NPZ

# We expects the data to be in npy format, return a dict of {group: items}, each file is a different group
function load_data(path::String; prefix::String="", swapDimension::Bool = true)
    groups_dict = Dict()
    arr = npzread(path * prefix * ".npy")
    for (index, value) in enumerate(arr)
        if isnan(value)
            arr[index] = 0.0
        end
    end
    return swapDimension ? transpose(arr) : arr
end

# #Note that we expect the log_likelihood_array to be in rows (samples) x columns (clusters) , this is due to making it more efficent that way.
function sample_log_cat_array!(labels::AbstractArray{Int64,1}, log_likelihood_array::AbstractArray{Float32,2})
    # println("lsample log cat" * string(log_likelihood_array))
    log_likelihood_array[isnan.(log_likelihood_array)] .= -Inf #Numerical errors arent fun
    max_log_prob_arr = maximum(log_likelihood_array, dims = 2)
    log_likelihood_array .-= max_log_prob_arr
    map!(exp,log_likelihood_array,log_likelihood_array)
    # println("lsample log cat2" * string(log_likelihood_array))
    sum_prob_arr = sum(log_likelihood_array, dims =[2])
    log_likelihood_array ./=  sum_prob_arr
    for i=1:length(labels)
        labels[i] = sample(1:size(log_likelihood_array,2), ProbabilityWeights(log_likelihood_array[i,:]))
    end
end

function create_sufficient_statistics(dist::distribution_hyper_params, pts::Array{Any,1})
    return create_sufficient_statistics(dist,dist, Array{Float32}(undef, 0, 0))
end

function get_labels_histogram(labels)
    hist_dict = Dict()
    for v in labels
        if haskey(hist_dict,v) == false
            hist_dict[v] = 0
        end
        hist_dict[v] += 1
    end
    return sort(collect(hist_dict), by=x->x[1])
end

function get_node_leaders_dict()
    leader_dict = Dict()
    cur_leader = (nworkers()== 0 ? procs() : workers())[1]
    leader_dict[cur_leader] = []
    for i in (nworkers()== 0 ? procs() : workers())
        if i in procs(cur_leader)
            push!(leader_dict[cur_leader], i)
        else
            cur_leader = i
            leader_dict[cur_leader] = [i]
        end
    end
    return leader_dict
end

function log_multivariate_gamma(x::Number, D::Number)
    res::Float32 = D*(D-1)/4*log(pi)
    for d = 1:D
        res += logabsgamma(x+(1-d)/2)[1]
    end
    return res
end

function dcolwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    n = length(r)
    @simd for j = 1:n
        v = zero(promote_type(eltype(a), eltype(b)))
        @simd for i = 1:size(a, 1)
            @fastmath @inbounds v += a[i, j]*b[i, j]
        end
        r[j] = v
    end
end

function suff_stats_aggregation(suff_statistics_l::Vector{Tuple{sufficient_statistics,Number}},suff_statistics_r::Vector{Tuple{sufficient_statistics,Number}})
    max_age = maximum(vcat([x[2] for x in suff_statistics_l],[x[2] for x in suff_statistics_r]))
    index_l,index_r = 1,1
    cur_time = 0
    suff_stats = Vector{Tuple{sufficient_statistics,Number}}()
    if length(suff_statistics_l) != length(suff_statistics_r)
        println("OH NOES THE DEVIL IS UPON US")
    end
    while cur_time <= max_age
        ss_vector = []
        if index_l <= length(suff_statistics_l) && suff_statistics_l[index_l][2] == cur_time
            push!(ss_vector,suff_statistics_l[index_l][1])
            index_l += 1
        end
        if index_r <= length(suff_statistics_r) && suff_statistics_r[index_r][2] == cur_time
            push!(ss_vector,suff_statistics_r[index_r][1])
            index_r += 1
        end
        if length(ss_vector) > 0
            push!(suff_stats,(reduce(aggregate_suff_stats, ss_vector),cur_time))
        end
        cur_time = min(index_l <= length(suff_statistics_l) ? suff_statistics_l[index_l][2] : max_age+1,
                        index_r <= length(suff_statistics_r) ? suff_statistics_r[index_r][2] : max_age+1)
    end
    if length(suff_stats) == 0 && length(suff_statistics_l) != 0 && length(suff_statistics_r) != 0
        println([x[2] for x in suff_statistics_l],[x[2] for x in suff_statistics_r])
    end
    return suff_stats
end

# function suff_stats_aggregation(suff_statistics_l::Vector{Tuple{sufficient_statistics,Number}},suff_statistics_r::Vector{Tuple{sufficient_statistics,Number}})
#     return vcat(suff_statistics_l,suff_statistics_r)
# end

# function suff_stats_aggregation(suff_statistics_l::Vector{Tuple{sufficient_statistics,Number}},suff_statistics_r::Vector{Tuple{sufficient_statistics,Number}})
#     max_age = maximum(vcat([x[2] for x in suff_statistics_l],[x[2] for x in suff_statistics_r]))
#     index_l,index_r = 1,1
#     cur_time = 0
#     suff_stats = Vector{Tuple{sufficient_statistics,Number}}()

#     for i=1:length(suff_statistics_l)
#         push!(suff_stats,(aggregate_suff_stats(suff_statistics_l[i][1],suff_statistics_r[i][1]),suff_statistics_l[i][2]))
#     end

#     return suff_stats
# end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

type IntSet <: AbstractSet{Int}
    bits::Array{UInt32,1}
    limit::Int
    fill1s::Bool

    IntSet() = new(zeros(UInt32,256>>>5), 256, false)
end
IntSet(itr) = (s=IntSet(); for a in itr; push!(s,a); end; s)

eltype(::Type{IntSet}) = Int64
similar(s::IntSet) = IntSet()

function show(io::IO, s::IntSet)
    print(io, "IntSet([")
    first = true
    for n in s
        if n > s.limit
            break
        end
        if !first
            print(io, ", ")
        end
        print(io, n)
        first = false
    end
    if s.fill1s
        print(io, ", ..., ", typemax(Int)-1)
    end
    print(io, "])")
end

copy(s::IntSet) = union!(IntSet(), s)

function sizehint!(s::IntSet, top::Integer)
    if top >= s.limit
        lim = ((top+31) & -32)>>>5
        olsz = length(s.bits)
        if olsz < lim
            resize!(s.bits, lim)
            fill = s.fill1s ? UInt32(-1) : UInt32(0)
            for i=(olsz+1):lim; s.bits[i] = fill; end
        end
        s.limit = top
    end
    s
end

function push!(s::IntSet, n::Integer)
    if n >= s.limit
        if s.fill1s
            return s
        else
            lim = Int(n + div(n,2))
            sizehint!(s, lim)
        end
    elseif n <= 0
        if n < 0
            throw(ArgumentError("IntSet elements cannot be negative"))
        else
            depwarn("storing zero in IntSets is deprecated", :push!)
        end
    end
    s.bits[n>>5 + 1] |= (UInt32(1)<<(n&31))
    return s
end

function union!(s::IntSet, ns)
    for n in ns
        push!(s, n)
    end
    return s
end

function pop!(s::IntSet, n::Integer, deflt)
    if n >= s.limit
        if s.fill1s
            lim = Int(n + div(n,2))
            sizehint!(s, lim)
        else
            return deflt
        end
    end
    if n <= 0
        if n < 0
            return deflt
        else
            depwarn("stored zeros in IntSet is deprecated", :pop!)
        end
    end
    mask = UInt32(1)<<(n&31)
    idx = n>>5 + 1
    b = s.bits[idx]
    if (b&mask)==0; return deflt; end
    s.bits[idx] = b&~mask
    return n
end

function pop!(s::IntSet, n::Integer)
    if pop!(s, n, n+1) == n+1
        throw(KeyError(n))
    end
    return n
end

# TODO: what should happen when fill1s == true?
pop!(s::IntSet) = pop!(s, last(s))

function delete!(s::IntSet, n::Integer)
    pop!(s, n, n)
    return s
end

function setdiff!(s::IntSet, ns)
    for n in ns
        delete!(s, n)
    end
    return s
end

setdiff(a::IntSet, b::IntSet) = setdiff!(copy(a),b)
symdiff(s1::IntSet, s2::IntSet) =
    (s1.limit >= s2.limit ? symdiff!(copy(s1), s2) : symdiff!(copy(s2), s1))

function empty!(s::IntSet)
    fill!(s.bits, 0)
    return s
end

"""
    symdiff!(s, n)

The set `s` is destructively modified to toggle the inclusion of integer `n`.
"""
function symdiff!(s::IntSet, n::Integer)
    if n >= s.limit
        lim = Int(n + div(n,2))
        sizehint!(s, lim)
    elseif n < 0
        throw(ArgumentError("IntSet elements cannot be negative"))
    end
    s.bits[n>>5 + 1] ⊻= UInt32(1)<<(n&31)
    return s
end

"""
    symdiff!(s, itr)

For each element in `itr`, destructively toggle its inclusion in set `s`.
"""
function symdiff!(s::IntSet, ns)
   for n in ns
       symdiff!(s, n)
   end
   return s
end

function copy!(to::IntSet, from::IntSet)
    if to === from
        return to
    else
        empty!(to)
        return union!(to, from)
    end
end

in(n, s::IntSet) = n < 0 ? false : (n > typemax(Int) ? s.fill1s : in(convert(Int, n), s))
function in(n::Int, s::IntSet)
    if n >= s.limit
        # max IntSet length is typemax(Int), so highest possible element is
        # typemax(Int)-1
        return s.fill1s && n >= 0 && n < typemax(Int)
    elseif n <= 0
        if n < 0
            return false
        else
            depwarn("stored zeros in IntSet is deprecated", :in)
        end
    end
    (s.bits[n>>5 + 1] & (UInt32(1)<<(n&31))) != 0
end

start(s::IntSet) = Int64(0)
done(s::IntSet, i) = (!s.fill1s && next(s,i)[1] >= s.limit) || i == typemax(Int)
function next(s::IntSet, i)
    if i >= s.limit
        n = Int64(i)
    else
        n = Int64(ccall(:bitvector_next, UInt64, (Ptr{UInt32}, UInt64, UInt64), s.bits, i, s.limit))
    end
    (n, n+1)
end

isempty(s::IntSet) =
    !s.fill1s && ccall(:bitvector_any1, UInt32, (Ptr{UInt32}, UInt64, UInt64), s.bits, 0, s.limit)==0

function first(s::IntSet)
    n = next(s,0)[1]
    if n >= s.limit
        throw(ArgumentError("set must be non-empty"))
    end
    return n
end

shift!(s::IntSet) = pop!(s, first(s))

function last(s::IntSet)
    if !s.fill1s
        for i = length(s.bits):-1:1
            w = s.bits[i]
            if w != 0
                return (i-1)<<5 + (31-leading_zeros(w))
            end
        end
    end
    throw(ArgumentError("set has no last element"))
end

length(s::IntSet) = Int(ccall(:bitvector_count, UInt64, (Ptr{UInt32}, UInt64, UInt64), s.bits, 0, s.limit)) +
    (s.fill1s ? typemax(Int) - s.limit : 0)

# Math functions
function union!(s::IntSet, s2::IntSet)
    if s2.limit > s.limit
        sizehint!(s, s2.limit)
    end
    lim = length(s2.bits)
    for n = 1:lim
        s.bits[n] |= s2.bits[n]
    end
    if s2.fill1s
        for n=lim+1:length(s.bits)
            s.bits[n] = UInt32(-1)
        end
    end
    s.fill1s |= s2.fill1s
    s
end

union(s1::IntSet) = copy(s1)
union(s1::IntSet, s2::IntSet) = (s1.limit >= s2.limit ? union!(copy(s1), s2) : union!(copy(s2), s1))
union(s1::IntSet, ss::IntSet...) = union(s1, union(ss...))

"""
    intersect!(s1::IntSet, s2::IntSet)

Intersects sets `s1` and `s2` and overwrites the set `s1` with the result. If needed, `s1`
will be expanded to the size of `s2`.
"""
function intersect!(s::IntSet, s2::IntSet)
    if s2.limit > s.limit
        sizehint!(s, s2.limit)
    end
    lim = length(s2.bits)
    for n = 1:lim
        s.bits[n] &= s2.bits[n]
    end
    if !s2.fill1s
        for n=lim+1:length(s.bits)
            s.bits[n] = UInt32(0)
        end
    end
    s.fill1s &= s2.fill1s
    s
end

intersect(s1::IntSet) = copy(s1)
intersect(s1::IntSet, s2::IntSet) =
    (s1.limit >= s2.limit ? intersect!(copy(s1), s2) : intersect!(copy(s2), s1))
intersect(s1::IntSet, ss::IntSet...) = intersect(s1, intersect(ss...))

"""
    symdiff!(s1, s2)

Construct the symmetric difference of sets `s1` and `s2`, storing the result in `s1`.
"""
function symdiff!(s::IntSet, s2::IntSet)
    if s2.limit > s.limit
        sizehint!(s, s2.limit)
    end
    lim = length(s2.bits)
    for n = 1:lim
        s.bits[n] ⊻= s2.bits[n]
    end
    if s2.fill1s
        for n=lim+1:length(s.bits)
            s.bits[n] = ~s.bits[n]
        end
    end
    s.fill1s ⊻= s2.fill1s
    s
end

function ==(s1::IntSet, s2::IntSet)
    if s1.fill1s != s2.fill1s
        return false
    end
    lim1 = length(s1.bits)
    lim2 = length(s2.bits)
    for i = 1:min(lim1,lim2)
        if s1.bits[i] != s2.bits[i]
            return false
        end
    end
    filln = s1.fill1s ? reinterpret(UInt32, Int32(-1)) : UInt32(0)
    if lim1 > lim2
        for i = lim2:lim1
            if s1.bits[i] != filln
                return false
            end
        end
    else
        for i = lim1+1:lim2
            if s2.bits[i] != filln
                return false
            end
        end
    end
    return true
end

const hashis_seed = UInt === UInt64 ? 0x88989f1fc7dea67d : 0xc7dea67d
function hash(s::IntSet, h::UInt)
    h += hashis_seed
    h += hash(s.fill1s)
    filln = s.fill1s ? ~zero(eltype(s.bits)) : zero(eltype(s.bits))
    for x in s.bits
        if x != filln
            h = hash(x, h)
        end
    end
    return h
end

issubset(a::IntSet, b::IntSet) = isequal(a, intersect(a,b))
<(a::IntSet, b::IntSet) = (a<=b) && !isequal(a,b)
<=(a::IntSet, b::IntSet) = issubset(a, b)

max_stars_repo_name: oxinabox/julia

max_issues_repo_name: oxinabox/julia

max_forks_repo_name: oxinabox/julia