text:
export logsumexp!, softmax3!, logsumexp_and_softmax!

"""
    sumprod(f, x, y, u, v) = sum(x .* y .* f.(u, v))
"""
function sumprod(f::Function, x::AbstractArray{T}, y::AbstractArray{T}, u::AbstractArray{T}, v::AbstractArray{T}) where {T<:Real}
    length(x)==length(y)==length(u)==length(v) || throw(DimensionMismatch())
    s = zero(T)
    @inbounds @simd for i in eachindex(x)
        s += x[i] * y[i] * f(u[i], v[i])
    end
    return s
end

"""
    sumprod3(x, y, z)

sumprod3 = sum( x .* y .* z)
"""
function sumprod3(x::AbstractArray{T}, y::AbstractArray{T}, z::AbstractArray{T}) where {T<:Real}
    length(x)==length(y)==length(z) || throw(DimensionMismatch())
    s = zero(T)
    @inbounds @simd for i in eachindex(x)
        s += x[i] * y[i] * z[i]
    end
    return s
end

function sumprod!(red::AbstractArray3, big::AbstractArray4, small::AbstractArray3)
     nz, nψ, nk, nd = size(big)
    (nz, nψ,     nd) == size(small) || throw(DimensionMismatch())
    (nz, nψ, nk,   ) == size(red)   || throw(DimensionMismatch())

    # second set w/ plus equals
    @inbounds for d in 1:nd, k in 1:nk
        if d == 1
            @views red[:,:,k] .=  small[:,:,d] .* big[:,:,k,d]
        else
            @views red[:,:,k] .+= small[:,:,d] .* big[:,:,k,d]
        end
    end
end

function sumprod!(red::AbstractMatrix, big::AbstractArray3, small::AbstractArray3)
     nz, nψ, nd  = size(big)
    (nz, nψ, nd) == size(small) || throw(DimensionMismatch())
    (nz, nψ)     == size(red) || throw(DimensionMismatch())

    fill!(red, 0)
    @inbounds for d in 1:nd
        @views red .+= small[:,:,d] .* big[:,:,d]
    end
end

# TODO would be good to evaluate whether I can do better with tmapreduce
# https://github.com/jw3126/ThreadingTools.jl
# https://discourse.julialang.org/t/how-to-speed-up-this-simple-code-multithreading-simd-inbounds/19681/39
# https://discourse.julialang.org/t/parallelizing-for-loop-in-the-computation-of-a-gradient/9154/7?u=tkoolen
# https://discourse.julialang.org/t/innefficient-paralellization-need-some-help-optimizing-a-simple-dot-product/9723/31
# https://discourse.julialang.org/t/for-loops-acceleration/24011/7

# tmapreduce((q,c,u,v) -> q*c*dψdρ(u,v), +, qm,cm,um,vm; init=zero(eltype(qm))) * βψ

# https://discourse.julialang.org/t/poor-performance-on-cluster-multithreading/12248/39

