text:
struct SymTen{T<:Number} <: AbstractTen{T}
    xx::T
    xy::T
    xz::T
    yy::T
    yz::T
    zz::T
end

Base.IndexStyle(a::Type{<:SymTen}) = Base.IndexLinear()

@inline function Base.getindex(a::SymTen,I::Integer)
    I <= 3 && return getfield(a,I)
    I == 4 && return a.xy
    I == 5 && return a.yy
    I == 6 && return a.yz
    I == 7 && return a.xz
    I == 8 && return a.yz
    I == 9 && return a.zz
end

@inline Base.size(a::SymTen) =
    (3,3)

@inline Base.length(a::SymTen) =
    9

@inline Base.zero(a::Type{SymTen{T}}) where {T} =
    SymTen{T}(zero(T),zero(T),zero(T),zero(T),zero(T),zero(T))

@inline SymTen(xx::T,xy::T,xz::T,yy::T,yz::T,zz::T) where {T} =
    SymTen{T}(xx,xy,xz,yy,yz,zz)

@inline SymTen(q,w,e,r,t,y) =
    SymTen(promote(q,w,e,r,t,y)...)

@inline SymTen(s::SymTen) = SymTen(s.xx,s.xy,s.xz,s.yy,s.yz,s.zz)

@inline Base.:+(a::SymTen{T},b::SymTen{T2}) where {T,T2} =
    SymTen{promote_type(T,T2)}(a.xx+b.xx, a.xy+b.xy, a.xz+b.xz, a.yy+b.yy, a.yz + b.yz, a.zz+b.zz)

@inline Base.:-(A::SymTen{T}) where T = SymTen{T}(-A.xx,-A.xy,-A.xz,-A.yy,-A.yz,-A.zz)

@inline Base.:-(a::SymTen{T},b::SymTen{T2}) where {T,T2} =
    SymTen{promote_type(T,T2)}(a.xx-b.xx, a.xy-b.xy, a.xz-b.xz, a.yy-b.yy, a.yz - b.yz, a.zz - b.zz)

@inline Base.:*(a::T,b::SymTen{T2}) where {T<:Number,T2<:Number} =
    SymTen{promote_type(T,T2)}(a*b.xx, a*b.xy, a*b.xz, a*b.yy, a*b.yz, a*b.zz)

@inline Base.:*(b::SymTen{T},a::T2) where {T<:Number,T2<:Number} =
    a*b

@inline Base.:/(b::SymTen{T},a::T2) where {T<:Number,T2<:Number} =
    inv(a)*b

@inline Base.:(:)(a::SymTen,b::SymTen) =
    muladd(a.xx, b.xx, muladd(a.yy, b.yy, a.zz*b.zz)) + 2muladd(a.xy, b.xy, muladd(a.xz, b.xz, a.yz*b.yz))

@inline Base.:(+)(L::LinearAlgebra.UniformScaling,A::SymTen) =
    SymTen(A.xx + L.λ, A.xy,A.xz,A.yy + L.λ, A.yz, A.zz + L.λ)

@inline Base.:(+)(A::SymTen, L::LinearAlgebra.UniformScaling) = L+A

@inline Base.:(-)(A::SymTen, L::LinearAlgebra.UniformScaling) =
    SymTen(A.xx - L.λ, A.xy,A.xz,A.yy - L.λ, A.yz, A.zz - L.λ)

@inline Base.:(-)(L::LinearAlgebra.UniformScaling,A::SymTen) = -(A-L)

@inline LinearAlgebra.tr(a::SymTen) =
    a.xx + a.yy + a.zz

@inline LinearAlgebra.norm(a::SymTen) =
    @fastmath sqrt(2(a:a))

@inline LinearAlgebra.det(t::SymTen) =
    #-txx*tyz^2 - txz^2*tyy - tzz*txy^2 + 2*txz*txy*tyz + tzz*txx*tyy
    t.xx*t.yy*t.zz + 2*t.xz*t.xy*t.yz - (t.xx*t.yz*t.yz + t.yy*t.xz*t.xz + t.zz*t.xy*t.xy)

@inline LinearAlgebra.dot(a::SymTen{T},b::Vec{T2}) where {T<:Number,T2<:Number} =
    Vec{promote_type(T,T2)}(muladd(a.xx, b.x, muladd(a.xy, b.y, a.xz*b.z)),
        muladd(a.xy, b.x, muladd(a.yy, b.y, a.yz*b.z)),
        muladd(a.xz, b.x, muladd(a.yz, b.y, a.zz*b.z)))

@inline LinearAlgebra.dot(b::Vec,a::SymTen) =
    LinearAlgebra.dot(a,b)

@inline LinearAlgebra.dot(A::AbstractMatrix,v::Vec) = A*v

@inline symouter(a::Vec,b::Vec) = #symmetric part of the outer product of two vectors
    SymTen(a.x*b.x,
           muladd(a.x,b.y,a.y*b.x)/2,
           muladd(a.x,b.z,a.z*b.x)/2,
           a.y*b.y,
           muladd(a.y,b.z,a.z*b.y)/2,
           a.z*b.z)

@inline antisymouter(a::Vec,b::Vec) =
    AntiSymTen((a.x*b.y - a.y*b.x)/2,
               (a.x*b.z - a.z*b.x)/2,
               (a.y*b.z - a.z*b.y)/2)

@inline traceless(S::SymTen) =
    S - inv(3)*tr(S)*I

@inline function Base.:^(t::SymTen{T},n::Integer) where {T}
    λ,e = eigvec(t)
    r = (λ[1]^n)*(symouter(e[1],e[1])) + (λ[2]^n)*(symouter(e[2],e[2])) + (λ[3]^n)*(symouter(e[3],e[3]))
end

@inline square(t::SymTen) =
    SymTen(muladd(t.xx,t.xx,muladd(t.xy,t.xy,t.xz^2)),
              muladd(t.xx,t.xy,muladd(t.yy,t.xy,t.xz*t.yz)),
              muladd(t.xx,t.xz,muladd(t.xy,t.yz,t.zz*t.xz)),
              muladd(t.xy,t.xy,muladd(t.yy,t.yy,t.yz^2)),
              muladd(t.xy,t.xz,muladd(t.yy,t.yz,t.zz*t.yz)),
              muladd(t.xz,t.xz,muladd(t.yz,t.yz,t.zz^2)))

@inline function Base.literal_pow(::typeof(Base.:^),t::SymTen{T},::Val{2}) where {T<:AbstractFloat}
    return square(t)
end

@inline Base.adjoint(a::SymTen{T}) where T<:Complex = SymTen{T}(conj(a.xx),conj(a.xy),conj(a.xz),conj(a.yy),conj(a.yz),conj(a.zz))
@inline Base.adjoint(a::SymTen) = a

####### Full Tensors

struct Ten{T<:Number} <: AbstractTen{T}
    xx::T
    yx::T
    zx::T
    xy::T
    yy::T
    zy::T
    xz::T
    yz::T
    zz::T
end

Base.IndexStyle(a::Type{<:Ten}) = Base.IndexLinear()

@inline Base.@propagate_inbounds function Base.getindex(a::Ten,I::Integer)
    @boundscheck checkbounds(a,I)
    @inbounds getfield(a,I)
end

@inline Base.size(a::Ten) =
    (3,3)

@inline Base.length(a::Ten) =
    9

@inline Base.zero(a::Type{Ten{T}}) where {T} =
    Ten{T}(zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T),zero(T))

@inline Ten(xx::T,xy::T,xz::T,yy::T,yz::T,zz::T) where {T} =
    Ten{T}(xx,xy,xz,yy,yz,zz)

@inline Ten(q,w,e,r,t,y,a,d,f) =
    Ten(promote(q,w,e,r,t,y,a,d,f)...)

