text:
#!/usr/bin/env julia
using ObjectDetectionStats
using Test

@testset "Basic Data Ops" begin
    z = DataFrame( (:a => [1,2], :B => [120,13], :c => [111,111]))
    @test highest_vote!(z, [:a,:B,:c])[:Highest_Vote] .== [ :B, :c ]

    barycenter = DataFrame( Dict( [   :notneeded    => ["a", "b"],
                                      :center_x     => [15, 15],
                                      :center_y     => [10, 10],
                                      :width        => [5, 10],
                                      :height       => [5, 10]
                                  ]) )

    @test has_required_columns(barycenter, [ :center_x, :center_y, :width, :height ]  )
    @test !has_required_columns(barycenter, [ :center_x, :center_y, :width, :height, :NOPE ]  )
    #Test DF for comparing all others with
    KeyDF = DataFrame( Dict( [    :notneeded        => ["a", "b"],
                                  :upper_left_x     => [12.5, 10.0],
                                  :lower_right_x    => [17.5, 20.0],
                                  :upper_left_y     => [7.50, 5.00],
                                  :lower_right_y    => [12.5, 15.0]
                             ] ) )
    #Convert barycenter boxes to ULLR
    barycenter_boxes!( barycenter )
    @test all( Matrix(KeyDF) .== Matrix(barycenter[!, names(KeyDF)]) )

    area!(barycenter)
    @test all( barycenter[ :Area ] .== [ 25, 100 ] )

end

@testset "Box Translate and Clamp Tests" begin
    a            = Box(5,5,15,15)
    trans_result = translate( a, -5, -5 )
    @test trans_result                          == Box(0,0,10,10)
    @test ObjectDetectionStats.clamp(trans_result, Box(5,5,15,15))   == Box(5,5,10,10)
    @test ObjectDetectionStats.clamp(trans_result, Box(0,0, 5, 5))   == Box(0,0, 5, 5)
end

@testset "Box Area Tests" begin
    @test area( Box(5,5,15,15) ) == 100
    @test area( Box(5,5,11,11) ) == 36
end

@testset "Box Intersection Tests" begin
    a = Box(5,5,15,15)
    b = Box(5,5,15,15)
    @test intersection_area( a, b ) == 100
    @test intersection_over_union( a, b ) == 1.0

    a = Box(10,10,15,15)
    b = Box(5,5,15,15)
    @test intersection_area( a, b ) == 25
    @test intersection_over_union( a, b ) == 0.25

    b = Box(10,10,15,15)
    a = Box(5,5,15,15)
    @test intersection_area( a, b ) == 25
    @test intersection_over_union( a, b ) == 0.25

    a = Box(1,1,15,15)
    b = Box(25,25,35,35)
    @test intersection_area( a, b ) == 0
    @test intersection_over_union( a, b ) == 0.0
end

@testset "3 TP test" begin
    ods_machine     = ObjectDetectionScore( 3 )
    pred_scores     = [ 0.2 0.9 0.5 ; #3
                        0.2 0.9 0.5;  #2
                        0.5 0.9 0.2]  #1
    pred_locations  = [ Box( 1,     1,  10,     10 ),
                        Box( 15,    1,  25,     20 ),
                        Box( 1,    15,  10,     25 ) ]
    GT_cold_encodings = [ 2, 2, 2 ]# all correct
    #prepare inputs for evaluation...
    hcl = HotClassLocalization( pred_scores, pred_locations )
    ccl = ColdClassLocalization( GT_cold_encodings, pred_locations )
    ods_machine( hcl, ccl )
    @test all(ods_machine.TP .== [0,3,0])
    @test all(ods_machine.FP .== [0,0,0])
    @test all(ods_machine.FN .== [0,0,0])
end

@testset "3 wrong class test" begin
    ods_machine     = ObjectDetectionScore( 3 )
    pred_scores     = [ 0.2 0.9 0.5 ; #3
                        0.2 0.9 0.5; #2
                        0.5 0.9 0.2]  #1
    pred_locations  = [ Box( 1,     1,  10,     10 ),
                        Box( 15,    1,  25,     20 ),
                        Box( 1,    15,  10,     25 ) ]
    GT_cold_encodings = [ 1, 1, 1 ]# all incorrect class but perfect overlap
    #prepare inputs for evaluation...
    hcl = HotClassLocalization( pred_scores, pred_locations )
    ccl = ColdClassLocalization( GT_cold_encodings, pred_locations )

    ods_machine( hcl, ccl )
    @test all(ods_machine.TP .== [0,0,0])
    @test all(ods_machine.FP .== [0,3,0])
    @test all(ods_machine.FN .== [3,0,0])
end

@testset "3 wrong class locations with correct predictions" begin
    ods_machine     = ObjectDetectionScore( 3 )
    pred_scores     = [ 0.2 0.9 0.5 ; #2
                        0.2 0.9 0.5;  #2
                        0.5 0.9 0.2 ] #2
    pred_locations  = [ Box( 1,     1,  10,     10 ),
                        Box( 15,    1,  25,     20 ),
                        Box( 1,    15,  10,     25 ) ]
    GT_cold_encodings = [ 2, 2, 2 ]# all correct class
    GT_locations  = [   Box( 100,    100,  100,     100 ),
                        Box( 150,    100,  250,     200 ),
                        Box( 100,    150,  100,     250 ) ]
    #prepare inputs for evaluation...
    hcl = HotClassLocalization( pred_scores, pred_locations )
    ccl = ColdClassLocalization( GT_cold_encodings, GT_locations )

    ods_machine( hcl, ccl )

    @test all(ods_machine.TP .== [0,0,0])
    @test all(ods_machine.FP .== [0,3,0])
    @test all(ods_machine.FN .== [0,3,0])
end

@testset "3 wrong class locations with wrong predictions" begin
    ods_machine     = ObjectDetectionScore( 3 )
    pred_scores     = [ 0.2 0.5 0.9; #3
                        0.2 0.9 0.5; #2
                        0.9 0.5 0.2] #1
    pred_locations  = [ Box( 1,     1,  10,     10 ),
                        Box( 15,    1,  25,     20 ),
                        Box( 1,    15,  10,     25 ) ]
    GT_cold_encodings     = [   3, #correct
                                1, #incorrect
                                1 ]#correct
    GT_locations  = [   Box( 1,     1,  10,     10 ), #TP class 3
                        Box( 15,    1,  25,     20 ), #FP class 1
                        Box( 30,    30, 42,     42 )] #FN class 1
    #prepare inputs for evaluation...
    hcl = HotClassLocalization( pred_scores, pred_locations )
    ccl = ColdClassLocalization( GT_cold_encodings, GT_locations )
    ods_machine( hcl, ccl )
    @test all(ods_machine.TP .== [0,0,1])
    @test all(ods_machine.FP .== [1,1,0])
    @test all(ods_machine.FN .== [2,0,0])
end

@testset "2 correct overlapping a prediction and 1 nonprediction" begin
    ods_machine     = ObjectDetectionScore( 3 )
    pred_scores     = [ 0.2 0.9 0.5 ; #2
                        0.2 0.9 0.5;  #2
                        0.0 0.0 0.0;] #0 not predicted
    pred_locations  = [ Box( 1,     1,  10,     10 ),
                        Box( 1,     1,  10,     10 ),
                        Box( 15,   15,  20,     20 )  ]
    GT_cold_encodings = [ 2 ]# all correct class
    GT_locations  = [   Box( 1,     1,  10,     10 ) ]
    #prepare inputs for evaluation...
    hcl = HotClassLocalization( pred_scores, pred_locations )
    ccl = ColdClassLocalization( GT_cold_encodings, GT_locations )

    ods_machine( hcl, ccl )
    @test all(ods_machine.TP .== [0,1,0])
    @test all(ods_machine.FP .== [0,1,0])
    @test all(ods_machine.FN .== [0,0,0])
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

@doc raw"""
`struct BasicPSH1pBasis <: OneParticleBasis`

One-particle basis of the form
```math
\phi_{nlm}({\bm r}, z_1, z_0) = J_{n}(r) Y_l^m(\hat{\br r})
```
where ``J_{n}`` denotes a radial basis.
"""
struct BasicPSH1pBasis{T, NZ, TJ <: ScalarBasis{T}} <: OneParticleBasis{T}
   J::TJ
   SH::SHBasis{T}
   zlist::SZList{NZ}
   spec::Vector{PSH1pBasisFcn}
   Aindices::Matrix{UnitRange{Int}}
end

function BasicPSH1pBasis(J::ScalarBasis{T}, zlist::SZList,
                         spec::Vector{PSH1pBasisFcn}) where {T}
   # now get the maximum L-degree to generate the SH basis
   maxL = maximum(b.l for b in spec)
   SH = SHBasis(maxL, T)
   # construct a preliminary Aindices array to get an incorrect basis
   NZ = length(zlist)
   P = BasicPSH1pBasis(J, SH, zlist, spec,
                       Matrix{UnitRange{Int}}(undef, NZ, NZ))
   # ... now fix the Aindices array
   set_Aindices!(P)
   return P
end

function BasicPSH1pBasis(J::ScalarBasis;
                         species = :X,
                         D::AbstractDegree = SparsePSHDegree())
   # get a generic basis spec
   spec = _build_PSH_1p_spec(length(J), D, species)
   # construct the basis
   zlist = ZList(species; static=true)
   return BasicPSH1pBasis(J, zlist, spec)
end

cutoff(basis::BasicPSH1pBasis) = cutoff(basis.J)

Base.length(basis::BasicPSH1pBasis, z0::AtomicNumber) =
      numz(basis) *  length(basis.spec)

Base.length(basis::BasicPSH1pBasis, iz0::Integer) =
      numz(basis) *  length(basis.spec)

Base.length(basis::BasicPSH1pBasis, iz::Integer, iz0::Integer) =
      length(basis.spec)

function get_basis_spec(basis::BasicPSH1pBasis, z0::AtomicNumber)
   iz0 = z2i(basis, z0)
   len_iz0 = sum(length(basis, iz, iz0) for iz = 1:numz(basis))
   spec = Vector{PSH1pBasisFcn}(undef, len_iz0)
   for iz = 1:numz(basis)
      z = i2z(basis, iz)
      spec[basis.Aindices[iz, iz0]] =
               [ PSH1pBasisFcn(b.n, b.l, b.m, z)   for b in basis.spec ]
   end
   return spec
end

get_basis_spec(basis::BasicPSH1pBasis, iz0::Integer, i::Integer) =
   get_basis_spec(basis, i2z(basis, iz0), i)

function get_basis_spec(basis::BasicPSH1pBasis, z0::AtomicNumber, i::Integer)
   iz = ((i-1) ÷ length(basis.spec)) + 1
   inew = mod1(i, length(basis.spec))
   b = basis.spec[inew]
   return PSH1pBasisFcn(b.n, b.l, b.m, i2z(basis, iz))
end

_build_specnl(maxdeg, D, z, z0) =
         gensparse(2, maxdeg;               #     n     l    m  z
                    tup2b = t -> PSH1pBasisFcn(t[1]+1, t[2], 0, z),
                    degfun = b -> D(b, z0),
                    ordered = false)

@doc raw"""
`function _build_PSH_1p_spec`

Construct the specification for a ``P \otimes Y`` type 1-particle basis.
These must be treated differently because of the requirements that complete
``l``-blocks are represented in the basis.

See also: `_get_1p_spec`.
"""
function _build_PSH_1p_spec(maxn::Integer, D::AbstractDegree, species)
   # find out what the largest degree is that we can allow:
   specnl = []
   for s1 in species, s2 in species
      z, z0 = AtomicNumber(s1), AtomicNumber(s2)
      maxdeg = maximum( D(PSH1pBasisFcn(n, 0, 0, z), z0) for n = 1:maxn )
      # generate the `spec::Vector{PSH1pBasisFcn}` using maxn
      specnl_zz0 = _build_specnl(maxdeg, D, z, z0)
      append!(specnl, specnl_zz0)
   end
   specnl = unique([ PSH1pBasisFcn(b.n, b.l, 0, 0) for b in specnl ])
   # add the m-parameters
   return [ PSH1pBasisFcn(b.n, b.l, m, 0)
              for b in specnl for m = -b.l:b.l ]
end

_build_PSH_1p_spec(maxn::Integer, D::AbstractDegree, species::Symbol) =
         _build_PSH_1p_spec(maxn::Integer, D, (species,))

# ------------------------------------------------------
#  FIO code

==(P1::BasicPSH1pBasis, P2::BasicPSH1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::BasicPSH1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_BasicPSH1pBasis",
           "J" => write_dict(basis.J),
          "SH" => write_dict(basis.SH),
        "spec" => write_dict.(basis.spec),
       "zlist" => write_dict(basis.zlist),
   )

read_dict(::Val{:SHIPs_BasicPSH1pBasis}, D::Dict) =
   read_dict(Val{:ACE_BasicPSH1pBasis}(), D::Dict)

function read_dict(::Val{:ACE_BasicPSH1pBasis}, D::Dict)
   J = read_dict(D["J"])
   SH = read_dict(D["SH"])
   zlist = read_dict(D["zlist"])
   spec = read_dict.(D["spec"])
   P = BasicPSH1pBasis(J, SH, zlist, spec,
                 Matrix{UnitRange{Int}}(undef, length(zlist), length(zlist)))
   set_Aindices!(P)
   return P
