text:
# inertial pose3

export
  PreintegralCompensationGradients,
  InertialPose3Container,
  oplus,
  ⊕,
  getSample,
  InertialPose3,
  PackedInertialPose3,
  PriorInertialPose3,
  PackedPriorInertialPose3,
  compare

abstract type PreintContainer end

"""
$(TYPEDEF)
"""
mutable struct PreintegralCompensationGradients <: PreintContainer
  # First order terms
  dPdDa::Array{Float64,2}
  dVdDa::Array{Float64,2}
  dPdDw::Array{Float64,2}
  dVdDw::Array{Float64,2}
  dRdDw::Array{Float64,2}
  # second order
  # ...
  PreintegralCompensationGradients(;
        dPdDa::Array{Float64,2}=zeros(3,3),
        dVdDa::Array{Float64,2}=zeros(3,3),
        dPdDw::Array{Float64,2}=zeros(3,3),
        dVdDw::Array{Float64,2}=zeros(3,3),
        dRdDw::Array{Float64,2}=zeros(3,3)  ) = new(dPdDa,dVdDa,dPdDw,dVdDw,dRdDw)
end
"""
$(TYPEDEF)
"""
mutable struct InertialPose3Container <: PreintContainer
  rRp::Array{Float64,2}
  rPosp::Vector{Float64}
  rVelp::Vector{Float64}
  pBw::Vector{Float64}
  pBa::Vector{Float64}
  rnTime::Int
  InertialPose3Container(;
        rRp::Array{Float64,2}=Matrix{Float64}(LinearAlgebra.I, 3,3),
        rPosp::Vector{Float64}=zeros(3),
        rVelp::Vector{Float64}=zeros(3),
        pBw::Vector{Float64}=zeros(3),
        pBa::Vector{Float64}=zeros(3),
        rnTime::Int=0  ) = new( rRp,rPosp,rVelp,pBw,pBa,rnTime )
end

function veeQuaternion(ip3::InertialPose3Container)
  ret = zeros(16)
  ret[1:7] = veeQuaternion(  SE3(ip3.rPosp, SO3(rRp))  )

end

function oplus(xi::InertialPose3Container, Dx::InertialPose3Container)
  return InertialPose3Container(
        rRp = xi.rRp*Dx.rRp,
        rPosp = xi.rPosp + xi.rRp*Dx.rPosp,
        rVelp = xi.rVelp + xi.rRp*Dx.rVelp*(Dx.rnTime*1e-9),
        pBw = xi.pBw + Dx.pBw,
        pBa = xi.pBa + Dx.pBa,
        rnTime = xi.rnTime+Dx.rnTime  )
end
⊕(xi::InertialPose3Container, Dx::InertialPose3Container) = oplus(xi,Dx)

