text:
include("../../HamiltonODE.jl")

using Plots
pyplot()

d = 2
m = (1)
q0 = [1.,0.]
p0 = [0.,1.]

F(q)= -q
U(q)= q⋅q
sol(q)=[cos(q),sin(q)]

HS = HamiltonSystem(d,m,q0,p0,F,U,sol)

T=2*pi
dt = 1e-1

t,Q,P=integrate(HS,Integrator(verlet_step,dt,T));

p1 = plot(t,Q')
p2 = plot(Q[1,:],Q[2,:])
plot(p1,p2)

max_stars_repo_name: BoergeSt/HamiltonODE

max_issues_repo_name: BoergeSt/HamiltonODE

max_forks_repo_name: BoergeSt/HamiltonODE


text:
"""
	$(SIGNATURES)

Step function for admission probabilities.

There are `n` endowment groups (e.g. GPA quartiles) defined by upper bounds for a score.
There are `nc` colleges.
For each college / endowment pair, some entry probs are fixed (e.g. close to 0 or 1); others are calibrated.
"""
mutable struct AdmProbFctStepSwitches{F1} <: AbstractAdmProbFctSwitches{F1}
    pvec :: ParamVector
    nColleges :: Int
    nOpenColleges :: Int
    # nGroups :: Int
    cutoffV :: Vector{F1}
    # caSwitches :: CalibratedArraySwitches{F1, 2}
end

mutable struct AdmProbFctStep{F1} <: AbstractAdmProbFct{F1}
    objId :: ObjectId
    switches :: AdmProbFctStepSwitches{F1}
    # calArray :: CalibratedArray{F1, 2}
    probMatrix :: Matrix{F1}
end

Lazy.@forward  AdmProbFctStepSwitches.pvec (
    ModelParams.get_object_id
    );
Lazy.@forward AdmProbFctStep.switches (
    ModelParams.get_pvector, ModelParams.has_pvector
    );

# ModelParams.get_object_id(switches :: AdmProbFctStepSwitches{F1}) where F1 =
#     get_object_id(switches.pvec);
ModelParams.has_pvector(switches :: AdmProbFctStepSwitches{F1}) where F1 = true;

param_names(switches :: AdmProbFctStepSwitches{F1}) where F1 = nothing;

# Cutoffs: group upper bounds (there is another group with cutoff Inf)
function step_adm_prob(x :: F1, cutoffV, probV) where F1 <: Real
    return probV[find_class(x, cutoffV)]
end

step_adm_prob(xV, cutoffV, probV) = [step_adm_prob(x, cutoffV, probV)  for x in xV];

## ----------  Constructors

# defaultProbM by [group, college]
function init_admprob_fct_step_switches(
    objId :: ObjectId, cutoffV :: AbstractVector{F1},
    defaultProbM :: Matrix{F1},  isCalM :: Matrix{Bool}
    ) where F1

    @assert length(cutoffV) == (size(defaultProbM, 1) - 1);
    @assert all_at_least(defaultProbM, zero(F1));
    @assert all_at_most(defaultProbM, one(F1));

    lbM = fill(0.01, size(defaultProbM));
    ubM = fill(0.99, size(defaultProbM));
    pCA = CalArray(:probMatrix, "Admission probabilities", "admProb",
        copy(defaultProbM), defaultProbM, lbM, ubM, isCalM);
    pvec = ParamVector(objId, [pCA]);
    # caSwitches = CalibratedArraySwitches(
    #     make_child_id(objId, :probMatrix),
    #     defaultProbM, lbM, ubM, isCalM);

    nc = size(defaultProbM, 2);
    nOpen = 1;
    return AdmProbFctStepSwitches(pvec, nc, nOpen, cutoffV)
    # return AdmProbFctStepSwitches(objId, nc, nOpen, cutoffV, caSwitches)
end

function make_test_admprob_fct_step_switches(nc)
    rng = MersenneTwister(43);
    ng = nc + 1;
    cutoffV = collect(LinRange(0.2, 0.9, ng-1));
    isCalM = Matrix{Bool}(rand(rng, ng, nc) .> 0.5);
    defaultProbM = LinRange(0.1, 0.8, ng) .+ LinRange(0.0, 0.1, nc)';
    # Open college
    isCalM[:, 1] .= false;
    defaultProbM[:, 1] .= 0.99;
    switches =
        init_admprob_fct_step_switches(ObjectId(:test), cutoffV, defaultProbM, isCalM);
    @assert validate_admprob_fct_switches(switches);
    return switches
end

function init_admprob_fct(switches :: AdmProbFctStepSwitches{F1}) where F1
    # calArray = CalibratedArray(switches.caSwitches);
    probM = copy(param_value(switches, :probMatrix));
    af = AdmProbFctStep(get_object_id(switches), switches, probM);
    @assert validate_admprob_fct(af);
    return af
end

function make_test_admprob_fct_step(nc)
    switches = make_test_admprob_fct_step_switches(nc);
    return init_admprob_fct(switches)
end

function validate_admprob_fct_switches(switches :: AdmProbFctStepSwitches{F1}) where F1
    isValid = true;
    return isValid;
end

function validate_admprob_fct(af :: AdmProbFctStep{F1}) where F1
    isValid = true;
    return isValid;
end

## ------------  Adm prob for one college

function make_admprob_function(af :: AdmProbFctStep{F1},
    iCollege :: Integer) where F1

    if iCollege <= n_open_colleges(af)
        fct = x -> 1.0;
    else
        # prob_gcM = ModelParams.values(af.calArray);
        prob_gcM = af.probMatrix;
        fct = x -> step_adm_prob(x, af.switches.cutoffV, prob_gcM[:, iCollege]);
    end
    return fct
end

function prob_admit(af :: AdmProbFctStep{F1}, iCollege :: Integer, hsGpa) where F1
    if iCollege <= n_open_colleges(af)
        probV = prob_admit_open(hsGpa);
    else
        prob_gcM = af.probMatrix; #  ModelParams.values(af.calArray);
        probV = step_adm_prob(hsGpa, af.switches.cutoffV, prob_gcM[:, iCollege]);
    end
    return probV
end

# ----------------

ext: jl


text:
"""
    BCType

Abstract supertype for boundary condition types.
"""
abstract type BCType end

"""
    Periodic

A type specifying a periodic boundary condition.

A condition may not be specified with a `Periodic` boundary condition.
"""
struct Periodic <: BCType end

"""
    Flux

A type specifying a boundary condition on the flux of a field.
"""
struct Flux <: BCType end

"""
    Gradient

A type specifying a boundary condition on the derivative or gradient of a field. Also
called a Neumann boundary condition.
"""
struct Gradient <: BCType end

"""
    Value

A type specifying a boundary condition on the value of a field. Also called a Dirchlet
boundary condition.
"""
struct Value <: BCType end

"""
    NoPenetration

A type specifying a no-penetration boundary condition for a velocity component that is normal to a wall.

Thus `NoPenetration` can only be applied to `u` along x, `v` along y, or `w` along z.
For all other cases --- fields located at (Cell, Cell, Cell), or `u`, `v`,
and `w` in (y, z), (x, z), and (x, y), respectively, either `Value`,
`Gradient`, or `Flux` conditions must be used.

A condition may not be specified with a `NoPenetration` boundary condition.

Note that this differs from a zero `Value` boundary condition as `Value` imposes values at the cell centers
(and could apply to tracers) while a no-penetration boundary condition only applies
to normal velocity components at a wall, where the velocity at the cell face collocated
at the wall is known and set to zero.
"""
struct NoPenetration <: BCType end


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

function GitBlame(repo::GitRepo, path::AbstractString; options::BlameOptions=BlameOptions())
    blame_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_blame_file, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Ptr{BlameOptions}),
                   blame_ptr_ptr, repo.ptr, path, Ref(options))
    return GitBlame(repo, blame_ptr_ptr[])
end

function counthunks(blame::GitBlame)
    return ccall((:git_blame_get_hunk_count, :libgit2), Int32, (Ptr{Void},), blame.ptr)
end

function Base.getindex(blame::GitBlame, i::Integer)
    if !(1 <= i <= counthunks(blame))
        throw(BoundsError(blame, (i,)))
    end
    hunk_ptr = ccall((:git_blame_get_hunk_byindex, :libgit2),
                      Ptr{BlameHunk},
                      (Ptr{Void}, Csize_t), blame.ptr, i-1)
    return unsafe_load(hunk_ptr)
end

function Base.show(io::IO, blame_hunk::BlameHunk)
    println(io, "GitBlameHunk:")
    println(io, "Original path: ", blame_hunk.orig.path)
    println(io, "Lines in hunk: ", blame_hunk.lines_in_hunk)
    println(io, "Final commit oid: ", blame_hunk.final_commit_id)
    print(io, "Final signature: ")
    show(io, blame_hunk.final_signature)
    println()
    println(io, "Original commit oid: ", blame_hunk.orig_commit_id)
    print(io, "Original signature: ")
    show(io, blame_hunk.orig_signature)
    println()
end

max_stars_repo_name: nassarhuda/julia

max_issues_repo_name: nassarhuda/julia

max_forks_repo_name: nassarhuda/julia


text:
module NHST
    export t_test

    export test_statistic, p_value, ci
    export BinomialTest, ChisqTest, TTest

    using Distributions

	# :twosided, :greater, :less

    include("types.jl")
    include("binomial.jl")
    include("chisq.jl")
    include("ttest.jl")
end


text:
module TablesTests

using Test
using Gridap
using Gridap.Arrays
using Gridap.Io
using JSON
using FillArrays