@inline Base.:+(a::Ten{T},b::Ten{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx+b.xx, a.yx+b.yx, a.zx+b.zx, a.xy+b.xy, a.yy+b.yy, a.zy+b.zy, a.xz+b.xz, a.yz + b.yz, a.zz+b.zz)

@inline Base.:-(A::Ten{T}) where T = Ten{T}(-A.xx,-A.yx,-A.zx,-A.xy,-A.yy,-A.zy,-A.xz,-A.yz,-A.zz)

@inline Base.:-(a::Ten{T},b::Ten{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx-b.xx, a.yx-b.yx, a.zx-b.zx, a.xy-b.xy, a.yy - b.yy, a.zy - b.zy, a.xz - b.xz, a.yz - b.yz, a.zz - b.zz)

@inline Base.:*(a::T,b::Ten{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(a*b.xx, a*b.yx, a*b.zx, a*b.xy, a*b.yy, a*b.zy, a*b.xz, a*b.yz, a*b.zz)

@inline Base.:*(b::Ten{T},a::T2) where {T<:Number,T2<:Number} =
    a*b

@inline Base.:/(b::Ten{T},a::T2) where {T<:Number,T2<:Number} =
    inv(a)*b

@inline Base.:(:)(a::Ten,b::Ten) =
    muladd(a.xx, b.xx, muladd(a.yx, b.yx, muladd(a.zx, b.zx, muladd(a.xy, b.xy, muladd(a.yy, b.yy, muladd(a.zy, b.zy, muladd(a.xz, b.xz, muladd(a.yz, b.yz, a.zz*b.zz))))))))

@inline Base.:(:)(a::SymTen,b::Ten) =
    muladd(a.xx, b.xx, muladd(a.xy, b.yx, muladd(a.xz, b.zx, muladd(a.xy, b.xy, muladd(a.yy, b.yy, muladd(a.yz, b.zy, muladd(a.xz, b.xz, muladd(a.yz, b.yz, a.zz*b.zz))))))))

@inline Base.:(:)(b::Ten,a::SymTen) = a:b

@inline Base.:(+)(L::LinearAlgebra.UniformScaling,A::Ten) =
    Ten(A.xx + L.λ, A.yx, A.zx, A.xy, A.yy + L.λ, A.zy, A.xz, A.yz, A.zz + L.λ)

@inline Base.:(+)(A::Ten, L::LinearAlgebra.UniformScaling) = L+A

@inline Base.:(-)(A::Ten, L::LinearAlgebra.UniformScaling) =
    Ten(A.xx - L.λ, A.yx, A.zx, A.xy, A.yy - L.λ, A.zy, A.xz, A.yz, A.zz - L.λ)

@inline Base.:(-)(L::LinearAlgebra.UniformScaling,A::Ten) = -(A-L)

@inline LinearAlgebra.tr(a::Ten) =
    a.xx + a.yy + a.zz

@inline LinearAlgebra.norm(a::Ten) =
    @fastmath sqrt(2(a:a))

@inline LinearAlgebra.det(t::Ten) =
    (t.xy*t.yz - t.xz*t.yy)*t.zx + (t.zy*t.xz - t.zz*t.xy)*t.yx + (t.zz*t.yy - t.zy*t.yz )*t.xx

@inline LinearAlgebra.dot(a::Ten{T},b::Ten{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(muladd(a.xx, b.xx, muladd(a.xy, b.yx, a.xz*b.zx)), muladd(a.yx, b.xx, muladd(a.yy, b.yx, a.yz*b.zx)), muladd(a.zx, b.xx, muladd(a.zy, b.yx, a.zz*b.zx)),
                                      muladd(a.xx, b.xy, muladd(a.xy, b.yy, a.xz*b.zy)), muladd(a.yx, b.xy, muladd(a.yy, b.yy, a.yz*b.zy)), muladd(a.zx, b.xy, muladd(a.zy, b.yy, a.zz*b.zy)),
                                      muladd(a.xx, b.xz, muladd(a.xy, b.yz, a.xz*b.zz)), muladd(a.yx, b.xz, muladd(a.yy, b.yz, a.yz*b.zz)), muladd(a.zx, b.xz, muladd(a.zy, b.yz, a.zz*b.zz)))

@inline LinearAlgebra.dot(a::Ten{T},b::Vec{T2}) where {T<:Number,T2<:Number} =
    Vec{promote_type(T,T2)}(muladd(a.xx, b.x, muladd(a.xy, b.y, a.xz*b.z)),
                                      muladd(a.yx, b.x, muladd(a.yy, b.y, a.yz*b.z)),
                                      muladd(a.zx, b.x, muladd(a.zy, b.y, a.zz*b.z)))

@inline LinearAlgebra.dot(b::Vec{T},a::Ten{T2}) where {T,T2} =
    Vec{promote_type(T,T2)}(muladd(a.xx, b.x, muladd(a.yx, b.y, a.zx*b.z)),
                                      muladd(a.xy, b.x, muladd(a.yy, b.y, a.zy*b.z)),
                                      muladd(a.xz, b.x, muladd(a.yz, b.y, a.zz*b.z)))

@inline LinearAlgebra.dot(a::SymTen{T},b::SymTen{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(muladd(a.xx, b.xx, muladd(a.xy, b.xy, a.xz*b.xz)), muladd(a.xy, b.xx, muladd(a.yy, b.xy, a.yz*b.xz)), muladd(a.xz, b.xx, muladd(a.yz, b.xy, a.zz*b.xz)),
    muladd(a.xx, b.xy, muladd(a.xy, b.yy, a.xz*b.yz)), muladd(a.xy, b.xy, muladd(a.yy, b.yy, a.yz*b.yz)), muladd(a.xz, b.xy, muladd(a.yz, b.yy, a.zz*b.yz)),
    muladd(a.xx, b.xz, muladd(a.xy, b.yz, a.xz*b.zz)), muladd(a.xy, b.xz, muladd(a.yy, b.yz, a.yz*b.zz)), muladd(a.xz, b.xz, muladd(a.yz, b.yz, a.zz*b.zz)))

@inline Base.adjoint(a::Ten{T}) where T = Ten{T}(conj(a.xx),conj(a.xy),conj(a.xz),conj(a.yx),conj(a.yy),conj(a.yz),conj(a.zx),conj(a.zy),conj(a.zz))

@inline outer(a::Vec,b::Vec) = #symmetric part of the outer product of two vectors
    Ten(a.x*b.x, a.y*b.x, a.z*b.x,
                  a.x*b.y, a.y*b.y, a.z*b.y,
                  a.x*b.z, a.y*b.z, a.z*b.z)

const ⊗ = outer

@inline traceless(S::Ten) =
    S - inv(3)*LinearAlgebra.tr(S)*LinearAlgebra.I

function Base.:^(t::Ten{T},n::Integer) where {T}
    n == 1 && return t
    n == 2 && return t⋅t
    return t⋅(t^(n-1))
end

@inline square(t::Ten) = t⋅t

@inline function Base.literal_pow(::typeof(Base.:^),t::Ten{T},::Val{2}) where {T<:AbstractFloat}
    return t⋅t
end

@inline function SymTen(t::Ten)
    if (t.yx == t.xy && t.zx == t.xz && t.zy == t.yz)
       return SymTen(t.xx, t.yx, t.zx, t.yy, t.zy, t.zz)
    else
        throw(InexactError(:SymTen,SymTen,t))
    end
end
@inline Ten(a::SymTen) = Ten(a.xx,a.xy,a.xz,a.xy,a.yy,a.yz,a.xz,a.yz,a.zz)

####### AntiSym Tensors

struct AntiSymTen{T<:Number} <: AbstractTen{T}
    xy::T
    xz::T
    yz::T
end

@inline Ten(a::AntiSymTen) = Ten(0,-a.xy,-a.xz,a.xy,0,-a.yz,a.xz,a.yz,0)

Base.IndexStyle(a::Type{<:AntiSymTen}) = Base.IndexLinear()

@inline Base.@propagate_inbounds function Base.getindex(a::AntiSymTen{T},I::Integer) where T
    @boundscheck checkbounds(a,I)
    (I == 1 || I == 5 || I == 9) && return zero(T)
    I == 2 && return -a.xy
    I == 3 && return -a.xz
    I == 4 && return a.xy
    I == 6 && return -a.yz
    I == 7 && return a.xz
    I == 8 && return a.yz
end

@inline Base.size(a::AntiSymTen) =
    (3,3)

@inline Base.length(a::AntiSymTen) =
    9

@inline Base.zero(a::Type{AntiSymTen{T}}) where {T} =
    AntiSymTen{T}(zero(T),zero(T),zero(T))

@inline AntiSymTen(xx::T,xy::T,xz::T) where {T} =
    AntiSymTen{T}(xx,xy,xz)

@inline AntiSymTen(q,w,e) =
    AntiSymTen(promote(q,w,e)...)

@inline Base.:+(a::AntiSymTen{T},b::AntiSymTen{T2}) where {T,T2} =
    AntiSymTen{promote_type(T,T2)}(a.xy+b.xy, a.xz+b.xz, a.yz + b.yz)

@inline Base.:-(A::AntiSymTen{T}) where T = AntiSymTen{T}(-A.xy,-A.xz,-A.yz)

@inline Base.:-(a::AntiSymTen{T},b::AntiSymTen{T2}) where {T,T2} =
    AntiSymTen{promote_type(T,T2)}(a.xy-b.xy, a.xz - b.xz, a.yz - b.yz)

@inline Base.:*(a::T,b::AntiSymTen{T2}) where {T<:Number,T2<:Number} =
    AntiSymTen{promote_type(T,T2)}(a*b.xy, a*b.xz, a*b.yz)

@inline Base.:*(b::AntiSymTen{T},a::T2) where {T<:Number,T2<:Number} =
    a*b

@inline Base.:/(b::AntiSymTen{T},a::T2) where {T<:Number,T2<:Number} =
    inv(a)*b

@inline Base.:(:)(a::AntiSymTen,b::AntiSymTen) =
    2muladd(a.xy, b.xy, muladd(a.xz, b.xz, a.yz*b.yz))

@inline Base.:(:)(a::Ten,b::AntiSymTen) =
    muladd(a.xy,b.xy,muladd(a.xz,b.xz,a.yz*b.yz)) - muladd(a.yx,b.xy,muladd(a.zx,b.xz, a.zy*b.yz))

@inline Base.:(:)(b::AntiSymTen,a::Ten) = a:b

@inline Base.:(:)(a::AntiSymTen,b::SymTen) = false

@inline Base.:(:)(a::SymTen,b::AntiSymTen) = false

@inline Base.:(+)(L::LinearAlgebra.UniformScaling,A::AntiSymTen) =
    Ten(L.λ, -A.xy, -A.xz, A.xy, L.λ, -A.yz, A.xz, A.yz, L.λ)

@inline Base.:(+)(A::AntiSymTen, L::LinearAlgebra.UniformScaling) = L+A

@inline Base.:(-)(A::AntiSymTen, L::LinearAlgebra.UniformScaling) =
    Ten(-L.λ, -A.xy, -A.xz, A.xy, -L.λ, -A.yz, A.xz, A.yz, -L.λ)

@inline Base.:(-)(L::LinearAlgebra.UniformScaling,A::AntiSymTen) = -(A-L)

@inline LinearAlgebra.tr(a::AntiSymTen{T}) where T = zero(T)

@inline LinearAlgebra.norm(a::AntiSymTen) =
    @fastmath sqrt(2(a:a))

@inline LinearAlgebra.det(a::AntiSymTen) = false

@inline LinearAlgebra.dot(a::AntiSymTen{T},b::AntiSymTen{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(-muladd(a.xy, b.xy, a.xz*b.xz), -a.yz*b.xz, a.yz*b.xy,
                                      -a.xz*b.yz, -muladd(a.xy,b.xy,a.yz*b.yz), -a.xz*b.xy,
                                      a.xy*b.yz, -a.xy*b.xz, -muladd(a.xz,b.xz,a.yz*b.yz))

@inline LinearAlgebra.dot(a::AntiSymTen{T},b::Vec{T2}) where {T<:Number,T2<:Number} =
    Vec{promote_type(T,T2)}(muladd(a.xy, b.y, a.xz*b.z),
                                      muladd(a.yz, b.z, -a.xy*b.x),
                                      -muladd(a.xz, b.x, a.yz*b.y))

@inline LinearAlgebra.dot(b::Vec,a::AntiSymTen) = -a⋅b

@inline Base.adjoint(a::AntiSymTen{T}) where T = AntiSymTen{T}(conj(-a.xy),conj(-a.xz),conj(-a.yz))

@inline traceless(S::AntiSymTen) = S

Base.:^(t::AntiSymTen{T},n::Integer) where {T} = Ten(t)^n

@inline square(a::AntiSymTen) = SymTen(-muladd(a.xy,a.xy,a.xz^2), -a.xz*a.yz, a.xy*a.yz, -muladd(a.xy,a.xy,a.yz^2), -a.xy*a.xz, -muladd(a.xz,a.xz,a.yz^2))

@inline function Base.literal_pow(::typeof(Base.:^),a::AntiSymTen{T},::Val{2}) where {T<:AbstractFloat}
    return square(a)
end

@inline Ten(a::AntiSymTen{T}) where {T} = Ten{T}(zero(T),-a.xy,-a.xz,a.xy,zero(T),-a.yz,a.xz,a.yz,zero(T))
@inline Vec(a::AntiSymTen) = Vec(-a.yz,a.xz,-a.xy)
@inline AntiSymTen(a::Vec) = AntiSymTen(-a.z,a.y,-a.x)

@inline Lie(s::SymTen,w::AntiSymTen) =
    SymTen(-2*(w.xy*s.xy + w.xz*s.xz),
           w.xy*s.xx - w.yz*s.xz - w.xy*s.yy - w.xz*s.yz,
           w.xz*s.xx + w.yz*s.xy - w.xy*s.yz - w.xz*s.zz,
           2*(w.xy*s.xy - w.yz*s.yz),
           w.xy*s.xz + w.xz*s.xy + w.yz*s.yy - w.yz*s.zz,
           2*(w.xz*s.xz + w.yz*s.yz))

@inline Lie(w::AntiSymTen,s::SymTen) = Lie(s,-w)
@inline Lie(t::AbstractTen,l::AbstractTen) = t⋅l - l⋅t

############################## Interoperation

@inline Base.:+(a::AntiSymTen{T},b::SymTen{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(b.xx, b.xy-a.xy, b.xz-a.xz, a.xy+b.xy, b.yy, b.yz - a.yz, a.xz+b.xz, a.yz+b.yz, b.zz)

@inline Base.:+(a::AntiSymTen{T},b::Ten{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(b.xx, b.yx-a.xy, b.zx-a.xz, a.xy+b.xy, b.yy, b.zy - a.yz, a.xz+b.xz, a.yz+b.yz, b.zz)

@inline Base.:+(a::SymTen{T},b::Ten{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx+b.xx, a.xy+b.yx, a.xz+b.zx, a.xy+b.xy, a.yy+b.yy, a.yz+b.zy, a.xz+b.xz, a.yz + b.yz, a.zz+b.zz)

@inline Base.:+(b::SymTen,a::AntiSymTen) = a+b

@inline Base.:+(b::Ten,a::AntiSymTen) = a+b

@inline Base.:+(b::Ten,a::SymTen) = a+b

@inline Base.:-(a::AntiSymTen{T},b::SymTen{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(-b.xx, -a.xy-b.xy, -a.xz-b.xz, a.xy-b.xy, -b.yy, -a.yz - b.yz, a.xz - b.xz, a.yz - b.yz, -b.zz)

@inline Base.:-(a::SymTen{T},b::AntiSymTen{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx, a.xy+b.xy, a.xz+b.xz, a.xy-b.xy, a.yy, a.yz + b.yz, a.xz - b.xz, a.yz - b.yz, a.zz)

@inline Base.:-(a::Ten{T},b::AntiSymTen{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx, a.yx+b.xy, a.zx+b.xz, a.xy-b.xy, a.yy, a.zy + b.yz, a.xz - b.xz, a.yz - b.yz, a.zz)

@inline Base.:-(a::AntiSymTen{T},b::Ten{T2}) where {T,T2} = -(b-a)

@inline Base.:-(a::Ten{T},b::SymTen{T2}) where {T,T2} =
    Ten{promote_type(T,T2)}(a.xx-b.xx, a.yx-b.xy, a.zx-b.xz, a.xy-b.xy, a.yy - b.yy, a.zy - b.yz, a.xz - b.xz, a.yz - b.yz, a.zz - b.zz)

@inline Base.:-(a::SymTen{T},b::Ten{T2}) where {T,T2} = -(b-a)

@inline LinearAlgebra.dot(a::Ten{T},b::SymTen{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(muladd(a.xx, b.xx, muladd(a.xy, b.xy, a.xz*b.xz)), muladd(a.yx, b.xx, muladd(a.yy, b.xy, a.yz*b.xz)), muladd(a.zx, b.xx, muladd(a.zy, b.xy, a.zz*b.xz)),
                                      muladd(a.xx, b.xy, muladd(a.xy, b.yy, a.xz*b.yz)), muladd(a.yx, b.xy, muladd(a.yy, b.yy, a.yz*b.yz)), muladd(a.zx, b.xy, muladd(a.zy, b.yy, a.zz*b.yz)),
                                      muladd(a.xx, b.xz, muladd(a.xy, b.yz, a.xz*b.zz)), muladd(a.yx, b.xz, muladd(a.yy, b.yz, a.yz*b.zz)), muladd(a.zx, b.xz, muladd(a.zy, b.yz, a.zz*b.zz)))

@inline LinearAlgebra.dot(a::SymTen{T},b::Ten{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(muladd(a.xx, b.xx, muladd(a.xy, b.yx, a.xz*b.zx)), muladd(a.xy, b.xx, muladd(a.yy, b.yx, a.yz*b.zx)), muladd(a.xz, b.xx, muladd(a.yz, b.yx, a.zz*b.zx)),
    muladd(a.xx, b.xy, muladd(a.xy, b.yy, a.xz*b.zy)), muladd(a.xy, b.xy, muladd(a.yy, b.yy, a.yz*b.zy)), muladd(a.xz, b.xy, muladd(a.yz, b.yy, a.zz*b.zy)),
    muladd(a.xx, b.xz, muladd(a.xy, b.yz, a.xz*b.zz)), muladd(a.xy, b.xz, muladd(a.yy, b.yz, a.yz*b.zz)), muladd(a.xz, b.xz, muladd(a.yz, b.yz, a.zz*b.zz)))

@inline LinearAlgebra.dot(a::Ten{T},b::AntiSymTen{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(-muladd(a.xy,b.xy,a.xz*b.xz), -muladd(a.yy,b.xy, a.yz*b.xz), -muladd(a.zy,b.xy, a.zz*b.xz),
                                      a.xx*b.xy - a.xz*b.yz, a.yx*b.xy - a.yz*b.yz, a.zx*b.xy - a.zz*b.yz,
                                      muladd(a.xx,b.xz,a.xy*b.yz), muladd(a.yx,b.xz,a.yy*b.yz), muladd(a.zx,b.xz, a.zy*b.yz))

@inline LinearAlgebra.dot(b::AntiSymTen{T2},a::Ten{T}) where {T,T2} =
    Ten{promote_type(T,T2)}(muladd(b.xy, a.yx, b.xz*a.zx), -a.xx*b.xy+a.zx*b.yz, -muladd(a.xx,b.xz,a.yx*b.yz),
                                      muladd(a.yy,b.xy,a.zy*b.xz), -a.xy*b.xy+a.zy*b.yz, -muladd(a.xy,b.xz,a.yy*b.yz),
                                      muladd(a.yz,b.xy,a.zz*b.xz), -a.xz*b.xy+a.zz*b.yz, -muladd(a.xz,b.xz,a.yz*b.yz))

@inline LinearAlgebra.dot(a::SymTen{T},b::AntiSymTen{T2}) where {T<:Number,T2<:Number} =
    Ten{promote_type(T,T2)}(-muladd(a.xy,b.xy,a.xz*b.xz), -muladd(a.yy,b.xy, a.yz*b.xz), -muladd(a.yz,b.xy, a.zz*b.xz),
                                      a.xx*b.xy - a.xz*b.yz, a.xy*b.xy - a.yz*b.yz, a.xz*b.xy - a.zz*b.yz,
                                      muladd(a.xx,b.xz,a.xy*b.yz), muladd(a.xy,b.xz,a.yy*b.yz), muladd(a.xz,b.xz, a.yz*b.yz))

@inline LinearAlgebra.dot(b::AntiSymTen{T2},a::SymTen{T}) where {T,T2} =
    Ten{promote_type(T,T2)}(muladd(a.xy,b.xy,a.xz*b.xz), -a.xx*b.xy+a.xz*b.yz, -muladd(a.xx,b.xz,a.xy*b.yz),
                                      muladd(a.yy,b.xy,a.yz*b.xz), -a.xy*b.xy+a.yz*b.yz, -muladd(a.xy,b.xz,a.yy*b.yz),
                                      muladd(a.yz,b.xy,a.zz*b.xz), -a.xz*b.xy+a.zz*b.yz, -muladd(a.xz,b.xz,a.yz*b.yz))

@inline symmetric(t::SymTen) = t
@inline symmetric(t::AntiSymTen{T}) where {T} = zero(SymTen{T})
@inline symmetric(t::Ten) = SymTen(t.xx, (t.yx+t.xy)/2, (t.xz+t.zx)/2, t.yy, (t.zy+t.yz)/2, t.zz)

@inline antisymmetric(t::SymTen{T}) where {T} = zero(AntiSymTen{T})
@inline antisymmetric(t::AntiSymTen) = t
@inline antisymmetric(t::Ten) = AntiSymTen((t.xy-t.yx)/2, (t.xz-t.zx)/2, (t.yz-t.zy)/2)

include("tensorarrays.jl")


text:
using SymbolicUtils: Sym, FnType, Term, symtype
using SymbolicUtils
using Test

@testset "@syms" begin
    let
        @syms a b::Float64 f(::Real) g(p, h(q::Real))::Int

        @test a isa Sym{Number}
        @test a.name === :a

        @test b isa Sym{Float64}
        @test b.name === :b

        @test f isa Sym{FnType{Tuple{Real}, Number}}
        @test f.name === :f

        @test g isa Sym{FnType{Tuple{Number, FnType{Tuple{Real}, Number}}, Int}}
        @test g.name === :g

        @test f(b) isa Term
        @test symtype(f(b)) === Number
        @test_throws ErrorException f(a)

        @test g(b, f) isa Term
        @test_throws ErrorException g(b, a)

        @test symtype(g(b, f)) === Int
    end
end

@testset "methods test" begin
    @syms w::Complex z::Complex a::Real b::Real x

    @test w + z == Term{Complex}(+, [w, z])
    @test z + a == Term{Number}(+, [z, a])
    @test a + b == Term{Real}(+, [a, b])
    @test a + x == Term{Number}(+, [a, x])
    @test a + z == Term{Number}(+, [a, z])

    # promote_symtype of identity
    @test Term(identity, [w]) == Term{Complex}(identity, [w])
    @test +(w) == w
    @test +(a) == a

    @test rem2pi(a, RoundNearest) == Term{Real}(rem2pi, [a, RoundNearest])
end

@testset "err test" begin
    @syms t()
    @test_throws ErrorException t(2)
end

ext: jl


text:
include("ssa.jl")
using .SSA

n = 101
D = 2
J = 10
M = 20

x = zeros(n, D, J)

y = hcat([sin.(a:0.3:(a+30)) for a=1:10]...)
x[:, 1, :] = y
x[:, 2, :] = -2*y

# M-SSA with multiple subseries
ssa_info = ssa_decompose(x, M)
r = SSA.ssa_reconstruct(ssa_info, 1:2, sum_modes=true)
@assert isapprox(x, r)
@assert size(x) == size(r)

# Regular M-SSA
x_mssa = x[:, :, 1]
ssa_info = ssa_decompose(x_mssa, M)
r = SSA.ssa_reconstruct(ssa_info, 1:2, sum_modes=true)
@assert isapprox(x_mssa, r)
@assert size(x_mssa) == size(r)

# Single-channel SSA
x_ssa = x[:, 1, 1]
ssa_info = ssa_decompose(x_ssa, M)
r = SSA.ssa_reconstruct(ssa_info, 1:2, sum_modes=true)
@assert isapprox(x_ssa, r)
@assert size(x_ssa) == size(r)

# Single-channel SSA with multiple subseries
x_ssa2 = x[:, 1:1, :]
ssa_info = ssa_decompose(x_ssa2, M)
r = SSA.ssa_reconstruct(ssa_info, 1:2, sum_modes=true)
@assert isapprox(x_ssa2, r)
@assert size(x_ssa2) == size(r)


text:
"""
    gradient(vi::VarInfo, model::Function, spl::Union{Nothing, Sampler})

Function to generate the gradient dictionary, with each prior map to its derivative of the logjoint probibilioty. This function uses chunk-wise forward AD with a chunk of size $(CHUNKSIZE) as default.

Example:

```julia
grad = gradient(vi, model, spl)
end
```
"""
gradient(vi::VarInfo, model::Function) = gradient(vi, model, nothing)
gradient(vi::VarInfo, model::Function, spl::Union{Nothing, Sampler}) = begin

  θ_hash = hash(vi[spl])

  if spl != nothing && haskey(spl.info, :grad_cache)
    if haskey(spl.info[:grad_cache], θ_hash)
      return spl.info[:grad_cache][θ_hash]
    end
  end

  # Initialisation
  grad = Vector{Float64}()

  # Split keys(vi) into chunks,
  @debug "making chunks..."
  vn_chunk = Set{VarName}(); vn_chunks = []; chunk_dim = 0;

  vns = getvns(vi, spl); vn_num = length(vns)

  for i = 1:vn_num
    l = length(getrange(vi, vns[i]))           # dimension for the current variable
    if chunk_dim + l > CHUNKSIZE
      push!(vn_chunks,        # store the previous chunk
            (vn_chunk, chunk_dim))
      vn_chunk = []           # initialise a new chunk
      chunk_dim = 0           # reset dimension counter
    end
    push!(vn_chunk, vns[i])       # put the current variable into the current chunk
    chunk_dim += l            # update dimension counter
  end
  push!(vn_chunks,            # push the last chunk
        (vn_chunk, chunk_dim))

  # Chunk-wise forward AD
  for (vn_chunk, chunk_dim) in vn_chunks
    # 1. Set dual part correspondingly
    @debug "set dual..."
    dim_count = 1
    for i = 1:vn_num
      range = getrange(vi, vns[i])
      l = length(range)
      vals = getval(vi, vns[i])
      if vns[i] in vn_chunk        # for each variable to compute gradient in this round
        for i = 1:l
          vi[range[i]] = ForwardDiff.Dual{Nothing, Float64, CHUNKSIZE}(realpart(vals[i]), SEEDS[dim_count])
          dim_count += 1      # count
        end
      else                    # for other varilables (no gradient in this round)
        for i = 1:l
          vi[range[i]] = ForwardDiff.Dual{Nothing, Float64, CHUNKSIZE}(realpart(vals[i]))
        end
      end
    end
    @debug "set dual done"

    # 2. Run model
    @debug "run model..."
    vi = runmodel(model, vi, spl)

    # 3. Collect gradient
    @debug "collect gradients from logp..."
    append!(grad, collect(dualpart(-getlogp(vi)))[1:chunk_dim])
  end

  if spl != nothing && haskey(spl.info, :grad_cache)
    spl.info[:grad_cache][θ_hash] = grad
  end

  vi.logp = realpart(vi.logp)

  grad
end

verifygrad(grad::Vector{Float64}) = begin
  if any(isnan.(grad)) || any(isinf.(grad))
    @warn("Numerical error has been found in gradients.")
    @warn("grad = $(grad)")
    false
  else
    true
  end
end

# Direct call of ForwardDiff.gradient; this is slow
gradient_slow(_vi::VarInfo, model::Function, spl::Union{Nothing, Sampler}) = begin

  vi = deepcopy(_vi)

  f(x::Vector) = begin
    vi[spl] = x
    -getlogp(runmodel(model, vi, spl))
  end

  g = x -> ForwardDiff.gradient(f, x::Vector, ForwardDiff.GradientConfig{min(length(x),CHUNKSIZE)}(x::Vector))

  g(vi[spl])
end

gradient_r(theta::AbstractVector{<:Real}, vi::VarInfo, model::Function) =
  gradient_r(theta, vi, model, nothing)
gradient_r(theta::AbstractVector{<:Real}, vi::Turing.VarInfo, model::Function, spl::Union{Nothing, Sampler}) = begin
  # Use Flux.Tracker to get gradient
  grad = Tracker.gradient(x -> (vi[spl] = x; -runmodel(model, vi, spl).logp), theta)
  # Clean tracked numbers
  # Numbers do not need to be tracked between two gradient calls
  vi.logp = vi.logp.data
  vi_spl = vi[spl]
  for i = 1:length(theta)
    vi_spl[i] = vi_spl[i].data
  end
  # Return non-tracked graident value
  return first(grad).data
end


text:
"""

Computing GMD contributions for the parts of the solute

"""

push!(LOAD_PATH,"../../../")

using M3GTools
using DelimitedFiles
using LaTeXStrings

function smooth(data)
  # smooth data
  n = length(data)
  y = copy(data)
  nav = 10
  for i in nav+1:n-nav
    av = 0.
    for j in i-nav:i+nav
      av = av + data[j]
    end
    y[i] = av / (2*nav)
  end
  for i in 1:nav
    y[i] = data[i]
    y[n-i] = data[n-i]
  end
  return y
end

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Data for my simulation (we only need the topoloy file here, might be PSF or GRO
topology = "./structure.psf"

# Read GMD file corresponding to solute contributions
println(" Reading GMD solute contribution file... ")
gmd_solute = readdlm("./gmd-GMD_ATOM_SOLUTE_CONTRIB.dat",comments=true,comment_char='#')

# The distance is the first column of the gmd file, and the total gmd is the second column:
d = gmd_solute[:,1]
gmd_total = gmd_solute[:,2]

gmd_total = smooth(gmd_total)

# Define the solute selection used for computation of the GMD here:
gmd_backbone = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and (not sidechain)")
gmd_aliphatic = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and aliphatic")
gmd_aromatic = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and aromatic")
gmd_polar = M3GTools.gmdget(topology,gmd_solute,data="protein",get="protein and sidechain and polar")

gmd_backbone = smooth(gmd_backbone)
gmd_aliphatic = smooth(gmd_aliphatic)
gmd_aromatic = smooth(gmd_aromatic)
gmd_polar = smooth(gmd_polar)

# Plot the results
plot(d,gmd_total,label="Total",xlim=[0,8],xlabel=L"\textrm{Distance / \AA}",ylabel=L"\textrm{MDDF}",linewidth=2)
plot!(d,gmd_backbone,label="Backbone",linewidth=2)
plot!(d,gmd_aliphatic,label="Aliphatic",linewidth=2)
plot!(d,gmd_aromatic,label="Aromatic",linewidth=2)
plot!(d,gmd_polar,label="Polar",linewidth=2)
plot!(fontsize=16)
savefig("./gmd_total.png")
println(" Created plot: ./gmd_total.png")


text:
module SGDOptim

using Compat
using Reexport
using ArrayViews
@reexport using EmpiricalRisks

import Base: start, next, done
import Base.LinAlg: axpy!

export
    # from ArrayViews
    views,

    # streams.jl
    SampleStream,
    GenericSampleSeq,
    sample_seq,
    minibatch_seq,

    # callback.jl
    simple_trace,
    gtcompare_trace,

    # sgd_std.jl
    sgd!,
    sgd

# source files

include("utils.jl")
include("streams.jl")
include("callback.jl")

include("sgdstd.jl")

end # module


text: Base.getindex(x::Array, ::Nothing) = x # allows for writing shorter, cleaner code at some places


text:
using Revise, LForaging
##

ispath("/home/beatriz/Documents/") ? (main_path ="/home/beatriz/Documents/") : (main_path = "/Users/dariosarra/Documents/Lab/Walton/LForaging")
fig_dir = joinpath(main_path,"Figures")
# pokes = DataFrame(CSV.read(joinpath(main_path,"Miscellaneous/Poke_data.csv")))#, DataFrame)
# pokes2 = DataFrame(CSV.read(joinpath(main_path,"Miscellaneous/New_Poke_data.csv")))#, DataFrame)
pokes = DataFrame(CSV.read(joinpath(main_path,"Miscellaneous/States_Poke_data.csv")))#, DataFrame)
pokes = preprocess_pokes(pokes)[:,[:MOUSE,:DATE,:SIDE,:KIND,:TRAVEL,:BOUT,:TRIAL,:BOUT_TRIAL,:POKE_TRIAL,:REWARD,:LEAVE,
    :IN, :OUT, :IN_TRIAL, :OUT_TRIAL, :IN_BOUT, :OUT_BOUT, :DURATION,:TIME]];
##
pokes[!,:BIN] = Int64.(round.(pokes.DURATION ./ 0.1))
filter!(r-> r.BIN >= 1, pokes)
pokes[!,:PatchRewRate] = get_rew_rate(pokes.REWARD,pokes.BIN, 0.1, 2)
transform!(groupby(pokes,[:MOUSE,:DATE]), [:REWARD,:BIN] => ((R,B) -> get_rew_rate(R,B, 0.1, 2)) => :PatchRewRate)
open_html_table(pokes[1:2000,:])

##
@df pokes density(:DURATION, group = :KIND, xrotation = 45)
savefig(joinpath(fig_dir,"Poke_Duration_density.png"))

@df pokes histogram(:DURATION, xrotation = 45, bins = 100)
savefig(joinpath(fig_dir,"Poke_Duration_histogram.png"))
##
fdf = filter(r-> r.SIDE != "travel", pokes)
fdf.MOUSE = categorical(fdf.MOUSE)
fdf.KIND = levels!(categorical(fdf.KIND),["poor", "medium", "rich"])
fdf.SIDE = categorical(fdf.SIDE)
fdf.TRAVEL = levels!(categorical(fdf.TRAVEL),["short", "long"])
fdf.REWARD = categorical(Bool.(fdf.REWARD))
f1 =  @formula(LEAVE ~ 1 + OUT_TRIAL+OUT_BOUT+DURATION+REWARD+TRAVEL+KIND + (1|MOUSE))
gm = fit(MixedModel, f1,fdf, Bernoulli())

##
pltdf = filter(r->r.DURATION <=2,fdf)
poke_rich, poke_rich_df =
    function_analysis(pltdf,:DURATION, cumulative_algorythm; grouping = :KIND, calc = :bootstrapping, xaxis = 0:0.1:2)
plot(poke_rich, ylabel = "Cumulative", xlabel = "Poke duration")
savefig(joinpath(fig_dir,"CumPokeRich.png"))
plt_trav = filter(r-> r.SIDE == "travel" && r.DURATION <=2, pokes)
dropmissing!(plt_trav)
poke_trav, poke_trav_df =
    function_analysis(plt_trav,:DURATION, cumulative_algorythm; grouping = :TRAVEL, calc = :bootstrapping, xaxis = 0:0.1:2)
plot(poke_trav, ylabel = "Cumulative", xlabel = "Poke duration")
savefig(joinpath(fig_dir,"CumPokeTrav.png"))
##
gd = groupby(filter(r->r.SIDE != "travel",pokes),[:MOUSE,:DATE,:TRIAL])
pre_trials = combine(gd) do dd
    DataFrame(
        DURATION_FULL = dd[end, :OUT_TRIAL],
        DURATION_SUM = round.(sum(dd.DURATION),digits =1),
        LAST_BOUT_FULL = dd[dd.LEAVE,:OUT_BOUT],
        LAST_BOUT_SUM = round.(sum(dd[dd.BOUT .== dd[dd.LEAVE,:BOUT],:DURATION]),digits =1),
        REWARDS = maximum(dd.BOUT),
        LAST_REWARD = Bool(dd[end,:REWARD]),
        SIDE = dd[1,:SIDE],
        KIND = dd[1,:KIND]
    )
end
open_html_table(pre_trials)

gd_tr = groupby(filter(r->r.SIDE == "travel",pokes),[:MOUSE,:DATE,:TRIAL])
trials_tr = combine(gd_tr) do dd
    DataFrame(
        TRAVEL_FULL = dd[end, :OUT_TRIAL],
        TRAVEL_SUM = round.(sum(dd.DURATION),digits = 1),
        TRAVEL_KIND = dd[1,:KIND]
    )
end
trials_tr.TRIAL = trials_tr.TRIAL .+1

trials = leftjoin(pre_trials, trials_tr, on = [:MOUSE,:DATE,:TRIAL])
open_html_table(pre_trials)
##
bout_sum_rich, bout_sum_rich_df =
    function_analysis(trials,:LAST_BOUT_SUM, cumulative_algorythm; grouping = :KIND, calc = :bootstrapping, xaxis = 0:0.5:10)
plot(bout_sum_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Last Bout duration (SUM)")
savefig(joinpath(fig_dir,"CumBoutSumRich.png"))
bout_full_rich, bout_full_rich_df =
    function_analysis(trials,:LAST_BOUT_FULL, cumulative_algorythm; grouping = :KIND, calc = :bootstrapping, xaxis = 0:1:40)
plot(bout_full_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Last Bout duration (FULL)")
savefig(joinpath(fig_dir,"CumBoutFullRich.png"))
##
trv_pltdf = filter(r->!ismissing(r.TRAVEL_KIND),trials)
bout_sum_rich, bout_sum_rich_df =
    function_analysis(trv_pltdf,:LAST_BOUT_SUM, cumulative_algorythm; grouping = :TRAVEL_KIND, calc = :bootstrapping, xaxis = 0:0.5:10)
plot(bout_sum_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Last Bout duration (SUM)")
savefig(joinpath(fig_dir,"CumBoutSumTrav.png"))
bout_full_rich, bout_full_rich_df =
    function_analysis(trv_pltdf,:LAST_BOUT_FULL, cumulative_algorythm; grouping = :TRAVEL_KIND, calc = :bootstrapping, xaxis = 0:1:40)
plot(bout_full_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Last Bout duration (FULL)")
savefig(joinpath(fig_dir,"CumBoutFullTrav.png"))
##
trial_sum_rich, trial_sum_rich_df =
    function_analysis(trials,:DURATION_SUM, cumulative_algorythm; grouping = :KIND, calc = :bootstrapping, xaxis = 0:1:30)
plot(trial_sum_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Trial duration (SUM)")
savefig(joinpath(fig_dir,"CumTrialSumRich.png"))
trial_full_rich, trial_full_rich_df =
    function_analysis(trials,:DURATION_FULL, cumulative_algorythm; grouping = :KIND, calc = :bootstrapping, xaxis = 0:5:100)
plot(trial_full_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Trial duration (FULL)")
savefig(joinpath(fig_dir,"CumTrialFullRich.png"))
##
trv_pltdf = filter(r->!ismissing(r.TRAVEL_KIND),trials)
trial_sum_rich, trial_sum_rich_df =
    function_analysis(trv_pltdf,:DURATION_SUM, cumulative_algorythm; grouping = :TRAVEL_KIND, calc = :bootstrapping, xaxis = 0:1:30)
plot(trial_sum_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Trial duration (SUM)")
savefig(joinpath(fig_dir,"CumTrialSumTrav.png"))
trial_full_rich, trial_full_rich_df =
    function_analysis(trv_pltdf,:DURATION_FULL, cumulative_algorythm; grouping = :TRAVEL_KIND, calc = :bootstrapping, xaxis = 0:5:100)
plot(trial_full_rich, legend = :bottomright, fillalpha = 0.3, ylabel = "Cumulative", xlabel = "Trial duration (FULL)")
savefig(joinpath(fig_dir,"CumTrialFullTrav.png"))
##
gs0 = trials
gs0.DURATION_FULL = round.(gs0.DURATION_FULL)
filter!(r->r.DURATION_FULL <=60, gs0)
gs1 = combine(groupby(gs0,[:MOUSE,:DURATION_FULL,:KIND]),
    :LAST_BOUT_FULL	.=> [mean,sem]
)

gs2 = combine(groupby(gs1,[:DURATION_FULL,:KIND]),
    :LAST_BOUT_FULL_mean .=> [mean,sem] .=> [:LAST_BOUT_FULL_mean,:LAST_BOUT_FULL_sem]
)

sort!(gs2,:DURATION_FULL)
@df gs2 plot(:DURATION_FULL, :LAST_BOUT_FULL_mean, ribbon = :LAST_BOUT_FULL_sem, group = :KIND,
    xlabel = "Trial Duration (Full)", ylabel = "Last Bout Duration (Full)")
savefig(joinpath(fig_dir,"Full(BoutxTrav)_Kind.png"))
##
gs0 = trials
gs0.DURATION_FULL = round.(gs0.DURATION_FULL)
filter!(r->r.DURATION_FULL <=60, gs0)
gs1 = combine(groupby(gs0,[:MOUSE,:DURATION_FULL,:TRAVEL_KIND]),
    :LAST_BOUT_FULL	.=> [mean,sem]
)

gs2 = combine(groupby(gs1,[:DURATION_FULL,:TRAVEL_KIND]),
    :LAST_BOUT_FULL_mean .=> [mean,sem] .=> [:LAST_BOUT_FULL_mean,:LAST_BOUT_FULL_sem]
)

sort!(gs2,:DURATION_FULL)
filter!(r->!ismissing(r.TRAVEL_KIND), gs2)
@df gs2 plot(:DURATION_FULL, :LAST_BOUT_FULL_mean, ribbon = :LAST_BOUT_FULL_sem, group = :TRAVEL_KIND,
    xlabel = "Trial Duration (Full)", ylabel = "Last Bout Duration (Full)")
savefig(joinpath(fig_dir,"Full(BoutxTrav)_Travel.png"))
##
sum(trials.LAST_REWARD)/nrow(trials)


text:
# Automatically generated using Clang.jl

# Skipping MacroDefinition: SC_CALC_VERSION ( major , minor , patchlevel ) ( ( ( major ) * 1000 + ( minor ) ) * 1000 + ( patchlevel ) )
# Skipping MacroDefinition: SC_GCC_VERSION SC_CALC_VERSION ( __GNUC__ , __GNUC_MINOR__ , __GNUC_PATCHLEVEL__ )
# Skipping MacroDefinition: SC_EXTERN_C_BEGIN void sc_extern_c_hack_3 ( void )
# Skipping MacroDefinition: SC_EXTERN_C_END void sc_extern_c_hack_4 ( void )

const SC_EPS = 2.220446049250313e-16
const SC_1000_EPS = "1000.*2.220446049250313e-16"

# Skipping MacroDefinition: SC_NOOP ( ) ( ( void ) ( 0 ) )
# Skipping MacroDefinition: SC_ABORT ( s ) sc_abort_verbose ( __FILE__ , __LINE__ , ( s ) )
# Skipping MacroDefinition: SC_ABORT_NOT_REACHED ( ) SC_ABORT ( "Unreachable code" )
# Skipping MacroDefinition: SC_CHECK_ABORT ( q , s ) ( ( q ) ? ( void ) 0 : SC_ABORT ( s ) )
# Skipping MacroDefinition: SC_CHECK_MPI ( r ) SC_CHECK_ABORT ( ( r ) == sc_MPI_SUCCESS , "MPI error" )
# Skipping MacroDefinition: SC_CHECK_ZLIB ( r ) SC_CHECK_ABORT ( ( r ) == Z_OK , "zlib error" )
# Skipping MacroDefinition: SC_ABORTF ( fmt , ... ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_CHECK_ABORTF ( q , fmt , ... ) ( ( q ) ? ( void ) 0 : SC_ABORTF ( fmt , __VA_ARGS__ ) )
# Skipping MacroDefinition: SC_ABORT1 ( fmt , a ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) )
# Skipping MacroDefinition: SC_ABORT2 ( fmt , a , b ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) , ( b ) )
# Skipping MacroDefinition: SC_ABORT3 ( fmt , a , b , c ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) , ( b ) , ( c ) )
# Skipping MacroDefinition: SC_ABORT4 ( fmt , a , b , c , d ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) )
# Skipping MacroDefinition: SC_ABORT5 ( fmt , a , b , c , d , e ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) , ( e ) )
# Skipping MacroDefinition: SC_ABORT6 ( fmt , a , b , c , d , e , f ) sc_abort_verbosef ( __FILE__ , __LINE__ , ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) , ( e ) , ( f ) )
# Skipping MacroDefinition: SC_CHECK_ABORT1 ( q , fmt , a ) ( ( q ) ? ( void ) 0 : SC_ABORT1 ( ( fmt ) , ( a ) ) )
# Skipping MacroDefinition: SC_CHECK_ABORT2 ( q , fmt , a , b ) ( ( q ) ? ( void ) 0 : SC_ABORT2 ( ( fmt ) , ( a ) , ( b ) ) )
# Skipping MacroDefinition: SC_CHECK_ABORT3 ( q , fmt , a , b , c ) ( ( q ) ? ( void ) 0 : SC_ABORT3 ( ( fmt ) , ( a ) , ( b ) , ( c ) ) )
# Skipping MacroDefinition: SC_CHECK_ABORT4 ( q , fmt , a , b , c , d ) ( ( q ) ? ( void ) 0 : SC_ABORT4 ( ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) ) )
# Skipping MacroDefinition: SC_CHECK_ABORT5 ( q , fmt , a , b , c , d , e ) ( ( q ) ? ( void ) 0 : SC_ABORT5 ( ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) , ( e ) ) )
# Skipping MacroDefinition: SC_CHECK_ABORT6 ( q , fmt , a , b , c , d , e , f ) ( ( q ) ? ( void ) 0 : SC_ABORT6 ( ( fmt ) , ( a ) , ( b ) , ( c ) , ( d ) , ( e ) , ( f ) ) )
# Skipping MacroDefinition: SC_ASSERT ( c ) SC_NOOP ( )
# Skipping MacroDefinition: SC_EXECUTE_ASSERT_FALSE ( expression ) do { ( void ) ( expression ) ; } while ( 0 )
# Skipping MacroDefinition: SC_EXECUTE_ASSERT_TRUE ( expression ) do { ( void ) ( expression ) ; } while ( 0 )
# Skipping MacroDefinition: SC_ALLOC ( t , n ) ( t * ) sc_malloc ( sc_package_id , ( n ) * sizeof ( t ) )
# Skipping MacroDefinition: SC_ALLOC_ZERO ( t , n ) ( t * ) sc_calloc ( sc_package_id , ( size_t ) ( n ) , sizeof ( t ) )
# Skipping MacroDefinition: SC_REALLOC ( p , t , n ) ( t * ) sc_realloc ( sc_package_id , ( p ) , ( n ) * sizeof ( t ) )
# Skipping MacroDefinition: SC_STRDUP ( s ) sc_strdup ( sc_package_id , ( s ) )
# Skipping MacroDefinition: SC_FREE ( p ) sc_free ( sc_package_id , ( p ) )
# Skipping MacroDefinition: SC_ALIGN_UP ( x , n ) ( ( ( n ) <= 0 ) ? ( x ) : ( ( x ) + ( n ) - 1 ) / ( n ) * ( n ) )
# Skipping MacroDefinition: SC_ARG_ALIGN ( p , t , n ) SC_NOOP ( )
# Skipping MacroDefinition: SC_ATTR_ALIGN ( n ) __attribute__ ( ( aligned ( n ) ) )
# Skipping MacroDefinition: SC_BZERO ( p , n ) ( ( void ) memset ( ( p ) , 0 , ( n ) * sizeof ( * ( p ) ) ) )
# Skipping MacroDefinition: SC_MIN ( a , b ) ( ( ( a ) < ( b ) ) ? ( a ) : ( b ) )
# Skipping MacroDefinition: SC_MAX ( a , b ) ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) )
# Skipping MacroDefinition: SC_SQR ( a ) ( ( a ) * ( a ) )
# Skipping MacroDefinition: SC_LOG2_8 ( x ) ( sc_log2_lookup_table [ ( x ) ] )
# Skipping MacroDefinition: SC_LOG2_16 ( x ) ( ( ( x ) > 0xff ) ? ( SC_LOG2_8 ( ( x ) >> 8 ) + 8 ) : SC_LOG2_8 ( x ) )
# Skipping MacroDefinition: SC_LOG2_32 ( x ) ( ( ( x ) > 0xffff ) ? ( SC_LOG2_16 ( ( x ) >> 16 ) ) + 16 : SC_LOG2_16 ( x ) )
# Skipping MacroDefinition: SC_LOG2_64 ( x ) ( ( ( x ) > 0xffffffffLL ) ? ( SC_LOG2_32 ( ( x ) >> 32 ) ) + 32 : SC_LOG2_32 ( x ) )
# Skipping MacroDefinition: SC_ROUNDUP2_32 ( x ) ( ( ( x ) <= 0 ) ? 0 : ( 1 << ( SC_LOG2_32 ( ( x ) - 1 ) + 1 ) ) )
# Skipping MacroDefinition: SC_ROUNDUP2_64 ( x ) ( ( ( x ) <= 0 ) ? 0 : ( 1LL << ( SC_LOG2_64 ( ( x ) - 1LL ) + 1 ) ) )

const SC_LC_GLOBAL = 1
const SC_LC_NORMAL = 2
const SC_LP_DEFAULT = -1
const SC_LP_ALWAYS = 0
const SC_LP_TRACE = 1
const SC_LP_DEBUG = 2
const SC_LP_VERBOSE = 3
const SC_LP_INFO = 4
const SC_LP_STATISTICS = 5
const SC_LP_PRODUCTION = 6
const SC_LP_ESSENTIAL = 7
const SC_LP_ERROR = 8
const SC_LP_SILENT = 9
const SC_LP_THRESHOLD = SC_LP_INFO

# Skipping MacroDefinition: SC_GEN_LOG ( package , category , priority , s ) ( ( priority ) < SC_LP_THRESHOLD ? ( void ) 0 : sc_log ( __FILE__ , __LINE__ , ( package ) , ( category ) , ( priority ) , ( s ) ) )
# Skipping MacroDefinition: SC_GLOBAL_LOG ( p , s ) SC_GEN_LOG ( sc_package_id , SC_LC_GLOBAL , ( p ) , ( s ) )
# Skipping MacroDefinition: SC_LOG ( p , s ) SC_GEN_LOG ( sc_package_id , SC_LC_NORMAL , ( p ) , ( s ) )
# Skipping MacroDefinition: SC_GEN_LOGF ( package , category , priority , fmt , ... ) ( ( priority ) < SC_LP_THRESHOLD ? ( void ) 0 : sc_logf ( __FILE__ , __LINE__ , ( package ) , ( category ) , ( priority ) , ( fmt ) , __VA_ARGS__ ) )
# Skipping MacroDefinition: SC_GLOBAL_LOGF ( p , fmt , ... ) SC_GEN_LOGF ( sc_package_id , SC_LC_GLOBAL , ( p ) , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_LOGF ( p , fmt , ... ) SC_GEN_LOGF ( sc_package_id , SC_LC_NORMAL , ( p ) , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_TRACE ( s ) SC_GLOBAL_LOG ( SC_LP_TRACE , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_LDEBUG ( s ) SC_GLOBAL_LOG ( SC_LP_DEBUG , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_VERBOSE ( s ) SC_GLOBAL_LOG ( SC_LP_VERBOSE , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_INFO ( s ) SC_GLOBAL_LOG ( SC_LP_INFO , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_STATISTICS ( s ) SC_GLOBAL_LOG ( SC_LP_STATISTICS , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_PRODUCTION ( s ) SC_GLOBAL_LOG ( SC_LP_PRODUCTION , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_LERROR ( s ) SC_GLOBAL_LOG ( SC_LP_ERROR , ( s ) )
# Skipping MacroDefinition: SC_GLOBAL_TRACEF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_TRACE , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_LDEBUGF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_DEBUG , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_VERBOSEF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_VERBOSE , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_INFOF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_INFO , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_STATISTICSF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_STATISTICS , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_PRODUCTIONF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_PRODUCTION , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_ESSENTIALF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_ESSENTIAL , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_GLOBAL_LERRORF ( fmt , ... ) SC_GLOBAL_LOGF ( SC_LP_ERROR , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_TRACE ( s ) SC_LOG ( SC_LP_TRACE , ( s ) )
# Skipping MacroDefinition: SC_LDEBUG ( s ) SC_LOG ( SC_LP_DEBUG , ( s ) )
# Skipping MacroDefinition: SC_VERBOSE ( s ) SC_LOG ( SC_LP_VERBOSE , ( s ) )
# Skipping MacroDefinition: SC_INFO ( s ) SC_LOG ( SC_LP_INFO , ( s ) )
# Skipping MacroDefinition: SC_STATISTICS ( s ) SC_LOG ( SC_LP_STATISTICS , ( s ) )
# Skipping MacroDefinition: SC_PRODUCTION ( s ) SC_LOG ( SC_LP_PRODUCTION , ( s ) )
# Skipping MacroDefinition: SC_ESSENTIAL ( s ) SC_LOG ( SC_LP_ESSENTIAL , ( s ) )
# Skipping MacroDefinition: SC_LERROR ( s ) SC_LOG ( SC_LP_ERROR , ( s ) )
# Skipping MacroDefinition: SC_TRACEF ( fmt , ... ) SC_LOGF ( SC_LP_TRACE , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_LDEBUGF ( fmt , ... ) SC_LOGF ( SC_LP_DEBUG , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_VERBOSEF ( fmt , ... ) SC_LOGF ( SC_LP_VERBOSE , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_INFOF ( fmt , ... ) SC_LOGF ( SC_LP_INFO , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_STATISTICSF ( fmt , ... ) SC_LOGF ( SC_LP_STATISTICS , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_PRODUCTIONF ( fmt , ... ) SC_LOGF ( SC_LP_PRODUCTION , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_ESSENTIALF ( fmt , ... ) SC_LOGF ( SC_LP_ESSENTIAL , ( fmt ) , __VA_ARGS__ )
# Skipping MacroDefinition: SC_LERRORF ( fmt , ... ) SC_LOGF ( SC_LP_ERROR , ( fmt ) , __VA_ARGS__ )

const sc_handler_t = Ptr{Cvoid}
const sc_log_handler_t = Ptr{Cvoid}
const sc_abort_handler_t = Ptr{Cvoid}

# Skipping MacroDefinition: sc_hash_rot ( x , k ) ( ( ( x ) << ( k ) ) | ( ( x ) >> ( 32 - ( k ) ) ) )
# Skipping MacroDefinition: sc_hash_mix ( a , b , c ) ( ( void ) ( a -= c , a ^= sc_hash_rot ( c , 4 ) , c += b , b -= a , b ^= sc_hash_rot ( a , 6 ) , a += c , c -= b , c ^= sc_hash_rot ( b , 8 ) , b += a , a -= c , a ^= sc_hash_rot ( c , 16 ) , c += b , b -= a , b ^= sc_hash_rot ( a , 19 ) , a += c , c -= b , c ^= sc_hash_rot ( b , 4 ) , b += a ) )
# Skipping MacroDefinition: sc_hash_final ( a , b , c ) ( ( void ) ( c ^= b , c -= sc_hash_rot ( b , 14 ) , a ^= c , a -= sc_hash_rot ( c , 11 ) , b ^= a , b -= sc_hash_rot ( a , 25 ) , c ^= b , c -= sc_hash_rot ( b , 16 ) , a ^= c , a -= sc_hash_rot ( c , 4 ) , b ^= a , b -= sc_hash_rot ( a , 14 ) , c ^= b , c -= sc_hash_rot ( b , 24 ) ) )
# Skipping MacroDefinition: SC_ARRAY_IS_OWNER ( a ) ( ( a ) -> byte_alloc >= 0 )
# Skipping MacroDefinition: SC_ARRAY_BYTE_ALLOC ( a ) ( ( size_t ) ( SC_ARRAY_IS_OWNER ( a ) ? ( a ) -> byte_alloc : - ( ( a ) -> byte_alloc + 1 ) ) )
# Skipping MacroDefinition: sc_array_new_size ( s , c ) ( sc_array_new_count ( ( s ) , ( c ) ) )

const sc_hash_function_t = Ptr{Cvoid}
const sc_equal_function_t = Ptr{Cvoid}
const sc_hash_foreach_t = Ptr{Cvoid}
const sc_array_t = sc_array

# Skipping Typedef: CXType_FunctionProto size_t

struct sc_mstamp
    elem_size::Csize_t
    per_stamp::Csize_t
    stamp_size::Csize_t
    cur_snext::Csize_t
    current::Cstring
    remember::sc_array_t
end

const sc_mstamp_t = sc_mstamp

struct sc_mempool
    elem_size::Csize_t
    elem_count::Csize_t
    zero_and_persist::Cint
    mstamp::sc_mstamp_t
    freed::sc_array_t
end

const sc_mempool_t = sc_mempool

struct sc_link
    data::Ptr{Cvoid}
    next::Ptr{sc_link}
end

const sc_link_t = sc_link

struct sc_list
    elem_count::Csize_t
    first::Ptr{sc_link_t}
    last::Ptr{sc_link_t}
    allocator_owned::Cint
    allocator::Ptr{sc_mempool_t}
end

const sc_list_t = sc_list

struct sc_hash
    elem_count::Csize_t
    slots::Ptr{sc_array_t}
    user_data::Ptr{Cvoid}
    hash_fn::sc_hash_function_t
    equal_fn::sc_equal_function_t
    resize_checks::Csize_t
    resize_actions::Csize_t
    allocator_owned::Cint
    allocator::Ptr{sc_mempool_t}
end

const sc_hash_t = sc_hash

struct sc_hash_array_data
    pa::Ptr{sc_array_t}
    hash_fn::sc_hash_function_t
    equal_fn::sc_equal_function_t
    user_data::Ptr{Cvoid}
    current_item::Ptr{Cvoid}
end

const sc_hash_array_data_t = sc_hash_array_data

struct sc_hash_array
    a::sc_array_t
    internal_data::sc_hash_array_data_t
    h::Ptr{sc_hash_t}
end

const sc_hash_array_t = sc_hash_array

struct sc_recycle_array
    elem_count::Csize_t
    a::sc_array_t
    f::sc_array_t
end

const sc_recycle_array_t = sc_recycle_array

@cenum sc_io_error_t::Int32 begin
    SC_IO_ERROR_NONE = 0
    SC_IO_ERROR_FATAL = -1
    SC_IO_ERROR_AGAIN = -2
end

@cenum sc_io_mode_t::UInt32 begin
    SC_IO_MODE_WRITE = 0
    SC_IO_MODE_APPEND = 1
    SC_IO_MODE_LAST = 2
end

@cenum sc_io_encode_t::UInt32 begin
    SC_IO_ENCODE_NONE = 0
    SC_IO_ENCODE_LAST = 1
end

@cenum sc_io_type_t::UInt32 begin
    SC_IO_TYPE_BUFFER = 0
    SC_IO_TYPE_FILENAME = 1
    SC_IO_TYPE_FILEFILE = 2
    SC_IO_TYPE_LAST = 3
end

struct sc_io_sink
    iotype::sc_io_type_t
    mode::sc_io_mode_t
    encode::sc_io_encode_t
    buffer::Ptr{sc_array_t}
    buffer_bytes::Csize_t
    file::Ptr{FILE}
    bytes_in::Csize_t
    bytes_out::Csize_t
end

const sc_io_sink_t = sc_io_sink

struct sc_io_source
    iotype::sc_io_type_t
    encode::sc_io_encode_t
    buffer::Ptr{sc_array_t}
    buffer_bytes::Csize_t
    file::Ptr{FILE}
    bytes_in::Csize_t
    bytes_out::Csize_t
    mirror::Ptr{sc_io_sink_t}
    mirror_buffer::Ptr{sc_array_t}
end

const sc_io_source_t = sc_io_source


text:
"""
    Code to loading mask shapes for CCF
Author: Eric Ford
Created: September 2020
"""

# constants

default_gaussian_ccf_σ = 5000.0
default_gaussian_ccf_truncation_scale_factor = 2*(sqrt(2*log(2)))
default_supergaussian_ccf_fwhm = 448.0*2*sqrt(2*log(2)) #*4.5
default_supergaussian_ccf_exponent = 1.3
default_supergaussian_ccf_truncation_scale_factor = 2
default_gaussian_mixture_ccf_truncation_Δv = 20000.0

# Can be specialized for specific instruments
default_ccf_mask_v_width(inst::AbstractInstrument) = 1e3

"A struct implementing a specific mask shapes should be a subtype of AbstractCCFMaskShape."
abstract type AbstractCCFMaskShape  end
export AbstractCCFMaskShape

include("tophat.jl")
export TopHatCCFMask
include("gaussian.jl")
export GaussianCCFMask
include("halfcos.jl")
export CosCCFMask
include("supergaussian.jl")
export SuperGaussianCCFMask

# TODO: Fix before adding
#include("gaussian_mixture.jl")
#export GaussianMixtureCCFMask

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
import Phylo.API: _src, _dst, _setsrc!, _setdst!, _getlength

"""
    Branch

    A directed branch connecting two AbstractNodes of phylogenetic tree
"""
mutable struct Branch{T} <: AbstractBranch
    source::T
    destination::T
    length::Float64

    function Branch(source::T, destination::T, length::Float64) where T
        length >= 0.0 || isnan(length) ||
            error("Branch length must be positive or NaN (no recorded length)")
        new{T}(source, destination, length)
    end
end

const SimpleBranch = Branch{Int}

_src(branch::Branch) = branch.source
_dst(branch::Branch) = branch.destination
_setsrc!(branch::Branch{T}, source::T) where T = branch.source = source
_setdst!(branch::Branch{T}, destination::T) where T = branch.destination = destination
_getlength(branch::Branch) = branch.length

function checkbranch(id::Int, branch::Branch, tree::AbstractTree)
    return id > 0 &&
        src(branch) != dst(branch) &&
        !haskey(getbranches(tree), id) &&
        haskey(getnodes(tree), src(branch)) &&
        haskey(getnodes(tree), dst(branch)) &&
        !hasinbound(getnodes(tree)[dst(branch)]) &&
        outboundspace(getnodes(tree)[src(branch)])
end


text:
include("./basis_types.jl")
include("./GLMFit.jl")
using .GLMFit
using .GLMFit.GLMMetrics
import .GLMFit.GLMPriors
import .RelayUtils
using LinearAlgebra, Statistics, Random
const VecVec{T} = AbstractVector{<:AbstractVector{T}}
# ============================================================================ #
struct Predictor{T<:AbstractBasis}
    ts::Vector{Float64}
    evt::Vector{Float64}
    basis::T
end

Predictor(x::Vector{Float64}) = Predictor(Float64[], x, NullBasis())

function Predictor(p::Predictor{L}, basis::T) where {L<:AbstractBasis, T<:AbstractBasis}
    T == L && return p
    return Predictor(p.ts, p.evt, basis)
end

@inline nrow(p::Predictor) = length(p.evt)
@inline ncol(p::Predictor) = nbasis(p.basis)
@inline Base.size(p::Predictor) = (nrow(p), ncol(p))

@inline generate(p::Predictor) = project(p.basis, p.ts, p.evt)

@inline get_basis(p::Predictor) = p.basis

# ============================================================================ #
struct PredictorSet
    keys::Vector{Symbol}
    d::Dict{Symbol, Predictor}
end

PredictorSet() = PredictorSet(Vector{Symbol}(), Dict{Symbol, Predictor}())

Base.IteratorSize(::PredictorSet) = Base.HasLength()
Base.IteratorEltype(::PredictorSet) = Base.HasEltype()
Base.eltype(p::PredictorSet) = Predictor
Base.length(p::PredictorSet) = length(p.keys)
Base.iterate(p::PredictorSet) = iterate(p, 1)
function Base.iterate(p::PredictorSet, k::Integer)
    k > length(p) && return nothing
    return p.d[p.keys[k]], k+1
end

Base.getindex(ps::PredictorSet, k::Symbol) = ps.d[k]
function Base.setindex!(ps::PredictorSet, p::Predictor, k::Symbol)
    !in(k, ps.keys) && push!(ps.keys, k)
    ps.d[k] = p
    return ps
end

@inline ncol(ps::PredictorSet) = sum(ncol, ps) + 1

function generate(ps::PredictorSet)
    nc = ncol(ps)
    nr = nrow(first(ps))
    @assert(all(x -> nrow(x) == nr, ps), "Predictors do not have a common number of observations (rows)!")

    # NOTE: the y-intercept is addeded here as the first column
    dm = ones(nr, nc)
    kl = 2
    for p in ps
        ke = kl + ncol(p) - 1
        dm[:,kl:ke] .= generate(p)
        kl = ke + 1
    end
    return dm
end

@inline get_basis(p::PredictorSet, k::Symbol) = get_basis(p[k])
# ============================================================================ #
abstract type AbstractPrior end
struct NullPrior <: AbstractPrior end
struct RidgePrior <: AbstractPrior end
struct SmoothingPrior <: AbstractPrior end
struct LassoPrior <: AbstractPrior end

generate(::Type{NullPrior}, flen::NTuple{N,Int}, lmb::NTuple{N,T}) where {N,T<:AbstractVector{<:Real}} = GLMPriors.null_prior(flen)
generate(::Type{RidgePrior}, flen::NTuple{N,Int}, lmb::NTuple{N,T}) where {N,T<:AbstractVector{<:Real}} = GLMPriors.ridge_prior(flen, lmb)
generate(::Type{SmoothingPrior}, flen::NTuple{N,Int}, lmb::NTuple{N,T}) where {N,T<:AbstractVector{<:Real}} = GLMPriors.smoothing_prior(flen, lmb)
generate(::Type{LassoPrior}, flen::NTuple{N,Int}, lmb::NTuple{N,T}) where {N,T<:AbstractVector{<:Real}} = GLMPriors.lasso_prior(flen, lmb)

# ============================================================================ #
abstract type AbstractGLM end

struct GLM <: AbstractGLM #{T<:AbstractPrior,R<:Real,L<:Real}
    predictors::PredictorSet
    response::Vector{Bool}
    # lambda::Vector{Vector{L}}
end
struct RegGLM{T<:AbstractPrior} <: AbstractGLM
    glm::GLM
    lambda::Vector{Vector{Float64}}
end
"""
lambda = [1:10,[]] # -> apply prior to only the first filter, ignore 2nd
lambda = [1:10] # -> apply lambda to all filters
"""
# function GLM(ps::PredictorSet, y::Vector{<:Real})
#     return GLM(ps, y)
# end
function GLM(ps::PredictorSet, y::Vector{Bool}, prior::Type{P}, lambda::VecVec{N}) where {P<:AbstractPrior,N<:Real}
    return RegGLM{P}(GLM(ps, y), map(collect, lambda))
end
function GLM(glm::GLM)
    return GLM(glm.predictors, glm.response)
end
function GLM(glm::GLM, prior::Type{P}, lambda::VecVec{N}) where {P<:AbstractPrior,N<:Real}
    return RegGLM{P}(glm, map(collect, lambda))
end
@inline get_basis(glm::GLM, name::Symbol) = get_basis(glm.predictors, name)
@inline get_predictor(glm::GLM, name::Symbol) = glm.predictors[name]
@inline get_keys(glm::GLM) = glm.predictors.keys
@inline predictors(glm::GLM) = glm.predictors
@inline response(glm::GLM) = glm.response

@inline get_basis(rg::RegGLM, name::Symbol) = get_basis(rg.glm.predictors, name)
@inline get_predictor(rg::RegGLM, name::Symbol) = rg.glm.predictors[name]
@inline get_keys(rg::RegGLM) = rg.glm.predictors.keys
@inline predictors(rg::RegGLM) = rg.glm.predictors
@inline response(rg::RegGLM) = rg.glm.response

@inline function generate_prior(rg::RegGLM{P}) where P
    # lambda = [1:10,[]] # -> apply prior to only the first filter, ignore 2nd
    # lambda = [1:10] # -> apply lambda to all filters
    # if lambda contains a single set apply it to all filters
    flen = length(rg.lambda) == 1 ? sum : map
    return generate(P, tuple(flen(ncol, rg.glm.predictors)...), tuple(rg.lambda...))
end
# ============================================================================ #
struct GLMFilter
    coef::Vector{Float64}
    error::Vector{Float64}
end
struct GLMResult{T<:PerformanceMetric}
    filters::Dict{Symbol,GLMFilter}
    lambda::Dict{Symbol,Float64}
    coef::Vector{Float64}
    metric::T
    nlli::Float64
    null_nlli::Float64
    converged::Bool
end
function GLMResult(coef::Vector{Float64}, metric::T, li::Real, nl::Real, convg::Bool) where T <: PerformanceMetric
    return GLMResult{T}(Dict{Symbol,GLMFilter}(), Dict{Symbol,Float64}(), coef, metric, li, nl, convg)
end

@inline function Base.setindex!(r::GLMResult, f::GLMFilter, k::Symbol)
    r.filters[k] = f
    return r
end
@inline get_coef(r::GLMResult, k::Symbol) = r.filters[k].coef
@inline get_error(r::GLMResult, k::Symbol) = r.filters[k].error
@inline metric(r::GLMResult) = r.metric
@inline function set_lambda!(r::GLMResult, k::Symbol, v::Float64)
    r.lambda[k] = v
    return r
end
@inline get_lambda(r::GLMResult, k::Symbol) = r.lambda[k]
# ============================================================================ #
wasrelayed(pre::Vector{Float64}, post::Vector{Float64}) = RelayUtils.relay_status(pre, post)
# ============================================================================ #
function cross_validate_impl(::Type{T}, ::Type{D}, ::Type{A}, glm::GLM, nfold::Integer, shfl::Bool) where {T<:PerformanceMetric,D<:Distribution,A<:ActivationFunction}
    bs = map(x -> getfield(get_basis(x), :bin_size), predictors(glm))
    @assert(all(x -> x==bs[1], bs), "Bin sizes *MUST* be consistent across predictors")

    # construct GLM data container from the old interface
    data = GLMFit.GLM(D, A, generate(predictors(glm)), glm.response, bs[1])

    # initial guess for model parameters is just the "relayed-spike triggered average"
    x0 = data.x' * data.y ./ sum(data.y)

    # old cross_validate() interface does the actual fitting
    return GLMFit.cross_validate(T, data, x0, nfold, shfl), data
end
#----------------------------------------------------------------------------- #
function cross_validate_impl(::Type{T}, ::Type{D}, ::Type{A}, glm::RegGLM, nfold::Integer, shfl::Bool) where {T<:PerformanceMetric,D<:Distribution,A<:ActivationFunction}
    bs = map(x -> getfield(get_basis(x), :bin_size), predictors(glm))
    @assert(all(x -> x==bs[1], bs), "Bin sizes *MUST* be consistent across predictors")

    prior = generate_prior(glm)
    pr = GLMPriors.allocate_prior(prior)

    # construct GLM data container from the old interface
    data = GLMFit.RegularizedGLM(D, A, generate(predictors(glm)), glm.glm.response, bs[1], pr)

    # initial guess for model parameters is just the "relayed-spike triggered average"
    x0 = data.glm.x' * data.glm.y ./ sum(data.glm.y)

    # old cross_validate() interface does the actual fitting
    return GLMFit.cross_validate(T, data, prior, x0, nfold, shfl), data.glm
end
#----------------------------------------------------------------------------- #
function cross_validate(::Type{T}, ::Type{D}, ::Type{A}, glm::GLMType; nfold=10, shuffle_design=true) where {T<:PerformanceMetric,D<:Distribution,A<:ActivationFunction,GLMType<:AbstractGLM}

    res, data = cross_validate_impl(T, D, A, glm, nfold, shuffle_design)

    coef = coefficients(res)

    se = try
        # standard error of the optimization
        sqrt.(diag(pinv(GLMFit.hessian(data, coef))))
    catch err
        @warn("Failed to calculate coefficient std-error: \"$(typeof(err))\"")
        fill(NaN, length(coef))
    end

    # cor(coef, x0) > 0.9 && @warn("Fitting may have failed: cor(coef, x0) = $(cor(coef, x0))")

    result = GLMResult(
        coef,
        GLMFit.metric(res),
        GLMFit.nlli(res),
        GLMModels.null_nlli(data) / length(data.y),
        res.nconvg == res.ntotal
        )

    klambda = 1
    kl = 2
    lm = lambda(res)
    for k in get_keys(glm)
        basis = get_basis(glm, k)
        ke = kl + nbasis(basis) - 1
        result[k] = GLMFilter(basis * coef[kl:ke], basis * se[kl:ke])
        if klambda > length(lm)
            set_lambda!(result, k, lm[end])
        else
            set_lambda!(result, k, lm[klambda])
        end
        kl = ke + 1
        klambda += 1
    end

    return result
end
# ============================================================================ #
function jensen_shannon(::Type{A}, glm::GLM, result::GLMResult) where A <: ActivationFunction
    bs = get_basis(glm.predictors[:ff]).bin_size
    x = generate(glm.predictors)

    l = x * result.coef

    pred, _ = GLMModels.activate(A, l, bs)

    krel = findall(isequal(1), glm.response)
    kfail = findall(iszero, glm.response)

    d = RelayUtils.jensen_shannon(pred, glm.response)
    d2 = RelayUtils.jensen_shannon(l, glm.response)

    return d, d2
end
# ============================================================================ #


text:
try
    @eval using Revise
catch err
    @warn "Error starting Revise" exception=err
end

max_stars_repo_name: oxinabox/dotfiles

max_issues_repo_name: oxinabox/dotfiles

max_forks_repo_name: oxinabox/dotfiles


text:
using OptimalEstimationProject
using LinearAlgebra
using AstroTime
using MATLAB

# Spacecraft simulator
t0      = TAIEpoch("2021-08-08T12:00:00")
scsim   = SpacecraftSim(t0.second)

# Initial state error variance
σr      = 0.1   # [km]
σv      = 0.01  # [km / s]
σm      = 1e-6  # [kg]
sqrtP0  = Diagonal([σr, σr, σr, σv, σv, σv, σm])
P0      = sqrtP0.^2

# Initial state estimate
(y0true, us) = OptimalEstimationProject.GetStateAndControl(scsim, t0.second)
xhat0   = y0true[1:7] + sqrtP0*randn(7)

# Measurement time stamps
Δt      = 60.0        # [sec]
gpsΔt   = 15*60     # [sec]
ts      = range(t0.second + Δt; step = Δt, stop = t0.second + scsim.ts[2]*86400)

# Measurement statistics
σρ      = 10.0e-3 # [km]     Pseudorange noise standard deviation
σr      = 5.0e-3  # [km]     GPS broadcast ephemeris standard deviation
σa      = 1.0e-6  # [km/s^2] Accelerometer noise standard deviation

# Process noise covariance
Q       = Diagonal([1.0e-16, 1.0e-16, 1.0e-16, 1.0e-12, 1.0e-12, 1.0e-12, 1.0e-3])

# GPS Simulation Span
startWeek   = 2170
startDay    = 0
endWeek     = 2174
endDay      = 6

# Create GPS simulation object
gpssim = GPSSim(startWeek, startDay, endWeek, endDay; σρ = σρ, σr = σr)

# Create IMU simulation object
imusim = IMUSim(σa)

# Create EKF
ukf = UKF(xhat0, P0, Q, (σρ^2 + 3*σr^2), σa^2, ts, gpsΔt, gpssim, imusim, scsim;
    lunaPerts = true, α = 1.0, β = 0.0, κ = 3.0+14.0);

@time runFilter!(ukf)

# Get true trajectory for plotting
xtrue = zeros(length(ukf.txp), 3)
for i in 1:length(ukf.txp)
    (y0l, us) = OptimalEstimationProject.GetStateAndControl(scsim, ukf.txp[i])
    xtrue[i, :] .= y0l[1:3]
end

# Plotting
plot(ukf, xtrue, length(ukf.txp))

max_stars_repo_name: GrantHecht/OptimalEstimationProject


text:
"""
    RowProjectionVector{T} <: AbstractArray{T, 1}

Implicit vector projecting a sequence of rows of a given vector with elements of type `T`.

# Constructors
    RowProjectionVector(base, rows)
    RowProjectionVector(base, rows...)
    RowProjectionVector{T}(base, rows)

Create a new row-projected vector from the given `base` vector and a list or sequence of
`rows` corresponding to it. Rows can be used multiple times and in any order.

# Examples

```jldoctest; setup = :(using ImplicitArrays)
julia> v = [10, 20, 30, 40, 50];

julia> RowProjectionVector(v, [1, 3, 5])
3-element RowProjectionVector{Int64}:
 10
 30
 50

julia> RowProjectionVector(v, 2, 3, 2)
3-element RowProjectionVector{Int64}:
 20
 30
 20
```

!!! note
    Changing an element through `RowProjectionVector` actually changes the element in the
    underlying array. The changes are still reflected in the `RowProjectionVector`. Thus,
    caution is advised when changing elements of rows that are contained multiple times
    in the projection, in order to avoid unwanted side effects.
"""
struct RowProjectionVector{T} <: AbstractArray{T, 1}
    "Base vector from which the projection is created"
    base::AbstractArray{T, 1}
    "Rows used for the projection"
    rows::Vector{Int}

    function RowProjectionVector{T}(
        base::AbstractArray{T, 1},
        rows::Vector{Int}
    ) where T
        @boundscheck _checkrowbounds(base, rows)
        new(base, rows)
    end
end

RowProjectionVector(
    base::AbstractArray{T, 1},
    rows::Vector{Int}
) where T = RowProjectionVector{T}(base, rows)

RowProjectionVector(
    base::AbstractArray{T, 1},
    rows::Vararg{Int, N}
) where {T, N} = RowProjectionVector(base, collect(rows))

@inline Base.size(
    M::RowProjectionVector{T}
) where T = (length(M.rows),)

@inline Base.getindex(
    M::RowProjectionVector{T},
    i::Int
) where T = Base.getindex(M.base, M.rows[i])

@inline Base.setindex!(
    M::RowProjectionVector{T},
    v::Any,
    i::Int
) where T = Base.setindex!(M.base, v, M.rows[i])

"""
    RowProjectionMatrix{T} <: AbstractArray{T, 2}

Implicit matrix projecting a sequence of rows of a given matrix with elements of type `T`.

# Constructors
    RowProjectionMatrix(base, rows)
    RowProjectionMatrix(base, rows...)
    RowProjectionMatrix{T}(base, rows)

Create a new row-projected matrix from the given `base` matrix and a list or sequence of
`rows` corresponding to it. Rows can be used multiple times and in any order.

# Examples

```jldoctest; setup = :(using ImplicitArrays)
julia> M = [10 20 30; 40 50 60; 70 80 90];

julia> RowProjectionMatrix(M, [1, 3])
2×3 RowProjectionMatrix{Int64}:
 10  20  30
 70  80  90

julia> RowProjectionMatrix(M, 2, 3, 2)
3×3 RowProjectionMatrix{Int64}:
 40  50  60
 70  80  90
 40  50  60
```

!!! note
    Changing an element through `RowProjectionMatrix` actually changes the element in the
    underlying array. The changes are still reflected in the `RowProjectionMatrix`. Thus,
    caution is advised when changing elements of rows that are contained multiple times
    in the projection, in order to avoid unwanted side effects.
"""
struct RowProjectionMatrix{T} <: AbstractArray{T, 2}
    "Base matrix from which the projection is created"
    base::AbstractArray{T, 2}
    "Rows used for the projection"
    rows::Vector{Int}

    function RowProjectionMatrix{T}(
        base::AbstractArray{T, 2},
        rows::Vector{Int}
    ) where T
        @boundscheck _checkrowbounds(base, rows)
        new(base, rows)
    end
end

RowProjectionMatrix(
    base::AbstractArray{T, 2},
    rows::Vector{Int}
) where T = RowProjectionMatrix{T}(base, rows)

RowProjectionMatrix(
    base::AbstractArray{T, 2},
    rows::Vararg{Int, N}
) where {T, N} = RowProjectionMatrix(base, collect(rows))

@inline Base.size(
    M::RowProjectionMatrix{T}
) where T = (length(M.rows), size(M.base)[2])

@inline Base.getindex(
    M::RowProjectionMatrix{T},
    I::Vararg{Int, 2}
) where T = Base.getindex(M.base, M.rows[I[1]], I[2])

@inline Base.setindex!(
    M::RowProjectionMatrix{T},
    v::Any,
    I::Vararg{Int, 2}
) where T = Base.setindex!(M.base, v, M.rows[I[1]], I[2])

function _checkrowbounds(
    base::Union{AbstractArray{T, 1}, AbstractArray{T, 2}},
    rows::Vector{Int}
) where T
    for row ∈ rows
        checkbounds(base, row, :)
    end
end


text:
using DataStructures

struct Graph
    vertices::Vector{Int}
    edges::Set{Tuple{Tuple, Int}}
end

function getNeighbors(arr, (i, j))
    maxrow = size(arr, 1)
    maxcol = size(arr, 2)
    neighboridxs = Tuple[]
    if i == 1 || i == maxrow || j == 1 || j == maxcol
        if (i,j) == (1,1)
            neighboridxs = push!(neighboridxs, (1,2))
            neighboridxs = push!(neighboridxs, (2,1))
        elseif (i,j) == (1, maxcol)
            neighboridxs = push!(neighboridxs, (1,maxcol-1))
            neighboridxs = push!(neighboridxs, (2,maxcol))
        elseif (i,j) == (maxrow, 1)
            neighboridxs = push!(neighboridxs, (maxrow-1,1))
            neighboridxs = push!(neighboridxs, (maxrow,2))
        elseif (i,j) == (maxrow, maxcol)
            neighboridxs = push!(neighboridxs, (maxrow-1,maxcol))
            neighboridxs = push!(neighboridxs, (maxrow,maxcol-1))
        elseif i == 1
            neighboridxs = push!(neighboridxs, (i, j-1))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
        elseif j == 1
            neighboridxs = push!(neighboridxs, (i-1, j))
            neighboridxs = push!(neighboridxs, (i, j+1))
            neighboridxs = push!(neighboridxs, (i+1, j))
        elseif j == maxcol
            neighboridxs = push!(neighboridxs, (i-1, maxcol))
            neighboridxs = push!(neighboridxs, (i+1, maxcol))
            neighboridxs = push!(neighboridxs, (i, maxcol - 1))
        elseif i == maxrow
            neighboridxs = push!(neighboridxs, (maxrow, j-1))
            neighboridxs = push!(neighboridxs, (maxrow, j+1))
            neighboridxs = push!(neighboridxs, (maxrow-1, j))
        end
    else
        neighboridxs = push!(neighboridxs, (i, j-1))
        neighboridxs = push!(neighboridxs, (i, j+1))
        neighboridxs = push!(neighboridxs, (i+1, j))
        neighboridxs = push!(neighboridxs, (i-1, j))
    end
    return neighboridxs
end

function toGraph(filename, part2 = false)
    values = readlines(filename*".txt")
    values = split.(values, "")
    intvalues = Vector{Int}[]
    for line in values
        convline = parse.(Int, line)
        intvalues = push!(intvalues, convline)
    end
    valarr = zeros(Int, size(intvalues)[1],size(intvalues[1])[1])
    for i in 1:length(intvalues)
        for j in 1:length(intvalues[1])
            valarr[i,j] = intvalues[i][j]
        end
    end
    if part2
        origvalarr = copy(valarr)
        arrs = []
        for i in 1:5
            line = []
            for j in 1:5
                arr = copy(origvalarr)
                toadd_x = i - 1
                toadd_y = j - 1
                toadd = toadd_x + toadd_y
                for k in 1:size(origvalarr,1)
                    for l in 1:size(origvalarr,2)
                        if origvalarr[k,l] + toadd != 9
                            arr[k,l] = (origvalarr[k,l] + toadd) % 9
                            if arr[k,l] == 0
                                arr[k,l] = 1
                            end
                        else
                            arr[k,l] = 9
                        end
                    end
                end
                line = push!(line, arr)
            end
            arrs = push!(arrs, line)
        end
        larrs = []
        for line in arrs
            larr = line[1]
            for i in 1:4
                larr = hcat(larr,line[i+1])
            end
            larrs = push!(larrs, larr)
        end
        valarr = larrs[1]
        for i in 1:4
            valarr = vcat(valarr,larrs[i+1])
        end
    end

    vertices = Int[]
    edges = Set{Tuple{Tuple,Int}}()
    for i in 1:size(valarr)[1]
        for j in 1:size(valarr)[2]
            vidx = ((i-1) * size(valarr)[2]) + j
            vertices = push!(vertices,vidx)
            ngb = getNeighbors(valarr, (i,j))
            for n in ngb
                org = vidx
                nidx = ((n[1]-1) * size(valarr)[2]) + n[2]
                weight = valarr[n[1], n[2]]
                edges = push!(edges, ((org,nidx),weight))
            end
        end
    end
    graph = Graph(vertices, edges)
    return graph, valarr
end

function get2Dcoords(x,valarr)
    row = Int(floor(x/size(valarr, 2))) + 1
    col = Int(x % (size(valarr, 2)))
    if col == 0
        col = size(valarr, 2)
        row = row - 1
    end
    return (row,col)
end

function dijkstra(graph, start, stop, valarr)
    dist = zeros(size(graph.vertices))
    previous = zeros(Int, size(graph.vertices))
    for (i,el) in enumerate(dist)
        dist[i] = Inf
        previous[i] = -1
    end
    dist[start] = 0

    G = PriorityQueue()
    for v in graph.vertices
        G = enqueue!(G, v, dist[v])
    end

    while isempty(G) == false
        u = dequeue!(G)
        if u == stop
            break
        end
        ngb = getNeighbors(valarr, get2Dcoords(u, valarr))
        for n in ngb
            d = dist[u] + valarr[n[1], n[2]]
            n = (n[1]-1) * size(valarr, 2) + n[2]
            if d < dist[n]
                dist[n] = d
                G[n] = d
                previous[n] = u
            end
        end
    end
    return dist, previous
end

function solve(filename, part2 = false)
    graph, valarr = toGraph(filename, part2)
    stopcoord = length(graph.vertices)
    dist, prev = dijkstra(graph, 1, stopcoord, valarr)

    println(dist[stopcoord])
end

solve("example")
solve("example", true)
solve("input")
solve("input", true)
#solve("test")
#solve("test2")


text:
using DataFrames, CSV
using Profile
include("../src/sdistance.jl") # 必ずダブルクオーテーション
include("../test/draw.jl")
include("../src/SignedDistanceFunction.jl")
include("../src/environments.jl")
include("../test/APT.jl")
import .Draw: parformance_graphs
using .SignedDistanceFunction
using .APT

# ====== Application product testing ======
p = plot_for_debug(parse(Int, ARGS[2]), "./test/mock_csv_data/interface.csv")
# p = plot_for_debug(parse(Int, ARGS[2]), "./test/mock_csv_data/circle.csv")
# p = plot_for_debug(parse(Int, ARGS[2]), "./test/mock_csv_data/interface.csv", "multi")
# p = plot_for_debug(parse(Int, ARGS[2]), "./test/mock_csv_data/multiple_curves.csv", "multi")
# p = plot_for_debug(parse(Int, ARGS[2]), "./test/mock_csv_data/double_circle.csv", "multi")
# plots_contours([i for i = 50:50:300], "./test/mock_csv_data/interface.csv", "multi")
# plots_contours([i for i = 50:50:300], "./test/mock_csv_data/interface.csv")
# plots_wireframe([i for i = 50:50:300], "./test/mock_csv_data/interface.csv", "multi")

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
module GradientStarvation

# activate environment
using Pkg
if isfile("Project.toml") && isfile("Manifest.toml")
    Pkg.activate(".")
end

using DataStructures

include("data.jl")
include("experiments.jl")
include("figures.jl")
include("neuralnetwork.jl")
include("plotlib.jl")
include("plotutils.jl")
include("regularization.jl")
include("structs.jl")

# includes & export needed for documentation
export Data, Figures, Experiments, NeuralNetwork, PlotLib, PlotUtils, Regularization

# figure 1
Figures.create_figure_1("medium")

# figures 2-13
# run different experiments on the 2D classification task (moons)
experiments = OrderedDict(
    "GD Δ0.5" => Experiment("GD", 0.5, false),
    "GD Δ1.0" => Experiment("GD", 1.0, false),
    "GD Δ0.5 + SD" => Experiment("GD", 0.5, true),
    "GD Δ1.0 + SD" => Experiment("GD", 1.0, true),

    "SGD Δ0.5" => Experiment("SGD", 0.5, false),
    "SGD Δ1.0" => Experiment("SGD", 1.0, false),
    "SGD Δ0.5 + SD" => Experiment("SGD", 0.5, true),
    "SGD Δ1.0 + SD" => Experiment("SGD", 1.0, true),

    "WD Δ0.5" => Experiment("WD", 0.5, false),
    "WD Δ1.0" => Experiment("WD", 1.0, false),
    "WD Δ0.5 + SD" => Experiment("WD", 0.5, true),
    "WD Δ1.0 + SD" => Experiment("WD", 1.0, true),

    "ADAM Δ1.0" => Experiment("ADAM", 1.0, false),
    "ADAM Δ0.5" => Experiment("ADAM", 0.5, false),
    "ADAM Δ0.5 + SD" => Experiment("ADAM", 0.5, true),
    "ADAM Δ1.0 + SD" => Experiment("ADAM", 1.0, true),
)
Experiments.run_experiments(experiments, show_plots = false, random_moons = false)

# additional experiments with different parameters

# figure 14
Experiments.run_experiments(
    OrderedDict("GD Δ0.5 + SD 10000" => Experiment("GD", 0.5, true)),
    show_plots = false, random_moons = false, epochs = 10000
)
# figure 15, 16
add_exp_adam = OrderedDict("ADAM Δ0.5 + SD" => Experiment("ADAM", 0.5, true))
Experiments.run_experiments(add_exp_adam, show_plots = false, random_moons = false, learning_rate = 1e-4, prefix = "1e-4")
Experiments.run_experiments(add_exp_adam, show_plots = false, random_moons = false, learning_rate = 1e-3, prefix = "1e-3")

end


text:
using Documenter, FinEconResearchUtils

makedocs(;
    modules=[FinEconResearchUtils],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tbeason/FinEconResearchUtils.jl/blob/{commit}{path}#L{line}",
    sitename="FinEconResearchUtils.jl",
    authors="Tyler Beason",
    assets=[],
)

deploydocs(;
    repo="github.com/tbeason/FinEconResearchUtils.jl",
)


text:
# define constants whose name changed with different versions
mf = Pkg.Operations.Context().env.manifest
ver = mf[findfirst(v -> v.name == "GLPK", mf)].version
if ver >= v"1"
    const GLPK_ON = GLPK.GLP_ON
else
    const GLPK_ON = GLPK.ON
end


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct MonitoredLine <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    r::Float64  # System per-unit value
    x::Float64  # System per-unit value
    b::NamedTuple{(:from, :to), Tuple{Float64, Float64}}  # System per-unit value
    flowlimits::NamedTuple{(:from_to, :to_from), Tuple{Float64, Float64}}  # TODO: throw warning above max SIL
    rate::Float64  # TODO: compare to SIL (warn) (theoretical limit)
    anglelimits::Min_Max
    internal::InfrastructureSystemsInternal
end

function MonitoredLine(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, flowlimits, rate, anglelimits, )
    MonitoredLine(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, flowlimits, rate, anglelimits, InfrastructureSystemsInternal())
end

function MonitoredLine(; name, available, activepower_flow, reactivepower_flow, arc, r, x, b, flowlimits, rate, anglelimits, )
    MonitoredLine(name, available, activepower_flow, reactivepower_flow, arc, r, x, b, flowlimits, rate, anglelimits, )
end

# Constructor for demo purposes; non-functional.

function MonitoredLine(::Nothing)
    MonitoredLine(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        b=(from=0.0, to=0.0),
        flowlimits=(from_to=0.0, to_from=0.0),
        rate=0.0,
        anglelimits=(min=-1.571, max=1.571),
    )
end

"""Get MonitoredLine name."""
get_name(value::MonitoredLine) = value.name
"""Get MonitoredLine available."""
get_available(value::MonitoredLine) = value.available
"""Get MonitoredLine activepower_flow."""
get_activepower_flow(value::MonitoredLine) = value.activepower_flow
"""Get MonitoredLine reactivepower_flow."""
get_reactivepower_flow(value::MonitoredLine) = value.reactivepower_flow
"""Get MonitoredLine arc."""
get_arc(value::MonitoredLine) = value.arc
"""Get MonitoredLine r."""
get_r(value::MonitoredLine) = value.r
"""Get MonitoredLine x."""
get_x(value::MonitoredLine) = value.x
"""Get MonitoredLine b."""
get_b(value::MonitoredLine) = value.b
"""Get MonitoredLine flowlimits."""
get_flowlimits(value::MonitoredLine) = value.flowlimits
"""Get MonitoredLine rate."""
get_rate(value::MonitoredLine) = value.rate
"""Get MonitoredLine anglelimits."""
get_anglelimits(value::MonitoredLine) = value.anglelimits
"""Get MonitoredLine internal."""
get_internal(value::MonitoredLine) = value.internal


text:
using Test
using DeconvolveDistribution
using Distributions
import Random

@testset "Fhat" begin
    n = 100
    F_X = MixtureModel([Normal(-1, 1), Normal(2, 0.8)], [0.6, 0.4])
    σ_distr = Gamma(1.0, 4.0)
    Random.seed!(1)
    X = rand(F_X, n)
    σ = rand(σ_distr, n)
    U = Normal.(0.0, σ)
    ϵ = rand.(U)
    W = X .+ ϵ

    n_xx = 120
    F_xx = collect(range(minimum(W), stop=maximum(W), length=n_xx))
    num_t = 50
    h = 0.5
    Fhat_xx = Fhat(F_xx, W, num_t, h, U)

    @test length(Fhat_xx) == n_xx
    # CDF starts at 0
    @test -0.01 < Fhat_xx[1] < 0.01
    # CDF ends at 1
    @test -0.99 < Fhat_xx[end] < 1.01
    # CDF is monotonic
    @test all(diff(Fhat_xx) .>= -0.005)
end # testset

@testset "decon" begin
    # simulate some data
    n = 400
    F_X = MixtureModel([Normal(-1, 1), Normal(2, 0.8)], [0.6, 0.4])
    σ_distr = Gamma(1.0, 4.0)
    Random.seed!(1)
    X = rand(F_X, n)
    σ = rand(σ_distr, n)
    U = Normal.(0.0, σ)
    ϵ = rand.(U)
    W = X .+ ϵ

    # estimate F_X from simulated data
    n_xx = 500
    F_xx = collect(range(-8.0, stop=8.0, length=n_xx))
    num_t = 50
    h = 0.3

    decon_settings = FourierDeconv(h, F_xx, num_t)
    decon_distr = decon(decon_settings, W, U)
end

@testset "fix_CDF" begin
    n = 100
    F_X = MixtureModel([Normal(-1, 1), Normal(2, 0.8)], [0.6, 0.4])
    σ_distr = Gamma(1.0, 4.0)
    Random.seed!(1)
    X = rand(F_X, n)
    σ = rand(σ_distr, n)
    U = Normal.(0.0, σ)
    ϵ = rand.(U)
    W = X .+ ϵ

    n_xx = 120
    F_xx = collect(range(minimum(W), stop=maximum(W), length=n_xx))
    num_t = 20
    h = 0.5
    Fhat_xx = Fhat(F_xx, W, num_t, h, U)

    fix_CDF!(Fhat_xx)

    @test length(Fhat_xx) == n_xx
    # CDF starts at 0
    @test Fhat_xx[1] == 0.0
    # CDF ends at 1
    @test Fhat_xx[end] == 1.0
    # CDF doesn't go outside of [0, 1]
    for Fx in Fhat_xx
        @test 0.0 <= Fx <= 1.0
    end
    # CDF is monotonic
    for dFx in diff(Fhat_xx)
        @test dFx >= 0.0
    end
end


text:
export get_data!, get_data

@doc """
    S = get_data(method, chans [, keywords])
    get_data!(S, method, chans [, keywords])

Wrapper to web requests for time-series data. Request data using `method` from channels `chans` using keywords `KWs`, storing the output in `S`.

* Methods: FDSN, IRIS, PH5
* Channels: See `?web_chanspec`
* Keywords: autoname, demean, detrend, fmt, msr, nd, opts, rad, reg, rr, s, si, src, t, taper, to, ungap, unscale, v, w, xf, y

This function is fully described in the official documentation at https://seisio.readthedocs.io/ under subheading **Web Requests**.

See also: `web_chanspec`, `parsetimewin`, `seis_www`, `SeisIO.KW`
""" get_data
function get_data(method_in::String, C::ChanOpts="*";
   autoname::Bool              = false         , # Auto-generate file names?
     demean::Bool              = false         , # Demean after download?
    detrend::Bool              = false         , # Detrend after download?
        fmt::String            = KW.fmt        , # File format
        msr::Bool              = false         , # Get multi-stage response?
         nd::Real              = KW.nd         , # Number of days per request
       opts::String            = KW.opts       , # Options string
        rad::Array{Float64, 1} = KW.rad        , # Query radius
        reg::Array{Float64, 1} = KW.reg        , # Query region
      prune::Bool              = KW.prune      , # Prune empty channels?
         rr::Bool              = false         , # Remove instrument response?
          s::TimeSpec          = 0             , # Start
         si::Bool              = KW.si         , # Fill station info?
        src::String            = KW.src        , # Data source
      taper::Bool              = false         , # Taper after download?
          t::TimeSpec          = (-600)        , # End or Length (s)
         to::Int64             = KW.to         , # Timeout (s)
      ungap::Bool              = false         , # Remove time gaps?
    unscale::Bool              = false         , # Unscale (divide by gain)?
          v::Integer           = KW.v          , # Verbosity
          w::Bool              = KW.w          , # Write to disc?
         xf::String            = "FDSNsta.xml" , # XML save file
          y::Bool              = KW.y          , # Sync
          )

  # Parse time window
  α, ω = parsetimewin(s, t)

  # KWs that ovewrite other KWs (rare; keep this behavior to an absolute minimum!)
  (autoname == true) && (w = true)

  # Generate SeisData
  S = SeisData()

  # Condense requests as much as possible
  if method_in == "FDSN"
    if isa(C, String)
      C = parse_chstr(C, ',', true, false)
    elseif isa(C, Array{String,1})
      C = parse_charr(C, '.', true)
    end
    R = minreq(C)
    parse_err = FDSNget!(S, R, α, ω,
      autoname, fmt, msr, nd, opts, rad, reg, si, src, to, v, w, xf, y)
  elseif method_in == "PH5"
    if isa(C, String)
      C = parse_chstr(C, ',', true, false)
    elseif isa(C, Array{String,1})
      C = parse_charr(C, '.', true)
    end
    R = minreq(C)
    parse_err = FDSNget!(S, R, α, ω,
      autoname, fmt, msr, nd, opts, rad, reg, si, src*"PH5", to, v, w, xf, y)
  elseif method_in == "IRIS"
    if isa(C, String)
      R = String[strip(String(j)) for j in split(C, ',')]
    elseif isa(C, Array{String,2})
      NC = size(C,1)
      R = Array{String,1}(undef, NC)
      for j = 1:NC
        R[j] = join(C[j,:],'.')
      end
    else
      R = deepcopy(C)
    end
    parse_err = IRISget!(S, R, α, ω, fmt, opts, to, v, w)
  end

  # DND DND DND
  # Wrapper to a generic "get" function -- leave as example code
  # getfield(SeisIO, Symbol(string(method_in, "get!")))(S, C, α, ω,
  #          f = f, opts = opts, si = si, src = src, to = to, v = v, w = w)
  # DND DND DND

  # ========================================================================
  # Viable operations for any data type
  if prune == true
    if parse_err == false
      v > 0 && @info(tnote("Removing empty channels."))
      prune!(S)
    else
      v > 0 && @info(tnote("Can't prune empty channels; web request wasn't fully parsed."))
    end
  end

  # Sync
  if y == true
    v > 0 && @info(tnote("Synchronizing data."))
    sync!(S)
  end

  # ========================================================================
  # Operations that only make sense for seismic or seismoacoustic data

  if any([demean, detrend, rr, taper, unscale])
    # Get list of channels with sane instrument codes
    CC = get_seis_channels(S)

    # Unscale
    if unscale == true
      v > 0 && @info(tnote("dividing out gain of seismic channels"))
      unscale!(S, chans=CC)
    end

    # Demean
    if demean == true
      v > 0 && @info(tnote("removing mean of seismic channels"))
      demean!(S, chans=CC)
    end

    # Detrend
    if detrend == true
      v > 0 && @info(tnote("detrending seismic channels"))
      detrend!(S, chans=CC)
    end

    # Taper
    if taper == true
      v > 0 && @info(tnote("tapering seismic channel data"))
      taper!(S, chans=CC)
    end

    # Ungap
    if ungap == true
      v > 0 && @info(tnote("ungapping seismic channel data"))
      ungap!(S, chans=CC)
    end

    # Remove response
    if rr == true
      v > 0 && @info(tnote("removing instrument response from seismic channels"))
      remove_resp!(S, chans=CC)
    end
  end

  return S
end

@doc (@doc get_data)
function get_data!(S::SeisData, method_in::String, C::ChanOpts="*";
   autoname::Bool              = false         , # Auto-generate file names?
     demean::Bool              = false         , # Demean after download?
    detrend::Bool              = false         , # Detrend after download?
        fmt::String            = KW.fmt        , # File format
        msr::Bool              = false         , # Get multi-stage response?
         nd::Real              = KW.nd         , # Number of days per request
       opts::String            = KW.opts       , # Options string
        rad::Array{Float64, 1} = KW.rad        , # Query radius
        reg::Array{Float64, 1} = KW.reg        , # Query region
      prune::Bool              = KW.prune      , # Prune empty channels?
         rr::Bool              = false         , # Remove instrument response?
          s::TimeSpec          = 0             , # Start
         si::Bool              = KW.si         , # Fill station info?
        src::String            = KW.src        , # Data source
      taper::Bool              = false         , # Taper after download?
          t::TimeSpec          = (-600)        , # End or Length (s)
         to::Int64             = KW.to         , # Timeout (s)
      ungap::Bool              = false         , # Remove time gaps?
    unscale::Bool              = false         , # Unscale (divide by gain)?
          v::Integer           = KW.v          , # Verbosity
          w::Bool              = KW.w          , # Write to disc?
         xf::String            = "FDSNsta.xml" , # XML save file
          y::Bool              = KW.y          , # Sync
          )

  U = get_data(method_in, C,
               autoname=autoname,
                 demean=demean,
                detrend=detrend,
                    fmt=fmt,
                    msr=msr,
                     nd=nd,
                   opts=opts,
                  prune=prune,
                    rad=rad,
                    reg=reg,
                     rr=rr,
                      s=s,
                     si=si,
                    src=src,
                      t=t,
                  taper=taper,
                     to=to,
                  ungap=ungap,
                unscale=unscale,
                      v=v,
                      w=w,
                     xf=xf,
                      y=y)
  v > 2 && println(stdout, "S = \n", U)
  append!(S, U)
  return nothing
end


text:
Base.@pure @generated function fndims(f)
    if hasmethod(ndims, Tuple{f})
        N = ndims(f)
        return :($N)
    else
        # Function is not type-stable if fndims_via_hasmethod
        # is copied here
        return quote
            (d = fndims_via_hasmethod(f, Union{})) > 0 && return d
            (d = fndims_via_hasmethod(f, Float64)) > 0 && return d
            (d = fndims_via_hasmethod(f, ComplexF64)) > 0 && return d
            (d = fndims_via_hasmethod(f, Float32)) > 0 && return d
            (d = fndims_via_hasmethod(f, ComplexF32)) > 0 && return d
            (d = fndims_via_hasmethod(f, Int)) > 0 && return d
            throw(ArgumentError("Could not determine the dimension of function $f"))
        end
    end
end
Base.@pure function fndims_via_hasmethod(f,T)
    Base.Cartesian.@nexprs 10 k->begin
        hasmethod(f, NTuple{k,T}) && return k
    end
    return 0
end

"""
    grideval(f, x::AbstractVector...)
    grideval(f, x::NTuple{N,AbstractVector})

Evaluate `f` on the grid spanned by `x`.

# Examples
```
julia> grideval(*, ([1,2],[3,4]))
2×2 Array{Int64,2}:
 3  4
 6  8
```
"""
grideval(f, x) = grideval(f,ntuple(i->x,Val(fndims(f))))
grideval(f, x...) = grideval(f,x)
grideval(f, x::NTuple{N,Any}) where {N} = f.(reshape4grideval(x)...)

reshape4grideval(x::NTuple{N,Any}) where {N} = reshape4grideval_.(x, Val(N), ntuple(identity, Val(N)))
reshape4grideval_(xi::Number,_,_) = xi
reshape4grideval_(xi::AbstractVector,::Val{N},i) where {N} = reshape(xi,ntuple(j -> i==j ? length(xi) : 1, Val(N)))

using MacroTools

struct GridFunction{N,F}
    fun::F
end
Base.ndims(::GridFunction{N}) where {N} = N
Base.ndims(::Type{<:GridFunction{N}}) where {N} = N
GridFunction{N}(f) where {N} = GridFunction{N,typeof(f)}(f)

(gf::GridFunction{N})(x::Vararg{Any,N}) where {N} = gf.fun(x...)
grideval(gf::GridFunction{N}, x::NTuple{N,Any}) where {N} = gf.fun(reshape4grideval(x)...)

cvec(x::Number) = x
cvec(x::AbstractArray) = vec(x)

"""
    @gridfun(f)

Assemble an anonymous function where nested function calls
annotated with `\$` get `grideval`ed separately.

# Examples

`grideval(@gridfun((x1,x2)->f(x1,x2) - \$g(x1,x2)),x)` is equivalent to

    grideval((x1,x2)->f(x1,x2) - grideval(g,x),x)
"""
macro gridfun(f)
    @capture(f, (x_) -> body_) || error("@gridfun must be used on an anonymous function")
    body = MacroTools.postwalk(body) do expr
        (@capture(expr, g_(y__)) && isa(g,Expr) && g.head == :$) || return expr
        return :($grideval($g, $cvec.(($(y...),))))
    end
    if isa(x, Expr) N = length(x.args)
    else N = 1; end
    return esc(:($GridFunction{$N}($x -> @. $body)))
end


text:
using NCDatasets, ColorSchemes, RPRMakie
using ImageShow

# Taken from https://lazarusa.github.io/BeautifulMakie/GeoPlots/topography/
cmap = dataset = Dataset(joinpath(@__DIR__, "ETOPO1_halfdegree.nc"))
lon = dataset["lon"][:]
lat = dataset["lat"][:]

data = Float32.(dataset["ETOPO1avg"][:, :])
function toCartesian(lon, lat; r=1, cxyz=(0, 0, 0))
    lat, lon = lat * π / 180, lon * π / 180
    x = cxyz[1] + (r + 80_000) * cos(lat) * cos(lon)
    y = cxyz[2] + (r + 80_000) * cos(lat) * sin(lon)
    z = cxyz[3] + (r + 80_000) * sin(lat)
    return (x, y, z) ./ 80_000
end

function lonlat3D(lon, lat, data; cxyz=(0, 0, 0))
    xyzw = zeros(size(data)..., 3)
    for (i, lon) in enumerate(lon), (j, lat) in enumerate(lat)
        x, y, z = toCartesian(lon, lat; r=data[i, j], cxyz=cxyz)
        xyzw[i, j, 1] = x
        xyzw[i, j, 2] = y
        xyzw[i, j, 3] = z
    end
    return xyzw[:, :, 1], xyzw[:, :, 2], xyzw[:, :, 3]
end
# this is needed in order to have a closed surface
lonext = cat(collect(lon), lon[1]; dims=1)
dataext = begin
    tmpdata = zeros(size(lon)[1] + 1, size(lat)[1])
    tmpdata[1:size(lon)[1], :] = data
    tmpdata[size(lon)[1] + 1, :] = data[1, :]
    tmpdata
end
#now let's plot that surface!
xetopo, yetopo, zetopo = lonlat3D(lonext, lat, dataext)

begin
    r = 30
    lights = [PointLight(Vec3f(2, 1, 3), RGBf(r, r, r))]
    fig = Figure(; resolution=(1200, 1200), backgroundcolor=:black)
    ax = LScene(fig[1, 1]; show_axis=false)#, scenekw=(lights=lights,))
    pltobj = surface!(ax, xetopo, yetopo, zetopo; color=dataext, colormap=:hot, colorrange=(-6000, 5000))
    cam = cameracontrols(ax.scene)
    cam.fov[] = 10
    cam.eyeposition[] = Vec3f(3, 1, 1)
    cam.lookat[] = Vec3f(0)
    cam.upvector[] = Vec3f(0, 0, 1)
    Makie.update_cam!(ax.scene, cam)
    RPRMakie.activate!(; iterations=100)
    display(ax.scene)
end