function zetaEmbedding(posei::InertialPose3Container, posej::InertialPose3Container; rGrav=[0.0;0.0;9.81])
  z = zeros(30)
  z[1:3] = logmap(SO3(posei.rRp'*posej.rRp))
  z[4:6] = posej.pBw
  z[7:9] = posej.rVelp
  z[10:12] = posej.rPosp
  z[13:15] = posej.pBa
  z[16:18] = posei.pBw
  z[19:21] = posei.rVelp
  z[22:24] = posei.rPosp
  z[25:27] = posei.pBa
  z[28:30] = rGrav
  return z
end

# predict the preintegral value
function constructL(posei::InertialPose3Container, Dt::Float64)
  biRw = posei.rRp'
  L = zeros(15,30)
  L[1:3,1:3] = Matrix{Float64}(LinearAlgebra.I, 3,3)
  L[7:9,7:9] = biRw
  L[10:12,10:12] = biRw
  L[7:9,19:21] = -biRw
  L[10:12,19:21] = -biRw*Dt # preintegral constant velocity component in delta position
  L[10:12,22:24] = -biRw
  return L
end

# First term in Taylor expansion
function constructC1(posei::InertialPose3Container, pido::T, Dt::Float64) where {T <: PreintContainer}
  biRw = posei.rRp'
  C1 = zeros(15,30)
  g1 = -biRw*Dt
  g2 = 0.5*g1*Dt
  C1[4:6,4:6] = Matrix{Float64}(LinearAlgebra.I, 3,3)
  C1[13:15,13:15] = Matrix{Float64}(LinearAlgebra.I, 3,3)
  C1[4:6,16:18] = -Matrix{Float64}(LinearAlgebra.I, 3,3)
  C1[13:15,25:27] = -Matrix{Float64}(LinearAlgebra.I, 3,3)
  C1[7:9,28:30] = g1
  C1[10:12,28:30] = g2
  C1[1:3,16:18] = pido.dRdDw
  C1[7:9,16:18] = pido.dVdDw
  C1[10:12,16:18] = pido.dPdDw
  C1[7:9,25:27] = pido.dVdDa
  C1[10:12,25:27] = pido.dPdDa
  return C1
end

function preintMeas(pido::T) where {T <: PreintContainer}
  return [logmap(SO3(pido.rRp)); pido.pBw; pido.rVelp; pido.rPosp; pido.pBa]
  # return [logmap(SO3(pido.iRj));zeros(3);pido.iDvj;pido.iDppj;zeros(3)] # temporarily suppressing bias updates
end

function predictDeltaXij(pido::T, posei::InertialPose3Container, posej::InertialPose3Container; rGrav=[0.0;0.0;9.81]) where {T <: PreintContainer}
  zet = zetaEmbedding(posei, posej, rGrav=rGrav)
  Dt = Float64(posej.rnTime - posei.rnTime)*1e-9
  L = constructL(posei,Dt)
  C1 = constructC1(posei, pido, Dt)
  return (L-C1)*zet #-0.5*C2*zet.^2
end

function residual!(res::Vector{Float64},
        pioc::InertialPose3Container,
        picg::PreintegralCompensationGradients,
        posei::InertialPose3Container,
        posej::InertialPose3Container, rGrav=[0.0;0.0;9.81]  )
  #
  res[1:15] = preintMeas(pioc) - predictDeltaXij(picg, posei, posej, rGrav=rGrav)
  nothing
end

"""
$(TYPEDEF)

Inertial Odometry version of preintegration procedure and used as a factor between InertialPose3 types for inertial navigation in factor graphs.
"""
mutable struct InertialPose3 <: FunctorPairwise #RoME.BetweenPoses
  # Zij is entropy of veeLie15, pioc is preintegral measurements, pido is compensation gradients.
  Zij::Distribution
  pioc::InertialPose3Container
  picg::PreintegralCompensationGradients
  reuse::Vector{Tuple{InertialPose3Container,InertialPose3Container, InertialPose3Container}} # number of threads

  InertialPose3(
    zij::Distributions.MvNormal,
    pioc::InertialPose3Container,
    picg::PreintegralCompensationGradients  ) =
      new( zij, pioc, picg,
        fill( (InertialPose3Container(),
               InertialPose3Container(),
               InertialPose3Container()),
        Threads.nthreads() )
      )
end

function getSample(ip3::InertialPose3, N::Int=1)
  return (rand( ip3.Zij, N ), )
end

function (ip3::InertialPose3)(
            res::Vector{Float64},
            userdata ,
            idx::Int,
            meas::Tuple,
            wIPi::Array{Float64,2},
            wIPj::Array{Float64,2}  )
  #
  # Function can be massively improved. Just getting it all wired at first.
  # get pointer to memory, local to this thread
  posei, posej, ENT = ip3.reuse[Threads.threadid()]

  ## this part must be moved in the general ApproxConv.jl code, since this reassignment
  ## does not have to happen every time
  # repoint to existing values for first pose
  posei.rRp[1:3,1:3] = convert(SO3, Euler(wIPi[4:6, idx]...)).R
  posei.rPosp = wIPi[1:3, idx]
  posei.rVelp = wIPi[7:9, idx]
  posei.pBw = wIPi[10:12, idx]
  posei.pBa = wIPi[13:15, idx]

  # repoint to existing load values for second pose
  posei.rRp[1:3,1:3] = convert(SO3, Euler(wIPj[4:6, idx]...)).R
  posei.rPosp = wIPj[1:3, idx]
  posei.rVelp = wIPj[7:9, idx]
  posei.pBw = wIPj[10:12, idx]
  posei.pBa = wIPj[13:15, idx]

  # repoint to existing load values for second pose
  ENT.rRp[1:3,1:3] = convert(SO3, so3(meas[1][4:6, idx])).R
  ENT.rPosp = meas[1][1:3, idx]
  ENT.rVelp = meas[1][7:9, idx]
  ENT.pBw = meas[1][10:12, idx]
  ENT.pBa = meas[1][13:15, idx]
  ENT.rnTime = ip3.pioc.rnTime # need time to enact velocity noise

  keeprntime = deepcopy(posej.rnTime)
  noisyposej = posej⊕ENT
  noisyposej.rnTime = keeprntime  # rest timestamp to correct value

  # function residual!(res::Vector{Float64},
  #         pioc::InertialPose3Container,
  #         picg::PreintegralCompensationGradients,
  #         posei::InertialPose3Container,
  #         posej::InertialPose3Container, rGrav=[0.0;0.0;9.81]  )
  residual!(res, ip3.pioc, ip3.picg, posei,  noisyposej)
  # (res'*(ip3.Zij.Σ.mat\res))[1]
  nothing
end

"""
$(TYPEDEF)
"""
mutable struct PackedInertialPose3 <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation, 3 velocities
  vecCov::Array{Float64,1}
  dimc::Int
  vecpioc::Vector{Float64}
  rnTime::Int
  picgvecdPdDa::Vector{Float64}
  picgvecdVdDa::Vector{Float64}
  picgvecdPdDw::Vector{Float64}
  picgvecdVdDw::Vector{Float64}
  picgvecdRdDw::Vector{Float64}
  PackedInertialPose3() = new()
  PackedInertialPose3(ip3::InertialPose3) = new( ip3.Zij.μ, ip3.Zij.Σ.mat[:], size(ip3.Zij.Σ.mat,1),
            veeQuaternion(ip3.pioc), ip3.pioc.rnTime,
            ip3.dPdDa[:] ,ip3.dVdDa[:] ,ip3.dPdDw[:] ,ip3.dVdDw[:] ,ip3.dRdDw[:] )
end

convert(::Type{PackedInertialPose3}, ip3::InertialPose3) = PackedInertialPose3(ip3)

function convert(::Type{InertialPose3}, pip3::PackedInertialPose3)
  pioc = InertialPose3Container(
      rPosp=pip3.vecpioc[1:3],
      rRp=Quaternion(vecpioc[4],vecpioc[5:7]),
      rVelp=vecpioc[8:10],
      pBw=vecpioc[11:13],
      pBa=vecpioc[14:16],
      rnTime=pip3.rnTime
  )
  picg = PreintegralCompensationGradients(
      reshapeVec2Mat(pip3.picgvecdPdDa, 3),
      reshapeVec2Mat(pip3.picgvecdVdDa, 3),
      reshapeVec2Mat(pip3.picgvecdPdDw, 3),
      reshapeVec2Mat(pip3.picgvecdVdDw, 3),
      reshapeVec2Mat(pip3.picgvecdRdDw, 3)
  )
  InertialPose3(Distributions.MvNormal(pip3.vecZij, reshapeVec2Mat(pip3.vecCov, pip3.dimc)), pioc, picg)
end

function compare(a::PreintegralCompensationGradients, b::PreintegralCompensationGradients)
  TP = true
  TP = TP && norm( a.dPdDa - b.dPdDa ) < 1e-10
  TP = TP && norm( a.dVdDa - b.dVdDa ) < 1e-10
  TP = TP && norm( a.dPdDw - b.dPdDw ) < 1e-10
  TP = TP && norm( a.dVdDw - b.dVdDw ) < 1e-10
  TP = TP && norm( a.dRdDw - b.dRdDw ) < 1e-10
  return TP
end
function compare(a::InertialPose3, b::InertialPose3)
  TP = true
  TP = TP && norm(a.Zij.μ - b.Zij.μ) < 1e-10
  TP = TP && norm(a.Zij.Σ.mat - b.Zij.Σ.mat) < 1e-10
  TP = TP && compare(a.pido, b.pido)
  return TP
end

# we also need a prior for InertialPose3

"""
$(TYPEDEF)
"""
mutable struct PriorInertialPose3 <: IncrementalInference.FunctorSingleton
  Zi::Distribution
end
function getSample(prip3::PriorInertialPose3, N::Int=1)
  return (rand( prip3.Zi, N )', )
end

"""
$(TYPEDEF)
"""
mutable struct PackedPriorInertialPose3 <: IncrementalInference.PackedInferenceType
  vecZi::Array{Float64,1} # 3translations, 3rotation, 3 velocities
  vecCov::Array{Float64,1}
  dimc::Int
  PackedPriorInertialPose3() = new()
  PackedPriorInertialPose3(prip3::PriorInertialPose3) = new( prip3.Zi.μ, prip3.Zi.Σ.mat[:], size(prip3.Zi.Σ.mat,1)  )
end

convert(::Type{PackedPriorInertialPose3}, prip3::PriorInertialPose3) = PackedPriorInertialPose3(prip3)

function convert(::Type{PriorInertialPose3}, pprip3::PackedPriorInertialPose3)
  PriorInertialPose3(Distributions.MvNormal(pprip3.vecZi, reshapeVec2Mat(pprip3.vecCov, pprip3.dimc)) )
end

#


text:
using Pkg
using Coverage

Pkg.test(; coverage=true)
coverage = process_folder()
LCOV.writefile("lcov.info", coverage)

ext: jl


text:
using Test
using PipeLayout

include("helpers.jl")

@testset "basics" begin
    include("gasphysics.jl")
    include("util.jl")
    include("pwl.jl")
    include("diameter.jl")
end

@testset "topology" begin
    include("topology/geosteiner.jl")
    include("topology/grid.jl")
    include("topology/mst.jl")
    include("topology/util.jl")
    include("topology/isomorph.jl")
end

@testset "flow" begin
    include("flow.jl")
end

@testset "pipe dimensioning" begin
    include("models/pipedimensioning.jl")
end

@testset "ground structure & discrete diameter" begin
    include("models/gndstruct_discdiam.jl")
end

@testset "junction location" begin
    include("models/junctionlocation.jl")
    include("models/junctionlocation_nlp.jl")
end


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

module Base64URL
import Base: read, write, close, eof, empty!
#export Base64EncodePipe, Base64DecodePipe, base64encode, base64decode

# Base64EncodePipe is a pipe-like IO object, which converts into base64 data sent
# to a stream. (You must close the pipe to complete the encode, separate from
# closing the target stream).  We also have a function base64encode(f,
# args...) which works like sprint except that it produces
# base64-encoded data, along with base64encode(args...)  which is equivalent
# to base64encode(write, args...), to return base64 strings.
# A Base64DecodePipe object can be used to decode base64-encoded data read from a stream
# , while function base64decode is useful for decoding strings
#############################################################################

mutable struct Base64EncodePipe <: IO
    io::IO
    # writing works in groups of 3, so we need to cache last two bytes written
    b0::UInt8
    b1::UInt8
    nb::UInt8 # number of bytes in cache: 0, 1, or 2

    function Base64EncodePipe(io::IO)
        b = new(io,0,0,0)
        finalizer(b, close)
        return b
    end
end

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

# Based on code by Stefan Karpinski from https://github.com/hackerschool/WebSockets.jl (distributed under the same MIT license as Julia)

# '-' and '_' instead of '+' and '/'
const b64chars = ['A':'Z';'a':'z';'0':'9';'-';'_']

const base64_pad = UInt8('=')

function b64(x::UInt8, y::UInt8, z::UInt8)
    n = Int(x)<<16 | Int(y)<<8 | Int(z)
    b64chars[(n >> 18) + 1],
    b64chars[(n >> 12) & 0b111111 + 1],
    b64chars[(n >> 6) & 0b111111 + 1],
    b64chars[(n ) & 0b111111 + 1]
end

function b64(x::UInt8, y::UInt8)
    a, b, c = b64(x, y, 0x0)
    a, b, c, base64_pad
end

function b64(x::UInt8)
    a, b = b64(x, 0x0, 0x0)
    a, b, base64_pad, base64_pad
end

const sentinel = typemax(UInt8)
const revb64chars = fill(sentinel, 256)
# Fill revb64chars
for (val, ch) in enumerate(b64chars)
    revb64chars[UInt8(ch)] = UInt8(val - 1)
end

# Decode a block of at least 2 and at most 4 bytes, received in encvec
# Returns the first decoded byte and stores up to two more in cache
function b64decode!(encvec::Vector{UInt8}, cache::Vector{UInt8})
    if length(encvec) < 2
        throw(ArgumentError("incorrect base64 format, block must be at least 2 and at most 4 bytes"))
    end
    @inbounds u = revb64chars[encvec[1]]
    @inbounds v = revb64chars[encvec[2]]
    empty!(cache)
    res = (u << 2) | (v >> 4)
    if length(encvec) > 2
        @inbounds w = revb64chars[encvec[3]]
        push!(cache, (v << 4) | (w >> 2))
    end
    if length(encvec) > 3
        @inbounds z = revb64chars[encvec[4]]
        push!(cache, (w << 6) | z)
    end
    res
end

function write(b::Base64EncodePipe, x::UInt8)
    if b.nb == 0
        b.b0 = x
        b.nb = 1
    elseif b.nb == 1
        b.b1 = x
        b.nb = 2
    else
        write(b.io, b64(b.b0,b.b1,x)...)
        b.nb = 0
    end
    1
end

function close(b::Base64EncodePipe)
    if b.nb > 0
        # write leftover bytes + padding
        if b.nb == 1
            write(b.io, b64(b.b0)...)
        else # b.nb == 2
            write(b.io, b64(b.b0, b.b1)...)
        end
        b.nb = 0
    end
    nothing
end

# like sprint, but returns base64 string
function base64urlencode(f::Function, args...)
    s = IOBuffer()
    b = Base64EncodePipe(s)
    f(b, args...)
    close(b)
    String(take!(s))
end
base64urlencode(x...) = base64urlencode(write, x...)

#############################################################################
mutable struct Base64DecodePipe <: IO
    io::IO
    # reading works in blocks of 4 characters that are decoded into 3 bytes and 2 of them cached
    cache::Vector{UInt8}
    encvec::Vector{UInt8}

    function Base64DecodePipe(io::IO)
        b = new(io,[],[])
        finalizer(b, close)
        return b
    end
end

function read(b::Base64DecodePipe, t::Type{UInt8})
    if !isempty(b.cache)
        return shift!(b.cache)
    else
        empty!(b.encvec)
        while !eof(b.io) && length(b.encvec) < 4
            c::UInt8 = read(b.io, t)
            @inbounds if revb64chars[c] != sentinel
                push!(b.encvec, c)
            end
        end
        return b64decode!(b.encvec,b.cache)
    end
end

eof(b::Base64DecodePipe) = isempty(b.cache) && eof(b.io)
close(b::Base64DecodePipe) = nothing

# Decodes a base64-encoded string
function base64urldecode(s)
    b = IOBuffer(s)
    try
        return read(Base64DecodePipe(b))
    finally
        close(b)
    end
end

end # module


text:
"""

#### References

An incremental extreme learning machine for online sequential learning problems.

Lu Guo, Jing-hua Hao, and Min Liu.

Neurocomputing, 2014 vol. 128 pp. 50-58.

http://linkinghub.elsevier.com/retrieve/pii/S0925231213010059

"""
mutable struct LSIELM{TN,TV} <: AbstractSLFN where TN <: AbstractNodeInput where TV AbstractArray{Float64}
    p::Int  # Number of training points
    q::Int  # Dimensionality of function domain
    s::Int  # number of neurons
    neuron_type::TN
    c::Float64  # regularization parameter
    μx::Vector{Float64}
    σx::Vector{Float64}
    Wt::Matrix{Float64}  # transpose of W matrix
    d::Vector{Float64}

    M::Matrix{Float64}  # internal
    p_tot::Int
    v::TV

    function LSIELM(p::Int, q::Int, s::Int, neuron_type::TN, c::Float64,
                    μx, σx) where TN <: AbstractNodeInput
        # NOTE: Wt and d are constant throughout all learning chunks
        Wt = 2*rand(q, s) .- 1
        d = rand(s)
        M = zeros(Float64, s, s)
        new{TN,typeof(d)}(p, q, s, neuron_type, c, μx, σx, Wt, d, M, 0)
    end
end

function LSIELM(x::AbstractArray, y::AbstractArray;
                c::Float64=1/(50*eps()),
                neuron_type::AbstractNodeInput=Linear(Tanh()),
                s::Int=size(x, 1))
    q = size(x, 2)  # dimensionality of function domain
    p = size(x, 1)  # number of training points
    s = min(p, s)   # Can't have more neurons than training points
    xn, μx, σx = standardize(x[:, :])
    out = LSIELM(p, q, s, neuron_type, c, μx, σx)
    fit!(out, xn, y)
end

## API methods
function fit!(elm::LSIELM, x::AbstractArray, y::AbstractArray)
    # check if we are boosting or not
    S = hidden_out(elm, x)
    if elm.p_tot == 0  # this is the initialization phase
        # elm.v = S \ y
        elm.M = inv(1/elm.c .+ S'S)
        elm.v = elm.M*S'y
    else
        K = I - elm.M * S' * inv(S*elm.M*S' + I)*S
        elm.v = K*(elm.v + elm.M*S'y)
        elm.M = K*elm.M
    end

    elm.p_tot += size(x, 1)
    elm

end

function Base.show(io::IO, elm::LSIELM{TA}) where TA
    s =
    """
    LSIELM with
      - $(TA) Activation function
      - $(elm.q) input dimension(s)
      - $(elm.s) neuron(s)
      - $(elm.p) initial training point(s)
      - $(elm.p_tot) total training points
    """
    print(io, s)
end

function test_me(::Type{LSIELM})
    x = range(0, stop=1, length=20)
    y = sin.(6*x[:])
    oselm = LSIELM(x, y)

    xt = range(0, stop=1, length=300)
    yt = sin.(6*xt)

    @show maximum(abs.(oselm(xt) - yt))

    for i in 1:50
        x2 = rand(15)
        y2 = sin.(6*x2)
        fit!(oselm, x2, y2)
        @printf "%i\t%2.4e\n" i maximum(abs.(oselm(xt) - yt))
    end
end


text:
import Base

mutable struct Broadcasting{F} <: DiffableFunction
    f :: F
    true_func :: DiffableFunction
    grad_field :: GradField
end

Broadcasting(f, true_func) = Broadcasting(f, true_func, GradField())

_func_to_jt_struct = Dict(
    (+) => Add,
    (-) => Sub,
    (*) => Mul,
    (/) => Div,
    (^) => Pow,
    (sin) => Sin,
    (cos) => Cos,
    (tan) => Tan,
    (log) => Log,
    (exp) => Exp
)

Base.broadcasted(::typeof(|>), x::Variable, f) = Base.broadcasted(f, x)

# fix https://github.com/abap34/JITrench.jl/issues/11
Base.broadcasted(::typeof(Base.literal_pow), ::typeof(^), x::Variable, ::Val{c}) where c = Broadcasting(^, Pow(GradField(), c))(x)
Base.broadcasted(::typeof(^), x::Variable, c) = Broadcasting(^, Pow(GradField(), c))(x)

for (func, jt_func) in _func_to_jt_struct
    (func == ^) && (continue)
    Base.broadcasted(::typeof(func), x::Variable...) = Broadcasting(func, jt_func(GradField()))(x...)
    Base.broadcasted(::typeof(func), x1::Variable, x2) = Broadcasting(func, jt_func(GradField()))(x1, Variable(x2))
    Base.broadcasted(::typeof(func), x1, x2::Variable) = Broadcasting(func, jt_func(GradField()))(Variable(x1), x2)
end

function forward(f::Broadcasting, x...)
    Base.materialize(Base.broadcasted(f.f, x...))
end

function forward(f::Broadcasting{typeof(^)}, x...)
    return Base.materialize(Base.broadcasted(f.f, x..., f.true_func.c))
end

function backward(f::Broadcasting, gy)
    f.true_func.grad_field = f.grad_field
    if length(f.grad_field.inputs) ==  1
        gx = backward(f.true_func, gy)
        return gx
    else
        x1 = f.grad_field.inputs[1]
        x2 = f.grad_field.inputs[2]
        gx1, gx2 = backward(f.true_func, gy)
        if size(x1.values) != size(x2.values)
            gx1 = sum_to(gx1, size(x1.values))
            gx2 = sum_to(gx2, size(x2.values))
        end
        return gx1, gx2
    end
end


text:
export @unionise

"""
    unionise_arg(arg)

Return an expression in which the argument expression `arg` is replaced with an argument
whos type admits `Node`s.
"""
unionise_arg(arg::Symbol) = arg

function unionise_arg(arg::Expr)
    if arg.head === :(::)
        Expr(:(::), arg.args[1:end-1]..., unionise_type(arg.args[end]))
    elseif arg.head === :...
        Expr(:..., unionise_arg(arg.args[1]))
    elseif arg.head === :kw
        Expr(:kw, unionise_arg(arg.args[1]), arg.args[2])
    elseif arg.head === :parameters
        arg  # Ignore keyword arguments and leave them untouched for now
    else
        throw(ArgumentError("Unrecognized argument in Symbol ($arg)."))
    end
end

"""
    unionise_subtype(arg::Union{Symbol, Expr})

Equivalent to `unionise_arg`, but replacing `::` with `<:`.
"""
unionise_subtype(arg::Symbol) = arg
unionise_subtype(arg::Expr) =
    arg.head == Symbol("<:") ?
        Expr(Symbol("<:"), arg.args[1:end-1]..., unionise_type(arg.args[end])) :
        throw(ArgumentError("Unrecognised argument in arg ($arg)."))

"""
    get_quote_body(code)

Return the body of a quoted expression or symbol. There is slightly different behaviour
depending upon whether a symbol or expression is quoted and how exactly it is quoted,
thus we dispatch on whether code is a `:quote` expression or a `QuoteNode` object.
"""
get_quote_body(code::Expr) = code.args[1]
get_quote_body(code::QuoteNode) = code.value

"""
    unionise_eval(code::Expr)

Unionise the code inside a call to `eval`, such that when the `eval` call actually occurs
the code inside will be unionised.
"""
function unionise_eval(code::Expr, linfo::LineNumberNode=LineNumberNode(0))
    body = Expr(:macrocall, Symbol("@unionise"), linfo, deepcopy(get_quote_body(code.args[end])))
    return length(code.args) == 3 ?
        Expr(:call, :eval, deepcopy(code.args[2]), quot(body)) :
        Expr(:call, :eval, quot(body))
end

"""
    unionise_macro_eval(code::Expr)

Unionise the code in a call to @eval, such that when the `eval` call actually occurs, the
code inside will be unionised.
"""
function unionise_macro_eval(code::Expr, linfo::LineNumberNode=LineNumberNode(0))
    body = Expr(:macrocall, Symbol("@unionise"), linfo, deepcopy(code.args[end]))
    return length(code.args) == 4 ?
        Expr(:macrocall, Symbol("@eval"), linfo, deepcopy(code.args[3]), body) :
        Expr(:macrocall, Symbol("@eval"), linfo, body)
end

"""
    unionise_sig(code)

`code` should be a `tuple`, `call` or `where` expression (containing a `tuple` or `call`).
Returns a signature which accepts `Node`s.
"""
unionise_sig(code::Symbol) = code
function unionise_sig(code::Expr)
    body = get_body(code)
    if body.head == :tuple
        new_body = Expr(:tuple, unionise_arg.(body.args)...)
    elseif body.head == :call
        new_body = Expr(:call, body.args[1], unionise_arg.(body.args[2:end])...)
    elseif body.head == Symbol("::")
        new_body = unionise_arg(body)
    end
    return replace_body(code, new_body)
end

"""
    unionise_struct(code)

`code` should be an `Expr` containing the definition of a type. The type will only be
changed if it is parametric. That is, if you wish to be able to differentiate through a
user-defined type, it must contain only `Any`s and parametric types.
"""
function unionise_struct(code::Expr)
    tmp = code.args[2]
    is_subtype_expr = tmp isa Expr && tmp.head == Symbol("<:")
    name = is_subtype_expr ? code.args[2].args[1] : code.args[2]
    if name isa Expr && name.head == :curly
        curly = Expr(:curly, name.args[1], unionise_subtype.(name.args[2:end])...)
        if is_subtype_expr
            return Expr(
                :struct,
                code.args[1],
                Expr(Symbol("<:"), curly, tmp.args[2]),
                code.args[3],
            )
        else
            return Expr(:struct, code.args[1], curly, code.args[3])
        end
    else
        return code
    end
end

"""
    unionise(code)

Return transformed code in which all function definitions are guaranteed to accept nodes as
arguments. This should not affect the existing functionality of the code.
"""
function unionise end

# If we get a symbol then we cannot have found a function definition, so ignore it.
unionise(code, linfo::LineNumberNode=LineNumberNode(0)) = code

# Recurse through an expression, bottoming out if we find a function definition or a
# quoted expression to be `eval`-ed.
function unionise(code::Expr, linfo::LineNumberNode=LineNumberNode(0))
    if code.head in (:function, Symbol("->"))
        return Expr(code.head, unionise_sig(code.args[1]), code.args[2])
    elseif code.head == Symbol("=") && !isa(code.args[1], Symbol) &&
        (get_body(code.args[1]).head == :tuple || get_body(code.args[1]).head isa Symbol)
        return Expr(code.head, unionise_sig(code.args[1]), code.args[2])
    elseif code.head == :call && code.args[1] == :eval
        return unionise_eval(code, linfo)
    elseif code.head == :macrocall && code.args[1] == Symbol("@eval")
        return unionise_macro_eval(code, linfo)
    elseif code.head == :struct
        return unionise_struct(code)
    else
        return Expr(code.head, [unionise(arg, linfo) for arg in code.args]...)
    end
end

"""
    @unionise code

Transform code such that each function definition accepts `Node` objects as arguments,
without effecting dispatch in other ways.
"""
macro unionise(code)
    return esc(unionise(code, __source__))
end


text:
#=
   The mag type used by Arb (fredrikj.net/arb/mag.html)
   see also (https://github.com/Nemocas/Nemo.jl/blob/master/src/arb/ArbTypes.jl)
   see also (https://github.com/thofma/Hecke.jl/blob/master/src/Misc/mag.jl)
=#

#=
type MagFloat
    radius_exponentOf2 :: Int        # exponentOf2
    radius_significand :: UInt64     # significand
end
=#

# define hash so other things work
const hash_arbmag_lo = (UInt === UInt64) ? 0x29f934c433d9a758 : 0x2578e2ce
const hash_0_arbmag_lo = hash(zero(UInt), hash_arbmag_lo)
if UInt===UInt64
   hash(z::MagFloat, h::UInt) = hash( reinterpret(UInt64, z.radius_exponentOf2), z.radius_significand )
else
   hash(z::MagFloat, h::UInt) = hash( reinterpret(UInt32, z.radius_exponentOf2) % UInt64, z.radius_significand )
end

# initialize and zero a variable of type MagFloat
function zero(::Type{MagFloat})
    z = MagFloat()
    # ccall(@libarb(mag_set_ui), Void, (Ptr{MagFloat}, UInt64), &z, zero(UInt64))
    return z
end
#zero(x::MagFloat) = zero(MagFloat)
function zero(x::MagFloat)
    return zero(MagFloat)
end

function one(::Type{MagFloat})
    z = MagFloat()
    ccall(@libarb(mag_set_ui), Void, (Ptr{MagFloat}, UInt64), &z, one(UInt64))
    return z
end
function one(x::MagFloat)
    return one(MagFloat)
end

#=
# how to write
for (T, postfix) in ((:Int, :int), (:Float64, :fp))
    @eval  ccall(  ___ )
end
# so Julia sees
ccall( (:convert_from_int, "myClib"), Void, (Ref{ Int     }, ), my_int_va
ccall( (:convert_from_fp,  "myClib"), Void, (Ref{ Float64 }, ), my_fp_var  )
# ------- from Yichao Yu, to get
ccall( (:convert_from_int, "myClib"), Void, (Ref{ Int     }, ), my_int_var )
# ------- the loop content is written this way
ccall( ($(QuoteNode(Symbol("convert_from", postfix))), "myClib"),
Void, (Ref{$T}, ), $(Symbol("my_", postfix, "_var"))
=#

for (T,M) in ((:UInt, :ui), (:Int, :si), (:Float64, :d))
  @eval begin
    function convert(::Type{MagFloat}, x::($T))
        z = MagFloat()
        ccall( :($(QuoteNode(Symbol("mag_set_", $M))), "libarb"), Void, (Ptr{$T}, ), &z )
        #ccall( ($(QuoteNode(Symbol("mag_set_", M))), "libarb"), Void, (Ref{$T}, ), z )
        return z
    end
  end
end

MagFloat(radius_exponentOf2::Int, radius_significand::Int64) =
    MagFloat(radius_exponentOf2, radius_significand % UInt64)

MagFloat(radius_exponentOf2::Int, radius_significand::Int32) =
    MagFloat(radius_exponentOf2, UInt64(radius_significand % UInt32) )

MagFloat(radius_exponentOf2::Int, radius_significand::Float64) =
    MagFloat(radius_exponentOf2, convert(UInt64, abs(radius_significand)) )

MagFloat(radius_exponentOf2::Int, radius_significand::Float32) =
    MagFloat(radius_exponentOf2, convert(UInt64, abs(radius_significand)) )

# conversions

# convert to MagFloat

Error_MagIsNegative() = throw(ErrorException("Magnitudes must be nonnegative."))

#=
function convert(::Type{MagFloat}, x::Float64)
    signbit(x) && Error_MagIsNegative()
    z = MagFloat()
    ccall(@libarb(mag_set_d), Void, (Ptr{MagFloat}, Ptr{Float64}), &z, &x)
    return z
end
convert(::Type{MagFloat}, x::Float32) = convert(MagFloat, convert(Float64, x))
convert(::Type{MagFloat}, x::Float16) = convert(MagFloat, convert(Float64, x))
=#

#=
   lowerbound returns upper bound of value
   upperbound returns lower bound of value
=#

if Int == Int64
  function upperbound(::Type{MagFloat}, x::UInt64)
    z = MagFloat()
    ccall(@libarb(mag_set_ui), Void, (Ptr{MagFloat}, Ptr{UInt64}), &z, &x)
    return z
  end
  function lowerbound(::Type{MagFloat}, x::UInt64)
    z = MagFloat()
    ccall(@libarb(mag_set_ui_lower), Void, (Ptr{MagFloat}, Ptr{UInt64}), &z, &x)
    return z
  end
  for T in (:UInt128, :UInt32, :UInt16, :UInt8)
    @eval upperbound(::Type{MagFloat}, x::($T)) = upperbound(MagFloat, convert(UInt64, x))
    @eval lowerbound(::Type{MagFloat}, x::($T)) = lowerbound(MagFloat, convert(UInt64, x))
  end
  for T in (:Int128, :Int32, :Int16, :Int8)
    @eval convert(::Type{MagFloat}, x::($T))  = convert(MagFloat, convert(Int64, x))
  end
else
  function upperbound(::Type{MagFloat}, x::UInt32)
    z = MagFloat()
    ccall(@libarb(mag_set_ui), Void, (Ptr{MagFloat}, Ptr{UInt32}), &z, &x)
    return z
  end
  function lowerbound(::Type{MagFloat}, x::UInt32)
    z = MagFloat()
    ccall(@libarb(mag_set_ui_lower), Void, (Ptr{MagFloat}, Ptr{UInt32}), &z, &x)
    return z
  end
  for T in (:UInt128, :UInt64, :UInt16, :UInt8)
    @eval upperbound(::Type{MagFloat}, x::($T)) = upperbound(MagFloat, convert(UInt32, x))
    @eval lowerbound(::Type{MagFloat}, x::($T)) = lowerbound(MagFloat, convert(UInt32, x))
  end
  for T in (:Int128, :Int64, :Int16, :Int8)
    @eval convert(::Type{MagFloat}, x::($T))  = convert(MagFloat, convert(Int32, x))
  end
end

# convert from MagFloat

function convert(::Type{Float64}, x::MagFloat)
    z = ccall(@libarb(mag_get_d), Float64, (Ptr{MagFloat}, ), &x)
    return z
end
function convert(::Type{Float32}, x::MagFloat)
    z = convert(Float64, x)
    convert(Float32, z)
    return z
end

# function convert(::Type{UInt}, x::MagFloat)

# promotions

for T in (:UInt, :Int, :Float32, :Float64)
    @eval promote_rule(::Type{MagFloat}, ::Type{$T}) = MagFloat
end

# string, show
#
function string(x::MagFloat)
    fp = convert(Float64, x)
    return string(fp)
end

function stringcompact(x::MagFloat)
    fp = convert(Float32, convert(Float64, x))
    return string(fp)
end

function show(io::IO, x::MagFloat)
    s = string(x)
    print(io, s)
    return nothing
end

function showcompact(io::IO, x::MagFloat)
    s = stringcompact(x)
    print(io, s)
    return nothing
end


text:
filenames = Dict{String, String}()

tests = ["lasertag",
         "simpleld",
         "subhunt",
         "vdpbarrier"]

for t in tests
    try
        println("\n$(uppercase(t))\n")
        include(t*"_table.jl")
        filenames[t] = filename
    catch ex
        showerror(STDERR, ex)
    end
end

# for dt in ["subhunt_discretization.jl", "ld_discretization.jl"]
#     try
#         println("\n$(uppercase(dt))\n")
#         include(dt)
#         filenames[dt] = filename
#     catch ex
#         showerror(STDERR, ex)
#     end
# end

println("\nFILENAMES\n")
for (k, v) in filenames
    println("$k => $v")
end

@show filenames


text:
import JSON
using Gadfly

adaptLog = JSON.parsefile("adapt3.log")
notAdaptLog = JSON.parsefile("notAdapt3.log")

p = plot(layer(x = map(x -> x["id"], adaptLog),
               y = map(x -> mean(x["results"]), adaptLog),
               xintercept = [10, 50],
               Geom.point,
               Geom.vline,
               Theme(highlight_width = 0pt)),

         layer(x = map(x -> x["id"], notAdaptLog),
               y = map(x -> mean(x["results"]), notAdaptLog),
               Geom.point,
               Theme(default_color=colorant"red", highlight_width = 0pt)),

         Scale.x_continuous(minvalue=0, maxvalue=60),

         Scale.y_log10,

         Guide.xlabel("Time (s)"),

         Guide.ylabel("Response time (ms)"),

         Guide.manual_color_key("Legend",
                                ["Self-adaptive", "Non-self-adaptive"],
                                ["deepskyblue", "red"]))

p2 = plot(
         layer(x = map(x -> x["id"], adaptLog),
               y = map(x->length(x["results"]), adaptLog),
               Geom.bar,
               Theme(bar_spacing=0.5mm)),

         Guide.xlabel("Time (s)"),

         Guide.ylabel("Number of requests"),

         Scale.x_continuous(minvalue=0, maxvalue=70),

         Scale.y_continuous(minvalue=0, maxvalue=30)
         )

draw(PDF("myplot1.pdf", 10inch, 5inch), p)

draw(PDF("myplot2.pdf", 10inch, 5inch), p2)

max_stars_repo_name: prl-tokyo/configuration-bigul

max_issues_repo_name: prl-tokyo/configuration-bigul

max_forks_repo_name: prl-tokyo/configuration-bigul


text:
using ConstrainedDynamics
using LinearAlgebra
using ForwardDiff
using Rotations

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1

p2 = [0.0;0.0;length1/2] # joint connection point

# Links
origin = Origin{Float64}()
link1 = Box(width, depth, length1, length1)
link2 = deepcopy(link1)

# Constraints
joint1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=p2))
joint2 = EqualityConstraint(Revolute(link1, link2, joint_axis; p1=-p2, p2=p2))

links = [link1;link2]
constraints = [joint1;joint2]

mech = Mechanism(origin, links, constraints, g=-9.81)
setPosition!(origin,link1,p2 = p2)
setPosition!(link1,link2,p1=-p2,p2 = p2)

xd=[-[p2];-[p2+p2+p2]]
qd=[[UnitQuaternion(RotX(0.0))];[UnitQuaternion(RotX(0.0))]]

A, Bu, Bλ, G = ConstrainedDynamics.linearsystem(mech, xd, [zeros(3) for i=1:2], qd, [zeros(3) for i=1:2], [[0.] for i=1:1], getid.(mech.bodies), [getid(constraints[2])])

steps = Base.OneTo(1000)
storage = Storage{Float64}(steps,length(mech.bodies))

ang01 = 0.1
ang02 = 0.0
setPosition!(origin,link1,p2 = p2,Δq = UnitQuaternion(RotX(ang01)))
setPosition!(link1,link2,p1=-p2,p2 = p2,Δq = UnitQuaternion(RotX(ang02)))
state1 = link1.state
state2 = link2.state
q01=[state1.xc-xd[1]; 0;0;0; ConstrainedDynamics.Vmat() * Rotations.params(state1.qc); 0;0;0]
q02=[state2.xc-xd[2]; 0;0;0; ConstrainedDynamics.Vmat() * Rotations.params(state2.qc); 0;0;0]
q0=[q01;q02]
simulate!(mech,storage,record = true)

function testfoo(vars)
    q0 = vars[1:24]
    q1 = vars[25:48]
    λ = vars[49:58]

    [
        A*q0 + Bλ*λ - q1
        G*(A*q0 + Bλ*λ)
    ]
end

function testnewton(q0)
    q1 = q0
    λ = zeros(10)
    for i=1:100
        vars = [q0;q1;λ]
        dqλ = -inv(ForwardDiff.jacobian(testfoo,vars)[:,25:58])*testfoo(vars)
        qλnew = [q0;λ] + dqλ
        varsnew = [q0;qλnew]
        for j=1:10
            if norm(testfoo(vars)) > norm(testfoo(varsnew))
                q0 = qλnew[1:24]
                λ = qλnew[25:34]
                break
            else
                qλnew = [q0;λ] + (1/2^(j))*dqλ
                varsnew = [q0;qλnew]
            end
        end

        if norm(testfoo(varsnew)) < 1e-7
            return q0, λ
        end
    end
    # display("newton not converged")
    return
end

Q = zeros(24,1000)
Q[:,1] = q0
for i=2:1000
    Q[:,i] = testnewton(Q[:,i-1])[1]
end

diff = getindex.(storage.x[1],2)-Q[2,:]
@test maximum(diff) < 0.01
diff = getindex.(storage.x[2],2)-Q[14,:]
@test maximum(diff) < 0.01


text:
using Chemfiles
using Test

TESTS = [
    "Atom.jl",
    "Residue.jl",
    "Topology.jl",
    "UnitCell.jl",
    "Frame.jl",
    "Property.jl",
    "Trajectory.jl",
    "Selection.jl",
    "misc.jl",
]

include("utils.jl")

function main()
    @testset "Generics" begin
        @test split(Chemfiles.version(), '-')[1] == "0.9.1"
    end
    root = dirname(@__FILE__)
    for test in TESTS
        include(test)
    end
end

main()


text:
# using Pkg
# Pkg.add("DelimitedFiles")

using DelimitedFiles

data = readdlm("../externals/Core.Math.Data/data/Pejcic_318.csv", ',')

print(data)

# println((data),"\n")


text:
mutable struct DeviceMenu <: TerminalMenus.AbstractMenu
    options::Vector{Schema.DeviceInfo}
    pagesize::Int
    indent::Int
    pageoffset::Int
    selected::Int
end

function DeviceMenu(devices::Vector{Schema.DeviceInfo}; pagesize::Int=5, warn::Bool=true)
    length(devices) < 1 && error("DeviceMenu must have at least one option")
    pagesize < 4 && error("minimum pagesize must be larger than 4")
    pagesize = pagesize == -1 ? length(devices) : pagesize
    # pagesize shouldn't be bigger than options
    pagesize = min(length(devices), pagesize)
    # after other checks, pagesize must be greater than 1
    pagesize < 1 && error("pagesize must be >= 1")

    pageoffset = 0
    selected = -1 # none
    indent = maximum(x->length(x.backend_name), devices) + 2
    return DeviceMenu(devices, pagesize, indent, pageoffset, selected)
end

function TerminalMenus.options(m::DeviceMenu)
    return m.options
end

TerminalMenus.cancel(m::DeviceMenu) = m.selected = -1

function TerminalMenus.pick(m::DeviceMenu, cursor::Int)
    m.selected = cursor
    return true #break out of the menu
end

function TerminalMenus.printmenu(out::IO, m::DeviceMenu, cursoridx::Int; oldstate=nothing, init::Bool=false)
    buf = IOBuffer()
    lastoption = length(m.options)
    ncleared = oldstate === nothing ? m.pagesize-1 : oldstate

    if init
        # like clamp, except this takes the min if max < min
        m.pageoffset = max(0, min(cursoridx - m.pagesize ÷ 2, lastoption - m.pagesize))
    else
        print(buf, "\x1b[999D\x1b[$(ncleared)A")   # move left 999 spaces and up `ncleared` lines
    end

    firstline = m.pageoffset+1
    lastline = min(m.pagesize+m.pageoffset, lastoption)
    curr_device = m.options[cursoridx]

    for i in firstline:lastline
        # clearline
        print(buf, "\x1b[2K")

        upscrollable = i == firstline && m.pageoffset > 0
        downscrollable = i == lastline && i != lastoption

        if upscrollable && downscrollable
            print(buf, TerminalMenus.updown_arrow(m)::Union{Char,String})
        elseif upscrollable
            print(buf, TerminalMenus.up_arrow(m)::Union{Char,String})
        elseif downscrollable
            print(buf, TerminalMenus.down_arrow(m)::Union{Char,String})
        else
            print(buf, ' ')
        end

        # TODO: use 1.6's implementation when we drop 1.5
        # TerminalMenus.printcursor(buf, m, i == cursoridx)
        print(buf, i == cursoridx ? '→' : ' ', ' ')

        name = m.options[i].backend_name
        indent = " "^(m.indent - length(name))
        device = m.options[cursoridx]

        line_idx = i - firstline + 1
        print(buf, GREEN_FG(name))

        if line_idx == 1
            print(buf, indent, LIGHT_BLUE_FG("nqubits: "), GREEN_FG(string(device.n_qubits)))
        elseif line_idx == 2
            print(buf, indent, LIGHT_BLUE_FG("max_shots: "), GREEN_FG(string(device.max_shots)))
        elseif line_idx == 3
            print(buf, indent, LIGHT_BLUE_FG("credits_required:: "), GREEN_FG(string(device.credits_required)))
        end

        (firstline == lastline || i != lastline) && print(buf, "\r\n")
    end

    newstate = lastline-firstline  # final line doesn't have `\n`
    if newstate < ncleared && oldstate !== nothing
        # we printed fewer lines than last time. Erase the leftovers.
        for i = newstate+1:ncleared
            print(buf, "\r\n\x1b[2K")
        end
        print(buf, "\x1b[$(ncleared-newstate)A")
    end

    print(out, String(take!(buf)))

    return newstate
end


text:
module onnx_utils

# Converter Utils
function UInt8toFloat32(val)
    dims = val.dims
    data = val.raw_data
    indices = collect(1:4:length(data))
    data_list = []
    for i in indices; push!(data_list, float_from_bitstring(Float32, generate_bitstring(data[i:i+3]))); end;
    if length(val.dims) != 0
        new_size = tuple(val.dims...)
        data_list = reshape(data_list, new_size)
    end
    #cast to Float32 (migth adjust this later)
    return Float32.(data_list)
end

function generate_bitstring(raw_4)
    bits = ""
    for i in reverse(raw_4); bits *= bitstring(i); end
    bits
end

function float_from_bitstring(::Type{T}, str::String) where {T<:Base.IEEEFloat}
    unsignedbits = Meta.parse(string("0b", str))
    thefloat  = reinterpret(T, unsignedbits)
    return thefloat
end

export UInt8toFloat32
end


text:
using CSV
using DataFrames
CSV.read(joinpath(@__DIR__,"sample.csv"), DataFrame)

max_stars_repo_name: jondea/juliab

max_issues_repo_name: jondea/juliab

max_forks_repo_name: jondea/juliab


text:
# These tests don't do much except ensure that the package loads,
# and does something sensible if it does.
using Loess
using Test
using Random
using Statistics

@testset " multivariate test case " begin

    Random.seed!(100)
    xs = 10 .* rand(100)
    ys = hcat(sin.(xs) .+ 0.5 * rand(100),
              cos.(xs) .+ 0.5 * rand(100))

    model = loess(xs, ys)

    us = collect(minimum(xs):0.1:maximum(xs))
    vs = predict(model, us)

    @test minimum(vs) >= -1.1
    @test maximum(vs) <= +1.1

end

@testset " 1d test case " begin

    Random.seed!(100)
    xs = 10 .* rand(100)
    ys = sin.(xs) .+ 0.5 * rand(100)

    model = loess(xs, ys)

    us = collect(minimum(xs):0.1:maximum(xs))
    vs = predict(model, us)

    @test minimum(vs) >= -1.1
    @test maximum(vs) <= +1.1

    x = [13.0,14.0,14.35,15.0,16.0]
    y = [0.369486,  0.355579, 0.3545, 0.356952, 0.36883]
    model = loess(x,y)
    @test Loess.predict(model,x) ≈ y

    let x = 1:10, y = sin.(1:10)
        model = loess(x, y)
        @test model.xs == reshape(collect(Float64, 1:10), 10, 1)
        @test model.ys == y
        pred = [1.02866, 0.798561, 0.533528, 0.253913, -0.0325918, -0.319578, -0.648763]
        @test predict(model, 1.0:0.5:4.0) ≈ pred atol=1e-5
    end

    @test_throws DimensionMismatch loess([1.0 2.0; 3.0 4.0], [1.0])

end


text:
@testset "Orbitals" begin
    @set basis cc-pvdz
    orbs = Fermi.Orbitals.GeneralRestrictedOrbitals(rand(5,5), sd_energy=-50.0)
    @test orbs.basis == "cc-pvdz" && orbs.sd_energy == -50.0 && size(orbs.C) == (5,5)
end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
abstract type Link <: Device end
abstract type Pipe <: Link end
abstract type Valve <: Link end


text:
(pwd() != @__DIR__) && cd(@__DIR__) # allow starting app from bin/ dir

using SwarmServer
push!(Base.modules_warned_for, Base.PkgId(SwarmServer))
SwarmServer.main()


text:
module GiftExchangeRandomizer

using Combinatorics: permutations

struct Exchanger
  id::Int
  name::String
end

struct ExchangeRelationship
  giver::Int
  receiver::Int
end

ExchangeRelationship(array) = ExchangeRelationship(array...)

exchangers() = [
                Exchanger(1, "Ed Baldwin"),
                Exchanger(2, "Karen Baldwin"),
                Exchanger(3, "Gordo Stevens"),
                Exchanger(4, "Tracy Stevens"),
                Exchanger(5, "Ellen Wilson"),
                Exchanger(6, "Larry Wilson"),
                Exchanger(7, "Molly Cobb"),
                Exchanger(8, "Wayne Cobb")
               ]

households() = [[1, 2], [3, 4], [5, 6], [7, 8]]

last_time() = [[1, 4], [2, 8], [3, 6], [4, 7], [5, 1], [6, 3], [7, 2], [8, 5]]
time_before_last() =
  [[1, 4], [2, 6], [3, 2], [4, 8], [5, 7], [6, 1], [7, 5], [8, 3]]

function relationships_from_array(exchanger_ids)
  relationships = Vector{ExchangeRelationship}(undef, length(exchanger_ids))
  for (i, id) = pairs(exchanger_ids)
    next_i = i + 1 <= length(exchanger_ids) ? i + 1 : 1
    relationships[i] = ExchangeRelationship(id, exchanger_ids[next_i])
  end
  relationships
end

function possible_relationships(exchanger_ids)
  relationships_from_array.(permutations(exchanger_ids)) |>
  Iterators.flatten |>
  Set
end

function household_relationships()
  union(possible_relationships.(households())...)
end

function possible_arrangements(exchanger_ids)
  relationships_from_array.(permutations(exchanger_ids))
end

function valid_arrangements(exchanger_ids)
  household_rels = household_relationships()
  has_household_rel(arrangement) =
    length(intersect(arrangement, household_rels)) != 0

  filter(!has_household_rel, possible_arrangements(exchanger_ids))
end

function scored_arrangements(exchanger_ids)
  arrangements = valid_arrangements(exchanger_ids)
  last_time_rels = Set(map(ExchangeRelationship, last_time()))
  time_before_last_rels = Set(map(ExchangeRelationship, time_before_last()))

  num_last_time_rels(arrangement) =
    length(intersect(arrangement, last_time_rels))
  num_time_before_last_rels(arrangement) =
    length(intersect(arrangement, time_before_last_rels))

  scored = Dict{Int, Vector{Vector{ExchangeRelationship}}}()
  for arrangement in arrangements
    score = (1000
             - (num_last_time_rels(arrangement) * 20)
             - (num_time_before_last_rels(arrangement) * 10))
    if !haskey(scored, score)
      scored[score] = []
    end
    push!(scored[score], arrangement)
  end
  scored
end

function next_arrangement(exchanger_ids)
  scored = scored_arrangements(exchanger_ids)
  max_score = max(keys(scored)...)

  arrangement = rand(scored[max_score])

  (max_score, arrangement)
end

end

max_stars_repo_name: rosston/GiftExchangeRandomizer

max_issues_repo_name: rosston/GiftExchangeRandomizer

max_forks_repo_name: rosston/GiftExchangeRandomizer


text:
using Turf

fc = GeoJSON.parsefile(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/getting-started/example.geojson")

centroid_point = centroid(fc)
center_point = center(fc)

push!(fc.features, Feature(centroid_point, Dict("marker-color" => "#fff000")))
push!(fc.features, Feature(center_point, Dict("marker-color" => "#eaa000")))

result = GeoJSON.geojson(fc)

open(joinpath(dirname(pathof(Turf)), "..", "docs") * "/examples/getting-started/example.result.geojson", "w") do file
    write(file, result)
end


text:
using OREnvironment
using Test

const testdir = dirname(@__FILE__)

tests = [
    "Status"
    "Solution"
    "Constraints"
    "Problem"
    "Movements"
]

@testset "OREnvironment.jl" begin
    for t in tests
        @info "Testing $t"
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
# import Base: getindex, setindex!, firstindex, lastindex, iterate, keys, isempty
struct BiDictMap{T <: Integer}
    int_sym::Dict{T,Symbol}
    sym_int::Dict{Symbol,T}
end

BiDictMap{T}(;sizehint=100) where T<:Integer = begin
    int_sym = Dict{T,Symbol}()
    sizehint!(int_sym, sizehint)
    sym_int = Dict{Symbol,T}()
    sizehint!(sym_int, sizehint)
    BiDictMap{T}(int_sym, sym_int)
end

BiDictMap(;sizehint=100) = BiDictMap{Int}(;sizehint=sizehint)

Base.getindex(b::BiDictMap, key::Int) = b.int_sym[key]
Base.getindex(b::BiDictMap, key::Symbol) = b.sym_int[key]

# setindex!(b, value, key) = b[key] = value
function Base.setindex!(b::BiDictMap, s::Symbol, i::Int)
    haskey(b.sym_int, s) && delete!(b.int_sym, b[s])
    haskey(b.int_sym, i) && delete!(b.sym_int, b[i])

    b.int_sym[i] = s
    b.sym_int[s] = i
end

function Base.setindex!(b::BiDictMap, i::Int, s::Symbol)
    haskey(b.int_sym, i) && delete!(b.sym_int, b[i])
    haskey(b.sym_int, s) && delete!(b.int_sym, b[s])

    b.int_sym[i] = s
    b.sym_int[s] = i
end

function Base.delete!(b::BiDictMap, i::Int)
    s = b[i]
    delete!(b.int_sym, i)
    delete!(b.sym_int, s)
    return b
end

Base.haskey(b::BiDictMap, s::Symbol) = haskey(b.sym_int, s)
Base.haskey(b::BiDictMap, i::Int) = haskey(b.int_sym, i)

Base.length(b::BiDictMap) = length(b.int_sym)
#NOTE to future self: This will work only with LightGraphs that always follows a 1:nv(g) index
Base.firstindex(v::BiDictMap) = 1
Base.lastindex(v::BiDictMap) = length(v.int_sym)
Base.iterate(v::BiDictMap, i=1) = (length(v.int_sym) < i ? nothing : (v.int_sym[i], i + 1))
Base.keys(v::BiDictMap) = Base.OneTo(length(v.int_sym))
Base.isempty(v::BiDictMap) = (length(v.int_sym) == 0)


text:
function calcmax(xa::Float64, p::Int64)
  λarr = collect(linspace(0.0, 4.0, p))
  const iterconst = 50000
  const ϵ = 1e-6
  xarr = Array{Array{Float64,1},1}(p)
  yarr = Array{Array{Float64,1},1}(p)
  max = 0
  for j = 1:p
    λ = λarr[j]
    x = xa
    xarr[j] = Array{Float64,1}(0)
    yarr[j] = Array{Float64,1}(0)
    for i = 1:iterconst
      x = λ*x*(1-x)
    end
    push!(yarr[j], x)
    found = false
    while !found
      x = λ*x*(1-x)
      for y in yarr[j]
        if almostequal(x, y, ϵ)
          found = true
          break
        end
      end
      if !found
        push!(yarr[j], x)
      end
    end
    push!(xarr[j], λ)
    max = (length(yarr[j]) > max ? length(yarr[j]) : max)
  end
  return max
end

function findeq(k::Int64)
  xa = 0.5
  arr = Vector{Int64}(k)
  for i = 1:k
    arr[i] = calcmax(xa, i*1000)
  end
  return arr
end

function timecalc(xa::Float64, p::Int64)
  _,_ = calcarr(xa, p)
  println("Array of arrays:")
  @time calcarr(xa, p)
  _,_ = calcnan(xa, p)
  println("NaN:")
  @time calcnan(xa, p)
  nothing
end

function calcarrp(xa::Float64, p::Int64)
  λarr = collect(linspace(0.0, 4.0, p))
  const iterconst = 50000
  const ϵ = 1e-6
  xarr = Array{Array{Float64, 1}, 1}(p)
  yarr = Array{Array{Float64, 1}, 1}(p)
  points = 0
  for j = 1:p
    λ = λarr[j]
    x = xa
    xarr[j] = Array{Float64, 1}(0)
    yarr[j] = Array{Float64, 1}(0)
    for i = 1:iterconst
      x = λ*x*(1-x)
    end
    push!(yarr[j], x)
    found = false
    while !found
      x = λ*x*(1-x)
      for y in yarr[j]
        if almostequal(x, y, ϵ)
          found = true
          break
        end
      end
      if !found
        push!(yarr[j], x)
      end
    end
    push!(xarr[j], λ)
    points += length(yarr[j])
  end
  println(points)
end

function calcnanp(xa::Float64, p::Int64, start = 0.0, fin = 4.0)
  λarr = collect(linspace(start, fin, p))
  const pconst = 2720 # calculated with findeq()
  const iterconst = 50000
  const ϵ = 1e-6
  yarr = Matrix{Float64}(p, pconst)
  fill!(yarr, NaN)
  points = 0
  for i = 1:p
    λ = λarr[i]
    x = xa
    for j = 1:iterconst
      x = λ*x*(1-x)
    end
    yarr[i, 1] = x # add first value
    m = 1
    found = false
    while !found # while we still have new values
      x = λ*x*(1-x) # iterate once more
      for k = 1:m # iterate through all found values
        if almostequal(x, yarr[i, k], ϵ) # and check if we already have the new one
          found = true
          break
        end
      end
      if !found
        m += 1
        yarr[i, m] = x # if not, add it to the array
      end
    end
    points += m
  end
  println(points)
end

function calc_numofpoints(xa::Float64, p::Int64)
  calcarrp(xa, p)
  calcnanp(xa, p)
end


text:
# Called in the OrdinaryDiffEQ.__init; All `OrdinaryDiffEqAlgorithm`s have one
function OrdinaryDiffEq.initialize!(integ, cache::GaussianODEFilterCache)
    if integ.opts.dense && !integ.alg.smooth
        error("To use `dense=true` you need to set `smooth=true`!")
    elseif !integ.opts.dense && integ.alg.smooth
        @warn "If you set dense=false for efficiency, you might also want to set smooth=false."
    end
    if !integ.opts.save_everystep && integ.alg.smooth
        error("If you do not save all values, you do not need to smooth!")
    end
    @assert integ.saveiter == 1

    integ.kshortsize = 1
    resize!(integ.k, integ.kshortsize)
    integ.k[1] = integ.u

    # Update the initial state to the known (given or computed with AD) initial values
    initial_update!(integ, cache, integ.alg.initialization)

    # These are necessary since the solution object is not 100% initialized by default
    OrdinaryDiffEq.copyat_or_push!(integ.sol.x_filt, integ.saveiter, cache.x)
    OrdinaryDiffEq.copyat_or_push!(
        integ.sol.pu,
        integ.saveiter,
        mul!(cache.pu_tmp, cache.SolProj, cache.x),
    )
    return nothing
end

"""Perform a step

Not necessarily successful! For that, see `step!(integ)`.

Basically consists of the following steps
- Coordinate change / Predonditioning
- Prediction step
- Measurement: Evaluate f and Jf; Build z, S, H
- Calibration; Adjust prediction / measurement covs if the diffusion model "dynamic"
- Update step
- Error estimation
- Undo the coordinate change / Predonditioning
"""
function OrdinaryDiffEq.perform_step!(
    integ,
    cache::GaussianODEFilterCache,
    repeat_step=false,
)
    @unpack t, dt = integ
    @unpack d, SolProj = integ.cache
    @unpack x, x_pred, u_pred, x_filt, u_filt, err_tmp = integ.cache
    @unpack x_tmp, x_tmp2 = integ.cache
    @unpack A, Q, Ah, Qh = integ.cache

    make_preconditioners!(cache, dt)
    @unpack P, PI = integ.cache

    tnew = t + dt

    # Build the correct matrices
    @. Ah = PI.diag .* A .* P.diag'
    X_A_Xt!(Qh, Q, PI)

    if isdynamic(cache.diffusionmodel)  # Calibrate, then predict cov

        # Predict
        predict_mean!(x_pred, x, Ah)
        mul!(view(u_pred, :), SolProj, x_pred.μ)

        # Measure
        evaluate_ode!(integ, x_pred, tnew)

        # Estimate diffusion
        cache.local_diffusion, cache.global_diffusion =
            estimate_diffusion(cache.diffusionmodel, integ)
        # Adjust prediction and measurement
        predict_cov!(x_pred, x, Ah, Qh, cache.C1, cache.global_diffusion)

        # Compute measurement covariance only now
        compute_measurement_covariance!(cache)

    else
        predict_mean!(x_pred, x, Ah)
        predict_cov!(x_pred, x, Ah, Qh, cache.C1)
        mul!(view(u_pred, :), SolProj, x_pred.μ)
        evaluate_ode!(integ, x_pred, tnew)
        compute_measurement_covariance!(cache)
        cache.local_diffusion, cache.global_diffusion =
            estimate_diffusion(cache.diffusionmodel, integ)
    end

    # Likelihood
    # cache.log_likelihood = logpdf(cache.measurement, zeros(d))

    # Update
    x_filt = update!(integ, x_pred)

    # Estimate error for adaptive steps - can already be done before filtering
    if integ.opts.adaptive
        err_est_unscaled = estimate_errors(cache)
        if integ.f isa DynamicalODEFunction # second-order ODE
            DiffEqBase.calculate_residuals!(
                err_tmp,
                dt * err_est_unscaled,
                integ.u[1, :],
                u_pred[1, :],
                integ.opts.abstol,
                integ.opts.reltol,
                integ.opts.internalnorm,
                t,
            )
        else # regular first-order ODE
            DiffEqBase.calculate_residuals!(
                err_tmp,
                dt * err_est_unscaled,
                integ.u,
                u_pred,
                integ.opts.abstol,
                integ.opts.reltol,
                integ.opts.internalnorm,
                t,
            )
        end
        integ.EEst = integ.opts.internalnorm(err_tmp, t) # scalar
    end

    # If the step gets rejected, we don't even need to perform an update!
    reject = integ.opts.adaptive && integ.EEst >= one(integ.EEst)
    if !reject

        # Save into u_filt and integ.u
        mul!(view(u_filt, :), SolProj, x_filt.μ)
        integ.u .= u_filt

        # Advance the state here
        copy!(integ.cache.x, integ.cache.x_filt)
        integ.sol.log_likelihood += integ.cache.log_likelihood
    end
end

function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    x_pred,
    t,
    second_order::Val{false},
)
    @unpack f, p, dt, alg = integ
    @unpack u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    # Mean
    _eval_f!(du, u_pred, p, t, f)
    integ.destats.nf += 1
    # z .= MM*E1*x_pred.μ .- du
    if f.mass_matrix == I
        H .= E1
    elseif f.mass_matrix isa UniformScaling
        H .= f.mass_matrix.λ .* E1
    else
        _matmul!(H, f.mass_matrix, E1)
    end

    _matmul!(z, H, x_pred.μ)
    z .-= du[:]

    # If EK1, evaluate the Jacobian and adjust H
    if alg isa EK1 || alg isa IEKS
        u_lin =
            (alg isa IEKS && !isnothing(alg.linearize_at)) ? alg.linearize_at(t).μ : u_pred

        # Jacobian is computed either with the given jac, or ForwardDiff
        if !isnothing(f.jac)
            _eval_f_jac!(ddu, u_lin, p, t, f)
        elseif isinplace(f)
            ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u_lin)
            integ.destats.nf += 1
        else
            ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u_lin)
            integ.destats.nf += 1
        end
        integ.destats.njacs += 1

        # _matmul!(H, f.mass_matrix, E1) # This is already the case (see above)
        _matmul!(H, ddu, E0, -1.0, 1.0)
    end

    return nothing
end

function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:EK1FDB},
    x_pred,
    t,
    second_order::Val{false},
)
    @unpack f, p, dt, alg = integ
    @unpack d, u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    (f.mass_matrix != I) && error("EK1FDB does not support mass-matrices right now")

    # Mean
    _eval_f!(du, u_pred, p, t, f)
    integ.destats.nf += 1
    # z .= MM*E1*x_pred.μ .- du
    H1, H2 = view(H, 1:d, :), view(H, d+1:2d, :)
    z1, z2 = view(z, 1:d), view(z, d+1:2d)

    H1 .= E1
    _matmul!(z1, H1, x_pred.μ)
    z1 .-= @view du[:]

    # If EK1, evaluate the Jacobian and adjust H
    u_lin = u_pred
    if !isnothing(f.jac)
        _eval_f_jac!(ddu, u_lin, p, t, f)
    elseif isinplace(f)
        ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u_lin)
        integ.destats.nf += 1
    else
        ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u_lin)
        integ.destats.nf += 1
    end
    integ.destats.njacs += 1
    _matmul!(H1, ddu, E0, -1.0, 1.0)

    z2 .= (E2 * x_pred.μ .- ddu * du)
    if integ.alg.jac_quality == 1
        # EK0-type approach
        H2 .= E2
    elseif integ.alg.jac_quality == 2
        H2 .= E2 - ddu * ddu * E0
    elseif integ.alg.jac_quality == 3
        _z2(m) = begin
            u_pred = E0 * m
            du = zeros(eltype(m), d)
            ddu = zeros(eltype(m), d, d)
            _eval_f!(du, u_pred, p, t, f)
            if !isnothing(f.jac)
                _eval_f_jac!(ddu, u_pred, p, t, f)
            elseif isinplace(f)
                ForwardDiff.jacobian!(ddu, (du, u) -> f(du, u, p, t), du, u_pred)
                # integ.destats.nf += 1
            else
                ddu .= ForwardDiff.jacobian(u -> f(u, p, t), u_pred)
                # integ.destats.nf += 1
            end
            # integ.destats.njacs += 1
            return (E2 * m .- ddu * du)
        end
        H2 .= ForwardDiff.jacobian(_z2, x_pred.μ)
    else
        error("EK1FDB's `jac_quality` has to be in [1,2,3]")
    end
    return nothing