"""
    logsumexp!(r, x)

Compute `r` = softmax(x) and return `logsumexp(x)`.

Based on code from
https://arxiv.org/pdf/1412.8695.pdf eq 3.8 for p(y)
https://discourse.julialang.org/t/fast-logsumexp/22827/7?u=baggepinnen for stable logsumexp
"""
function logsumexp!(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
    n = length(x)
    length(r) == n || throw(DimensionMismatch())
    isempty(x) && return -T(Inf)
    1 == stride1(r) == stride1(x) || throw(error("Arrays not strided"))

    u = maximum(x)                                       # max value used to re-center
    abs(u) == Inf && return any(isnan, x) ? T(NaN) : u   # check for non-finite values

    s = zero(T)
    for i = 1:n
        tmp = exp(x[i] - u)
        r[i] = tmp
        s += tmp
    end

    invs = inv(s)
    r .*= invs

    return log1p(s-1) + u
end

logsumexp!(x) = logsumexp!(x,x)

add_1_dim(x::AbstractArray) = reshape(x, size(x)..., 1)

"""
    softmax3!(q, lse, tmp, x, maxk)

set `q[i,j,k] = Pr( k | x[i,j,:])` for k in 1:maxk

Sets `lse = ∑_k exp(x[i,j,k])`

Uses temporary array `tmp`
"""
function softmax3!(q::AA, lse::A, tmpmax::A, x::AA, maxk=size(q, ndims(q)) ) where {T<:Real, A<:Array{T}, AA<:AbstractArray{T}}
        ndims(q) == 1+ndims(lse) || throw(DimensionMismatch())
        xsizes = size(x)
        xsizes == size(q) || throw(DimensionMismatch("size(x) = $(size(x)) but size(q) = $(size(q))"))
        nk = last(xsizes)
        for i = OneTo(ndims(lse))
            size(q,i) == size(lse,i) == size(tmpmax,i) || throw(DimensionMismatch("size(x) = $(size(x)),  size(lse) = $(size(lse)), and size(tmpmax) = $(size(tmpmax))"))
        end
        0 < maxk <= nk || throw(DomainError(maxk))
        1 == stride1(q) == stride1(x) || throw(error("Arrays not strided"))

        isempty(x) && throw(error("x empty"))
        check_finite(x)

        maximum!(add_1_dim(tmpmax), x)
        fill!(lse, zero(T))

        xx = reshape(x, :, nk)
        qq = reshape(q, :, nk)
        tmpmaxvec = vec(tmpmax)
        lsevec = vec(lse)

        for k in OneTo(maxk)
            qk = uview(qq, :, k)
            xk = uview(xx, :, k)
            for i in eachindex(lsevec)
                tmp = exp(xk[i] - tmpmaxvec[i])
                lsevec[i] += tmp
                qk[i] = tmp
            end
        end

        for k in maxk+1:nk
            qk = uview(qq, :, k)
            xk = uview(xx, :, k)
            for i in eachindex(lsevec)
                tmp = exp(xk[i] - tmpmaxvec[i])
                lsevec[i] += tmp
            end
        end

        qq[:,OneTo(maxk)] ./= vec(lse)
end

"use in vfit!"
function logsumexp_and_softmax!(lse, q, tmpmax, x, args...)
    softmax3!(q, lse, tmpmax, x, args...)
    lse .-= 1
    lse .= log1p.(lse) .+ tmpmax
end

function softmax3test!(q::AbstractArray3, x::AbstractArray3)
    ni, nj, nk = size(x)
    (ni, nj, nk) == size(q) || throw(DimensionMismatch())
    @inbounds for j = OneTo(nj)
        for i in OneTo(ni)
            @views softmax!(q[i,j,:], x[i,j,:])
        end
    end
end

function logsumexp3test!(x::AbstractArray3{T}) where {T}
    ni, nj, nk = size(x)
    lse = zeros(T,ni,nj)
    tmp = zeros(T,nk)
    # @inbounds
    for j = OneTo(nj)
        for i in OneTo(ni)
            tmp .= view(x, i, j, :)
            lse[i,j] = logsumexp!(tmp)
        end
    end
    return lse
end

@inline function sumsq(x, start::Integer, stop::Integer)
    @assert start >= 1
    @assert stop <= length(x)
    s = zero(eltype(x))
    @inbounds @simd for i in start:stop
        s += x[i]*x[i]
    end
    return s
end


text:
using DataFrames, Dates, JSON, ThreadsX

### Function to list all files with specific pattern in directory recursively
function listFiles(dir::String, pattern::String )
	f = String[]
	for (root, dirs, files) in walkdir(dir)
		for file in files
			if split(file, ".")[end] == pattern push!(f, joinpath(root,file)) end
		end
	end
return(f)
end

### get files
files = listFiles("./", "json")

### print options for read in's
println("Read In Options:
	heart_rate
	altitude
	calories
	demographicVO2Max
	distance
	lightly_active_minutes
	moderately_active_minutes
	very_active_minutes
	resting_heart_rate
	run_vo2_max
	sedentary_minutes
	steps
	swim_lengths_data
	time_in_heart_rate_zones
	exercise
	sleep")

##############################################################################################
function read_in(inp::String, files=files, ree=ree)
	f = files[occursin.(Regex(string("/",inp)),files)]
	h = reduce(vcat, ThreadsX.collect(ree(i) for i in f))
	h = DataFrame(h)
	h[!,:dateTime] = try
				DateTime.(h.dateTime, dateformat"mm/dd/yy HH:MM:SS") .+ Dates.Year(2000)
			catch
				 DateTime.(h.dateTime, dateformat"yyyy-mm-ddTHH:MM:SS")
		 	end
  return(h)
end

#################################################################################################
function ree(x::String)
        dat = JSON.parse(read(x, String))

	if occursin(r"/heart_rate", x) == true
		f= Array{hr}(undef,length(dat))
		@simd for i in 1:length(dat)
            		@inbounds f[i] = hr(dat[i]["dateTime"], dat[i]["value"]["bpm"], dat[i]["value"]["confidence"])
		end

	elseif occursin("altitude", x) == true
        	f= Array{alt}(undef,length(dat))
        	[f[i] = alt(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin("calories", x) == true
        	f= Array{cal}(undef,length(dat))
        	[f[i] = cal(dat[i]["dateTime"], parse(Float64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin(r"/demographic", x) == true
        	f= Array{demVM}(undef,length(dat))
        	for i in 1:length(dat)
                	f[i] = demVM(dat[i]["dateTime"],
				 dat[i]["value"]["demographicVO2MaxError"],
			 	dat[i]["value"]["demographicVO2Max"],
			 	dat[i]["value"]["filteredDemographicVO2MaxError"],
			 	dat[i]["value"]["filteredDemographicVO2Max"] )
        	end
	elseif occursin("distance", x) == true
        	f= Array{dist}(undef,length(dat))
        	[f[i] = dist(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin("lightly_active_minutes" , x) == true
        	f= Array{lam}(undef,length(dat))
        	[f[i] = lam(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin("moderately_active_minutes", x) == true
        	f= Array{mam}(undef,length(dat))
        	[f[i] = mam(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin("very_active_minutes", x) == true
        	f= Array{vam}(undef,length(dat))
        	[f[i] = vam(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin(r"/resting_heart_rate*", x) == true
		f= Array{rest_hr}(undef,length(dat))
       		[f[i] = rest_hr(dat[i]["dateTime"], dat[i]["value"]["value"], dat[i]["value"]["error"]) for i in 1:length(dat)]
	elseif occursin("run_vo2_max", x) == true
        	f= Array{run_vo2}(undef,length(dat))
        	for i in 1:length(dat)
                	f[i] = run_vo2(dat[i]["dateTime"],
                        	 dat[i]["value"]["exerciseId"],
                         	dat[i]["value"]["runVO2MaxError"],
                         	dat[i]["value"]["filteredRunVO2MaxError"],
                         	dat[i]["value"]["filteredRunVO2Max"],
				 dat[i]["value"]["runVO2Max"] )
        	end
	elseif occursin("sedentary_minutes", x) == true
        	f= Array{sed_min}(undef,length(dat))
        	[f[i] = sed_min(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]

	elseif occursin("steps", x) == true
       		f= Array{step}(undef,length(dat))
        	[f[i] = step(dat[i]["dateTime"], parse(Int64,values(dat[i]["value"]))) for i in 1:length(dat)]
	elseif occursin("swim_lengths_data", x) == true
		f= Array{swim}(undef,length(dat))
        	for i in 1:length(dat)
                	f[i] = swim(dat[i]["dateTime"],
                        	 dat[i]["value"]["strokeCount"],
                         	dat[i]["value"]["lapDurationSec"],
                         	dat[i]["value"]["swimAlgorithmType"],
                         	dat[i]["value"]["swimStrokeType"])
        	end
	elseif occursin("time_in_heart_rate_zones", x) == true
		f= Array{tihrz}(undef,length(dat))
        	for i in 1:length(dat)
                	f[i] = tihrz(dat[i]["dateTime"],
                        	dat[i]["value"]["valuesInZones"]["IN_DEFAULT_ZONE_2"],
			 	dat[i]["value"]["valuesInZones"]["BELOW_DEFAULT_ZONE_1"],
                         	dat[i]["value"]["valuesInZones"]["IN_DEFAULT_ZONE_1"],
                         	dat[i]["value"]["valuesInZones"]["IN_DEFAULT_ZONE_3"])
        	end
	elseif occursin("exercise", x) == true
        	f= Array{exercs}(undef,length(dat))
        	for i in 1:length(dat)
			f[i] = exercs(dat[i]["startTime"],
			      	 	dat[i]["calories"],
				 	dat[i]["logType"],
				 	try
				 		dat[i]["steps"]
			 	 	catch
				 		missing
			 	 	end,
				 		dat[i]["manualValuesSpecified"]["calories"],
				 		dat[i]["manualValuesSpecified"]["distance"],
				 		dat[i]["manualValuesSpecified"]["steps"],
				 	try
				 		dat[i]["distanceUnit"]
			 	 	catch
			  	 		missing
			 	 	end,
				 		dat[i]["hasActiveZoneMinutes"],
				 		dat[i]["activeDuration"],
				 		dat[i]["duration"],
				 	try
					 	dat[i]["vo2Max"]["vo2Max"]
				 	catch
						missing
				 	end,
				 	try
			   	 		dat[i]["elevationGain"]
			 	 	catch
				 		missing
			 	 	end,
				 	try
				 		dat[i]["distance"]
			 	 	catch
				 		missing
			 	 	end,
				 		dat[i]["activityName"],
					 try
				  		dat[i]["pace"]
				 	catch
				 		missing
				 	end,
				 		dat[i]["originalDuration"],
				 		dat[i]["shouldFetchDetails"],
				 	try
				 		dat[i]["heartRateZones"][1]["minutes"]
			 	 	catch
				 		missing
			 	 	end,
				 	try
				 		dat[i]["heartRateZones"][2]["minutes"]
			 	 	catch
				 		missing
			 	 	end,
				 	try
				 		dat[i]["heartRateZones"][3]["minutes"]
			 	 	catch
				 		missing
			 	 	end,
				 	try
				 		dat[i]["heartRateZones"][4]["minutes"]
			 	 	catch
				 		missing
			 	 	end,
				 		dat[i]["logId"],
				 		dat[i]["activityTypeId"],
				 	try
				 		dat[i]["averageHeartRate"]
			 	 	catch
				 		missing
			 	 	end,
				 		dat[i]["hasGps"],
				 		dat[i]["activityLevel"][1]["minutes"],
				 		dat[i]["activit,yLevel"][2]["minutes"],
				 	dat[i]["activityLevel"][3]["minutes"],
				 	dat[i]["activityLevel"][4]["minutes"],
				 	try
				 		dat[i]["speed"]
			 	 	catch
				 		missing
			 	 	end,
				 		dat[i]["lastModified"])
        	end
	elseif occursin("sleep", x) == true

		f= Array{sle}(undef,length(dat))
		for i in 1:length(dat)
			f[i] = sle(dat[i]["startTime"][1:19],
			   dat[i]["endTime"],
			   dat[i]["minutesAsleep"],
			   dat[i]["minutesToFallAsleep"],
			   dat[i]["minutesAwake"],
			   dat[i]["duration" ],
			   dat[i]["timeInBed"],
			   dat[i]["mainSleep" ],
			   dat[i]["infoCode"],
			   dat[i]["minutesAfterWakeup"],
			   dat[i]["efficiency"],
			   dat[i]["type"],
			   try
			   	dat[i]["levels"]["summary"]["awake"]["minutes"]
		   	   catch
				dat[i]["levels"]["summary"]["wake"]["minutes"]
		   	   end,
			   try
			   	dat[i]["levels"]["summary"]["asleep"]["minutes"]
			   catch
		           	missing
			   end,
			   try
			   	dat[i]["levels"]["summary"]["restless"]["minutes"]
			   catch
				missing
			   end,
			   try
				dat[i]["levels"]["summary"]["light"]["minutes"]
			   catch
				missing
			  end,
			  try
				dat[i]["levels"]["summary"]["deep"]["minutes"]
			 catch
			 	missing
			 end,
			 try
			       dat[i]["levels"]["summary"]["rem"]["minutes"]
			 catch
			      missing
		      end)
		end
	else
		println("not a metric")

	end

    return(f)
end

###################################################################################################
struct alt
	dateTime::String
	altitude::Int64
end

####################################################################################################
struct cal
	dateTime::String
	calories::Float64
end

###################################################################################################
struct demVM
	dateTime::String
	demographicVO2Error::Float64
	demographicVO2Max::Float64
	filteredDemVO2MaxError::Float64
	filteredDemVO2Max::Float64
end

####################################################################################################
struct dist
	dateTime::String
        distance::Int64
end

#####################################################################################################
struct lam
	dateTime::String
	light_active_min::Int64
end

#####################################################################################################
struct mam
        dateTime::String
        moderately_active_min::Int64
end

####################################################################################################
struct vam
	dateTime::String
	very_active_minutes::Int64
end
######################################################################################################
struct rest_hr
	dateTime::String
	resting_hr::Float64
	error::Float64
end
###################################################################################################
struct run_vo2
	dateTime::String
	exerciseId::Int64
	runVO2MaxError::Float64
	filteredRunVO2MaxError::Float64
	filteredRunVO2Max::Float64
	runVO2Max::Float64
end

######################################################################################################

struct sed_min
	dateTime::String
	sedentary_minutes::Int64
end

######################################################################################################

struct step
	dateTime::String
	steps::Int64
end
######################################################################################################
struct swim
	dateTime::String
	strokeCount::Int64
	lapDurationSec::Int64
	swimAlgorithmType::String
	swimStrokeType::String
end

#######################################################################################################
struct tihrz
	dateTime::String
	IN_DEFAULT_ZONE_2::Float64
	BELOW_DEFAULT_ZONE_1::Float64
	IN_DEFAULT_ZONE_1::Float64
	IN_DEFAULT_ZONE_3::Float64
end

########################################################################################################

struct exercs
	dateTime::String
	calories::Int64
	logType::String
	steps::Union{Int64,Missing}
	manualValueSpecified_calories::Bool
	manualValueSpecified_distance::Bool
	manualValueSpecified_steps::Bool
	distanceUnit::Union{String,Missing}
	hasActiveZoneMinutes::Bool
	activeDuration::Int64
	duration::Int64
	vo2Max::Union{Float64,Missing}
	elevationGain::Union{Float64,Missing}
	distance::Union{Float64,Missing}
	activityName::String
	pace::Union{Float64,Missing}
	originalDuration::Int64
	shouldFetchDetails::Bool
	heartRateZones_out_of_range::Union{Int64,Missing}
	heartRateZones_fat_burn::Union{Int64,Missing}
	heartRateZones_cardio::Union{Int64,Missing}
	heartRateZones_peak::Union{Int64,Missing}
	logId::Int64
	activityTypeId::Int64
	averageHeartRate::Union{Int64,Missing}
	hasGps::Bool
	activityLevel_sedentary::Int64
	activityLevel_lightly::Int64
	activityLevel_fairly::Int64
	activityLevel_very::Int64
	speed::Union{Float64,Missing}
	lastModified::String
end

#########################################################################################################

struct sle
	dateTime::String
	endTime::String
	minutesAsleep::Int64
	minutesToFallAsleep::Int64
	minutesAwake::Int64
	duration::Int64
	timeInBed::Int64
	mainSleep::Bool
	infoCode::Int64
	minutesAfterWakeup::Int64
	efficiency::Int64
	type::String
	awake::Int64
	asleep::Union{Int64,Missing}
	restless::Union{Int64,Missing}
	light::Union{Int64,Missing}
	deep::Union{Int64,Missing}
	rem::Union{Int64,Missing}
end
###################################################################################

struct hr
        dateTime::String
        bpm::Int64
	confidence::Int64
end

###################################################################################

ext: jl

max_stars_repo_name: atott/Fitbit-Julia

max_issues_repo_name: atott/Fitbit-Julia

max_forks_repo_name: atott/Fitbit-Julia


text:
function LinearOperator(in,operators,parameters;adj=true)
	if adj
		d = copy(in)
		m = [];
		for j = 1 : 1 : length(operators)
			op = operators[j]
			m = op(d,true;parameters[j]...)
			d = copy(m)
		end
		return m
	else
		m = copy(in)
		d = [];
		for j = length(operators) : -1 : 1
			op = operators[j]
			d = op(m,false;parameters[j]...)
			m = copy(d)
		end
		return d
	end

end

function LinearOperator(m::AbstractString,d::AbstractString,operators,parameters;adj=true)

    rand_string = string(round(Int,rand()*100000))
    tmp_m = join(["tmp_CG_m_",rand_string])
    tmp_d = join(["tmp_CG_d_",rand_string])
    if adj
        SeisCopy(d,tmp_d)
        for j = 1 : 1 : length(operators)
            op = operators[j]
            op(tmp_m,tmp_d,true;parameters[j]...)
            SeisCopy(tmp_m,tmp_d)
        end
        SeisCopy(tmp_m,m)
        SeisRemove(tmp_m)
        SeisRemove(tmp_d)
    else
        SeisCopy(m,tmp_m)
        for j = length(operators) : -1 : 1
            op = operators[j]
            op(tmp_m,tmp_d,false;parameters[j]...)
            SeisCopy(tmp_d,tmp_m)
        end
        SeisCopy(tmp_d,d)
        SeisRemove(tmp_m)
        SeisRemove(tmp_d)
    end

end

function LinearOperator(m::Array{AbstractString,1},d::Array{AbstractString,1},operators,parameters;adj=true)

    rand_string = string(round(Int,rand()*100000))
    tmp_m = [join(["tmp_CG_m1_",rand_string]);join(["tmp_CG_m2_",rand_string])]
    tmp_d = [join(["tmp_CG_d1_",rand_string]);join(["tmp_CG_d2_",rand_string])]
    if adj
        SeisCopy(d,tmp_d)
        for j = 1 : 1 : length(operators)
            op = operators[j]
            if (length(methods(op,(Array,Array,Dict{Any,Any}))) > 0)
                op(tmp_m,tmp_d,true;parameters[j]...)
                SeisCopy(tmp_m,tmp_d)
            else
                for k = 1 : length(m)
                    op(tmp_m[k],tmp_d[k],true;parameters[j]...)
                end
                SeisCopy(tmp_m,tmp_d)
            end
        end
        SeisCopy(tmp_m,m)
        SeisRemove(tmp_m)
        SeisRemove(tmp_d)
    else
        SeisCopy(m,tmp_m)
        for j = length(operators) : -1 : 1
            op = operators[j]
            if (length(methods(op,(Array,Array,Dict{Any,Any}))) > 0)
                op(tmp_m,tmp_d,false;parameters[j]...)
                SeisCopy(tmp_d,tmp_m)
            else
                for k = 1 : length(m)
                    op(tmp_m[k],tmp_d[k],false;parameters[j]...)
                end
                SeisCopy(tmp_d,tmp_m)
            end
        end
        SeisCopy(tmp_d,d)
        SeisRemove(tmp_m)
        SeisRemove(tmp_d)
    end

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
export prufer_code, is_tree, prufer_restore

"""
    is_tree(G)
Determines if the `SimpleGraph` is a tree.
"""
function is_tree(G::SimpleGraph)
    return is_connected(G) && (NE(G) == NV(G) - 1)
end

"""
    lowest_leaf(G)
Return the leaf of the `SimpleGraph` with the smallest label.
"""
function lowest_leaf(G::SimpleGraph)
    leaves = [v for v ∈ G.V if deg(G, v) == 1]
    return findmin(leaves)[1]
end

"""
    lowest_leaf_neighbor(G::SimpleGraph)
Return the unique neighbor of `lowest_leaf(G)`.
"""
function lowest_leaf_neighbor(G::SimpleGraph)
    v = lowest_leaf(G)
    return G[v][1]
end

"""
    prufer_code(G)
Return the Prufer code of the `SimpleGraph` which must be a tree.
The vertices must be `<`-comparable and preferrably are the integers
`{1,2,...,n}` (otherwise we cannot decode the sequence generated).
"""
function prufer_code(G::SimpleGraph{T})::Vector{T} where {T}
    if !is_tree(G)
        error("Graph must be a tree")
    end

    n = NV(G)

    if n < 2
        error("Graph must have at least two vertices")
    end

    if G.V != Set(1:n)
        @warn "The vertex set of this tree is not of the form {1,2,...,n}"
    end

    if n <= 2
        return T[]
    end
    GG = deepcopy(G)
    return prufer_work(GG)

end

function prufer_work(G::SimpleGraph{T})::Vector{T} where {T}
    if NV(G) <= 2
        return T[]
    end
    v = lowest_leaf(G)
    w = G[v][1]   # neighbor of lowest leaf

    delete!(G, v)

    code = prufer_work(G)
    prepend!(code, [w])
    return code
end

"""
    prufer_restore(code::Vector{Int})
Create a tree from its Prufer code.
"""
function prufer_restore(code::Vector{Int})
    n = length(code)
    G = IntGraph(n+2)

    # find the degree sequence of this graph
    ds = ones(Int,n+2)
    for v ∈ code
        ds[v] += 1
    end

    for v ∈ code
        # find first leaf in
        leaves = [w for w ∈ G.V if ds[w] == 1]
        w = findmin(leaves)[1]
        add!(G,v,w)
        ds[v] -= 1
        ds[w] -= 1
    end

    leaves = [w for w ∈ G.V if ds[w] == 1]
    add!(G, leaves[1], leaves[2])

    return G
end

code_to_tree = prufer_restore  # for backward compatability


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput03.txt")
function TobogganTrajectory(data; current_pos = 1, right = 3, down = 1)
    c = 0
    lenrow = length(data[1])
    for row in data[1:down:end]
        if current_pos > lenrow
            row = row ^(Int(ceil(current_pos/lenrow)))
            🌲 = row[current_pos]
            if 🌲 == '#'
                c += 1
            end
        else
            🌲 = row[current_pos]
            if 🌲 == '#'
                c += 1
            end
        end
        current_pos += right
    end
    c
end
TobogganTrajectory(data)

movesRight = [1,3,5,7,1]
movesDown = [1,1,1,1,2]
narboles = [TobogganTrajectory(data; right = movesRight[n], down = movesDown[n])
        for n in 1:length(movesRight)]
prod(narboles)


text:
using FFTW
export QFTCircuit

CRk(i::Int, j::Int, k::Int) = control([i, ], j=>shift(2π/(1<<k)))
CRot(n::Int, i::Int) = chain(n, i==j ? put(i=>H) : CRk(j, i, j-i+1) for j = i:n)
QFTCircuit(n::Int) = chain(n, CRot(n, i) for i = 1:n)


text:
using FFTW
using DSP
using ImageTransformations
using Dierckx
using Plots
using Colors
using PyCall
using LinearAlgebra
using CoordinateTransformations, Rotations, StaticArrays
using Dierckx

function rescale(x)
    mx = maximum(x)
    mi = minimum(x)
    return map(i-> (i-mi)/(mx-mi),x)
end

function ident(ϕ)
     t1(x) = x
     t2(x) = x
     return (t1,t2)
end

function nonlinear_transformation(ϕ)
    m1 = (sin(0.002*(ϕ*(257/(2*π))-1)*3/256))/4
    m2 = (sin(0.005*(ϕ*(257/(2*π))-1)*3/256))/4
    f(x) = x-128.5
    g1(x) = m1*x
    g2(x) = m2*x
    h(x) = x+1
    j(x) = x^5
    k(x) = x-1
    l1(x) = x/(5*m1)
    l2(x) = x/(5*m2)
    q(x) = x+128.5
    v1(x) = x |> f |> g1 |> h |> j |> k |> l1 |> q
    v2(x) = x |> f |> g2 |> h |> j |> k |> l2 |> q

    t(x) = SArray{Tuple{2},Float64,1,2}(v1(x[1]),v2(x[2]))
    return t
end

function nonlinear_transformation_inv(ϕ)
    m1 = (sin(0.002*(ϕ*(257/(2*π))-1)*3/256))/4
    m2 = (sin(0.005*(ϕ*(257/(2*π))-1)*3/256))/4
    f(x) = x+128.5
    g1(x) = x/m1
    g2(x) = x/m2
    h(x) = x-1
    j(x) = abs((x+0.0*im)^(1/5))*sign(x)
    k(x) = x+1
    l1(x) = x*(5*m1)
    l2(x) = x*(5*m2)
    q(x) = x-128.5
    v1(x) = x |> q |> l1 |> k |> j |> h |> g1 |> f
    v2(x) = x |> q |> l2 |> k |> j |> h |> g2 |> f

    return (v1,v2)
end

replace_nan(v) = map(x -> isnan(x) ? zero(x) : x, v)

function precompute_kernels(ρ)
    f = fftfreq(length(ρ))
    fourier_filter = 2 * abs.(f)
    return fourier_filter
end

function filtering(data,kernel)
    FT = fft(data)
    g = real.(ifft(kernel.*FT))
    return g
end

function dynamic_backprojection()
end

function Hahn(sinogram, reconstrcution, inverse_motion)
    return reconstruction
end

function backprojection(data_filtered, θ, inv)
    (Γ_inv_x, Γ_inv_y) = inv(θ)
    output_size = length(data_filtered)
    mid_index = (output_size+1) / 2
    X = Γ_inv_x.(repeat(collect(1:output_size), 1, output_size))
    Y = Γ_inv_y.(repeat(collect(1:output_size)', output_size, 1))
    xpr = X .- mid_index
    ypr = Y .- mid_index
    t = ypr * cos(θ) - xpr * sin(θ)
    x = collect(1:output_size) .- mid_index
    spl = Spline1D(x, data_filtered[1,:], k=1, bc="zero")
    backprojected = spl.(t)
    return backprojected
end

function rod_phantom()
    rod = zeros(Float64, (600,600))
    rod[292:308,300:600] .= 1.0
    return rod
end

function disc_phantom(img, a, b, r)
    height, width = size(img)
    #find the possibility which has most surrounding pixels
    disc = zeros(Float64, (height,width))
    for x = 1:width
        for y = 1:height
            disc[y,x] = (x-a)^2 + (y-b)^2 < r^2 ? 1.0 : 0.0
        end
    end
    return disc
end

function FBP(sinogram, θ, ρ, inv)
    N = length(ρ)
    kernel = precompute_kernels(ρ)
    reconstruction = zeros(Float64,N,N)
    for i =1:length(θ)
        g = filtering(sinogram[:,i],kernel)
        reconstruction += backprojection(g', θ[i], inv)
    end
    return reconstruction * (π / (2 * length(θ)))
end

function radon(image,angles)
    theta = rad2deg.(angles)
    trans = pyimport("skimage.transform")
    sino = trans.radon(image, theta, true)
    return sino
end

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)

c_1 = 128.5
c_2 = 128.5
rad = 50

H,W = 256, 256

img = zeros(Float64, (H, W))
img = disc_phantom(img, c_1, c_2, rad)
p1 = plot(Gray.(img), aspect_ratio=:equal)
angles = collect(range(0,2*π, length=201))
pop!(angles)
bins = collect(-127.5:127.5)

sino = zeros(Float64,length(bins), length(angles))

anim = @animate for i = 1:length(angles)
    trans =nonlinear_transformation(angles[i])
    imgw = replace_nan(warp(img, trans, axes(img)))
    plot(Gray.(imgw), aspect_ratio=:equal, framestyle=:none)
    sino[:,i] = radon(imgw,[angles[i]])
end

gif(anim, "deformation.gif", fps = 20)

display(plot(Gray.(rescale(sino))))

reconstruction = FBP(sino, angles, bins, nonlinear_transformation_inv)
display(plot(Gray.(rescale(reconstruction))))

reconstruction_no_comp = FBP(sino, angles, bins, ident)
display(plot(Gray.(rescale(reconstruction_no_comp))))

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
insertsort!(A) = reduce(1:length(A)) do _, k
    while k>1 && A[k-1] > A[k]
        A[k-1], A[k] = A[k], A[k-1]
        k-=1
    end
    return A
end


text:
export NoHandler, DataMissingKey

struct NoHandler <: Exception
    stream::Stream
end

struct MissingCookies{T} <: Exception
    t::Type{T}
    k::Vector{String}
end

struct MissingHeaders{T} <: Exception
    t::Type{T}
    k::Array{String}
end

struct DataMissingKey{T} <: Exception
    t::Type{T}
    struct_keys::Array{Symbol}
    data_keys::Array{Symbol}
end

function Base.show(io::IO, e::DataMissingKey)
    println(io, "DataMissingKey:")
    println(io, "  Expected - $(e.struct_keys)")
    print(io, "  Given    - $(e.data_keys)")
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
#   Unit "fftw.jl" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements a convenient interface to FFTW.jl and AbstractFFTs.jl

# ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤

# flags to be passed as `planflag` parameter in several functions
# see www.fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags for details
const plan_estimate        = UInt32(64)      # very small search =(1 << 6)
const plan_measure         = UInt32(0)       # small search
const plan_patient         = UInt32(32)      # large search =(1 << 5)
const plan_exhaustive      = UInt32(8)       # very large search =(1 << 3)
const plan_conserve_memory = UInt32(4)       # = (1 << 2)
const plan_wisdom_only     = UInt32(2097152) # = (1 << 21)
const plan_unaligned       = UInt32(2)       # = (1 << 1)
const plan_preserve_input  = UInt32(16)      # do not copy the input data = (1 << 4)

"""
FFTW plans to be used by all functions in **Fourier Analysis**
are incapsulated in this structure. This is the only *operator object*
created by this package, all others being *data objects*.

Spectral and cross-spectral computations (frequency domain objects)
need only a forward plan (real FFT),
while objects based on the analytic signal (time-frequency domain objects)
need both a forward (`.p `) and a backward (`.ip `) plan (complex iFFT).

Argument `flags` must be one of the constants here above.
Basic usage of the flags involves a trade-off between the time needed to compute
the planner and the efficacy for computing the FFTs. The following flags
are sorted in ascending order of time needed to be computed and efficacy:
`plan_estimate`, `plan_measure`, `plan_patient`,
`plan_exhaustive`.
By default *FourierAnalysis* adopts `plan_estimate`,
which allows the quickest search, but the least optimal choice.
Other flags may be worth if several FFT computations are to be done
with the same setting.
See [here](http://www.fftw.org/fftw3_doc/Planner-Flags.html) for details.

Argument `timelimit` is the maximum time (in seconds) allowed to FFTW
for optimizing the plans. Setting this to `-1.0` amounts to imposing
no time limits.

FFTW plans are computed for a given window length `wl` and data type `type`.

`Planners` objects may be passed as an argument to constructors of
[FDobjects](@ref) and [TFobjects](@ref) by using one of the `Planner` constuctor
here below.

**Constructors**

```julia
Planner(flags :: UInt32,
    timelimit :: Union{Int, Float64},
           wl :: Int,
         type :: Type,
           bw :: Bool = false)
```

Use this to create a Planner object, passing as argument a FFTW `flags` constant
and `timelimit`, the window length `wl` and the type of the input data `type`,
which must be real, e.g., Float64.

If `bw` is false (default), a dummy backward plan is created,
otherwise a backward plan is created for the complex data type corresponding
to `type` (e.g., if `type` is Float64, it will created for ComplexF64 data.)

For example, suppose 𝐗 is a vector of many matrices of multivariate time-series
sampled at 128 samples per second and that we want to compute the spectra
for all of them using a 256-point FFT. We first create a plan by

```julia
p=Planner(plan_exhaustive, 10.0, 256, eltype(𝐗[1]))
```

Then we invoke the [Spectra](@ref) function passing the plan as argument:

```julia
𝐒=spectra(𝐗, sr, t; planner=p)
```

A shorter construction syntax is available when only the
forward plan is needed and the type of the data is Float64:

```julia
Planner(flags :: UInt32,
    timelimit :: Union{Int, Float64},
           wl :: Int)
```

For example, the line above could have been written more shortly as

```julia
p=Planner(plan_exhaustive, 10.0, 256)
```

"""
struct Planner
    flags     ::  UInt32            # see
    timelimit ::  Float64           # =for no time limit set to -1.0
    wl        ::  Int               # window length
    type      ::  Type              # element type of forward planner (e.g., Float64)
    p         ::  FFTW.rFFTWPlan    # forward plan (FFT)
    ip        ::  FFTW.cFFTWPlan    # backward plan (inverse FFT)
end

Planner(flags::UInt32, timelimit::Union{Int, Float64}, wl::Int, type::Type, bw::Bool=false) =
    bw ? Planner(flags, Float64(timelimit), wl, type, plan_rfft(zeros(type, wl);
                 flags=flags, timelimit=timelimit), plan_bfft(zeros(eltype(Complex(type(0))), wl); flags=flags, timelimit=timelimit)) :
         Planner(flags, Float64(timelimit), wl, type, plan_rfft(zeros(type, wl);
                 flags=flags, timelimit=timelimit), plan_bfft(zeros(ComplexF64, 2)))

# create a planner with only a forward plan defaulting to the Float64 type
Planner(flags::UInt32, timelimit::Union{Int, Float64}, wl::Int) =
    Planner(flags, Float64(timelimit), wl, Float64, plan_rfft(zeros(Float64, wl), flags=flags, timelimit=timelimit), plan_bfft(zeros(ComplexF64, 2)))

# example of full constructor:
# plan=Planner( plan_exhaustive, 10, 256, Float64,
#               plan_rfft(zeros(Float64, 256), plan_exhaustive, timelimit=10),
#               plan_bfft(zeros(ComplexF64, 256), plan_exhaustive, timelimit=10))

# fast constructor writing only the forward planner passed as argument.
# all other fields are dummy. This is used internally by spectra and
# cross-spectra functions
_Planner(p::FFTW.rFFTWPlan)=Planner(plan_estimate, -1.0, 0, Float64, p, plan_bfft(zeros(ComplexF64, 2)))

# dummy FFTW planner used for default settings. Only for internal use
getplanner=_Planner(plan_rfft(zeros(2)))

_getflags(flags::UInt32) =
 if     flags==UInt32(64) return "plan_estimate"
 elseif flags==UInt32(0)  return "plan_measure"
 elseif flags==UInt32(32) return "plan_patient"
 elseif flags==UInt32(8)  return "plan_exhaustive"
 elseif flags==UInt32(4)  return "plan_conserve_memory"
 elseif flags==UInt32(2097152) return "plan_wisdom_only"
 elseif flags==UInt32(2)  return "plan_unaligned"
 elseif flags==UInt32(16) return "plan_preserve_input"
 elseif @error 📌*", the FFTW flags does not exist"
 end

 # ++++++++++++++++++++  Show override  +++++++++++++++++++ # (REPL output)
function Base.show(io::IO, ::MIME{Symbol("text/plain")}, p::Planner)
    println(io, titleFont, "∿ Planner type")
    #println(io, "□  □    □      □        □           □", defaultFont)
    println(io, separatorFont, "⭒  ⭒    ⭒      ⭒        ⭒           ⭒", defaultFont)
    println(io, "FFTW flags        (.flags): "*_getflags(p.flags)*" ($(p.flags))")
    p.timelimit==-1.0 ? println(io, "FFTW timelimit(.timelimit): no limit") :
                        println(io, "FFTW timelimit(.timelimit): $(p.timelimit)")
    println(io, "window length        (.wl): $(p.wl)")
    println(io, "element type       (.type): $(p.type)")
    println(io, "forward plan          (.p): ")
    println(io, greyFont, "$(p.p)", defaultFont)
    println(io, "backward plan        (.ip): ")
    println(io, greyFont, "$(p.ip)", defaultFont)
end
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ #
# ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤


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

using CrimsonSkyline
using Distributions: Beta, Geometric
using Logging
using StatsBase: mean, std
using Random: seed!
using Plots

seed!(2021)

num_constraint(x, num, bound) = (x > num - bound) && (x < num + bound)

function constrained(t::Trace, evidence::Vector{Int64}, bound::Int64)
    prob = sample(t, "prob", Beta(2.0, 2.0))
    num = sample(t, "num", Geometric(prob))
    plate(t, observe, "constraint", SoftDelta(x -> num_constraint(x, num, bound), -2.0), evidence)
end

function unconstrained(t::Trace, evidence::Vector{Int64})
    prob = sample(t, "prob", Beta(2.0, 2.0))
    plate(t, observe, "num", Geometric(prob), evidence)
end

function main()
    data = [1, 5, 2, 9, 3, 2, 1, 3, 2, 1, 1, 4, 3, 2]
    bound = 3
    @time constrained_results = mh(constrained; params = (data, bound))
    @time unconstrained_results = mh(unconstrained; params = (data,))
    aic_constrained = aic(constrained_results)
    aic_unconstrained = aic(unconstrained_results)
    @info "AIC of constrained model = $aic_constrained"
    @info "AIC of unconstrained model = $aic_unconstrained"
    hpdi_con = hpdi(constrained_results, 0.5, ["prob"])
    hpdi_uncon = hpdi(unconstrained_results, 0.5, ["prob"])
    @info "50% HPDI of prob under constrained model: $hpdi_con"
    @info "50% HPDI of prob under unconstrained model: $hpdi_uncon"
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
function update_hangartner_temp_vars!(u::Matrix{Int}, X::Matrix{Int}, t::Int)
    d = size(X, 2)

    for j in 1:d
        u[X[t, j], j] += 1
    end
end

function hangartner_inner(Y::AbstractMatrix, m::Int)
    ## setup temp vars
    n, d = size(Y)

    # Count for each category in each chain
    u = zeros(Int, m, d)
    v = zeros(Int, m, d)

    for t in 1:n
        # fill out temp vars
        update_hangartner_temp_vars!(u, Y, t)
    end
    phia, chi_stat, m_tot = weiss_sub(u, v, n)

    return (n * sum(chi_stat), m_tot)
end

"""
    weiss(X::AbstractMatrix)

Assess the convergence of the MCMC chains with the Weiss procedure.

It computes ``\\frac{X^2}{c}`` and evaluates a p-value from the ``\\chi^2`` distribution with ``(|R| − 1)(s − 1)`` degrees of freedom.
"""
function weiss(X::AbstractMatrix)
    ## number of iterations, number of chains
    n, d = size(X)

    ## mapping of values to integers
    v_dict = Dict{eltype(X),Int}()

    ## max unique categories
    mc = map(c -> length(unique(X[:, c])), 1:d)
    m = length(unique(X))

    ## counter of number of unique values in each chain
    r0 = 0

    ## Count for each category in each chain
    u = zeros(Int, m, d)

    ## Number of times a category did not transition in each chain
    v = zeros(Int, m, d)

    for t in 1:n
        for c in 1:d
            if !(X[t, c] in keys(v_dict))
                r0 += 1
                v_dict[X[t, c]] = r0
            end
            idx1 = v_dict[X[t, c]]
            u[idx1, c] += 1

            if t > 1
                if X[t - 1, c] == X[t, c]
                    v[idx1, c] += 1
                end
            end
        end
    end
    phia, chi_stat, m_tot = weiss_sub(u, v, n)
    ca = (1 + phia) / (1 - phia)
    stat = (n / ca) * sum(chi_stat)
    pval = NaN
    if ((m_tot - 1) * (d - 1)) >= 1
        pval = 1 - Distributions.cdf(Distributions.Chisq((m_tot - 1) * (d - 1)), stat)
    end

    return (stat, m_tot, pval, ca)
end

function weiss_sub(u::Matrix{Int}, v::Matrix{Int}, t::Int)
    m, d = size(u)
    nt = 0.0
    dt = 0.0
    m_tot = 0

    mp = zeros(Float64, m, d)
    ma = zeros(Float64, m)
    phia = 0.0
    ca = 0.0
    df = 0.0

    chi_stat = zeros(Float64, d)

    for j in 1:m
        p1 = 0.0
        p2 = 0.0
        for l in 1:d
            #aggregate
            p1 += v[j, l] / (d * (t - 1))
            p2 += u[j, l] / (d * t)

            #per chain
            mp[j, l] = u[j, l] / t
            ma[j] += u[j, l] / (d * t)
        end
        nt += p1
        dt += p2^2

        if ma[j] > 0
            m_tot += 1
            for l in 1:d
                chi_stat[l] += (mp[j, l] - ma[j])^2 / ma[j]
            end
        end
    end
    phia = 1.0 + (1.0 / t) - ((1 - nt) / (1 - dt))
    phia = min(max(phia, 0.0), 1.0 - eps())
    return (phia, chi_stat, m_tot)
end

function update_billingsley_temp_vars!(f::Array{Int,3}, X::Matrix{Int}, t::Int)
    d = size(X, 2)
    for j in 1:d
        if t > 1
            f[X[t - 1, j], X[t, j], j] += 1
        end
    end
end

function billingsley_sub(f::Array{Int,3})
    df = 0.0
    stat = 0.0

    m, d = size(f)[2:3]

    # marginal transitions, i.e.
    # number of transitions from each category
    mf = mapslices(sum, f; dims=[2])

    # For each category, number of chains for which
    # that category was present
    A = vec(mapslices((x) -> sum(x .> 0), mf; dims=[3]))

    # For each category, number of categories it
    # transitioned to
    B = vec(mapslices((x) -> sum(x .> 0), mapslices(sum, f; dims=[3]); dims=[2]))

    # transition probabilities in each chain
    P = f ./ mf

    # transition probabilities
    mP = (mapslices(sum, f; dims=[3]) ./ mapslices(sum, mf; dims=[3]))
    mP = reshape(mP, size(mP)[1:2])

    idx = findall((A .* B) .> 0)
    for j in idx
        #df for billingsley
        df += (A[j] - 1) * (B[j] - 1)

        #billingsley
        for k in idx
            if (mP[j, k] > 0.0)
                for l in 1:d
                    if mf[j, 1, l] > 0 && isfinite(P[j, k, l])
                        stat += mf[j, 1, l] * (P[j, k, l] - mP[j, k])^2 / mP[j, k]
                    end
                end
            end
        end
    end
    return (stat, df, mP)
end

function bd_inner(Y::AbstractMatrix, m::Int)
    num_iters, num_chains = size(Y)
    # Transition matrix for each chain
    f = zeros(Int, m, m, num_chains)

    for t in 1:num_iters
        # fill out temp vars
        update_billingsley_temp_vars!(f, Y, t)
    end
    return billingsley_sub(f)
end

function simulate_NDARMA(
    N::Int, p::Int, q::Int, prob::Vector{Float64}, phi::Vector{Float64}
)
    sampler1 = Distributions.sampler(Distributions.Multinomial(1, phi))
    sampler2 = Distributions.sampler(Distributions.Categorical(prob))

    alphabeta = Vector{Int}(undef, length(phi))
    alphabeta_p = @view(alphabeta[1:p])
    alphabeta_q = @view(alphabeta[(end - q):end])
    eps = Vector{Int}(undef, q + 1)
    X = Vector{Int}(undef, N)
    Random.rand!(sampler2, @view(X[1:p]))
    for t in (p + 1):N
        Random.rand!(sampler1, alphabeta)
        Random.rand!(sampler2, eps)
        X[t] =
            LinearAlgebra.dot(alphabeta_p, @view(X[(t - p):(t - 1)])) +
            LinearAlgebra.dot(alphabeta_q, eps)
    end

    return X
end

function simulate_MC(N::Int, P::Matrix{Float64})
    X = zeros(Int, N)
    n, m = size(P)
    X[1] = StatsBase.sample(1:n)
    for i in 2:N
        X[i] = StatsBase.wsample(1:n, vec(P[X[i - 1], :]))
    end
    return X
end

function diag_all(
    X::AbstractMatrix, method::Symbol, nsim::Int, start_iter::Int, step_size::Int
)

    ## number of iterations, number of chains
    n, d = size(X)

    ## mapping of values to integers
    v_dict = Dict{eltype(X),Int}()

    ## max unique categories
    mc = map(c -> length(unique(X[:, c])), 1:d)
    m = length(unique(X))

    ## counter of number of unique values in each chain
    r0 = 0

    ## Count for each category in each chain
    u = zeros(Int, m, d)

    ## Number of times a category did not transition in each chain
    v = zeros(Int, m, d)

    ## transition matrix for each chain
    f = zeros(Int, m, m, d)

    length_result = length(start_iter:step_size:n)
    result = (
        stat=Vector{Float64}(undef, length_result),
        df=Vector{Float64}(undef, length_result),
        pvalue=Vector{Float64}(undef, length_result),
    )
    result_iter = 1
    for t in 1:n
        for c in 1:d
            if !(X[t, c] in keys(v_dict))
                r0 += 1
                v_dict[X[t, c]] = r0
            end
            idx1 = v_dict[X[t, c]]
            u[idx1, c] += 1

            if t > 1
                idx2 = v_dict[X[t - 1, c]]
                f[idx1, idx2, c] += 1

                if X[t - 1, c] == X[t, c]
                    v[idx1, c] += 1
                end
            end
        end

        if ((t >= start_iter) && (((t - start_iter) % step_size) == 0))
            phia, chi_stat, m_tot = weiss_sub(u, v, t)
            hot_stat, df, mP = billingsley_sub(f)

            phat = mapslices(sum, u; dims=[2])[:, 1] / sum(mapslices(sum, u; dims=[2]))
            ca = (1 + phia) / (1 - phia)
            stat = NaN
            pval = NaN
            df0 = NaN

            if method == :hangartner
                stat = t * sum(chi_stat)
                df0 = (m - 1) * (d - 1)
                if m > 1 && !isnan(stat)
                    pval =
                        1 - Distributions.cdf(Distributions.Chisq((m - 1) * (d - 1)), stat)
                end
            elseif method == :weiss
                stat = (t / ca) * sum(chi_stat)
                df0 = (m - 1) * (d - 1)
                pval = NaN
                if m > 1 && !isnan(stat)
                    pval =
                        1 - Distributions.cdf(Distributions.Chisq((m - 1) * (d - 1)), stat)
                end
            elseif method == :DARBOOT
                stat = t * sum(chi_stat)
                bstats = zeros(Float64, nsim)
                for b in 1:nsim
                    Y = reduce(
                        hcat,
                        [simulate_NDARMA(t, 1, 0, phat, [phia, 1 - phia]) for j in 1:d],
                    )
                    s = hangartner_inner(Y, m)[1]
                    bstats[b] = s
                end
                non_nan_bstats = filter(!isnan, bstats)
                df0 = Statistics.mean(non_nan_bstats)
                pval = Statistics.mean(stat <= x for x in non_nan_bstats)
            elseif method == :MCBOOT
                bstats = zeros(Float64, nsim)
                for b in 1:nsim
                    Y = reduce(hcat, [simulate_MC(t, mP) for j in 1:d])
                    s = hangartner_inner(Y, m)[1]
                    bstats[b] = s
                end
                non_nan_bstats = filter(!isnan, bstats)
                df0 = Statistics.mean(non_nan_bstats)
                pval = Statistics.mean(stat <= x for x in non_nan_bstats)
            elseif method == :billingsley
                stat = hot_stat
                df0 = df
                if df > 0 && !isnan(hot_stat)
                    pval = 1 - Distributions.cdf(Distributions.Chisq(df), hot_stat)
                end
            elseif method == :billingsleyBOOT
                stat = hot_stat
                bstats = zeros(Float64, nsim)
                for b in 1:nsim
                    Y = reduce(hcat, [simulate_MC(t, mP) for j in 1:d])
                    (s, sd) = bd_inner(Y, m)[1:2]
                    bstats[b] = s / sd
                end
                non_nan_bstats = filter(!isnan, bstats)
                df0 = Statistics.mean(non_nan_bstats)
                statodf = stat / df
                pval = Statistics.mean(statodf <= x for x in non_nan_bstats)
            else
                error("Unexpected")
            end
            result.stat[result_iter] = stat
            result.df[result_iter] = df0
            result.pvalue[result_iter] = pval
            result_iter += 1
        end
    end
    return result
end

function discretediag_sub(
    c::AbstractArray{<:Real,3},
    frac::Real,
    method::Symbol,
    nsim::Int,
    start_iter::Int,
    step_size::Int,
)
    num_iters, num_vars, num_chains = size(c)

    ## Between-chain diagnostic
    length_results = length(start_iter:step_size:num_iters)
    plot_vals_stat = Matrix{Float64}(undef, length_results, num_vars)
    plot_vals_pval = Matrix{Float64}(undef, length_results, num_vars)
    between_chain_vals = (
        stat=Vector{Float64}(undef, num_vars),
        df=Vector{Float64}(undef, num_vars),
        pvalue=Vector{Float64}(undef, num_vars),
    )
    for j in 1:num_vars
        X = convert(AbstractMatrix{Int}, c[:, j, :])
        result = diag_all(X, method, nsim, start_iter, step_size)

        plot_vals_stat[:, j] .= result.stat ./ result.df
        plot_vals_pval[:, j] .= result.pvalue

        between_chain_vals.stat[j] = result.stat[end]
        between_chain_vals.df[j] = result.df[end]
        between_chain_vals.pvalue[j] = result.pvalue[end]
    end

    ## Within-chain diagnostic
    within_chain_vals = (
        stat=Matrix{Float64}(undef, num_vars, num_chains),
        df=Matrix{Float64}(undef, num_vars, num_chains),
        pvalue=Matrix{Float64}(undef, num_vars, num_chains),
    )
    for k in 1:num_chains
        for j in 1:num_vars
            x = convert(AbstractVector{Int}, c[:, j, k])

            idx1 = 1:round(Int, frac * num_iters)
            idx2 = round(Int, num_iters - frac * num_iters + 1):num_iters
            x1 = x[idx1]
            x2 = x[idx2]
            n_min = min(length(x1), length(x2))
            Y = [x1[1:n_min] x2[(end - n_min + 1):end]]

            result = diag_all(Y, method, nsim, n_min, step_size)
            within_chain_vals.stat[j, k] = result.stat[end]
            within_chain_vals.df[j, k] = result.df[end]
            within_chain_vals.pvalue[j, k] = result.pvalue[end]
        end
    end

    return between_chain_vals, within_chain_vals, plot_vals_stat, plot_vals_pval
end

"""
    discretediag(chains::AbstractArray{<:Real,3}; frac=0.3, method=:weiss, nsim=1_000)

Compute discrete diagnostic where `method` can be one of `:weiss`, `:hangartner`,
`:DARBOOT`, `:MCBOOT`, `:billinsgley`, and `:billingsleyBOOT`.
"""
function discretediag(
    chains::AbstractArray{<:Real,3}; frac::Real=0.3, method::Symbol=:weiss, nsim::Int=1000
)
    valid_methods = (:weiss, :hangartner, :DARBOOT, :MCBOOT, :billingsley, :billingsleyBOOT)
    method in valid_methods || throw(
        ArgumentError("`method` must be one of :" * join(valid_methods, ", :", " and :")),
    )
    0 < frac < 1 || throw(ArgumentError("`frac` must be in (0,1)"))

    num_iters = size(chains, 1)
    between_chain_vals, within_chain_vals, _, _ = discretediag_sub(
        chains, frac, method, nsim, num_iters, num_iters
    )

    return between_chain_vals, within_chain_vals
end


text:
# const ID = Int
const IDT = UInt64
const ID = Vector{IDT}
const uidcounter = Counter(IDT)
toid(id::Integer) = ID([id])

"Unique id, (scope 1), (scope, 2), ..."
uid(; scope = scope()) = append(scope, increment!(uidcounter))
@spec :nocheck (x = [uid() for i = 1:Inf]; unique(x) == x)

"Construct globally unique id for indices for ω"
macro uid()
  uid()
end

# Scope
const SCOPE = Box(base(ID))
scope()::ID = SCOPE.val

"""
Introduce a scope where `id` is prepended to any id created in this scope
```
withscope(21) do
  x =~ uniform(0, 1)
  y =~ uniform(0, 1)
end
```
"""
function withscope(f, id::ID)
  oldscope = scope()
  scope!(id)
  ret = f()
  scope!(oldscope)
  ret
end

withscope(f, ID) = withscope(f, base(ID, hash(ID)))

"Change the global scope to `id`"
function scope!(id::ID)::ID
  reset!(uidcounter)
  SCOPE.val = id
end

scope!(id) = scope!(base(ID, hash(id)))


text:
mutable struct TransportDirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Box, Ball}}
    agent_pos::CartesianIndex{2}
    agent_dir::AbstractDirection
    reward::T
    rng::R
    terminal_reward::T
    ball_pos::CartesianIndex{2}
    box_pos::CartesianIndex{2}
    has_ball::Bool
    done::Bool
end

@generate_getters(TransportDirected)
@generate_setters(TransportDirected)

mutable struct TransportUndirected{T, R} <: AbstractGridWorld
    world::GridWorldBase{Tuple{Agent, Wall, Box, Ball}}
    agent_pos::CartesianIndex{2}
    reward::T
    rng::R
    terminal_reward::T
    ball_pos::CartesianIndex{2}
    box_pos::CartesianIndex{2}
    has_ball::Bool
    done::Bool
end

@generate_getters(TransportUndirected)
@generate_setters(TransportUndirected)

#####
# Directed
#####

function TransportDirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, BOX, BALL)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    box_pos = CartesianIndex(height - 1, width - 1)
    world[BOX, box_pos] = true

    ball_pos = CartesianIndex(2, 3)
    world[BALL, ball_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    agent_dir = RIGHT
    reward = zero(T)
    terminal_reward = one(T)
    has_ball = false
    done = false

    env = TransportDirected(world, agent_pos, agent_dir, reward, rng, terminal_reward, ball_pos, box_pos, has_ball, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::TransportDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const TRANSPORT_DIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::TransportDirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_agent_dir(env), get_half_size(env), TRANSPORT_DIRECTED_LAYERS)

RLBase.state_space(env::TransportDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::TransportDirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = (copy(get_grid(env)), get_agent_dir(env))

RLBase.action_space(env::TransportDirected, ::RLBase.DefaultPlayer) = (DIRECTED_NAVIGATION_ACTIONS..., PICK_UP, DROP)
RLBase.reward(env::TransportDirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::TransportDirected) = get_done(env)

function RLBase.reset!(env::TransportDirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[BOX, get_box_pos(env)] = false
    world[BALL, get_ball_pos(env)] = false

    new_box_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_box_pos!(env, new_box_pos)
    world[BOX, new_box_pos] = true

    new_ball_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_ball_pos!(env, new_ball_pos)
    world[BALL, new_ball_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_agent_dir!(env, rand(rng, DIRECTIONS))

    set_reward!(env, zero(T))
    set_has_ball!(env, false)
    set_done!(env, false)

    return nothing
end

function (env::TransportDirected{T})(action::AbstractTurnAction) where {T}
    new_dir = turn(action, get_agent_dir(env))
    set_agent_dir!(env, new_dir)
    world = get_world(env)

    if world[BALL, get_box_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::TransportDirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, get_agent_dir(env), agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[BALL, get_box_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Undirected
#####

function TransportUndirected(; T = Float32, height = 8, width = 8, rng = Random.GLOBAL_RNG)
    objects = (AGENT, WALL, BOX, BALL)
    world = GridWorldBase(objects, height, width)

    room = Room(CartesianIndex(1, 1), height, width)
    place_room!(world, room)

    box_pos = CartesianIndex(height - 1, width - 1)
    world[BOX, box_pos] = true

    ball_pos = CartesianIndex(2, 3)
    world[BALL, ball_pos] = true

    agent_pos = CartesianIndex(2, 2)
    world[AGENT, agent_pos] = true

    reward = zero(T)
    terminal_reward = one(T)
    has_ball = false
    done = false

    env = TransportUndirected(world, agent_pos, reward, rng, terminal_reward, ball_pos, box_pos, has_ball, done)

    RLBase.reset!(env)

    return env
end

RLBase.state_space(env::TransportUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = nothing
const TRANSPORT_UNDIRECTED_LAYERS = SA.SVector(2, 3, 4)
RLBase.state(env::TransportUndirected, ::RLBase.Observation, ::RLBase.DefaultPlayer) = get_grid(get_world(env), get_agent_pos(env), get_half_size(env), TRANSPORT_UNDIRECTED_LAYERS)

RLBase.state_space(env::TransportUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = nothing
RLBase.state(env::TransportUndirected, ::RLBase.InternalState, ::RLBase.DefaultPlayer) = copy(get_grid(env))

RLBase.action_space(env::TransportUndirected, player::RLBase.DefaultPlayer) = (UNDIRECTED_NAVIGATION_ACTIONS..., PICK_UP, DROP)
RLBase.reward(env::TransportUndirected, ::RLBase.DefaultPlayer) = get_reward(env)
RLBase.is_terminated(env::TransportUndirected) = get_done(env)

function RLBase.reset!(env::TransportUndirected{T}) where {T}
    world = get_world(env)
    rng = get_rng(env)

    world[AGENT, get_agent_pos(env)] = false
    world[BOX, get_box_pos(env)] = false
    world[BALL, get_ball_pos(env)] = false

    new_box_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_box_pos!(env, new_box_pos)
    world[BOX, new_box_pos] = true

    new_ball_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_ball_pos!(env, new_ball_pos)
    world[BALL, new_ball_pos] = true

    new_agent_pos = rand(rng, pos -> !any(@view world[:, pos]), env)
    set_agent_pos!(env, new_agent_pos)
    world[AGENT, new_agent_pos] = true

    set_reward!(env, zero(T))
    set_has_ball!(env, false)
    set_done!(env, false)

    return nothing
end

function (env::TransportUndirected{T})(action::AbstractMoveAction) where {T}
    world = get_world(env)

    agent_pos = get_agent_pos(env)
    dest = move(action, agent_pos)
    if !world[WALL, dest]
        world[AGENT, agent_pos] = false
        set_agent_pos!(env, dest)
        world[AGENT, dest] = true
    end

    if world[BALL, get_box_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

#####
# Common
#####

function (env::Union{TransportDirected{T}, TransportUndirected{T}})(::PickUp) where {T}
    world = get_world(env)
    agent_pos = get_agent_pos(env)

    if world[BALL, agent_pos]
        world[BALL, agent_pos] = false
        set_has_ball!(env, true)
    end

    if world[BALL, get_box_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end

function (env::Union{TransportDirected{T}, TransportUndirected{T}})(::Drop) where {T}
    world = get_world(env)
    agent_pos = get_agent_pos(env)

    if get_has_ball(env)
        set_has_ball!(env, false)
        set_ball_pos!(env, agent_pos)
        world[BALL, agent_pos] = true
    end

    if world[BALL, get_box_pos(env)]
        set_done!(env, true)
        set_reward!(env, get_terminal_reward(env))
    else
        set_done!(env, false)
        set_reward!(env, zero(T))
    end

    return nothing
end


text:
module Helpers

using CSV
using DataFrames
using Formatting
using StatsBase
using StatsModels

export get_coefs, allequal, featurestrings
export ys1d, roundclip, f1, log2file, xdecim

get_coefs(models::Vector{<:RegressionModel})::Vector{Vector{Float64}} = coef.(models)
function allequal(x::Vector)
    length(x) < 2 && return true
    e1 = x[1]
    @inbounds for i in 2:length(x)
        x[i] == e1 || return false
    end
    return true
end
function featurestrings(models::Vector{<:RegressionModel})::Vector{String}
    coefs = get_coefs(models)
    selected = [c .!= 0 for c in coefs]
    if allequal(selected)
        return [join([b ? "1" : "0" for b in selected[1]])]
    else
        return [join([b ? "1" : "0" for b in s]) for s in selected]
    end
end

ys1d(ymat::VecOrMat{Int})::Vector{Int} = ymat
function ys1d(ymat::Matrix{Float64})::Vector{Int}
    if size(ymat, 2) > 1
        [argmax(row) for row in eachrow(ymat)]
    else
        vec(1 * (ymat .>= .5))
    end
end
roundclip(vec::Vector{Float64}, range::UnitRange{Int})::Vector{Int} = max.(first(range), min.(last(range), round.(Int, vec)))

xdecim(num::Real, dec::Int) = sprintf1("%.$(dec)f", num)

function log2file(filename::String, x::String)
    path = "log/$filename.txt"
    dir, _ = splitdir(path)
    mkpath(dir)
    open(path, "a") do io
        write(io, "$x\n")
    end
end

function log2file(filename::String, df::DataFrame, header::Bool = false)
    CSV.write("log/$filename.txt", df, append = true, writeheader = header)
end
log2file(filename::String, mat::Matrix, names::Vector, header::Bool = false) = log2file(filename, DataFrame(mat, names), header)

end

using .Helpers

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
module nc

using ..lib

add_lib(Library(:ncurses,"libncursestw"))
add_lib(Library(:form, "libformtw"))
add_lib(Library(:menu, "libmenutw"))
add_lib(Library(:panel, "libpaneltw"))
add_lib(Library(:cdk, "libcdk.dylib.6.0.4"))

const COLOR_BLACK   = 0
const COLOR_RED     = 1
const COLOR_GREEN   = 2
const COLOR_YELLOW  = 3
const COLOR_BLUE    = 4
const COLOR_MAGENTA = 5
const COLOR_CYAN    = 6
const COLOR_WHITE   = 7

export COLOR_BLACK, COLOR_RED, COLOR_GREEN, COLOR_YELLOW, COLOR_BLUE, COLOR_MAGENTA, COLOR_CYAN, COLOR_WHITE

chtype = UInt8

COLOR_PAIR( n ) = NCURSES_BITS(n, 0)
export COLOR_PAIR

const NCURSES_ATTR_SHIFT =8
function NCURSES_BITS( m, shf)
    m << (shf + NCURSES_ATTR_SHIFT)
end
export NCURSES_BITS

const A_NORMAL     = uint32(0)
const A_ATTRIBUTES = ~uint32(0)
const A_CHARTEXT   = (NCURSES_BITS(uint32(1),0) - uint32(1))
const A_COLOR      = NCURSES_BITS((uint32(1) << 8) - uint32(1),0)
const A_STANDOUT   = NCURSES_BITS(uint32(1),8)
const A_UNDERLINE  = NCURSES_BITS(uint32(1),9)
const A_REVERSE    = NCURSES_BITS(uint32(1),10)
const A_BLINK      = NCURSES_BITS(uint32(1),11)
const A_DIM        = NCURSES_BITS(uint32(1),12)
const A_BOLD       = NCURSES_BITS(uint32(1),13)
const A_ALTCHARSET = NCURSES_BITS(uint32(1),14)
const A_INVIS      = NCURSES_BITS(uint32(1),15)
const A_PROTECT    = NCURSES_BITS(uint32(1),16)
const A_HORIZONTAL = NCURSES_BITS(uint32(1),17)
const A_LEFT       = NCURSES_BITS(uint32(1),18)
const A_LOW        = NCURSES_BITS(uint32(1),19)
const A_RIGHT      = NCURSES_BITS(uint32(1),20)
const A_TOP        = NCURSES_BITS(uint32(1),21)
const A_VERTICAL   = NCURSES_BITS(uint32(1),22)

export A_NORMAL, A_ATTRIBUTES, A_CHARTEXT, A_COLOR, A_STANDOUT, A_UNDERLINE, A_REVERSE, A_BLINK, A_DIM, A_BOLD, A_ALTCHARSET, A_INVIS, A_PROTECT, A_HORIZONTAL, A_LEFT, A_LOW, A_TOP, A_VERTICAL

# not really used. See getmouse() hack in ccall.jl
NCURSES_MOUSE_MASK(b,m) = m<<((b-1)*6) # NCURSES_MOUSE_VERSION=1
const NCURSES_BUTTON_PRESSED = uint32(2)
const BUTTON1_PRESSED = NCURSES_MOUSE_MASK( 1, NCURSES_BUTTON_PRESSED )
const BUTTON2_PRESSED = NCURSES_MOUSE_MASK( 2, NCURSES_BUTTON_PRESSED )
const BUTTON3_PRESSED = NCURSES_MOUSE_MASK( 3, NCURSES_BUTTON_PRESSED )
const BUTTON4_PRESSED = NCURSES_MOUSE_MASK( 4, NCURSES_BUTTON_PRESSED )
const REPORT_MOUSE_POSITION = NCURSES_MOUSE_MASK( 5, 8 ) # NCURSES_MOUSE_VERSION=1

export NCURSES_MOUSE_MASK, NCURSES_BUTTON_PRESSED, BUTTON1_PRESSED, BUTTON2_PRESSED, BUTTON3_PRESSED, BUTTON4_PRESSED, REPORT_MOUSE_POSITION

ncurses = get_lib(:ncurses)
panel = get_lib(:panel)
form = get_lib(:form)
menu = get_lib(:menu)
cdk = get_lib(:cdk)

export ncurses, panel, form, menu

NCURSES_VERSION = ""

export NCURSES_VERSION

function get_nc_version()
    nc_ver = ccall( dlsym( ncurses, :curses_version ), Ptr{UInt8}, ())
    if nc_ver != C_NULL
        NCURSES_VERSION = bytestring(nc_ver)
    end
end

function initscr()
    ccall( dlsym( ncurses, :initscr), Ptr{Void}, () )
end

function endwin()
    ccall( dlsym( ncurses, :endwin), Int, () )
end

function isendwin()
    ccall( dlsym( ncurses, :isendwin), Bool, () )
end

function newwin( lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( dlsym( ncurses, :newwin ), Ptr{Void}, ( Int, Int, Int, Int ),
        lines, cols, origy, origx )
end

function subwin( win::Ptr{Void}, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( dlsym( ncurses, :subwin ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end

function derwin( win, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( dlsym( ncurses, :derwin ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end

function delwin( win::Ptr{Void} )
    ccall( dlsym( ncurses, :delwin ), Void, ( Ptr{Void}, ), win )
end

function mvwaddch( win, y::Int, x::Int, c )
    ccall( dlsym( ncurses, :mvwaddch), Void,
        ( Ptr{Void}, Int, Int, Int ), win, y, x, c )
end

function mvwprintw( win::Ptr{Void}, row::Int, height::Int, fmt::String, str::String )
    ccall( dlsym( ncurses, :mvwprintw), Void,
        ( Ptr{Void}, Int, Int, Ptr{Uint8}, Ptr{Uint8}),
        win, row, height, fmt, str )
end

function printw( str::String )
    ccall( dlsym( ncurses, :printw), Void, (Ptr{Uint8},), str )
end

function mvprintw( row::Int, height::Int, str::String )
    ccall( dlsym( ncurses, :printw), Void, (Int, Int, Ptr{Uint8}), row, height, str )
end

function wprintw( win::Ptr{Void} , str::String )
    ccall( dlsym( ncurses, :printw), Void, (Ptr{Void}, Ptr{Uint8}), win, str )
end

function wmove( win, y::Int, x::Int )
    ccall( dlsym( ncurses, :wmove), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end

function wrefresh( win::Ptr{Void} )
    ccall( dlsym( ncurses, :wrefresh ), Void, ( Ptr{Void}, ), win )
end

function touchwin( win::Ptr{Void} )
    ccall( dlsym( ncurses, :touchwin ), Void, ( Ptr{Void}, ), win )
end

function refresh()
    ccall( dlsym( ncurses, :refresh ), Void, ( ) )
end

function erase()
    ccall( dlsym( ncurses, :erase ), Void, () )
end

function werase( win::Ptr{Void} )
    ccall( dlsym( ncurses, :werase ), Void, (Ptr{Void},), win )
end

function wclear( win::Ptr{Void} )
    ccall( dlsym( ncurses, :wclear ), Void, (Ptr{Void},), win )
end

function box( win, vchr, hchr )
    ccall( dlsym( ncurses, :box ), Void, (Ptr{Void}, Char, Char), win, vchr, hchr )
end

function wgetch( win::Ptr{Void } )
    ccall( dlsym( ncurses, :wgetch ), Uint32, (Ptr{Void},), win )
end

function keypad( win, bf )
    ccall( dlsym( ncurses, :keypad ), Int, (Ptr{Void}, Bool), win, bf )
end

function cbreak()
    ccall( dlsym( ncurses, :cbreak), Void, ( ) )
end

function nocbreak()
    ccall( dlsym( ncurses, :nocbreak), Void, ( ) )
end

function echo()
    ccall( dlsym( ncurses, :echo), Void, ( ) )
end

function noecho()
    ccall( dlsym( ncurses, :noecho), Void, ( ) )
end

function nodelay( win, bf )
    ccall( dlsym( ncurses, :nodelay ), Int, (Ptr{Void}, Bool), win, bf )
end

function raw()
    ccall( dlsym( ncurses, :raw ), Int, () )
end

function noraw()
    ccall( dlsym( ncurses, :noraw ), Int, () )
end

function notimeout( win, bf )
    ccall( dlsym( ncurses, :notimeout ), Int, (Ptr{Void}, Bool), win, bf )
end

function timeout( delay::Int )
    ccall( dlsym( ncurses, :timeout ), Void, (Int,), delay )
end

function wtimeout( win, delay::Int )
    ccall( dlsym( ncurses, :wtimeout ), Void, (Ptr{Void}, Int), win, delay )
end

# not standard but convenient
function getwinmaxyx( win )
    maxy = ccall( dlsym(ncurses, :getmaxy), Int, ( Ptr{Void}, ), win )
    maxx = ccall( dlsym(ncurses, :getmaxx), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function getwinbegyx( win )
    maxy = ccall( dlsym(ncurses, :getbegy), Int, ( Ptr{Void}, ), win )
    maxx = ccall( dlsym(ncurses, :getbegx), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function mvwin( win, y::Int, x::Int )
    ccall( dlsym( ncurses, :mvwin), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end

function beep()
    ccall( dlsym( ncurses, :beep ), Void, () )
end

function flash()
    ccall( dlsym( ncurses, :flash), Void, () )
end

function is_term_resized( lines::Int, cols::Int )
    ccall( dlsym( ncurses, :is_term_resized ), Bool, (Int, Int), lines, cols )
end

function wresize( win, lines::Int, cols::Int )
    ccall( dlsym( ncurses, :wresize), Int, (Ptr{Void}, Int, Int), win, lines, cols )
end

function getcuryx( win )
    cury = ccall( dlsym(ncurses, :getcury), Cint, ( Ptr{Void}, ), win )
    curx = ccall( dlsym(ncurses, :getcurx), Cint, ( Ptr{Void}, ), win )
    ( cury, curx )
end

function start_color()
    ccall( dlsym( ncurses, :start_color), Void, () )
end

function init_pair( pair, f, b )
    ccall( dlsym( ncurses, :init_pair ), Int, ( Int16, Int16, Int16 ), pair, f, b )
end

function init_color( color, r,g,b )
    ccall( dlsym( ncurses, :init_color ), Int, ( Int16, Int16, Int16, Int16 ),
        color, r,g,b )
end

function has_colors()
    ccall( dlsym( ncurses, :has_colors), Bool, () )
end

function wattroff( win, attrs )
    ccall( dlsym(ncurses, :wattroff), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end

function wattron( win, attrs )
    ccall( dlsym(ncurses, :wattron), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end

function wattrset( win, attrs )
    ccall( dlsym(ncurses, :wattrset), Int, ( Ptr{Void}, Uint32 ), win, attrs )
end

function wbkgdset( win, ch )
    ccall( dlsym(ncurses, :wbkgdset ), Void, ( Ptr{Void}, Uint32 ), win, ch )
end

function wbkgd( win, ch )
    ccall( dlsym(ncurses, :wbkgd ), Void, ( Ptr{Void}, Uint32 ), win, ch )
end

function curs_set( vis )
    ccall( dlsym(ncurses, :curs_set), Int, ( Int, ), vis )
end

function has_mouse()
    ccall( dlsym(ncurses, :has_mouse), Bool, () )
end

function mousemask( mask )
    oldmm = Array( Uint64, 1 )
    resultmm = ccall( dlsym( ncurses, :mousemask), Uint64, (Uint64, Ptr{Uint64}), mask, oldmm )
    ( resultmm, oldmm[1])
end

function mouseinterval( n::Int )
    ccall( dlsym(ncurses, :mouseinterval), Int, (Int, ) , n)
end

#hack!
const mouseByteString = bytestring( Array( Uint8, 64 ) )
function getmouse()
    #=
    type Mouse_Event_t
        id::Int8 # short
        x::Int32 # int
        y::Int32 # int
        z::Int32 # int
        bstate::Uint64 # unsigned long
    end
    =#
    # the 5th byte is x, 9th byte is y
    # 19th byte is x08 if scroll up, 20th byte is x08 if scroll down
    # 17th byte is x02 if button 1 pressed
    # 17th byte is x01 if button 1 released
    # 17-18th is 0xfffd if mousewheel is pressed down
    ccall( dlsym( ncurses, :getmouse), Int, (Ptr{Uint8}, ), mouseByteString )
    bs = mouseByteString
    x = uint8(bs[5])
    y = uint8(bs[9])
    state=:unknown
    println(bs, x, y, state)
    # if bs[17] & 0x02 != 0
    #     state = :button1_pressed
    # elseif bs[19] & 0x08 != 0
    #     state = :scroll_up
    # elseif bs[20] & 0x08 != 0
    #     state = :scroll_down
    # end
    ( state, x, y, bs )
end

function baudrate()
    ccall( dlsym( ncurses, :baudrate), Int, () )
end

function clearok( win, bf )
    ccall( dlsym( ncurses, :clearok), Int, ( Ptr{Void}, Bool ), win, bf )
end

function immedok( win, bf )
    ccall( dlsym( ncurses, :immedok), Int, ( Ptr{Void}, Bool ), win, bf )
end

function napms( ms::Int )
    ccall( dlsym( ncurses, :napms), Int, (Int, ), ms )
end

#===== PANEL library ====#

function update_panels()
    ccall( dlsym( panel, :update_panels), Void, () )
end

function doupdate()
    ccall( dlsym( panel, :doupdate ), Void, () )
end

function new_panel( win::Ptr{Void} )
    ccall( dlsym( panel, :new_panel ), Ptr{Void}, ( Ptr{Void}, ), win )
end

function top_panel( pan::Ptr{Void} )
    ccall( dlsym( panel, :top_panel ), Int, ( Ptr{Void}, ), pan )
end

function bottom_panel( pan::Ptr{Void} )
    ccall( dlsym( panel, :bottom_panel ), Int, ( Ptr{Void}, ), pan )
end

function move_panel( pan::Ptr{Void}, starty::Int, startx::Int )
    ccall( dlsym( panel, :move_panel ), Ptr{Void}, ( Ptr{Void}, Int, Int ), pan, starty, startx )
end

function del_panel( panel::Ptr{Void} )
    ccall(dlsym( panel, :del_panel ), Void, (Ptr{Void}, ), panel )
end

function show_panel( panel::Ptr{Void} )
    ccall(dlsym( panel, :show_panel ), Int, (Ptr{Void}, ), panel )
end

function hide_panel( panel::Ptr{Void} )
    ccall(dlsym( panel, :hide_panel ), Int, (Ptr{Void}, ), panel )
end

function panel_hidden( panel::Ptr{Void } )
    ccall(dlsym( panel, :panel_hidden ), Int, (Ptr{Void}, ), panel ) != 0
end

function replace_panel( panel::Ptr{Void}, window::Ptr{Void} )
    ccall(dlsym( panel, :replace_panel ), Void, (Ptr{Void}, Ptr{Void} ), panel, window )
end

function set_panel_userptr( p1::Ptr{Void}, p2::Ptr{Void} )
    ccall(dlsym( panel, :set_panel_userptr), Void, (Ptr{Void}, Ptr{Void}), p1, p2 )
end

function wborder( win::Ptr{Void},ls::chtype,rs::chtype,ts::chtype,bs::chtype,tl::chtype,tr::chtype,bl::chtype,br::chtype)
    ccall(dlsym( panel, :wborder), Int, (Ptr{Void}, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8, Uint8), ls, rs,ts, bs, tl, tr, bl, br)
end

export initscr, endwin, isendwin, newwin, subwin, delwin, derwin, werase, erase, refresh, mvwaddch, wclear, mvwprintw, touchwin, wmove, box, wrefresh, wgetch, nocbreak, keypad, nodelay, noecho, cbreak, echo, raw, noraw, timeout, notimeout, mvwin, getwinbegyx, beep, wtimeout, flash, getwinmaxyx, is_term_resized, wresize, getcuryx, start_color, init_pair, init_color, has_colors, wattroff, wattron, wattrset, wbkgdset, wbkgd, curs_set, has_mouse, mousemask, mouseinterval, getmouse, mouseByteString, baudrate, clearok, immedok, napms, update_panels, doupdate, new_panel, top_panel, bottom_panel, hide_panel, move_panel, del_panel, show_panel, panel_hidden, replace_panel, set_panel_userptr

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
function plotquality(df)
    legendmarkers["feeder"] = (color = :black, marker = '▴', strokecolor = :black, markerstrokewidth = 0.5, strokewidth = 0.5, markersize = 15px)
    legendmarkers["pickup"] = (color = :black, marker = '↑', strokecolor = :black, markerstrokewidth = 0.5, strokewidth = 0.5, markersize = 15px)
    legendmarkers["dropoff"] = (color = :black, marker = '↓', strokecolor = :black, markerstrokewidth = 0.5, strokewidth = 0.5, markersize = 15px)
    delete!(legendmarkers, "mean ± FWHM")
    delete!(legendmarkers, "center of search")
    polys = OrderedDict(string(df.group[1]) => (color = df.groupcolor[1], strokecolor = :transparent))
    for r in eachrow(df)
        c = r.groupcolor
        scene, layout = layoutscene()#0, resolution = (2max_width, 900.0))
        ax = layout[1,1] = LAxis(scene, aspect = DataAspect())
        lines!(ax, homing(r.track))#; legendmarkers["track"]..., color = :red)
        lines!(ax, searching(r.track); legendmarkers["track"]..., color = c)
        scatter!(ax, [turningpoint(r.track)]; legendmarkers["turning point"]..., color = RGBA(c, 0.75))
        # scatter!(ax, [searchcenter(r.track)]; legendmarkers["center of search"]..., color = RGBA(c, 0.75))
        if !ismissing(r.feeder)
            scatter!(ax, [r.feeder]; legendmarkers["feeder"]..., color = RGBA(c, 0.75))
        end
        if !ismissing(r.nest)
            scatter!(ax, [r.nest]; legendmarkers["burrow"]..., color = RGBA(c, 0.75))
        end
        if !ismissing(r.pickup)
            scatter!(ax, [r.pickup]; legendmarkers["pickup"]..., color = RGBA(c, 0.75))
        end
        scatter!(ax, [r.dropoff]; legendmarkers["dropoff"]..., color = RGBA(c, 0.75))
        scatter!(ax, [r.fictive_nest]; legendmarkers["fictive burrow"]...)#, color = RGBA(c, 0.75))
        layout[2, 1] = LLegend(scene, apply_element(values(legendmarkers)), collect(keys(legendmarkers)), orientation = :horizontal, nbanks = 2, tellheight = true, height = Auto(), groupgap = 30);
        limits!(ax, (-130, 130), (-200, 20))
        FileIO.save("$(first(splitext(r.comment))).pdf", scene)
    end
end


text:
#Worked with Kasra#
                            #Question No.1#
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables

#part (l)#
function q1()
#Part (a)#
Random.seed!(1234)
A=rand(Uniform(-5,10),(10,7))
B=rand(Normal(-2,15),(10,7))
C=[A[1:5,1:5] B[1:5,6:7]]
D=(A.<=0).*A

#Part (b)#
length(A)

#Part (c)#
length(unique(D))

#Part (d)#
E=reshape(B,:,1)
E=vec(B) #Alternative: Trese command for vectorization

#Part (e)#
F=cat(A, B;dims=3)

#Part (f)#
F=permutedims(F,(3,1,2)) # 3 1 2 are dimentions

#Part (g)#
G=kron(B,C)
kron(C,F) #Kronecher product is not possible since dimentions of C & F do not match

#Part (h)#
cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\0 PS Solutions_Julia Directory")
JLD2.@save "matrixpractice.jld2" A B C D E F G

#Part (i)#
JLD2.@save "firstmatrix.jld2" A B C D

#Part (j)#
CC=convert(DataFrame, C) # df=DataFrame(C) --> Better code which Kasra told me
CSV.write("Cmatrix.csv", CC) # CSV.write("Cmatrix.csv", df)

#Part (k)#
DD=convert(DataFrame, D) # dfD=DataFrame(D) --> Better code by Kasra
CSV.write("Dmatrix.dat", DD, delim='\t') #CSV.write("Dmatrix.dat",df,delim='\t')

#Part (l) --> done
      return A,B,C,D
end

                            #Question No.2#
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables

#part (f)
function q2(a,b,c)
#Part (a)#
AB=zeros(size(A))
for i=1:length(A) #Alternatively define n=length(A)
       AB[i]=A[i]*B[i]
       end
AB2=A.*B

#Part (b)#
Cprime=[]
for i=1:length(C)
      if -5<=C[i]<=5
            append!(Cprime,C[i])    #Took help from Kasra to figure out the loops
      end
end
Cprime2=[i for i in C if -5<=i<=5]

#Part (c)#
X1=cat(ones(15169,5), [rand(Bernoulli(0.75*(6-t)/5)) for n in 1:15169,t in 1:5], [rand(Normal(15+t-1,5(t-1))) for n in 1:15169,t in 1:5], [rand(Normal(π*(6-t)/3,1/ℯ)) for n in 1:15169,t in 1:5], rand(Binomial(20,0.6),15169,5), rand(Binomial(20,0.5),15169,5);dims=3)
X=permutedims(X1,(1,3,2)) # Or replace X1 with the entire code of X1

#Part (d)#
β=hcat([[1+.25*(t-1), log(t), -sqrt(t), ℯ^t-ℯ^(t+1), t, t/3] for t in 1:5]...)
#Part (e)#
Y=hcat([X[:,:,i]*β[:,i].+rand(Normal(0,.36)) for i in 1:5]...)
#part(f)#
      return
end

                            #Question No.3#
#Cleared Julia Manually by restarting#
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables

#part (g)#
function q3()

#Part (a)#
cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\0 PS Solutions_Julia Directory")
df=CSV.read("nlsw88.csv", header=true, missingstring="")
JLD2.@save "nlsw88.jld" df

#Part (b)#
describe(df)
describe(df[:,:never_married])
# Since never_married is a binary variable which takes the value of 1 to represent
# people who never married; their precentage would be meean*100 = 10.4185% #
describe(df[:,:collgrad])
# By the same logic as above, the percentage of college grads would be 23.6866%#

#Part (c)#
tabulate(df,:race)
#Category 1 = 72.885% | Category 2 = 25.957% | Category 3 = 1.158%
#Discussed with Kasra and Waleed, and used TexTables package for tabulate command#

#Part (d)#
descriptives=describe(df, :mean, :median, :std, :min, :max, :nunique, :q75, :q25)
interquartile_range= descriptives[:,:q75] - descriptives[:,:q25]
#Cannot figure out how to find missing observations#

#Part (e)#
freqtable(df[:,:industry],df[:,:occupation])

#Part (f)#
df_f=df[:,[:industry, :occupation , :wage]]
gdf=groupby(df_f,[:industry, :occupation])
gdf=combine(gdf, valuecols(gdf) .=> mean)
freqtable(gdf,:wage_mean,:industry,:occupation);
#Part (g)#
      return
end
q1()=A,B,C,D
q2(A,B,C)
q3()

                              #Question No.4#
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables

#Part (h)#
function q4()
cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\0 PS Solutions_Julia Directory")

#Part (a)#
JLD2.@load "firstmatrix.jld2" A B C D

#Part (b+c+e)#
function matrixops(A,B)
      if if size(A)==size(B)
      AdotB=A.*B
      AtransB=A'*B
      sumAB=sum(A+B)
      return AdotB, AtransB, sumAB
else
      println("inputs must have the same size.")
end
end

#Part (d)#
matrixops(A,B)

#Part (f)#
matrixops(C,D)
#Does not work because inputs must have the same size.

#Part (g)#
JLD2.@load "nlsw88.jld"
matrixops(convert(Array,df.ttl_exp),convert(Array,df.wage));

#Part (h)#
      return
end
end
q4()
                              #THE END#


text:
julia_version = v"1.3.1"
gcc_version = v"8"
include("../common.jl")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Flux, Flux.Tracker
using Distributions
import LinearAlgebra.logabsdet
import SpecialFunctions.lgamma
# include("sbt.jl")

struct VP
  m::TrackedArray
  log_s::TrackedArray
  real_w::TrackedArray

  VP(K::Integer) = new(param(randn(K, 2)), param(randn(K, 2)), param(zeros(K - 1, 2)))
end

function lpdf_vp(p, x)
  return sum(logpdf.(Normal.(p[:, 1], exp.(p[:, 2])), x))
end

function rsample(p)
  return p[:, 1] .+ randn(size(p, 1)) .* exp.(p[:, 2])
end

function logsumexp(logx)
  mx = maximum(logx)
  return mx + log(sum(exp.(logx .- mx)))
end

function logsumexp(logx; dims::Integer)
  mx = maximum(logx, dims=dims)
  return mx .+ log.(sum(exp.(logx .- mx), dims=dims))
end

function logsumexp0p(logx)
  return log1p(sum(exp.(logx)))
end

function lpdf_logx(dist, logx)
  return logpdf(dist, exp(logx)) + logx
end

function softmax_safe(x)
  logsumexp_x = logsumexp(x)
  return exp.(x .- logsumexp_x)
end

"""
Transform a real vector to a simplex. The sum of the resulting simplex
is less than 1, and one minus that sum is the probability of the first
element in the K+1 dimensional simplex.
"""
function softmax_fullrank(x; complete::Bool=true)
  logsumexp0p_x = logsumexp0p(x)
  p_reduced = exp.(x .- logsumexp0p_x)
  if complete
    K = length(p_reduced) + 1
    g =  vcat(1 .- sum(p_reduced, dims=1), p_reduced)
    return g
  else
    return p_reduced
  end
end

function lpdf_Dirichlet_fullrank(a, p)
  out = lgamma(sum(a)) - sum(lgamma.(a)) + sum((a[2:end] .- 1.0) .* log.(p))
  out += (a[1] - 1.0) * log1p(-sum(p))
  return out
end

# TODO: Test
function lpdf_real_simplex(alpha, real_simplex)
  dim = length(real_simplex)
  p = softmax_fullrank(real_simplex, complete=false)
  jacobian = [i == j ? p[i] * (1 - p[i]) : -p[i] * p[j] for i in 1:dim, j in 1:dim]
  return lpdf_Dirichlet_fullrank(alpha, p) + logabsdet(jacobian)[1]

  # prob = sbt(real_simplex)
  # return logpdf(Dirichlet(alpha), prob) + sum(sbt_logabsdetJ(real_simplex, prob))
end

# loglike
function loglike(y, w, m, s)
  N = length(y)

  # NOTE: This is slow because of how tracker is implemented.
  #       Tensor expressions are still needed.
  # out = 0.0
  # for i in 1:N
  #   out += sum(logsumexp(log.(w) .+ logpdf.(Normal.(m, s), y[i]), dims=1))
  # end
  # return out

  # NOTE: Preferred Implementation
  K = length(m)
  log_w = reshape(log.(w), 1, K)
  lpdf = logpdf.(Normal.(reshape(m, 1, K), reshape(s, 1, K)), y)
  return sum(logsumexp(log_w .+ lpdf, dims=2))

  # NOTE: Tracker can't keep track of reshapes of Distribution objects
  #       i.e. REPLACING THE PREVIOUS TWO LINES WITH THE FOLLOWING WILL NOT
  #       WORK!!!
  # lpdf = logpdf.(reshape(Normal.(m, s), 1, K), reshape(y, N, 1))
  # return sum(logsumexp(log_w .+ lpdf, dims=2))
end

# log_p
function log_p(real_w, m, log_s)
  K = length(m)
  log_p_m = sum(logpdf.(Normal(0, 1), m))
  log_p_log_s = sum(lpdf_logx.(LogNormal(0, 1), log_s))
  log_p_real_w = lpdf_real_simplex(ones(K) / K, real_w)
  return log_p_m + log_p_log_s + log_p_real_w
end

# log_q
function log_q(real_w, m, log_s, vp::VP)
  log_q_m = lpdf_vp(vp.m, m)
  log_q_log_s = lpdf_vp(vp.log_s, log_s)
  log_q_real_w = lpdf_vp(vp.real_w, real_w)
  return log_q_m + log_q_log_s + log_q_real_w
end

# ELBO
function elbo(y, vp::VP)
  m = rsample(vp.m)

  log_s = rsample(vp.log_s)
  s = exp.(log_s)

  real_w = rsample(vp.real_w)
  w = softmax_fullrank(real_w)

  return loglike(y, w, m, s) + log_p(real_w, m, log_s) - log_q(real_w, m, log_s, vp)
end


text:
# Gumbel distribution

export Gumbel

@parameterized Gumbel(μ,σ) ≃ Lebesgue(ℝ)

@kwstruct Gumbel()

@μσ_methods Gumbel()

function logdensity(d::Gumbel{()} , x)
    return -exp(-x) - x
end

import Base

function Base.rand(rng::AbstractRNG, d::Gumbel{()})
    u = rand(rng)
    -log(-log(u))
end

TV.as(::Gumbel) = asℝ

≪(::Gumbel, ::Lebesgue{X}) where X <: Real = true

distproxy(::Gumbel{()}) = Dists.Gumbel()


text:
##

using ACE, JuLIP, ACEatoms, ACEbase, Test, LinearAlgebra
using ACE: evaluate, evaluate_d, SymmetricBasis, SimpleSparseBasis, PIBasis
using ACEbase.Testing: fdtest

##

# construct the 1p-basis
@info("Construcing a Linear ACE Model")
maxdeg = 8
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

species = [:Ti, :Al]
B1p = ACEatoms.ZμRnYlm_1pbasis(; species = species, maxdeg=maxdeg, Bsel = Bsel,
                                 rin = 1.2, rcut = 5.0)

ACE.init1pspec!(B1p, Bsel)
basis = SymmetricBasis(ACE.Invariant(), B1p, Bsel)
cTi = randn(length(basis))
cAl = randn(length(basis))
models = Dict(:Ti => ACE.LinearACEModel(basis, cTi; evaluator = :standard),
              :Al => ACE.LinearACEModel(basis, cAl; evaluator = :standard) )

@info("Convert LinearACEModel into an ACESitePotential")
V = ACEatoms.ACESitePotential(models)

## FIO

@info("Check FIO")
using ACEbase.Testing: test_fio
println(@test(all(test_fio(V; warntype = false))))

##

@info("Create random TiAl configuration")
zTi = AtomicNumber(:Ti)
zAl = AtomicNumber(:Al)
at = bulk(:Ti, cubic=true) * 3
at.Z[2:3:end] .= zAl
at = rattle!(at, 0.1)

##

@info("Test that evaluate(model) ≈ evaluate(site-potential)")
nlist = neighbourlist(at, cutoff(V))
env, _ = ACEatoms.environment(V, at, nlist, 2)
Js, Rs, Zs = JuLIP.Potentials.neigsz(nlist, at, 2)
z0 = at.Z[2]
sym = chemical_symbol(z0)
env1 = ACEatoms.environment(V, Rs, Zs, z0)
ACEbase._allfieldsequal(env1, env)
v1 = evaluate(models[sym], env).val
v2 = evaluate(V, Rs, Zs, z0)
println(@test v1 ≈ v2)
energy(V, at)

##
@info("Test that grad_config(model) ≈ evaluate_d(site-potential)")
dEs = zeros(JVecF, length(env))
dv1 = ACE.grad_config!(dEs, V.models[z0], env)
dv2 = JuLIP.evaluate_d!(dEs, nothing, V, Rs, Zs, z0)
dv3 = evaluate_d(V, Rs, Zs, z0)

println(@test(dv1 ≈ dv2))
println(@test(dv1 ≈ dv3))

forces(V, at)
##

@info("Finite-difference test on total energy")
println(@test JuLIP.Testing.fdtest(V, at))

##

@info("Check also the ACESitePotentialBasis interface")

ipbasis = ACEatoms.basis(V)
cc = [cTi; cAl]

# check get_params, set_params

@info("  ... energy")
val1 = energy(V, at)
val2 = sum(cc .* energy(ipbasis, at))
println(@test (val1 ≈ val2))

@info("  ... forces")
val1 = forces(V, at)
frcB = forces(ipbasis, at)
val2 = sum(cc .* frcB)
println(@test val1 ≈ val2)

@info("  ... virial")
val1 = virial(V, at)
virB = virial(ipbasis, at)
val2 = sum(cc .* virB)
println(@test val1 ≈ val2)