text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# ╔═╡ 61fd28c2-3ed0-11eb-083a-05e43d8bae12
begin #imports
	using Gadfly
	using DataFrames:DataFrame,select,select!,Not
	using CSV:File
end

# ╔═╡ 95b71286-3ed0-11eb-16a1-c1f5f0281a23
begin #running this cell sets our data and target
data=File("FluxTrain")|>DataFrame
#select!(data,Not(:Column1))
end

# ╔═╡ 87118aa0-4005-11eb-074a-951a3d43a444
plot(data,y="error(ms)")

# ╔═╡ bbb04e40-4005-11eb-0a5e-af47afd2c872
plot(data,x=:Column1,y="error(ms)",Geom.line,color=[colorant"red"])

# ╔═╡ cb220b46-3ed6-11eb-0e87-6f9fb238cc8c
begin
	err=layer(data,y="error(ms)",x=:Column1,Geom.line,color=[colorant"red"])
	left=layer(data,x=:Column1,y=:Ly,color=[colorant"silver",colorant"gold"])
	right=layer(data,y=:Ry,Geom.line,color=[colorant"green",colorant"blue"])
	plot(err,left,right,Guide.title("Everything Juxtaposed"))
end

# ╔═╡ 0839f3a0-4002-11eb-2d2b-2b8611a991a2
vstack(map(plot,[err,left]))

# ╔═╡ bf46a650-4003-11eb-0499-c7e69cb45f96
hstack(map(plot,[err,left]))

# ╔═╡ d0bec360-4056-11eb-1ca6-412c59828f38
plot(data,y="error(ms)",x=:Column1,Geom.line,color=[colorant"red"],Scale.y_sqrt(minvalue=1,maxvalue=100))

# ╔═╡ Cell order:
# ╠═61fd28c2-3ed0-11eb-083a-05e43d8bae12
# ╠═95b71286-3ed0-11eb-16a1-c1f5f0281a23
# ╠═87118aa0-4005-11eb-074a-951a3d43a444
# ╠═bbb04e40-4005-11eb-0a5e-af47afd2c872
# ╠═cb220b46-3ed6-11eb-0e87-6f9fb238cc8c
# ╠═0839f3a0-4002-11eb-2d2b-2b8611a991a2
# ╠═bf46a650-4003-11eb-0499-c7e69cb45f96
# ╠═d0bec360-4056-11eb-1ca6-412c59828f38

max_stars_repo_name: punkmath/Misc

max_issues_repo_name: punkmath/Misc

max_forks_repo_name: punkmath/Misc


text:
__precompile__()

"""
Provides the [`@suffix`](@ref) macro for generating shorthand numeric conversion suffixes.

**Examples**

```julia
x = 1i8        # Convert `1` to `Int8`.
y = 3u32       # Convert `3` to `UInt32`.
z = (x + y)r16 # Convert `4` to `Float16`.
```

**Default Suffixes**

$(length(suffixes)) suffixes are provided by default:

$(join(["- `$(rpad(k, 4)) --> $v`\n" for (k, v) in suffixes]))

Additional user defined suffixes can be made using the `@suffix` macro.
"""
module NumericSuffixes

using Base.Meta, Compat

export @suffix

abstract AbstractSuffix

gentype(n :: Symbol) = :(Base.@__doc__(immutable $n <: $AbstractSuffix end))
gentype(n)           = error("invalid `@suffix` usage. `name` must be a `Symbol`.")

function genfunc(n, ex :: Expr)
    isexpr(ex, :->, 2)      || error("expression must be an anonymous function.")
    isa(ex.args[1], Symbol) || error("anonymous function must have one argument.")
    :(@inline (*)($(ex.args[1]) :: Number, ::Type{$n}) = $(ex.args[2]))
end
genfunc(n, T :: Symbol) = genfunc(n, :(x -> $(T)(x)))
genfunc(n, other)       = error("invalid arguments for `@suffix` macro.")

"""
Generate a new numeric suffix named `n` which converts numbers to type `T`.

**Signature**

```julia
@suffix n T
```

**Examples**

```julia
@suffix i8 Int8

x = 1i8
typeof(x) == Int8
```

More complex conversions can be made using anonymous function syntax:

```julia
@suffix t16 x -> trunc(Int16, x)

x = 1.2t16
typeof(x) == Int16
```
"""
macro suffix(n, T)
    quote
        import Base.Operators: (*)
        $(gentype(n))
        $(genfunc(n, T))
        $n
    end |> esc
end

const suffixes = (
    # Signed integers.
    :i8   => :Int8,
    :i16  => :Int16,
    :i32  => :Int32,
    :i64  => :Int64,
    :i128 => :Int128,
    # Unsigned integers.
    :u8   => :UInt8,
    :u16  => :UInt16,
    :u32  => :UInt32,
    :u64  => :UInt64,
    :u128 => :UInt128,
    # Floating point.
    :r16  => :Float16,
    :r32  => :Float32,
    :r64  => :Float64,
    # Complex floating point.
    :c32  => :Complex32,
    :c64  => :Complex64,
    :c128 => :Complex128,
)

for (n, T) in suffixes
    @eval begin
        export $n
        @doc("Numeric suffix for `$($T)` values.", @suffix($n, $T))
    end
end

end # module

ext: jl


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ e2c752db-738b-429d-bf21-90fecfe0f783
using Images,TestImages

# ╔═╡ 849b2c54-9db4-4b8e-9132-7ebbd7fc4f66
using LinearAlgebra

# ╔═╡ 2fb1ff10-bac7-4f7b-a90d-e6a0cefe2da2
img = float.(testimage("mandrill"))

# ╔═╡ eb6124d9-e0b4-4c74-a7c7-a7f01ab32689
channels = channelview(img)

# ╔═╡ a60655fe-eea5-41ef-bd09-cffc5e1bfb08
begin
	function rank_approx(F::SVD, k)
    U, S, V = F
    M = U[:, 1:k] * Diagonal(S[1:k]) * V[:, 1:k]'
    clamp01!(M)
    end
	svdfactors = svd.(eachslice(channels; dims=1))
	imgs = map((10, 50, 100)) do k
	    colorview(RGB, rank_approx.(svdfactors, k)...)
	end

end

# ╔═╡ 1e3d7470-ae34-4912-95b0-e97de061641d
mosaicview(img, imgs...; nrow=1, npad=5)

# ╔═╡ Cell order:
# ╠═e2c752db-738b-429d-bf21-90fecfe0f783
# ╠═849b2c54-9db4-4b8e-9132-7ebbd7fc4f66
# ╠═2fb1ff10-bac7-4f7b-a90d-e6a0cefe2da2
# ╠═eb6124d9-e0b4-4c74-a7c7-a7f01ab32689
# ╠═a60655fe-eea5-41ef-bd09-cffc5e1bfb08
# ╠═1e3d7470-ae34-4912-95b0-e97de061641d

max_stars_repo_name: BioMedAI/Pluto-Note

max_issues_repo_name: BioMedAI/Pluto-Note

max_forks_repo_name: BioMedAI/Pluto-Note


text:
# functions related to gamma distribution

# R implementations
# For pdf and logpdf we use the Julia implementation
using .RFunctions:
    gammacdf,
    gammaccdf,
    gammalogcdf,
    gammalogccdf,
    gammainvcdf,
    gammainvccdf,
    gammainvlogcdf,
    gammainvlogccdf

# Julia implementations
gammapdf(k::Real, θ::Real, x::Real) = exp(gammalogpdf(k, θ, x))

gammalogpdf(k::Real, θ::Real, x::Real) = gammalogpdf(promote(k, θ, x)...)
function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
    # we ensure that `log(x)` does not error if `x < 0`
    xθ = max(x, 0) / θ
    val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
    return x < 0 ? oftype(val, -Inf) : val
end


text:
# mutable struct MvNormal{T}
#     mu::AbstractVector{T}
#     cov::AbstractMatrix{T}
# end


text:
##analytical gradient- now used mostly for test (too slow for regular use)
DgXki(ti, Xik, wLAsignXik, t, Gkt) = (ti>t)*(Xik/Gkt) + (t==ti)*wLAsignXik

function dfactor!(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T}, DgX; dinfo...) where {TM,TV,T}
  @unpack Gkt, k, t, wLAsignX, dtsM1, dims, X= dinfo

  DgX[:, k] .= DgXki.(dtsM1, X[:,k], wLAsignX[:,k], t, Gkt)

  return nothing
end