end

function evaluate_ode!(
    integ::OrdinaryDiffEq.ODEIntegrator{<:AbstractEK},
    x_pred,
    t,
    second_order::Val{true},
)
    @unpack f, p, dt, alg = integ
    @unpack d, u_pred, du, ddu, measurement, R, H = integ.cache
    @assert iszero(R)
    du2 = du

    @unpack E0, E1, E2 = integ.cache

    z, S = measurement.μ, measurement.Σ

    # Mean
    # _u_pred = E0 * x_pred.μ
    # _du_pred = E1 * x_pred.μ
    if isinplace(f)
        f.f1(du2, view(u_pred, 1:d), view(u_pred, d+1:2d), p, t)
    else
        du2 .= f.f1(view(u_pred, 1:d), view(u_pred, d+1:2d), p, t)
    end
    integ.destats.nf += 1
    _matmul!(z, E2, x_pred.μ)
    z .-= @view du2[:]

    # Cov
    if alg isa EK1
        (alg isa IEKS) && error("IEKS is currently not supported for SecondOrderODEProbems")

        if isinplace(f)
            H .= E2

            J = ddu
            ForwardDiff.jacobian!(
                J,
                (du2, du_u) -> f.f1(du2, view(du_u, 1:d), view(du_u, d+1:2d), p, t),
                du2,
                u_pred,
            )
            integ.destats.nf += 1
            integ.destats.njacs += 1
            _matmul!(H, J, integ.cache.SolProj, -1.0, 1.0)
        else
            J = ForwardDiff.jacobian(
                (du_u) -> f.f1(view(du_u, 1:d), view(du_u, d+1:2d), p, t),
                u_pred,
            )
            integ.destats.nf += 1
            integ.destats.njacs += 1
            H .= E2 .- J * integ.cache.SolProj
        end
    end

    return measurement