end

# ------------------------------------------------------
#  Evaluation code

fltype(basis::BasicPSH1pBasis{T}) where T = Complex{T}
rfltype(basis::BasicPSH1pBasis{T}) where T = T

alloc_temp(basis::BasicPSH1pBasis, args...) =
      (
        BJ = alloc_B(basis.J, args...),
        tmpJ = alloc_temp(basis.J, args...),
        BY = alloc_B(basis.SH, args...),
        tmpY = alloc_temp(basis.SH, args...),
       )

function add_into_A!(A, tmp, basis::BasicPSH1pBasis,
                     R, iz::Integer, iz0::Integer)
   # evaluate the r-basis and the R̂-basis for the current neighbour at R
   evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(R))
   evaluate!(tmp.BY, tmp.tmpY, basis.SH, R)
   # add the contributions to the A_zklm
   @inbounds for (i, nlm) in enumerate(basis.spec)
      A[i] += tmp.BJ[nlm.n] * tmp.BY[index_y(nlm.l, nlm.m)]
   end
   return nothing
end

# function add_into_A!(A, inds, tmp, basis::BasicPSH1pBasis,
#                      R, iz::Integer, iz0::Integer)
#    # evaluate the r-basis and the R̂-basis for the current neighbour at R
#    evaluate!(tmp.BJ, tmp.tmpJ, basis.J, norm(R))
#    evaluate!(tmp.BY, tmp.tmpY, basis.SH, R)
#    # add the contributions to the A_zklm
#    @inbounds for (i, nlm) in zip(inds, basis.spec)
#       A[i] += tmp.BJ[nlm.n] * tmp.BY[index_y(nlm.l, nlm.m)]
#    end
#    return nothing
# end

alloc_temp_d(basis::BasicPSH1pBasis, args...) =
      (
        BJ = alloc_B(basis.J, args...),
        tmpJ = alloc_temp(basis.J, args...),
        BY = alloc_B(basis.SH, args...),
        tmpY = alloc_temp(basis.SH, args...),
        dBJ = alloc_dB(basis.J, args...),
        tmpdJ = alloc_temp_d(basis.J, args...),
        dBY = alloc_dB(basis.SH, args...),
        tmpdY = alloc_temp_d(basis.SH, args...),
       )

function add_into_A_dA!(A, dA, tmpd, basis::BasicPSH1pBasis, R, iz::Integer, iz0::Integer)
   r = norm(R)
   R̂ = R / r
   # evaluate the r-basis and the R̂-basis for the current neighbour at R
   evaluate_d!(tmpd.BJ, tmpd.dBJ, tmpd.tmpdJ, basis.J, r)
   evaluate_d!(tmpd.BY, tmpd.dBY, tmpd.tmpdY, basis.SH, R)
   # add the contributions to the A_zklm, ∇A
   @inbounds for (i, nlm) in enumerate(basis.spec)
      iY = index_y(nlm.l, nlm.m)
      A[i] += tmpd.BJ[nlm.n] * tmpd.BY[iY]
      dA[i] = (tmpd.dBJ[nlm.n] * tmpd.BY[iY]) * R̂ + tmpd.BJ[nlm.n] * tmpd.dBY[iY]
   end
   return nothing
end

ext: jl


text:
abstract type AbstractWatson <: AbstractFloat end

for DrWatsonN in (:DrWatson16, :DrWatson32, :DrWatson64)

@eval FloatN = $(Symbol("Float",string(DrWatsonN)[end-1:end]))

@eval begin

struct $DrWatsonN{f} <: AbstractWatson
    val::$FloatN
end

# conversions back from DrWatson64
Base.Float64(x::$DrWatsonN) = Float64(x.val)
Base.Float32(x::$DrWatsonN) = Float32(x.val)
Base.Float16(x::$DrWatsonN) = Float16(x.val)
Base.Int64(x::$DrWatsonN) = Int64(x.val)
Base.Int32(x::$DrWatsonN) = Int32(x.val)
Base.Int16(x::$DrWatsonN) = Int16(x.val)

# conversion among DrWatsons
$DrWatsonN{f}(x::$DrWatsonN) where f = $DrWatsonN{f}(x.val)
$DrWatsonN(x::$DrWatsonN{f}) where {f} = $DrWatsonN{f}(x.val)

Base.promote_rule(::Type{<:Integer},::Type{$DrWatsonN{f}}) where f = $DrWatsonN{f}
Base.promote_rule(::Type{Float64},::Type{$DrWatsonN{f}}) where f = $DrWatsonN{f}
Base.promote_rule(::Type{Float32},::Type{$DrWatsonN{f}}) where f = $DrWatsonN{f}
Base.promote_rule(::Type{Float16},::Type{$DrWatsonN{f}}) where f = $DrWatsonN{f}

Base.bitstring(x::$DrWatsonN) = bitstring(x.val)

Base.eps(::Type{$DrWatsonN{f}}) where f = eps(Float64)
Base.eps(x::$DrWatsonN) = eps(x.val)

Base.typemin(::$DrWatsonN{f}) where f = $DrWatsonN{f}(typemin($FloatN))
Base.typemax(::$DrWatsonN{f}) where f = $DrWatsonN{f}(typemax($FloatN))
Base.floatmin(::$DrWatsonN{f}) where f = $DrWatsonN{f}(floatmin($FloatN))
Base.floatmax(::$DrWatsonN{f}) where f = $DrWatsonN{f}(floatmax($FloatN))
Base.precision(::$DrWatsonN{f}) where f = $DrWatsonN{f}(precision($FloatN))

end

for O in (:(+), :(-), :(*), :(/), :(^))
    @eval function Base.$O(x::$DrWatsonN{f},y::$DrWatsonN{f}) where f
        r = $O(x.val, y.val)
        scent(f,r)
        return $DrWatsonN{f}(r)
    end
end

for O in (:(-), :(+),
          :sign,
          :prevfloat, :nextfloat,
          :round, :trunc, :ceil, :floor,
          :inv, :abs, :sqrt, :cbrt,
          :exp, :expm1, :exp2, :exp10,
          :log, :log1p, :log2, :log10,
          :rad2deg, :deg2rad, :mod2pi, :rem2pi,
          :sin, :cos, :tan, :csc, :sec, :cot,
          :asin, :acos, :atan, :acsc, :asec, :acot,
          :sinh, :cosh, :tanh, :csch, :sech, :coth,
          :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
          :sinc, :sinpi, :cospi,
          :sind, :cosd, :tand, :cscd, :secd, :cotd,
          :asind, :acosd, :atand, :acscd, :asecd, :acotd
         )
    @eval begin
        function Base.$O(x::$DrWatsonN{f}) where f
            r = $O(x.val)
            scent(f,r)
            return $DrWatsonN{f}(r)
        end
    end
end

for O in (:(<), :(<=))     # do not trigger scent for those
    @eval Base.$O(x::$DrWatsonN, y::$DrWatsonN) = $O(x.val, y.val)
end

@eval Random.rand(rng::Random.AbstractRNG, ::Random.Sampler{T}) where T<:$DrWatsonN = convert(T, rand(rng, $FloatN))
@eval Random.randn(rng::Random.AbstractRNG, ::Type{T}) where T<:$DrWatsonN = convert(T, randn(rng, $FloatN))
@eval Random.randexp(rng::Random.AbstractRNG, ::Type{T}) where T<:$DrWatsonN = convert(T, randexp(rng, $FloatN))

@eval round(::Type{T}, x::$DrWatsonN) where {T<:Integer} = round(T, x.val)
@eval rem(x::$DrWatsonN, y::$DrWatsonN) = Core.Intrinsics.rem_float(x.val, y.val)

end


text:
# Convenient dictionary for mapping powers of ten to an SI prefix.
const prefixdict = Dict(
    -24 => "y",
    -21 => "z",
    -18 => "a",
    -15 => "f",
    -12 => "p",
    -9  => "n",
    -6  => "μ",     # tab-complete \mu, not option-m on a Mac!
    -3  => "m",
    -2  => "c",
    -1  => "d",
    0   => "",
    1   => "da",
    2   => "h",
    3   => "k",
    6   => "M",
    9   => "G",
    12  => "T",
    15  => "P",
    18  => "E",
    21  => "Z",
    24  => "Y"
)

"""
`abbr(x)` provides abbreviations for units or dimensions. Since a method should
always be defined for each unit and dimension type, absence of a method for a
specific unit or dimension type is likely an error. Consequently, we return ❓
for generic arguments to flag unexpected behavior.
"""
abbr(x) = "❓"     # Indicate missing abbreviations

"""
```
prefix(x::Unit)
```

Returns a string representing the SI prefix for the power-of-ten held by
this particular unit.
"""
function prefix(x::Unit)
    if haskey(prefixdict, tens(x))
        return prefixdict[tens(x)]
    else
        error("Invalid power-of-ten prefix.")
    end
end

"""
```
show(io::IO, x::Quantity)
```

Show a unitful quantity by calling `show` on the numeric value, appending a
space, and then calling `show` on a units object `U()`.
"""
function show(io::IO, x::Quantity)
    show(io,x.val)
    if unit(x) != NoUnits
        print(io," ")
        show(io, unit(x))
    end
    nothing
end

"""
```
show{T,D,U}(io::IO, ::Type{Quantity{T,D,U}})
```

Show the type of a unitful quantity in a succinct way. Otherwise,
array summaries are nearly unreadable.
"""
function show{T,D,U}(io::IO, ::Type{Quantity{T,D,U}})
    print(io, "Quantity{", string(T),
            ", Dimensions:{", string(D()),
            "}, Units:{", string(U()), "}}")
    nothing
end

"""
```
show(io::IO,x::Unitlike)
```

Call [`Unitful.showrep`](@ref) on each object in the tuple that is the type
variable of a [`Unitful.Units`](@ref) or [`Unitful.Dimensions`](@ref) object.
"""
function show(io::IO,x::Unitlike)
    first = ""
    tup = typeof(x).parameters[1]
    map(tup) do y
        print(io,first)
        showrep(io,y)
        first = " "
    end
    nothing
end