#try to compute the analytical gradient
function ∇lossag(G::Matrix{T}, Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterControlLevel;
  RegressionType::Val = Val(PARAM[:iterregressiontype])) where {TM<:CuMatrix,TV<:CuVector,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2),KW=size(Θ.A,1) + size(Xv.W,2),
    Nexpanded=size(Xv.ws,1), expandedinds = collect(1:size(Xv.ws,1)) |> CuVector{Int})
  @unpack v = Xv

  #pre-compute some common values
  X::TM = hcat(initializeÃandfactor(Θ, Xv), Xv.W)
  cuI = TM(I, dims.KW, dims.KW) #need an identity matrix
  XtXinv = cholesky!(X'*X)\cuI
  XXtXinv = X*XtXinv
  XtXinvXtv = XtXinv * (X'*v)
  εt = (X*XtXinvXtv-Xv.v)'

  wLAsignX = signabsμₖ.(Θ.expand.Ã', Θ.expand.LÃ', Xv.ws, Xv.Rtws, Xv.RLws) .* Xv.ws .* abs.(Θ.expand.LÃ')

  #preallocate this carefully high performance
  DgX::TM = CUDA.zeros(T, dims.Nexpanded, dims.KW)

  DG = similar(G)
  Pv = similar(Xv.v)
  dtsM1 = Θ.tsM1 |> CuVector{Int}
  for k ∈ 1:dims.K
    for t ∈ 1:(dims.T-1)

      dfactor!(Θ, Xv, DgX, Gkt=G[k,t]; X, wLAsignX, k, t, dims, dtsM1)
      copyto!(Pv, DgX*XtXinvXtv -
        XXtXinv*(DgX'*X+X'*DgX)*XtXinvXtv +
        XXtXinv*(DgX'*v))
      DG[k,t] = 2.0 * (εt * Pv)

    end
    DgX[:,k] .= T(0.0)
  end

  return DG
end

#try to compute the analytical gradient
function ∇lossag(G::Matrix{T}, Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterControlLevel
  ) where {TM,TV,T}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2), KW=size(Θ.A,1) + size(Xv.W,2),
    Nexpanded=size(Xv.ws,1), expandedinds = collect(1:size(Xv.ws,1)) |> CuVector{Int})
  @unpack v = Xv

  #pre-compute some common values
  #note the reversed intercept location
  X::TM = hcat(initializeÃandfactor(Θ, Xv), Xv.W)
  cuI = TM(I, dims.KW, dims.KW) #need an identity matrix
  XtXinv = cholesky!(X'*X)\cuI
  XXtXinv = X*XtXinv
  XtXinvXtv = XtXinv * (X'*v)
  εt = (X*XtXinvXtv-Xv.v)'

  wLAsignX = signabsμₖ.(Θ.expand.Ã', Θ.expand.LÃ', Xv.ws, Xv.RLws) .* Xv.ws .* abs.(Θ.expand.LÃ')

  #preallocate this carefully high performance
  DgXs::Vector{TM} = [zeros(T, dims.Nexpanded, dims.KW) for i ∈ 1:Threads.nthreads()]
  Pvs = [similar(Xv.v) for i ∈ 1:Threads.nthreads()]

  DG = similar(G)

  dtsM1 = Θ.tsM1
  for k ∈ 1:dims.K
    Threads.@threads for t ∈ 1:(dims.T-1)
      DgX = DgXs[Threads.threadid()]
      Pv = Pvs[Threads.threadid()]

      dfactor!(Θ, Xv, DgX, Gkt=G[k,t]; X, wLAsignX, k, t, dims, dtsM1)
      copyto!(Pv, DgX*XtXinvXtv -
        XXtXinv*(DgX'*X+X'*DgX)*XtXinvXtv +
        XXtXinv*(DgX'*v))
      DG[k,t] = 2.0 * (εt * Pv)

    end
    for DgX ∈ DgXs
      DgX[:,k] .= T(0.0)
    end
  end

  return DG
end

function testag(Θ, Xv::AbstractXYIter{TM,TV,T},
  PredictionType::Val=Val{:level}(),
  GrowthType::Val = Val(PARAM[:itergrowthtype])) where {TM,TV,T}
  CUDA.allowscalar(false)

  @info "Beginning X_g accuracy test. PredictionType=$PredictionType, GrowthType=$GrowthType"

  dims = (K=size(Θ.A,1), T = size(Θ.A,2),
    Nexpanded=size(Xv.ws,1), expandedinds = collect(1:size(Xv.ws,1)) |> CuVector{Int})

  Θtest = deepcopy(Θ)
  Θtest.G .= exp.(Θtest.G)

  local Ãffactor, Ãag
  #first check accuracy of d/dG(X)
  function ffactor(G)
    someones = ones(T, dims.K, 1)
    prodmatG = cumprodbyrow(G)
    Ã = hcat(someones, prodmatG)
    Ãffactor = deepcopy(Ã)
    Ãexpanded = CuMatrix{T}(Ã[Θtest.expand.index.Ã...])
    LÃexpanded = CuMatrix{T}(Ã[Θtest.expand.index.LÃ...])
    factored = genabsμₖ.(Ãexpanded, LÃexpanded, Xv.ws', Xv.Rtws', Xv.RLws')' #0.4ms

    return sum(factored)
  end

  function ∇Xag(G)
    @assert (G|>Matrix) == (Θtest.G|>Matrix)
    X = initializeÃandfactor(Θtest, Xv) |> deepcopy
    DgX = zeros(T, dims.Nexpanded, dims.K) |> TM
    wLAsignX = signabsμₖ.(Θtest.expand.Ã', Θtest.expand.LÃ', Xv.ws, Xv.Rtws, Xv.RLws) .* Xv.ws .* abs.(Θtest.expand.LÃ')
    DG = similar(G)
    Ãag = deepcopy(Θtest.Ã)
    dtsM1 = Θ.tsM1 |> CuVector{Int}
    for k ∈ 1:dims.K
      for t ∈ 1:(dims.T-1)
        dfactor!(Θtest, Xv, DgX, Gkt=G[k,t]; wLAsignX, k, t, dims, X, dtsM1)
        DG[k,t] = sum(DgX)
      end
      DgX[:,k] .= 0.0
    end

    return DG
  end

  ∇X = gradient(ffactor, Θtest.G |> Matrix)
  ∇Xans = ∇X[1]
  ∇agXans = ∇Xag(Θtest.G |> Matrix)
  @assert (Ãag |> Matrix) ≈ (Ãffactor |> Matrix) #check product matrix
  if !(∇agXans ≈ ∇Xans)
    error("!(∇Xag(Gtest) ≈ ∇X(Gtest))!!
      ∇Xag: $(∇agXans)
      ∇X: $(∇Xans)")
  end
  @info "ag X accuracy test passed"

  #acquire gradient functions
  ∇zygote = ∇lossforzygote(Θ, Xv, GrowthType, GradType=Val{:zygote}())
  ∇ag = ∇lossforzygote(Θ, Xv, GrowthType, GradType=Val{:ag}())

  #create the inputs and space for the outputs
  g = Θ.g |> TV |> deepcopy
  ∇gzygote = g |> TV |> deepcopy |> Vector
  ∇gag = g |> TV |> deepcopy |> Vector

  #check the accuracy
  @info "beginning ag accuracy test"

  ∇zygote(∇gzygote, g)
  ∇ag(∇gag, g)
  if !(∇gag ≈ ∇gzygote)
    error("!(∇gag ≈ ∇gzygote)!!
      ∇gag: $∇gag
      ∇gzygote: $∇gzygote")
  end
  @info "ag accuracy test passed"

  #check performance
  if PARAM[:benchgrad]
    print("Zygote performance: ")
    @btime CUDA.@sync $∇zygote($∇gzygote, $g)
    print("Analytic performance: ")
    @btime CUDA.@sync $∇ag($∇gag, $g)
  end

  return nothing
end

##
#this used to be part of VolumePartsIter, now it is here for testing purposes
function volumepartsxsection(Θ::VolumePartsIter)
  @unpack A,G,ts = Θ

  dims = (T=size(A,2), K=size(A,1))
  xsections = Vector{NamedTuple}(undef, dims.T-1)
  for t ∈ 2:dims.T
    #count the number of stocks in each window (t and t-1)
    Nt = sum(ts.==t)

    #repeat t for each stock
    Ainds = repeat([t],Nt)
    LA_inds = repeat([t-1], Nt)
    xsections[t-1] = (A = view(A, :, Ainds),
      LA = view(A, :, LA_inds),
      G = view(G, :, LA_inds))
  end
  xsection::NamedTuple = (
    A = (x->x.A).(xsections),
    LA = (x->x.LA).(xsections),
    G = (x->x.G).(xsections),
  )

  return xsection
end

###These are legacy code routines moved from VolumePartsIter, still sometimes used in testing
#creates a factored object and associated cross-secitons for building up the
#regression rhs
genfactoredLA(wsᵢ, Rtwsᵢ, RLwsᵢ,Gᵢt)= abs(Rtwsᵢ  - RLwsᵢ) .+ abs(Rtwsᵢ  - Gᵢt * wsᵢ)
#factors the |G*LA*w-LA*RLw| expression(pulls out the LA)
function factorLA(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  ) where {TM,TV,T}
#use the pre-allocation
  factored =  genfactoredLA.(Xv.ws, Xv.Rtws, Xv.RLws, Θ.expand.G')

  return factored
end

function unfactorLA(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  factoredLA::TM) where {TM, TV, T}
  absμₖ = factoredLA .* (abs).(Θ.expand.LA)'

  return absμₖ
end

#factors the |G*LA*w-LA*RLw| expression(pulls out the LA) for each crosssection
function xfactorLA(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},xsection::NamedTuple,
  ) where {TM,TV,T}

  @assert !((TM <: CuArray) || TV <: CuArray)
  Θxsection = volumepartsxsection(Θ)
  #use the pre-allocation
  factoredLA::Vector{TM} = deepcopy(xsection.ws)
  Threads.@threads for i ∈ 1:(size(Θ.A,2)-1)
    factoredLA[i] .= genfactoredLA.(xsection.ws[i],xsection.Rtws[i],xsection.RLws[i], Θxsection.G[i]')
  end

  return factoredLA
end

#gpu version of above
#primarily for testing since the xsection is no longer pre-allocated
function xfactorLA(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},xsection::NamedTuple,
  ) where {TM<:CuMatrix,TV<:CuVector,T}

  Θxsection = volumepartsxsection(Θ)

  #form the crosssection
  factoredLA::Vector{TM} = CUDA.deepcopy(xsection.ws)
  for i ∈ 1:(size(Θ.A,2)-1)
    factoredLA[i] .= genfactoredLA.(xsection.ws[i],xsection.Rtws[i],xsection.RLws[i], Θxsection.G[i]')
  end

  return factoredLA
end

function unfactorA(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  factoredA::TM) where {TM, TV, T}
  absμₖ = factoredA .* abs.(Θ.expand.A)'

  return absμₖ
end

#copies relevant data into global variables for testing purposes
function testglobal(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterLevel,
    ) where {TM,TV,T}
  global thetag, Xvg = deepcopy(Θ), deepcopy(Xv)
  global factoredg = factorLA(Θ, Xv) |> Matrix{T}
  global df = DataFrame(factoredg)
  df.ts = (i->"t$i").(thetag.ts)
  df.v = Xvg.v |> Vector{T}

  global v = df.v |> TV
  global X = factoredg |> TM
  return nothing
end

#copies relevant data into global variables for testing purposes
function testglobal(Θ::AbstractVolumePartsIter{TM,TV,T}, Xv::XYIterDemean,
    ) where {TM,TV,T}
  global thetag, Xvg = deepcopy(Θ), deepcopy(Xv)
  global factoredg = factorLA(Θ, Xv) |> Matrix{T}
  global df = DataFrame(factoredg)
  df.ts = (i->"t$i").(thetag.ts)
  df.v = Xvg.v |> Vector{T}
  #f = @formula(0 ~ x1 + x2 + ts + 0)
  global MM = hcat(factoredg, Xv.M.F|> collect |> Matrix{T})
  global F = MM[:,(1+size(Xv.ws,2)):end] |> TM
  global v = df.v |> TV
  global X = factoredg |> TM
  return nothing
end

function testfixed(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
    PredictionType::Val,
    RegressionType::Val=Val{PARAM[:iterregressiontype]}(),
    benchmeasure::Bool = true) where {TM,TV,T}

  CUDA.allowscalar(false)
  dims = (K=size(Θ.A,1), T=size(Θ.A,2))
  factoredg = factorLA(Θ, Xv) |> Matrix{T}

  #F = Xv.M.F |> TM
  v = Xv.v |> TV
  X = factoredg |> TM

  @assert Xv.M*v ≈ Xv.ṽ

  #now check the demeaned predictions by strategy
  absμₖdemean::Matrix{Float64} = similar(Xv.ws)
  D = (ws=Xv.ws, Rtws=Xv.Rtws, RLws=Xv.RLws)
  @unpack ws, Rtws, RLws = dataxsection((ws=Xv.ws, Rtws=Xv.Rtws, RLws=Xv.RLws),Θ.ts,TM,TV,T)
  @unpack xabsμₖ = dataxsection((xabsμₖ=absμₖdemean,),Θ.ts,Matrix{Float64},Vector{Float64},Float64)

  Θxsection=volumepartsxsection(Θ)
  for t ∈ 1:(dims.T-1)
    xabsμₖ[t] .= genabsμₖ.(Θxsection.A[t]', Θxsection.LA[t]', ws[t], Rtws[t], RLws[t])
    xabsμₖ[t] .-= mean(xabsμₖ[t], dims=1) #demeaning step
  end
  absμₖdemean .= reduce(vcat, xabsμₖ)
  absμₖproject::Matrix{Float64} = projectstrategyvolume(Θ, Xv)

  absμₖrmserror = sqrt(mean((absμₖdemean .- absμₖproject).^2))
  if T.(absμₖdemean) ≈ T.(absμₖproject)
    @info "testfixed absμₖ test passed with rmserror=$absμₖrmserror"
  else
    println("absμₖ test failed!!! absμₖdemean[1:10,:]")
    printmln(absμₖdemean[1:10,:])
    println("absμₖproject[1:10,:]")
    printmln(absμₖproject[1:10,:])
    error("absμₖ test failed!!! rmserror=$absμₖrmserror")
  end

  vhatdemean = sum(absμₖdemean, dims=2) |> vec
  vhatproject = sum(absμₖproject, dims=2) |> vec

  vrmserror = sqrt(maximum((vhatdemean .- vhatproject).^2))
  if T.(vhatdemean)  ≈ T.(vhatproject)
    @info "testfixed v test passed with rmserror=$vrmserror"
  else
    error("objective test failed!!! rmserror=$vrmserror")
  end

end

#compares the zygote loss/prediction against the standard loss/prediction
function testpredictionzygote(Θ::AbstractVolumePartsIter, Xv::AbstractXYIter,
  ::Type{T}=PARAM[:itertype];
  RegressionType::Val=Val{PARAM[:iterregressiontypezygote]}()) where T

  #now get the zygote prediction
  G::Matrix = Θ.G |> Matrix
  vhatzygotever = Θ(Xv, RegressionType)
  vhatzygote = Θ(G, Xv, RegressionType)  |> Vector
  (T.(vhatzygotever |> Vector) ≈ T.(vhatzygote)) || error("vhatzygote≠vhattest
    vhattest[1:10]: $(vhatzygote[1:10])
    vhatzygote[1:10]: $(vhatzygotever[1:10])")

  λzygote = lossforzygote(G, Θ, Xv)

  Δver = vhatzygotever .- Xv.ṽ
  λzygotever = sum(Δver .* Δver)

  #!(λzygote ≈ λzygotever) && @warn "λzygote !≈ λzygotever. If zygote fixes cumprod,
  #  replace the zygote cumprod and we can get much better with this."
  (λzygote ≈ λzygotever) || error("λzygote ≠ λzygotever
    λzygote: $λzygote
    λzygotever: $λzygotever
    vhattest[1:10]: $(vhatzygote[1:10])
    vhatzygote[1:10]: $(vhatzygotever[1:10])")
end

#compares the zygote loss/prediction against the standard loss/prediction
function testcudazygote(Θ::VolumePartsIter, Xv::XYIterControl,
  Θcpu::VolumePartsIter, Xvcpu::XYIterControl,
  ::Type{T}=PARAM[:itertype];
  RegressionType::Val=Val{PARAM[:iterregressiontypezygote]}(),
  GrowthType::Val = Val{PARAM[:itergrowthtype]}(),
  benchmeasure::Bool = PARAM[:benchitermeasure]) where T

  ###first compare standard predictions

  vhatcpu = Θcpu(Xvcpu, RegressionType) |> Vector
  vhat = Θ(Xv, RegressionType)  |> Vector
  (T.(vhatcpu) ≈ T.(vhat)) || error("vhat≠vhatcpu (cuda zygote test)
    vhat[1:10]: $(vhat[1:10])
    vhatcpu[1:10]: $(vhatcpu[1:10])")

  λ = iterloss(Θ, Xv, RegressionType)
  λcpu = iterloss(Θcpu, Xvcpu, RegressionType)

  @assert λ ≈ λcpu

  ###compare cpu and cuda zygote prediction
  vhatzygotecpu = Θcpu(Θcpu.G, Xvcpu, RegressionType) |> Vector
  vhatzygote = Θ(Θcpu.G, Xv, RegressionType)  |> Vector
  (T.(vhatzygotecpu |> Vector) ≈ T.(vhatzygote)) || error("vhatzygote≠vhattest (cuda zygote test)
    vhatzygote[1:10]: $(vhatzygote[1:10])
    vhatzygotecpu[1:10]: $(vhatzygotecpu[1:10])")

  λzygote = lossforzygote(Θcpu.G, Θ, Xv,)
  λzygotecpu = lossforzygote(Θcpu.G, Θcpu, Xvcpu,)

  Δcpu = vhatzygotecpu .- Xvcpu.ṽ
  λzygotecpu = sum(Δcpu .* Δcpu)
  @assert (λzygote ≈ λzygotecpu) || error("λzygote ≈/ λzygotecpu!!!
    λzygote: $λzygote ; λzygotecpu: $λzygotecpu")

  ###now compare the gradients
  x = deepcopy(Θcpu.g |> Vector)

  if typeof(GrowthType) <: Val{:log}
    x .= log.(abs.(x))
    Θcpu.g .= abs.(Θcpu.g)
    Θ.g .= abs.(Θ.g)
  end

  outcpu = similar(x)
  out = similar(x)

  ∇Θcpu = ∇lossforzygote(Θcpu, Xvcpu, GrowthType)
  ∇Θ = ∇lossforzygote(Θ, Xv, GrowthType)

  ∇Θcpu(outcpu, x)
  ∇Θ(out, x)

  @assert Θcpu.g ≈ (Θ.g |> Vector)
  @assert Θcpu.G[Θcpu.expand.index.G...] ≈ (Θ.G[Θ.expand.index.G...]  |> Matrix)

  (out ≈ outcpu) || error("(out ≈/ outcpu)!!!  out: $out outcpu: $outcpu")

  if benchmeasure
    print("Loss performance cpu:")
    @btime iterloss($Θ, $Xv, $RegressionType)
    print("Loss performance gpu:")
    @btime iterloss($Θcpu, $Xvcpu, $RegressionType)

    print("Gradient performance cpu:")
    @btime $∇Θcpu($outcpu, $x)
    print("Gradient performance gpu:")
    @btime $∇Θ($out, $x)
  end
end

#compares the zygote loss/prediction against the standard loss/prediction
function testpredictionzygote(Θ::VolumePartsIter, Xv::XYIterDemean,
  ::Type{T}=PARAM[:itertype];
  RegressionType::Val=Val{PARAM[:iterregressiontypezygote]}()) where T

  #generate the standard pre-alloc prediction
  factored = initializeÃandfactor(Θ, Xv)
  Θ.A[:, 1] .= reg(Xv.M * factored, Xv.ṽ,Val{:glm}())
  #@info "Θ.A[:, 1] (testpredictionzygote): $(Θ.A[:, 1])"

  #now get the zygote prediction
  vhatzygotever =  factored * (abs).(Θ.A[:,1])
  vhatzygote = Θ(Θ.G, Xv, Val{:choleskyzygote}())  |> Vector
  (vhatzygotever≈vhatzygote) || error("vhatzygote≠vhattest
    vhattest[1:10]: $(vhatzygote[1:10])
    vhatzygote[1:10]: $(vhatzygotever[1:10])")

  @assert lossforzygote(Θ.G, Θ, Xv) ≈ sum((Xv.ṽ .- (Float64).(vhatzygotever)).^2)
end

#some validation tests for xsection integrity
function validatexsection(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
    p::NamedTuple, xsection::NamedTuple) where {TM,TV,T}
  @unpack dims, tidx, cleanpanel = p

  #@unpack xsection = Xv
  Θxsection=volumepartsxsection(Θ)

  for i ∈ 1:(dims.T-1)
    xsectionsize = size(Θxsection.A[i]')
    (xsectionsize == size(Θxsection.LA[i]')) || error("unpected size of Θxsection.LA[i]'")
    (xsectionsize == size(Θxsection.G[i]')) || error("unpected length of Θxsection.G[i]'")
    (xsectionsize == size(xsection.ws[i]) == size(xsection.RLws[i])) || error(
      "Unexpected array sizing for xsection $i: xsectionsize=$(xsectionsize)
      xsection.ws[i]=$(size(xsection.ws[i]))
      xsection.RLws[i]=$(size(xsection.RLws[i]))")
    @assert xsectionsize[1] == length(xsection.ṽ[i])
    xsectiondate = cleanpanel.date[Θ.ts .== i+1]
    @assert sum(xsectiondate .== xsectiondate[1]) == length(xsectiondate) #all the same date
    @assert i+1 == tidx[xsectiondate[1]]
  end

  return nothing
end

validatexsection(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  p::NamedTuple) where {TM,TV,T} = validatexsection(
    Θ, Xv, p, dataxsection((ws=Xv.ws, RLws=Xv.RLws, ṽ=Xv.ṽ, v=Xv.v), Θ.ts, TM, TV, T))

function testG(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  xsection::NamedTuple) where {TM,TV,T}

  #@unpack xsection = Xv
  Θxsection=volumepartsxsection(Θ)
  #first, generate predicatable test data and run the update fuinction
  dims = (K=size(Θ.A,1), T=size(Θ.A,2))
  testrowA::TM = exp.([isodd(i) ? log(T(1.1)) : log(T(0.9)) for i ∈ 1:dims.T])' |> TM
  testrowA .*= (i->T(-1)^T(i+1)).(1:dims.T)' |> TM
  testrowA .= cumprod(testrowA, dims=2)#, dims=2)
  testrowG::TM = testrowA[:,2:end] ./ testrowA[:,1:(end-1)]

  #exp.(collect(1:dims.T) .* (i->isodd(i)*log(T(1.1)) + iseven(i))).(1:dims.T)'
  # testrowG
  Θ.G .= testrowG
  Θ.A .= hcat(ones(T, dims.K,1) |> TM, cumprodbyrow(Θ.G))
  Θ.G .= Θ.A[:, 2:end] ./ Θ.A[:, 1:(end-1)]

  #test updateAfromGAₜ!
  Θcpy=deepcopy(Θ)
  Θcpy.A[:,25] .= T(25)
  updateAfromGAₜ!(Θcpy, Xv,25)
  ΠG = hcat(ones(T, dims.K,1) |> TM, cumprodbyrow(Θ.G))
  @assert Θcpy.A ≈ ΠG .* T(25) ./  ΠG[:, 25]

  #create test factored flows
  factoredLA = deepcopy(factorLA(Θ, Xv))
  xfactoredLA = deepcopy(xfactorLA(Θ, Xv, xsection))

  @assert factoredLA == reduce(vcat, xfactoredLA)

  #make sure the values of G check out
  testG::TM = reduce(vcat, [testrowG for i ∈1:dims.K])
  testA =hcat(ones(T,dims.K,1)|>TM, cumprodbyrow(testG))
  testG .= testA[:,2:end] ./ testA[:,1:(end-1)]

  (Θ.G |> Matrix! ≈ testG |> Matrix!) || error(
    "Differences found between G and test
    G[:,1:5]: $(Θ.G[:, 1:5]) expected $(testG[:, 1:5])")

  #verify the overall data structure is as expected
  expected = Matrix(abs.(Θ.expand.LA' .* Xv.Rtws .- Θ.expand.LA' .* Xv.RLws) .+
    abs.(Θ.expand.LA' .* Xv.Rtws .- Θ.expand.A' .* Xv.ws))
  @assert expected ≈ Matrix!(factoredLA .* abs.(Θ.expand.LA)')
  @assert expected ≈ Matrix!(unfactorLA(Θ,Xv,factoredLA))
  @assert expected ≈ Matrix!(genabsμₖ.(Θ.expand.A', Θ.expand.LA', Xv.ws, Xv.Rtws, Xv.RLws))

  #verify each crossection is as expected
  expectedprojectt = Vector{Matrix{Float64}}()
  for i ∈ 1:(dims.T-1)
    Nᵢ::Int = size(xsection.ws[i],1)
    expected = Matrix!(abs.(Θxsection.LA[i]' .* xsection.Rtws[i] .- Θxsection.LA[i]' .* xsection.RLws[i]
      ) .+ abs.(Θxsection.LA[i]' .* xsection.Rtws[i] .- Θxsection.A[i]' .* xsection.ws[i]))
    xpredLA = Matrix!(xfactoredLA[i]) .* Matrix!(abs.(Θxsection.LA[i]))'

    @assert (expected ≈ xpredLA) "
      size expected:  $(size(expected))
      size Matrix(xpredLA):  $(size(xpredLA))
      expected[:,1:5]: $(expected[1:5,:])
      predLA[:,1:5]: $(xpredLA[1:5,:])"

    #the below is a simple component test of the projectt methods and
    #genabsμk
    if Xv.xsection.xM ≡ nothing
      expected = Matrix!(abs.(Θxsection.LA[i]' .* Xv.xsection.xRtws[i] .- Θxsection.LA[i]' .* Xv.xsection.xRLws[i]
        ) .+ abs.(Θxsection.LA[i]' .* Xv.xsection.xRtws[i] .- Θxsection.A[i]' .* Xv.xsection.xws[i]))
      @assert expected ≈ Matrix!(projectt(Θxsection.A[i], Θxsection.LA[i],
        Xv.xsection.xws[i], Xv.xsection.xRtws[i], Xv.xsection.xRLws[i], nothing))
      push!(expectedprojectt, expected)
    else
      expected = Matrix!(Xv.xsection.xM[i] * (abs.(Θxsection.LA[i]' .* Xv.xsection.xRtws[i] .- Θxsection.LA[i]' .* Xv.xsection.xRLws[i]
        ) .+ abs.(Θxsection.LA[i]' .* Xv.xsection.xRtws[i] .- Θxsection.A[i]' .* Xv.xsection.xws[i])))
      @assert expected ≈ Matrix!(projectt(Θxsection.A[i], Θxsection.LA[i],
        Xv.xsection.xws[i], Xv.xsection.xRtws[i], Xv.xsection.xRLws[i], Xv.xsection.xM[i]))
      push!(expectedprojectt, expected)
    end
  end

  combinedexpectedproj = reduce(vcat, expectedprojectt)
  @assert combinedexpectedproj ≈ Matrix!(projectbyt(testA, Xv.xsection))
  @info "completed test of G"

  return nothing
end

#verify that standardized predictions are equivelent to the unstandardized
#(after the scaling adjustment)
function teststandardization(panel::AbstractDataFrame,  zs::ZSpec,
  ::Type{TM}, ::Type{TV}, ::Type{T}) where {TM, TV, T}
  @unpack Zms = zs
  ms = zs.originalms

  ZFw::Vector{Symbol} = Fweights(Zms, :Fw)
  ZFRtw::Vector{Symbol} = Fweights(Zms, :FRtw)
  ZFRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  Fw::Vector{Symbol} = Fweights(ms, :Fw)
  FRtw::Vector{Symbol} = Fweights(ms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(ms, :FRLw)

  p = modelparts(panel, T, TV, TM,
    weightdata = (;Fw, FRtw, FRLw),
    Fvol=ms.Fvol, zs=zs,)
  Zp = modelparts(panel, T, TV, TM,
    weightdata = (;Fw=ZFw, FRtw=ZFRtw, FRLw=ZFRLw),
    Fvol=Zms.Fvol, ms=Zms, zs=zs)

  dims = p.dims
  ws::TM = p.dat[:Fw]
  Rtws::TM = p.dat[:FRtw]
  RLws::TM = p.dat[:FRLw]
  v::TV = p.dat[:Fvol]
  ts::Vector{Int} = p.ts
  Nexpanded::Int = length(ts)

  #use the standardized verions
  Xv::XYIterLevel = XYIterLevel(ws,Rtws, RLws, v, ts)
  Θ = VolumePartsIter(dims.T, dims.K, ts, T)

  Zws::TM = Zp.dat[:Fw]
  ZRtws::TM = Zp.dat[:FRtw]
  ZRLws::TM = Zp.dat[:FRLw]
  Zv::TV = Zp.dat[:Fvol]

  #now test the level versions
  ZXv::XYIterLevel = XYIterLevel(Zws,ZRtws,ZRLws, Zv, ts)
  ZΘ = VolumePartsIter(dims.T, dims.K, ts, T)

  #placeholder values for G
  Θ.G .= (T(0.75) .+ T(0.5) .* rand(Uniform(),size(Θ.G))) .* rand((-1,1,1,1), size(Θ.G)) |> TM
  ZΘ.G .= Θ.G

  #now make the predictions
  iterloss(Θ, Xv, Val(:cholesky))
  iterloss(ZΘ, ZXv, Val(:cholesky))

  #sync up A
  updateAfromGAₜ!(Θ,Xv, 1)
  updateAfromGAₜ!(ZΘ,ZXv, 1)

  ZA = ZΘ.A |> deepcopy |> Matrix
  A = Θ.A |> deepcopy |> Matrix

  Fcoefs::Vector{Symbol} = (ξ->ξ.X[:Fw]).(ms.ξs)
  FZcoefs::Vector{Symbol} = (ξ->ξ.X[:Fw]).(Zms.ξs)

  #de-scale ZA
  ZAscale = (s->zs.scale[s]).(FZcoefs)
  σᵥ::Float64 = zs.scale[zs.Zms.Fvol]
  println("mean scaling ratio (pre-descaling): $(mean(A ./ ZA))")
  AfromZA = (ZA .* σᵥ) ./ ZAscale

  if !( AfromZA ≈ A)
    println("σᵥ: $σᵥ")
    println("ZAscale: $ZAscale")

    println("AfromZA: ")
    printmln(AfromZA[:,1:10])
    println("A:")
    printmln(A[:,1:10])
    println("AfromZA/A:")
    printmln(AfromZA[:,1:10] ./ A[:,1:10])
    error("!( AfromZA ≈ A) !!!")
  end
  @info "Passed standardization test"

  return nothing
end

#this fucntion tests the equivelency of level, intercept, demean specs against the control spec
#note we may want to get rid of this or at least put it in a seperate file
#along with all the legacy specs
function testcontrolequivelency(panel::AbstractDataFrame,  zs::ZSpec,
  ::Type{TM}, ::Type{TV}, ::Type{T};) where {TM,TV,T}

  @info "Beginning control equivlency test"

  CUDA.allowscalar(false)
  Fw::Vector{Symbol} = Fweights(zs.Zms, :Fw)
  FRtw::Vector{Symbol} = Fweights(zs.Zms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(zs.Zms, :FRLw)
  p = modelparts(panel, T, TV, TM,
    weightdata = (;Fw, FRtw, FRLw),
    Fvol=zs.Zms.Fvol, zs=zs)

  dims = p.dims
  ws::TM = p.dat[:Fw]
  Rtws::TM = p.dat[:FRtw]
  RLws::TM = p.dat[:FRLw]
  v::TV = p.dat[:Fvol]
  ts::Vector{Int} = p.ts
  Nexpanded::Int = length(ts)

  Xvverlevel::XYIterLevel = XYIterLevel(ws, Rtws, RLws, v, ts)
  Xvverintercept::XYIterLevel = XYIterLevel(ws ,Rtws, RLws, v, ts)
  Xvverdemean::XYIterDemean = XYIterDemean(ws|>Array ,Rtws|>Array, RLws|>Array, v|>Array, ts) #no gpu version here
  Xvtestlevel::XYIterControl= XYIterControl(ws ,Rtws, RLws, v, ts)
  Xvtestintercept::XYIterControl= XYIterControl(ws ,Rtws, RLws, v, ts, W=ones(Nexpanded,1)|>TM)
  Xvtestdemean::XYIterControl= XYIterControl(ws ,Rtws, RLws, v, ts, Wgroup=ones(Nexpanded,1)|>TM)

  #placeholder values for G

  Θverlevel = VolumePartsIter(dims.T, dims.K, ts, T)
  G = (T(0.75) .+ T(0.5) .* rand(Uniform(),size(Θverlevel.G))) |> TM
  Θverlevel.G .= G
  g = vec(log.(G)) |> Array #use the log form since this is what I use in practice

  #form the other prediction objects
  Θverintercept= Θverlevel |> deepcopy
  Θverdemean = VolumePartsIter(dims.T, dims.K, ts, T, Matrix{T}, Vector{T})
  Θverdemean.G .= G |> Array
  Θtestlevel = Θverlevel |> deepcopy
  Θtestintercept = Θverlevel |> deepcopy
  Θtestdemean = Θverlevel|> deepcopy

  #compute losses
  lossverlevel = iterloss(Θverlevel, Xvverlevel, Val(:cholesky))
  lossverintercept = iterloss(Θverintercept, Xvverintercept, Val(:intercept))
  lossverdemean = iterloss(Θverdemean, Xvverdemean, Val(:cholesky))
  losstestlevel = iterloss(Θtestlevel, Xvtestlevel, Val(:cholesky))
  losstestintercept = iterloss(Θtestintercept, Xvtestintercept, Val(:cholesky))
  losstestdemean = iterloss(Θtestdemean, Xvtestdemean, Val(:cholesky))

  #computepredictions
  μverlevel = Θverlevel(Xvverlevel, Val{:cholesky}()) |> Vector
  μverintercept = Θverintercept(Xvverintercept, Val{:intercept}()) |> Vector
  μverdemean = Θverdemean(Xvverdemean, Val{:cholesky}()) |> Vector
  μtestlevel = Θtestlevel(Xvtestlevel, Val{:cholesky}()) |> Vector
  μtestintercept = Θtestintercept(Xvtestintercept, Val{:cholesky}()) |> Vector
  μtestdemean = Θtestdemean(Xvtestdemean, Val{:cholesky}()) |> Vector

  #check the pre-conditioned volume
  @assert (Xvtestdemean.ṽ |> Array )≈ (Xvverdemean.ṽ |> Array)

  #check the predictions
  @assert μverlevel ≈ μtestlevel
  @assert μverintercept ≈ μtestintercept
  #@assert μverdemean ≈ μtestdemean
  (μverdemean ≈ μtestdemean) || error("!(μverdemean ≈ μtestdemean)!!
    μverdemean[1:10] $(μverdemean[1:10])
    μtestdemean[1:10] $(μtestdemean[1:10])")

  #check the aggregate loss
  @assert lossverlevel ≈ losstestlevel
  @assert lossverintercept ≈ losstestintercept
  @assert lossverdemean ≈ losstestdemean

  #acquire gradient functions
  ∇verlevel! = ∇lossforzygote(Θverlevel, Xvverlevel,
    Val(:log), RegressionType=Val(:choleskyzygote))
  ∇verintercept! = ∇lossforzygote(Θverintercept, Xvverintercept,
    Val(:log), RegressionType=Val(:interceptzygote))
  ∇verdemean! = ∇lossforzygote(Θverdemean, Xvverdemean,
    Val(:log), RegressionType=Val(:choleskyzygote))
  ∇testlevel! = ∇lossforzygote(Θtestlevel, Xvtestlevel,
    Val(:log), RegressionType=Val(:choleskyzygote))
  ∇testintercept! = ∇lossforzygote(Θtestintercept, Xvtestintercept,
    Val(:log), RegressionType=Val(:choleskyzygote))
  ∇testdemean! = ∇lossforzygote(Θtestdemean, Xvtestdemean,
    Val(:log), RegressionType=Val(:choleskyzygote))

  #these will hold the gradient output
  gverlevel  = similar(g)
  gverintercept = similar(g)
  gverdemean = similar(g)
  gtestlevel = similar(g)
  gtestintercept = similar(g)
  gtestdemean = similar(g)

  global g = g
  global gverlevel = gverlevel
  global verlevel = ∇verlevel!

  ∇verlevel!(gverlevel, g)
  ∇verintercept!(gverintercept, g)
  ∇verdemean!(gverdemean, g)
  ∇testlevel!(gtestlevel, g)
  ∇testintercept!(gtestintercept, g)
  ∇testdemean!(gtestdemean, g)

  if PARAM[:benchcontrol]
    print("Benchmark ∇verlevel!(gverlevel, g): ")
    @btime CUDA.@sync $∇verlevel!($gverlevel, $g)
    print("Benchmark ∇verintercept!(gverintercept, g): ")
    @btime CUDA.@sync $∇verintercept!($gverintercept, $g)
    print("Benchmark ∇verdemean!(gverdemean, g): ")
    @btime CUDA.@sync $∇verdemean!($gverdemean, $g)
    print("Benchmark ∇testlevel!(gtestlevel, g): ")
    @btime CUDA.@sync $∇testlevel!($gtestlevel, $g)
    print("Benchmark ∇testintercept!(gtestintercept, g): ")
    @btime CUDA.@sync $∇testintercept!($gtestintercept, $g)
    print("Benchmark ∇testdemean!(gtestdemean, g): ")
    @btime CUDA.@sync $∇testdemean!($gtestdemean, $g)
  end
  @assert gverlevel ≈ gtestlevel
  @assert gverintercept ≈ gtestintercept
  @assert gverdemean ≈ gtestdemean

  @info "Test equivelancy completed successfully."
end

####main entry point for testing
function testiter(panel::AbstractDataFrame,  zs::ZSpec,
    ::Type{T}=PARAM[:itertype],
    ::Type{TV}=PARAM[:itergpu] ? CUDA.CuVector{T} : Vector{T},
    ::Type{TM}=PARAM[:itergpu] ? CUDA.CuMatrix{T} : Matrix{T};
    tol::Float64 = PARAM[:testtol]
    ) where {
      T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  @unpack Zms = zs
  #NOTE: disabled as experiment if PARAM[:fluxmodel] ≠ :levellsq
  #  error("PARAM[:fluxmodeltype]≠:levellsq is not coded for use with Jacobian")
  #end

  testcumprodbyrow()
  choltest()
  testXYIterControl()
  teststandardization(panel, zs, TM, TV, T)

  global gtestcontrolequivelency = ()->testcontrolequivelency(panel, zs, TM, TV, T)
  testcontrolequivelency(panel, zs, TM, TV, T)

  Fw::Vector{Symbol} = Fweights(Zms, :Fw)
  FRtw::Vector{Symbol} = Fweights(Zms, :FRtw)
  FRLw::Vector{Symbol} = Fweights(Zms, :FRLw)
  essentialfields = [Fw; FRtw; FRLw; Zms.Fvol; ]

  #form the data into a more usable format and expand the arrays
  p = modelparts(panel, T, TV, TM,
    weightdata = (;Fw, FRtw, FRLw),
    Fvol=Zms.Fvol, zs=zs)

  @warn "test gpu commented out, should work after dependencies of packages are sorted"
  #testgpu(p)

  dims = p.dims
  ws::TM = p.dat[:Fw]
  Rtws::TM = p.dat[:FRtw]
  RLws::TM = p.dat[:FRLw]
  v::TV = p.dat[:Fvol]
  ts::Vector{Int} = p.ts
  Nexpanded::Int = length(ts)
  #this creates the factored objects- which are accurate given that initially G=1
  #absΔ(wsᵢ::CuArray,RLwsᵢ)=CUDA.abs(wsᵢ-RLwsᵢ)
  absΔ(wsᵢ,Rtwsᵢ,RLwsᵢ)=abs(Rtwsᵢ-RLwsᵢ)

  #now test the level versions
  Xvlevel = XYIterControl(ws , Rtws,RLws, v, ts, W=ones(T, Nexpanded, 1) |> TM) #includes intercept
  Xvdemean = XYIterControl(ws , Rtws,RLws, v, ts, Wgroup=ones(T, Nexpanded, 1) |> TM)

  Θ = VolumePartsIter(dims.T, dims.K, ts, T)

  Dlevel = (ws=Xvlevel.ws, Rtws=Xvlevel.Rtws, RLws=Xvlevel.RLws,)
  xsectionlevel = dataxsection(Dlevel,Θ.ts,TM, TV, T)

  testG(Θ, Xvlevel, xsectionlevel)
  @assert Θ.G |> Matrix! ≈ Θ.A[:,2:end] ./ Θ.A[:,1:(end-1)] |> Matrix!

  #WARNING WARNING WARNING uncomment below!!!
  testpredictionzygote(deepcopy(Θ),Xvlevel)
  testiterloss(deepcopy(Θ),Xvlevel)

  Ddemean = (ws=Xvdemean.ws, Rtws=Xvdemean.Rtws, RLws=Xvdemean.RLws, ṽ = Xvdemean.ṽ,)
  xsectiondemean = dataxsection(Ddemean,Θ.ts,TM, TV, T)
  validatexsection(Θ, Xvdemean, p, xsectiondemean)

  testG(Θ, Xvdemean, xsectiondemean)
  @assert Θ.G |> Matrix! ≈ Θ.A[:,2:end] ./ Θ.A[:,1:(end-1)] |> Matrix!

  testXv(Θ, Xvlevel, Xvdemean)

  testpredictionzygote(deepcopy(Θ),Xvdemean,)
  testiterloss(deepcopy(Θ),Xvdemean,)

  if (TM<:CuArray)
    pcpu = modelparts(panel, T, Vector{T}, Matrix{T},
      weightdata = (;Fw, FRtw, FRLw),
      Fvol=Zms.Fvol,
      zs=zs)

    wscpu::Matrix{T} = pcpu.dat[:Fw]
    Rtwscpu::Matrix{T} = pcpu.dat[:FRtw]
    RLwscpu::Matrix{T} = pcpu.dat[:FRLw]
    vcpu::Vector{T} = pcpu.dat[:Fvol]

    Xvlevelcpu::XYIterControl = XYIterControl(wscpu, Rtwscpu, RLwscpu, vcpu, ts, W= Xvlevel.W|>Matrix)
    Θcpu = VolumePartsIter(dims.T, dims.K, ts, T, Matrix{T}, Vector{T})
    Θcpu.G .= Θ.G |> Matrix
    Θcpu.A .= Θ.A |> Matrix

    Xvdemeancpu::XYIterControl = XYIterControl(wscpu, Rtwscpu, RLwscpu, vcpu, ts, Wgroup= ones(T, Nexpanded, 1))

    @info "beginning tests on Xvlevel gpu v cpu"
    testcudazygote(deepcopy(Θ),Xvlevel,Θcpu|>deepcopy, Xvlevelcpu)

    @info "beginning tests on Xvdemean gpu v cpu"
    testcudazygote(deepcopy(Θ),Xvdemean,Θcpu|>deepcopy, Xvdemeancpu)
    testag(Θ, Xvlevel)
  end

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using Sched

sched = Scheduler()

# Time as Float64
# global _time = time

# Time as DateTime
global _time = UTCDateTimeFunc

function print_time_noparam()
    println("From print_time_noparam $(_time())")
end

function print_time_args(x)
    println("From print_time_args $(_time()) $x")
end

function print_time_kwargs(; a="default")
    println("From print_time_kwargs $(_time()) $a")
end

function print_some_times()
    println(_time())
    enter(sched, Dates.Second(10), 1, print_time_noparam)
    enter(sched, Dates.Second(5), 2, print_time_args, ("positional, argument"))
    enter(sched, Dates.Second(5), 1, print_time_kwargs; Dict(:a=>"keyword")...)
    run(sched)
    println(_time())
end

print_some_times()


text:
# Julia GARCH package
# Copyright 2013 Andrey Kolev
# Distributed under MIT license (see LICENSE.md)

"Generalized Autoregressive Conditional Heteroskedastic (GARCH) models for Julia."
module GARCH

using NLopt, Distributions, Printf, LinearAlgebra, SpecialFunctions, TimeSeries

export garchFit, predict

include("stattests.jl")

"Fitted GARCH model object."
struct GarchFit
    data::Vector
    params::Vector
    llh::Float64
    status::Symbol
    converged::Bool
    sigma::Vector
    hessian::Array{Float64,2}
    cvar::Array{Float64,2}
    secoef::Vector
    tval::Vector
end

function Base.show(io::IO ,fit::GarchFit)
    pnorm(x) = 0.5 * (1 + erf(x / sqrt(2)))
    prt(x) = 2 * (1 - pnorm(abs(x)))
    jbstat, jbp = jbtest(fit.data./fit.sigma)

    @printf io "Fitted garch model \n"
    @printf io " * Coefficient(s):    %-15s%-15s%-15s\n" "ω" "α" "β"
    @printf io "%-22s%-15.5g%-15.5g%-15.5g\n" "" fit.params[1] fit.params[2] fit.params[3]
    @printf io " * Log Likelihood: %.5g\n" fit.llh
    @printf io " * Converged: %s\n" fit.converged
    @printf io " * Solver status: %s\n\n" fit.status
    @printf io " * Standardised Residuals Tests:\n"
    @printf io "   %-26s%-15s%-15s\n" "" "Statistic" "p-Value"
    @printf io "   %-21s%-5s%-15.5g%-15.5g\n\n" "Jarque-Bera Test" "χ²" jbstat jbp
    @printf io " * Error Analysis:\n"
    @printf io "   %-7s%-15s%-15s%-15s%-15s\n" "" "Estimate" "Std.Error" "t value" "Pr(>|t|)"
    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "ω" fit.params[1] fit.secoef[1] fit.tval[1] prt(fit.tval[1])
    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "α" fit.params[2] fit.secoef[2] fit.tval[2] prt(fit.tval[2])
    @printf io "   %-7s%-15.5g%-15.5g%-15.5g%-15.5g\n" "β" fit.params[3] fit.secoef[3] fit.tval[3] prt(fit.tval[3])
end

"Estimate Hessian using central difference approximation."
function cdHessian(params, f)
    eps = 1e-4 * params
    n = length(params)
    H = zeros(n, n)
    function step(x, i1, i2, d1, d2)
        xc = copy(x)
        xc[i1] += d1
        xc[i2] += d2
        f(xc)
    end
    for i in 1:n
        for j in 1:n
            H[i,j] = (step(params, i, j, eps[i], eps[j]) -
                      step(params, i, j, eps[i], -eps[j]) -
                      step(params, i, j, -eps[i], eps[j]) +
                      step(params, i, j, -eps[i], -eps[j])) / (4.0*eps[i]*eps[j])
        end
    end
    H
end

"Simulate GARCH process."
function garchSim(ɛ²::Vector, ω, α, β)
    h = similar(ɛ²)
    h[1] = mean(ɛ²)
    for i = 2:length(ɛ²)
        h[i] = ω + α*ɛ²[i-1] + β*h[i-1]
    end
    h
end

"Normal GARCH log likelihood function."
function garchLLH(y::Vector, params::Vector)
    ɛ² = y.^2
    T = length(y)
    h = garchSim(ɛ², params...)
    -0.5*(T-1)*log(2π) - 0.5*sum(log.(h) + (y./sqrt.(h)).^2)
end

"""
    predict(fit::GarchFit, n::Integer=1)

Make n-step prediction using fitted object returned by garchFit (default step=1).

# Arguments
* `fit::GarchFit` : fitted model object returned by garchFit.
* `n::Integer` : the number of time-steps to be forecasted, by default 1 (returns scalar for n=1 and array for n>1).

# Examples

```
fit = garchFit(ret)
predict(fit, n=2)
```

"""
function predict(fit::GarchFit, n::Integer=1)
    if n < 1
        throw(ArgumentError("n shoud be >= 1 !"))
    end
    ω, α, β = fit.params
    y = fit.data
    ɛ² = y.^2
    h = garchSim(ɛ², ω, α, β)
    pred = ω + α*ɛ²[end] + β*h[end]
    if n == 1
        return sqrt(pred)
    end
    pred = [pred]
    for i in 2:n
        push!(pred, ω + (α + β)*pred[end])
    end
    sqrt.(pred)
end

"""
    garchFit(y::Vector)

Estimate parameters of the univariate normal GARCH process.

# Arguments
* `y::Vector`: univariate time-series array

# Examples

```
filename = Pkg.dir("GARCH", "test", "data", "price.csv")
price = Array{Float64}(readdlm(filename, ',')[:,2])
ret = diff(log.(price))
ret = ret - mean(ret)
fit = garchFit(ret)
```

"""
function garchFit(y::Vector)
    ɛ² = y.^2
    T = length(y)
    h = zeros(T)

    #opt = Opt(Sys.isapple() ? (:LN_PRAXIS) : (:LN_SBPLX), 3)  # LN_SBPLX has a problem on mac currently
    opt = Opt(:LN_SBPLX, 3)
    lower_bounds!(opt, [1e-10, 0.0, 0.0])
    upper_bounds!(opt, [1, 0.3, 0.99])
    min_objective!(opt, (params, grad) -> -garchLLH(y, params))
    (minf, minx, ret) = optimize(opt, [1e-5, 0.09, 0.89])
    h = garchSim(ɛ², minx...)

    converged = minx[1] > 0 && all(minx[2:3] .>= 0) && sum(minx[2:3]) < 1.0
    H = cdHessian(minx, x -> garchLLH(y, x))
    cvar = -inv(H)
    secoef = sqrt.(diag(cvar))
    tval = minx ./ secoef
    GarchFit(y, minx, -minf, ret, converged, sqrt.(h), H, cvar, secoef, tval)
end

garchFit(y::TimeArray) = garchFit(values(y))

end  #module


text:
module Makie

using AbstractPlotting, GLMakie
import FileIO
using GLMakie
using GLMakie: assetpath, loadasset

for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
    @eval export $(name)
end

function logo()
    FileIO.load(joinpath(@__DIR__, "..", "assets", "logo.png"))
end

end


text:
module Topologies

import ..Geometry
import ..Domains: Domains, coordinate_type
import ..Meshes: Meshes, domain

"""
   AbstractTopology

Subtypes of `AbstractHorizontalTopology` define connectiveness of a
mesh in the horizontal domain.

# Interfaces

- [`nlocalelems`](@ref)
- [`vertex_coordinates`](@ref)
- [`interior_faces`](@ref)
- [`vertices`](@ref)
- [`boundary_tags`](@ref)
- [`boundary_tag`](@ref)
- [`boundary_faces`](@ref)

"""
abstract type AbstractTopology end

coordinate_type(topology::AbstractTopology) = coordinate_type(domain(topology))

"""
    nlocalelems(topology)

The number of local elements in `topology`.
"""
function nlocalelems end

"""
    (c1,c2,c3,c4) = vertex_coordinates(topology, elem)

The coordinates of the 4 vertices of element `elem`.
"""
function vertex_coordinates end

"""
    (opelem, opface, reversed) = opposing_face(topology, elem, face)

The opposing face of face number `face` of element `elem` in `topology`.

- `opelem` is the opposing element number, 0 for a boundary, negative for a ghost element
- `opface` is the opposite face number, or boundary face number if a boundary
- `reversed` indicates whether the opposing face has the opposite orientation.
"""
function opposing_face end

"""
    i,j = face_node_index(face, Nq, q, reversed=false)

The node indices of the `q`th node on face `face`, where `Nq` is the number of
face nodes in each direction.
"""
@inline function face_node_index(face, Nq, q, reversed = false)
    if reversed
        q = Nq - q + 1
    end
    if face == 1
        return q, 1
    elseif face == 2
        return Nq, q
    elseif face == 3
        return Nq - q + 1, Nq
    else
        return 1, Nq - q + 1
    end
end

"""
    i,j = vertex_node_index(vertex_num, Nq)

The node indices of `vertex_num`, where `Nq` is the number of face nodes in
each direction.
"""
function vertex_node_index(vertex_num, Nq)
    if vertex_num == 1
        return 1, 1
    elseif vertex_num == 2
        return Nq, 1
    elseif vertex_num == 3
        return Nq, Nq
    else
        return 1, Nq
    end
end

"""
    interior_faces(topology::AbstractTopology)

An iterator over the interior faces of `topology`. Each element of the iterator
is a 5-tuple the form

    (elem1, face1, elem2, face2, reversed)

where `elemX, faceX` are the element and face numbers, and `reversed` indicates
whether they have opposing orientations.
"""
function interior_faces(topology)
    InteriorFaceIterator(topology)
end

struct InteriorFaceIterator{T <: AbstractTopology}
    topology::T
end

"""
    boundary_tags(topology)

A `Tuple` or `NamedTuple` of the boundary tags of the topology. A boundary tag
is an integer that uniquely identifies a boundary.
"""
function boundary_tags end

"""
    boundary_tag(topology, name::Symbol)

The boundary tag of the topology for boundary name `name`. A boundary tag
is an integer that uniquely identifies a boundary.
"""
function boundary_tag end

"""
    boundary_faces(topology, boundarytag)

An iterator over the faces of `topology` which face the boundary with tag
`boundarytag`. Each element of the iterator is an `(elem, face)` pair.
"""
function boundary_faces end

"""
    vertices(topology)

An iterator over the unique (shared) vertices of the topology `topology`.
Each vertex is an iterator over `(element, vertex_number)` pairs.
"""
function vertices(topology)
    VertexIterator(topology)
end
struct VertexIterator{T <: AbstractTopology}
    topology::T
end
struct Vertex{T <: AbstractTopology, V}
    topology::T
    num::V
end
Base.eltype(::Type{<:Vertex}) = Tuple{Int, Int}

include("interval.jl")
include("topology2d.jl")

# deprecate
@deprecate boundaries(topology::AbstractTopology) boundary_tags(topology)
@deprecate GridTopology(mesh) Topology2D(mesh)
@deprecate Topology2D(mesh) Topology2D(mesh)

end # module


text:
# Generate the Frequency Array of a String

const NUCLEOTIDE_ORDER = Dict{Char, Int64}('A' => 0, 'C' => 1, 'G' => 2, 'T' => 3)

pattern_to_number(pattern) = foldl((acc, nucleotide) -> 4 * acc + NUCLEOTIDE_ORDER[nucleotide], pattern; init=0)

function computing_frequencies(text, k)
    frequencies = zeros(Int32, 4 ^ k)
    for ix in 1:(length(text)-k+1)
        pattern = text[ix:ix+k-1]
        code = pattern_to_number(pattern)
        frequencies[code+1] += 1
    end
    frequencies
end

function main()
    dna = readline()
    k = parse(Int64, readline())
    result = computing_frequencies(dna, k)
    println(join(result, " "))
end

main()


text:
# ============================================================
# Theories
# ============================================================

export prim_conserve, conserve_prim, mixture_prim_conserve, mixture_conserve_prim
export advection_flux, burgers_flux, euler_flux, euler_jacobi, em_coefficients
export gauss_moments, mixture_gauss_moments, discrete_moments
export moments_conserve, diatomic_moments_conserve, mixture_moments_conserve
export pdf_slope, mixture_pdf_slope, moments_conserve_slope, mixture_moments_conserve_slope
export flux_conserve!, pressure, stress, heat_flux
export maxwellian, energy_maxwellian, maxwellian!, mixture_maxwellian, mixture_maxwellian!
export shakhov, shakhov!, rykov!
export reduce_distribution, full_distribution
export ref_vhs_vis, vhs_collision_time, rykov_zr
export aap_hs_collision_time, aap_hs_prim, aap_hs_diffeq!
export shift_pdf!
export hs_boltz_kn, kernel_mode, fsm_kernel, boltzmann_fft, boltzmann_fft!
export boltzmann_ode!, boltzmann_nuode!, bgk_ode!
export chapman_enskog
export heat_capacity_ratio, sound_speed

include("theory_continuum.jl")
include("theory_plasma.jl")
include("theory_maxwellian.jl")
include("theory_shakhov.jl")
include("theory_diatomic.jl")
include("theory_atom.jl")
include("theory_pdf.jl")
include("theory_fsm.jl")
include("theory_moments_pure.jl")
include("theory_moments_mixture.jl")
include("theory_thermo.jl")
include("theory_particle.jl")


text:
## Preimage refinement by breadth first search
## ===========================================

immutable BFSPartition <: PartitionAlgorithm end

# Preimage of Y under F, unioned with X
function pre_partition{D <: Domain}(
    Y::RandVar{Bool},
    init_box::D,
    ::Type{BFSPartition};
    precision::Float64 = default_precision(),
    dontstop::Function = loose_bounds,
    args...)

  under = Deque{D}()     # Partition of under approximation of Y-1({true})
  rest = Deque{D}()      # (Partition of over approximation of Y-1({true})) \ under
  push!(under, init_box)

  i = 0
  while dontstop(under, rest, i) && !isempty(under)
    lens(:partition_loop, i, under, rest)
    box = shift!(under)
    image = call(Y, box)

    # If all of the box is within preimage keep it
    if isequal(image,t) || issmall(box, precision)
      push!(rest, box)

    # Otherwise split it into disjoint subsets and repeat for each part
    elseif isequal(image,tf)
      for child in mid_split(box)
        push!(under, child)
      end
    end
    i += 1
  end

  ApproxPartition{D}(collect(rest), collect(under))
end

## Stop Functons
## =============

memory_left(under, rest, i; limit::Int = 1000) = length(under) + length(rest) < limit
iters_left(under, rest, i; limit::Int = 1000) = i < limit

"returns true (i.e. continue) when lower and upper probability bounds are far"
function loose_bounds(under, rest, i; delta = 0.05, do_every_i::Int = 100)
  # For speed check only every do_every_i-th iteration
  if i % do_every_i == 0
    partition = ApproxPartition(collect(rest), collect(under))
    prob_bounds = measure(partition)
    return measure(prob_bounds) > delta
  else
    return true
  end
end

"Stops when precision is reached or time or memory constraints exhausted"
loose_bounds_limited(under, rest, i) =
  iters_left(under,rest,i) && memory_left(under, rest, i) && loose_bounds(under, rest, i)


text:
const ρw = 1000.0 # Density of water, kg/m³
const golovin_b = (1500.0)*1e-3 # Golovin collision coefficient, input cm³/g/s to -> m³/kg/s

"""
    mass_from_r(r; ρ=ρw)

Compute mass (kg) of a droplet with given radius (m) and density (kg/m³)
"""
mass_from_r(r; ρ=ρw) = (4*π/3)*ρ*(r^3)

"""
    r_from_mass(x; ρ=ρw)

Compute radius (m) of a droplet with given mass (kg) and density (kg/m³)
"""
r_from_mass(x; ρ=ρw) = (3*x/4/π/ρ)^(1/3)

# Initial cloud droplet distribution

## Collision Kernels

"""
    golovin_kernel(xᵢ, xⱼ)

Compute the Golovin collision kernel for droplets of two given masses (in kg)
"""
golovin_kernel(xᵢ, xⱼ) = golovin_b * (xᵢ + xⱼ)

function _interior_hydro_kernel(E_coal, E_coll, r_sum, tv_diff)
    E_coal * E_coll * π * r_sum*r_sum * abs(tv_diff)
end

"""
    hydrodynamic_kernel(xᵢ, xⱼ)

Compute the hydrodynamic collision kernel for droplets of the two given masses
(in kg) assuming unity collision and coalescence efficiency.
"""
function hydrodynamic_kernel(xᵢ, xⱼ)
    tvᵢ = terminal_v(xᵢ)
    tvⱼ = terminal_v(xⱼ)
    tv_diff = tvⱼ - tvᵢ
    r_sum = r_from_mass(xᵢ) + r_from_mass(xⱼ)
    _interior_hydro_kernel(1.0, 1.0, r_sum, tv_diff)
end

"""
    long_kernel(xᵢ, xⱼ)

Compute the hydrodynamic collision kernel with parameterization of collision
efficiency from Long (1974) for two droplets of given masses (in kg)
"""
function long_kernel(xᵢ, xⱼ)
    rᵢ = r_from_mass(xᵢ)
    rⱼ = r_from_mass(xⱼ)
    tvᵢ = terminal_v(xᵢ)
    tvⱼ = terminal_v(xⱼ)

    tv_diff = tvⱼ - tvᵢ
    r_sum = rᵢ + rⱼ

    r_small = min(rᵢ, rⱼ) * 1e6 # m -> micron
    r_large = max(rᵢ, rⱼ) * 1e6

    # Collection efficiency cut-off in limit of very large drops
    if r_large >= 50
        E_coll = 1.
    else
        E_coll = 4.5e-4 * r_large^2 * (1.0 - 3.0/(max(3.0, r_small) + 1e-2))
    end

    # Limit collection efficiency to 0 <= E_coll <= 1.0
    E_coll = min(E_coll, 1.)
    E_coll = max(0., E_coll)

    _interior_hydro_kernel(1.0, E_coll, r_sum, tv_diff)
end

"""
    terminal_v(x)

Compute terminal velocity of a droplet of given mass (in kg) following the
parameterization of Beard (1976)
"""
function terminal_v(x)
    r = r_from_mass(x)
    d = 2 * r * 1e6 # diameter, m -> μm
    x = x * 1e3 # mass, kg -> g

    if d ≤ 134.43
        α = 4.5795e5
        x_to_beta = x^(2 / 3)
    elseif 134.43 < d && d ≤ 1511.64
        α = 4962.0
        x_to_beta = x^(1 / 3)
    elseif 1511.64 < d && d ≤ 3477.84
        α = 1732.0
        x_to_beta = x^(1 / 6)
    else
        α = 917.0
        x_to_beta = 1.0
    end

    tv = 1e-2 * α * x_to_beta # cm/s -> m/s
end

"""
    kernels(x::AbstractArray{FT, 1}, kernel::Symbol) where {FT <: AbstractFloat}

Pre-compute and cache pair-wise collision kernels for a 1D binned discretization
of mass space.
"""
function kernels(x::AbstractArray{FT, 1}, kernel::Symbol) where {FT <: AbstractFloat}

  m = length(x)
  ck = zeros(FT, m, m)
  cck = zeros(FT, m, m)

  # Compute collision kernel for all potential bin interactions
  for j ∈ 1:m
    for i ∈ 1:j

      if kernel == :golovin
        cck[j, i] = golovin_kernel(x[i], x[j])
      elseif kernel == :hydro
        cck[j, i] = hydrodynamic_kernel(x[i], x[j])
      elseif kernel == :long
        cck[j, i] = long_kernel(x[i], x[j])
      else # default to Golovin
        cck[j, i] = golovin_kernel(x[i], x[j])
      end # kernel

      # Copy over diagonal
      cck[i, j] = cck[j, i]
    end # i
  end # j

  # cache 2d interpolation on kernel vals
  for i ∈ 1:m
    for j ∈ 1:m
      jm = max(j - 1, 1)
      im = max(i - 1, 1)
      jp = min(j + 1, m)
      ip = min(i + 1, m)
      ck[i, j] = 0.125 * (
        cck[i, jm] + cck[im, j] + cck[ip, j] + cck[i, jp]
      )  + 0.5 * cck[i, j]
      if i == j
        ck[i, j] = 0.5 * ck[i, j]
      end
    end
  end

  return ck

end


text:
"""
    $TYPEDEF

$TYPEDFIELDS
"""
@with_kw mutable struct History
    # Number of completed training epochs
    epochs::Int = 0
    # total number of completed steps
    nsteps::Int = 0
    # number of completed steps in current epoch
    nstepsepoch::Int = 0
    # number of seen samples during training, usually `nsteps * batchsize`
    nsamples::Int = 0
    # stores metrics of epochs grouped by phase
    epochmetrics::DefaultDict =
        DefaultDict{AbstractFittingPhase, MVHistory}(() -> MVHistory())
    # stores values of all metrics for every step
    stepmetrics::MVHistory = MVHistory()
    # stores values of all hyperparameters for every step
    stephyperparams::MVHistory = MVHistory()
end

"""
    Recorder(epoch, step, epochstats, stepstats)
"""
struct Recorder <: Callback end

stateaccess(::Recorder) = (cbstate = (history = Write()), batch = Read(), callbacks = Read())
runafter(::Recorder) = (AbstractMetric,)

function on(::EpochBegin, ::AbstractTrainingPhase, recorder::Recorder, learner)
    if isnothing(get(learner.cbstate, :history, nothing))
        learner.cbstate[:history] = History()
    end
    learner.cbstate[:history].nstepsepoch = 0
end

function on(::BatchEnd, phase::AbstractTrainingPhase, recorder::Recorder, learner)
    history = learner.cbstate[:history]

    # increment counters
    history.nsteps += 1
    history.nstepsepoch += 1
    history.nsamples += size(learner.batch.xs)[end]

    # save metrics
    for metric in getmetrics(learner.callbacks)
        push!(history.stepmetrics, Symbol(metric), history.nsteps, stepvalue(metric))
    end

    # TODO: refactor together with hyperparameters
    # save hyper parameters
    # for (OptimParam, value) in getoptimparams(learner.opt)
    #     push!(history.stephyperparams, Symbol(OptimParam.name), history.nsteps, value)
    # end
end

function on(::EpochEnd, phase::AbstractFittingPhase, recorder::Recorder, learner)
    history = learner.cbstate[:history]

    if phase isa AbstractTrainingPhase
        history.epochs += 1
    end

    for metric in getmetrics(learner.callbacks)
        push!(history.epochmetrics[phase], Symbol(string(metric)), history.epochs, epochvalue(metric))
    end
end


text:
if Sys.iswindows()
    cd("C:\\Users\\Federico\\Documents\\GitHub\\pysical_locomotion\\analysis\\behavior")
else
    cd("/Users/federicoclaudi/Documents/Github/LocomotionControl/analysis/behavior")
end

import Pkg;
Pkg.activate(".")

using Plots, Statistics, Colors, StatsBase
import KernelDensity: kde
using MultiKDE
using DataFrames: DataFrame
using CSV
using Glob, NaturalSort
using Interpolations

using Term
import Term: install_term_logger
import Term: install_stacktrace
install_term_logger()
install_stacktrace()

import MyterialColors: salmon, green_dark, grey, grey_dark, grey_darker, black, blue_grey_darker, blue_dark

import jcontrol.comparisons: ComparisonPoints, ComparisonPoint
import jcontrol: FULLTRACK, Solution
import jcontrol.io: load_cached_trials
import jcontrol.trial: Trial
using jcontrol.visuals
import jcontrol: interpolate_wrt_to

"""
Useful things for behavior analysis. Mostly Used in Thesis/Chpt3 analyses.
"""

if Sys.iswindows()
    PLOTS_FOLDER = "D:\\Dropbox (UCL)\\Rotation_vte\\Writings\\THESIS\\Chpt3\\Plots"
else
    PLOTS_FOLDER = "/Users/federicoclaudi/Dropbox (UCL)/Rotation_vte/Writings/THESIS/Chpt3/Plots"
end

# ---------------------------------------------------------------------------- #
#                                    CURVES                                    #
# ---------------------------------------------------------------------------- #

"""
Stores information about a curve's position in the track
"""
struct Curve
    s0::Float64    # position of "start"
    s::Float64     # position of apex
    sf::Float64    # position of end
    maxdur::Float64
    name::String
    tpos::Float64   # where the angular velocity should be at 0 after previous turn
    direction::Int     # 1 for left turns, -1 for right turns

    cp_s0::ComparisonPoint
    cp_s::ComparisonPoint
    cp_sf::ComparisonPoint
end

function Curve(;
    s0::Float64,
    s::Float64 = 0.0,
    sf::Float64 = 0.0,
    maxdur::Float64 = 0.0,
    name::String = "",
    tpos::Float64 = 0.0,
    direction::Int = 0  ,
    )

    # get comparison points
    cpt = ComparisonPoints(FULLTRACK, [s0, s, sf])

    return Curve(s0, s, sf, maxdur, name, tpos, direction, cpt.points...)
end

curves = [
    Curve(;
        name="first", s0 = 0.0, s=32.0, sf=40.0, tpos=20.0, maxdur=1.5, direction=-1
    ),
    Curve(;
        name="second", s0 = 45.0, s=84.0, sf=96.0, tpos=60.0, maxdur=1.5, direction=-1
    ),
    Curve(;
        name="third", s0 = 98.0, s=135.0, sf=144.0, tpos=115.0, maxdur=1.5, direction=1
    ),
    Curve(;
        name="fourth", s0 = 150.0, s=187.0, sf=205.0, tpos=165.0, maxdur=1.5, direction=1
    ),
]

curves_colors  = range(HSL(360, .6, .7), stop=HSL(200, .5, .7), length=length(curves))

"""
Get the frames in a trial in which the mouse is between s0 and s
to get when in a `Curve`
"""
function get_curve_indices(trial, s0, s1)
    start = findlast(trial.s .< s0)
    start = isnothing(start) ? 1 : start
    start += 2

    stop = findfirst(trial.s[start:end] .>= s1) + start
    stop = isnothing(stop) ? length(trial.s) : stop
    stop -= 2

    idxs = repeat([false], length(trial.s))
    idxs[start:stop] .= true
    idx = findfirst(idxs)
    return idxs, idx, start, stop
end

# ---------------------------------------------------------------------------- #
#                                    TRIALS                                    #
# ---------------------------------------------------------------------------- #
all_trials = load_cached_trials(; keep_n = nothing);
@info "Loaded $(length(all_trials)) trials"

# ---------------------------------------------------------------------------- #
#                              TORTUOSITY ANALYSIS                             #
# ---------------------------------------------------------------------------- #
function pathlength(x, y)
    return sum(
        sqrt.(diff(x).^2 .+ diff(y).^2)
    )
end

function linelength(x, y)
    return sqrt((x[end]-x[1])^2 + (y[end]-y[1])^2)
end

"""
Compute tortuosity along an XY, by taking the trace
in Δt intervals and computing the ratio of the path length
vs the line length between the start and the end.

Assumes the trace is at 60fps
"""
function get_max_tortuosity(x, y; Δt=.25)
    T = 0:(1/60):(length(x)/60 - Δt)
    if length(T) == 0
        return 0
    end
    snapshots = T[1]:Δt:(T[end])

    tortuosity = []
    for tval in snapshots
        tf = tval + Δt

        i0 = findfirst(T .>= tval)
        i1 = max(length(x), findlast(T .<= tf))

        _x = x[i0:i1]
        _y = y[i0:i1]

        push!(tortuosity, pathlength(_x, _y) / linelength(_x, _y))
    end
    return max(tortuosity...)
end

# get the tortuosity at each curve for each trial
tortuosity = []
for (i, curve) in enumerate(curves)
    _tortuosity = []
    for trial in all_trials
        _, _, start, stop = get_curve_indices(trial, curve.s0, curve.sf)

        x = trial.x[start:stop]
        y = trial.y[start:stop]

        push!(_tortuosity, get_max_tortuosity(x, y; Δt=.1))
    end
    push!(tortuosity, _tortuosity)
end

# get which trials have low torosity
clean_trials = Dict{Any, Vector}(1=>[], 2=>[], 3=>[], 4=>[])
for (i, curve) in enumerate(curves)
    th = percentile(tortuosity[i], 97.5)
    for t in tortuosity[i]
        push!(clean_trials[i], t < th)
    end

end

# mark trials that are clean at each curve
clean_trials["wholetrial"] = .*(values(clean_trials)...)

tot = length(all_trials)
discarded = length(filter(i->clean_trials["wholetrial"][i] == 0, 1:tot))
@info "After tortuosity analysis, discarded $(round(discarded/tot * 100; digits=3))% of trials | $(tot - discarded) trials left"

# filter trials
trials = [t for (i,t) in enumerate(all_trials) if clean_trials["wholetrial"][i]]
# trials = all_trials
high_torosity_trials = [t for (i,t) in enumerate(all_trials) if clean_trials["wholetrial"][i] == 0]

# compute quantities on clean trials
S = getfield.(trials, :s)
X = getfield.(trials, :x)
Y = getfield.(trials, :y)
Θ = getfield.(trials, :θ)
U = getfield.(trials, :u)
Ω = getfield.(trials, :ω)

# ---------------------------------------------------------------------------- #
#                               COMPARISON POINTS                              #
# ---------------------------------------------------------------------------- #
function get_comparison_points(trials; δs=2.0)
    # load comparison points
    cpoints = ComparisonPoints(FULLTRACK; δs=δs, trials=trials, s₀=1.0, s₁=250.0);
    cpoints_valid = cpoints.points[2:end-2]
    return cpoints, cpoints_valid
end

# ---------------------------------------------------------------------------- #
#                                 MTM SOLUTIONS                                #
# ---------------------------------------------------------------------------- #
function fix_solution_dtype(sol)
    sol.Fu = map(v-> v isa Number ? v : parse(Float64, v), sol.Fu)
    sol
end

fixvec(vec) = map(v-> v isa Number ? v : parse(Float64, v), vec)

function df2sol(df)
    _keys = filter!(k -> Symbol(k) != :δt, names(df))
    dt = "δt" ∈ names(df) ? df[1, :δt] : 0.0

    return Solution(; Dict(map(k -> Symbol(k)=>fixvec(df[:, k]), _keys))..., δt=dt)
end

function load_global_solution()

    if Sys.iswindows()
        fld = "D:\\Dropbox (UCL)\\Rotation_vte\\Locomotion\\analysis\\behavior"
    else
        fld = "/Users/federicoclaudi/Dropbox (UCL)/Rotation_vte/Locomotion/analysis/behavior"
    end

    gsol_path = joinpath(fld, "globalsolution.csv")
    globalsolution = DataFrame(CSV.File(gsol_path))
    return Solution((df2sol ∘ fix_solution_dtype)(globalsolution));
end

function load_global_solution(gsol_path)
    globalsolution = DataFrame(CSV.File(gsol_path))
    return Solution((df2sol ∘ fix_solution_dtype)(globalsolution));
end

function load_mtm_solutions(; folder=nothing, name="multiple_horizons_mtm_horizon_length")
    # load individual solutions
    folder = isnothing(folder) ? PATHS["horizons_sims_cache"] : folder
    files = sort(glob("$(name)*.csv", folder), lt=natural)

    solutions_df = map(file->DataFrame(CSV.File(file)), files)
    solutions_df = map(sol -> fix_solution_dtype(sol), solutions_df)

    solutions = df2sol.(solutions_df)
    # _names = map(file -> split(file, "_")[end][1:end-4], files);
    _names = map(file -> basename(file), files)
    return solutions, _names
end

"""
Get a solution's variable value at equaly spaced s values
between s0 and sf.

Effectively it resamples the selected variable at equally spaced S values.
"""
function get_solution_variable_at_s_values(solution::Union{Trial, Solution}, s0::Number, sf::Number, δs::Number, var::Symbol;  stop=nothing)
    S = collect(s0:δs:sf)

    stop = isnothing(stop) ? length(solution.s) : stop
    _s = solution.s[1:stop]
    _s = Interpolations.deduplicate_knots!(_s; move_knots=true)
    _x = getfield(solution, var)[1:stop]

    # x = interpolate_wrt_to(_s, _x)
    itp = LinearInterpolation(_s, _x, extrapolation_bc=Flat())
    return itp.(S)
end
# ---------------------------------------------------------------------------- #
#                                   UTILITIES                                  #
# ---------------------------------------------------------------------------- #
calcspeed(solution, idx) = sqrt.(solution.u[idx] .^ 2 .+ solution.v[idx] .^ 2)

euclidean(x::Number, y::Number, x2, y2) = sqrt((x2-x)^2 + (y-y2)^2)
euclidean(x::Vector, y::Vector, x2, y2) = sqrt.((x .- x2) .^ 2 .+ (y .- y2) .^ 2)

filterdf(df, col, val) = df[df[:, col] .== val, :]
filterdf(df, col1, val1, col2, val2) = df[(df[:, col1] .== val1) .* (df[:, col2] .== val2), :]

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
struct DictTable{I, T <: NamedTuple, Data <: NamedTuple{<:Any, <:Tuple{Vararg{AbstractDictionary}}}, Inds <: AbstractIndices{I}} <: AbstractDictionary{I, T}
    data::Data
    indices::Inds

    function DictTable{I, T, Data, Inds}(data::NamedTuple{<:Any, <:Tuple{Vararg{AbstractDictionary}}}, indices::AbstractIndices) where {I, T, Data, Inds}
        # All the data columns have to share tokens - copy if we have to
        data = map(data) do col
            if sharetokens(keys(col), indices)
                return col
            else
                out = similar(indices, eltype(col))
                copyto!(out, col)
                return out
            end
        end
        return new(data, indices)
    end
end

@generated function DictTable{I, T, NamedTuple{Names, Types}, Inds}() where {I, T, Inds, Names, Types}
    return quote
        inds = Inds()
        data = NamedTuple{Names, Types}(tuple($((:($(p)()) for p in Types.parameters)...)))
        DictTable{I, T, NamedTuple{Names, Types}, Inds}(data, inds)
    end
end

function DictTable(ts...; kwargs...)
    cols = removenothings_dict(merge(_columns(ts...), kwargs.data))
    inds = keys(first(cols))
    return DictTable{eltype(inds), _eltypes(cols), typeof(cols), typeof(inds)}(cols, inds)
end

_columns(t::DictTable) = columns(t)

@generated function removenothings_dict(nt::NamedTuple{names, T}) where {names, T}
    exprs = []
    newnames = []
    params = T.parameters
    first = true
    preamble = :()
    for i in 1:length(params)
        if params[i] === Nothing
            continue
        end

        push!(newnames, names[i])
        if first
            # If the first column is not a dictionary, then assume it is the primary key
            first = false
            if params[i] <: AbstractDictionary
                preamble = :(indices = keys(getfield(nt, $(QuoteNode(names[i])))))
                push!(exprs, :(getfield(nt, $(QuoteNode(names[i])))))
            else
                preamble = :(indices = Indices(getfield(nt, $(QuoteNode(names[i])))))
                push!(exprs, :indices)
            end
        else
            if params[i] <: AbstractDictionary
                push!(exprs, :(getfield(nt, $(QuoteNode(names[i])))))
            else
                push!(exprs, :(Dictionary(indices, getfield(nt, $(QuoteNode(names[i]))))))
            end
        end
    end

    return quote
        $preamble
        return NamedTuple{$(tuple(newnames...))}(tuple($(exprs...)))
    end
end

function DictTable(cols::NamedTuple{names, <:Tuple{Vararg{AbstractDictionary}}}) where names
    if isempty(names)
        throw(MethodError(DictTable, Tuple{NamedTuple{(), Tuple{}}}))
    end
    indices = keys(cols[names[1]])
    I = eltype(indices)
    T = _eltypes(cols)
    return DictTable{I, T, typeof(cols), typeof(indices)}(cols, indices)
end

function DictTable{I, T}() where {I, T <: NamedTuple}
    indices = Indices{I}()
    cols = _similar(indices, T)
    return DictTable{I, T, typeof(cols), typeof(indices)}(cols, indices)
end

@generated function _similar(indices, ::Type{T}) where {T <: NamedTuple}
    exprs = (:(similar(indices, $V)) for V in T.parameters[2].parameters)
    return :(NamedTuple{$(T.parameters[1])}(tuple($(exprs...))))
end

Base.similar(t::DictTable) = similar(t, eltype(t))
function Base.similar(t::DictTable, ::Type{T}) where {T <: NamedTuple}
    return DictTable(_similar(keys(t), T))
end

Dictionaries.empty_type(::Type{DictTable{<:Any, <:Any, <:Any, Inds}}, ::Type{I}) where {I, Inds} = Dictionaries.empty_type(Inds, I)

@generated function Dictionaries.empty_type(::Type{<:DictTable{<:Any, <:Any, <:Any, Inds}}, ::Type{I}, ::Type{NamedTuple{Names, T}}) where {I, Names, T, Inds}
    return quote
        DictTable{I, NamedTuple{Names, T}, NamedTuple{Names, Tuple{$((:(Dictionaries.empty_type(Inds, I, $p)) for p in T.parameters)...)}}, Dictionaries.empty_type(Inds, I)}
    end
end

function Base.empty(t::DictTable{I, <:NamedTuple{names, <:Any}}) where {I, names}
    # Check if first element contains the keys (TODO expand search to any element?)
    if length(names) > 0
        first_name = names[1]
        cols = columns(t)
        first_col = cols[first_name]
        if first_col isa AbstractIndices{I}
            first_col = empty(first_col)
            cols = merge(NamedTuple{(first_name,)}((first_col,)), _similar(first_col, typeof(Base.tail(cols))))
            return DictTable(cols)
        end
    end
    return empty(t, I, eltype(t))
end

function Base.empty(t::DictTable, ::Type{I}, ::Type{T}) where {I, T <: NamedTuple}
    indices = empty(keys(t), I)
    return DictTable(_similar(indices, T))
end

Base.keys(d::DictTable) = getfield(d, :indices)
Tables.columns(d::DictTable) = getfield(d, :data)
@inline Tables.rows(t::DictTable) = t

Base.getproperty(d::DictTable, s::Symbol) = getproperty(columns(d), s)

# Private fields are never exposed since they can conflict with column names
Base.propertynames(t::DictTable, private::Bool=false) = columnnames(t)

Tables.istable(::Type{<:DictTable}) = true
Tables.rowaccess(::Type{<:DictTable}) = true
Tables.columnaccess(::Type{<:DictTable}) = true
Tables.schema(::DictTable{<:Any, T}) where {T} = Tables.Schema(T)
Tables.materializer(::DictTable) = DictTable

function Dictionaries.istokenassigned(d::DictTable, t)
    return all(col -> istokenassigned(col, t), columns(d))
end

function Dictionaries.gettokenvalue(d::DictTable, t)
    return map(col -> gettokenvalue(col, t), columns(d))
end

function Dictionaries.issettable(d::DictTable)
    return all(col -> issettable(col), columns(d))
end

function Dictionaries.settokenvalue!(d::DictTable{I, V}, t, val::V) where {I, V}
    foreach((col, v) -> settokenvalue!(col, t, v), columns(d), val)
    return d
end

function Dictionaries.isinsertable(d::DictTable)
    return all(col -> isinsertable(col), columns(d))
end

function Dictionaries.gettoken!(d::DictTable{I}, i::I) where {I}
    return gettoken!(keys(d), i, map(tokenized, columns(d)))
end

function Dictionaries.deletetoken!(d::DictTable, t)
    return deletetoken!(keys(d), t, map(tokenized, columns(d)))
end

function Base.empty!(d::DictTable)
    empty!(keys(d), map(tokenized, columns(d)))
    return d
end

columnnames(::AbstractDictionary{<:Any, <:NamedTuple{names}}) where {names} = names

@generated function keyname(::DictTable{<:Any, <:Any, Data, Inds}) where { Data, Inds }
    params = Data.parameters[2].parameters
    for i in 1:length(params)
        if params[i] === Inds
            return QuoteNode(Data.parameters[1][i])
        end
    end
    return nothing
end

# show
function Base.show(io::IO, ::MIME"text/plain", t::DictTable)
    if get(io, :compact, false)
        invoke(show, Tuple{typeof(io), MIME"text/plain", AbstractDictionary}, io, MIME"text/plain"(), t)
    else
        showtable(io, t, keyname(t))
    end
end

# function Base.show(io::IO, t::DictTable)
#     if get(io, :compact, false)
#         invoke(show, Tuple{typeof(io), AbstractDictionary}, io, t)
#     else
#         showtable(io, t, keyname(t))
#     end
# end

propertytype(::DictTable) = DictTable
(::GetProperties{()})(t::DictTable) = FillDictionary(keys(t), (;))

# filter!

function Base.filter!(pred, t::DictTable)
    pred2 = token -> pred(map(col -> gettokenvalue(col, token), columns(t)))
    Dictionaries._filter!(pred2, keys(t), map(tokenized, columns(t)))
    return t
end

# Column-based operations: Some operations on rows are faster when considering columns

Base.map(::typeof(identity), t::DictTable) = copy(t)

Base.map(::typeof(merge), t::DictTable) = copy(t)

function Base.map(::typeof(merge), t1::DictTable, t2::DictTable)
    return copy(DictTable(merge(columns(t1), columns(t2))))
end

function Base.map(f::GetProperty, t::DictTable)
    return copy(f(t))
end

@inline function Base.map(f::GetProperties, t::DictTable)
    return copy(f(t))
end

@inline function Base.map(f::Compute{names}, t::DictTable) where {names}
    # minimize number of columns before iterating over the rows
    map(f, GetProperties(names)(t))
end

@inline function Base.map(f::Compute{names}, t::DictTable{<:Any, <:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return map(f.f, getproperty(names[1])(t))
    elseif length(names) == 2
        return map(f.f, getproperty(names[1])(t), getproperty(names[2])(t))
    end

    invoke(map, Tuple{Function, typeof(t)}, f, t)
end

@generated function Base.map(s::Select{names}, t::DictTable) where {names}
    exprs = [:($(names[i]) = map(s.fs[$i], t)) for i in 1:length(names)]

    return :(DictTable($(Expr(:tuple, exprs...))))
end

SplitApplyCombine.mapview(::typeof(merge), t::DictTable) = t

function SplitApplyCombine.mapview(::typeof(merge), t1::DictTable, t2::DictTable)
    return DictTable(merge(columns(t1), columns(t2)))
end

@inline function SplitApplyCombine.mapview(f::GetProperty, t::DictTable)
    return f(t)
end

@inline function SplitApplyCombine.mapview(f::GetProperties, t::DictTable)
    return f(t)
end

@inline function SplitApplyCombine.mapview(f::Compute{names}, t::DictTable) where {names}
    # minimize number of columns before iterating over the rows
    mapview(f, GetProperties(names)(t))
end

@inline function SplitApplyCombine.mapview(f::Compute{names}, t::DictTable{<:Any, <:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases (consider 2-argument version)
        return mapview(f.f, getproperty(names[1])(t))
    end

    invoke(mapview, Tuple{Function, typeof(t)}, f, t)
end

@generated function SplitApplyCombine.mapview(s::Select{names}, t::DictTable) where {names}
    exprs = [:($(names[i]) = mapview(s.fs[$i], t)) for i in 1:length(names)]

    return :(DictTable($(Expr(:tuple, exprs...))))
end

function Indexing.getindices(t::DictTable, inds::AbstractIndices)
    return DictTable(map(col -> getindices(col, inds), columns(t)))
end

function Base.view(t::DictTable, inds::AbstractIndices)
    return DictTable(map(col -> view(col, inds), columns(t)))
end

@inline function Broadcast.broadcasted(::Dictionaries.DictionaryStyle, ::typeof(merge), ts::DictTable...) where {N}
	return DictTable(merge(map(columns, ts)...))
end

@inline function Broadcast.broadcasted(::Dictionaries.DictionaryStyle, f::GetProperty, t::DictTable) where {N}
    return f(t)
end

@inline function Broadcast.broadcasted(::Dictionaries.DictionaryStyle, f::GetProperties, t::DictTable) where {N}
    return DictTable(f(t))
end

@inline function Broadcast.broadcasted(::Dictionaries.DictionaryStyle, f::Compute, t::DictTable) where {N}
    return map(f, t)
end

@inline function Broadcast.broadcasted(::Dictionaries.DictionaryStyle, f::Select, t::DictTable) where {N}
    return map(f, t)
end

# findall

function Base.findall(f::GetProperty, t::DictTable)
    return findall(identity, f(t))
end

function Base.findall(f::Compute{names}, t::DictTable) where {names}
    # minimize number of columns before iterating over the rows
    return findall(f, GetProperties(names)(t))
end

function Base.findall(f::Compute{names}, t::DictTable{<:NamedTuple{names}}) where {names}
    # efficient to iterate over rows with a minimal number of columns
    if length(names) == 1 # unwrap in the simple cases
        return findall(f.f, getproperty(names[1])(t))
    end

    invoke(findall, Tuple{Function, typeof(t)}, f, t)
end

#
# grouping
#

function SplitApplyCombine.group(groups::AbstractDictionary, values::DictTable)
    I = eltype(groups)
    T = eltype(values)
    D = Dictionaries.empty_type(typeof(values), keytype(values), T)

    indices = keys(values)
    out = empty(groups, I, D)

    if sharetokens(groups, values)
        @inbounds for token in tokens(groups)
            group = gettokenvalue(groups, token)
            key = gettokenvalue(indices, token)
            value = gettokenvalue(values, token)
            insert!(get!(D, out, group), key, value)
        end
    else
        if length(groups) != length(values)
            throw(KeyError("Indices must match"))
        end

        for (i, group) in pairs(groups)
            value = values[i]
            insert!(get!(D, out, group), i, value)
        end
    end

    return out
end


text:
export Emoji

"""
An emoji.
More details [here](https://discordapp.com/developers/docs/resources/emoji#emoji-object).
"""
struct Emoji
    id::Nullable{Snowflake}
    name::String
    roles::Optional{Vector{Snowflake}}
    user::Optional{User}
    require_colons::Optional{Bool}
    managed::Optional{Bool}
    animated::Optional{Bool}
end
@boilerplate Emoji :constructors :docs :lower :merge :mock


text:
function _canonicalize_period(p::Dates.CompoundPeriod)
    return Dates.canonicalize(p)
end

function _canonicalize_period(p::Dates.Period)
    return _canonicalize_period(Dates.CompoundPeriod(p))
end


text:
using EffectiveWaves
using LinearAlgebra, Statistics, Test

@testset "Pair-correlation tests" begin

    # import EffectiveWaves: gls_pair_radial_fun

    particle_medium = Acoustic(3; ρ=0.0, c=0.0);

    exclusion_distance = 1.01

    s1 = Specie(
        particle_medium, Sphere(0.8), exclusion_distance = exclusion_distance
    );

    # minimal distance between particle centres
    a12 = 2.0 * s1.exclusion_distance * outer_radius(s1)
    R = 10.0
    polynomial_order = 40

    pair_corr_inf(z) = hole_correction_pair_correlation([0.0,0.0,0.0],s1, [0.0,0.0,z],s1)

    # import EffectiveWaves: smooth_pair_corr_distance

    pair_corr_inf_smooth = smooth_pair_corr_distance(pair_corr_inf, a12; smoothing = 1.0, max_distance = 2R, polynomial_order = polynomial_order)

    # using Plots
    # rs = 0.0:0.1:(2R)
    # plot(rs,pair_corr_inf_smooth.(rs))
    # plot!(rs,pair_corr_inf.(rs))

    pair_radial = pair_radial_fun(pair_corr_inf_smooth, a12;
        sigma_approximation = false,
        polynomial_order = polynomial_order
    )

    r1 = 4.0; r2 = 4.2;
    cosθs = LinRange(0.6,1.0,200);

    pair_corr = pair_radial_to_pair_corr(pair_radial)

    p_corrs = map(eachindex(cosθs)) do i
        sinθ = sqrt(1.0 - cosθs[i]^2)
        x1 = [0.0, 0.0, r1]
        x2 = [r2 * sinθ, 0.0, r2 * cosθs[i]]
        pair_corr(x1,s1,x2,s1)
    end;

    p_rads = pair_radial.(r1,r2,cosθs);
    p_infs = pair_corr_inf_smooth.(sqrt.(r1^2 .+ r2.^2 .- 2r1 .* r2 .* cosθs))

    @test maximum(abs.(p_corrs - p_infs)) < 1e-3
    @test maximum(abs.(p_rads - p_infs)) < 1e-3

    polynomial_order = 90;
    pair_corr_inf(z) = hole_correction_pair_correlation([0.0,0.0,0.0],s1, [0.0,0.0,z],s1) #* (1 + sin(5*abs(z-a12)) * exp(-z))

    pair_radial = pair_radial_fun(pair_corr_inf, a12;
        sigma_approximation = true,
        # sigma = true,
        polynomial_order = polynomial_order,
        mesh_size = 12polynomial_order + 1
    )

    # polynomial_order = 10
    # sigmas = [one(T); sin.(pi .* ls[2:end] ./ (polynomial_order+1)) ./ (pi .* ls[2:end] ./ (polynomial_order+1))]
    # plot(sigmas)

    p_rads = pair_radial.(r1,r2,cosθs);
    p_infs = pair_corr_inf.(sqrt.(r1^2 .+ r2.^2 .- 2r1 .* r2 .* cosθs))
    # using Plots
    #
    # plot(cosθs,p_rads)
    # plot!(cosθs,p_infs)

    r1s = (20 * a12) .* rand(1000);
    r2s = (20 * a12) .* rand(1000);
    cosθs = LinRange(-1.0,1.0,1000);

    p_rads = pair_radial.(r1s,r2s,cosθs);
    p_infs = pair_corr_inf.(sqrt.(r1s.^2 .+ r2s.^2 .- 2r1s .* r2s .* cosθs));

    @test mean(abs.(p_rads - p_infs)) < 1e-3
    @test maximum(abs.(p_rads - p_infs)) < 0.6 # infinite sharp function may have point wise error of 0.5, but exceptionally rare.
end

max_stars_repo_name: arturgower/EffectiveWaves

max_issues_repo_name: arturgower/EffectiveWaves

max_forks_repo_name: arturgower/EffectiveWaves


text:
using SFML

function handle_client(client)
    while true
        packet = Packet()

        while receive(client, packet) != SOCKET_DONE
            sleep(0)
        end

        msg = read_string(packet)

        println("Received: $msg")

        out_packet = Packet()
        write(out_packet, msg)
        send(client, out_packet)
    end
end

server = TcpListener()
set_blocking(server, false)

if listen(server, 8080) != SOCKET_DONE
    println("Failed")
    exit(1)
end

println("Listening on port 8080")

while true
    client = SocketTCP()
    while accept(server, client) != SOCKET_DONE
        sleep(0)
    end
    println("Client connected from $(get_remoteaddress(client))")

    @async begin
        handle_client(client)
    end
end