end
evaluate_ode!(integ, x_pred, t) =
    evaluate_ode!(integ, x_pred, t, Val(integ.f isa DynamicalODEFunction))

# The following functions are just there to handle both IIP and OOP easily
_eval_f!(du, u, p, t, f::AbstractODEFunction{true}) = f(du, u, p, t)
_eval_f!(du, u, p, t, f::AbstractODEFunction{false}) = (du .= f(u, p, t))
_eval_f_jac!(ddu, u, p, t, f::AbstractODEFunction{true}) = f.jac(ddu, u, p, t)
_eval_f_jac!(ddu, u, p, t, f::AbstractODEFunction{false}) = (ddu .= f.jac(u, p, t))

compute_measurement_covariance!(cache) =
    X_A_Xt!(cache.measurement.Σ, cache.x_pred.Σ, cache.H)

function update!(integ, prediction)
    @unpack measurement, H, R, x_filt = integ.cache
    @unpack K1, K2, x_tmp2, m_tmp = integ.cache
    update!(x_filt, prediction, measurement, H, K1, x_tmp2.Σ.mat, m_tmp)
    return x_filt
end

function smooth_all!(integ)
    integ.sol.x_smooth = copy(integ.sol.x_filt)

    @unpack A, Q = integ.cache
    @unpack x_smooth, t, diffusions = integ.sol
    @unpack x_tmp, x_tmp2 = integ.cache
    x = x_smooth

    for i in length(x)-1:-1:1
        dt = t[i+1] - t[i]
        if iszero(dt)
            copy!(x[i], x[i+1])
            continue
        end

        make_preconditioners!(integ.cache, dt)
        P, PI = integ.cache.P, integ.cache.PI

        mul!(x_tmp, P, x[i])
        mul!(x_tmp2, P, x[i+1])
        smooth!(x_tmp, x_tmp2, A, Q, integ, diffusions[i])
        mul!(x[i], PI, x_tmp)
    end