data = Float64[2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = [1,4,4,7,13]
a = Table(data,ptrs)
b = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(a,b)

c = convert(Table{Float64,Vector{Float64},Vector{Int32}},a)
test_array(c,b)

data = Fill(1.3,12)
ptrs = [1,4,4,7,13]
d = Table(data,ptrs)
e = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(d,e)

data = Int[2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = [1,4,4,7,13]
a = Table(data,ptrs)

data  = reinterpret(Int,Vector{Float64}(undef,12))
data[1:6] .= a.data[7:12]
data[7:9] .= a.data[4:6]
data[10:12] .= a.data[1:3]

perm = Vector{Int}(undef,12)
perm[1:3]  .= 10:12
perm[4:6]  .= 7:9
perm[7:12] .= 1:6
data = lazy_map(Reindex(data),perm)
b = Table(data,ptrs)
test_array(a,b)

vv = Array{Array{Int,2},2}(undef,2,2)
vv[1,1] = rand(1:10,2,2)
vv[2,1] = rand(1:10,2,2)
vv[2,2] = rand(1:10,2,2)
vv[1,2] = rand(1:10,2,2)
b = [[ vv[i][j] for j in 1:length(vv[i]) ] for i in 1:length(vv)]
a = Table(vv)
test_array(a,b)

data = Float64[2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = Int32[1,4,4,7,13]
a = Table(data,ptrs)
b = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(a,b)

data = Float64[]
ptrs = [1,]
a = Table(data,ptrs)
b = [ data[ptrs[i]:ptrs[i+1]-1] for i in 1:length(ptrs)-1]
test_array(a,b)

vv = [[1,2,3],[2,3],[5,8],Int[],[1,2,4]]
a = Table(vv)
test_array(a,vv)

_data, _ptrs = generate_data_and_ptrs(vv)

data = [1, 2, 3, 2, 3, 5, 8, 1, 2, 4]
ptrs = [1, 4, 6, 8, 8, 11]

@test data == _data
@test ptrs == _ptrs

a = [9,2,1,2,4,7,4]
b = [1,9,2,1,2,4,7]

rewind_ptrs!(a)
@test a == b

a = [3,2,4,2]
b = [1,3,7,9]

length_to_ptrs!(a)
@test a == b

pa = [1,3,5,7,9]
pb = [1,3,5,7]

pc = append_ptrs(pa,pb)

@test pa == [1,3,5,7,9]
@test pb == [1,3,5,7]
@test pc == [1, 3, 5, 7, 9, 11, 13, 15]

vv1 = [[1,2,3],[2,3],[5,8],Int[],[1,2,4]]
table1 = Table(vv1)
vv2 = [[1,3],[4,2,3],Int[],Int[],[1,2,4]]
table2 = Table(vv2)

table3 = append_tables_globally(table1,table2)
@test table3 == vcat(table1,table2)

table4 = append_tables_locally((0,5),(table1,table2))
@test table4 == [[1, 2, 3, 6, 8], [2, 3, 9, 7, 8], [5, 8], Int[], [1, 2, 4, 6, 7, 9]]

a_to_lb_to_b = [[1,2,3],[2,3],[5,8],[2],[1,2,4]]
a_to_lb_to_b = Table(a_to_lb_to_b)
a_to_lb = fill(1,length(a_to_lb_to_b))
a_to_b = lazy_map(getindex,a_to_lb_to_b,a_to_lb)
r = [ lb_to_b[a_to_lb[a]] for (a,lb_to_b) in enumerate(a_to_lb_to_b) ]
test_array(a_to_b,r)

b_to_la_to_a = [[5,1,4,2,3],[1,2,3],[5,4],[2,4,5,3],[5,1,2,4]]
b_to_la_to_a = Table(b_to_la_to_a)
         #1,2,3,4,5
a_to_b = [2,5,4,3,1]
a_to_la = find_local_index(a_to_b,b_to_la_to_a)
r = [1,3,4,2,1]
test_array(a_to_la,r)

a = [[5,1,4,2,3],[1,2,3],[5,4],[2,4,5,3],[5,1,2,4]]
t = Table(a)
r = 2:4
u = t[r]
b = a[r]
@test isa(u,Table)
test_array(u,b)

r = [2,1,3]
u = t[r]
b = a[r]
@test isa(u,Table)
test_array(u,b)

t = empty_table(3)
r = Vector{Int}[[], [], []]
test_array(t,r)

a_to_bs = Table([[1,2,3],[7,8],[4,5,6]])
b_to_a = flatten_partition(a_to_bs)
r = [1, 1, 1, 3, 3, 3, 2, 2]
test_array(b_to_a,r)

data = Float64[2,3,1,3,6,7,3,2,5,6,3,4]
ptrs = [1,4,4,7,13]
a = Table(data,ptrs)

dict = to_dict(a)
b = from_dict(Table{Float64,Vector{Float64},Vector{Int32}},dict)
@test a == b

s = to_json(a)
b = from_json(Table{Float64,Vector{Float64},Vector{Int32}},s)
@test a == b

d = mktempdir()
f = joinpath(d,"a.jld2")

to_jld2_file(a,f)
@test a == from_jld2_file(typeof(a),f)

rm(d,recursive=true)

end # module


text:
function create(c::Client, ::Type{DiscordChannel}, g::AbstractGuild; kwargs...)
    return create_guild_channel(c, g.id; kwargs...)
end
function create(c::Client, ::Type{DiscordChannel}; kwargs...)
    return if haskey(kwargs, :recipient_id)
        create_dm(c; kwargs...)
    else
        create_group_dm(c; kwargs...)
    end
end

retrieve(c::Client, ::Type{DiscordChannel}, channel::Integer) = get_channel(c, channel)
retrieve(c::Client, ::Type{DiscordChannel}, g::AbstractGuild) = get_guild_channels(c, g.id)

update(c::Client, ch::DiscordChannel; kwargs...) = modify_channel(c, ch.id; kwargs...)

delete(c::Client, ch::DiscordChannel) = delete_channel(c, ch.id)


text:
include("typedefs.jl")

function calc_D_MN(providers::Vector{Node}, clients::Vector{Node}, avg_speed::Float64)::Matrix{Int}
    m = length(providers)
    n = length(clients)

    D_MN = zeros(Int, (m, n))

    for k in range(1, stop=m)
        for i in range(1, stop=n)
            p = providers[k]
            c = clients[i]
            d_ki = round(hypot(p.cx-c.cx, p.cy-c.cy) / avg_speed)
            if (d_ki == 0) && (k != i)
                d_ki += 1
            end
            D_MN[k, i] = d_ki
        end
    end

    return D_MN
end

function calc_D_NN(clients::Vector{Node}, avg_speed::Float64)::Matrix{Int}
    n = length(clients)

    D_NN = zeros((n, n))

    for i in range(1, stop=n)
        for j in range(1, stop=n)
            c_outer = clients[i]
            c_inner = clients[j]
            d_ij = round(hypot(c_outer.cx-c_inner.cx, c_outer.cy-c_inner.cy) / avg_speed)
            if (d_ij == 0) && (i != j)
                d_ij += 1
            end
            D_NN[i, j] = d_ij
        end
    end

    return D_NN
end

max_stars_repo_name: dwirya/nodegen

max_issues_repo_name: dwirya/nodegen

max_forks_repo_name: dwirya/nodegen


text:
while (true)
  print("a")
end


text:
## PulseDesigner
module PulseDesigner
using ..Koma
using ..Koma: γ, Scanner, get_bvalue, get_max_grad

###############
## RF Pulses ##
###############
RF_hard(B1, T, sys::Scanner; G=[0,0,0]) = begin
	ζ = sum(G) / sys.Smax
	EX = Sequence([	Grad(G[1],T,ζ);	 #Gx
					Grad(G[2],T,ζ);  #Gy
					Grad(G[3],T,ζ);;], #Gz
					 [RF(B1,T,0,ζ);;]	 #RF
					)
	EX
end
##################
## Gradient SEQ ##
##################
# # Without moment nulling
# DIF_base(G0,Δ,δ;verbose=false) = begin
# 	DIF = Sequence([Grad(G0,δ) Delay(Δ-δ) -Grad(G0,δ); #Gx
# 	         		Grad(0, δ) Delay(Δ-δ)  Grad(0 ,δ)]) #Gy
# 	if verbose
# 	b_value = get_bvalue(DIF)[1]
# 	GDmax = get_max_grad(DIF)
# 	println("Sequence with diffusion gradient: "*string(round(GDmax*1e3,digits=2))*" mT/m") # s/mm2
# 	println("Sequence with b-value: "*string(round(b_value*1e-6,digits=2))*" s/mm2") # s/mm2
# 	end
# 	DIF
# end
# # With moment nulling: M0,M1,M2
# DIF_null(G0,Δ,δ,verbose=false) = begin
# 	if abs(Δ-4δ/3)<1e-6 #Necessary condition
# 	#Stoecker
# 	DIF = Sequence([Grad(G0,δ/3) -Grad(G0,2δ/3) Delay(Δ-δ) Grad(G0,2δ/3) -Grad(G0,δ/3); #Gx
#              		Delay(δ/3) Delay(2δ/3) Delay(Δ-δ) Delay(2δ/3) Delay(δ/3)]) #Gy
# 	else
# 	#Castillo
# 	DIF = Sequence([ Grad(G0,δ/4) -Grad(G0,δ/2) Grad(G0,δ/4) Delay(Δ-δ)
# 					-Grad(G0,δ/4)  Grad(G0,δ/2)-Grad(G0,δ/4); #Gx
#              		 Delay(δ/4) Delay(δ/2) Delay(δ/4) Delay(Δ-δ)
# 					 Delay(δ/4) Delay(δ/2) Delay(δ/4)])
# 	end

# 	if verbose
# 	b_value = get_bvalue(DIF)[1]
# 	GDmax = get_max_grad(DIF)
# 	println("Sequence with diffusion gradient: "*string(round(GDmax*1e3,digits=2))*" mT/m") # s/mm2
# 	println("Sequence with b-value: "*string(round(b_value*1e-6,digits=2))*" s/mm2") # s/mm2
# 	end
# 	DIF
# end
# EPI
EPI(FOV::Float64, N::Int, sys::Scanner) = begin
    #TODO: consider when N is even
	Δt = sys.ADC_Δt
	Gmax = sys.Gmax
	Nx = Ny = N #Square acquisition
	Δx = FOV/(Nx-1)
	Ta = Δt*(Nx-1) #4-8 us
    Δτ = Ta/(Ny-1)
	Ga = 1/(γ*Δt*FOV)
	ζ = Ga / sys.Smax
	Ga ≥ sys.Gmax ? error("Ga=$(Ga*1e3) mT/m exceeds Gmax=$(Gmax*1e3) mT/m, increase Δt to at least Δt_min="
	*string(round(1/(γ*Gmax*FOV),digits=2))*" us.") : 0
	ϵ1 = Δτ/(Δτ+ζ)
	#EPI base
	EPI = Sequence(vcat(
	    [mod(i,2)==0 ? Grad(Ga*(-1)^(i/2),Ta,ζ) : Grad(0.,Δτ,ζ) for i=0:2*Ny-2],  #Gx
	 	[mod(i,2)==1 ? ϵ1*Grad(Ga,Δτ,ζ) :         Grad(0.,Ta,ζ) for i=0:2*Ny-2])) #Gy
	EPI.ADC = [mod(i,2)==1 ? ADC(0,Δτ,ζ) : ADC(N,Ta,ζ) for i=0:2*Ny-2]
	# Relevant parameters
	Δfx_pix = 1/Ta
	Δt_phase = (Ta+2ζ)*Ny + (Δτ+2ζ)*Ny
	Δfx_pix_phase = 1/Δt_phase
	println("## EPI parameters ##")
	println("Δx = $(round(Δx*1e3,digits=2)) mm")
	println("Pixel Δf in freq. direction $(round(Δfx_pix,digits=2)) Hz")
	println("Pixel Δf in phase direction $(round(Δfx_pix_phase,digits=2)) Hz")
	#Pre-wind and wind gradients
	ϵ2 = Ta/(Ta+ζ)
    PHASE =   Sequence(1/2*[Grad(-Ga, Ta, ζ)      ; ϵ2*Grad(-Ga, Ta, ζ);;]) #This needs to be calculated differently
	DEPHASE = Sequence(1/2*[Grad((-1)^N*Ga, Ta, ζ); ϵ2*Grad(-Ga, Ta, ζ);;]) #for even N
	seq = PHASE+EPI+DEPHASE
	#Saving parameters
	seq.DEF = Dict("Nx"=>Nx,"Ny"=>Ny,"Nz"=>1,"Name"=>"epi")
	seq
end

# Radial
radial_base(FOV::Float64, Nr::Int, sys::Scanner) = begin
	Δt = sys.ADC_Δt
	Gmax = sys.Gmax
	Δx = FOV/(Nr-1)
	Ta = Δt*(Nr-1)
	Ga = 1/(γ*Δt*FOV)
	ζ = Ga / sys.Smax
	Ga ≥ sys.Gmax ? error("Ga=$(Ga*1e3) mT/m exceeds Gmax=$(Gmax*1e3) mT/m, increase Δt to at least Δt_min="
	*string(round(1/(γ*Gmax*FOV),digits=2))*" us.") : 0
	#Radial base
	rad = Sequence([Grad(Ga,Ta,ζ)]) #Gx
	rad.ADC = [ADC(Nr,Ta,ζ)]
	# Acq/Recon parameters
	Nspokes = ceil(Int64, π/2 * Nr ) #Nyquist in the radial direction
	Δθ = π / Nspokes
	println("## Radial parameters ##")
	println("FOVr = $(round(FOV*1e2,digits=2)) cm; Δr = $(round(Δx*1e3,digits=2)) mm")
	println("Nspokes = $Nspokes, to satisfy the Nyquist criterion")
    PHASE = Sequence([Grad(-Ga/2, Ta, ζ)])
	seq = PHASE+rad+PHASE
	#Saving parameters
	seq.DEF = Dict("Nx"=>Nr,"Ny"=>Nspokes,"Nz"=>1,"Δθ"=>Δθ,"Name"=>"radial","FOV"=>[FOV, FOV, 0])

	seq
end

export EPI, radial_base
end


text:
# Test TLE checksum
let
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    csum1 = tle_checksum(line1[1:end-1])
    csum2 = tle_checksum(line2[1:end-1])

    @test csum1 == parse(Int, line1[end])
    @test csum2 == parse(Int, line2[end])
end

# Test TLE initialization
let
    # Bad line 1 checksum
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4752"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    @test_throws ArgumentError TLE(line1, line2)

    # Bad line 2 checksum
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413668"

    @test_throws ArgumentError TLE(line1, line2)

    # Bad line 1 length
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  759"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    @test_throws ArgumentError TLE(line1, line2)

    # Bad line 1 length
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.8241915741361"

    @test_throws ArgumentError TLE(line1, line2)
end

let
    # Test Spacecraft Report 3 TLE 1
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    tle = TLE(line1, line2)
    rv  = sgp4(tle, tle.epoch)
    array_isapprox(rv[1:3], [7022.46647249, -1400.06656182, 0.05106558], atol=1e-8)
    # array_isapprox(rv[4:6], [1.893831081, 6.405894873, 4.534806701], atol=1e-8)

    # Test Spacecraft Report 3 TLE 2
    line1 = "1 04632U 70093B   04031.91070959 -.00000084  00000-0  10000-3 0  9955"
    line2 = "2 04632  11.4628 273.1101 1450506 207.6000 143.9350  1.20231981 44145"

    tle = TLE(line1, line2)
    rv  = sgp4(tle, tle.epoch)
    # array_isapprox(rv[1:3], [2334.12240168, -41920.42940755, -0.03695311], atol=1e-8)
    # array_isapprox(rv[4:6], [2.826320327, -0.065091320, 0.570935915], atol=1e-8)

    # Test Spacecraft Report 3 TLE 3
    line1 = "1 06251U 62025E   06176.82412014  .00008885  00000-0  12808-3 0  3985"
    line2 = "2 06251  58.0579  54.0425 0030035 139.1568 221.1854 15.56387291  6774"

    tle = TLE(line1, line2)
    rv  = sgp4(tle, tle.epoch)
    array_isapprox(rv[1:3], [3988.29488121, 5498.97562026, 0.92898335], atol=1e-8)
    # array_isapprox(rv[4:6], [-3.290042995, 2.357636871, 6.496621836], atol=1e-8)
end

let
    # Test Spacecraft Report 3 TLE 1
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    tle = TLE(line1, line2)
    rv  = ecef(tle, tle.epoch)
end

let
    # Test Spacecraft Report 3 TLE 1
    line1 = "1 00005U 58002B   00179.78495062  .00000023  00000-0  28098-4 0  4753"
    line2 = "2 00005  34.2682 348.7242 1859667 331.7664  19.3264 10.82419157413667"

    tle = TLE(line1, line2)
    rv  = eci(tle, tle.epoch)
end


text:
mutable struct Multi
    lock   :: ReentrantLock
    handle :: Ptr{Cvoid}
    timer  :: Ptr{Cvoid}
    easies :: Vector{Easy}
    grace  :: UInt64

    function Multi(grace::Integer = typemax(UInt64))
        timer = jl_malloc(Base._sizeof_uv_timer)
        uv_timer_init(timer)
        multi = new(ReentrantLock(), C_NULL, timer, Easy[], grace)
        finalizer(multi) do multi
            uv_timer_stop(multi.timer)
            uv_close(multi.timer, cglobal(:jl_free))
            done!(multi)
        end
    end
end

function init!(multi::Multi)
    multi.handle != C_NULL && return
    multi.handle = curl_multi_init()
    add_callbacks(multi)
    set_defaults(multi)
end

function done!(multi::Multi)
    multi.handle == C_NULL && return
    curl_multi_cleanup(multi.handle)
    multi.handle = C_NULL
end

# adding & removing easy handles

function cleanup_callback(uv_timer_p::Ptr{Cvoid})::Cvoid
    ## TODO: use a member access API
    multi_p = unsafe_load(convert(Ptr{Ptr{Cvoid}}, uv_timer_p))
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    done!(multi)
    return
end

function add_handle(multi::Multi, easy::Easy)
    lock(multi.lock) do
        if isempty(multi.easies)
            preserve_handle(multi)
            uv_timer_stop(multi.timer) # stop grace timer
        end
        push!(multi.easies, easy)
        init!(multi)
        @check curl_multi_add_handle(multi.handle, easy.handle)
    end
end

function remove_handle(multi::Multi, easy::Easy)
    lock(multi.lock) do
        @check curl_multi_remove_handle(multi.handle, easy.handle)
        deleteat!(multi.easies, findlast(==(easy), multi.easies))
        !isempty(multi.easies) && return
        cleanup_cb = @cfunction(cleanup_callback, Cvoid, (Ptr{Cvoid},))
        if multi.grace <= 0
            done!(multi)
        elseif 0 < multi.grace < typemax(multi.grace)
            uv_timer_start(multi.timer, cleanup_cb, multi.grace, 0)
        end
        unpreserve_handle(multi)
    end
end

# multi-socket options

function set_defaults(multi::Multi)
    # currently no defaults
end

# libuv callbacks

struct CURLMsg
   msg  :: CURLMSG
   easy :: Ptr{Cvoid}
   code :: CURLcode
end

# should already be locked
function check_multi_info(multi::Multi)
    while true
        p = curl_multi_info_read(multi.handle, Ref{Cint}())
        p == C_NULL && return
        message = unsafe_load(convert(Ptr{CURLMsg}, p))
        if message.msg == CURLMSG_DONE
            easy_handle = message.easy
            easy_p_ref = Ref{Ptr{Cvoid}}()
            @check curl_easy_getinfo(easy_handle, CURLINFO_PRIVATE, easy_p_ref)
            easy = unsafe_pointer_to_objref(easy_p_ref[])::Easy
            @assert easy_handle == easy.handle
            easy.code = message.code
            close(easy.progress)
            close(easy.output)
            easy.input = nothing
            notify(easy.ready)
        else
            @async @error("curl_multi_info_read: unknown message", message)
        end
    end
end

function event_callback(
    uv_poll_p :: Ptr{Cvoid},
    status    :: Cint,
    events    :: Cint,
)::Cvoid
    ## TODO: use a member access API
    multi_p = unsafe_load(convert(Ptr{Ptr{Cvoid}}, uv_poll_p))
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    sock_p = uv_poll_p + Base._sizeof_uv_poll
    sock = unsafe_load(convert(Ptr{curl_socket_t}, sock_p))
    flags = 0
    events & UV_READABLE != 0 && (flags |= CURL_CSELECT_IN)
    events & UV_WRITABLE != 0 && (flags |= CURL_CSELECT_OUT)
    lock(multi.lock) do
        @check curl_multi_socket_action(multi.handle, sock, flags)
        check_multi_info(multi)
    end
end

function timeout_callback(uv_timer_p::Ptr{Cvoid})::Cvoid
    ## TODO: use a member access API
    multi_p = unsafe_load(convert(Ptr{Ptr{Cvoid}}, uv_timer_p))
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    lock(multi.lock) do
        @check curl_multi_socket_action(multi.handle, CURL_SOCKET_TIMEOUT, 0)
        check_multi_info(multi)
    end
end

# curl callbacks

function timer_callback(
    multi_h    :: Ptr{Cvoid},
    timeout_ms :: Clong,
    multi_p    :: Ptr{Cvoid},
)::Cint
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    @assert multi_h == multi.handle
    if timeout_ms == 0
        lock(multi.lock) do
            @check curl_multi_socket_action(multi.handle, CURL_SOCKET_TIMEOUT, 0)
            check_multi_info(multi)
        end
    elseif timeout_ms >= 0
        timeout_cb = @cfunction(timeout_callback, Cvoid, (Ptr{Cvoid},))
        uv_timer_start(multi.timer, timeout_cb, max(1, timeout_ms), 0)
    elseif timeout_ms == -1
        uv_timer_stop(multi.timer)
    else
        @async @error("timer_callback: invalid timeout value", timeout_ms)
        return -1
    end
    return 0
end

function socket_callback(
    easy_h    :: Ptr{Cvoid},
    sock      :: curl_socket_t,
    action    :: Cint,
    multi_p   :: Ptr{Cvoid},
    uv_poll_p :: Ptr{Cvoid},
)::Cint
    multi = unsafe_pointer_to_objref(multi_p)::Multi
    if action in (CURL_POLL_IN, CURL_POLL_OUT, CURL_POLL_INOUT)
        if uv_poll_p == C_NULL
            uv_poll_p = uv_poll_alloc()
            uv_poll_init(uv_poll_p, sock)
            ## TODO: use a member access API
            unsafe_store!(convert(Ptr{Ptr{Cvoid}}, uv_poll_p), multi_p)
            sock_p = uv_poll_p + Base._sizeof_uv_poll
            unsafe_store!(convert(Ptr{curl_socket_t}, sock_p), sock)
            lock(multi.lock) do
                @check curl_multi_assign(multi.handle, sock, uv_poll_p)
            end
        end
        events = 0
        action != CURL_POLL_IN  && (events |= UV_WRITABLE)
        action != CURL_POLL_OUT && (events |= UV_READABLE)
        event_cb = @cfunction(event_callback, Cvoid, (Ptr{Cvoid}, Cint, Cint))
        uv_poll_start(uv_poll_p, events, event_cb)
    elseif action == CURL_POLL_REMOVE
        if uv_poll_p != C_NULL
            uv_poll_stop(uv_poll_p)
            uv_close(uv_poll_p, cglobal(:jl_free))
            lock(multi.lock) do
                @check curl_multi_assign(multi.handle, sock, C_NULL)
            end
        end
    else
        @async @error("socket_callback: unexpected action", action)
        return -1
    end
    return 0
end

function add_callbacks(multi::Multi)
    # stash multi handle pointer in timer
    multi_p = pointer_from_objref(multi)
    ## TODO: use a member access API
    unsafe_store!(convert(Ptr{Ptr{Cvoid}}, multi.timer), multi_p)

    # set timer callback
    timer_cb = @cfunction(timer_callback, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
    setopt(multi, CURLMOPT_TIMERFUNCTION, timer_cb)
    setopt(multi, CURLMOPT_TIMERDATA, multi_p)

    # set socket callback
    socket_cb = @cfunction(socket_callback,
        Cint, (Ptr{Cvoid}, curl_socket_t, Cint, Ptr{Cvoid}, Ptr{Cvoid}))
    setopt(multi, CURLMOPT_SOCKETFUNCTION, socket_cb)
    setopt(multi, CURLMOPT_SOCKETDATA, multi_p)
end


text:
using Interpolations
using kRPC
using CSV
using LinearAlgebra
using DataFrames
conn = kRPCConnect("Aerodata")
using kRPC.Remote.SpaceCenter
totable = 1
outp_a = 1
outp_b = 1

vessel = ActiveVessel()
flight = Flight(vessel)
orbit = Orbit(vessel)
body = Body(orbit)
body_frame = ReferenceFrame(body)
control = Control(vessel)

function compute_reference_frame(body, lat, lon, alt::Float64)
    landing_position = SurfacePosition(body, lat, lon, body_frame)
    q_long = (
        0.0,
        sin(-lon * 0.5 * pi / 180),
        0.0,
        cos(-lon * 0.5 * pi / 180))
    q_lat = (
        0.0,
        0.0,
        sin(lat * 0.5 * pi / 180),
        cos(lat * 0.5 * pi / 180))
    cr1 = CreateRelative(body_frame, position=landing_position, rotation=q_long)
    cr2 = CreateRelative(cr1,rotation=q_lat)
    return CreateRelative(cr2,position=(alt,0.0,0.0))
end
guidance_frame = compute_reference_frame(body, Latitude(flight), Longitude(flight), 0.0)

guidance_flight = Flight(vessel)#, referenceFrame=guidance_frame)
direction = Float64[-1,0,0]
speed_of_sound = SpeedOfSound(flight)

function compute_aero_forces(mach, cos_aoa)
	velocity = mach*speed_of_sound
	velocity_vec = broadcast(*, Float64[cos_aoa,sqrt(1-cos_aoa^2),0], velocity)
	return kRPC.Remote.SpaceCenter.Delayed.SimulateAerodynamicForceAt(guidance_flight, body, (0.0,0.0,0.0), ((velocity_vec)...,)),
		   kRPC.Remote.SpaceCenter.Delayed.SimulateAerodynamicTorqueAt(guidance_flight, body, (0.0,0.0,0.0), ((velocity_vec)...,)), velocity_vec
end

function sweep_aero()
	cforce = kRPC.kPC[]
	ctorque = kRPC.kPC[]
	drg = Vector{Float64}[]
	lft = Vector{Float64}[]
	trq = Vector{Float64}[]
	aoal = Float64[]
	machl = Float64[]
	for mach=0.0:0.025:1.5
		for cos_aoa=cosd(180):1/90:cosd(0)
			cll,ctr,bv = compute_aero_forces(mach, cos_aoa)
			push!(cforce, cll)
			push!(ctorque, ctr)
			normbv = bv/norm(bv)
			push!(drg, normbv)
			torqd = cross(normbv, [1.0,0.0,0.0])
			liftd = cross(torqd,normbv)
			liftd = liftd/norm(liftd)
			push!(lft, liftd)
			push!(trq, torqd/norm(torqd))
			push!(aoal, cos_aoa)
			push!(machl, mach)
		end
	end
	res = kRPC.SendMessage(conn, cforce)
	tres = kRPC.SendMessage(conn, ctorque)
	hcat(([aoal[i], machl[i],
			if !isnan(dot(res[i], drg[i])) dot(res[i], drg[i]) else 0.0 end,
			if !isnan(dot(res[i], lft[i])) dot(res[i], lft[i]) else 0.0 end,
			if !isnan(dot(tres[i], trq[i])) dot(tres[i], trq[i]) else 0.0 end] for i=1:length(res))...)
end

function compute_aero_table()
	aero_force = sweep_aero()'
	aero_df = DataFrame(aoa=aero_force[:,1], mach=aero_force[:,2], drag=aero_force[:,3], lift=aero_force[:,4], torque=aero_force[:,5])
	CSV.write("lift_drag.csv", aero_df);
end

pts = Parts(vessel)
finsOn = WithTag(pts, "nosefins2")
fin1 = Modules(finsOn[1])[1]
fin2 = Modules(finsOn[2])[1]
finsOff = WithTag(pts, "nosefins")
fin3 = Modules(finsOff[1])[1]
fin4 = Modules(finsOff[2])[1]

function compute_fin_force_table()
	SetFieldFloat.([fin3,fin4], "Flp/Splr Dflct", convert(Float32, 0))
	SetFieldFloat.([fin1,fin2], "Flp/Splr Dflct", convert(Float32, 0))
	sleep(4)

	function aero_profile()
		return map(x -> [x...], kRPC.SendMessage(conn, [begin frcc,trqc,vv = compute_aero_forces(mach, 1.0); frcc end for mach=0.01:0.025:1.5]))
	end
	bl = aero_profile()
	forces = Vector{Vector{Float64}}[]
	for ang=0.0:0.1:90.0
		SetFieldFloat.([fin1,fin2], "Flp/Splr Dflct", convert(Float32, ang))
		sleep(0.1)
		push!(forces, aero_profile())
	end
	aerod = hcat(map((force,aoa) -> vcat(hcat((force .- bl)...), transpose(collect(0.01:0.025:1.5)), fill(aoa, 1, length(bl))), forces, 0.0:0.1:90.0)...)
	df = DataFrame(lift=aerod[2,:], drag=aerod[1,:], mach=aerod[4,:], aoa=aerod[5,:])
	CSV.write("fin.csv", df);
end

#=

function compute_aero_forces_tform(bv,vv)
	if (vv>0)
		alpha = acos(bv)
	else
		alpha = 0 #zero velocity anyways
	end
	mach = sqrt(vv)
	return compute_aero_forces(mach, alpha)
end

ctrls = kRPC.Remote.FerramControlSurface.FlightControls(vessel)
map(x->kRPC.Remote.FerramControlSurface.ControlMaxDeflection!(x, 80.0f0), ctrls)

function compute_control_force_defl(mach, aoamin, aoamax, minangle, maxangle, func)
	kRPC.Remote.FerramControlSurface.ControlTimeConstant!(0.00000001)
	func(control, minangle)
	sleep(0.1)

	forces = []
	for aoa in aoamin:pi/64:aoamax
		func(control, minangle)
		sleep(0.2)
		basf = kRPC.SendMessage(conn, compute_aero_forces(mach, aoa)[1])
		for angle in minangle:0.1f0:maxangle
			func(control, angle)
			println(angle)
			sleep(0.1)
			frce = [[(basf .- kRPC.SendMessage(conn, compute_aero_forces(mach, aoa)[1]))...] for i=1:10]
			mfrce = mean(frce)
			unshift!(mfrce, aoa)
			push!(forces, mfrce)
		end
	end
	func(control, 0.0f0)
	sleep(0.1)
	return hcat(forces...)
end

function compute_control_force_drag(mach, angle, delta)
	start_a = angle - delta/2
	end_a = angle + delta/2
	stepsize = delta/40
	force = []
	kRPC.Remote.FerramControlSurface.ControlTimeConstant!(0.00000001)
	Pitch!(control, convert(Float32, angle/deg2rad(80.0)))
	sleep(0.3)
	basef = kRPC.SendMessage(conn, compute_aero_forces(mach, pi+angle)[1])
	Pitch!(control, convert(Float32, start_a/deg2rad(80.0)))
	sleep(0.3)
	for ang in start_a:stepsize:end_a
		Pitch!(control, convert(Float32, ang/deg2rad(80.0)))
		sleep(0.05)
		push!(force, [(kRPC.SendMessage(conn, compute_aero_forces(mach, pi+angle)[1]) .- basef)..., ang])
	end
	return basef,force
end

function find_apex(mach, angle, delta)
	start_a = angle - delta/2
	end_a = angle + delta/2
	stepsize = delta/10
	angles = Any[]
	for ang in start_a:stepsize:end_a
		bf,fr = compute_control_force_drag(mach, ang, delta)
		fr = hcat(fr...)
		sz = size(fr)[2]
		which_angle = indmin(abs.((fr[1,1:sz-1] - fr[1,2:sz])./(fr[2,1:sz-1] - fr[2,2:sz])))
		push!(angles, [bf,fr])
	end
	return angles
end

function fit_parabola(t,b)
	inp = hcat(t.^2,t.^4)
	println(inp, " ", b)
	return inp\b
end

#=
bv_step = 0.005
vv_step = 0.1
max_mach = 2
reqs = [(bv,vv,compute_aero_forces_tform(bv,vv))
             for vv in 0:vv_step:max_mach^2
             for bv in 1:-bv_step:-1]
tosend = map(x->x[3][1], reqs)
res = kRPC.SendMessage(conn, tosend)

density = AtmosphereDensity(flight)
function decompose(aero, velocity, bv, vv)
	iv = velocity
	ia = aero
	velocity = velocity[1:2] # hey this is really a 2d problem
	aero = aero[1:2]/speed_of_sound
	velocity_dir = velocity/norm(velocity)
	drag = dot(aero,velocity_dir)
	lift = dot([-velocity_dir[2],velocity_dir[1]], aero - drag*velocity_dir)

	bx = bv
	by = sqrt(1-bx^2)

	if abs(by) < 0.0001
		res= (0.0, drag/sqrt(vv))
	else
		res= (lift/by, (drag-bx*lift/by)/sqrt(vv))
	end
	if isnan(res[1]) res = (0.0, res[2]) end
	if isnan(res[2]) res = (res[1], 0.0) end
	return res
end
totable = map((inp,res)->(inp[1],inp[2])=>decompose([res...], inp[3][2], inp[1], inp[2]), reqs, res)

bvs = unique(map(x->x[1][1], totable))
machsq = unique(map(x->x[1][2], totable))
map_bv = Dict([bvs[i] => i for i = 1:length(bvs)])
map_machsq = Dict([machsq[i] => i for i = 1:length(machsq)])
outp_a = Array{Float64,2}(length(bvs), length(machsq))
outp_b = Array{Float64,2}(length(bvs), length(machsq))
for result in totable
	outp_a[map_bv[result[1][1]], map_machsq[result[1][2]]] = result[2][1]
	outp_b[map_bv[result[1][1]], map_machsq[result[1][2]]] = result[2][2]
end
=#
=#

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


text:
"""
    All learning algorithms assume a single value function!!!
"""

using Flux
import LinearAlgebra: dot
import Statistics: mean

_square(x::AbstractArray) = x.*x
_prod(x::Array{T, 1}, y::Array{T, 1}) where {T<:Number} = x.*y
_prod(x::Array{T, 1}, y::AbstractArray) where {T<:Number} = x.*y

tderror(v_t::Array{<:AbstractFloat, 1}, c::Array{<:AbstractFloat, 1}, γ_tp1::Array{<:AbstractFloat, 1}, ṽ_tp1::Array{<:AbstractFloat, 1}) where {A<:AbstractArray} = v_t .- (c .+ γ_tp1.*ṽ_tp1)
tderror(v_t, c, γ_tp1, ṽ_tp1) = v_t - (c + γ_tp1*ṽ_tp1)

function offpolicy_tdloss(ρ_t::Array{Array{T, 1}, 1},
                          v_t::AbstractArray,
                          c::Array{Array{T, 1}, 1},
                          γ_tp1::Array{Array{T,1}, 1},
                          ṽ_tp1::Array{Array{T,1}, 1}) where {T<:AbstractFloat}
    target = c .+ _prod.(γ_tp1, ṽ_tp1)
    error = _square.(v_t .- target)
    return (T(0.5))*sum(mean.(_prod.(ρ_t, error))) * (1//length(ρ_t))
end

tdloss(v_t, c, γ_tp1, ṽ_tp1) =
    0.5*Flux.mse(v_t, Flux.data(c .+ γ_tp1.*ṽ_tp1))

"""
    LearningUpdate
WIP - Currently LearningUpdate and Optimizer are haphazardly similar....
"""
abstract type LearningUpdate end

update!(model, opt, lu::LearningUpdate, ρ, s_t, s_tp1, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0) =
    update!(model, opt, lu::LearningUpdate, ρ, s_t, s_tp1, r, γ, terminal)

abstract type AbstractTDUpdate <: LearningUpdate end

mutable struct BatchTD <: AbstractTDUpdate end

function update!(model, opt, lu::BatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0f0)
    v_t = model.(s_t)
    v_tp1 = model.(s_tp1)
    loss = offpolicy_tdloss(ρ.*corr_term, v_t, r, γ, Flux.data.(v_tp1))
    grads = Flux.gradient(()->loss, params(model))
    for weights in params(model)
        Flux.Optimise.update!(opt, weights, grads[weights])
    end
end

function update!(model::SingleLayer, opt, lu::BatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    v_t = model.(s_t)
    v_tp1 = model.(s_tp1)
    dvdt = deriv.(model, s_t)
    δ = ρ.*tderror(v_t, r, γ, v_tp1)
    Δ = mean(δ.*dvdt.*s_t)
    model.W .+= -Flux.Optimise.apply!(opt, model.W, corr_term*Δ)
end

function update!(model::SparseLayer, opt::Descent, lu::BatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    v_t = model.(s_t)
    v_tp1 = model.(s_tp1)
    dvdt = [deriv(model, s) for s in s_t]
    δ = ρ.*tderror(v_t, r, γ, v_tp1)
    Δ = δ.*dvdt.*1//length(ρ)

    for i in 1:length(ρ)
        model.W[s_t[i]] .-= opt.eta*corr_term*Δ[i]
    end
end

function update!(model::SparseLayer, opt::RMSProp, lu::BatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    v_t = model.(s_t)
    v_tp1 = model.(s_tp1)

    dvdt = [deriv(model, s) for s in s_t]
    δ = ρ.*tderror(v_t, r, γ, v_tp1)
    fill!(model.ϕ, 0.0f0)
    Δ = δ.*dvdt.*1//length(ρ)
    for i in 1:length(ρ)
        model.ϕ[s_t[i]] .+= corr_term*Δ[i]
    end
    feats = unique(collect(Iterators.flatten(s_t)))

    acc = get!(opt.acc, model.W, zero(model.W))::typeof(model.W)
    acc .*= opt.rho
    acc[feats] .+= (1-opt.rho) .* model.ϕ[feats].^2
    model.W[feats] .-= model.ϕ[feats].*(opt.eta./sqrt.(acc[feats] .+ 1e-8))

end

function update!(model::TabularLayer, opt::Descent, lu::BatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    # println(model, s_t)
    v_t = model.(s_t)
    v_tp1 = model.(s_tp1)
    δ = ρ.*tderror(v_t, r, γ, v_tp1)
    Δ = corr_term.*δ.*(1.0/length(v_t))
    for (s_idx, s) in enumerate(s_t)
        model.W[s...] -= opt.eta*Δ[s_idx]
    end
end

struct IncTD <: AbstractTDUpdate end

function update!(model::SparseLayer,
                 opt::Descent,
                 lu::IncTD,
                 ρ, s_t, s_tp1,
                 r, γ, terminal;
                 corr_term=1.0)

    for i in 1:length(ρ)
        v_t = model(s_t[i])
        v_tp1 = model(s_tp1[i])
        dvdt = deriv(model, s_t[i])
        δ = ρ[i]*tderror(v_t, r[i], γ[i], v_tp1)
        model.W[s_t[i]] .-= opt.eta*corr_term*δ*dvdt
    end

end

mutable struct WISBatchTD{T<:AbstractTDUpdate} <: LearningUpdate
    batch_td::T
end

function update!(model, opt, lu::WISBatchTD, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    wis_sum = sum(ρ) + 1e-8
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=(1.0/wis_sum)*corr_term*length(ρ))
end

function update!(model, opt, lu::WISBatchTD, ρ::Array{T,1}, s_t, s_tp1, r, γ, terminal; corr_term=1.0) where {T<:AbstractArray}
    wis_sum = (sum(ρ) .+ 1e-8)
    # println(typeof(wis_sum))
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=T(length(ρ)*corr_term./wis_sum))
end

mutable struct WISBatchTD_Rupam <: LearningUpdate
    z::IdDict
    u::IdDict
    d::IdDict
    v::IdDict
    prev_model::IdDict
    u_0::Float64
    normalize_eta::Bool
    WISBatchTD_Rupam(u_0, normalize_eta=false) = new(IdDict(), IdDict(), IdDict(), IdDict(), IdDict(), u_0, normalize_eta)
end

_many_hot(type::Type, size, x::Array{Int, 1}) =
    begin; t = zeros(type, size); t[x] .= one(type); return t; end;

_many_hot(size, x) = _many_hot(Float64, size, x)

function update!(model::SparseLayer, opt::Descent, lu::WISBatchTD_Rupam, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)

    # println(s_t)
    # z_vec = get!(lu.z, model, zeros(length(model.W)))
    # u_vec = get!(lu.u, model, zeros(length(model.W)))
    d_vec = get!(lu.d, model, lu.u_0.*length(s_t[1]).*ones(length(model.W)))::Array{Float64, 1}
    # v_vec = get!(lu.v, model, zeros(length(model.W)))

    # println("Step Update")
    for i in 1:length(ρ)
        ϕ = _many_hot(length(model.W), s_t[i])
        # @show ϕ
        ϕϕ = ϕ
        ϕnext = _many_hot(length(model.W), s_tp1[i])
        R = r[i]
        g = 0
        λ = 0
        λnext = 0
        gnext = γ[i]
        rho = ρ[i]
        η = opt.eta
        if lu.normalize_eta
            η = opt.eta/(lu.u_0*length(s_t[1]))
        end

        d_vec .= d_vec - η.*(ϕϕ).*d_vec + rho.*ϕϕ
        dtemp = copy(d_vec)
        dtemp[dtemp.==0.0] .= 1
        alpha = 1 ./ dtemp
        αϕ = alpha[s_t[i]].*ϕ[s_t[i]]
        prednext = model(s_tp1[i])

        v_old = 0.0
        if model ∈ keys(lu.prev_model)
            v_old = lu.prev_model[model](s_t[i])
        end

        model.W[1,s_t[i]] .= model.W[1,s_t[i]] .+ rho.*((R + gnext*prednext - v_old) + (v_old - model(s_t[i]))).*αϕ

        if model ∈ keys(lu.prev_model)
            lu.prev_model[model].W .= model.W
        else
            lu.prev_model[model] = deepcopy(model)
        end
    end

    # @show sum(model.W)

end

function update!(model::TabularLayer, opt::Descent, lu::WISBatchTD_Rupam, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)

    # println(s_t)
    z_vec = get!(lu.z, model, zeros(length(model.W)))
    u_vec = get!(lu.u, model, zeros(length(model.W)))
    d_vec = get!(lu.d, model, lu.u_0.*ones(length(model.W)))
    v_vec = get!(lu.v, model, zeros(length(model.W)))

    # println("Step Update")
    for i in 1:length(ρ)
        ϕ = _many_hot(length(model.W), [s_t[i]])
        # @show ϕ
        ϕϕ = ϕ.*ϕ
        ϕnext = _many_hot(length(model.W), [s_tp1[i]])
        R = r[i]
        g = 0
        λ = 0
        λnext = 0
        gnext = γ[i]
        rho = ρ[i]
        η = opt.eta

        d_vec .= d_vec - η.*(ϕϕ).*d_vec + rho*ϕϕ
        dtemp = copy(d_vec)
        dtemp[dtemp.==0.0] .= 1
        alpha = 1 ./ dtemp
        αϕ = alpha.*ϕ
        # println(minimum(αϕ), " ", maximum(αϕ))
        v_vec .= rho.*ϕϕ
        z_vec .= rho.*αϕ
        prednext = model(s_tp1[i])

        v_old = 0.0
        if model ∈ keys(lu.prev_model)
            v_old = lu.prev_model[model](s_t[i])
        end
        # @show (rho, gnext, R, η)
        # @show maximum(αϕ)
        # @show maximum((R + gnext*prednext - v_old)*z_vec + rho*(v_old - model(s_t[i]))*αϕ)

        # @show (R + gnext*prednext - v_old)
        # @show maximum(rho*(v_old - model(s_t[i]))*αϕ)
        # @show size(model.W)
        res_W = reshape(model.W, length(model.W))

        model.W .= model.W + ((R + gnext*prednext - v_old)*z_vec + rho*(v_old - model(s_t[i]))*αϕ)
        lu.prev_model[model] = deepcopy(model)
    end

    # @show sum(model.W)

end

# function update!(model::SparseLayer, opt::Descent, lu::WISBatchTD_Rupam, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
#     v_t = model.(s_t)
#     v_tp1 = model.(s_tp1)
#     dvdt = [deriv(model, s) for s in s_t]
#     # δ = ρ.*tderror(v_t, r, γ, v_tp1)
#     # Δ = δ.*dvdt.*1//length(ρ)

#     u_1 = zeros(length(lu.u_vec))
#     u_2 = zeros(length(lu.u_vec))
#     for i in 1:length(ρ)
#         fill!(u_1, 0.0)
#         fill!(u_2, 0.0)

#         v_old = zero(v_t)
#         if model ∈ keys(lu.prev_model)
#             v_old = lu.prev_model[model].(s_t)
#         end

#         δ = ρ.*(tderror(v_old, r, γ, v_tp1) + v_old - v_t)
#         Δ = -δ
#         # u_1[s_t[i]] .+= 1 - opt.eta
#         u_1[s_t[i]] .+= opt.eta
#         u_2[s_t[i]] .+= ρ[i]

#         lu.u_vec .= (1 .- u_1).*lu.u_vec .+ u_2
#         # lu.u_vec .= (1 .- u_1.*(1//length(ρ))).*lu.u_vec .+ u_2.*(1//length(ρ))

#         # println(minimum(lu.u_vec[s_t[i]]), " ", maximum(lu.u_vec[s_t[i]]))
#         # println(minimum(1 ./ (lu.u_vec[s_t[i]])), " ", maximum(1 ./ (lu.u_vec[s_t[i]])))

#         # @show v_t
#         # @show r
#         # for i in 1:length(ρ)
#             # println(length(corr_term*Δ[i] ./ lu.u_vec[s_t[i]]))
#         model.W[s_t[i]] .+= corr_term*Δ[i] ./ (lu.u_vec[s_t[i]] .+ 1e-8)
#         # println(Δ[i])
#         # println(minimum(Δ[i] ./ lu.u_vec[s_t[i]]), " ", maximum(Δ[i] ./ lu.u_vec[s_t[i]]))
#         # end
#         lu.prev_model[model] = deepcopy(model)
#     end
#     # lu.u_vec .= (u_1.*(1//length(ρ))).*lu.u_vec .+ u_2.*(1//length(ρ))

#     println(minimum(model.W), " ", maximum(model.W))

# end

mutable struct VTrace{T<:AbstractTDUpdate} <: LearningUpdate
    ρ_bar::Float64
    batch_td::T
end

# VTrace(ρ_bar, tdupdate) = new(ρ_bar, tdupdate)

function update!(model, opt, lu::VTrace, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    update!(model, opt, lu.batch_td, clamp.(ρ, 0.0, lu.ρ_bar), s_t, s_tp1, r, γ, terminal; corr_term=corr_term)
end

function update!(model, opt, lu::VTrace, ρ::Array{Array{T, 1}, 1}, s_t, s_tp1, r, γ, terminal; corr_term=1.0)  where {T<:AbstractFloat}
    clamp_ρ = [T.(clamp.(_ρ, T(0.0f0), T(lu.ρ_bar))) for _ρ in ρ]
    update!(model, opt, lu.batch_td, clamp_ρ, s_t, s_tp1, r, γ, terminal; corr_term=T(corr_term))
end

mutable struct IncNormIS{T<:AbstractTDUpdate} <: LearningUpdate
    max_is::IdDict
    batch_td::T
end

IncNormIS(tdupdate) = IncNormIS(IdDict(), tdupdate)

function update!(model, opt, lu::IncNormIS, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    if !(model ∈ keys(lu.max_is))
        lu.max_is[model] = 1.0f0
    end
    max_is = max(lu.max_is[model], maximum(ρ))
    lu.max_is[model] = max_is
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term/max_is)
end

function update!(model, opt, lu::IncNormIS, ρ::Array{Array{T, 1}, 1}, s_t, s_tp1, r, γ, terminal; corr_term=1.0f0)  where {T<:AbstractFloat}
    if !(model ∈ keys(lu.max_is))
        lu.max_is[model] = ones(T, length(ρ[1]))
    end
    max_is = lu.max_is[model]
    for idx in 1:length(max_is)
        max_is[idx] = max(max_is[idx], maximum(getindex.(ρ, idx)))
    end
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term./max_is)
end

mutable struct WSNormIS{T<:AbstractTDUpdate} <: LearningUpdate
    beta::Float32
    weighted_sum_is::IdDict
    batch_td::T
end

WSNormIS(beta::Float32, tdupdate) = WSNormIS(beta, IdDict(), tdupdate)

function update!(model, opt, lu::WSNormIS, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    get!(lu.weighted_sum_is, model, 0.0f0)::Float32
    lu.weighted_sum_is[model] = Float32(lu.beta*maximum(ρ)^2 + (1.0f0-lu.beta)*lu.weighted_sum_is[model])
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term/sqrt(lu.weighted_sum_is[model]))
end

function update!(model, opt, lu::WSNormIS, ρ::Array{Array{T, 1}, 1}, s_t, s_tp1, r, γ, terminal; corr_term=1.0f0)  where {T<:AbstractFloat}
    weighted_sum_is = get!(lu.weighted_sum_is, model, zeros(T, length(ρ[1])))
    for idx in 1:length(weighted_sum_is)
        weighted_sum_is[idx] = lu.beta*maximum(getindex.(ρ, idx))^2 + (1.0f0-lu.beta)*weighted_sum_is[idx]
    end
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term./sqrt.(weighted_sum_is))
end

mutable struct WSAvgNormIS{TD<:AbstractTDUpdate} <: LearningUpdate
    beta::Float32
    weighted_sum_is::IdDict
    batch_td::TD
end

WSAvgNormIS(beta::Float32, tdupdate) = WSAvgNormIS(beta, IdDict(), BatchTD())

function update!(model, opt, lu::WSAvgNormIS, ρ, s_t, s_tp1, r, γ, terminal; corr_term=1.0)
    get!(lu.weighted_sum_is, model, 0.0f0)::Float32
    lu.weighted_sum_is[model] = Float32(lu.beta*mean(ρ)^2 + (1.0f0-lu.beta)*lu.weighted_sum_is[model])
    update!(model, opt, lu.batch_td, ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term/sqrt(lu.weighted_sum_is[model]))
end

function update!(model, opt, lu::WSAvgNormIS, ρ::Array{Array{T, 1}, 1}, s_t, s_tp1, r, γ, terminal; corr_term=1.0f0)  where {T<:AbstractFloat}
    weighted_sum_is = get!(lu.weighted_sum_is, model, zeros(T, length(ρ[1])))::type(Array{T, 1})
    for idx in 1:length(max_is)
        weighted_sum_is[idx] = lu.beta*mean(getindex.(ρ, idx))^2 + (1.0f0-lu.beta)*weighted_sum_is[idx]
    end
    update!(model, opt, lu.batch_td, clamp_ρ, s_t, s_tp1, r, γ, terminal; corr_term=corr_term./sqrt.(weighted_sum_is))
end

"""
    All Q Learning algorithms are defined with multiple outputs, there is no cont action cases yet.
"""

# Q Learning

TabularQLayer(num_actions, dims...) = TabularLayer(num_actions, dims...)

mutable struct BatchSarsa <: LearningUpdate end

function update!(model, opt, lu::BatchSarsa, ρ, s_t, s_tp1, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)
    preds_t = [model(s_t[idx])[act] for (idx, act) in enumerate(a_t)]
    preds_tp1 = [model(s_tp1[idx])[act] for (idx, act) in enumerate(a_tp1)]
    grads = Flux.gradient(()->offpolicy_tdloss(preds, r, γ, Flux.data.(preds_tp1)), params(model))
    for weights in params(model)
        update!(opt, weights, -grads[weights])
    end
end

function update!(model::SingleLayer, opt, lu::BatchSarsa, ρ, s_t, s_tp1, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)
    q_t = [model(s_t[idx])[act] for (idx, act) in enumerate(action_t)]
    q_tp1 = [model(s_tp1[idx])[act] for (idx, act) in enumerate(a_tp1)]
    dvdt = [deriv(model, s_t[idx]) for (idx, act) in enumerate(action_t)]
    δ = tderror(v_t, r, γ, v_tp1)
    Δ = mean(δ.*dvdt.*s_t)
    model.W .+= -apply!(opt, model.W, corr_term*Δ)
end

function update!(model::TabularLayer, opt::Descent, lu::BatchSarsa, ρ, s_t::Array{Array{Int64, 1}, 1}, s_tp1::Array{Array{Int64, 1}, 1}, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)

    cis_t = [CartesianIndex(act, s_t[idx]...) for (idx, act) in enumerate(a_t)]
    q_t = model.(cis_t)
    q_tp1 = [model(CartesianIndex(act, s_tp1[idx]...)) for (idx, act) in enumerate(a_tp1)]
    δ = tderror(q_t, r, γ, q_tp1)
    Δ = corr_term.*δ.*(1.0/length(q_t))
    for (ci_idx, ci) in enumerate(cis_t)
        model[ci] -= opt.eta*Δ[ci_idx]
    end
end

function update!(model::TabularLayer, opt::Descent, lu::BatchSarsa, ρ, s_t::Array{Int64, 1}, s_tp1::Array{Int64, 1}, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)

    cis_t = [CartesianIndex(act, s_t[idx]...) for (idx, act) in enumerate(a_t)]
    q_t = model.(cis_t)
    q_tp1 = [model(CartesianIndex(act, s_tp1[idx]...)) for (idx, act) in enumerate(a_tp1)]
    δ = tderror(q_t, r, γ, q_tp1)
    Δ = corr_term.*δ.*(1.0/length(q_t))
    for (ci_idx, ci) in enumerate(cis_t)
        model[ci] -= opt.eta*Δ[ci_idx]
    end
end

mutable struct BatchExpectedSarsa <: LearningUpdate end

function update!(model, opt, lu::BatchExpectedSarsa, ρ, s_t, s_tp1, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)
    preds_t = [model(s_t[idx])[act] for (idx, act) in enumerate(a_t)]
    exp_q = [dot(target_policy, model(s_tp1[idx])) for (idx, act) in enumerate(a_tp1)]
    grads = Flux.gradient(()->offpolicy_tdloss(preds, r, γ, Flux.data.(exp_q)), params(model))
    for weights in params(model)
        update!(opt, weights, -grads[weights])
    end
end

function update!(model::SingleLayer, opt, lu::BatchExpectedSarsa, ρ, s_t, s_tp1, r, γ, terminal, a_t, a_tp1, target_policy; corr_term=1.0)
    q_t = [model(s_t[idx])[act] for (idx, act) in enumerate(action_t)]
    exp_q = [dot(target_policy, model(s_tp1[idx])) for (idx, act) in enumerate(a_tp1)]
    dvdt = [deriv(model, s_t[idx]) for (idx, act) in enumerate(action_t)]
    δ = tderror(v_t, r, γ, exp_q)
    Δ = mean(δ.*dvdt.*s_t)
    model.W .+= -apply!(opt, model.W, corr_term*Δ)
end

function update!(model::TabularLayer, opt::Descent, lu::BatchExpectedSarsa, ρ, s_t::Array{Array{Int64, 1}, 1}, s_tp1::Array{Array{Int64, 1}, 1}, r, γ, terminal, a_t, a_tp1, target_policy::Array{Float64, 1}; corr_term=1.0)
    q_t = [model(CartesianIndex(act, s_t[idx]...)) for (idx, act) in enumerate(a_t)]
    exp_q = [dot(target_policy, (model[:, s...])) for s in s_tp1]
    δ = tderror(q_t, r, γ, exp_q)
    Δ = corr_term.*δ.*(1.0/length(q_t))
    for (s_idx, s) in enumerate(s_t)
        model[a_t[s_idx], s...] -= opt.eta*Δ[s_idx]
    end
end

function update!(model::TabularLayer, opt::Descent, lu::BatchExpectedSarsa, ρ, s_t::Array{Int64, 1}, s_tp1::Array{Int64, 1}, r, γ, terminal, a_t, a_tp1, target_policy::Array{Float64, 1}; corr_term=1.0)
    q_t = [model(CartesianIndex(act, s_t[idx]...)) for (idx, act) in enumerate(a_t)]
    exp_q = [dot(target_policy, (model[:, s...])) for s in s_tp1]
    δ = tderror(q_t, r, γ, exp_q)
    Δ = corr_term.*δ.*(1.0/length(q_t))
    for (s_idx, s) in enumerate(s_t)
        model[a_t[s_idx], s...] -= opt.eta*Δ[s_idx]
    end
end

function update!(model::TabularLayer, opt::Descent, lu::BatchExpectedSarsa, ρ, s_t::Array{Array{Int64, 1}, 1}, s_tp1::Array{Array{Int64, 1}, 1}, r, γ, terminal, a_t, a_tp1, target_policy::P; corr_term=1.0) where {P<:AbstractPolicy}
    q_t = [model(CartesianIndex(act, s_t[idx]...)) for (idx, act) in enumerate(a_t)]
    # println([get(target_policy, s, a, nothing, nothing, nothing) for a in 1:size(model.W)[1]])
    exp_q = [dot([get(target_policy, s, a, nothing, nothing, nothing) for a in 1:size(model.W)[1]], (model[:, s...])) for s in s_tp1]
    δ = tderror(q_t, r, γ, exp_q)
    Δ = corr_term.*δ.*(1.0/length(q_t))
    for (s_idx, s) in enumerate(s_t)
        model[a_t[s_idx], s...] -= opt.eta*Δ[s_idx]
    end
end


text:
# ----------------------------------------------------
# --- Import modules
# ----------------------------------------------------
using DigitalComm
using DSP
using Test
# ----------------------------------------------------
# --- Tests
# ----------------------------------------------------
println("Tests for filters");
@testset "Raised Cosine filter test" begin
    # Create the vector
    h = raisedCosine(12,0.5,16)
    # Check type
    @test h isa Vector{Float64}
    # Check size
    sH = 12*16*2+1
    @test length(h)== sH
    # Check it is 1 in the middle
    @test h[1+ (sH-1)÷2] == 1
    # Check Nyquist criterion
    for n ∈ 0 : 11
        @test ≈(h[1 + n*16],0,atol=1e-8)
    end
end

@testset "Raised Cosine filter test" begin
    # Create the vector
    h = sqrtRaisedCosine(12,0.5,16)
    # Check type
    @test h isa Vector{Float64}
    # Check size
    sH = 12*16*2+1
    @test length(h)== sH
    # Check it is 1 in the middle
    @test h[1+ (sH-1)÷2] == 1
    # Check Nyquist criterion
    # sqrt is not a Nyquist filter but h*h is !
    p = conv(h,h)
    p = p /maximum(p)
    for n ∈ 0 : 2*11
        @test ≈(p[1 + n*16],0,atol=1e-4)
    end
end


text:
# Tests for conditional entropy function

# Test with 1 values array conditioned on itself
@test get_conditional_entropy(arr1, arr1) ≈ 0 atol = 0.05
println("Conditional entropy with 1 array passed.")

# Test conditional entropy formula
@test get_conditional_entropy(arr1, arr2) ≈ get_entropy(arr1, arr2) - get_entropy(arr2) atol = 0.05
println("Conditional entropy formula passed.")


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

"""
The `ChangePrecision` module exports a macro `@changeprecision T expression`
that changes the "default" floating-point precision in a given `expression`
to a new floating-point type `T`.
"""
module ChangePrecision

import Random, Statistics, LinearAlgebra
using Random: AbstractRNG

## Note: code in this module must be very careful with math functions,
#        because we've defined module-specific versions of very basic
#        functions like + and *.   Call Base.:+ etcetera if needed.

export @changeprecision

############################################################################
# The @changeprecision(T, expr) macro, below, takes calls to
# functions f that default to producing Float64 (e.g. from integer args)
# and converts them to calls to ChangePrecision.f(T, args...).  Then
# we implement our f(T, args...) to default to T instead.  The following
# are a list of function calls to transform in this way.

const randfuncs = (:rand, :randn, :randexp) # random-number generators
const matfuncs = (:ones, :zeros) # functions to construct arrays
const complexfuncs = (:abs, :angle) # functions that give Float64 for Complex{Int}
const binaryfuncs = (:*, :+, :-, :^) # binary functions on irrationals that make Float64

# math functions that convert integer-like arguments to floating-point results
# (from https://docs.julialang.org/en/release-0.6/manual/mathematical-operations/, up to date as of 0.6)
const intfuncs = (:/, :\, :inv, :float,
                  # powers logs and roots
                  :√,:∛,:sqrt,:cbrt,:hypot,:exp,:exp2,:exp10,:expm1,:log,:log2,:log10,:log1p,:cis,
                  # trig
                  :sin,    :cos,    :tan,    :cot,    :sec,    :csc,
                  :sinh,   :cosh,   :tanh,   :coth,   :sech,   :csch,
                  :asin,   :acos,   :atan,   :acot,   :asec,   :acsc,
                  :asinh,  :acosh,  :atanh,  :acoth,  :asech,  :acsch,
                  :sinc,   :cosc,   :atan2,
                  :cospi,  :sinpi,
                  # trig in degrees
                  :deg2rad,:rad2deg,
                  :sind,   :cosd,   :tand,   :cotd,   :secd,   :cscd,
                  :asind,  :acosd,  :atand,  :acotd,  :asecd,  :acscd,
                  )

# functions that convert integer arrays to floating-point results
const statfuncs = (:mean, :std, :stdm, :var, :varm, :median, :cov, :cor)
const linalgfuncs = (:opnorm, :norm, :normalize,
                     :factorize, :cholesky, :bunchkaufman, :ldlt, :lu, :qr, :lq,
                     :eigen, :eigvals, :eigfact, :eigmax, :eigmin, :eigvecs,
                     :hessenberg, :schur, :svd, :svdvals,
                     :cond, :condskeel, :det, :logdet, :logabsdet,
                     :pinv, :nullspace, :lyap, :sylvester)

# functions to change to ChangePrecision.func(T, ...) calls:
const changefuncs = Set([randfuncs..., matfuncs...,
                         intfuncs..., complexfuncs...,
                         statfuncs..., linalgfuncs...,
                         binaryfuncs..., :include])

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

changeprecision(T, x) = x
changeprecision(T::Type, x::Float64) = parse(T, string(x)) # change float literals
function changeprecision(T, x::Symbol)
    if x ∈ (:Inf, :NaN)
        return :(convert($T, $x))
    else
        return x
    end
end
function changeprecision(T, x::Float64)
    if T === :Float16
        return Float16(x)
    elseif T === :Float32
        return Float32(x)
    elseif T === :Float64
        return x
    else
        return :(parse($T, $(string(x))))
    end
end
function changeprecision(T, ex::Expr)
    if Meta.isexpr(ex, :call, 3) && ex.args[1] == :^ && ex.args[3] isa Int
        # mimic Julia 0.6/0.7's lowering to literal_pow
        return Expr(:call, ChangePrecision.literal_pow, T, :^, changeprecision(T, ex.args[2]), Val{ex.args[3]}())
    elseif Meta.isexpr(ex, :call, 2) && ex.args[1] == :include
        return :($include($T, @__MODULE__, $(ex.args[2])))
    elseif Meta.isexpr(ex, :call) && ex.args[1] in changefuncs
        return Expr(:call, Core.eval(ChangePrecision, ex.args[1]), T, changeprecision.(T, ex.args[2:end])...)
    elseif Meta.isexpr(ex, :., 2) && ex.args[1] in changefuncs && Meta.isexpr(ex.args[2], :tuple)
        return Expr(:., Core.eval(ChangePrecision, ex.args[1]), Expr(:tuple, T, changeprecision.(T, ex.args[2].args)...))
    elseif Meta.isexpr(ex, :macrocall) && ex.args[1] == Symbol("@__dot__")
        return changeprecision(T, Base.Broadcast.__dot__(ex.args[end]))
    else
        return Expr(ex.head, changeprecision.(T, ex.args)...)
    end
end

# calls to include(f) are changed to include(T, mod, f) so that
# @changeprecision can apply recursively to included files.
function include(T, mod, filename::AbstractString)
    # use the undocumented parse_input_line function so that we preserve
    # the filename and line-number information.
    s = string("begin; ", read(filename, String), "\nend\n")
    expr = Base.parse_input_line(s, filename=filename)
    Core.eval(mod, changeprecision(T, expr))
end

"""
    @changeprecision T expression

Change the "default" precision in the given `expression` to the floating-point
type `T`.

This changes floating-point literals, integer expressions like `1/3`,
random-number functions like `rand`, and matrix constructors like `ones`
to default to the new type `T`.

For example,
```
@changeprecision Float32 begin
    x = 7.3
    y = 1/3
    z = rand() .+ ones(3,4)
end
```
uses `Float32` precision for all of the expressions in `begin ... end`.
"""
macro changeprecision(T, expr)
    esc(changeprecision(T, expr))
end

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

# integer-like types that get converted to Float64 by various functions
const HWInt = Union{Bool,Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128}
const RatLike = Union{Rational{<:HWInt}, Complex{<:Rational{<:HWInt}}}
const IntLike = Union{HWInt, Complex{<:HWInt}}
const IntRatLike = Union{IntLike,RatLike}
const Promotable = Union{IntLike, RatLike, Irrational}
const PromotableNoRat = Union{IntLike, Irrational}

@inline tofloat(T, x) = T(x)
@inline tofloat(::Type{T}, x::Complex) where {T<:Real} = Complex{T}(x)
@inline tofloat(T, x::AbstractArray) = copyto!(similar(x, T), x)
@inline tofloat(::Type{T}, x::AbstractArray{<:Complex}) where {T<:Real} = copyto!(similar(x, Complex{T}), x)

###########################################################################
# ChangePrecision.f(T, args...) versions of Base.f(args...) functions.

# define our own versions of rand etc. that override the default type,
# which which still respect a type argument if it is explicitly provided
for f in randfuncs
    @eval begin
        $f(T) = Random.$f(T)
        $f(T, dims::Integer...) = Random.$f(T, dims...)
        $f(T, dims::Tuple{<:Integer}) = Random.$f(T, dims)
        $f(T, rng::AbstractRNG, dims::Integer...) = Random.$f(rng, T, dims...)
        $f(T, rng::AbstractRNG, dims::Tuple{<:Integer}) = Random.$f(rng, T, dims)
        $f(T, args...) = Random.$f(args...)
    end
end

# similarly for array constructors like ones
for f in matfuncs
    @eval begin
        $f(T) = Base.$f(T)
        $f(T, dims::Integer...) = Base.$f(T, dims...)
        $f(T, dims::Tuple{<:Integer}) = Base.$f(T, dims)
        $f(T, args...) = Base.$f(args...)
    end
end

# we want to change expressions like 1/2 to produce the new floating-point type
for f in intfuncs
    @eval begin
        $f(T, n::Promotable) = Base.$f(tofloat(T, n))
        $f(T, m::Promotable, n::Promotable) = Base.$f(tofloat(T, m), tofloat(T, n))
        $f(T, args...) = Base.$f(args...)
    end
end

# exception to intfuncs above: division on rationals produces an exact rational
inv(T, x::RatLike) = Base.inv(x)
/(T, y::IntLike, x::RatLike) = Base.:/(y, x)
\(T, x::RatLike, y::IntLike) = Base.:\(x, y)
/(T, y::RatLike, x::IntLike) = Base.:/(y, x)
\(T, x::IntLike, y::RatLike) = Base.:\(x, y)
/(T, y::RatLike, x::RatLike) = Base.:/(y, x)
\(T, x::RatLike, y::RatLike) = Base.:\(x, y)

for f in complexfuncs
    @eval begin
        $f(T, z::Union{Complex{<:HWInt},Complex{<:Rational{<:HWInt}}}) = Base.$f(tofloat(T, z))
        $f(T, args...) = Base.$f(args...)
    end
end

for f in binaryfuncs
    @eval begin
        $f(T, x::Irrational, y::Promotable) = Base.$f(tofloat(T, x), tofloat(T, y))
        $f(T, x::Promotable, y::Irrational) = Base.$f(tofloat(T, x), tofloat(T, y))
        $f(T, x::Irrational, y::Irrational) = Base.$f(tofloat(T, x), tofloat(T, y))
        $f(T, args...) = Base.$f(args...)
    end
end
-(T::Type, x::Irrational) = Base.:-(tofloat(T, x))
for f in (:+, :*) # these functions can accept 3+ arguments
    # FIXME: these methods may be slow compared to the built-in + or *
    #        because they do less inlining?
    @eval begin
        @inline $f(T, x::Promotable, y, z, args...) = $f(T, x, $f(T, y, z, args...))
        @inline $f(T, x::IntRatLike, y::IntRatLike, z::IntRatLike, args::IntRatLike...) = Base.$f(x, y, z, args...)
    end
end

^(T, x::Union{AbstractMatrix{<:Promotable},Promotable}, y::Union{RatLike,Complex{<:HWInt}}) = Base.:^(tofloat(T, x), y)

# e^x is handled specially
const esym = VERSION < v"0.7.0-DEV.1592" ? :e : :ℯ # changed in JuliaLang/julia#23427
^(T, x::Irrational{esym}, y::Promotable) = Base.exp(tofloat(T, y))
literal_pow(T, op, x::Irrational{esym}, ::Val{n}) where {n} = Base.exp(tofloat(T, n))

# literal integer powers are specially handled in Julia
if VERSION < v"0.7.0-DEV.843" # JuliaLang/julia#22475
    literal_pow(T, op, x::Irrational, ::Val{n}) where {n} = Base.literal_pow(op, tofloat(T, x), Val{n})
    @inline literal_pow(T, op, x, ::Val{n}) where {n} = Base.literal_pow(op, x, Val{n})
else
    literal_pow(T, op, x::Irrational, p) = Base.literal_pow(op, tofloat(T, x), p)
    @inline literal_pow(T, op, x, p) = Base.literal_pow(op, x, p)
end

for f in (statfuncs...,linalgfuncs...)
    m = f ∈ statfuncs ? :Statistics : :LinearAlgebra
    # for functions like factorize, if we are converting the matrix to floating-point
    # anyway then we might as well call factorize! instead to overwrite our temp array:
    if f ∈ (:factorize, :cholesky, :bunchkaufman, :ldlt, :lu, :qr, :lq, :eigen, :svd, :eigvals!, :svdvals!, :median)
        f! = Symbol(f, :!)
        @eval begin
            $f(T, x::AbstractArray{<:Promotable}, args...; kws...) = $m.$f!(tofloat(T, x), args...; kws...)
            $f(T, x::AbstractArray{<:Promotable}, y::AbstractArray{<:Promotable}, args...; kws...) = $m.$f!(tofloat(T, x), tofloat(T, y), args...; kws...)
        end
    else
        @eval begin
            $f(T, x::AbstractArray{<:Promotable}, args...; kws...) = $m.$f(tofloat(T, x), args...; kws...)
            $f(T, x::AbstractArray{<:Promotable}, y::AbstractArray{<:Promotable}, args...; kws...) = $m.$f(tofloat(T, x), tofloat(T, y), args...; kws...)
        end
    end
    @eval begin
        $f(T, x::AbstractArray{<:Promotable}, y::AbstractArray, args...; kws...) = $m.$f(x, y, args...; kws...)
        $f(T, args...; kws...) = $m.$f(args...; kws...)
    end
end
for f in (:varm, :stdm) # look at type of second (scalar) argument
    @eval begin
        $f(T, x::AbstractArray{<:Promotable}, m::Union{AbstractFloat,Complex{<:AbstractFloat}}, args...; kws...) = Statistics.$f(x, m, args...; kws...)
        $f(T, x::AbstractArray{<:PromotableNoRat}, m::PromotableNoRat, args...; kws...) = Statistics.$f(tofloat(T, x), tofloat(T, m), args...; kws...)
    end
end
inv(T, x::AbstractArray{<:PromotableNoRat}) = Base.inv(tofloat(T, x))
/(T, x::AbstractArray{<:Promotable}, y::Union{PromotableNoRat,AbstractArray{<:PromotableNoRat}}) = Base.:/(tofloat(T, x), tofloat(T, y))
\(T, y::Union{PromotableNoRat,AbstractArray{<:PromotableNoRat}}, x::AbstractArray{<:Promotable}) = Base.:\(tofloat(T, y), tofloat(T, x))

# more array functions that are exact for rationals: don't convert
for f in (:mean, :median, :var, :std, :cor, :cov, :ldlt, :lu)
    m = f ∈ statfuncs ? :Statistics : :LinearAlgebra
    @eval begin
        $f(T, x::AbstractArray{<:RatLike}, y::AbstractArray{<:Promotable}, args...; kws...) = $m.$f(tofloat(T, x), tofloat(T, y), args...; kws...)
        $f(T, x::AbstractArray{<:RatLike}, y::AbstractArray{<:RatLike}, args...; kws...) = $m.$f(x, y, args...; kws...)
        $f(T, x::AbstractArray{<:RatLike}, args...; kws...) = $m.$f(x, args...; kws...)
    end
end

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

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
"""https://github.com/moruten/julia-code/blob/99b50bf329195ba102a9496602a22138f07d186b/Test-2021-8-12-no6.ipynb"""
module Original

using Plots
using Random

function any_width_rand!(n_rand,width_l,width_r)
    rands=zeros(n_rand)
    L=width_r - width_l
    for i=1:n_rand
        rands[i] = width_l + L*rand()
    end
    return rands
end

function action_cal(x)
    exp(0.5*x^2)
end

function MP_test(x_ini,x_pro)
    S_ini = action_cal(x_ini)
    S_pro = action_cal(x_pro)

    if exp(S_ini - S_pro) > rand()
        return x_pro
    else
        return x_ini
    end
end

n_rand=1_000_000

MP_rands=zeros(n_rand)

x_propose  = 0.0
x_ini      = 0.0
width_rand = 2.0

@time for i=2:n_rand
    x_propose   = x_ini + any_width_rand!(1,-width_rand,width_rand)[1]
    MP_rands[i] = MP_test(x_ini,x_propose)
    global x_ini       = MP_rands[i]
end

f(x) = exp(-0.5*x^2)/sqrt(2*pi)
xr   = range(-3,3,length=10001)

plot(MP_rands,st=:histogram,norm=:true,nbins=100,label="Metropolis method")
plot!(xr,f.(xr),lw=3,label="Gaussian") |> display

end

# %%
"""https://github.com/moruten/julia-code/blob/99b50bf329195ba102a9496602a22138f07d186b/Test-2021-8-12-no6.ipynb"""
module Original_plus

using Plots
using Random

function any_width_rand!(n_rand,width_l,width_r)
    rands=zeros(n_rand)
    L=width_r - width_l
    for i=1:n_rand
        rands[i] = width_l + L*rand()
    end
    return rands
end

function action_cal(x)
    exp(0.5*x^2)
end

function MP_test(x_ini,x_pro)
    S_ini = action_cal(x_ini)
    S_pro = action_cal(x_pro)

    if exp(S_ini - S_pro) > rand()
        return x_pro
    else
        return x_ini
    end
end

n_rand=1_000_000

MP_rands=zeros(n_rand)

x_propose  = 0.0
x_ini      = 0.0
width_rand = 2.0

@time for i=2:n_rand
    x_propose   = x_ini + any_width_rand!(1,-width_rand,width_rand)[1]
    MP_rands[i] = MP_test(x_ini,x_propose)
    global x_ini       = MP_rands[i]
end

f(x) = exp(-0.5*x^2)/sqrt(2*pi)
xr   = range(-3,3,length=10001)

using QuadGK
G(x) = exp(-exp(0.5*x^2))
Z = quadgk(G, -Inf, Inf)[1]
g(x) = G(x)/Z

histogram(MP_rands, norm=true, nbins=100, label="Metropolis method", alpha=0.3)
plot!(xr, f.(xr), lw=3, label="Gaussian")
plot!(xr, g.(xr), lw=3, label="exp(-exp(0.5*x^2))/Z", ls=:dash) |> display

end

# %%
module Corrected

using Plots
using Random

function any_width_rand!(n_rand,width_l,width_r)
    rands=zeros(n_rand)
    L=width_r - width_l
    for i=1:n_rand
        rands[i] = width_l + L*rand()
    end
    return rands
end

function action_cal(x)
    exp(0.5*x^2)
end

function MP_test(x_ini,x_pro)
    S_ini = action_cal(x_ini)
    S_pro = action_cal(x_pro)

    if S_ini / S_pro > rand() # <-- 変更はこの行だけ！
        return x_pro
    else
        return x_ini
    end
end

n_rand=1_000_000

MP_rands=zeros(n_rand)

x_propose  = 0.0
x_ini      = 0.0
width_rand = 2.0

@time for i=2:n_rand
    x_propose   = x_ini + any_width_rand!(1,-width_rand,width_rand)[1]
    MP_rands[i] = MP_test(x_ini,x_propose)
    global x_ini       = MP_rands[i]
end

f(x) = exp(-0.5*x^2)/sqrt(2*pi)
xr   = range(-5,5,length=10001)

plot(MP_rands,st=:histogram,norm=:true,nbins=100,label="Metropolis method", alpha=0.3)
plot!(xr,f.(xr),lw=3,label="Gaussian") |> display

end

# %%
module MH

using Plots
using Random

function any_width_rand!(n_rand,width_l,width_r)
    rands=zeros(n_rand)
    L=width_r - width_l
    for i=1:n_rand
        rands[i] = width_l + L*rand()
    end
    return rands
end

function action_cal(x)
    exp(0.5*x^2)
end

function MP_test(x_ini,x_pro)
    S_ini = action_cal(x_ini)
    S_pro = action_cal(x_pro)
    a = S_ini/S_pro
    if a > 1 || rand() ≤ a
        return x_pro
    else
        return x_ini
    end
end

n_rand=1_000_000

MP_rands=zeros(n_rand)

x_propose  = 0.0
x_ini      = 0.0
width_rand = 2.0

@time for i=2:n_rand
    x_propose   = x_ini + any_width_rand!(1,-width_rand,width_rand)[1]
    MP_rands[i] = MP_test(x_ini,x_propose)
    global x_ini       = MP_rands[i]
end

f(x) = exp(-0.5*x^2)/sqrt(2*pi)
xr   = range(-5,5,length=10001)

plot(MP_rands,st=:histogram,norm=:true,nbins=100,label="MH method", alpha=0.3)
plot!(xr,f.(xr),lw=3,label="Gaussian") |> display

end

# %%
using Plots
using Distributions

function update(H, prop, x)
    y = prop(x)
    a = exp(H(x) - H(y))
    a > 1 || rand() ≤ a ? y : x
end

function iter_update(H, prop, x0, N; warmup=min(N ÷ 10, 10^3))
    x = x0
    for _ in 1:warmup
        x = update(H, prop, x)
    end
    X = similar([x], N)
    X[1] = x
    for i in 2:N
        X[i] = update(H, prop, X[i-1])
    end
    X
end

H(x) = x^2/2
prop(x) = x + 0.2randn()
x0 = randn()
N = 10^6
X = iter_update(H, prop, x0, N)

histogram(X; norm=true, alpha=0.3, bin=100, label="MH sample")
plot!(x -> pdf(Normal(), x), -5, 5; ls=:dash, lw=2, label="normal dist.")

# %%
H(x) = x ≤ 0 ? Inf : x
prop(x) = x + 0.1randn()
x0 = abs(randn())
N = 10^6
X = iter_update(H, prop, x0, N)

histogram(X; norm=true, alpha=0.3, bin=100, label="MH sample")
plot!(x -> pdf(Exponential(), x), eps(), 10; label="exp. dist.")
plot!(xlim=(-0.1, 8))

# %%
H(x) = x ≤ 0 ? Inf : x/3 - (4 - 1)*log(x)
prop(x) = x + randn()
x0 = abs(randn())
N = 10^6
X = iter_update(H, prop, x0, N)

histogram(X; norm=true, alpha=0.3, bin=100, label="MH sample")
plot!(x -> pdf(Gamma(4, 3), x), eps(), 50; label="Gamma(4, 3)")

# %%

max_stars_repo_path: 0015/Metropolis method.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/Metropolis method.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/Metropolis method.jl

max_forks_repo_name: genkuroki/public


text:
import Random.seed!
seed!(model::Model) = Random.seed!(model.rng)
seed!(model::Model, seed...) = Random.seed!(model.rng, seed...)

export Model, Ising, XY, Potts, Clock, AshkinTeller
export QuantumXXZ

include("common/union_find.jl")

## Classical
include("Ising/Ising.jl")
include("Potts/Potts.jl")
include("Clock/Clock.jl")
include("XY/XY.jl")
include("AshkinTeller/AshkinTeller.jl")

## Quantum
include("common/LoopOperator.jl")
include("QuantumXXZ/QuantumXXZ.jl")


text:
using FPAging

dy = 0.1
y_range = (-3.0:dy:10.0)[2:end-1]

# this does not work
# My problem is related to https://discourse.julialang.org/t/exporting-a-runtimegeneratedfunction/63628
# but I can't figure out how to get this to work

d = FPDist(0.35,0.0048,1.06,0.147,0.33,0.29,-0.5,0.1,40.0,60.0,y_range)
rand(d,3)
logpdf(d,[-1.0,0.0,1.0,2.0])

max_stars_repo_name: hstrey/FPAging

max_issues_repo_name: hstrey/FPAging

max_forks_repo_name: hstrey/FPAging


text:
export filter_reads, filter_muts

function _plot_stat(dir)
    dir == nothing ? nothing : function(args...)
        open(joinpath(dir, string(join(map(x->x.name, args[1:length(args)÷2]), '_'), ".html")), "w") do f
            plot_stat(f, args...)
        end
    end
end

function filter_reads(rules, reads; statdir=nothing, debug=false)
    out = Read[]
    build_rule_func(false, rules, function(read, pass)
        if !pass
            read.flag |= 0x0200
        end
        push!(out, read)
    end, _plot_stat(statdir), debug=debug)(reads)
    out
end

function filter_muts(f, rules, bam, reads; statdir=nothing, debug=false)
    build_rule_func(true, rules, function(reads, chr, mut, info, filters)
        dp = length(reads)
        ad = count(reads) do r
            map_to_read(mut, r) in r.muts
        end
        rd = dp - ad
        write_vcf_line(f, bam, chr, mut, info, filters, (dp, rd, ad))
    end, _plot_stat(statdir), debug=debug)(@task pileup(reads))
end


text:
# Autogenerated wrapper script for librdkafka_jll for i686-w64-mingw32-cxx03
export librdkafka

using Lz4_jll
using Zstd_jll
using Zlib_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("librdkafka")
JLLWrappers.@declare_library_product(librdkafka, "librdkafka.dll")
function __init__()
    JLLWrappers.@generate_init_header(Lz4_jll, Zstd_jll, Zlib_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        librdkafka,
        "bin\\librdkafka.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
  Tug of war
   From Hakaru example documentation/tugofwar_rejection.hk

  Compare with true_skill.jl

  cf ~/webppl/tug_of_war.wppl
     ~/blog/tug_of_war.blog
=#
using Turing, StatsPlots, DataFrames

#=
# For num_people = 5
strength (latent):
strength[1] > strength[2]: 0.67203
strength[1] > strength[3]: 0.641555
strength[1] > strength[4]: 0.76379
strength[1] > strength[5]: 0.73293
strength[2] > strength[3]: 0.60609
strength[2] > strength[4]: 0.69859
strength[2] > strength[5]: 0.688345
strength[3] > strength[4]: 0.54991
strength[3] > strength[5]: 0.640705
strength[4] > strength[5]: 0.49402
pulls (observed):
pulls[1] > pulls[2]: 1.0
pulls[1] > pulls[3]: 1.0
pulls[1] > pulls[4]: 1.0
pulls[1] > pulls[5]: 1.0
pulls[2] > pulls[3]: 1.0
pulls[2] > pulls[4]: 1.0
pulls[2] > pulls[5]: 1.0
pulls[3] > pulls[4]: 0.99546
pulls[3] > pulls[5]: 0.99546
pulls[4] > pulls[5]: 1.0

=#
@model function tug_of_war(n)

    strength = TArray{Real}(undef, n)
    pulls = TArray{Real}(undef, n)
    # strength = Vector{Real}(undef, n)
    # pulls = Vector{Real}(undef, n)

    for p in 1:n
        strength[p] ~ Uniform(0,1) # truncated(Normal(0,1),0,10) # latent
        pulls[p]    ~ truncated(Normal(strength[p],1),0,10) # observed # truncated(Normal(strength[p],1),0,10) # observed
    end

    # The Hakary/BLOG model had "pulls(a) < pulls(b) "
    # which is a little strange so I change to return the winner.
    function winner(a,b)
        return pulls[a] > pulls[b] ? a : b
    end

    # match1 = winner(alice, bob)
    # match2 = winner(bob, carol)
    # match3 = winner(alice, carol)

    # match1 = winner(1, 2)
    # match2 = winner(2, 3)
    # match3 = winner(1, 3)

    # A little more general approach
    for p1 in 1:n-1, p2 in p1+1:n
        true ~ Dirac(winner(p1,p2) == p1)
    end

end

num_people = 3
model = tug_of_war(num_people)

num_chains = 4
num_samples = num_people * 10_000 # 10_000

# High rhat
chains = sample(model,  MH(), MCMCThreads(), num_samples, num_chains)

# PG(20) doesn't differentiate the people's pull/strength
# chains = sample(model,  PG(20), MCMCThreads(), num_samples, num_chains)

# chains = sample(model,  HMC(0.1, 5), MCMCThreads(), num_samples, num_chains)

# IS() doesn't differentiate the people's pull/strength
# chains = sample(model,  IS(), MCMCThreads(), num_samples, num_chains)

# High rhat. Alot of rejected proposals
# chains = sample(model,  HMC(0.1, 5), MCMCThreads(), num_samples, num_chains)

# chains = sample(model,  HMCDA(200, 0.65, 0.3), MCMCThreads(), num_samples, num_chains)

# chains = sample(model,  NUTS(1000, 0.65), MCMCThreads(), num_samples, num_chains)

# Too slow
# chains = sample(model,  SMC(), MCMCThreads(), num_samples, num_chains)

# Good: slightly high rhat
# chains = sample(model,  NUTS(1000,0.65), MCMCThreads(), num_samples, num_chains)

display(chains)
display(gelmandiag(chains))
# display(plot(chains))
df = DataFrame(chains)

# Show the latent strength
println("strength (latent):")
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("strength[$p1] > strength[$p2]: ", mean(df[!,"strength[$p1]"] .>df[!,"strength[$p2]"]))
end

println("pulls (observed):")
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("pulls[$p1] > pulls[$p2]: ", mean(df[!,"pulls[$p1]"] .>df[!,"pulls[$p2]"]))
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"""
    mese! = MESEBlochSim(TR, TE, nechoes, [rfex, rfref])
    mese!(spin, [workspace])

Simulate a multi-echo spin echo (MESE) scan on `spin`, overwriting the spin's
magnetization vector. Returns a `Vector` with the magnetization vectors at each
echo.

# Arguments
- `TR::Real`: Repetition time (ms)
- `TE::Real`: First echo time, and echo spacing (ms)
- `nechoes::Integer`: Number of echoes to readout
- `rfex::AbstractRF = InstantaneousRF(π/2)`: Excitation RF pulse
- `rfref::AbstractRF = InstantaneousRF(π, -π/2)`: Refocussing RF pulse

`workspace isa MESEBlochSimWorkspace`.
"""
struct MESEBlochSim{T1<:AbstractRF,T2<:AbstractRF}
    TR::Float64
    TE::Float64
    nechoes::Int
    rfex::T1
    rfref::T2

    function MESEBlochSim(TR, TE, nechoes, rfex::T1, rfref::T2) where {T1,T2}

        TR >= TE * nechoes + duration(rfex) / 2 ||
            error("TR must be greater than or equal to TE * nechoes + duration(rfex) / 2")
        (TE - duration(rfref)) / 2 >= duration(rfex) / 2 ||
            error("the refocussing pulse must not overlap the excitation pulse")
        TE >= duration(rfref) || error("TE must be longer than the refocussing pulse")
        new{T1,T2}(TR, TE, nechoes, rfex, rfref)

    end
end

MESEBlochSim(TR, TE, nechoes) = MESEBlochSim(TR, TE, nechoes, InstantaneousRF(π/2), InstantaneousRF(π, -π/2))

Base.show(io::IO, mese::MESEBlochSim) =
    print(io, "MESEBlochSim(", mese.TR, ", ", mese.TE, ", ", mese.nechoes, ", ", mese.rfex, ", ", mese.rfref, ")")

function Base.show(io::IO, ::MIME"text/plain", mese::MESEBlochSim)

    print(io, "Multi-Echo Spin Echo (MESE) Bloch Simulation:")
    print(io, "\n TR = ", mese.TR, " ms")
    print(io, "\n TE (and echo spacing) = ", mese.TE, " ms")
    print(io, "\n nechoes = ", mese.nechoes)
    print(io, "\n rfex (excitation pulse) = ")
    show(io, "text/plain", mese.rfex)
    print(io, "\n rfref (refocussing pulses) = ")
    show(io, "text/plain", mese.rfref)

end

struct MESEBlochSimWorkspace{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12}
    Aex::T1
    Bex::T2
    Aref::T3
    Bref::T4
    Ate1::T5
    Bte1::T6
    Ate::T5
    Bte::T6
    Atr::T5
    Btr::T6
    Aecho1::T7
    Becho1::T6
    Aecho::T7
    Becho::T6
    tmpA1::T7
    tmpB1::T6
    tmpA2::T7
    tmpB2::T6
    mat::T8
    vec::T9
    bm_workspace::T10
    ex_workspace::T11
    ref_workspace::T12
end

function MESEBlochSimWorkspace(
    spin::AbstractSpin,
    scan::MESEBlochSim,
    bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    MESEBlochSimWorkspace(typeof(spin), typeof(scan), bm_workspace)

end

function MESEBlochSimWorkspace(
    spin::Union{Type{Spin{T}},Type{SpinMC{T,N}}},
    scan::Type{MESEBlochSim{T1,T2}},
    bm_workspace = spin <: Spin ? nothing : BlochMcConnellWorkspace(spin)
) where {T,N,T1,T2}

    if T1 <: InstantaneousRF
        Aex = ExcitationMatrix{T}()
        Bex = nothing
        ex_workspace = nothing
    else
        if spin <: Spin
            Aex = BlochMatrix{T}()
            Bex = Magnetization{T}()
        else
            Aex = BlochMcConnellMatrix{T}(N)
            Bex = MagnetizationMC{T}(N)
        end
        ex_workspace = ExcitationWorkspace(spin, bm_workspace)
    end
    if T2 <: InstantaneousRF
        Aref = ExcitationMatrix{T}()
        Bref = nothing
        ref_workspace = nothing
    else
        if spin <: Spin
            Aref = BlochMatrix{T}()
            Bref = Magnetization{T}()
        else
            Aref = BlochMcConnellMatrix{T}(N)
            Bref = MagnetizationMC{T}(N)
        end
        ref_workspace = ExcitationWorkspace(spin, bm_workspace)
    end
    if spin <: Spin
        Ate1 = FreePrecessionMatrix{T}()
        Bte1 = Magnetization{T}()
        Ate = FreePrecessionMatrix{T}()
        Bte = Magnetization{T}()
        Atr = FreePrecessionMatrix{T}()
        Btr = Magnetization{T}()
        Aecho1 = BlochMatrix{T}()
        Becho1 = Magnetization{T}()
        Aecho = BlochMatrix{T}()
        Becho = Magnetization{T}()
        tmpA1 = BlochMatrix{T}()
        tmpB1 = Magnetization{T}()
        tmpA2 = BlochMatrix{T}()
        tmpB2 = Magnetization{T}()
        mat = Matrix{T}(undef, 3, 3)
        vec = Vector{T}(undef, 3)
    else
        Ate1 = BlochMcConnellMatrix{T}(N)
        Bte1 = MagnetizationMC{T}(N)
        Ate = BlochMcConnellMatrix{T}(N)
        Bte = MagnetizationMC{T}(N)
        Atr = BlochMcConnellMatrix{T}(N)
        Btr = MagnetizationMC{T}(N)
        Aecho1 = BlochMcConnellMatrix{T}(N)
        Becho1 = MagnetizationMC{T}(N)
        Aecho = BlochMcConnellMatrix{T}(N)
        Becho = MagnetizationMC{T}(N)
        tmpA1 = BlochMcConnellMatrix{T}(N)
        tmpB1 = MagnetizationMC{T}(N)
        tmpA2 = BlochMcConnellMatrix{T}(N)
        tmpB2 = MagnetizationMC{T}(N)
        mat = Matrix{T}(undef, 3N, 3N)
        vec = Vector{T}(undef, 3N)
    end
    MESEBlochSimWorkspace(Aex, Bex, Aref, Bref, Ate1, Bte1, Ate, Bte, Atr, Btr,
        Aecho1, Becho1, Aecho, Becho, tmpA1, tmpB1, tmpA2, tmpB2, mat, vec,
        bm_workspace, ex_workspace, ref_workspace)

end

function (scan::MESEBlochSim)(spin::AbstractSpin, workspace::MESEBlochSimWorkspace = MESEBlochSimWorkspace(spin, scan))

    excite!(workspace.Aex, workspace.Bex, spin, scan.rfex, workspace.ex_workspace)
    excite!(workspace.Aref, workspace.Bref, spin, scan.rfref, workspace.ex_workspace)
    freeprecess!(workspace.Ate1, workspace.Bte1, spin, (scan.TE - duration(scan.rfex) - duration(scan.rfref)) / 2, workspace.bm_workspace)
    freeprecess!(workspace.Ate, workspace.Bte, spin, (scan.TE - duration(scan.rfref)) / 2, workspace.bm_workspace)
    freeprecess!(workspace.Atr, workspace.Btr, spin, scan.TR - scan.TE * scan.nechoes, workspace.bm_workspace)
    (S,) = spoil(spin)

    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Ate, workspace.Bte, workspace.Aref, workspace.Bref)
    combine!(workspace.Aecho, workspace.Becho, workspace.tmpA1, workspace.tmpB1, workspace.Ate, workspace.Bte)

    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Ate1, workspace.Bte1, workspace.Aref, workspace.Bref)
    combine!(workspace.Aecho1, workspace.Becho1, workspace.tmpA1, workspace.tmpB1, workspace.Ate, workspace.Bte)
    copyto!(workspace.tmpA1, workspace.Aecho1)
    copyto!(workspace.tmpB1, workspace.Becho1)
    for e = 2:scan.nechoes
        combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Aecho, workspace.Becho)
        copyto!(workspace.tmpA1, workspace.tmpA2)
        copyto!(workspace.tmpB1, workspace.tmpB2)
    end
    combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Atr, workspace.Btr)
    combine!(workspace.tmpA1, workspace.tmpB1, workspace.tmpA2, workspace.tmpB2, S, nothing)
    combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Aex, workspace.Bex)
    subtract!(workspace.mat, I, workspace.tmpA2)
    copyto!(workspace.vec, workspace.tmpB2)
    F = lu!(workspace.mat)
    ldiv!(F, workspace.vec)
    copyto!(spin.M, workspace.vec)

    Mout = Vector{typeof(spin.M)}(undef, scan.nechoes)
    applydynamics!(spin, workspace.tmpB1, workspace.Aecho1, workspace.Becho1)
    Mout[1] = copy(spin.M)
    for e = 2:scan.nechoes
        applydynamics!(spin, workspace.tmpB1, workspace.Aecho, workspace.Becho)
        Mout[e] = copy(spin.M)
    end

    return Mout

end


text:
abstract type AbstractAnalysis end

# TODO document AbstractAnalysis

# modify!(analysis::Type{<:AbstractAnalysis}, eclass::EClass) =
#     error("Analysis does not implement modify!")
modify!(analysis::Type{<:AbstractAnalysis}, g, id) = nothing
join(analysis::Type{<:AbstractAnalysis}, a, b) =
    error("Analysis does not implement join")
make(analysis::Type{<:AbstractAnalysis}, g, a) =
    error("Analysis does not implement make")

islazy(an::Type{<:AbstractAnalysis})::Bool = false


text:
# Autogenerated wrapper script for FLAC_jll for x86_64-w64-mingw32-cxx03
export flac, libflac, libflacpp, metaflac

using Ogg_jll
JLLWrappers.@generate_wrapper_header("FLAC")
JLLWrappers.@declare_library_product(libflac, "libFLAC-8.dll")
JLLWrappers.@declare_library_product(libflacpp, "libFLAC++-6.dll")
JLLWrappers.@declare_executable_product(flac)
JLLWrappers.@declare_executable_product(metaflac)
function __init__()
    JLLWrappers.@generate_init_header(Ogg_jll)
    JLLWrappers.@init_library_product(
        libflac,
        "bin\\libFLAC-8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libflacpp,
        "bin\\libFLAC++-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        flac,
        "bin\\flac.exe",
    )

    JLLWrappers.@init_executable_product(
        metaflac,
        "bin\\metaflac.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export encode, encode_message
export encode_integer

#------------------------------------------------------------------------------
# External API
#------------------------------------------------------------------------------

function encode_message (x)
    encoded, nbytes = encode(x)
    # Add 1 to arity due to magic number
    length_header = encode_integer(uint32(nbytes + 1))
    pkt = [length_header, 131, encoded]
    uint8(pkt), length(pkt)
end

#------------------------------------------------------------------------------
# Encoding Functions
#------------------------------------------------------------------------------

function encode (n::Integer)
    if n < 256
        pkt = [97, n]
    else
        pkt = [98, encode_integer(n)]
    end
    pkt, length(pkt)
end

function encode (n::Real)
    bytes = hex2bytes(num2hex(n))
    pkt = [70, bytes]
    pkt, length(pkt)
end

function encode (n::BigInt)
    sign = n > 0 ? 0 : 1
    size = sizeof(n)
    bytes = zeros(Uint8, sizeof(n))
    for i = 1:size
        bytes[i] = n & 0xFF
        n >>= 8
    end
    arity = length(bytes)
    pkt = [110, arity, sign, bytes]
    pkt, length(pkt)
end

function encode (sym::Symbol)
    bytes = map(uint8, collect(string(sym)))
    arity = encode_integer(uint16(length(bytes)))
    pkt = [100, arity, bytes]
    pkt, length(pkt)
end

function encode (tup::Tuple)
    encoded = Array(Uint8, 0)
    for elt in tup
        pkt, nbytes = encode(elt)
        encoded = vcat(encoded, pkt)
    end

    if length(tup) < 256
        arity = length(tup)
        pkt = [104, arity, encoded]
        pkt, length(pkt)
    else
        arity = encode_integer(uint32(length(tup)))
        pkt = [105, arity, encoded]
        pkt, length(pkt)
    end
end

function encode (str::String)
    if isempty(str)
        pkt = [106]
        pkt, length(pkt)
    else
        arity = encode_integer(uint16(length(str)))
        pkt = [107, arity, collect(str)]
        pkt, length(pkt)
    end
end

function encode (arr::Array)
    if isempty(arr)
        pkt = [106]
        pkt, length(pkt)
    else
        xs = Array(Uint8, 0)
        for i = 1:length(arr)
            bytes, nbytes = encode(arr[i])
            xs = vcat(xs, bytes)
        end

        arity = encode_integer(uint32(length(arr)))
        pkt = [108, arity, xs, 106]
        pkt, length(pkt)
    end
end

function encode (bin::Array{Uint8})
    arity = encode_integer(uint32(length(bin)))
    pkt = [109, arity, bin]
    pkt, length(pkt)
end

#------------------------------------------------------------------------------
# Internal Encoding Functions
#------------------------------------------------------------------------------

function encode_integer (n::Integer)
    sz = sizeof(n)
    bs = Array(Uint8, sz)
    ms = zero(n) | 0xFF
    for i = sz:-1:1
        bs[i] = n & ms
        n >>= 8
    end
    bs
end

max_stars_repo_name: esl/FsBert

max_issues_repo_name: esl/FsBert

max_forks_repo_name: esl/FsBert


text:
#=

# Setting options

To set [Package options](@ref), you can either specify key-value pairs via a `LocalPreferences.toml` file in your environment or by calling `Vulkan.set_preferences!`. Here is an example for setting `LOG_DESTRUCTION` to "true" for debugging purposes:

=#

using Vulkan

Vulkan.set_preferences!("LOG_DESTRUCTION" => "true")

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
"""
   Delegates loading functions & traits for the EXPRES spectrograph
   http://exoplanets.astro.yale.edu/expresBlog/
   https://ui.adsabs.harvard.edu/abs/2016SPIE.9908E..6TJ/abstract
Author: Eric Ford and collaborators
Created: August 2020
"""

import RvSpectMLBase: min_order, max_order, min_pixel_in_order, max_pixel_in_order, min_pixel, max_pixel
min_order(::EXPRES2D) = 1
max_order(::EXPRES2D) = 86
min_pixel_in_order(::EXPRES2D) = 1
max_pixel_in_order(::EXPRES2D) = 7920

min_pixel(::EXPRES1D) = 1
#max_pixel(::EXPRES1D) = orders_to_use_default(EXPRES2D())*(max_col_default(EXPRES2D())-min_col_default(EXPRES2D()+1))

# Values hard coded based on HD 101501.  Someone will need to make automatic or ask if these will stay fixed if these are to be used for anything other than exploring.
const minmax_col_nonnan = [2149:5973, 1920:6121, 1787:6333, 1634:6532, 1485:6744, 1327:6934, 1179:7096, 1061:7242, 956:7282, 853:7292, 777:7291, 711:7292, 678:7292, 639:7292, 581:7292, 547:7292, 518:7296, 456:7297, 382:7296, 365:7293, 356:7289, 348:7289, 336:7292, 318:7295, 295:7297, 272:7296, 258:7296, 247:7297, 238:7295, 229:7294, 222:7292, 216:7291, 211:7288, 206:7285, 202:7282, 199:7279, 197:7278, 193:7275, 189:7275, 183:7273, 180:7273, 173:7273, 167:7274, 162:7274, 157:7273, 154:7272, 150:7271, 146:7269, 143:7268, 140:7267, 137:7266, 133:7265, 130:7264, 127:7262, 124:7261, 121:7260, 118:7258, 115:7258, 112:7256, 110:7256, 107:7254, 104:7253, 101:7252, 99:7250, 96:7248, 94:7247, 91:7245, 89:7243, 87:7241, 85:7239, 83:7237, 81:7234, 79:7230, 78:7227, 76:7223, 75:7218, 75:7211, 75:7204, 76:7194, 79:7181, 84:7167, 93:7142, 109:7097, 157:6352, 500:5119, 3961:3961 ]
min_col_nonnan(::EXPRES2D, ord::Integer) = first(minmax_col_nonnan[ord])
max_col_nonnan(::EXPRES2D, ord::Integer) = last(minmax_col_nonnan[ord])

const minmax_col_excalibur_avail = vcat(fill(1:1,31),
    [ 2658:5073, 1962:5830, 1650:6030, 1550:6364, 1322:6405, 1156:6476, 1063:6848, 843:6927, 779:7053, 672:7136, 469:7172, 487:7177, 399:7148, 331:7203, 337:7191, 257:7196, 264:7201, 226:7207, 221:7212, 228:7218, 189:7206, 184:7212, 168:7200, 163:7206, 158:7212, 153:7200, 148:7206, 143:7194, 150:7200, 145:7207, 140:7194, 135:7201, 130:7188, 138:7195, 133:7182, 128:7189, 123:7176, 132:7183, 127:7170, 122:7177, 131:7164, 126:7150, 121:7113, 145:7098, 155:5401 ],
    fill(1:1, max_order(EXPRES2D())-77+1) )

min_col_excalibur(inst::EXPRES2D, ord::Integer) = first(minmax_col_excalibur_avail[ord])
max_col_excalibur(inst::EXPRES2D, ord::Integer) = last(minmax_col_excalibur_avail[ord])

import RvSpectMLBase: orders_to_use_default, min_col_default, max_col_default

orders_to_use_default(::EXPRES2D) = 43:75    # Based on methods for inital RV described at http://exoplanets.astro.yale.edu/science/activity.php
min_col_default(::EXPRES2D, ord::Integer) = 770   # Seems like they should be good
max_col_default(::EXPRES2D, ord::Integer) = 6650
#min_col_default(::EXPRES2D, ord::Integer) = 1825   # Narrower range for DCPCA
#max_col_default(::EXPRES2D, ord::Integer) = 7000

import RvSpectMLBase: get_pixel_range
function get_pixel_range(inst::EXPRES2D, ord::Integer)
    min_pixels_to_use_in_order = 60
    if max_col_excalibur(inst, ord)-min_col_excalibur(inst, ord) > min_pixels_to_use_in_order
        minc = max(min_col_default(inst, ord), min_col_excalibur(inst,ord), min_col_nonnan(inst,ord))
        maxc = min(max_col_default(inst, ord), max_col_excalibur(inst,ord), max_col_nonnan(inst,ord))
    else
        minc = max(min_col_default(inst, ord), min_col_nonnan(inst,ord))
        maxc = min(max_col_default(inst, ord), max_col_nonnan(inst,ord))
    end
    return minc:maxc
end

import RvSpectMLBase: metadata_symbols_default, metadata_strings_default
metadata_symbols_default(::AnyEXPRES) = Symbol[:midpoint, :bjd, :target, :exposure_time, :airmass, :pwv, :moondist, :sundist, :snr_prelim, :expcount, :wavecal]
metadata_strings_default(::AnyEXPRES) = String["MIDPOINT", "BARYMJD", "OBJECT", "AEXPTIME", "AIRMASS", "PWV", "MOONDIST", "SUNDIST", "channel5_SNR", "EXPCOUNT", "WAVE_CAL"]

#metadata_hdu2_symbols_default(::AnyEXPRES) = Symbol[:S_indicator, :Hα_indicator, :Hα_width, :ccf_width, :σ_ccf_width, :bis_indicator]
#metadata_hdu2_strings_default(::AnyEXPRES) = String["S-VALUE", "HALPHA", "HWIDTH", "CCFFWHM", "CCFFWHME", "BIS"]
# Left out S and BIS since their type is to be standardized in next update of input files
metadata_hdu2_symbols_default(::AnyEXPRES) = Symbol[:Hα_indicator, :Hα_width, :ccf_width, :σ_ccf_width]
metadata_hdu2_strings_default(::AnyEXPRES) = String["HALPHA", "HWIDTH", "CCFFWHM", "CCFFWHME"]

import RvSpectMLBase: default_ccf_mask_v_width
default_ccf_mask_v_width(::AnyEXPRES) = 448.0   # TODO: Update value for EXPRES

import RvSpectMLBase: get_inst_module
get_inst_module(::AnyEXPRES) = EXPRES

import RvSpectMLBase: get_λ_range
function get_λ_range(data::CLT) where { T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
                                       IT<:AnyEXPRES, CLT<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT} }
   mask = data.metadata[:excalibur_mask]
   (λmin, λmax) = extrema(data.λ[mask])
   return (min=λmin, max=λmax)
end

default_filter_line_list_λmin = 4738.0  # Based on HD 101501, should generalize
default_filter_line_list_λmax = 7227.0  #

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
function addxy(x::Int64, y::Int64)
    return sum([x,y])
end


text:
using KomePackage
using Test

@testset "KomePackage.jl group 1" begin
    # Write your own tests here.
    @test KomePackage.greet() == "Hello Test"
    @test KomePackage.greet("Package") == "Hello Package"

end

@testset "KomePackage.jl group 2" begin
    # Write your own tests here.
    @test KomePackage.greet() == "Hello Test"
    @test KomePackage.greet("PackageTest") == "Hello PackageTest"

end


text:
# psi is assumed to be already orthonormalized elsewhere
# Potentials are not updated
function calc_Energies( PW::PWGrid, Potentials, Focc::Array{Float64},
                        psi::Array{ComplexF64,2}, E_NN )

    Ω = PW.Ω
    Ns = PW.Ns
    Npoints = prod(Ns)

    Ngwx = size(psi)[1]
    Nstates = size(psi)[2]
    #
    Kpsi = op_K( PW, psi )
    E_kin = 0.0
    for is = 1:Nstates
        E_kin = E_kin + Focc[is] * real( dot( psi[:,is], Kpsi[:,is] ) )
    end

    # Compute rho
    rho = calc_rho( PW, Focc, psi )

    V_Hartree = real(G_to_R( Ns, Poisson_solve( PW, rho ) )) # 4*pi factor is handled in Poisson_solve
    E_Hartree = 0.5*dot( V_Hartree, rho ) * Ω/Npoints

    E_xc = dot( excVWN(rho), rho ) * Ω/Npoints

    E_ionic = dot( Potentials.Ionic, rho ) * Ω/Npoints

    Energies = EnergiesT(0.0, 0.0, 0.0, 0.0, 0.0, E_NN)
    Energies.Kinetic = E_kin
    Energies.Ionic   = E_ionic
    Energies.Hartree = E_Hartree
    Energies.XC      = E_xc
    Energies.Total   = E_kin + E_ionic + E_Hartree + E_xc + E_NN

    return Energies
end