"""
```
showrep(io::IO, x::Unit)
```

Show the unit, prefixing with any decimal prefix and appending the exponent as
formatted by [`Unitful.superscript`](@ref).
"""
function showrep(io::IO, x::Unit)
    print(io, prefix(x))
    print(io, abbr(x))
    print(io, (power(x) == 1//1 ? "" : superscript(power(x))))
    nothing
end

"""
```
showrep(io::IO, x::Dimension)
```

Show the dimension, appending any exponent as formatted by
[`Unitful.superscript`](@ref).
"""
function showrep(io::IO, x::Dimension)
    print(io, abbr(x))
    print(io, (power(x) == 1//1 ? "" : superscript(power(x))))
end

"""
```
superscript(i::Rational)
```

Prints exponents.
"""
function superscript(i::Rational)
    i.den == 1 ? "^"*string(i.num) : "^"*replace(string(i),"//","/")
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractDomain{T,N}

Spatial domain in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type AbstractDomain{T,N} <: AbstractSpatialObject{T,N} end

"""
    view(domain, locations)

Return a view of `domain` with all points in `locations`.
"""
Base.view(domain::AbstractDomain,
          locations::AbstractVector{Int}) = DomainView(domain, locations)

# ------------
# IO methods
# ------------
function Base.show(io::IO, domain::AbstractDomain{T,N}) where {N,T}
  npts = npoints(domain)
  print(io, "$npts SpatialDomain{$T,$N}")
end

#------------------
# IMPLEMENTATIONS
#------------------
include("domains/curve.jl")
include("domains/point_set.jl")
include("domains/regular_grid.jl")
include("domains/structured_grid.jl")


text:
"Parses a regex into a string"
function regtostr(r::Regex)
  regx = string(r)
  regx = regx[3:length(regx)-1]
  regx = strip(regx,'\$')
  return regx
end

oprp = Regex("[+|-]")
oprm = Regex("[/|*]")
oprln = Regex("not")
type_str = Regex("(\").*(\")\$")
type_int = Regex(regtostr(oprp)*"?[0-9]+\$")
type_float = Regex("("*regtostr(type_int)*")[\".\"']{1}[0-9]*\$")
id = Regex("^[A-Z|a-z]+(\w)*")

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
module TestClustering

# using Revise
using MLJBase
using Test
using Random:seed!
import LinearAlgebra: norm
import Distances: evaluate

# load code to be tested:
import MLJModels
import Clustering # MLJModels.Clustering_ now available for loading
using MLJModels.Clustering_

seed!(132442)

task = load_crabs()

X, y = X_and_y(task)

####
#### KMEANS
####

barekm = KMeans()

fitresult, cache, report = MLJBase.fit(barekm, 1, X)

R = MLJBase.matrix(MLJBase.transform(barekm, fitresult, X))

X_array = MLJBase.matrix(X)

# distance from first point to second center
@test R[1, 2] ≈ norm(view(X_array, 1, :) .- view(fitresult, :, 2))^2
@test R[10, 3] ≈ norm(view(X_array, 10, :) .- view(fitresult, :, 3))^2

p = MLJBase.predict(barekm, fitresult, X)

@test argmin(R[1, :]) == p[1]
@test argmin(R[10, :]) == p[10]

# km = machine(barekm, X)
# fit!(km)

info(barekm)

####
#### KMEDOIDS
####

barekm = KMedoids()

fitresult, cache, report = MLJBase.fit(barekm, 1, X)

R = MLJBase.matrix(MLJBase.transform(barekm, fitresult, X))

@test R[1, 2] ≈ evaluate(barekm.metric, view(X_array, 1, :), view(fitresult, :, 2))
@test R[10, 3] ≈ evaluate(barekm.metric, view(X_array, 10, :), view(fitresult, :, 3))

p = MLJBase.predict(barekm, fitresult, X)

@test all(report.assignments .== p)

# km = machine(barekm, X)
# fit!(km)

info(barekm)

end # module
true


text:
"""
    simplecycles_limited_length(g, n, ceiling=10^6)

Compute and return at most `ceiling` cycles of length at most `n` of
the given graph. Both directed and undirected graphs are supported.

### Performance
The number of cycles grows very fast with the number of vertices and
the allowed length of the cycles. This function is intended for
finding short cycles. If you want to find cycles of any length in a
directed graph, `simplecycles` or `simplecycles_iter` may be more
efficient.
"""
function simplecycles_limited_length(graph::AbstractGraph{T}, n::Int,
                                     ceiling = 10^6) where {T}
    cycles = Vector{Vector{T}}()
    n < 1 && return cycles
    cycle = Vector{T}(undef, n)
    @inbounds for v in vertices(graph)
        cycle[1] = v
        simplecycles_limited_length!(graph, n, ceiling, cycles, cycle, 1)
        length(cycles) >= ceiling && break
    end
    return cycles
end

function simplecycles_limited_length!(graph, n, ceiling, cycles, cycle, i)
    length(cycles) >= ceiling && return
    for v in outneighbors(graph, cycle[i])
        if v == cycle[1]
            push!(cycles, cycle[1:i])
        elseif (i < n
                && v > cycle[1]
                && !repeated_vertex(v, cycle, 2, i))
            cycle[i + 1] = v
            simplecycles_limited_length!(graph, n, ceiling, cycles, cycle, i + 1)
        end
    end
end

# Note: Doing the same thing as `@views v in cycle[n1:n2]`, but until
# views are completely allocation free this can be expected to be
# faster.
function repeated_vertex(v, cycle, n1, n2)
    for k = n1:n2
        cycle[k] == v && return true
    end
    return false
end


text:
"""
    band(x, ylower, yupper; kwargs...)
    band(lower, upper; kwargs...)

Plots a band from `ylower` to `yupper` along `x`. The form `band(lower, upper)` plots a [ruled surface](https://en.wikipedia.org/wiki/Ruled_surface)
between the points in `lower` and `upper`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Band, lowerpoints, upperpoints) do scene
    attr = Attributes(;
        default_theme(scene, Mesh)...,
        colorrange = automatic,
    )
    attr[:shading][] = false
    attr
end

convert_arguments(::Type{<: Band}, x, ylower, yupper) = (Point2f.(x, ylower), Point2f.(x, yupper))

function band_connect(n)
    ns = 1:n-1
    ns2 = n+1:2n-1
    [GLTriangleFace.(ns, ns .+ 1, ns2); GLTriangleFace.(ns .+ 1, ns2 .+ 1, ns2)]
end

function plot!(plot::Band)
    @extract plot (lowerpoints, upperpoints)
    @lift(@assert length($lowerpoints) == length($upperpoints) "length of lower band is not equal to length of upper band!")
    coordinates = @lift([$lowerpoints; $upperpoints])
    connectivity = lift(x -> band_connect(length(x)), plot[1])

    meshcolor = lift(Any, plot.color) do c
        if c isa AbstractArray
            # if the same number of colors is given as there are
            # points on one side of the band, the colors are mirrored to the other
            # side to make an even band
            if length(c) == length(lowerpoints[])
                return repeat(c, 2)
            # if there's one color for each band vertex, the colors are used directly
            elseif length(c) == 2 * length(lowerpoints[])
                return c
            else
                error("Wrong number of colors. Must be $(length(lowerpoints[])) or double.")
            end
        else
            return c
        end
    end

    attr = Attributes(plot)
    attr[:color] = meshcolor
    mesh!(plot, attr, coordinates, connectivity)
end

function fill_view(x, y1, y2, where::Nothing)
    x, y1, y2
end
function fill_view(x, y1, y2, where::Function)
    fill_view(x, y1, y2, where.(x, y1, y2))
end
function fill_view(x, y1, y2, bools::AbstractVector{<: Union{Integer, Bool}})
    view(x, bools), view(y1, bools), view(y2, bools)
end

"""
    fill_between!(scenelike, x, y1, y2; where = nothing, kw_args...)

fill the section between 2 lines with the condition `where`
"""
function fill_between!(scenelike, x, y1, y2; where = nothing, kw_args...)
    xv, ylow, yhigh = fill_view(x, y1, y2, where)
    band!(scenelike, xv, ylow, yhigh; kw_args...)
end

export fill_between!


text:
#cd(homedir()*"/GitHubRepos/Working-directory/HyperModularity")
# using Revise
# using Pkg; Pkg.activate(".")
using HyperModularity
using Test
using StatsBase
using Combinatorics
using SparseArrays

include("../src/test_funs.jl")

@testset "utils" begin
    include("utils_tests.jl")
end

@testset "Ω" begin
    include("omega_tests.jl")
end

@testset "hsbm" begin
    n = 20
    Z = rand(1:5, n)
    Z = Int64.(Z)
    ϑ = dropdims(ones(1,n) + rand(1,n), dims = 1)
    μ = mean(ϑ)
    kmax = Int64(4)

    ω(x, α) = (x[1]+1)^(-α[x[2]]*x[2])

    α = repeat([2.0], kmax)

    Ω = sumOfExteriorDegreesIntensityFunction(ω, kmax)
    H = sampleSBM(Z, ϑ, Ω;α=α, kmax=kmax, kmin = 1)

    @test length(keys(H.E)) == kmax
end

@testset "louvain_utils" begin

    include("louvain_utils_tests.jl")
end

# @testset "vol" begin
#     # some tests here are erroring
#     # include("vol_tests.jl")
# end

# @testset "inference" begin
#
# end
#
# @testset "analysis_helpers" begin
#
# end
#
# @testset "louvain_utils" begin
#
# end
#
# @testset "graph_louvain" begin
#
# end
#
# @testset "general_louvain" begin
#
# end
#
# @testset "aon_louvain" begin
#
# end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_translated)" --startup-file=no %s 2>&1 | FileCheck %s

# has the terminator unreachable
hasunreachable(x::Float64) = sqrt(x)
emit(hasunreachable, Float64)

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.hasunreachable), Float64}, svec(), hasunreachable(x::Float64) in Main.Main at /{{.*}}/test/Codegen/translate/hasunreachable.jl:4, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.hasunreachable), Float64}"(%arg0: !jlir<"typeof(Main.hasunreachable)">, %arg1: !jlir.Float64) -> !jlir.Float64 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"#<intrinsic #33 lt_float>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Float64
# CHECK-NEXT:     %2 = "jlir.call"(%0, %arg1, %1) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Float64, !jlir.Float64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%2)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir.Base.Math.throw_complex_domainerror} : () -> !jlir<"typeof(Base.Math.throw_complex_domainerror)">
# CHECK-NEXT:     %4 = "jlir.constant"() {value = #jlir<":sqrt">} : () -> !jlir.Symbol
# CHECK-NEXT:     %5 = "jlir.invoke"(%3, %4, %arg1) {methodInstance = #jlir<"throw_complex_domainerror(Symbol, Float64) from throw_complex_domainerror(Symbol, Any)">} : (!jlir<"typeof(Base.Math.throw_complex_domainerror)">, !jlir.Symbol, !jlir.Float64) -> !jlir<"Union{}">
# CHECK-NEXT:     %6 = "jlir.undef"() : () -> !jlir.Float64
# CHECK-NEXT:     "jlir.return"(%6) : (!jlir.Float64) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %7 = "jlir.constant"() {value = #jlir<"#<intrinsic #78 sqrt_llvm>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %8 = "jlir.call"(%7, %arg1) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     "jlir.goto"()[^bb4] : () -> ()
# CHECK-NEXT:   ^bb4:  // pred: ^bb3
# CHECK-NEXT:     "jlir.return"(%8) : (!jlir.Float64) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
"""
    crossTable(x, y; row=true)

Make cross-tabulation table between arrays `x` and `y`.

Optional `row` report row percentage where `row=false` report column percentage, `digits` specifies the digits use in `round` function.

Return table of complete cases only, i.e., missing values are removed. This table can be used as the input to `chisquaretest`, `oddsratio`, and `riskratio` functions.

## Examples
```julia
julia> using Epi

julia> gender = ["Female", "Male"][rand(1:2, 100)];

julia> age = ["0-5", "6-10", "11-15"][rand(1:3, 100)];

julia> tm_ = crossTable(age, gender);
Frequency (row percentage)
4×3 Named Array{Text{String},2}
A ╲ B │      Male     Female      Total
──────┼────────────────────────────────
11-15 │ 16 (0.53)  14 (0.47)   30 (1.0)
0-5   │ 15 (0.43)  20 (0.57)   35 (1.0)
6-10  │ 18 (0.51)  17 (0.49)   35 (1.0)
Total │ 49 (0.49)  51 (0.51)  100 (1.0)

julia> tm_ = crossTable(age, gender, row=false);
Frequency (column percentage)
4×3 Named Array{Text{String},2}
A ╲ B │      Male     Female      Total
──────┼────────────────────────────────
11-15 │ 16 (0.33)  14 (0.27)   30 (0.3)
0-5   │ 15 (0.31)  20 (0.39)  35 (0.35)
6-10  │ 18 (0.37)  17 (0.33)  35 (0.35)
Total │  49 (1.0)   51 (1.0)  100 (1.0)
```
"""
function crossTable(x::AbstractArray, y::AbstractArray; row=true, digits=2)
	if length(x)!=length(y) error("Different length") end
	data = removeMissing(hcat(x, y))
	xx, yy = (unique(data[:, 1]), unique(data[:, 2]))
	nr, nc = (length(xx), length(yy))
	if nc | nr > 7 print("Consider grouping the variables") end
	A = NamedArray{Int}(nr, nc)
	for i in 1:nr, j in 1:nc
		A[i, j] = crossCount(xx[i], yy[j], data)
	end
 	A = hcat(A, sum(A, dims=2))
 	A = vcat(A, sum(A, dims=1))
	setnames!(A, push!(string.(xx), string("Total")), 1)
	setnames!(A, push!(string.(yy), string("Total")), 2)
	if row
		pc = round.(A ./ A[:, end], digits=digits)
		println("Frequency (row percentage)")
		println(NamedArray(Text.(string.(Array(A), " (", pc, ")")),
		                   (names(A, 1), names(A, 2))))
	else
		println("Frequency (column percentage)")
		pr = round.(A ./ A[end, :]', digits=digits)
		println(NamedArray(Text.(string.(Array(A), " (", pr, ")")),
		                   (names(A, 1), names(A, 2))))
	end
	A # return Frequency table, no %
end


text:
using Qaintessent: measure
using Test
using TestSetExtensions
using LinearAlgebra
using Qaintessent
using StatsBase

##==----------------------------------------------------------------------------------------------------------------------
@testset ExtendedTestSet "measurement" begin
    @testset "Measure single instance" begin
        state_a = ComplexF64[1;0;0;0]
        state_b = ComplexF64[1;0;0]
        state_c = ComplexF64[0.5;0;0;0]
        @test measure(state_a) == Dict(0=>1)
        @test_throws ErrorException("Shots has to be a natural integer above 0.")  measure(state_a, 0)
        @test_throws ErrorException("Statevector is not a viable quantum state (length $(length(state_b))")  measure(state_b)
        @test_throws ErrorException("Statevector is not a viable statevector (norm is $(norm(state_c))")  measure(state_c)
    end

    @testset "Measure multiple instances" begin
        state_a = ComplexF64[1;0;0;0]
        circ = [circuit_gate(1, HadamardGate())]
        state_b = apply(state_a, circ)
        samples = measure(state_b, 20000)
        # Measurement result should follow binomial distribution, this should return false ~1 out of a trillion
        @test (samples[0]<10500 && samples[0]>9500)
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

#
# those functions realize the Msg protocol
# this is the 2nd actor layer
#

# return the behavior function
function bhvf(A::_ACT)
    f = A.bhv isa Bhv ? A.bhv.f : A.bhv
    return f isa Function ? f : typeof(f)
end

# Become
onmessage(A::_ACT, msg::Become) = A.bhv = _current(msg.x)
function onmessage(A::_ACT, msg::Call)
    A.res = A.bhv(msg.x...)
    send(msg.from, Response(A.res, A.self))
end

# Cast
onmessage(A::_ACT, msg::Cast) = A.res = A.bhv(msg.x...)

# Diag
function onmessage(A::_ACT, msg::Diag)
    x = first(msg.x)
    res = x == :act   ? A :
          x == :task  ?
            myid() == msg.from.pid ?
                current_task() :
                repr(current_task()) :
          x == :tid   ? pqtid() :
          x == :pid   ? myid() :
          x == :err   ? errored() :
          x == :info  ? Info(A.mode,bhvf(A),myid(),threadid(),tid(),pqtid(),A.name) :
          x == :state ? :ok : :unknown_request
    send(msg.from, Response(res, A.self))
end

# Connect
function onmessage(A::_ACT, msg::Connect)
    if !msg.remove
        push!(A.conn, msg.x)
        unique!(A.conn)
    else
        filter!(c->c.lk!=msg.x, A.conn)
    end
end

# Exec
onmessage(A::_ACT, msg::Exec) = send(msg.from, Response(_current(msg.func)(), A.self))

# Init
function onmessage(A::_ACT, msg::Init)
    A.init = _current(msg.x)
    if A.self.chn isa RemoteChannel
        foreach(A.conn) do c
            c isa Super && send(c.lk, ChildInit(A.self, A.init))
        end
    end
end

# Query
function onmessage(A::_ACT, msg::Query)
    msg.x in (:mode,:bhv,:res,:sta,:usr) ?
        send(msg.from, Response(getfield(A, msg.x), A.self)) :
        send(msg.from, Response("$(msg.x) not available", A.self))
end

# Update
function onmessage(A::_ACT, msg::Update)
    if msg.s in (:name,:self,:sta,:usr)
        setfield!(A, msg.s, msg.x)
    elseif msg.s == :mode
        A.mode = msg.x
        A.self.mode = msg.x
    elseif msg.s == :arg
        A.bhv = Bhv(A.bhv.f, msg.x.args...;
            pairs((; merge(A.bhv.kw, msg.x.kwargs)...))...)
    end
end

# Term
onmessage(A::_ACT, msg::Term) = A.term = _current(msg.x)

# dispatch on Request or user defined Msg
onmessage(A::_ACT, msg::Msg) = A.res = A.bhv(msg)


text:
# T is the type of the elements produced by this iterator
struct EnumerableFilter{T,S,Q<:Function} <: Enumerable
    source::S
    filter::Q
end

Base.eltype(iter::Type{EnumerableFilter{T,S,Q}}) where {T,S,Q} = T

function filter(source::Enumerable, filter_func::Function, filter_expr::Expr)
    T = eltype(source)
    S = typeof(source)
    Q = typeof(filter_func)
    return EnumerableFilter{T,S,Q}(source, filter_func)
end

function Base.iterate(iter::EnumerableFilter{T,S,Q}, state...) where {T,S,Q}
    ret = iterate(iter.source, state...)

    while ret!==nothing
        if iter.filter(ret[1])
            return ret
        end

        ret = iterate(iter.source, ret[2])
    end

    return nothing
end


text:
using OpenQuantumTools, OrdinaryDiffEq, Plots, Printf, LaTeXStrings

β = 4
T = β_2_temperature(β)
η = 0.1
fc= 10/(2π)
bath = Ohmic(η, fc, T)

plot(bath, :γ, range(0,10,length=100), linewidth=2, label="")

τsb, err_τsb = τ_SB((x)->correlation(x, bath))
@printf("τ_sb of the Ohmic bath is %.6f with error estimation %.2e \n", τsb, err_τsb)
τb, err_τb = τ_B((x)->correlation(x, bath), 100, τsb)
@printf("τ_b of the Ohmic bath is %.6f with error estimation %.2e \n", τb, err_τb)

Hp = 0.5*σz⊗σi - 0.7*σi⊗σz + 0.3*σz⊗σz
Hd = standard_driver(2)
H = DenseHamiltonian([(s)->1-s, (s)->s], [-Hd, Hp], unit=:ħ)

plot(H, range(0,1,length=100), 4, linewidth=2)
xlabel!("s")
ylabel!(L"P(s)")

tf = 20
ρ0 = (σi+σx)⊗(σi+σx)/4
coupling = ConstantCouplings([σz⊗σi, σi⊗σz], unit=:ħ)
annealing = Annealing(H, ρ0, bath=bath, coupling=coupling)
close_sol = solve_von_neumann(annealing, tf, alg = Tsit5(), abstol=1e-6, reltol=1e-6);

plot(close_sol, H, 1, range(0,tf,length=100), linewidth=2)
xlabel!("t")
ylabel!(L"P_G(s)")

t_axis = range(0,tf,length=100)
p_computational_basis = [real(diag(close_sol(s))) for s in t_axis]
p_computational_basis = hcat(p_computational_basis...)
plot(t_axis, p_computational_basis', linewidth=2, label=[L"\rho_{00}" L"\rho_{11}" L"\rho_{22}" L"\rho_{33}"])
xlabel!("t")
ylabel!(L"\rho")

tf = 20
U = solve_unitary(annealing, tf, alg = Tsit5(), abstol=1e-7, reltol=1e-7);
redfield_sol = solve_redfield(annealing, tf, U, alg = Tsit5(), abstol=1e-7, reltol=1e-7);

t_axis = range(0,tf,length=100)
p_computational_basis = [real(diag(redfield_sol(s))) for s in t_axis]
p_computational_basis = hcat(p_computational_basis...)
plot(t_axis, p_computational_basis', linewidth=2, label=[L"\rho_{00}" L"\rho_{11}" L"\rho_{22}" L"\rho_{33}"])
xlabel!("t")
ylabel!(L"\rho")

redfield_sol = solve_redfield(annealing, tf, U, alg = Tsit5(), abstol=1e-7, reltol=1e-7, callback=PositivityCheckCallback())


text:
include("centrality.jl")
using Combinatorics
using Random
using Test

function test3()
    dim = 11
    Random.seed!(1234)

    entries = collect(combinations(1:dim, 3))
    vals = rand(length(entries))
    T = SymTensor3([e[1] for e in entries],
                   [e[2] for e in entries],
                   [e[3] for e in entries],
                   vals)
    entry_vals = Dict{NTuple{3,Int64},Float64}()
    for (e, v) in zip(entries, vals)
        entry_vals[NTuple{3, Int64}(e)] = v
    end

    y = zeros(Float64, dim)
    x = rand(dim)
    for i in 1:dim, j in 1:dim, k in 1:dim
        key = NTuple{3,Int64}(sort([i, j, k]))
        if haskey(entry_vals, key)
            Tijk = entry_vals[key]
            y[i] += Tijk * x[j] * x[k]
        end
    end
    @test y ≈ apply(T, x)

    Y = zeros(Float64, dim, dim)
    for i in 1:dim, j in 1:dim, k in 1:dim
        key = NTuple{3,Int64}(sort([i, j, k]))
        if haskey(entry_vals, key)
            Tijk = entry_vals[key]
            Y[i, j] += Tijk * x[k]
        end
    end
    @test Y ≈ Matrix(reduce(T, x))
end

function test4()
    dim = 10
    Random.seed!(1234)

    entries = collect(combinations(1:dim, 4))
    vals = rand(length(entries))
    T = SymTensor4([e[1] for e in entries],
                   [e[2] for e in entries],
                   [e[3] for e in entries],
                   [e[4] for e in entries],
                   vals)
    entry_vals = Dict{NTuple{4,Int64},Float64}()
    for (e, v) in zip(entries, vals)
        entry_vals[NTuple{4,Int64}(e)] = v
    end

    y = zeros(Float64, dim)
    x = rand(dim)
    for i in 1:dim, j in 1:dim, k in 1:dim, l in 1:dim
        key = NTuple{4,Int64}(sort([i, j, k, l]))
        if haskey(entry_vals, key)
            Tijkl = entry_vals[key]
            y[i] += Tijkl * x[j] * x[k] * x[l]
        end
    end
    @test y ≈ apply(T, x)

    Y = zeros(Float64, dim, dim)
    for i in 1:dim, j in 1:dim, k in 1:dim, l in 1:dim
        key = NTuple{4,Int64}(sort([i, j, k, l]))
        if haskey(entry_vals, key)
            Tijkl = entry_vals[key]
            Y[i, j] += Tijkl * x[k] * x[l]
        end
    end
    @test Y ≈ Matrix(reduce(T, x))
end

function test5()
    dim = 9
    Random.seed!(1234)

    entries = collect(combinations(1:dim, 5))
    vals = rand(length(entries))
    T = SymTensor5([e[1] for e in entries],
                   [e[2] for e in entries],
                   [e[3] for e in entries],
                   [e[4] for e in entries],
                   [e[5] for e in entries],
                   vals)
    entry_vals = Dict{NTuple{5,Int64},Float64}()
    for (e, v) in zip(entries, vals)
        entry_vals[NTuple{5,Int64}(e)] = v
    end

    y = zeros(Float64, dim)
    x = rand(dim)
    for i in 1:dim, j in 1:dim, k in 1:dim, l in 1:dim, m in 1:dim
        key = NTuple{5,Int64}(sort([i, j, k, l, m]))
        if haskey(entry_vals, key)
            Tijklm = entry_vals[key]
            y[i] += Tijklm * x[j] * x[k] * x[l] * x[m]
        end
    end
    @test y ≈ apply(T, x)

    Y = zeros(Float64, dim, dim)
    for i in 1:dim, j in 1:dim, k in 1:dim, l in 1:dim, m in 1:dim
        key = NTuple{5,Int64}(sort([i, j, k, l, m]))
        if haskey(entry_vals, key)
            Tijklm = entry_vals[key]
            Y[i, j] += Tijklm * x[k] * x[l] * x[m]
        end
    end
    @test Y ≈ Matrix(reduce(T, x))
end

function main()
    test3()
    test4()
    test5()
end
;

max_stars_repo_name: joshua-matt/HO-HyperRank

max_issues_repo_name: joshua-matt/HO-HyperRank

max_forks_repo_name: joshua-matt/HO-HyperRank


text:
"""
This package facilitates working with probability distributions by means of Monte-Carlo methods, in a way that allows for propagation of probability distributions through functions. This is useful for, e.g.,  nonlinear [uncertainty propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty). A variable or parameter might be associated with uncertainty if it is measured or otherwise estimated from data. We provide two core types to represent probability distributions: `Particles` and `StaticParticles`, both `<: Real`. (The name "Particles" comes from the [particle-filtering](https://en.wikipedia.org/wiki/Particle_filter) literature.) These types all form a Monte-Carlo approximation of the distribution of a floating point number, i.e., the distribution is represented by samples/particles. **Correlated quantities** are handled as well, see [multivariate particles](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/#Multivariate-particles-1) below.

A number of type `Particles` behaves just as any other `Number` while partaking in calculations. After a calculation, an approximation to the **complete distribution** of the output is captured and represented by the output particles. `mean`, `std` etc. can be extracted from the particles using the corresponding functions. `Particles` also interact with [Distributions.jl](https://github.com/JuliaStats/Distributions.jl), so that you can call, e.g., `Normal(p)` and get back a `Normal` type from distributions or `fit(Gamma, p)` to get a `Gamma`distribution. Particles can also be iterated, asked for `maximum/minimum`, `quantile` etc. If particles are plotted with `plot(p)`, a histogram is displayed. This requires Plots.jl. A kernel-density estimate can be obtained by `density(p)` is StatsPlots.jl is loaded.

## Quick start
```julia
julia> using MonteCarloMeasurements, Plots

julia> a = π ± 0.1 # Construct Gaussian uncertain parameters using ± (\\pm)
Particles{Float64,2000}
 3.14159 ± 0.1

julia> b = 2 ∓ 0.1 # ∓ (\\mp) creates StaticParticles (with StaticArrays)
StaticParticles{Float64,100}
 2.0 ± 0.0999

julia> std(a)      # Ask about statistical properties
0.09999231528930486

julia> sin(a)      # Use them like any real number
Particles{Float64,2000}
 1.2168e-16 ± 0.0995

julia> plot(a)     # Plot them

julia> b = sin.(1:0.1:5) .± 0.1; # Create multivariate uncertain numbers

julia> plot(b)                   # Vectors of particles can be plotted

julia> using Distributions

julia> c = Particles(500, Poisson(3.)) # Create uncertain numbers distributed according to a given distribution
Particles{Int64,500}
 2.882 ± 1.7
```

For further help, see the [documentation](https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable), the [examples folder](https://github.com/baggepinnen/MonteCarloMeasurements.jl/tree/master/examples) or the [arXiv paper](https://arxiv.org/abs/2001.07625).
"""
module MonteCarloMeasurements
using LinearAlgebra, Statistics, Random, StaticArrays, RecipesBase, MacroTools, SLEEFPirates
using Distributed: pmap
import Base: add_sum

using Distributions, StatsBase, Requires

const DEFAULT_NUM_PARTICLES = 2000
const DEFAULT_STATIC_NUM_PARTICLES = 100

function pmean end

"""
The function used to reduce particles to a number for comparison. Defaults to `mean`. Change using `unsafe_comparisons`.
"""
const COMPARISON_FUNCTION = Ref{Function}(pmean)
const COMPARISON_MODE = Ref(:safe)

"""
    unsafe_comparisons(onoff=true; verbose=true)
Toggle the use of a comparison function without warning. By default `mean` is used to reduce particles to a floating point number for comparisons. This function can be changed, example: `set_comparison_function(median)`

    unsafe_comparisons(mode=:reduction; verbose=true)
One can also specify a comparison mode, `mode` can take the values `:safe, :montecarlo, :reduction`. `:safe` is the same as calling `unsafe_comparisons(false)` and `:reduction` corresponds to `true`.
If
"""
function unsafe_comparisons(mode=true; verbose=true)
    mode == false && (mode = :safe)
    mode == true && (mode = :reduction)
    COMPARISON_MODE[] = mode
    if mode != :safe && verbose
        if mode === :reduction
            @info "Unsafe comparisons using the function `$(COMPARISON_FUNCTION[])` has been enabled globally. Use `@unsafe` to enable in a local expression only or `unsafe_comparisons(false)` to turn off unsafe comparisons"
        elseif mode === :montecarlo
            @info "Comparisons using the monte carlo has been enabled globally. Call `unsafe_comparisons(false)` to turn off unsafe comparisons"
        end
    end
    mode ∉ (:safe, :montecarlo, :reduction) && error("Got unsupported comparison model")
end
"""
    set_comparison_function(f)

Change the Function used to reduce particles to a number for comparison operators
Toggle the use of a comparison Function without warning using the Function `unsafe_comparisons`.
"""
function set_comparison_function(f)
    if f in (mean, median, maximum, minimum)
        @warn "This comparison function ($(f)) is probably not the right choice, consider if you want the particle version (p$(f)) instead."
    end
    COMPARISON_FUNCTION[] = f
end

"""
    @unsafe expression
Activates unsafe comparisons for the provided expression only. The expression is surrounded by a try/catch block to robustly restore unsafe comparisons in case of exception.
"""
macro unsafe(ex)
    ex2 = if @capture(ex, assigned_vars__ = y_)
        if length(assigned_vars) == 1
            esc(assigned_vars[1])
        else
            esc.(assigned_vars[1].args)
        end
    else
        :(res)
    end
    quote
        previous_state = COMPARISON_MODE[]
        unsafe_comparisons(true, verbose=false)
        local res
        try
            res = ($(esc(ex)))
        finally
            unsafe_comparisons(previous_state, verbose=false)
        end
        $ex2 = res
    end
end

export ±, ∓, .., ⊠, ⊞, AbstractParticles,Particles,StaticParticles, MvParticles, sigmapoints, transform_moments, ≲,≳, systematic_sample, ess, outer_product, meanstd, meanvar, register_primitive, register_primitive_multi, register_primitive_single, ℝⁿ2ℝⁿ_function, ℝⁿ2ℂⁿ_function, ℂ2ℂ_function, ℂ2ℂ_function!, resample!, bootstrap, sqrt!, exp!, sin!, cos!, wasserstein, with_nominal, nominal, nparticles, particleeltype
# Plot exports
export errorbarplot, mcplot, ribbonplot

# Statistics reexport
export mean, std, cov, var, quantile, median
export pmean, pstd, pcov, pcor, pvar, pquantile, pmedian, pmiddle, piterate, pextrema, pminimum, pmaximum
# Distributions reexport
export Normal, MvNormal, Cauchy, Beta, Exponential, Gamma, Laplace, Uniform, fit, logpdf

export unsafe_comparisons, @unsafe, set_comparison_function

export bymap, bypmap, @bymap, @bypmap, @prob, Workspace, with_workspace, has_particles, mean_object

include("types.jl")
include("register_primitive.jl")
include("sampling.jl")
include("particles.jl")
include("distances.jl")
include("complex.jl")
include("sigmapoints.jl")
include("resampling.jl")
include("bymap.jl")
include("deconstruct.jl")
include("diff.jl")
include("plotting.jl")
include("optimize.jl")
include("sleefpirates.jl")
include("nominal.jl")

# This is defined here so that @bymap is loaded
LinearAlgebra.norm2(p::AbstractArray{<:AbstractParticles}) = bymap(LinearAlgebra.norm2,p)
Base.:\(x::AbstractVecOrMat{<:AbstractParticles}, y::AbstractVecOrMat{<:AbstractParticles}) = bymap(\, x, y)
Base.:\(x::Diagonal{<:AbstractParticles}, y::Vector{<:AbstractParticles}) = bymap(\, x, y) # required for ambiguity

function __init__()
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
    @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" include("measurements.jl")
    @require Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" include("unitful.jl")
end

end


text:
# TODO: use elements for any kind of object and use cell for the other!

"""
    Mesh

Holds mesh related attributes that are useful for the ray tracing algorithm, such as a
`kdtree` for nodes searches, the mesh bounding box and a mapping between nodes and cells
(`node_cells` and `cell_nodes`).
"""
struct Mesh{M,K,N,C,B}
    model::M
    kdtree::K
    node_cells::N
    cell_nodes::C
    bbmin::B
    bbmax::B
end

"""
    Mesh(model::UnstructuredDiscreteModel)

Builds a Mesh using data from a [`UnstructuredDiscreteModel`](@ref).
"""
function Mesh(model::UnstructuredDiscreteModel)
    grid = get_grid(model)
    kdtree = KDTree(grid)
    node_cells = get_faces(get_grid_topology(model), 0, num_cell_dims(model))
    cell_nodes = get_cell_node_ids(grid)
    bbmin, bbmax = bounding_box(grid)
    return Mesh(model, kdtree, node_cells, cell_nodes, bbmin, bbmax)
end

"""
    KDTree(grid::UnstructuredGrid)

Builds a [`KDTree`](@ref) using node data from a `grid`.
"""
function KDTree(grid::UnstructuredGrid{Dc,Dp,Tp}) where {Dc,Dp,Tp}
    nodes = get_node_coordinates(grid)
    snodes = convert.(SVector{Dp,Tp}, nodes) # Dc or Dp?
    return KDTree(snodes)
end

num_dims(mesh::Mesh) = num_dims(mesh.model)
num_cells(mesh::Mesh) = num_cells(mesh.model)
num_nodes(mesh::Mesh) = num_nodes(mesh.model)

"""
    bounding_box(grid::UnstructuredGrid)

Returns the bounding box of an unstructured `grid`.
"""
function bounding_box(grid::UnstructuredGrid{Dc,Dp,Tp}) where {Dc,Dp,Tp}
    nodes = get_node_coordinates(grid)

    xmin = MVector{Dp,Tp}(zeros(Tp, Dp)) # Dc or Dp?
    xmax = MVector{Dp,Tp}(zeros(Tp, Dp))

    for i in 1:Dp
        xs = getindex.(nodes, i)
        xmin[i] = min(xs...)
        xmax[i] = max(xs...)
    end

    bottomleft = convert(Point2D{Tp}, xmin)
    topright = convert(Point2D{Tp}, xmax)

    return bottomleft, topright
end

"""
    width(mesh::Mesh)

Returns the width of the rectangular mesh.
"""
@inline width(mesh::Mesh) = mesh.bbmax[1] - mesh.bbmin[1]

"""
    width(mesh::Mesh)

Returns the height of the rectangular mesh.
"""
@inline height(mesh::Mesh) = mesh.bbmax[2] - mesh.bbmin[2]

"""
    inboundary(mesh::Mesh, x::Point2D, [atol::Real=0]) -> Bool

Checks if a point `x` lies in the boundary of the `mesh` with certain absolute tolerance
`atol`.
"""
function inboundary(mesh::Mesh, x::Point2D, atol::Real=0)
    @unpack bbmin, bbmax = mesh
    return isapprox(x[1], bbmax[1], atol=atol) || isapprox(x[1], bbmin[1], atol=atol) ||
           isapprox(x[2], bbmax[2], atol=atol) || isapprox(x[2], bbmin[2], atol=atol)
end

"""
    find_element(mesh::Mesh, x::Point, k::Int=2)

Finds the cell or element of the `mesh` that contains a given point `x` by searching nodes
and elements using nearest neighbor searches.
"""
function find_element(mesh::Mesh, x::Point2D, k::Int=2)
    @unpack model, kdtree, node_cells, cell_nodes = mesh

    # get the nearest node id closest to `x`
    nn_id, nn_distance = nn(kdtree, x)

    # get the associated cell ids that contain the nearest node
    cell_ids = node_cells[nn_id]

    # loop over those cells until we find the one that contains the point `x`
    for cell_id in cell_ids
        node_ids = cell_nodes[cell_id]
        if point_in_element(mesh, node_ids, x)
            return cell_id
        end
    end

    # the mesh might be deformed, i.e. the cells that contain the nearest node do not
    # contain the point `x`. In that case, we need to search for more nodes. We could either
    # use `inrange` or `knn` for this purpose.
    nn_ids, nn_distances = knn(kdtree, x, k, true, i -> isequal(i, nn_id)) # TODO: cache!
    for node_id in nn_ids
        cell_ids = node_cells[node_id]
        for cell_id in cell_ids
            node_ids = cell_nodes[cell_id]
            if point_in_element(mesh, node_ids, x)
                return cell_id
            end
        end
    end

    # nn_ids = inrange(kdtree, x, factor * nn_distance, true)
    # for node_id in nn_ids
    #     cell_ids = node_cells[node_id]
    #     for cell_id in cell_ids
    #         node_ids = cell_nodes[cell_id]
    #         if point_in_element(grid, node_ids, x)
    #             return cell_id
    #         end
    #     end
    # end

    return -1
end

# Use dispatch once I get the info about the element type using Gridap topology.
point_in_element(mesh::Mesh, node_ids::AbstractVector{<:Int32}, x::Point2D) =
    point_in_triangle(mesh, node_ids, x)

"""
    point_in_triangle(mesh::Mesh, node_ids, x) -> Bool

Checks if whether a given point `x` lies inside, the edge or corner of the triangle given by
its node coordinates ids `node_ids`.
"""
function point_in_triangle(mesh::Mesh, node_ids::AbstractVector{<:Int32}, x::Point2D)
    @unpack model = mesh
    node_coordinates = get_node_coordinates(get_grid(model))

    x1, y1 = node_coordinates[node_ids[1]]
    x2, y2 = node_coordinates[node_ids[2]]
    x3, y3 = node_coordinates[node_ids[3]]

    R = @SMatrix [x1 x2 x3; y1 y2 y3; 1 1 1]
    r = @SVector [x[1], x[2], 1]
    λ = R \ r

    # return true if it lies in or on the triangle
    T = eltype(λ)
    tol = sqrt(eps(T))
    domain = ClosedInterval{T}(zero(T) - tol, one(T) + tol)
    return λ[1] in domain && λ[2] in domain && λ[3] in domain
    # return all(in.(λ, Ref(domain))) # allocates
end

"""
    point_in_quadrangle(mesh::Mesh, node_ids, x) -> Bool

Checks if whether a given point `x` lies inside the quadrangle given by its node coordinates
ids `node_ids`.
"""
function point_in_quadrangle(mesh::Mesh, node_ids::AbstractVector{<:Int32}, x::Point2D)

    # triangle node ids
    t_node_ids = @MVector zeros(3)

    # look on 4 triangles because we do not know the order of the nodes
    for i in 1:4
        for j in 1:3
            k = mod1(i + j - 1, 4) # k = (i + j - 2) % 4 + 1, k = mod(i + j - 1, 1:4)
            t_node_ids[j] = node_ids[k]
        end
        if point_in_triangle(mesh, t_node_ids, x)
            return true
        end
    end

    return false
end


text:
include("test_coverage.jl")
include("test_throughput.jl")
include("test_rpc.jl")

AMQPTestCoverage.runtests()
AMQPTestThroughput.runtests()
AMQPTestRPC.runtests()


text:
# references
#https://github.com/sctplab/usrsctp/blob/master/Manual.md
#https://tools.ietf.org/html/rfc6458#section-4.1.2
#https://github.com/JuliaLang/julia/blob/master/stdlib/Sockets/src/Sockets.jl
#https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/#Closure-cfunctions
#https://man7.org/linux/man-pages/man2/bind.2.html
module Sctp

include("UsrSctp.jl")
using .UsrSctp
include("defines.jl")
# using Sockets
import Sockets: IPv4, IPv6, @ip_str
export IPv4, @ip_str
export SockaddrIn
export sctp_init, SctpSocket
# import Base.close

"""
    Every application has to start with usrsctp_init(). This function calls sctp_init() and reserves the memory necessary to administer the data transfer.
    The flag `use_udp` determines whether to use SCTP over UDP. Defaults to false
"""
function sctp_init(use_udp=false)
    udp_port = use_udp ? 9899 : 0
    UsrSctp.usrsctp_init(udp_port, Ptr{Cvoid}(), Ptr{Cvoid}())
end

function sctp_finish()
  # @show

  ccall((:usrsctp_finish, "libusrsctp"), Cint, ())
end

struct SctpSocket
  ptr :: Ptr{UsrSctp.Socket}
end

function SctpSocket end

# function SCTPSocket(
#     domain :: Cint,
#     type :: Cint,
#     protocol :: Cint,
#     receive_cb :: Ptr{Cvoid},
#     send_cb :: Ptr{Cvoid},
#     sb_threshold :: Cuint,
function SctpSocket()
    ptr = UsrSctp.usrsctp_socket(
        AF_INET6,
        SOCK_SEQPACKET,
        SOL_SCTP,
        Ptr{Cvoid}(), #receive_cb
        Ptr{Cvoid}(), #send_cb
        0, # sb_threshold
        Ptr{Cvoid}()#ulp_info
    )
    if ptr === nothing # TODO if ptr == 0
        error("socket was not created")
    end
    return SctpSocket(ptr)
end

function SctpSocket(ptr :: Union{Ptr{}, Nothing})
  if ptr === nothing
    error("socket was not created")
  else
    return SctpSocket(ptr)
  end

end
function Base.bind(so :: SctpSocket,
                  addr :: IPv4,
                  port :: Integer = 0) # default 0, allow OS to choose ephemeral port
    p :: Cushort = UInt16(port)
    # sock_addr = UsrSctp.in_addr()
    sockaddr = UsrSctp.sockaddr_in(AF_INET,
                                hton(p),
                                UsrSctp.in_addr(hton(addr.host)),
                                (Cuchar(0) for _ in 1:8) |> collect |> Tuple)
    addrlen = sizeof(sockaddr)
    UsrSctp.usrsctp_bind(so.ptr, Ref(sockaddr), addrlen)
end

#=
function connect(so :: SCTPSocket,
                sockaddr :: SockaddrIn)
  addrlen = sizeof(sockaddr)
  if ccall((:usrsctp_connect, "libusrsctp"), Cint, (Ptr{Cvoid}, Ref{SockaddrIn}, Cint), so.ptr, sockaddr, addrlen) < 0
    printcerr()
  end
end

"""
    An application uses listen() to mark a socket as being able to accept new associations.
    if `backlog` is non-zero, enable listening, else disable listening.
"""
function listen(so :: SCTPSocket,
                backlog :: Int)
  ccall((:usrsctp_listen, "libusrsctp"), Cint, (Ptr{Cvoid}, Cint), so.ptr, backlog)
end

function accept(so :: SCTPSocket,
                addr :: IPv4,
                port :: Cushort) :: SCTPSocket
  sockaddr = SockaddrIn(addr, port)
  addrlen = sizeof(sockaddr)
  ptr = ccall((:usrsctp_bind, "libusrsctp"), Cint, (Ptr{Cvoid}, Ref{SockaddrIn}, Cint), so.ptr, sockaddr, addrlen)
  return SCTPSocket(ptr)
end

function Base.close(so :: SCTPSocket)
  ccall((:usrsctp_socket, "libusrsctp"), Cvoid, (Ptr{Cvoid},), so.ptr)
end

# return # of bytes sent, or -1 if error
function sendv(so :: SCTPSocket,
                  # can be nothing in the case of a connected socket
                  addr :: Union{SockaddrIn, Nothing},
                  data :: Vector{Cint})
  # info = nothing
  info = SctpSndinfo(
                     0,
  # snd_sid :: UInt16
                     SCTP_UNORDERED,
  # snd_flags :: UInt16
                     0,
  # snd_ppid :: UInt32
                    0,
  # snd_context :: UInt32
                    0
  # snd_assoc_id :: SctpAssocT
 )
  res = ccall((:usrsctp_sendv, "libusrsctp"), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{SockaddrIn}, Cint, Ptr{Cvoid}, Cuint, Cint, Cint),
        so.ptr,
        # "gather buffer", treated as user message
        Ref(data),
        sizeof(data), # number of elements in iov
        Ref(addr),
        # number of addresses (max 1)
        # 1
        addr === nothing ? 0 : 1,
        # Ref(info),
        Ref{Cvoid}(nothing),
        # sizeof(info), # length of info
        0,
        SCTP_SENDV_NOINFO,
        # SCTP_SENDV_SNDINFO,
        0 # flags
       )

  if res < 0
    printcerr()
  end
end

function printcerr()
  println("Error $(Libc.errno()): $(Libc.strerror())")
end
function recvv(so :: SCTPSocket)
  dbuf = Array{Cint}(undef, 2048)
  @show from = Ref{SockaddrIn}()
  fromlen = Ref(Cint(sizeof(from)))
  info = nothing
  infolen = Ref(Cint(sizeof(info)))
  bytes_recvd = ccall((:usrsctp_recvv, "libusrsctp"), Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Cint, Ptr{SockaddrIn}, Ptr{Cint}, Ptr{Cvoid}, Ptr{Cint}, Cuint, Cint),
        so.ptr,
        # "gather buffer", treated as user message
        dbuf,
        sizeof(dbuf), # number of elements in iov
        from,
        fromlen,
        Ref(info),
        infolen,
        SCTP_RECVV_NOINFO,
        0 # flags
       )
  @show bytes_recvd
  @show from
  return dbuf
end

function remote_udp_encaps_port(so :: SCTPSocket)
  getsockopt(so, SCTP_REMOTE_UDP_ENCAPS_PORT)
# SCTP_REMOTE_UDP_ENCAPS_PORT
end

function remote_udp_encaps_port!(so :: SCTPSocket, port :: Integer)
  udp_encaps = SctpUdpencaps(
                             SockaddrIn(ip"127.0.0.1", 7),
                             0,
                             Cushort(port)
                            )
  setsockopt(so, SCTP_REMOTE_UDP_ENCAPS_PORT, udp_encaps)
# SCTP_REMOTE_UDP_ENCAPS_PORT
end

function getsockopt(so :: SCTPSocket,
                    optname :: Integer)
  optval = Ref(Array{Cuint}(undef, 2048))
  @show optlen = Ref{Cuint}(sizeof(optval.x))
  funcall = ccall((:usrsctp_getsockopt, "libusrsctp"), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}, Ptr{Cuint}),
        so.ptr,
        SOL_SCTP,
        Cint(optname),
        optval,
        optlen
       )
  if funcall < 0
    printcerr()
    return
  end
  @show optlen
  result = optval.x[1:optlen.x]
  return result
end

function setsockopt(so :: SCTPSocket,
                    optname :: Integer,
                    optval)
  # optval = Ref(Array{Cuint}(undef, 2048))
  @show optlen = sizeof(optval)
  funcall = ccall((:usrsctp_setsockopt, "libusrsctp"), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cvoid}, Cint),
        so.ptr,
        SOL_SCTP,
        Cint(optname),
        Ref(optval),
        optlen
       )
  if funcall < 0
    printcerr()
    return
  end
  # result = optval.x[1:optlen.x]
end

=#

end


text:
using LabView0mqJlExamples
using Test


text:
# Autogenerated wrapper script for PPL_jll for aarch64-linux-musl-cxx11
export libppl, libppl_c, ppl_config, ppl_lcdd, ppl_pips

using GMP_jll
JLLWrappers.@generate_wrapper_header("PPL")
JLLWrappers.@declare_library_product(libppl, "libppl.so.14")
JLLWrappers.@declare_library_product(libppl_c, "libppl_c.so.4")
JLLWrappers.@declare_executable_product(ppl_config)
JLLWrappers.@declare_executable_product(ppl_lcdd)
JLLWrappers.@declare_executable_product(ppl_pips)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libppl,
        "lib/libppl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libppl_c,
        "lib/libppl_c.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        ppl_config,
        "bin/ppl-config",
    )

    JLLWrappers.@init_executable_product(
        ppl_lcdd,
        "bin/ppl_lcdd",
    )

    JLLWrappers.@init_executable_product(
        ppl_pips,
        "bin/ppl_pips",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using NumCME
using StaticArrays
using Test

states = [[@MVector([0, i]) for i in 0:10];[@MVector([1, i]) for i in 0:10]]
vals = [[0.0 for i in 0:10]; [1.0 for i in 0:10]]
v = FspVectorSparse(states, vals)
v1 = sum(v, [1])
v2 = sum(v, [2])
@test typeof(v) <: AbstractFspVector
@test typeof(v) <: FspVectorSparse{2, <:Integer, <:Real}
@test sum(v) == 11.0
@test nnz(v1) == 11
@test nnz(v2) == 2

@test_throws ArgumentError FspVectorSparse([@MVector([0, 1])], [1.0,2.0])
@test_throws ArgumentError FspVectorSparse([@MVector([0,1]), @MVector([1,2])], [1.0])


text:
export check_dΠψ, _dψ1dθρ, _ρ, _ψ2, _ψ1

@inline _ρ(θρ::Real) = logistic(θρ)
@inline _dρdθρ(θρ::Real) = (z = logistic(θρ); z*(1-z) )

@inline _ρ2(θρ::Real) = _ρ(θρ)^2
_dρdσ = _dρdθρ # alias

@inline _ψ1(u::Real,v::Real,ρ::Real) = ρ*u + sqrt(1-ρ^2)*v
@inline _ψ2(u::Real,v::Real, ρ::Real) = _ψ2(u,v)
@inline _ψ2(u::Real,v::Real) = u

@inline _dψ1dρ(u::Real,v::Real,ρ::Real) = u - ρ/sqrt(1-ρ^2)*v

@inline _dψ1dθρ(u::Real, v::Real, ρ::Real, θρ::Real) = _dψ1dρ(u,v,ρ)*_dρdθρ(θρ::Real)

@inline _z(x2::Real, x1::Real, Δ::Real, ρ::Real) = (x2 - ρ*x1 + Δ)/sqrt(1-ρ^2)

# derivatives
@inline _dzdρ(x2::Real, x1::Real, ρ::Real, z::Real) = -x1/sqrt(1-ρ^2) + ρ*z/(1-ρ^2)
@inline _dπdρ(x2::Real, x1::Real, Δ::Real, ρ::Real) = (z = _z(x2,x1,Δ,ρ);  normpdf(z) * _dzdρ(x2,x1,ρ,z))

# finite difference versions
@inline _ρ(x::Real, h::Real) = _ρ(x+h)
@inline _z(x2::Real, x1::Real, Δ::Real, ρ::Real, h::Real) = _z(x2, x1+h, Δ, ρ)

# ------------------------------ matrix updates -------------------------

function _βΠψdθρ!(P::AbstractMatrix, y::StepRangeLen, θρ::Real, β::Real)
    n = length(y)
    (n,n) == size(P) || throw(DimensionMismatch())
    ρ = _ρ(θρ)
    dρdθρ = _dρdθρ(θρ)
    Δ = 0.5 * step(y)

    @inbounds for (j,yj) in enumerate(y)
        if j == 1
            @. P[ :,j] = β * _dπdρ(yj,y,Δ,ρ) * dρdθρ
        elseif j == n
            @. P[ :,j] = -β * _dπdρ(yj,y,-Δ,ρ) * dρdθρ
        else
            @. P[ :,j] = β * ( _dπdρ(yj,y,Δ,ρ) - _dπdρ(yj,y,-Δ,ρ) ) * dρdθρ
        end
    end
end

function _βΠψ!(P::AbstractMatrix, y1::StepRangeLen, y2::StepRangeLen, θp::Real, β::Real)
    n = length(y2)
    n == length(y1) || throw(DimensionMismatch())
    (n,n) == size(P) || throw(DimensionMismatch())
    ρ = _ρ(θp)
    Δ = 0.5 * step(y2)

    @inbounds for (j,yj) in enumerate(y2)
        if j == 1
            @. P[ :,j] = β * normcdf(_z(yj, y1, Δ, ρ))
        elseif j == n
            @. P[ :,j] = β * normccdf(_z(yj, y1, -Δ, ρ))
        else
            @. P[ :,j] = β * ( normcdf(_z(yj, y1, Δ, ρ)) - normcdf(_z(yj, y1, -Δ, ρ) ))
        end
    end
end

_βΠψ!(P::AbstractMatrix, y1::StepRangeLen, θp::Real, β::Real) = _βΠψ!(P, y1, y1, θp, β)

function _βΠψ(y1::StepRangeLen{T}, θp::Real, β::Real) where {T}
    P = Matrix{T}(undef, length(y1), length(y1))
    _βΠψ!(P, y1, y1, θp, β)
    return P
end

# ------------------------------ derivative check -------------------------

function check_dΠψ(θρ::Real, ψspace::StepRangeLen)

    Δ = 0.5 * step(ψspace)

    for yj in ψspace
        for y in ψspace
            fdσ = Calculus.derivative((sig) -> normcdf(_z(yj, y, Δ, _ρ(sig))), θρ)
            dσ  = _dρdθρ(θρ) * _dπdρ(yj, y, Δ, _ρ(θρ))
            abs(fdσ - dσ ) < 1e-7 || throw(error("bad σ grad at σ = $σ, ψ2 = $yj, ψ1 = $y"))
        end
    end
    return true
end


text:
# ============================================================================ #
# Point
# ============================================================================ #
abstract type AbstractPoint end

struct Point{T} <: AbstractPoint
    x::T
    y::T
    function Point{T}(x::T, y::T) where {T}
        new(x, y)
    end
end

Point(x::T, y::T) where {T} = Point{T}(x, y)
Point(x, y) = Point( promote(x, y)... )
Point{T}(x, y) where {T} = Point( promote( convert(T, x), convert(T, y))... )
Point(a::A) where {A<:AbstractArray} = Point(a[1], a[2])
Point() = Point(0.0, 0.0)

mutable struct MutablePoint{T} <: AbstractPoint
    x::T
    y::T
    function MutablePoint{T}(x::T, y::T) where {T}
        new(x, y)
    end
end

MutablePoint(x::T, y::T) where {T} = MutablePoint{T}(x, y)
MutablePoint(x, y) = MutablePoint( promote(x, y)... )
MutablePoint{T}(x, y) where {T} = MutablePoint( promote( convert(T, x), convert(T, y))... )
MutablePoint(a::A) where {A<:AbstractArray} = MutablePoint(a[1], a[2])
MutablePoint() = MutablePoint(0.0, 0.0)

# ---------------------------------------------------------------------------- #
# Methods
# ---------------------------------------------------------------------------- #

function Base.getindex(p::T, index::Int) where {T <: AbstractPoint}
    if index == 1
        return p.x
    elseif index == 2
        return p.y
    else
        return NaN
    end
end

function Base.size(p::T) where {T <: AbstractPoint}
    return (2,)
end

function Base.size(p::T, dim = 1) where {T <: AbstractPoint}
    if dim == 1
        return 2
    elseif dim == 2
        return 1
    else
        return 0
    end
end

Base.:(==)(p1::AbstractPoint, p2::AbstractPoint) = p1.x == p2.x && p1.y == p2.y

Base.:+(p1::AbstractPoint, p2::AbstractPoint) = [p1.x + p2.x, p1.y + p2.y]
Base.:-(p1::AbstractPoint, p2::AbstractPoint) = [p1.x - p2.x, p1.y - p2.y]

Base.:+(p1::AbstractPoint, p2::AbstractArray) = [p1.x + p2[1], p1.y + p2[2]]
Base.:-(p1::AbstractPoint, p2::AbstractArray) = [p1.x - p2[1], p1.y - p2[2]]

Base.:+(p1::AbstractArray, p2::AbstractPoint) = [p1[1] + p2.x, p1[2] + p2.y]
Base.:-(p1::AbstractArray, p2::AbstractPoint) = [p1[1] - p2.x, p1[2] - p2.y]

Base.:+(p1::AbstractArray, p2::T) where {T <: Number} = [p1[1] + p2, p1[2] + p2]
Base.:-(p1::AbstractArray, p2::T) where {T <: Number} = [p1[1] - p2, p1[2] - p2]

Base.:+(p1::AbstractPoint, p2::T) where {T <: Number} = [p1[1] + p2, p1[2] + p2]
Base.:-(p1::AbstractPoint, p2::T) where {T <: Number} = [p1[1] - p2, p1[2] - p2]

Base.:*(f::N, point::P) where {N <: Number, P <: AbstractPoint} = P(f * point.x, f * point.y)
Base.:*(point::P, f::N) where {N <: Number, P <: AbstractPoint} = f * point;

Base.:/(point::P, f::N) where {N <: Number, P <: AbstractPoint} =  P( point.x / f, point.y / f);

function shift(p::T, v) where {T <: AbstractPoint}
    return shift(p::T, v[1], v[2])
end

function shift(p::T, x::S, y::S) where {S, T <: AbstractPoint}
    return T(p.x + x, p.y + y)
end
function scale(p::T, factor) where {T <: AbstractPoint}
    return factor * p
end
function rotate(p::T, rot_mat) where {T <: AbstractPoint}
    x = rot_mat[1,1] * p.x + rot_mat[1,2] * p.y
    y = rot_mat[2,1] * p.x + rot_mat[2,2] * p.y
    return T(x, y)
end

function random_point(n::Int, t::Type{PointType})  where {PointType <: AbstractPoint}
    # point x, y in the range of [0,1) x [0,1)

    point_array = Vector{t}(undef, n)
    x = rand(n)
    y = rand(n)
    for kk = 1:n
        point_array[kk] = t(x[kk], y[kk])
    end
    return point_array
end

random_point(n::Int) = random_point(n, Point)

function discretize_for_plot(point::AbstractPoint)
    return (point.x, point.y)
end

function discretize_for_plot(array::Vector{<:AbstractPoint}, closed_curve::Bool = false)
    L = length(array)
    if L <= 0
        return ([],[])
    end
    Lc = closed_curve ? 1 : 0

    ElementType = typeof(array[1].x)
    x_ = Vector{ElementType}(undef, L + Lc)
    y_ = Vector{ElementType}(undef, L + Lc)
    for kk = 1:L
        x_[kk] = array[kk].x
        y_[kk] = array[kk].y
    end
    if closed_curve
        x_[end] = array[1].x
        y_[end] = array[1].y
    end
    return (x_, y_)
end

function x(array::Vector{<:AbstractPoint})
    return map(p->p.x, array)
end

function y(array::Vector{<:AbstractPoint})
    return map(p->p.y, array)
end

function point_array_to_matrix(point_array)

    L = length(point_array)
    mat = zeros(L,2)

    for kk = 1:L
        mat[kk, 1] = point_array[kk].x
        mat[kk, 2] = point_array[kk].y
    end

    return mat
end

function matrix_to_point_array(mat)
    L = size(mat, 1)
    pa = Vector{MutablePoint}(undef, L)
    for kk = 1:L
        pa[kk] = MutablePoint(mat[kk,1], mat[kk,2])
    end
    return pa
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
Pilot response object. Limits pilot compliance with CAS instructions.
"""
Base.@kwdef mutable struct PilotResponse
    compliant::Bool = true       # whether pilot complies with CAS advisories
	turn_rate::R    = 3.0        # complying turn rate
	vert_rate::R    = 10.0       # complying vertical rate
	cmd::Command    = Command()  # resultant effective command
end

function initialize(pr::PilotResponse)
	pr.cmd = Command()
end

function update(pr::PilotResponse, cmd::Command, cas::CAS)
	pr.cmd = deepcopy(cmd)
    pr.compliant || return

	# horizontal
	if cas.hra == :cw
		pr.cmd.psi_d = pr.turn_rate
	elseif cas.hra == :ccw
		pr.cmd.psi_d = -pr.turn_rate
	end

	# vertical
	if cas.vra == :asc
		pr.cmd.h_d = pr.vert_rate
	elseif cas.vra == :desc
		pr.cmd.h_d = -pr.vert_rate
	end
end


text:
using GLVisualize, Colors, GeometryTypes, GLAbstraction, Reactive

if !isdefined(:runtests)
	window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end
const N = 87
const range = linspace(-5f0, 5f0, N)

function contour_inner(i, x, y)
    Float32(sin(1.3*x*i)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x))
end
const data = zeros(Float32, N, N)

function contourdata(t)
    for i=1:size(data, 1)
        @simd for j=1:size(data, 2)
            @inbounds data[i,j] = contour_inner(t, range[i], range[j])
        end
    end
    data
end

heightfield = map(contourdata, timesignal)
mini = Vec3f0(first(range), first(range), minimum(value(heightfield)))
maxi = Vec3f0(last(range), last(range), maximum(value(heightfield)))
barsvis = visualize(
    heightfield,
    scale_x = 0.07,
    scale_y = 0.07,
    color_map=map(RGBA{U8}, colormap("Blues")),
    color_norm=Vec2f0(0,1),
    ranges=(range, range),
    boundingbox=Signal(AABB{Float32}(mini, maxi))
)
_view(barsvis, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# Main functions concerning "doubt coefficient" for work on MCL
# Berkeley 2015 - Maxime Legrand

#Pkg.add("Gadfly")
using Gadfly
include("../plot_tools.jl")
include("../matrix_tools.jl")
include("mcl.jl")

#####################################
# Step by step (customed algorithm) #
#####################################

"""Apply MCL algorithm with expansion parameter e, inflation parameter r and p iterations
Customed version for doubt coefficient"""
function mcl(A;e=2,r=3,p=20,dc="none",de=(if e==2 0 else 1 end),dr=(r*0.2))
    A = mcl_norm(A)
    if dc=="none"
        for i=1:p
            A = mcl_norm(mcl_inflate(mcl_expand(A,e),r))
        end
        return A
    elseif dc=="parallel"
        B1,B2,B3,B4 = copy(A),copy(A),copy(A),copy(A)
        for i=1:p
            A = mcl_norm(mcl_inflate(mcl_expand(A,e),r))
            B1 = mcl_norm(mcl_inflate(mcl_expand(B1,e+de),r+dr))
            B2 = mcl_norm(mcl_inflate(mcl_expand(B2,e-de),r+dr))
            B3 = mcl_norm(mcl_inflate(mcl_expand(B3,e-de),r-dr))
            B4 = mcl_norm(mcl_inflate(mcl_expand(B4,e+de),r-dr))
        end
        dc = maximum([norm(A-B1), norm(A-B2), norm(A-B3), norm(A-B4)])
        return (A,dc)
    elseif dc=="forks"
        dc = zeros(p)
        for i=1:p
            A0 = copy(A)
            A = mcl_norm(mcl_inflate(mcl_expand(A0,e),r))
            B1 = mcl_norm(mcl_inflate(mcl_expand(A0,e+de),r+dr))
            B2 = mcl_norm(mcl_inflate(mcl_expand(A0,e-de),r+dr))
            B3 = mcl_norm(mcl_inflate(mcl_expand(A0,e-de),r-dr))
            B4 = mcl_norm(mcl_inflate(mcl_expand(A0,e+de),r-dr))
            dc[i] = maximum([norm(A-B1), norm(A-B2), norm(A-B3), norm(A-B4)])
        end
        return (A,dc)
    end
end

#################
# Input-focused #
#################

##################
# Output-focused #
##################

# Tools

"""SBM form
Input : A, adjacency matrix
    C, vector of clusters (represented as vectors as well)
Output : B, reorganised adjacency matrix from A to fit SBM output common outfit"""
function sbm_form(A,C)
    # Vertices permutation
    n = size(A,1)
    K = size(C,1)
    v = zeros(n)
    s = 1
    for k = 1:K
        for i = 1:(size(C[k],1))
            v[s] = C[k][i]
            s += 1
        end
    end
    # Output making
    B = zeros(n,n)
    for i = 1:n
        for j = 1:n
            B[i,j] = A[v[i],v[j]]
        end
    end
    return B
end

# SBM Prototype

"""SBM ML parameters
Input : A, adjacency matrix in SBM form
    l, vector of the sizes of the clusters
Output : p, vector of their intra-cluster connection probability
    p0, inter-cluster connection probability"""
function sbm_parameters(A, l)
    n = size(A,1)
    K = size(l,1)
    p = zeros(K)
    s1 = 0
    p0 = 0
    for k1 = 1:K
        s2 = 0
        for k2 = 1:K
            for i = 1:l[k1]
                for j = 1:l[k2]
                    if k1 == k2
                        p[k1] += A[s1+i,s2+j]
                    else
                        p0 += A[s1+i,s2+j]
                    end
                end
            end
            s2 += l[k2]
        end
        s1 += l[k1]
    end
    l2 = map(x -> x^2, l)
    p = p./l2
    n0 = n^2-sum(l2)
    if n0 != 0
        p0 = p0/n0
    end
    return (p,p0)
end

# DC-SBM Prototype

################################
# Partition of the input space #
################################


text:
# compute R²
function rsq(Y, X)
    β = X\Y
    errors = Y - X*β
    R² = 1.0 - sum(errors.^2.0)/sum((Y .- mean(Y)).^2.0)
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
using Base.Test, StatsFuns, DataFrames, GLM

function test_show(x)
    io = IOBuffer()
    show(io, x)
end

const glm_datadir = joinpath(dirname(@__FILE__), "..", "data")

## Formaldehyde data from the R Datasets package
form = DataFrame(Carb=[0.1,0.3,0.5,0.6,0.7,0.9],OptDen=[0.086,0.269,0.446,0.538,0.626,0.782])
lm1 = fit(LinearModel, OptDen ~ Carb, form)
test_show(lm1)
@test_approx_eq coef(lm1) collect(linreg(convert(Array, form[:Carb]), convert(Array, form[:OptDen])))
Σ = [6.136653061224592e-05 -9.464489795918525e-05
    -9.464489795918525e-05 1.831836734693908e-04]
@test_approx_eq vcov(lm1) Σ
@test_approx_eq cor(lm1.model) diagm(diag(Σ))^(-1/2)*Σ*diagm(diag(Σ))^(-1/2)
@test df(lm1) == 3
@test_approx_eq deviance(lm1) 0.0002992000000000012
@test_approx_eq loglikelihood(lm1) 21.204842144047973
@test_approx_eq nulldeviance(lm1) 0.3138488333333334
@test_approx_eq nullloglikelihood(lm1) 0.33817870295676444
@test r²(lm1) == r2(lm1)
@test_approx_eq r²(lm1) 0.9990466748057584
@test adjr²(lm1) == adjr2(lm1)
@test_approx_eq adjr²(lm1) 0.998808343507198
@test_approx_eq aic(lm1) -36.409684288095946
@test_approx_eq aicc(lm1) -24.409684288095946
@test_approx_eq bic(lm1) -37.03440588041178

dobson = DataFrame(Counts=[18.,17,15,20,10,20,25,13,12], Outcome=gl(3,1,9), Treatment=gl(3,3))
gm1 = fit(GeneralizedLinearModel, Counts ~ Outcome + Treatment, dobson, Poisson())
test_show(gm1)
@test df(gm1) == 5
@test_approx_eq deviance(gm1) 5.12914107700115
@test_approx_eq loglikelihood(gm1) -23.380659200978837
@test_approx_eq aic(gm1) 56.76131840195767
@test_approx_eq aicc(gm1) 76.76131840195768
@test_approx_eq bic(gm1) 57.74744128863877
@test_approx_eq coef(gm1)[1:3] [3.044522437723423,-0.45425527227759555,-0.29298712468147375]

## Example from http://www.ats.ucla.edu/stat/r/dae/logit.htm
admit = readtable(joinpath(glm_datadir, "admit.csv.gz"))
admit[:rank] = pool(admit[:rank])

for distr in (Binomial, Bernoulli)
    gm2 = fit(GeneralizedLinearModel, admit ~ gre + gpa + rank, admit, distr())
    test_show(gm2)
    @test df(gm2) == 6
    @test_approx_eq deviance(gm2) 458.5174924758994
    @test_approx_eq loglikelihood(gm2) -229.25874623794968
    @test_approx_eq aic(gm2) 470.51749247589936
    @test_approx_eq aicc(gm2) 470.7312329339146
    @test_approx_eq bic(gm2) 494.4662797585473
    @test_approx_eq coef(gm2) [-3.9899786606380734,0.0022644256521549043,0.8040374535155766,-0.6754428594116577,-1.3402038117481079,-1.5514636444657492]
end

gm3 = fit(GeneralizedLinearModel, admit ~ gre + gpa + rank, admit, Binomial(), ProbitLink())
test_show(gm3)
@test df(gm3) == 6
@test_approx_eq deviance(gm3) 458.4131713833386
@test_approx_eq loglikelihood(gm3) -229.20658569166932
@test_approx_eq aic(gm3) 470.41317138333864
@test_approx_eq aicc(gm3) 470.6269118413539
@test_approx_eq bic(gm3) 494.36195866598655
@test_approx_eq coef(gm3) [-2.3867922998680786,0.0013755394922972369,0.47772908362647015,-0.4154125854823675,-0.8121458010130356,-0.9359047862425298]

gm4 = fit(GeneralizedLinearModel, admit ~ gre + gpa + rank, admit, Binomial(), CauchitLink())
test_show(gm4)
@test df(gm4) == 6
@test_approx_eq deviance(gm4) 459.3401112751141
@test_approx_eq loglikelihood(gm4) -229.6700556375571
@test_approx_eq aic(gm4) 471.3401112751142
@test_approx_eq aicc(gm4) 471.5538517331295
@test_approx_eq bic(gm4) 495.28889855776214

gm5 = fit(GeneralizedLinearModel, admit ~ gre + gpa + rank, admit, Binomial(), CloglogLink())
test_show(gm5)
@test df(gm5) == 6
@test_approx_eq deviance(gm5) 458.89439629612616
@test_approx_eq loglikelihood(gm5) -229.44719814806314
@test_approx_eq aic(gm5) 470.8943962961263
@test_approx_eq aicc(gm5) 471.1081367541415
@test_approx_eq bic(gm5) 494.8431835787742

## Example with offsets from Venables & Ripley (2002, p.189)
anorexia = readtable(joinpath(glm_datadir, "anorexia.csv.gz"))
anorexia[:Treat] = pool(anorexia[:Treat])

gm6 = fit(GeneralizedLinearModel, Postwt ~ Prewt + Treat, anorexia, Normal(), IdentityLink(), offset=convert(Array, anorexia[:Prewt]))
test_show(gm6)
@test df(gm6) == 5
@test_approx_eq deviance(gm6) 3311.262619919613
@test_approx_eq loglikelihood(gm6) -239.9866487711122
@test_approx_eq aic(gm6) 489.9732975422244
@test_approx_eq aicc(gm6) 490.8823884513153
@test_approx_eq bic(gm6) 501.35662813730465
@test_approx_eq coef(gm6) [49.7711090149846,-0.5655388496391,-4.0970655280729,4.5630626529188]
@test_approx_eq GLM.dispersion(gm6.model, true) 48.6950385282296
@test_approx_eq stderr(gm6) [13.3909581420259,0.1611823618518,1.8934926069669,2.1333359226431]

gm7 = fit(GeneralizedLinearModel, Postwt ~ Prewt + Treat, anorexia, Normal(), LogLink(), offset=convert(Array, anorexia[:Prewt]),
          convTol=1e-8)
test_show(gm7)
@test_approx_eq deviance(gm7) 3265.207242977156
@test_approx_eq coef(gm7) [3.992326787835955,-0.994452693131178,-0.050698258703974,0.051494029957641]
@test_approx_eq GLM.dispersion(gm7.model, true) 48.01787789178518
@test_approx_eq stderr(gm7) [0.157167944259695,0.001886285986164,0.022584069426311,0.023882826190166]

## Gamma example from McCullagh & Nelder (1989, pp. 300-2)
clotting = DataFrame(u = log.([5,10,15,20,30,40,60,80,100]),
                     lot1 = [118,58,42,35,27,25,21,19,18])
gm8 = fit(GeneralizedLinearModel, lot1 ~ u, clotting, Gamma())
test_show(gm8)
@test df(gm8) == 3
@test_approx_eq deviance(gm8) 0.016729715178484157
@test_approx_eq loglikelihood(gm8) -15.994961974777247
@test_approx_eq aic(gm8) 37.989923949554495
@test_approx_eq aicc(gm8) 42.78992394955449
@test_approx_eq bic(gm8) 38.58159768156315
@test_approx_eq coef(gm8) [-0.01655438172784895,0.01534311491072141]
@test_approx_eq GLM.dispersion(gm8.model, true) 0.002446059333495581
@test_approx_eq stderr(gm8) [0.0009275466067257,0.0004149596425600]

gm9 = fit(GeneralizedLinearModel, lot1 ~ u, clotting, Gamma(), LogLink(), convTol=1e-8)
test_show(gm9)
@test df(gm9) == 3
@test_approx_eq deviance(gm9) 0.16260829451739
@test_approx_eq loglikelihood(gm9) -26.24082810384911
@test_approx_eq aic(gm9) 58.48165620769822
@test_approx_eq aicc(gm9) 63.28165620769822
@test_approx_eq bic(gm9) 59.07332993970688
@test_approx_eq coef(gm9) [5.50322528458221,-0.60191617825971]
@test_approx_eq GLM.dispersion(gm9.model, true) 0.02435442293561081
@test_approx_eq stderr(gm9) [0.19030107482720,0.05530784660144]

gm10 = fit(GeneralizedLinearModel, lot1 ~ u, clotting, Gamma(), IdentityLink(), convTol=1e-8)
test_show(gm10)
@test df(gm10) == 3
@test_approx_eq deviance(gm10) 0.60845414895344
@test_approx_eq loglikelihood(gm10) -32.216072437284176
@test_approx_eq aic(gm10) 70.43214487456835
@test_approx_eq aicc(gm10) 75.23214487456835
@test_approx_eq bic(gm10) 71.02381860657701
@test_approx_eq coef(gm10) [99.250446880986,-18.374324929002]
@test_approx_eq GLM.dispersion(gm10.model, true) 0.1041772704067886
@test_approx_eq stderr(gm10) [17.864388462865,4.297968703823]

# Logistic regression using aggregated data and weights
admit_agr = DataFrame(count=[28, 97, 93, 55, 33, 54, 28, 12],
                      admit=repeat([false, true], inner=[4]),
                      rank=PooledDataArray(repeat([1, 2, 3, 4], outer=[2])))
for distr in (Binomial, Bernoulli)
    gm14 = fit(GeneralizedLinearModel, admit ~ rank, admit_agr, distr(), wts=Vector{Float64}(admit_agr[:count]))
    @test df(gm14) == 4
    @test nobs(gm14) == 400
    @test_approx_eq deviance(gm14) 474.9667184280627
    @test_approx_eq loglikelihood(gm14) -237.48335921403134
    @test_approx_eq aic(gm14) 482.96671842822883
    @test_approx_eq aicc(gm14) 483.0679842510136
    @test_approx_eq bic(gm14) 498.9325766164946
    @test_approx_eq_eps coef(gm14) [0.16430305129127593,-0.7500299832244263,-1.364697929944679,-1.6867286645732025] 1e-5
end

# Logistic regression using aggregated data with proportions of successes and weights
admit_agr2 = DataFrame(count=[61, 151, 121, 67],
                       admit=[33, 54, 28, 12],
                       rank=PooledDataArray([1, 2, 3, 4]))
admit_agr2[:p] = Array(admit_agr2[:admit]) ./ Array(admit_agr2[:count])

gm15 = fit(GeneralizedLinearModel, p ~ rank, admit_agr2, Binomial(), wts=Vector{Float64}(admit_agr2[:count]))
test_show(gm15)
@test df(gm15) == 4
@test nobs(gm15) == 400
@test_approx_eq deviance(gm15) -2.4424906541753456e-15
@test_approx_eq loglikelihood(gm15) -9.50254433604239
@test_approx_eq aic(gm15) 27.00508867208478
@test_approx_eq aicc(gm15) 27.106354494869592
@test_approx_eq bic(gm15) 42.970946860516705
@test_approx_eq coef(gm15) [0.1643030512912767,-0.7500299832303851,-1.3646980342693287,-1.6867295867357475]

# Weighted Gamma example (weights are totally made up)
gm16 = fit(GeneralizedLinearModel, lot1 ~ u, clotting, Gamma(), wts=[1.5,2.0,1.1,4.5,2.4,3.5,5.6,5.4,6.7])
test_show(gm16)
@test df(gm16) == 3
@test nobs(gm16) == 32.7
@test_approx_eq deviance(gm16) 0.03933389380881689
@test_approx_eq loglikelihood(gm16) -43.35907878769152
@test_approx_eq aic(gm16) 92.71815757538305
@test_approx_eq aicc(gm16) 93.55439450918095
@test_approx_eq bic(gm16) 97.18028280909267
@test_approx_eq coef(gm16) [-0.017217012615523237,0.015649040411276433]

# Weighted Poisson example (weights are totally made up)
gm17 = fit(GeneralizedLinearModel, Counts ~ Outcome + Treatment, dobson, Poisson(),
           wts=[1.5,2.0,1.1,4.5,2.4,3.5,5.6,5.4,6.7])
test_show(gm17)
@test df(gm17) == 5
@test_approx_eq deviance(gm17) 17.699857821414266
@test_approx_eq loglikelihood(gm17) -84.57429468506352
@test_approx_eq aic(gm17) 179.14858937012704
@test_approx_eq aicc(gm17) 181.39578038136298
@test_approx_eq bic(gm17) 186.5854647596431
@test_approx_eq coef(gm17) [3.1218557035404793,  -0.5270435906931427,-0.40300384148562746,
                           -0.017850203824417415,-0.03507851122782909]

## Fitting GLMs with sparse matrices
srand(1)
X = sprand(1000, 10, 0.01)
β = randn(10)
y = Bool[rand() < logistic(x) for x in X * β]

gmsparse = fit(GeneralizedLinearModel, X, y, Binomial())
gmdense = fit(GeneralizedLinearModel, full(X), y, Binomial())

@test_approx_eq deviance(gmsparse) deviance(gmdense)
@test_approx_eq coef(gmsparse) coef(gmdense)
@test_approx_eq vcov(gmsparse) vcov(gmdense)

## Prediction for GLMs
srand(1)
X = rand(10, 2)
Y = Vector{Float64}(logistic.(X * [3; -3])) # Julia 0.4 loses type information so Vector{Float64} can be dropped when we don't support 0.4

gm11 = fit(GeneralizedLinearModel, X, Y, Binomial())
@test_approx_eq predict(gm11) Y

newX = rand(5, 2)
newY = Vector{Float64}(logistic.(newX * coef(gm11))) # Julia 0.4 loses type information so Vector{Float64} can be dropped when we don't support 0.4
@test_approx_eq predict(gm11, newX) newY

off = rand(10)
newoff = rand(5)

@test_throws ArgumentError predict(gm11, newX, offset=newoff)

gm12 = fit(GeneralizedLinearModel, X, Y, Binomial(), offset=off)
@test_throws ArgumentError predict(gm12, newX)
@test_approx_eq predict(gm12, newX, offset=newoff) logistic.(newX * coef(gm12) .+ newoff)

## Prediction from DataFrames
d = convert(DataFrame, X)
d[:y] = Y

gm13 = fit(GeneralizedLinearModel, y ~ 0 + x1 + x2, d, Binomial())
@test predict(gm13) == predict(gm13, d[[:x1, :x2]])
@test predict(gm13) == predict(gm13, d)

newd = convert(DataFrame, newX)
predict(gm13, newd)

# Issue 118
@inferred nobs(lm(randn(10, 2), randn(10)))

# Issue 84
let
    X=[1 1; 2 4; 3 9]
    Xf=[1 1; 2 4; 3 9.]
    y = [2, 6, 12]
    yf = [2, 6, 12.]
    @test_approx_eq(lm(X, y).pp.beta0, ones(2))
    @test_approx_eq(lm(Xf, y).pp.beta0, ones(2))
    @test_approx_eq(lm(X, yf).pp.beta0, ones(2))
end