end

function estimate_errors(cache::GaussianODEFilterCache)
    @unpack local_diffusion, Qh, H, d = cache

    if local_diffusion isa Real && isinf(local_diffusion)
        return Inf
    end

    L = cache.m_tmp.Σ.squareroot

    if local_diffusion isa Diagonal
        _matmul!(L, H, sqrt.(local_diffusion) * Qh.squareroot)
        error_estimate = sqrt.(diag(L * L'))
        return view(error_estimate, 1:d)

    elseif local_diffusion isa Number
        _matmul!(L, H, Qh.squareroot)
        # error_estimate = local_diffusion .* diag(L*L')
        @tullio error_estimate[i] := L[i, j] * L[i, j]
        error_estimate .*= local_diffusion

        # @info "it's small anyways I guess?" error_estimate cache.measurement.μ .^ 2
        # error_estimate .+= cache.measurement.μ .^ 2

        error_estimate .= sqrt.(error_estimate)
        return view(error_estimate, 1:d)
    end
end


text:
export TriMesh,
    GBMesh,
    gbmeshes,
    load_trimesh,
    save_trimesh,
    get_verts_list,
    get_verts_packed,
    get_verts_padded,
    get_faces_list,
    get_faces_packed,
    get_faces_padded,
    get_edges_packed,
    get_edges_to_key,
    get_faces_to_edges_packed,
    get_laplacian_packed,
    get_laplacian_sparse,
    compute_verts_normals_list,
    compute_verts_normals_packed,
    compute_verts_normals_padded,
    compute_faces_normals_list,
    compute_faces_normals_packed,
    compute_faces_normals_padded,
    compute_faces_areas_list,
    compute_faces_areas_packed,
    compute_faces_areas_padded

import GeometryBasics, Printf, MeshIO
import GeometryBasics:
    Point3f, GLTriangleFace, NgonFace, convert_simplex, meta, triangle_mesh, value

import Zygote: @ignore

"""
    TriMesh

Initialize Triangle Mesh representation.

### Fields:

- `N`                       - Batch size of TriMesh.
- `V`                       - Maximum vertices per mesh in TriMesh.
- `F`                       - Maximum faces per mesh in TriMesh.
- `equalised`               - Bool, indicates all mesh have same verts and faces size.
- `valid`                   -
- `offset`                  - Offset indicating number to be added to migrate to 0-indexed system.

- `_verts_len`              - Number of vertices in each mesh of TriMesh.
- `_verts_list`             - Vertices in list format.
- `_verts_packed`           - Vertices in packed format.
- `_verts_padded`           - Vertices in padded format.

- `_faces_len`              - Number of faces in each mesh of TriMesh.
- `_faces_list`             - Vertices in list format.
- `_faces_packed`           - Vertices in packed format.
- `_faces_padded`           - Vertices in padded format.

- `_edges_packed`           - Edges in packed format (according to packed vertices).
- `_faces_to_edges_packed`  - Faces formed by edges in packed format (according to packed edges).
- `_laplacian_packed`       - Laplacian sparce matrix in packed format.
- `_edges_to_key`           - Dict mapping edges tuple to unique key.

### Available Contructor:

- `TriMesh(verts_list, faces_list; offset::Number = -1)`
- `TriMesh(m::Vector{<:GeometryBasics.Mesh})`
- `TriMesh(m::GeometryBasics.Mesh)`
- `TriMesh(m::TriMesh)`

"""
mutable struct TriMesh{T<:AbstractFloat,R<:Integer,S} <: AbstractObject
    N::Int64
    V::Int64
    F::Int64
    equalised::Bool
    valid::BitArray{1}
    offset::Int8
    _verts_len::Vector{Int64}
    _faces_len::Vector{Int64}

    _verts_packed::S
    _verts_padded::S
    _verts_list::Vector{<:S}
    _verts_packed_valid::Bool
    _verts_padded_valid::Bool
    _verts_list_valid::Bool

    _faces_packed::Array{R,2}
    _faces_padded::Array{R,3}
    _faces_list::Vector{Array{R,2}}
    _faces_packed_valid::Bool
    _faces_padded_valid::Bool

    _edges_packed::Union{Nothing,Array{R,2}}
    _faces_to_edges_packed::Union{Nothing,Array{R,2}}
    _laplacian_packed::Union{Nothing,AbstractSparseArray{T,R,2}}

    _edges_to_key::Union{Nothing,Dict{Tuple{R,R},R}}
end

# TriMesh(
#     verts::Vector{<:AbstractArray{T,2}},
#     faces::Vector{<:AbstractArray{R,2}};
#     offset::Number = -1,
# ) where {T<:Number,R<:Number} =
#     TriMesh([Float32.(v) for v in verts], [UInt32.(f) for f in faces]; offset = offset)

TriMesh(
    verts::Vector{<:CuArray{T,2}},
    faces::Vector{<:AbstractArray{R,2}};
    offset::Number = -1,
) where {T<:AbstractFloat,R<:Integer} = TriMesh(CuArray, verts, faces; offset = offset)

TriMesh(
    verts::Vector{<:Array{T,2}},
    faces::Vector{<:AbstractArray{R,2}};
    offset::Number = -1,
) where {T<:AbstractFloat,R<:Integer} = TriMesh(Array, verts, faces; offset = offset)

function TriMesh(
    ::Type{S},
    verts::Vector{<:AbstractArray{T,2}},
    faces::Vector{<:AbstractArray{R,2}};
    offset::Number = -1,
) where {S,T<:AbstractFloat,R<:Integer}

    length(verts) == length(faces) || error(
        "batch size of verts and faces should match, $(length(verts)) != $(length(faces))",
    )

    # To remove lazy wrappers on verts and faces list
    verts = [T.(v) for v in verts]
    faces = [R.(f) for f in faces]

    _verts_len = size.(verts, 2)
    _faces_len = size.(faces, 2)

    N = length(verts)
    V = maximum(_verts_len)
    F = maximum(_faces_len)
    equalised = all(_verts_len .== V) && all(_faces_len .== F)
    valid = _faces_len .> 0
    offset = Int8(offset)

    _verts_list = verts::Vector{<:S{T,2}}
    _faces_list = faces::Vector{Array{R,2}}

    return TriMesh{T,R,S}(
        N,
        V,
        F,
        equalised,
        valid,
        offset,
        _verts_len,
        _faces_len,
        S{T,2}(undef, 3, sum(_verts_len)),
        S{T,3}(undef, 3, V, N),
        _verts_list,
        false,
        false,
        true,
        Array{R,2}(undef, 3, sum(_faces_len)),
        Array{R,3}(undef, 3, F, N),
        _faces_list,
        false,
        false,
        nothing,
        nothing,
        nothing,
        nothing,
    )
end

function TriMesh(ms::Vector{<:GeometryBasics.Mesh})
    verts_list = Array{Float32,2}[]
    faces_list = Array{UInt32,2}[]
    for (i, m) in enumerate(ms)
        (verts, faces) = _load_meta(m)
        push!(verts_list, verts)
        push!(faces_list, faces)
    end
    return TriMesh(verts_list, faces_list)
end

TriMesh(m::GeometryBasics.Mesh) = TriMesh([m])
TriMesh(m::TriMesh) = TriMesh(get_verts_list(m), get_faces_list(m))
TriMesh(v::AbstractArray{2}, f::AbstractArray{2}) = TriMesh([v], [f])

functor(x::TriMesh) = (_verts_list = x._verts_list,),
xs -> TriMesh(xs._verts_list, x._faces_list; offset = x.offset)

function Base.show(io::IO, m::TriMesh{T,R,S}) where {T,R,S}
    print(
        io,
        "TriMesh{$(T), $(R), $(S)} Structure:\n    Batch size: ",
        m.N,
        "\n    Max verts: ",
        m.V,
        "\n    Max faces: ",
        m.F,
        "\n    offset: ",
        m.offset,
        "\n    Storage type: ",
        S,
    )
end

function Base.setproperty!(m::TriMesh, f::Symbol, v)
    if (f == :_verts_packed) || (f == :_verts_padded) || (f == :_verts_list)
        # only attempt to change if given object is not same as object already
        # present inside struct, which can be checked with !== operator.
        # List is always valid, so to make sure that we refresh list always
        if (f == :_verts_packed) && (getproperty(m, f) !== v)
            setfield!(m, f, convert(fieldtype(typeof(m), f), v))
            setfield!(m, :_verts_padded_valid, false)
            setfield!(m, :_verts_list_valid, false)
            _compute_verts_list(m, true)
        elseif (f == :_verts_padded) && (getproperty(m, f) !== v)
            setfield!(m, f, convert(fieldtype(typeof(m), f), v))
            setfield!(m, :_verts_packed_valid, false)
            setfield!(m, :_verts_list_valid, false)
            _compute_verts_list(m, true)
        elseif (f == :_verts_list) && (getproperty(m, f) !== v)
            setfield!(m, f, convert(fieldtype(typeof(m), f), v))
            setfield!(m, :_verts_packed_valid, false)
            setfield!(m, :_verts_padded_valid, false)
        end
    else
        setfield!(m, f, convert(fieldtype(typeof(m), f), v))
    end
end

Base.getindex(m::TriMesh, inds...) =
    (get_verts_list(m)[inds...], get_faces_list(m)[inds...])

"""
    GBMesh(m::TriMesh; index::Int = 1)
    GBMesh(verts::AbstractArray{T,2}, faces::AbstractArray{R,2}) where {T,R}

Initialize GeometryBasics.Mesh from triangle mesh in TriMesh `m` at `index`.

See also: [`gbmeshes`](@ref)
"""
function GBMesh(verts::AbstractArray{T,2}, faces::AbstractArray{R,2}) where {T,R}
    points = Point3f[GeometryBasics.Point{3,Float32}(verts[:, i]) for i = 1:size(verts, 2)]
    verts_dim = size(faces, 1)
    poly_face = NgonFace{verts_dim,UInt32}[
        NgonFace{verts_dim,UInt32}(faces[:, i]) for i = 1:size(faces, 2)
    ]
    # faces = convert_simplex.(GLTriangleFace, poly_face)
    faces = GLTriangleFace.(poly_face)
    return GeometryBasics.Mesh(meta(points), faces)
end

GBMesh(m::TriMesh, index::Int = 1) = GBMesh(m[index]...)

"""
    gbmeshes(m::TriMesh)

Initialize list of GeometryBasics.Mesh from TriMesh `m`

See also: [`gbmeshes`](@ref)
"""
function gbmeshes(m::TriMesh)
    gbmeshes = [GBMesh(m, i) for i = 1:m.N]
    return gbmeshes
end

"""
    _load_meta(m::GeometryBasics.Mesh)

Returns `vertices` and `faces` in Array format.

"""
function _load_meta(m::GeometryBasics.Mesh)
    if !(m isa GeometryBasics.Mesh{3,Float32,<:GeometryBasics.Triangle})
        m = triangle_mesh(m)
    end
    vs = m.position
    vertices = [Array(v) for v in vs]
    vertices = reduce(hcat, vertices)
    fs = getfield(getfield(m, :simplices), :faces)
    faces = [Array(UInt32.(value.(f))) for f in fs]
    faces = reduce(hcat, faces)
    return (vertices, faces)
end

"""
    load_trimesh(fn::String...)

Load TriMesh from file(s).
It will load TriMesh with multiple meshes, if multiple files are given.

Supported formats are `obj`, `stl`, `ply`, `off` and `2DM`.

"""
function load_trimesh(fns::String...; elements_types...)
    verts_list = Array{Float32,2}[]
    faces_list = Array{UInt32,2}[]
    for (i, fn) in enumerate(fns)
        mesh = load(fn; elements_types...)
        verts, faces = _load_meta(mesh)
        push!(verts_list, verts)
        push!(faces_list, faces)
    end
    return TriMesh(verts_list, faces_list)
end

"""
    save_trimesh(fn::String, mesh::TriMesh, index::Int = 1)
    save_trimesh(fn::String, mesh::GeometryBasics.Mesh)

Save mesh in given `fn`.
`index` is an optional argument specifing the index of mesh,
incase of multiple meshes in TriMesh `mesh`.

Supported formats are `obj`, `stl`, `ply`, `off` and `2DM`.

"""
save_trimesh(fn::String, mesh::TriMesh, index::Int = 1) =
    save_trimesh(fn, GBMesh(mesh, index))

function save_trimesh(fn::String, mesh::GeometryBasics.Mesh)
    save(fn, mesh)
end

"""
    get_verts_packed(m::TriMesh; refresh::Bool = false)

Returns vertices of TriMesh `m` in packed format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_verts_padded`](@ref), [`get_verts_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_verts_packed(m)
```
"""
function get_verts_packed(m::TriMesh{T,R,S}; refresh::Bool = false)::S{T,2} where {T,R,S}
    _compute_verts_packed(m, refresh)
    return m._verts_packed
end

"""
    get_verts_padded(m::TriMesh; refresh::Bool = false)

Returns vertices of TriMesh `m` in padded format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_verts_packed`](@ref), [`get_verts_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_verts_padded(m)
```
"""
function get_verts_padded(m::TriMesh{T,R,S}; refresh::Bool = false)::S{T,3} where {T,R,S}
    _compute_verts_padded(m, refresh)
    return m._verts_padded
end

"""
    get_verts_list(m::TriMesh; refresh::Bool = false)

Returns vertices of TriMesh `m` in list format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_verts_padded`](@ref), [`get_verts_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_verts_list(m)
```
"""
function get_verts_list(
    m::TriMesh{T,R,S};
    refresh::Bool = false,
)::Vector{<:S{T,2}} where {T,R,S}
    _compute_verts_list(m, refresh)
    return m._verts_list
end

"""
    get_faces_packed(m::TriMesh; refresh::Bool = false)

Returns faces of TriMesh `m` in packed format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_faces_padded`](@ref), [`get_faces_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_faces_packed(m)
```
"""
function get_faces_packed(m::TriMesh{T,R}; refresh::Bool = false)::Array{R,2} where {T,R}
    @ignore _compute_faces_packed(m, refresh)
    return m._faces_packed
end

"""
    get_faces_padded(m::TriMesh; refresh::Bool = false)

Returns faces of TriMesh `m` in padded format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_faces_padded`](@ref), [`get_faces_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_faces_padded(m)
```
"""
function get_faces_padded(m::TriMesh{T,R}; refresh::Bool = false)::Array{R,3} where {T,R}
    @ignore _compute_faces_padded(m, refresh)
    return m._faces_padded
end

"""
    get_faces_list(m::TriMesh; refresh::Bool = false)

Returns faces of TriMesh `m` in list format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_faces_padded`](@ref), [`get_faces_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_faces_list(m)
```
"""
function get_faces_list(
    m::TriMesh{T,R};
    refresh::Bool = false,
)::Vector{Array{R,2}} where {T,R}
    return m._faces_list
end

"""
    get_edges_packed(m::TriMesh; refresh::Bool = false)

Returns edges of TriMesh `m` in packed format. Edges are according
to the indices of corresponding vertices in `get_verts_packed(m)`

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_verts_packed`](@ref), [`get_faces_to_edges_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_edges_packed(m)
```
"""
function get_edges_packed(m::TriMesh{T,R}; refresh::Bool = false)::Array{R,2} where {T,R}
    @ignore _compute_edges_packed(m, refresh)
    return m._edges_packed
end

"""
    get_edges_to_key(m::TriMesh; refresh::Bool = false)

Returns dict mapping edges (tuple) of TriMesh `m` to unique key.
Edges are according to the indices of corresponding vertices in`get_verts_packed(m)`

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`get_verts_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_edges_to_key(m)
```
"""
function get_edges_to_key(
    m::TriMesh{T,R};
    refresh::Bool = false,
)::Dict{Tuple{R,R},R} where {T,R}
    @ignore _compute_edges_packed(m, refresh)
    return m._edges_to_key
end

"""
    get_faces_to_edges_packed(m::TriMesh; refresh::Bool = false)

Returns faces of TriMesh `m` in form of edges.
Each edge corresponds to the indices of `get_edges_packed(m)`.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`. Assuming face f consists of
(v1,v2,v3) vertices, and e1 = {v2,v3}, e2 = {v3,v1}, e3 = {v1,v2},
so face f in form of edges would be (e1,e2,e3).

See also: [`get_edges_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_faces_to_edges_packed(m)
```
"""
function get_faces_to_edges_packed(
    m::TriMesh{T,R};
    refresh::Bool = false,
)::Array{R,2} where {T,R}
    @ignore _compute_edges_packed(m, refresh)
    return m._faces_to_edges_packed
end

"""
    get_laplacian_packed(m::TriMesh; refresh::Bool = false)
    get_laplacian_sparse(m::TriMesh; refresh::Bool = false)

Returns Laplacian sparce matrix of TriMesh `m`.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`laplacian_loss`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> get_laplacian_packed(m)
```
"""
function get_laplacian_packed(
    m::TriMesh{T,R};
    refresh::Bool = false,
)::SparseMatrixCSC{T,R} where {T,R}
    _compute_laplacian_packed(m, refresh)
    return m._laplacian_packed
end

const get_laplacian_sparse = get_laplacian_packed

"""
    compute_verts_normals_packed(m::TriMesh)

Computes Unit normal of vertices of TriMesh `m` in packed format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.
Normal vec of each vertex is weighted sum of normals
of adjacent faces, weighted by corresponding faces areas
and then normalize to unit vector.

See also: [`compute_verts_normals_padded`](@ref), [`compute_verts_normals_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_verts_normals_packed(m)
```
"""
function compute_verts_normals_packed(m::TriMesh{T,R,S})::S{T,2} where {T,R,S}

    verts = get_verts_packed(m)
    faces = get_faces_packed(m)

    vert_faces = verts[:, faces]
    vertex_normals = similar(verts, size(verts))
    vertex_normals = @ignore fill!(vertex_normals, 0.0)
    vertex_normals = Zygote.bufferfrom(vertex_normals)

    # normal of each vertex is sum of normals of its faces weighted by
    # corresponding areas, ie. (A1 .* fn1) + (A2 .* fn2) + (A2 .* fn2)
    # where fn1, fn2, fn3 are normals of faces shared by that vertex
    # and A1, A2, A3 are corresponding areas.

    vertex_normals[:, faces[1, :]] += _lg_cross(
        vert_faces[:, 2, :] - vert_faces[:, 1, :],
        vert_faces[:, 3, :] - vert_faces[:, 1, :],
    )
    vertex_normals[:, faces[2, :]] += _lg_cross(
        vert_faces[:, 3, :] - vert_faces[:, 2, :],
        vert_faces[:, 1, :] - vert_faces[:, 2, :],
    )
    vertex_normals[:, faces[3, :]] += _lg_cross(
        vert_faces[:, 1, :] - vert_faces[:, 3, :],
        vert_faces[:, 2, :] - vert_faces[:, 3, :],
    )

    return _normalize(copy(vertex_normals), dims = 1)
end

"""
    compute_verts_normals_padded(m::TriMesh)

Computes Unit normal of vertices of TriMesh `m` in padded format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.
Normal vec of each vertex is weighted sum of normals
of adjacent faces, weighted by corresponding faces areas
and then normalize to unit vector.

See also: [`compute_verts_normals_packed`](@ref), [`compute_verts_normals_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_verts_normals_padded(m)
```
"""
function compute_verts_normals_padded(m::TriMesh{T,R,S})::S{T,3} where {T,R,S}
    normals_packed = compute_verts_normals_packed(m)
    normals_padded = _packed_to_padded(normals_packed, m._verts_len, 0.0)
    return normals_padded
end

"""
    compute_verts_normals_list(m::TriMesh)

Computes Unit normal of vertices of TriMesh `m` in list format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.
Normal vec of each vertex is weighted sum of normals
of adjacent faces, weighted by corresponding faces areas
and then normalize to unit vector.

See also: [`compute_verts_normals_padded`](@ref), [`compute_verts_normals_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_verts_normals_list(m)
```
"""
function compute_verts_normals_list(m::TriMesh{T,R,S})::Vector{<:S{T,2}} where {T,R,S}
    normals_packed = compute_verts_normals_packed(m)
    normals_list = _packed_to_list(normals_packed, m._verts_len)
    return normals_list
end

"""
    compute_faces_normals_packed(m::TriMesh)

Computes Unit normal of faces of TriMesh `m` in packed format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_normals_padded`](@ref), [`compute_faces_normals_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_normals_packed(m)
```
"""
function compute_faces_normals_packed(m::TriMesh{T,R,S})::S{T,2} where {T,R,S}
    verts = get_verts_packed(m)
    faces = get_faces_packed(m)

    vert_faces = verts[:, faces]
    # normal vec of face [f1, f2, f3] is (f2-f1)X(f3-f1)
    face_normals = _lg_cross(
        vert_faces[:, 2, :] - vert_faces[:, 1, :],
        vert_faces[:, 3, :] - vert_faces[:, 1, :],
    )
    return _normalize(face_normals, dims = 1)
end

"""
    compute_faces_normals_padded(m::TriMesh)

Computes Unit normal of faces of TriMesh `m` in padded format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_normals_packed`](@ref), [`compute_faces_normals_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_normals_padded(m)
```
"""
function compute_faces_normals_padded(m::TriMesh{T,R,S})::S{T,3} where {T,R,S}
    normals_packed = compute_faces_normals_packed(m)
    normals_padded = _packed_to_padded(normals_packed, m._faces_len, 0.0)
    return normals_padded
end

"""
    compute_faces_normals_list(m::TriMesh)

Computes Unit normal of faces of TriMesh `m` in list format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_normals_padded`](@ref), [`compute_faces_normals_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_normals_list(m)
```
"""
function compute_faces_normals_list(m::TriMesh{T,R,S})::Vector{<:S{T,2}} where {T,R,S}
    normals_packed = compute_faces_normals_packed(m)
    normals_list = _packed_to_list(normals_packed, m._faces_len)
    return normals_list
end

"""
    compute_faces_areas_packed(m::TriMesh)

Computes area of faces of TriMesh `m` in packed format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_areas_padded`](@ref), [`compute_faces_areas_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_areas_packed(m)
```
"""
function compute_faces_areas_packed(
    m::TriMesh{T,R,S};
    eps::Number = 1e-6,
)::S{T,1} where {T,R,S}
    verts = get_verts_packed(m)
    faces = get_faces_packed(m)

    vert_faces = verts[:, faces]
    face_normals_vec = _lg_cross(
        vert_faces[:, 2, :] - vert_faces[:, 1, :],
        vert_faces[:, 3, :] - vert_faces[:, 1, :],
    )
    face_norm = _norm(face_normals_vec; dims = 1)
    face_areas = dropdims(face_norm ./ 2; dims = 1)
    return face_areas
end

"""
    compute_faces_areas_padded(m::TriMesh)

Computes area of faces of TriMesh `m` in padded format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_areas_packed`](@ref), [`compute_faces_areas_list`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_areas_padded(m)
```
"""
function compute_faces_areas_padded(
    m::TriMesh{T,R,S};
    eps::Number = 1e-6,
)::S{T,3} where {T,R,S}
    areas_packed = compute_faces_areas_packed(m; eps = eps)
    # increasing ndims to 2, to support _packed_to_list
    areas_packed = reshape(areas_packed, 1, :)
    areas_padded = _packed_to_padded(areas_packed, m._faces_len, 0.0)
    return areas_padded
end

"""
    compute_faces_areas_list(m::TriMesh)

Computes area of faces of TriMesh `m` in list format.

`refresh` is an optional keyword argument to
recompute this format for TriMesh `m`.

See also: [`compute_faces_areas_padded`](@ref), [`compute_faces_areas_packed`](@ref)

### Examples:

```jldoctest
julia> m = load_trimesh("teapot.obj")
julia> compute_faces_areas_list(m)
```
"""
function compute_faces_areas_list(
    m::TriMesh{T,R,S};
    eps::Number = 1e-6,
)::Vector{<:S{T,2}} where {T,R,S}
    areas_packed = compute_faces_areas_packed(m; eps = eps)
    # increasing ndims to 2, to support _packed_to_list
    areas_packed = reshape(areas_packed, 1, :)
    areas_list = _packed_to_list(areas_packed, m._faces_len)
    return areas_list
end

function _compute_verts_packed(m::TriMesh, refresh::Bool = false)
    if refresh || !(m._verts_packed_valid)
        _compute_verts_list(m)
        verts_packed = _list_to_packed(m._verts_list)
        # avoiding setproperty!, as we are building packed
        # from list and list is always valid
        setfield!(m, :_verts_packed, verts_packed)
        setfield!(m, :_verts_packed_valid, true)
        return nothing
    end
end

function _compute_verts_padded(m::TriMesh, refresh::Bool = false)
    if refresh || !(m._verts_padded_valid)

        # FIXME: #38 chamfer_loss gradient doesn't reach packed, so using temp fix
        _compute_verts_packed(m)
        verts_padded = _packed_to_padded(m._verts_packed, m._verts_len, 0)
        # m._verts_padded = _padded

        # _list_to_padded!(m._verts_padded, get_verts_list(m), 0, (3, m.V))

        # avoiding setproperty!, as we are building padded
        # from list and list is always valid
        setfield!(m, :_verts_padded, verts_padded)
        setfield!(m, :_verts_padded_valid, true)
        return nothing
    end
end

function _compute_verts_list(m::TriMesh, refresh::Bool = false)
    if refresh || !(m._verts_list_valid)
        if m._verts_packed_valid
            verts_list = _packed_to_list(m._verts_packed, m._verts_len)
        elseif m._verts_padded_valid
            verts_list = _padded_to_list(m._verts_padded, m._verts_len)
        else
            error("not possible to contruct list without padded and packed")
        end
        # avoiding setproperty, cause verts_list is always valid.
        setfield!(m, :_verts_list, verts_list)
        setfield!(m, :_verts_list_valid, true)
        return nothing
    end
end

function _compute_faces_packed(m::TriMesh{T,R}, refresh::Bool = false) where {T,R}
    if refresh || !(m._faces_packed_valid)
        faces_packed = _list_to_packed(m._faces_list)
        _, verts_packed_first_idx, _ = _auxiliary_mesh(m._verts_list)
        _, _, faces_packed_list_idx = _auxiliary_mesh(m._faces_list)
        # offset will be equal to first idx in verts_packed of corresponding faces_idx-1
        faces_packed_offset = verts_packed_first_idx[faces_packed_list_idx] .- 1
        faces_packed = faces_packed .+ reshape(faces_packed_offset, 1, :)
        setfield!(m, :_faces_packed, R.(faces_packed))
        setfield!(m, :_faces_packed_valid, true)
        return nothing
    end
end

function _compute_faces_padded(m::TriMesh, refresh::Bool = false)
    if refresh || !(m._faces_padded_valid)
        _list_to_padded!(m._faces_padded, m._faces_list, 0, (3, m.F))
        # setfield!(m, :_faces_padded, faces_padded)
        setfield!(m, :_faces_padded_valid, true)
        return nothing
    end
end

function _compute_edges_packed(m::TriMesh{T,R}, refresh::Bool = false) where {T,R}
    if refresh ||
       (any([m._edges_packed, m._edges_to_key, m._faces_to_edges_packed] .=== nothing))

        faces = get_faces_packed(m)
        verts = get_verts_packed(m)

        e12 = cat(faces[1, :], faces[2, :], dims = 2)
        e23 = cat(faces[2, :], faces[3, :], dims = 2)
        e31 = cat(faces[3, :], faces[1, :], dims = 2)

        # Sort edges (v0, v1) such that v0 <= v1
        e12 = sort(e12; dims = 2)
        e23 = sort(e23; dims = 2)
        e31 = sort(e31; dims = 2)

        # Edges including duplicates
        edges = cat(e12, e23, e31, dims = 1)

        # Converting edge (v0, v1) into integer hash, ie. (V+1)*v0 + v1.
        # There will be no collision, which is asserted by (V+1), as 1<=v0<=V.
        V_hash = size(verts, 2) + 1
        edges_hash = (V_hash .* edges[:, 1]) .+ edges[:, 2]

        # Sort and remove duplicate edges_hash
        sort!(edges_hash)
        unique!(edges_hash)

        # Convert edges_hash to edges
        edges = cat((edges_hash .÷ V_hash), (edges_hash .% V_hash); dims = 2)

        # Edges to key
        edges_to_key =
            Dict{Tuple{R,R},R}([(Tuple(edges[i, :]), i) for i = 1:size(edges, 1)])

        # e12 -> tuple -> get
        e12_tup = [Tuple(e12[i, :]) for i = 1:size(e12, 1)]
        e23_tup = [Tuple(e23[i, :]) for i = 1:size(e23, 1)]
        e31_tup = [Tuple(e31[i, :]) for i = 1:size(e31, 1)]
        faces_to_edges_tuple = cat(e23_tup, e31_tup, e12_tup; dims = 2)

        faces_to_edges = map(x -> get(edges_to_key, x, -1), faces_to_edges_tuple)

        m._edges_packed = edges
        m._edges_to_key = edges_to_key
        m._faces_to_edges_packed = faces_to_edges
        return nothing
    end
end

function _compute_laplacian_packed(m::TriMesh{T,R}, refresh::Bool = false) where {T,R}
    if refresh || (m._laplacian_packed isa Nothing)
        # laplacian sparse `L`:
        # where i, j are vertices indices, then
        # L[i, j] =    -1       , if i == j
        # L[i, j] = 1 / deg(i)  , if (i, j) is an edge
        # L[i, j] =    0        , otherwise
        verts = get_verts_packed(m)
        edges = get_edges_packed(m)

        # computing all i, j of edges (i,j), as (j,i) is also an edge
        e1 = edges[:, 1]
        e2 = edges[:, 2]

        # computing all edges, considering (i,j) and (j,i)
        idx12 = cat(e1, e2, dims = 2)
        idx21 = cat(e2, e1, dims = 2)
        idx = cat(idx12, idx21, dims = 1)

        A = sparse(
            idx[:, 1],
            idx[:, 2],
            fill(1, size(idx, 1)),
            size(verts, 2),
            size(verts, 2),
        )

        # computing degree of each vertices
        deg = sum(A, dims = 2)

        # if degree >= 1, then there exist edge with that vertex,
        # so computing 1/deg[i]
        deg1 = map(x -> T.(x > 0 ? 1 / x : x), deg[e1])
        deg2 = map(x -> T.(x > 0 ? 1 / x : x), deg[e2])
        diag = fill(T.(-1.0), size(verts, 2))

        # assuming (i,j) is an edge, so we computing 1/deg[i],
        # and store at L[i,j] and same for (j,i).
        # We also store -1 for L[i,i] for all vertices.
        Is = cat(e1, e2, R.(1:size(verts, 2)); dims = 1)
        Js = cat(e2, e1, R.(1:size(verts, 2)); dims = 1)
        Vs = cat(deg1, deg2, diag; dims = 1)
        m._laplacian_packed = sparse(Is, Js, Vs, size(verts, 2), size(verts, 2))
        return nothing
    end
end