text:
function x0_volume_liquid(model::SAFTgammaMieModel,T,z)
    v_lb = lb_volume(model,z)
    return v_lb*2.0
end
function lb_volume(model::SAFTgammaMieModel, z = SA[1.0])
    vk  = model.groups.n_flattenedgroups
    seg = model.params.segment.values
    S   = model.params.shapefactor.values
    σᵢᵢ = model.params.sigma.diagvalues
    val = π/6*N_A*sum(z[i]*sum(vk[i][k]*seg[k]*S[k]*σᵢᵢ[k]^3 for k in @groups(i)) for i in @comps)
    return val
end

function T_scale(model::SAFTgammaMieModel,z=SA[1.0])
    return T_scale(model.vrmodel,z)
end

function p_scale(model::SAFTgammaMieModel,z=SA[1.0])
    V = zero(first(z))
    T = zero(first(z))
    σ̄3 = @f(σ3x)
    ϵ̄ = T_scale(model,z)
    val    = σ̄3*N_A/R̄/ϵ̄
    return 1/val
end

function a_res(model::SAFTgammaMieModel, V, T, z)
    _data = @f(data)
    return @f(a_hs,_data) + @f(a_disp,_data) + @f(a_chain,_data) + @f(a_assoc,_data)
end

function a_mono(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    return @f(a_hs,_data) + @f(a_disp,_data)
end

function data(model::SAFTgammaMieModel, V, T, z)
    _d_gc = @f(d_gc)
    _d = @f(d,_d_gc)
    ζi = @f(ζ0123,_d_gc)
    _ζ_X,σ3x = @f(ζ_X_σ3,_d_gc)
    _ρ_S = @f(ρ_S)
    _ζst = _ζst = σ3x*_ρ_S*π/6
    return (_d_gc,(_d,_ρ_S,ζi,_ζ_X,_ζst,σ3x))
end

function d_gc(model::SAFTgammaMieModel, V, T, z)
    _ϵ = model.params.epsilon.diagvalues
    _σ = model.params.sigma.diagvalues
    _λa = model.params.lambda_a.diagvalues
    _λr = model.params.lambda_r.diagvalues
    u = SAFTγMieconsts.u
    w = SAFTγMieconsts.w
    _0 = zero(T+first(z))
    n = length(_ϵ)
    _d_gc = zeros(typeof(_0),n)
    for k ∈ 1:n
        ϵ = _ϵ[k]
        σ = _σ[k]
        λa = _λa[k]
        λr = _λr[k]
        θ = Cλ(model.vrmodel,V,T,z,λa,λr)*ϵ/T
        di = _0
        λrinv = 1/λr
        λaλr = λa/λr
        for j ∈ 1:5
            di += w[j]*(θ/(θ+u[j]))^λrinv*(exp(θ*(1/(θ/(θ+u[j]))^λaλr-1))/(u[j]+θ)/λr)
        end
        _d_gc[k] = σ*(1-di)
    end
    return _d_gc
end

function d(model::SAFTgammaMieModel,V,T,z,_d_gc = @f(d_gc))
    _d = zeros(eltype(_d_gc),length(z))
    _0 = zero(eltype(_d))
    _zz = model.mie_zfractions.z.values
    for i ∈ @comps
        _z = _zz[i]
        m = length(_z)
        di = _0
        for k ∈ 1:m
            zk = _z[k]
            iszero(zk) && continue
            dk = _d_gc[k]
            di += zk*zk*dk^3
            for l ∈ 1:k - 1
                di += 0.25*zk*_z[l]*(dk + _d_gc[l])^3
            end
        end
        _d[i] = cbrt(di)
    end
    #d̂[i] = cbrt(∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*@f(d,k,l)^3 for l ∈ @groups) for k ∈ @groups))
    return _d
end

function ζ0123(model::SAFTgammaMieModel, V, T, z,_d_gc=@f(d_gc))
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    _0 = zero(V+T+first(z))
    ζ0,ζ1,ζ2,ζ3 = _0,_0,_0,_0
    _ϵ = model.params.epsilon.diagvalues
    σ = model.params.sigma.diagvalues
    _λr = model.params.lambda_r.diagvalues
    _λa = model.params.lambda_a.diagvalues
    u = SAFTVRMieconsts.u
    _mi = model.params.mixedsegment.values
    w = SAFTVRMieconsts.w
    for i ∈ @groups
        λa = _λa[i]
        λr = _λr[i]
        ϵ = _ϵ[i]
        di =_d_gc[i]
        mi = _mi[i]
        xS = dot(mi,z)*m̄inv
        ζ0 += xS
        ζ1 += xS*di
        ζ2 += xS*di*di
        ζ3 += xS*di*di*di
    end
    c = π/6*N_A*m̄/V
    ζ0,ζ1,ζ2,ζ3 = c*ζ0,c*ζ1,c*ζ2,c*ζ3
    return ζ0,ζ1,ζ2,ζ3
end

#return π/6*@f(ρ_S)*∑(@f(x_S,k)*@f(d,k)^m for k ∈ @groups)

function ρ_S(model::SAFTgammaMieModel, V, T, z)
    return ρ_S(model.vrmodel,V,T,z)
end

function ζ_X_σ3(model::SAFTgammaMieModel, V, T, z,_d = @f(d_gc))
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    σ = model.params.sigma.values
    _mi = model.params.mixedsegment.values
    ρS = N_A/V*m̄
    _ζ_X = zero(first(_d))
    kρS = ρS* π/6
    σ3_x = _ζ_X
    r1 =_ζ_X
    for i ∈ @groups
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        σ3_x += x_Si*x_Si*(σ[i,i]^3)
        _ζ_X += x_Si*x_Si*(_d[i])^3
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            σ3_x += 2*x_Si*x_Sj*(σ[i,j]^3)
            dij = 0.5*(_d[i] + _d[j])
            r1 = 2*x_Si*x_Sj*dij^3
            _ζ_X += r1
        end
    end

    #return π/6*@f(ρ_S)*∑(@f(x_S,i)*@f(x_S,j)*(@f(d,i)+@f(d,j))^3/8 for i ∈ comps for j ∈ comps)
    return kρS*_ζ_X,σ3_x
end

function σ3x(model::SAFTgammaMieModel, V, T, z)
    vrmodel = model.vrmodel
    m = vrmodel.params.segment.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    σ = model.params.sigma.values
    _mi = model.params.mixedsegment.values
    σ3_x =  zero(first(z))
    for i ∈ @groups
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        σ3_x += x_Si*x_Si*(σ[i,i]^3)
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            σ3_x += 2*x_Si*x_Sj*(σ[i,j]^3)
        end
    end
    return σ3_x
end

function  a_hs(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_hs(model.vrmodel,V,T,z,vrdata)
end

function a_chain(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_chain(model.vrmodel,V,T,z,vrdata)
end

function a_assoc(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    _,vrdata = _data
    return a_assoc(model.vrmodel,V,T,z,vrdata)
end

function a_disp(model::SAFTgammaMieModel, V, T, z,_data = @f(data))
    groups = @groups
    gcmodel = model
    model = gcmodel.vrmodel
    _d,vrdata = _data
    _,ρS,ζi,_ζ_X,_ζst,_ = vrdata
    comps = @comps
    l = length(comps)
    ∑z = ∑(z)
    m = model.params.segment.values
    _mi = gcmodel.params.mixedsegment.values
    _ϵ = gcmodel.params.epsilon.values
    _λr = gcmodel.params.lambda_r.values
    _λa = gcmodel.params.lambda_a.values
    _σ = gcmodel.params.sigma.values
    m̄ = dot(z, m)
    m̄inv = 1/m̄
    a₁ = zero(V+T+first(z))
    a₂ = a₁
    a₃ = a₁
    _ζst5 = _ζst^5
    _ζst8 = _ζst^8
    _KHS = @f(KHS,_ζ_X,ρS)
    for i ∈ groups
        j = i
        mi = _mi[i]
        x_Si = dot(mi,z)*m̄inv
        ϵ = _ϵ[i,j]
        λa = _λa[i,j]
        λr = _λr[i,j]
        σ = _σ[i,j]
        _C = @f(Cλ,λa,λr)
        dij = _d[i]
        x_0ij = σ/dij
        dij3 = dij^3
        x_0ij = σ/dij
        #calculations for a1 - diagonal
        aS_1_a = @f(aS_1,λa,_ζ_X)
        aS_1_r = @f(aS_1,λr,_ζ_X)
        B_a = @f(B,λa,x_0ij,_ζ_X)
        B_r = @f(B,λr,x_0ij,_ζ_X)
        a1_ij = (2*π*ϵ*dij3)*_C*ρS*
        (x_0ij^λa*(aS_1_a+B_a) - x_0ij^λr*(aS_1_r+B_r))

        #calculations for a2 - diagonal
        aS_1_2a = @f(aS_1,2*λa,_ζ_X)
        aS_1_2r = @f(aS_1,2*λr,_ζ_X)
        aS_1_ar = @f(aS_1,λa+λr,_ζ_X)
        B_2a = @f(B,2*λa,x_0ij,_ζ_X)
        B_2r = @f(B,2*λr,x_0ij,_ζ_X)
        B_ar = @f(B,λr+λa,x_0ij,_ζ_X)
        α = _C*(1/(λa-3)-1/(λr-3))
        f1,f2,f3,f4,f5,f6 = @f(f123456,α)
        _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
        a2_ij = π*_KHS*(1+_χ)*ρS*ϵ^2*dij3*_C^2 *
        (x_0ij^(2*λa)*(aS_1_2a+B_2a)
        - 2*x_0ij^(λa+λr)*(aS_1_ar+B_ar)
        + x_0ij^(2*λr)*(aS_1_2r+B_2r))

        #calculations for a3 - diagonal
        a3_ij = -ϵ^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
        #adding - diagonal
        a₁ += a1_ij*x_Si*x_Si
        a₂ += a2_ij*x_Si*x_Si
        a₃ += a3_ij*x_Si*x_Si
        for j ∈ 1:i-1
            mj = _mi[j]
            x_Sj = dot(mj,z)*m̄inv
            ϵ = _ϵ[i,j]
            λa = _λa[i,j]
            λr = _λr[i,j]
            σ = _σ[i,j]
            _C = @f(Cλ,λa,λr)
            dij = 0.5*(_d[i]+_d[j])
            x_0ij = σ/dij
            dij3 = dij^3
            x_0ij = σ/dij
            #calculations for a1
            a1_ij = (2*π*ϵ*dij3)*_C*ρS*
            (x_0ij^λa*(@f(aS_1,λa,_ζ_X)+@f(B,λa,x_0ij,_ζ_X)) - x_0ij^λr*(@f(aS_1,λr,_ζ_X)+@f(B,λr,x_0ij,_ζ_X)))

            #calculations for a2
            α = _C*(1/(λa-3)-1/(λr-3))
            f1,f2,f3,f4,f5,f6 = @f(f123456,α)
            _χ = f1*_ζst+f2*_ζst5+f3*_ζst8
            a2_ij = π*_KHS*(1+_χ)*ρS*ϵ^2*dij3*_C^2 *
            (x_0ij^(2*λa)*(@f(aS_1,2*λa,_ζ_X)+@f(B,2*λa,x_0ij,_ζ_X))
            - 2*x_0ij^(λa+λr)*(@f(aS_1,λa+λr,_ζ_X)+@f(B,λa+λr,x_0ij,_ζ_X))
            + x_0ij^(2*λr)*(@f(aS_1,2λr,_ζ_X)+@f(B,2*λr,x_0ij,_ζ_X)))

            #calculations for a3
            a3_ij = -ϵ^3*f4*_ζst * exp(f5*_ζst+f6*_ζst^2)
            #adding
            a₁ += 2*a1_ij*x_Si*x_Sj
            a₂ += 2*a2_ij*x_Si*x_Sj
            a₃ += 2*a3_ij*x_Si*x_Sj
        end
    end
    a₁ = a₁*m̄/T/∑z
    a₂ = a₂*m̄/(T*T)/∑z
    a₃ = a₃*m̄/(T*T*T)/∑z
    #@show (a₁,a₂,a₃)
    adisp =  a₁ + a₂ + a₃
    return adisp
end

const SAFTγMieconsts =(
    A = SA[ 0.81096    1.7888   -37.578   92.284;
          1.02050  -19.341    151.26  -463.50 ;
         -1.90570   22.845   -228.14   973.92 ;
          1.08850   -6.1962   106.98  -677.64  ],

    ϕ = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9  ],
        [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430 ],
        [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230  ],
        [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530 ],
        [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2 ],
        [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2 ],
        [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6  ]],

    u = [0.26356031971814109102031,1.41340305910651679221800,3.59642577104072208122300,7.08581000585883755692200,12.6408008442757826594300],
    w = [0.5217556105828086524759,0.3986668110831759274500,7.5942449681707595390e-2,3.6117586799220484545e-3,2.3369972385776227891e-5],

         c  = [0.0756425183020431	-0.128667137050961	 0.128350632316055	-0.0725321780970292	   0.0257782547511452  -0.00601170055221687	  0.000933363147191978  -9.55607377143667e-05  6.19576039900837e-06 -2.30466608213628e-07 3.74605718435540e-09
                0.134228218276565	    -0.182682168504886 	 0.0771662412959262	-0.000717458641164565 -0.00872427344283170	0.00297971836051287	 -0.000484863997651451	 4.35262491516424e-05 -2.07789181640066e-06	4.13749349344802e-08 0
               -0.565116428942893	     1.00930692226792   -0.660166945915607	 0.214492212294301	  -0.0388462990166792	0.00406016982985030	 -0.000239515566373142	 7.25488368831468e-06 -8.58904640281928e-08	0	                 0
               -0.387336382687019	    -0.211614570109503	 0.450442894490509	-0.176931752538907	   0.0317171522104923  -0.00291368915845693	  0.000130193710011706  -2.14505500786531e-06  0	                0	                 0
                2.13713180911797	    -2.02798460133021 	 0.336709255682693	 0.00118106507393722  -0.00600058423301506	0.000626343952584415 -2.03636395699819e-05	 0	                   0	                0	                 0
               -0.300527494795524	     2.89920714512243   -0.567134839686498	 0.0518085125423494	  -0.00239326776760414	4.15107362643844e-05  0	                     0	                   0	                0                    0
               -6.21028065719194	    -1.92883360342573	 0.284109761066570	-0.0157606767372364	   0.000368599073256615	0 	                  0	                     0	                   0	                0	                 0
                11.6083532818029	     0.742215544511197  -0.0823976531246117	 0.00186167650098254   0	                0	                  0	                     0	                   0	                0	                 0
               -10.2632535542427	    -0.125035689035085	 0.0114299144831867	 0	                   0	                0	                  0	                     0	                   0	                0	                 0
                4.65297446837297	    -0.00192518067137033 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0
               -0.867296219639940	     0	                 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0],
)


text:
println("https://docs.julialang.org/en/latest/manual/methods/")

f(x::Float64, y::Float64) = 2x + y
f(2., 3.)

f(x::Number, y::Number) = 2x + y
f(2, 3.)

methods(f)
methods(+)

# parametric methods
same_type(x::T, y::T) where {T} = true
same_type(x, y) = false
same_type(1, 2)
same_type(1., 2)

myappend(v::Vector{T}, x::T) where {T <: Number} = [v..., x]
myappend([1,2,3], 4)

function matmul(A::AbstractMatrix, B::AbstractMatrix)
    op = (ai, bi) -> ai * bi + ai * bi
    R = promote_op(op, eltype(A), eltype(B))
    output = similar(B, R, (size(A, 1), size(B, 2)))
    if size(A, 2) > 0
        for j in 1:size(B, 2)
            for i in 1:size(B, 1)
                ab::R = a[i, 1] * b[1, j]
                for k in 2:size(A, 2)
                    ab += a[i, k] * b[k, j]
                end
                output[i, j] = ab
            end
        end
    end
    return output
end

# function-like objects
struct Poly{R}
    coef::Vector{R}
end

function (p::Poly)(x)
    v = p.coef[end]
    for i = (length(p.coef)-1):-1:1
        v = v*x + p.coef[i]
    end
    return v
end

p = Poly([1,10,100])
println(p(3))

function emptyfunction
end

ext: jl

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
struct Demand
    itemId::String
    uom::String
    nodeId::String
    type::String
    quantity::Int128
end # struct

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
using Distributions, LinearAlgebra, Manifolds, PDMats, Random, Statistics, StatsBase, Test

@testset "DiffusionNormal" begin
    M = SpecialOrthogonal(3)
    @testset "constructors" begin
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        @inferred DiffusionNormal(M, μ, Σ)
        d = DiffusionNormal(M, μ, Σ)
        @test d isa DiffusionNormal
        @test d.Σ isa PDMats.PDMat
        @test d.μ == μ
        @test isapprox(M, d.e, identity(M, μ))
        @test d.direction === Manifolds.LeftAction()

        Σ = rand() * I
        @inferred DiffusionNormal(M, μ, Σ)
        d = DiffusionNormal(M, μ, Σ; direction=Manifolds.RightAction())
        @test d.Σ isa PDMats.ScalMat
        @test d.direction === Manifolds.RightAction()

        Σ = Diagonal(rand(3))
        @inferred DiffusionNormal(M, μ, Σ)
        d = DiffusionNormal(M, μ, Σ; e=Matrix{Float32}(Matrix(Diagonal(ones(3)))))
        @test d.Σ isa PDMats.PDiagMat
        @test d.e isa Matrix{Float32}
    end

    @testset "eltype" begin
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        @test eltype(d) === typeof(μ)
    end

    @testset "inversion" begin
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        dinv = inversion(d)
        @test dinv.μ == inv(M, μ)
        @test dinv.Σ == d.Σ
        @test dinv.basis === d.basis
        @test dinv.e == d.e
        @test dinv.direction === Manifolds.RightAction()
    end

    @testset "rand! SpecialOrthogonal" begin
        M = SpecialOrthogonal(3)
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        q = Matrix{Float64}(undef, 3, 3)
        d = DiffusionNormal(M, μ, Σ)
        @test rand!(d, q) === q
        @test is_manifold_point(M, q, true)
    end

    @testset "rand SpecialOrthogonal" begin
        M = SpecialOrthogonal(3)
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        @inferred rand(d)
        q = rand(d)
        @test is_manifold_point(M, q, true)

        ps = rand(d, 2)
        @test typeof(ps) == Vector{typeof(μ)}
        @test length(ps) == 2
        @test is_manifold_point(M, ps[1], true)
        @test is_manifold_point(M, ps[2], true)
    end

    @testset "insupport" begin
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        p = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        @test insupport(d, p)
        @test !insupport(d, randn(3))
        @test !insupport(d, randn(4, 4))
        @test !insupport(d, zeros(3, 3))
        @test !insupport(d, Matrix(Diagonal([-1.0, 1.0, 1.0])))
    end

    @testset "mode" begin
        μ = group_exp(M, hat(M, Matrix(Diagonal(ones(3))), randn(3)))
        Σ = exp(Symmetric(randn(3, 3)))
        d = DiffusionNormal(M, μ, Σ)
        @test mode(d) == μ
    end
end


text:
"""\
 C version: Dmitri Laikov
 F77 version: Christoph van Wuellen, http://www.ccl.net
 Python version: Richard P. Muller, 2002.
 Julia version: John C. Snyder, 2015.

 This subroutine is part of a set of subroutines that generate
 Lebedev grids [1-6] for integration on a sphere. The original
 C-code [1] was kindly provided by Dr. Dmitri N. Laikov and
 translated into fortran by Dr. Christoph van Wuellen.
 This subroutine was translated from C to fortran77 by hand.

 Users of this code are asked to include reference [1] in their
 publications, and in the user- and programmers-manuals
 describing their codes.

   [1] V.I. Lebedev, and D.N. Laikov
       'A quadrature formula for the sphere of the 131st
        algebraic order of accuracy'
       Doklady Mathematics, Vol. 59, No. 3, 1999, pp. 477-481.

   [2] V.I. Lebedev
       'A quadrature formula for the sphere of 59th algebraic
        order of accuracy'
       Russian Acad. Sci. Dokl. Math., Vol. 50, 1995, pp. 283-286.

   [3] V.I. Lebedev, and A.L. Skorokhodov
       'Quadrature formulas of orders 41, 47, and 53 for the sphere'
       Russian Acad. Sci. Dokl. Math., Vol. 45, 1992, pp. 587-592.

   [4] V.I. Lebedev
       'Spherical quadrature formulas exact to orders 25-29'
       Siberian Mathematical Journal, Vol. 18, 1977, pp. 99-107.

   [5] V.I. Lebedev
       'Quadratures on a sphere'
       Computational Mathematics and Mathematical Physics, Vol. 16,
       1976, pp. 10-24.

   [6] V.I. Lebedev
       'Values of the nodes and weights of ninth to seventeenth
        order Gauss-Markov quadrature formulae invariant under the
        octahedron group with inversion'
       Computational Mathematics and Mathematical Physics, Vol. 15,
       1975, pp. 44-51.
"""

function genOh_a00(v)
    "(0,0,a) etc. (6 points)"
    a=1.0
    hcat([a,0,0,v],[-a,0,0,v],[0,a,0,v],[0,-a,0,v],[0,0,a,v],[0,0,-a,v])'
end

function genOh_aa0(v)
    "(0,a,a) etc, a=1/sqrt(2) (12 points)"
    a=sqrt(0.5)
    hcat([0,a,a,v],[0,-a,a,v],[0,a,-a,v],[0,-a,-a,v],
         [a,0,a,v],[-a,0,a,v],[a,0,-a,v],[-a,0,-a,v],
         [a,a,0,v],[-a,a,0,v],[a,-a,0,v],[-a,-a,0,v])'
 end

function genOh_aaa(v)
    "(a,a,a) etc, a=1/sqrt(3) (8 points)"
    a = sqrt(1./3.)
    hcat([a,a,a,v],[-a,a,a,v],[a,-a,a,v],[-a,-a,a,v],
         [a,a,-a,v],[-a,a,-a,v],[a,-a,-a,v],[-a,-a,-a,v])'
 end

function genOh_aab(v,a)
    "(a,a,b) etc, b=sqrt(1-2 a^2), a input (24 points)"
    b = sqrt(1.0 - 2.0*a*a)
    hcat([a,a,b,v],[-a,a,b,v],[a,-a,b,v],[-a,-a,b,v],
         [a,a,-b,v],[-a,a,-b,v],[a,-a,-b,v],[-a,-a,-b,v],
         [a,b,a,v],[-a,b,a,v],[a,-b,a,v],[-a,-b,a,v],
         [a,b,-a,v],[-a,b,-a,v],[a,-b,-a,v],[-a,-b,-a,v],
         [b,a,a,v],[-b,a,a,v],[b,-a,a,v],[-b,-a,a,v],
         [b,a,-a,v],[-b,a,-a,v],[b,-a,-a,v],[-b,-a,-a,v])'
 end

function genOh_ab0(v,a)
    "(a,b,0) etc, b=sqrt(1-a^2), a input (24 points)"
    b=sqrt(1.0-a*a)
    hcat([a,b,0,v],[-a,b,0,v],[a,-b,0,v],[-a,-b,0,v],
         [b,a,0,v],[-b,a,0,v],[b,-a,0,v],[-b,-a,0,v],
         [a,0,b,v],[-a,0,b,v],[a,0,-b,v],[-a,0,-b,v],
         [b,0,a,v],[-b,0,a,v],[b,0,-a,v],[-b,0,-a,v],
         [0,a,b,v],[0,-a,b,v],[0,a,-b,v],[0,-a,-b,v],
         [0,b,a,v],[0,-b,a,v],[0,b,-a,v],[0,-b,-a,v])'
end

function genOh_abc(v,a,b)
    "(a,b,c) etc, c=sqrt(1-a^2-b^2), a,b input  (48 points)"
    c=sqrt(1.0 - a*a - b*b)
    hcat([a,b,c,v],[-a,b,c,v],[a,-b,c,v],[-a,-b,c,v],
         [a,b,-c,v],[-a,b,-c,v],[a,-b,-c,v],[-a,-b,-c,v],
         [a,c,b,v],[-a,c,b,v],[a,-c,b,v],[-a,-c,b,v],
         [a,c,-b,v],[-a,c,-b,v],[a,-c,-b,v],[-a,-c,-b,v],
         [b,a,c,v],[-b,a,c,v],[b,-a,c,v],[-b,-a,c,v],
         [b,a,-c,v],[-b,a,-c,v],[b,-a,-c,v],[-b,-a,-c,v],
         [b,c,a,v],[-b,c,a,v],[b,-c,a,v],[-b,-c,a,v],
         [b,c,-a,v],[-b,c,-a,v],[b,-c,-a,v],[-b,-c,-a,v],
         [c,a,b,v],[-c,a,b,v],[c,-a,b,v],[-c,-a,b,v],
         [c,a,-b,v],[-c,a,-b,v],[c,-a,-b,v],[-c,-a,-b,v],
         [c,b,a,v],[-c,b,a,v],[c,-b,a,v],[-c,-b,a,v],
         [c,b,-a,v],[-c,b,-a,v],[c,-b,-a,v],[-c,-b,-a,v])'
end

function leb6()
    genOh_a00(0.1666666666666667)
end

function leb14()
    [genOh_a00(0.06666666666666667);
     genOh_aaa(0.07500000000000000)]
end

function leb26()
    [genOh_a00(0.04761904761904762);
     genOh_aa0(0.03809523809523810);
     genOh_aaa(0.03214285714285714)]
end

function leb38()
    [genOh_a00(0.009523809523809524);
     genOh_aaa(0.3214285714285714E-1);
     genOh_ab0(0.2857142857142857E-1,0.4597008433809831E+0)]
end

function leb50()
    [genOh_a00(0.1269841269841270E-1);
     genOh_aa0(0.2257495590828924E-1);
     genOh_aaa(0.2109375000000000E-1);
     genOh_aab(0.2017333553791887E-1,0.3015113445777636E+0)]
end

function leb74()
    [genOh_a00(0.5130671797338464E-3);
     genOh_aa0(0.1660406956574204E-1);
     genOh_aaa(-0.2958603896103896E-1);
     genOh_aab(0.2657620708215946E-1,0.4803844614152614E+0);
     genOh_ab0(0.1652217099371571E-1,0.3207726489807764E+0)]
end

function leb86()
    [genOh_a00(0.1154401154401154E-1);
     genOh_aaa(0.1194390908585628E-1);
     genOh_aab(0.1111055571060340E-1,0.3696028464541502E+0);
     genOh_aab(0.1187650129453714E-1,0.6943540066026664E+0);
     genOh_ab0(0.1181230374690448E-1,0.3742430390903412E+0)]
end

function leb110()
    [genOh_a00(0.3828270494937162E-2);
     genOh_aaa(0.9793737512487512E-2);
     genOh_aab(0.8211737283191111E-2,0.1851156353447362E+0);
     genOh_aab(0.9942814891178103E-2,0.6904210483822922E+0);
     genOh_aab(0.9595471336070963E-2,0.3956894730559419E+0);
     genOh_ab0(0.9694996361663028E-2,0.4783690288121502E+0)]
end

function leb146()
    [genOh_a00(0.5996313688621381E-3);
     genOh_aa0(0.7372999718620756E-2);
     genOh_aaa(0.7210515360144488E-2);
     genOh_aab(0.7116355493117555E-2,0.6764410400114264E+0);
     genOh_aab(0.6753829486314477E-2,0.4174961227965453E+0);
     genOh_aab(0.7574394159054034E-2,0.1574676672039082E+0);
     genOh_abc(0.6991087353303262E-2,0.1403553811713183E+0,
                     0.4493328323269557E+0)]
end

function leb170()
    [genOh_a00(0.5544842902037365E-2);
     genOh_aa0(0.6071332770670752E-2);
     genOh_aaa(0.6383674773515093E-2);
     genOh_aab(0.5183387587747790E-2,0.2551252621114134E+0);
     genOh_aab(0.6317929009813725E-2,0.6743601460362766E+0);
     genOh_aab(0.6201670006589077E-2,0.4318910696719410E+0);
     genOh_ab0(0.5477143385137348E-2,0.2613931360335988E+0);
     genOh_abc(0.5968383987681156E-2,0.4990453161796037E+0,
                     0.1446630744325115E+0)]
end

function leb194()
    [genOh_a00(0.1782340447244611E-2);
     genOh_aa0(0.5716905949977102E-2);
     genOh_aaa(0.5573383178848738E-2);
     genOh_aab(0.5608704082587997E-2,0.6712973442695226E+0);
     genOh_aab(0.5158237711805383E-2,0.2892465627575439E+0);
     genOh_aab(0.5518771467273614E-2,0.4446933178717437E+0);
     genOh_aab(0.4106777028169394E-2,0.1299335447650067E+0);
     genOh_ab0(0.5051846064614808E-2,0.3457702197611283E+0);
     genOh_abc(0.5530248916233094E-2,0.1590417105383530E+0,
                     0.8360360154824589E+0)]
end

leb = Dict(
    6   => leb6,
    14  => leb14,
    26  => leb26,
    38  => leb38,
    50  => leb50,
    74  => leb74,
    86  => leb86,
    110 => leb110,
    146 => leb146,
    170 => leb170,
    194 => leb194
)

# generate lebedev grids
function generate_lebedev()

    function lebgrid(n::Integer)
        data = leb[n]()'
        w = vec(4pi*data[end,:])
        p = data[1:3,:]
        p = Point3{Float64}[tuple(p[:,i]...) for i in 1:size(p,2)]
        AngularGrid(p,w)
    end

    lebedev = Dict{Int,AngularGridf}([k=>lebgrid(k) for k in keys(leb)])

    open(joinpath(Pkg.dir("Element"), "deps/lebedev.jls"), "w") do io serialize(io, lebedev) end

end

max_stars_repo_name: johncsnyder/Element


text:
# views that provide efficient construction and element access

abstract AbstractUnsafeView{T, N} <: AbstractArray{T, N}

typealias ContiguousArray{T<:Number, N} Union(Array{T, N}, AbstractUnsafeView{T, N})
typealias ContiguousVector{T<:Number} ContiguousArray{T, 1}
typealias ContiguousMatrix{T<:Number} ContiguousArray{T, 2}
typealias ContiguousCube{T<:Number} ContiguousArray{T, 3}

typealias ContiguousArrOrNum{T<:Number} Union(ContiguousArray{T}, T)

typealias ContiguousRealArray{T<:Real,N} ContiguousArray{T,N}
typealias ContiguousRealArrOrNum{T<:Real} Union(ContiguousArray{T}, T)

similar(a::AbstractUnsafeView, T::DataType, shape::NTuple) = Array(T, shape)
similar{T}(a::AbstractUnsafeView{T}) = Array(T, size(a))

function copy{T}(a::AbstractUnsafeView{T})
    r = similar(a)
    if isbits(T)
        unsafe_copy!(pointer(r), pointer(a), length(a))
    else
        for i = 1 : length(a)
            r[i] = a[i]
        end
    end
    r
end

immutable UnsafeVectorView{T} <: AbstractUnsafeView{T, 1}
    ptr::Ptr{T}
    len::Int
end

pointer(a::UnsafeVectorView) = a.ptr
pointer{T}(a::UnsafeVectorView{T}, i::Int) = a.ptr + (i - 1) * sizeof(T)
size(a::UnsafeVectorView) = (a.len,)
size(a::UnsafeVectorView, d::Int) = d == 1 ? a.len : 1
length(a::UnsafeVectorView) = a.len

getindex(a::UnsafeVectorView, i::Int) = unsafe_load(a.ptr, i)
setindex!(a::UnsafeVectorView, v, i::Int) = unsafe_store!(a.ptr, v, i)

immutable UnsafeMatrixView{T} <: AbstractUnsafeView{T, 2}
    ptr::Ptr{T}
    dim1::Int
    dim2::Int
    len::Int

    UnsafeMatrixView(ptr::Ptr{T}, d1::Int, d2::Int) = new(ptr, d1, d2, d1 * d2)
end

pointer(a::UnsafeMatrixView) = a.ptr
pointer{T}(a::UnsafeMatrixView{T}, i::Int) = a.ptr + (i - 1) * sizeof(T)
size(a::UnsafeMatrixView) = (a.dim1, a.dim2)
size(a::UnsafeMatrixView, d::Int) = d == 1 ? a.dim1 : d == 2 ? a.dim2 : 1
length(a::UnsafeMatrixView) = a.len

getindex(a::UnsafeMatrixView, i::Int) = unsafe_load(a.ptr, i)
getindex(a::UnsafeMatrixView, i::Int, j::Int) = unsafe_load(a.ptr, i + (j-1) * a.dim1)

setindex!(a::UnsafeMatrixView, v, i::Int) = unsafe_store!(a.ptr, v, i)
setindex!(a::UnsafeMatrixView, v, i::Int, j::Int) = unsafe_store!(a.ptr, v, i + (j-1) * a.dim1)

immutable UnsafeCubeView{T} <: AbstractUnsafeView{T, 3}
    ptr::Ptr{T}
    dim1::Int
    dim2::Int
    dim3::Int
    len::Int
    plen::Int

    UnsafeCubeView(ptr::Ptr{T}, d1::Int, d2::Int, d3::Int) = new(ptr, d1, d2, d3, d1 * d2 * d3, d1 * d2)
end

pointer(a::UnsafeCubeView) = a.ptr
pointer{T}(a::UnsafeCubeView{T}, i::Int) = a.ptr + (i - 1) * sizeof(T)
size(a::UnsafeCubeView) = (a.dim1, a.dim2, a.dim3)
size(a::UnsafeCubeView, d::Int) = d == 1 ? a.dim1 : d == 2 ? a.dim2 : d == 3 ? a.dim3 : 1
length(a::UnsafeCubeView) = a.len

getindex(a::UnsafeCubeView, i::Int) = unsafe_load(a.ptr, i)
getindex(a::UnsafeCubeView, i::Int, j::Int) = unsafe_load(a.ptr, i + (j-1) * a.dim1)
getindex(a::UnsafeCubeView, i::Int, j::Int, k::Int) = unsafe_load(a.ptr, i + (j-1) * a.dim1 + (k-1) * a.plen)

setindex!(a::UnsafeCubeView, v, i::Int) = unsafe_store!(a.ptr, v, i)
setindex!(a::UnsafeCubeView, v, i::Int, j::Int) = unsafe_store!(a.ptr, v, i + (j-1) * a.dim1)
setindex!(a::UnsafeCubeView, v, i::Int, j::Int, k::Int) = unsafe_store!(a.ptr, v, i + (j-1) * a.dim1 + (k-1) * a.plen)

# calculate offset of an element

offset(a::ContiguousArray, i::Int) = (i - 1)
offset(a::ContiguousArray, i::Int, j::Int) = (i - 1) + size(a, 1) * (j - 1)
offset(a::ContiguousArray, i::Int, j::Int, k::Int) = (i - 1) + size(a, 1) * ((j - 1) + size(a, 2) * (k - 1))

offset(a::UnsafeMatrixView, i::Int, j::Int) = (i - 1) + a.dim1 * (j - 1)
offset(a::UnsafeCubeView, i::Int, j::Int, k::Int) = (i - 1) + a.dim1 * ((j - 1) + a.dim2 * (k - 1))

# view constructors

offset_view{T}(a::ContiguousArray{T}, offset::Int, d1::Int) = UnsafeVectorView(pointer(a) + sizeof(T) * offset, d1)
offset_view{T}(a::ContiguousArray{T}, offset::Int, d1::Int, d2::Int) = UnsafeMatrixView{T}(pointer(a) + sizeof(T) * offset, d1, d2)

function offset_view{T}(a::ContiguousArray{T}, offset::Int, d1::Int, d2::Int, d3::Int)
    UnsafeCubeView{T}(pointer(a) + sizeof(T) * offset, d1, d2, d3)
end

unsafe_view(v) = v  # fallback

# 1D

unsafe_view(a::ContiguousVector) = UnsafeVectorView(pointer(a), length(a))

unsafe_view(a::ContiguousArray, r::Colon) = UnsafeVectorView(pointer(a), length(a))
unsafe_view(a::ContiguousArray, r::Colon, j::Int) = offset_view(a, offset(a, 1, j), size(a, 1))
unsafe_view(a::ContiguousArray, r::Colon, j::Int, k::Int) = offset_view(a, offset(a, 1, j, k), size(a, 1))

unsafe_view(a::ContiguousArray, r::Range1{Int}) = offset_view(a, r[1] - 1, length(r))
unsafe_view(a::ContiguousArray, r::Range1{Int}, j::Int) = offset_view(a, offset(a, r[1], j), length(r))
unsafe_view(a::ContiguousArray, r::Range1{Int}, j::Int, k::Int) = offset_view(a, offset(a, r[1], j, k), length(r))

# 2D

unsafe_view{T}(a::ContiguousMatrix{T}) = UnsafeMatrixView{T}(pointer(a), size(a, 1), size(a, 2))

unsafe_view{T}(a::ContiguousVector{T}, r1::Colon, r2::Colon) = UnsafeMatrixView{T}(pointer(a), length(a), 1)
unsafe_view{T}(a::ContiguousMatrix{T}, r1::Colon, r2::Colon) = UnsafeMatrixView{T}(pointer(a), size(a, 1), size(a, 2))

function unsafe_view{T}(a::ContiguousArray{T}, r1::Colon, r2::Colon)
    siz = size(a)
    UnsafeMatrixView{T}(pointer(a), siz[1], succ_length(siz, 1))
end

unsafe_view(a::ContiguousArray, r::Colon, r2::Range1{Int}) = offset_view(a, offset(a, 1, r2[1]), size(a, 1), length(r2))
unsafe_view(a::ContiguousArray, r::Colon, r2::Range1{Int}, k::Int) = offset_view(a, offset(a, 1, r2[1], k), size(a, 1), length(r2))
unsafe_view(a::ContiguousArray, r::Colon, r2::Colon, k::Int) = offset_view(a, offset(a, 1, 1, k), size(a, 1), size(a, 2))

# 3D

unsafe_view{T}(a::ContiguousCube{T}) = UnsafeCubeView{T}(pointer(a), size(a, 1), size(a, 2), size(a, 3))

unsafe_view{T}(a::ContiguousVector{T}, r::Colon, r2::Colon, r3::Colon) = UnsafeCubeView{T}(pointer(a), size(a, 1), 1, 1)
unsafe_view{T}(a::ContiguousMatrix{T}, r::Colon, r2::Colon, r3::Colon) = UnsafeCubeView{T}(pointer(a), size(a, 1), size(a, 2), 1)

function unsafe_view{T}(a::ContiguousArray{T}, r::Colon, r2::Colon, r3::Colon)
    siz = size(a)
    UnsafeCubeView{T}(pointer(a), siz[1], siz[2], succ_length(siz, 2))
end

function unsafe_view(a::ContiguousArray, r::Colon, r2::Colon, r3::Range1{Int})
    offset_view(a, offset(a, 1, 1, r3[1]), size(a, 1), size(a, 2), length(r3))
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module ZumkellerNumbers

using NumberTheory, IterTools, CombinationsIterator

export ModuleZumkellerNumbers
export isZumkeller, is083207, I083207, F083207, L083207, V083207

"""
A Zumkeller number ``n`` is an integer whose divisors can be partitioned into two disjoint sets whose sums are both ``σ(n)/2``.

* isZumkeller, is083207, I083207, F083207, L083207, V083207
"""
const ModuleZumkellerNumbers = ""

"""
Is ``n`` a Zumkeller number? A Zumkeller number ``n`` is an integer whose divisors can be partitioned into two disjoint sets whose sums are both ``σ(n)/2``.
"""
function isZumkeller(n::Int)
    n == 0 && return false
    s = σ(n)
    ((s % 2 ≠ 0) || (s < 2n)) && return false
    S = s >> 1 - n
    D = [d for d in Divisors(n) if d ≤ S]
    D == [] && return true
    for c in Combinations(D)
        S == sum(c) && return true
    end
    return false
end

"""
Is ``n`` a Zumkeller number?
"""
is083207(n) = isZumkeller(n)

"""
Iterate over the first ``n`` Zumkeller numbers.
"""
I083207(n) = Iterators.take(Iterators.filter(isZumkeller, Iterators.countfrom(1)), n)

"""
Iterate over the Zumkeller numbers ``z`` which are below ``n, (1 ≤ z ≤ n)``.
"""
F083207(n) = Iterators.filter(isZumkeller, 1:n)

"""
List the first ``n`` Zumkeller numbers.
"""
L083207(n) = collect(I083207(n))

"""
Return the ``n``-th Zumkeller number.
"""
V083207(n) = nth(I083207(n), n)

#START-TEST-########################################################

using Test

function test()
    @testset "Zumkeller" begin
        @test isZumkeller(17000) == true
        @test isZumkeller(27472) == true
        @test isZumkeller(29062) == false
        @test isZumkeller(43464) == true
    end
end

function demo()
    println(L083207(10))
    for a in I083207(30) print(a, ", ") end; println("...")
    for a in F083207(30) print(a, ", ") end; println("...")

    for n in 20:30
        println(n, " ↦ ",  isZumkeller(n))
    end

    for (index, value) in enumerate(I083207(10))
        println("$index -> $value")
    end

    println(V083207(10))

    for n in 1:6
        println(n, " ↦ ", L083207(n))
    end
end

"""
@time (for n in 1:2000 isZumkeller(n) end) ::
    0.311896 seconds (1.86 M allocations: 63.875 MiB, 31.28% gc time)
@time L083207(500) ::
    0.295718 seconds (2.21 M allocations: 75.920 MiB, 40.89% gc time)
"""
function perf()
    GC.gc()
    @time (for n in 1:2000 isZumkeller(n) end)
    @time L083207(500)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
[6, 12, 20, 24, 28, 30, 40, 42, 48, 54]
6, 12, 20, 24, 28, 30, 40, 42, 48, 54, 56, 60, 66, 70, 78, 80, 84,
88, 90, 96, 102, 104, 108, 112, 114, 120, 126, 132, 138, 140, ...
6, 12, 20, 24, 28, 30, ...

20 ↦ true
21 ↦ false
22 ↦ false
23 ↦ false
24 ↦ true
25 ↦ false
26 ↦ false
27 ↦ false
28 ↦ true
29 ↦ false
30 ↦ true

1 -> 6
2 -> 12
3 -> 20
4 -> 24
5 -> 28
6 -> 30
7 -> 40
8 -> 42
9 -> 48
10 -> 54
54

1 ↦ [6]
2 ↦ [6, 12]
3 ↦ [6, 12, 20]
4 ↦ [6, 12, 20, 24]
5 ↦ [6, 12, 20, 24, 28]
6 ↦ [6, 12, 20, 24, 28, 30]

  0.311551 seconds (1.86 M allocations: 63.820 MiB, 60.62% gc time)
  0.294909 seconds (2.21 M allocations: 75.858 MiB, 41.31% gc time)
=#


text:
using JuMP
using MathOptInterface

optimizer_solver = optimizer_with_attributes(
    ActiveSetMethods.Optimizer,
    "external_optimizer" => GLPK.Optimizer,
    "OutputFlag" => 0,
)

model = Model(optimizer_solver)

@variable(model, X);
@variable(model, Y);
@objective(model, Min, X^2 + X);
@NLconstraint(model, X^2 - X == 2);
@NLconstraint(model, X * Y == 1);
@NLconstraint(model, X * Y >= 0);
@constraint(model, X >= -2);

JuMP.optimize!(model);

xsol = JuMP.value.(X)
ysol = JuMP.value.(Y)
status = termination_status(model)

max_stars_repo_name: exanauts/ActiveSetMethods

max_issues_repo_name: exanauts/ActiveSetMethods

max_forks_repo_name: exanauts/ActiveSetMethods


text:
using Fst

text = "0 1 a b 1.0\n1 2 d c 1.0\n2 1.0"

create_pdf(read_wfst(text), "ioread.pdf")


text:
# special_sequences.jl
# A collection of special sequences.

"The ZeroSequence is a sequence of all zeros."
immutable ZeroSequence{T} <: Sequence{T}
end

ZeroSequence() = ZeroSequence{Float64}()

getindex(s::ZeroSequence, k) = zero(eltype(s))

ztransform(s::ZeroSequence, z) = zero(eltype(s))

"The DiracSequence is a sequence of all zeros, except for `h[0] = 1`."
immutable DiracSequence{T} <: Sequence{T}
end

DiracSequence() = DiracSequence{Float64}()

getindex(s::DiracSequence, k) = k == 0 ? one(eltype(s)) : zero(eltype(s))

ztransform(s::DiracSequence, z) = one(eltype(s))


text:
using BenchmarkTools
using Hydroelastics
using LinearAlgebra
using StaticArrays
import Hydroelastics: intersect_halfplanes, intersect_halfplanes_slow, HalfPlane, Point

macro run_bench(name, expr, args...)
    quote
        printstyled("Benchmark: ", $name, "\n"; color = :yellow)
        display(@benchmark esc($expr) $(args...))
        println("\n\n")
    end
end

tet1 = Object(
    Mesh(
        [
            0.0 1.0 0.0 0.0
            0.0 0.0 1.0 0.0
            0.0 0.0 0.0 1.0
        ],
        reshape([1, 2, 3, 4], :, 1),
        [0.0, 0.0, 0.0, 1.0],
    ),
)
tet2 = Object(
    Mesh(
        [
            0.0 1.0 0.0 0.0
            0.0 0.0 1.0 0.0
            1.8 1.8 1.8 0.8
        ],
        reshape([1, 2, 3, 4], :, 1),
        [0.0, 0.0, 0.0, 1.0],
    ),
)

@run_bench "Tet force" begin
    Hydroelastics.tet_force(tet1, tet2, 1, 1)
end

sphere = make_icosphere(2)
@run_bench "Icosphere force" begin
    local object1 = translate(sphere, @SVector [0.031, -0.5, 0.052])
    local object2 = translate(sphere, @SVector [0, 0.5, 0])
    compute_force(object1, object2)
end

halfplanes = [
    HalfPlane(Point(-0.16469227056401645, -0.2664776914569725), Point(-0.8506508083520399, 0.5257311121191335)),
    HalfPlane(Point(0.028287259696359788, -0.040547756753354335), Point(0.8201447309019166, 0.5721561153861965)),
    HalfPlane(Point(-0.16620458460952442, 0.09547326678562677), Point(-0.49810107091297723, -0.8671189786617204)),
    HalfPlane(Point(0.2970629430715192, -0.23999299562154927), Point(0.6284282897093981, 0.77786752387082)),
    HalfPlane(Point(-0.05197402131154934, 0.06002998491806123), Point(-0.7560128937297079, -0.6545567236797993)),
    HalfPlane(Point(0.1879473775300143, 0.30410524493997143), Point(0.85065080835204, -0.5257311121191336)),
    HalfPlane(Point(0.1619918586413753, -0.0851641599977757), Point(0.46534112719498627, 0.8851314226378474)),
    HalfPlane(Point(0.2970629430715192, -0.23999299562154938), Point(0.6284282897093982, 0.7778675238708199)),
]

@run_bench "Original halfplane algorithm using Polyhedra" begin
    intersect_halfplanes_slow(halfplanes)
end

@run_bench "New halfplane algorithm" begin
    intersect_halfplanes(halfplanes)
end

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
"""
    interpolate_and_bin(
            left_bin_edges::AbstractRange{T2}, xs, ys,
            intp::Linear, intp_grid::AbstractRange{T1},
             extrapolation_bc::Union{<:Real, Interpolations.BoundaryCondition}
            f::Function, args...; kwargs...) where {T1, T2}

Linearly interpolate the `xs` and `ys` to `intp_grid`, then bin the interpolated data
on the grid defined by the `N`-element `left_bin_edges` grid. After binning, apply
`f` with arguments `args` and keyword arguments `kwargs` element-wise to each of the bins.
Returns a vector of `N-1` bin summaries.

## Details

Binning is performed *after* interpolating `xs` and `ys` to `intp_grid`. Binning is done
by distributing the `intp_grid` values among the bins defined by `left_bin_edges`. That is,
if `intp_grid[i]` falls in the `k`-th bin, then `intp(xs)[i]` is assigned to the `k`-th bin.
If `left_bin_edges` contain `N` grid points, then there will be `N - 1` bins in total.

If `xs` contain `NaN` values, then these are filled internally by linear interpolation
*before* interpolating to the provided `intp_grid`.
"""
function interpolate_and_bin(f::Function,
            left_bin_edges::AbstractRange{T2}, xs, ys,
            intp::Linear, intp_grid::AbstractRange{T1}, extrapolation_bc::Union{<:Real, Interpolations.BoundaryCondition},
            args...; kwargs...) where {T1, T2}

    if length(xs) != length(ys)
        msg = "`length(xs)` not equal to `length(ys)`"
        throw(ArgumentError(msg))
    end

    #=
    Handle `NaN` values initially present in `xs` by linear interpolation and
    extrapolating using the provided boundary condition. Only internal `NaN`s are
    interpolated; if there are `NaN`s present at the edges of `xs` initially, then
    the corresponding `ys` are ignored.
    =#
    if any(isnan.(xs))
        println("There were `NaN`s in `xs`. Interpolating them by introducing an internal grid in `xs`.")
        xs = fill_nans(xs, Linear())
    end

    if any(isnan.(ys))
        println("There were `NaN`s in `ys`. Interpolating them by introducing an internal grid in `ys`.")
        ys = fill_nans(ys, Linear())
    end

    isnotnans = findall(isfinite.(xs) .& isfinite.(ys))
    tmp_intp = LinearInterpolation(xs[isnotnans], ys[isnotnans],
        extrapolation_bc = extrapolation_bc)

    # Interpolate `ys` to `intp_grid`
    intp_ys = tmp_intp(intp_grid)

    # Each bin is represented by (potentially differently sized) vectors of type `eltype(ys)`
    YT = eltype(intp_ys)
    n_bins = length(left_bin_edges) - 1
    bin_vals = [Array{YT}(undef, 0) for i = 1:n_bins]

    # Assign values to bins. Can be done faster if `xs` is guaranteed to be ordered.
    g_min = minimum(left_bin_edges)
    g_step = step(left_bin_edges)

    for (x, y) in zip(intp_grid, intp_ys)
        arr_idx = ceil(Int, (x - g_min)/g_step)
        # Check that the value falls inside the grid; if not, ignore it
        if 0 < arr_idx <= n_bins
            push!(bin_vals[arr_idx], y)
        end
    end

    # Get bin summaries by applying function to non-empty bins.
    bin_summaries = fill(NaN, n_bins)
    inds_nonempty_bins = findall(length.(bin_vals) .> 0)
    bin_summaries[inds_nonempty_bins] = f.(bin_vals[inds_nonempty_bins], args...; kwargs...)

    return bin_summaries
end

export interpolate_and_bin


text:
module LAPACK2

    using LinearAlgebra
    using LinearAlgebra: BlasInt, chkstride1, LAPACKException
    using LinearAlgebra.BLAS: @blasfunc
    using LinearAlgebra.LAPACK: chkdiag, chkside, chkuplo

    ## Standard QR/QL
    function steqr!(compz::Char,
                    d::StridedVector{Float64},
                    e::StridedVector{Float64},
                    Z::StridedMatrix{Float64} = compz == 'N' ?
                        Matrix{Float64}(undef, 0, 0) :
                        Matrix{Float64}(undef, length(d), length(d)),
                    work::StridedVector{Float64} = compz == 'N' ?
                        Vector{Float64}(undef, 0) :
                        Vector{Float64}(undef, max(1, 2*length(d) - 2)))

        # Extract sizes
        n = length(d)
        ldz = stride(Z, 2)

        # Checks
        length(e) + 1 == n || throw(DimensionMismatch(""))
        if compz == 'N'
        elseif compz in ('V', 'I')
            chkstride1(Z)
            n <= size(Z, 1) || throw(DimensionMismatch("Z matrix has too few rows"))
            n <= size(Z, 2) || throw(DimensionMismatch("Z matrix has too few colums"))
        else
            throw(ArgumentError("compz must be either 'N', 'V', or 'I'"))
        end

        # Allocations
        info = Vector{BlasInt}(undef, 1)

        ccall((@blasfunc("dsteqr_"), Base.liblapack_name),Cvoid,
            (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
             Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ptr{BlasInt}),
             compz, n, d, e,
             Z, max(1, ldz), work, info)
        info[1] == 0 || throw(LAPACKException(info[1]))

        return d, Z
    end

    ## Square root free QR/QL (values only, but fast)
    function sterf!(d::StridedVector{Float64}, e::StridedVector{Float64})

        # Extract sizes
        n = length(d)

        # Checks
        length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short"))

        # Allocations
        info = BlasInt[0]

        ccall((@blasfunc("dsterf_"), Base.liblapack_name), Cvoid,
            (Ref{BlasInt}, Ptr{Float64}, Ptr{Float64}, Ptr{BlasInt}),
            n, d, e, info)

        info[1] == 0 || throw(LAPACKException(info[1]))

        d
    end

    ## Divide and Conquer
    function stedc!(compz::Char,
                    d::StridedVector{Float64},
                    e::StridedVector{Float64},
                    Z::StridedMatrix{Float64},
                    work::StridedVector{Float64},
                    lwork::BlasInt,
                    iwork::StridedVector{BlasInt},
                    liwork::BlasInt)

        # Extract sizes
        n = length(d)
        ldz = stride(Z, 2)

        # Checks
        length(e) + 1 == n || throw(DimensionMismatch(""))
        # length(work)
        # length(lwork)

        # Allocations
        info = BlasInt[0]

        ccall((@blasfunc("dstedc_"), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
                 Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
                 Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}),
                compz, n, d, e,
                Z, max(1, ldz), work, lwork,
                iwork, liwork, info)

        info[1] == 0 || throw(LAPACKException(info[1]))

        return d, Z
    end

    function stedc!(compz::Char,
                    d::StridedVector{Float64},
                    e::StridedVector{Float64},
                    Z::StridedMatrix{Float64} = compz == 'N' ?
                        Matrix{Float64}(undef, 0, 0) :
                        Matrix{Float64}(undef, length(d), length(d)))

        work::Vector{Float64} = Float64[0]
        iwork::Vector{BlasInt} = BlasInt[0]
        stedc!(compz, d, e, Z, work, -1, iwork, -1)

        lwork::BlasInt = work[1]
        work = Vector{Float64}(undef, lwork)
        liwork = iwork[1]
        iwork = Vector{BlasInt}(undef, liwork)

        return stedc!(compz, d, e, Z, work, lwork, iwork, liwork)
    end

    ## RRR
    for (lsymb, elty) in ((:dstemr_, :Float64), (:sstemr_, :Float32))
        @eval begin
            function stemr!(jobz::Char,
                            range::Char,
                            dv::StridedVector{$elty},
                            ev::StridedVector{$elty},
                            vl::$elty,
                            vu::$elty,
                            il::BlasInt,
                            iu::BlasInt,
                            w::StridedVector{$elty},
                            Z::StridedMatrix{$elty},
                            nzc::BlasInt,
                            isuppz::StridedVector{BlasInt},
                            work::StridedVector{$elty},
                            lwork::BlasInt,
                            iwork::StridedVector{BlasInt},
                            liwork::BlasInt)

                # Extract sizes
                n = length(dv)
                ldz = stride(Z, 2)

                # Checks
                length(ev) >= n - 1 || throw(DimensionMismatch("subdiagonal is too short"))

                # Allocations
                eev::Vector{$elty} = length(ev) == n - 1 ? [ev; zero($elty)] : copy(ev)
                abstol = Vector{$elty}(undef, 1)
                m      = Vector{BlasInt}(undef, 1)
                tryrac = BlasInt[1]
                info   = Vector{BlasInt}(undef, 1)

                ccall((@blasfunc($lsymb), Base.liblapack_name), Cvoid,
                    (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                    Ptr{$elty}, Ref{$elty}, Ref{$elty}, Ref{BlasInt},
                    Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
                    Ref{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
                    Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                    Ptr{BlasInt}),
                    jobz, range, n, dv,
                    eev, vl, vu, il,
                    iu, m, w, Z,
                    max(1, ldz), nzc, isuppz, tryrac,
                    work, lwork, iwork, liwork,
                    info)

                info[1] == 0 || throw(LAPACKException(info[1]))

                w, Z, tryrac[1]
            end

            function stemr!(jobz::Char,
                            range::Char,
                            dv::StridedVector{$elty},
                            ev::StridedVector{$elty},
                            vl::$elty = typemin($elty),
                            vu::$elty = typemax($elty),
                            il::BlasInt = 1,
                            iu::BlasInt = length(dv))

                n = length(dv)
                w = Vector{$elty}(undef, n)
                if jobz == 'N'
                    # Z = Matrix{$elty}(0, 0)
                    nzc::BlasInt = 0
                    isuppz = Vector{BlasInt}(undef, 0)
                elseif jobz == 'V'
                    nzc = -1
                    isuppz = similar(dv, BlasInt, 2n)
                else
                    throw(ArgumentError("jobz must be either 'N' or 'V'"))
                end

                work = Vector{$elty}(undef, 1)
                lwork::BlasInt = -1
                iwork = BlasInt[0]
                liwork::BlasInt = -1
                Z = Matrix{$elty}(undef, jobz == 'N' ? 1 : n, 1)
                nzc = -1

                stemr!(jobz, range, dv, ev, vl, vu, il, iu, w, Z, nzc, isuppz, work, lwork, iwork, liwork)

                lwork = work[1]
                work = Vector{$elty}(undef, lwork)
                liwork = iwork[1]
                iwork = Vector{BlasInt}(undef, liwork)
                nzc = Z[1]
                Z = similar(dv, $elty, n, nzc)

                return stemr!(jobz, range, dv, ev, vl, vu, il, iu, w, Z, nzc, isuppz, work, lwork, iwork, liwork)
            end
        end
    end

    # Non-symmetric eigenvalue problem
    function lahqr!(wantt::Bool,
                    wantz::Bool,
                    H::StridedMatrix{Float64},
                    ilo::Integer,
                    ihi::Integer,
                    wr::StridedVector{Float64},
                    wi::StridedVector{Float64},
                    iloz::Integer,
                    ihiz::Integer,
                    Z::StridedMatrix{Float64})
        # SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
        #                         ILOZ, IHIZ, Z, LDZ, INFO )

        #      .. Scalar Arguments ..
        #      INTEGER            IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
        #      LOGICAL            WANTT, WANTZ
        #      ..
        #      .. Array Arguments ..
        #      DOUBLE PRECISION   H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )

        n   = LinearAlgebra.checksquare(H)
        ldh = stride(H, 2)
        ldz = stride(Z, 2)

        info = Ref{BlasInt}(0)

        ccall((@blasfunc("dlahqr_"), Base.liblapack_name), Cvoid,
            (Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
             Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
             Ptr{Float64}, Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
             Ref{BlasInt}, Ptr{BlasInt}),
            wantt, wantz, n, ilo,
            ihi, H, max(1, ldh), wr,
            wi, iloz, ihiz, Z,
            max(1, ldz), info)

        if info[] != 0
            throw(LAPACKException(info[]))
        end

        return wr, wi, Z
    end
    function lahqr!(H::StridedMatrix{Float64})
        n = size(H, 1)
        return lahqr!(false,
                      false,
                      H,
                      1,
                      n,
                      Vector{Float64}(undef, n),
                      Vector{Float64}(undef, n),
                      1,
                      1,
                      Matrix{Float64}(undef, 0, 0))
    end

    ## Cholesky + singular values
    function spteqr!(
        compz::Char,
        d::StridedVector{Float64},
        e::StridedVector{Float64},
        Z::StridedMatrix{Float64},
        work::StridedVector{Float64} = Vector{Float64}(undef, 4length(d)))

        n = length(d)
        ldz = stride(Z, 2)

        # Checks
        length(e) >= n - 1 || throw(DimensionMismatch("subdiagonal vector is too short"))
        if compz == 'N'
        elseif compz == in('V', 'I')
            size(Z, 1) >= n     || throw(DimensionMismatch("Z does not have enough rows"))
            size(Z, 2) >= ldz   || throw(DimensionMismatch("Z does not have enough columns"))
        else
            throw(ArgumentError("compz must be either 'N', 'V', or 'I'"))
        end

        info = Vector{BlasInt}(undef, 1)

        ccall((@blasfunc(:dpteqr_), Base.liblapack_name), Cvoid,
            (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
             Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ptr{BlasInt}),
            compz, n, d, e,
            Z, max(1, ldz), work, info)

        info[1] == 0 || throw(LAPACKException(info[1]))

        d, Z
    end

    # Gu's dnc eigensolver
    for (f, elty) in ((:dsyevd_, :Float64),
                      (:ssyevd_, :Float32))

        @eval begin
            function syevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty})
                n      = LinearAlgebra.checksquare(A)
                lda    = stride(A, 2)
                w      = Vector{$elty}(n)
                work   = Vector{$elty}(undef, 1)
                lwork  = BlasInt(-1)
                iwork  = Vector{BlasInt}(undef, 1)
                liwork = BlasInt(-1)
                info = BlasInt[0]
                for i = 1:2
                    ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                        (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                         Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                         Ptr{BlasInt}, Ref{BlasInt}, Ptr{BlasInt}),
                        jobz, uplo, n, A,
                        lda, w, work, lwork,
                        iwork, liwork, info)

                    if info[1] != 0
                        return LinearAlgebra.LAPACKException(info[1])
                    end
                    if i == 1
                        lwork  = BlasInt(work[1])
                        resize!(work, lwork)
                        liwork = iwork[1]
                        resize!(iwork, liwork)
                    end
                end
                return w, A
            end
        end
    end
    for (f, elty, relty) in ((:zheevd_, :ComplexF64, :Float64),
                             (:cheevd_, :ComplexF32, :Float32))

        @eval begin
            function heevd!(jobz::Char, uplo::Char, A::StridedMatrix{$elty})
                n      = LinearAlgebra.checksquare(A)
                lda    = stride(A, 2)
                w      = Vector{$relty}(n)
                work   = Vector{$elty}(undef, 1)
                lwork  = BlasInt(-1)
                rwork  = Vector{$relty}(undef, 1)
                lrwork = BlasInt(-1)
                iwork  = Vector{BlasInt}(undef, 1)
                liwork = BlasInt(-1)
                info = BlasInt[0]
                for i = 1:2
                    ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                        (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                         Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt},
                         Ptr{$relty}, Ref{BlasInt}, Ptr{BlasInt}, Ref{BlasInt},
                         Ptr{BlasInt}),
                        jobz, uplo, n, A,
                        lda, w, work, lwork,
                        rwork, lrwork, iwork, liwork,
                        info)

                    if info[1] != 0
                        return LinearAlgebra.LAPACKException(info[1])
                    end

                    if i == 1
                        lwork = BlasInt(work[1])
                        resize!(work, lwork)
                        lrwork = BlasInt(rwork[1])
                        resize!(rwork, lrwork)
                        liwork = iwork[1]
                        resize!(iwork, liwork)
                    end
                end
                return w, A
            end
        end
    end

# Generalized Eigen
## Eigenvectors of (pseudo) triangular matrices
for (f, elty) in ((:dtgevc_, :Float64),
                   (:stgevc_, :Float32))

    @eval begin
        function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty}, VL::StridedMatrix{$elty}, VR::StridedMatrix{$elty}, work::Vector{$elty})

            n = LinearAlgebra.checksquare(S)
            if LinearAlgebra.checksquare(P) != n
                throw(DimensionMismatch("the two matrices must have same size"))
            end

            lds, ldp, ldvl, ldvr = stride(S, 2), stride(P, 2), stride(VL, 2), stride(VR, 2)
            if side == 'L'
                mm = size(VL, 2)
                if size(VL, 1) != n
                    throw(DimensionMismatch("first dimension of output is $size(VL, 1) but should be $n"))
                end
            elseif side == 'R'
                mm = size(VR, 2)
                if size(VR, 1) != n
                    throw(DimensionMismatch("first dimension of output is $size(VR, 1) but should be $n"))
                end
            elseif side == 'B'
                mm = size(VL, 2)
                if mm != size(VR, 2)
                    throw(DimensionMismatch("the matrices for storing the left and right eigenvectors must have the same number of columns when side == 'B'"))
                end
                if size(VL, 1) != n
                    throw(DimensionMismatch("first dimension of output is $size(VL, 1) but should be $n"))
                end            else
                if size(VR, 1) != n
                    throw(DimensionMismatch("first dimension of output is $size(VR, 1) but should be $n"))
                end
                throw(ArgumentError("side must be either 'L', 'R', or 'B'"))
            end
            if howmny == 'A' || howmny == 'B'
                if mm != n
                    throw(DimensionMismatch("the number of columne in the output matrices are $mm, but should be $n"))
                end
            elseif howmny == 'S'
                mx, mn = extrama(select)
                if mx > 1 || nm < 0
                    throw(ArgumentError("the elements of select must be either 0 or 1"))
                end
                if sum(howmny) != mm
                    throw(DimensionMismatch("the number of columns in the output arrays is $mm, but you have selected $(sum(howmny)) vectors"))
                end
            else
                throw(ArgumentError("howmny must be either A, B, or S"))
            end

            m = BlasInt[0]
            info = BlasInt[0]

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ptr{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ref{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                side, howmny, select, n,
                S, lds, P, ldp,
                VL, ldvl, VR, ldvr,
                mm, m, work, info)

            if info[1] != 0
                throw(LAPACKException(info[1]))
            end

            return VL, VR, m[1]
        end

        function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty}, VL::StridedMatrix{$elty}, VR::StridedMatrix{$elty})

            n = LinearAlgebra.checksquare(S)
            work = Vector{$elty}(6n)

            return tgevc!(side, howmny, select, S, P, VL, VR, work)
        end

        function tgevc!(side::Char, howmny::Char, select::Vector{BlasInt}, S::StridedMatrix{$elty}, P::StridedMatrix{$elty})
            # No checks here as they are done in method above
            n = LinearAlgebra.checksquare(S)
            if side == 'L'
                VR = Matrix{$elty}(n, 0)
                if howmny == 'A' || howmny == 'B'
                    VL = Matrix{$elty}(n, n)
                else
                    VL = Matrix{$elty}(n, sum(select))
                end
            elseif side == 'R'
                VL = Matrix{$elty}(n, 0)
                if howmny == 'A' || howmny == 'B'
                    VR = Matrix{$elty}(n, n)
                else
                    VR = Matrix{$elty}(n, sum(select))
                end
            else
                if howmny == 'A' || howmny == 'B'
                    VL = Matrix{$elty}(n, n)
                    VR = Matrix{$elty}(n, n)
                else
                    VL = Matrix{$elty}(n, sum(select))
                    VR = Matrix{$elty}(n, sum(select))
                end
            end

            return tgevc!(side, howmny, select, S, P, VL, VR)
        end
    end
end

# Rectangular full packed format

## Symmetric rank-k operation for matrix in RFP format.
for (f, elty, relty) in ((:dsfrk_, :Float64, :Float64),
                   (:ssfrk_, :Float32, :Float32),
                   (:zhfrk_, :ComplexF64, :Float64),
                   (:chfrk_, :ComplexF32, :Float32))

    @eval begin
        function sfrk!(transr::Char, uplo::Char, trans::Char, alpha::Real, A::StridedMatrix{$elty}, beta::Real, C::StridedVector{$elty})
            chkuplo(uplo)
            chkstride1(A)
            if trans in ('N', 'n')
                n, k = size(A)
            elseif trans in ('T', 't', 'C', 'c')
                k, n = size(A)
            end
            lda = max(1, stride(A, 2))

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                 Ref{BlasInt}, Ref{$relty}, Ptr{$elty}, Ref{BlasInt},
                 Ref{$relty}, Ptr{$elty}),
                transr, uplo, trans, n,
                k, alpha, A, lda,
                beta, C)
            C
        end
    end
end

# Cholesky factorization of a real symmetric positive definite matrix A
for (f, elty) in ((:dpftrf_, :Float64),
                   (:spftrf_, :Float32),
                   (:zpftrf_, :ComplexF64),
                   (:cpftrf_, :ComplexF32))

    @eval begin
        function pftrf!(transr::Char, uplo::Char, A::StridedVector{$elty})
            chkuplo(uplo)
            n = round(Int,div(sqrt(8length(A)), 2))
            info = Vector{BlasInt}(undef, 1)

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                 Ptr{BlasInt}),
                transr, uplo, n, A,
                info)
            A
        end
    end
end

# Computes the inverse of a (real) symmetric positive definite matrix A using the Cholesky factorization
for (f, elty) in ((:dpftri_, :Float64),
                   (:spftri_, :Float32),
                   (:zpftri_, :ComplexF64),
                   (:cpftri_, :ComplexF32))

    @eval begin
        function pftri!(transr::Char, uplo::Char, A::StridedVector{$elty})
            chkuplo(uplo)
            n = round(Int,div(sqrt(8length(A)), 2))
            info = Vector{BlasInt}(undef, 1)

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                 Ptr{BlasInt}),
                transr, uplo, n, A,
                info)

            A
        end
    end
end

# DPFTRS solves a system of linear equations A*X = B with a symmetric positive definite matrix A using the Cholesky factorization
for (f, elty) in ((:dpftrs_, :Float64),
                   (:spftrs_, :Float32),
                   (:zpftrs_, :ComplexF64),
                   (:cpftrs_, :ComplexF32))

    @eval begin
        function pftrs!(transr::Char, uplo::Char, A::StridedVector{$elty}, B::StridedVecOrMat{$elty})
            chkuplo(uplo)
            chkstride1(B)
            n = round(Int,div(sqrt(8length(A)), 2))
            if n != size(B, 1)
                throw(DimensionMismatch("B has first dimension $(size(B,1)) but needs $n"))
            end
            nhrs = size(B, 2)
            ldb = max(1, stride(B, 2))
            info = Vector{BlasInt}(undef, 1)

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                 transr, uplo, n, nhrs,
                 A, B, ldb, info)

            B
        end
    end
end

# Solves a matrix equation (one operand is a triangular matrix in RFP format)
for (f, elty) in ((:dtfsm_, :Float64),
                   (:stfsm_, :Float32),
                   (:ztfsm_, :ComplexF64),
                   (:ctfsm_, :ComplexF32))

    @eval begin
        function tfsm!(transr::Char,
                       side::Char,
                       uplo::Char,
                       trans::Char,
                       diag::Char,
                       alpha::$elty,
                       A::StridedVector{$elty},
                       B::StridedVecOrMat{$elty})
            chkuplo(uplo)
            chkside(side)
            chkdiag(diag)
            chkstride1(B)
            m, n = size(B, 1), size(B, 2)
            if round(Int, div(sqrt(8length(A)), 2)) != m
                throw(DimensionMismatch("First dimension of B must equal $(round(Int, div(sqrt(8length(A)), 2))), got $m"))
            end
            ldb = max(1, stride(B, 2))

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{UInt8},
                 Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{$elty},
                 Ptr{$elty}, Ptr{$elty}, Ref{BlasInt}),
                transr, side, uplo, trans,
                diag, m, n, alpha,
                A, B, ldb)

            return B
        end
    end
end

# Computes the inverse of a triangular matrix A stored in RFP format.
for (f, elty) in ((:dtftri_, :Float64),
                   (:stftri_, :Float32),
                   (:ztftri_, :ComplexF64),
                   (:ctftri_, :ComplexF32))

    @eval begin
        function tftri!(transr::Char, uplo::Char, diag::Char, A::StridedVector{$elty})
            chkuplo(uplo)
            chkdiag(diag)
            n = round(Int,div(sqrt(8length(A)), 2))
            info = Vector{BlasInt}(undef, 1)

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{UInt8}, Ref{BlasInt},
                 Ptr{$elty}, Ptr{BlasInt}),
                transr, uplo, diag, n,
                A, info)

            A
        end
    end
end

# Copies a triangular matrix from the rectangular full packed format (TF) to the standard full format (TR)
for (f, elty) in ((:dtfttr_, :Float64),
                   (:stfttr_, :Float32),
                   (:ztfttr_, :ComplexF64),
                   (:ctfttr_, :ComplexF32))

    @eval begin
        function tfttr!(transr::Char, uplo::Char, Arf::StridedVector{$elty})
            chkuplo(uplo)
            n = round(Int,div(sqrt(8length(Arf)), 2))
            info = Vector{BlasInt}(undef, 1)
            A = similar(Arf, $elty, n, n)

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}),
                transr, uplo, n, Arf,
                A, n, info)

            A
        end
    end
end

# Copies a triangular matrix from the standard full format (TR) to the rectangular full packed format (TF).
for (f, elty) in ((:dtrttf_, :Float64),
                   (:strttf_, :Float32),
                   (:ztrttf_, :ComplexF64),
                   (:ctrttf_, :ComplexF32))

    @eval begin
        function trttf!(transr::Char, uplo::Char, A::StridedMatrix{$elty})
            chkuplo(uplo)
            chkstride1(A)
            n = size(A, 1)
            lda = max(1, stride(A, 2))
            info = Vector{BlasInt}(undef, 1)
            Arf = similar(A, $elty, div(n*(n+1), 2))

            ccall((@blasfunc($f), Base.liblapack_name), Cvoid,
                (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ptr{$elty},
                 Ref{BlasInt}, Ptr{$elty}, Ptr{BlasInt}),
                transr, uplo, n, A,
                lda, Arf, info)

            Arf
        end
    end
end

end


text:
using ArgParse, CSV, TimerOutputs, FileIO, Colors, DataFrames, Images
using LAP_julia

println("running file...")

function save_string_to_file(str, file_path)
    open(file_path, "w") do io
        write(io, str)
    end
end

function load_image_gray(img_path)
    img = Float64.(Gray.(load(img_path)))
end

function load_image(img_path)
    img = load(img_path)
end

function save_shift_landmarks(source_landmarks_path, flow, args_dict)
    source_landmarks_df = CSV.read(source_landmarks_path) |> DataFrame
    locations = [source_landmarks_df[:, "Y"] source_landmarks_df[:, "X"]]
    shifted_landmarks = LAP_julia.move_landmarks(locations, flow)
    shifted_landmarks_df = DataFrame(Column1 = 0:(size(locations,1)-1), X = shifted_landmarks[:,2], Y = shifted_landmarks[:,1])
    save_path = joinpath(args_dict["output_path"], args_dict["land_warped_fname"])
    CSV.write(save_path, shifted_landmarks_df)
end

function save_time(time_in_secs, args_dict)
    save_string_to_file(string(time_in_secs), joinpath(args_dict["output_path"], args_dict["time_fname"]))
end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--opt1"
            help = "an option with an argument"
        "--time_fname"
            help = "name of the file to save execution time"
            arg_type = String
            default = "exec_time_in_secs.txt"
        "--land_warped_fname"
            help = "name of the file to save warped landmarks"
            arg_type = String
            default = "warped_landmarks.csv"
        "--diag_pix"
            help = "an option without argument, i.e. a flag"
            arg_type = Int
            default = 500
        "img_path"
            help = "target image path"
            required = true
        "imgw_path"
            help = "source image path"
            required = true
        "moving_landmark_path"
            help = "source landmark path"
            required = true
        "output_path"
            help = "output path"
            required = true
        "reg_alg"
            help = "registration algorithm"
            required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Running with the arguments:")
    args_dict = Dict()
    println("******ARGS:********")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
        args_dict[arg] = val
    end
    args_dict["reg_alg"] = getfield(LAP_julia, Symbol(args_dict["reg_alg"]))

    target = load_image_gray(args_dict["img_path"])
    source = load_image_gray(args_dict["imgw_path"])

    resized_source, source_resize_ratio = resize_to_diag_size(source, args_dict["diag_pix"])
    resized_target = imresize(target, size(resized_source))

    time_in_secs = 0
    # Take the 3 run of the reg alg, allowing for julia compilation.
    for _ in 1:3
        timer = TimerOutput(string(args_dict["reg_alg"]))
        flow_est, source_reg, timer, time_in_secs = time_reg_alg(args_dict["reg_alg"],
                                                                 ready_target,
                                                                 ready_source,
                                                                 timer=timer,
                                                                 method_kwargs=Dict(:display => false, :timer => timer, :match_source_histogram => true))
        global flow_est, source_reg, timer, time_in_secs
        println()

    end

    # save time
    save_time(time_in_secs, args_dict)

    # rescale flow_est
    flow_est = imresize(flow_est.*(1/source_resize_ratio), ratio=1/source_resize_ratio)
    @assert size(flow_est) == size(source)

    # move landmarks
    save_shift_landmarks(args_dict["moving_landmark_path"], flow_est, args_dict)

    # save source_reg
    colored_source = load_image(args_dict["imgw_path"])
    source_reg = warp_img(colored_source, real.(flow_est), imag.(flow_est))
    save_path = joinpath(args_dict["output_path"], "source_reg.jpg")
    save(save_path, source_reg)
end

main()


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# Wrapper type for solutions from Trixi's own time integrators, partially mimicking
# DiffEqBase.ODESolution
struct TimeIntegratorSolution{tType, uType, P}
  t::tType
  u::uType
  prob::P
end

include("methods_2N.jl")
include("methods_3Sstar.jl")

end # @muladd


text:
@benchgroup "max-flow" begin
  for n in 9:5:29
    srand(1)
    p = 8.0 / n
    A = sprand(n,n,p)
    g = DiGraph(A)
    cap = round.(A*100)
    @bench "n = $n" LightGraphs.maximum_flow($g, 1, $n, $cap)
  end
end # max-flow

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
export  w_legendre,
        build_w_jacobi,
        w_jacobi,
        w_laguerre,
        w_hermite,
        build_w_genhermite,
        build_w_genlaguerre,
        w_meixner_pollaczek,
        build_w_meixner_pollaczek,

        w_gaussian,
        w_uniform01,
        w_uniform_11,
        w_logistic,
        w_genhermite,
        build_w_beta,
        build_w_gamma

_throwError(t) = throw(DomainError(t, "not in support"))

function w_legendre(t)
    -1. <= t <= 1. ? 1. : _throwError(t)
end

function build_w_jacobi(a,b)
    return t->w_jacobi(t,a,b)
end

function w_jacobi(t,a,b)
    -1. <= t <= 1. ? (1-t)^a*(1+t)^b : _throwError(t)
end

function w_hermite(t)
    exp(-t^2)
end

function build_w_genhermite(mu)
    return t->w_genhermite(t,mu)
end

function w_genhermite(t,μ)
    abs(t)^(2*μ)*exp(-t^2)
end

function build_w_genlaguerre(a)
    return t -> w_genlaguerre(t,a)
end

function w_laguerre(t)
    t >= 0. ? exp(-t) : _throwError(t)
end

function w_meixner_pollaczek(t,lambda,phi)
    1 / (2pi) * exp((2*phi - pi)*t) * abs(gamma(lambda+im*t))^2
end

function build_w_meixner_pollaczek(lambda,phi)
    t->w_meixner_pollaczek(t,lambda,phi)
end

##################################################
# probability density functions
function w_gaussian(t)
    1 / (sqrt(2*pi))*exp(-0.5*t^2)
end

function build_w_beta(α,β)
    return t->w_beta(t,α,β)
end

function w_beta(t,α,β)
    -1 <= t <= 1 ? t^(α-1)*(1-t)^(β-1)/beta(α,β) : _throwError(t)
end

function build_w_gamma(α)
    return t->w_gamma(t,α)
end

function w_gamma(t,α)
    t >= 0. ? (1/gamma(α)*Float64(t)^(α-1)*exp(-t)) : _throwError(t)
end

function w_uniform01(t)
    0. <= t <= 1. ? 1. : _throwError(t)
end

function w_uniform_11(t)
    -1. <= t <= 1. ? 0.5 : _throwError(t)
end

function w_logistic(t)
    0.25*sech(0.5t)^2
end


text:
#
# SIRD model trajectory structure
#

struct SIRDTraj
  nsteps :: Int64
  U :: Vector{Float64}
  S :: Vector{Float64}
  D :: Vector{Float64}
  I :: Vector{Float64}
  time :: Vector{Float64}
end

function SIRDTraj( input :: SIRDInput )
  nsteps = input.nsave
  U = Vector{Float64}(undef,nsteps)
  S = Vector{Float64}(undef,nsteps)
  D = Vector{Float64}(undef,nsteps)
  I = Vector{Float64}(undef,nsteps)
  time = Vector{Float64}(undef,nsteps)
  return SIRDTraj(nsteps,U,S,D,I,time)
end


text:
using Unitful, UnitfulRecipes, UnitfulAstro
using Dierckx

#
# The file containing the various differential equation systems to solve.
#

#--------------------------
# Definition of equations
#--------------------------

"The differential equation system derived by Draine 2011"
function d2011system!(dv, v, p, y)
    u, ϕ, τ = v
    γ, β = p

    dv[1] = du = -1 - γ*(β*exp(-τ) + ϕ)*u/y^2
    dv[2] = dϕ = -(y/u)^2 - γ * ϕ/u
    dv[3] = dτ = γ/u
end

"""
    draine2011model(γ, β; u0=454, ystart=10, yend=100)

Solve the Draine (2011) model with model parameters γ and β.
This formulates the ODEProblem and solves it

"""
function draine2011model(γ, β; u0=454.0, ystart=10.0, yend=100.0)

    starting_point = [u0, 1.0, 0.0]
    yspan = [ystart, yend]

    p = [γ, β]

    # The system is only defined for ϕ>=0 so we interrupt on ϕ=u[2]=0
    condition(u, y, integrator) = u[2]
    affect!(integrator) = terminate!(integrator)
    cb = ContinuousCallback(condition, affect!);
    prob = ODEProblem(d2011system!, starting_point, yspan, p)

    sol = solve(prob, Tsit5(), callback=cb);

    result = Draine2011(sol, Draine2011Parameters(γ, β), u0,
                        starting_point[2], starting_point[3],
                        ystart, ymax(sol))

    return result

end

draine2011model(d::NamedTuple; kwargs...) = draine2011model(d.γ, d.β; u0=d.u0, kwargs...)

function d2011wgravity!(dv, v, h, p, y)
    u, ϕ, τ = v
    γ, β, η, ymin = p # Ymin is indirectly a parameter of the model

    rint(x) = x^2/h(p, x)[1]
    Ry, err = quadgk(rint, ymin, y)

    dv[1] = du = u^2*η*Ry/y^2 -1 - γ*(β*exp(-τ) + ϕ)*u/y^2
    dv[2] = dϕ = -(y/u)^2 - γ * ϕ/u
    dv[3] = dτ = γ/u
end

function d2011wgravitymodel(γ, β, η; u0=454, ystart=10, yend=100)

    starting_point = [u0, 1, 0]
    yspan = [ystart, yend]

    p = [γ, β, η, ystart]

    # The system is only defined for ϕ>=0 so we interrupt on ϕ=u[2]=0
    condition(u, y, integrator) = u[2]
    affect!(integrator) = terminate!(integrator)
    cb = ContinuousCallback(condition, affect!);

    # This specifies the original starting point. I think it is probably
    # necessary to rethink this  but haven't done this yet
    history(p, t) = [u0, 1, 0]
    prob = DDEProblem(d2011wgravity!, starting_point, history, yspan, p)

    # The choice for stepper is following https://diffeq.sciml.ai/stable/solvers/dde_solve/
    # see also https://github.com/SciML/DifferentialEquations.jl/issues/383
    sol = solve(prob, MethodOfSteps(RK4()), callback=cb);

    return sol

end

#
# Handling of solutions themselves
#

ymax(solution) = maximum(solution.t)
ymin(solution) = minimum(solution.t)

# Return the functions of the differential equation
tau(y, model::AbstractDustyHII) = model.solution(y; idxs=3)
phi(y, model::AbstractDustyHII) = model.solution(y; idxs=2)
ufunc(y, model::AbstractDustyHII) = model.solution(y; idxs=1)

#
# To-from initial conditions.
#
# In this case we have two options:
#
#    1. The user might give physical parameters of a model and want
#       the parameters to use for initial conditions and model parameters γ and β.
#
#    2. The user might have a solution for a given γ and β and wants physical
#       parameters that are consistent with this. This is the more complex case
#       because of the various extra parameters that might have to be passed.
#

"""
    physical_to_model(; Teff=45000u"K", Tgas=1e4u"K", ninner=1e3u"1/cm^3")

Calculates β, Q0, hν_i for a star with the given Teff and calculates γ and
η for the other properties.

"""
function physical_to_model(;Teff=45000u"K", Tgas=1e4u"K", ninner=1e3u"1/cm^3",
                           σ=1e-21*u"1/cm^2")

    # Calculate first a Black-body spectrum - I fix wavelength ranges to
    # be reasonable for a hot star.
    lrange = collect(LinRange(1., 6e4, 10000))
    flux = [PlanckAstro(l*u"angstrom", Teff) for l in lrange];
    itp = interpolate_spectrum(lrange, flux)
    resStar = specQuantities(lrange, itp)

    # Then calculate γ
    γ = gamma(;T=Tgas, Q0=resStar.Q0, hnui=resStar.hnui, σ=σ)

    #  η
    η = upreferred(Eg(Tgas; Q0=resStar.Q0, hnui=resStar.hnui, σ=σ)/Et(Tgas))

    # u0
    u0 = ustrip(uval(ninner; T=Tgas, Q0=resStar.Q0, hnui=resStar.hnui))

    return (Ln=resStar.Ln, Li=resStar.Li, hnui=resStar.hnui,
            Q0=resStar.Q0, β=resStar.β, γ=γ, η=η, u0=u0, Teff=Teff,
            Tgas=Tgas, ninner=ninner, σ=σ)

end

"""
    model_to_physical(m)

Takes a model which has been calculated in unitless form and converts quantities
to physical units.

NOT YET FUNCTIONAL
"""
function model_to_physical(m, T, Q0, hnui, σ)
    return Draine2011Physical(m, T, Q0, hnui, σ)

end


text:
using ThreadPinning
using BandwidthBenchmark
using UnicodePlots

# Compact Pinning
pinthreads(:compact)
println("Check compact pinning: ", getcpuids(), "\n")

# bwbench(; verbose=true);
println("bwbench():")
bwbench(; verbose=true, write_allocate=true);

println("bwscaling():")
data = bwscaling()
p = lineplot(data[:,1], data[:,2], title = "Bandwidth Scaling, Compact Pinning", xlabel = "# cores", ylabel = "MB/s", border = :ascii, canvas = AsciiCanvas)
println(p)

println("flopsscaling():")
data = BandwidthBenchmark.flopscaling()
p = lineplot(data[:,1], data[:,2], title = "Flops Scaling, Compact Pinning", xlabel = "# cores", ylabel = "MFlops/s", border = :ascii, canvas = AsciiCanvas)
println(p)

# Scatter Pinning
pinthreads(:scatter)
println("Check scatter pinning: ", getcpuids(), "\n")

# bwbench(; verbose=true);
println("bwbench():")
bwbench(; verbose=true, write_allocate=true);

println("bwscaling():")
data = bwscaling()
p = lineplot(data[:,1], data[:,2], title = "Bandwidth Scaling, Scattered Pinning", xlabel = "# cores", ylabel = "MB/s", border = :ascii, canvas = AsciiCanvas)
println(p)

println("flopsscaling():")
data = BandwidthBenchmark.flopscaling()
p = lineplot(data[:,1], data[:,2], title = "Flops Scaling, Scattered Pinning", xlabel = "# cores", ylabel = "MFlops/s", border = :ascii, canvas = AsciiCanvas)
println(p)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
####################################################################
function gpmodelHMC(x, flux, σ2, idx; maxiterml = 10000, n_samples=2_000, n_adapts =1_000, maxrandom=500, seed = 1, show_trace = true)
####################################################################

    #--------------------------------------------------------
    # Auxiliary
    #--------------------------------------------------------

    rg = MersenneTwister(seed)

    JITTER = 1e-8

    Σdiag = Diagonal(σ2) + JITTER*I

    prvlogs = NaN

    numF = length(unique(idx))

    # @assert(numF == length(galaxyline(3000.0)))

    numparam = 2*numF + 1

    @assert(all(unique(idx) .== collect(1:numF)))

    N = length(x)

    @assert(N == length(flux) == length(idx) == length(σ2))

    @printf("Running gpmodelHMC with %d and %d filters\n", N, numF)

    @printf("\tNumber of free parameters is %d\n", numparam)

    #--------------------------------------------------------
    # Pre-allocate arrays and pre-calculate distances
    #--------------------------------------------------------

    A    = Diagonal(zeros(N))

    bvec = zeros(N)

    K    = zeros(N, N)

    C    = zeros(N, N)

    prvlogs = NaN

    dist2 = pairwise(SqEuclidean(), reshape(x, 1, N), dims=2)

    #--------------------------------------------------------
    function unpack(param)
    #--------------------------------------------------------

        @assert(length(param) == numparam)

        local a     = param[1:numF]
        local b     = param[numF+1:2*numF]
        local logs  = param[end]

        a, b, logs

    end

    #--------------------------------------------------------
    function objective(param)
    #--------------------------------------------------------

        local a, b, logs = unpack(param)

        -1.0 * logl(a = a, logs = logs, b = b)

    end

    #--------------------------------------------------------
    function blockvector!(bvec, idx, b)
    #--------------------------------------------------------

        for i in 1:N
            bvec[i] = b[idx[i]]
        end

    end

    #--------------------------------------------------------
    function diagonalmatrix!(A, idx, a)
    #--------------------------------------------------------

        for i in 1:N
            A.diag[i] = a[idx[i]]
        end

    end

    #--------------------------------------------------------
    function logl(; a = a, b = b, logs = logs)
    #--------------------------------------------------------

        local K = rbf(dist2, [0.0; logs])
        local A = Diagonal(zeros(eltype(a), N))
        local bvec = zeros(eltype(a), N)

        for i in 1:N
                bvec[i] = b[idx[i]]

                A.diag[i] = a[idx[i]]
        end

        C = (A * K * A) .+ Σdiag ; C .= (C .+ C').*0.5

        local L = try

            cholesky(C).L

        catch err

            if isa(err, PosDefException)

                @warn("Cholesky failed, matrix not posdef. Return -Inf.")

                return -Inf

            else

                throw(err)

            end

        end

        local logdetC = 2.0*sum(log.(diag(L)))

        bvec .= L\(flux .- bvec)

        local yTCinvy = sum(bvec.^2)

        return - 0.5*logdetC - 0.5*yTCinvy - 0.5*N*log(2.0π)

    end

    # a = randn(rg, 6).+1; b = randn(rg, 6).+1
    # @show logl(a=a,b=b,logs=0.0) logl0(a=a,b=b,logs=0.0)
    # # @btime  $logl(a=$a,b=$b,logs=0.0)
    # # @btime $logl0(a=$a,b=$b,logs=0.0)
    # return 0,0

    #--------------------------------------------------------
    # Initialise GP hyperparameters with random search
    # and few iterations
    #--------------------------------------------------------

    initsolutions       = [randn(rg, numparam) for i=1:maxrandom]
    initialfitness      = @showprogress "Initial random search " map(objective, initsolutions)
    bestinitialsolution = initsolutions[argmin(initialfitness)]

    #--------------------------------------------------------
    # Call optimiser
    #--------------------------------------------------------

    opt    = Optim.Options(iterations = maxiterml, show_trace = show_trace, show_every=100)
    result = Optim.optimize(objective, bestinitialsolution, NelderMead(), opt)

    #--------------------------------------------------------
    # Retrieve optimised parameters
    #--------------------------------------------------------

    a, b, logs = unpack(result.minimizer)

    #--------------------------------------------------------
    # Inference via HMC
    #--------------------------------------------------------

    newobjective(x) = -1.0* objective([x; logs])

    # Choose parameter dimensionality and initial parameter value
    initial_θ = [a;b]; D = length(initial_θ)

    # Define the target distribution
    ℓπ(θ) = newobjective(θ)

    # Define a Hamiltonian system
    metric = DiagEuclideanMetric(D)
    hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)

    # Define a leapfrog solver, with initial step size chosen heuristically
    initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)
    integrator = Leapfrog(initial_ϵ)

    # Define an HMC sampler, with the following components
    #   - multinomial sampling scheme,
    #   - generalised No-U-Turn criteria, and
    #   - windowed adaption for step-size and diagonal mass matrix
    proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    # Run the sampler to draw samples from the specified Gaussian, where
    #   - `samples` will store the samples
    #   - `stats` will store diagnostic statistics for each sample
    return samples, stats = sample(hamiltonian, proposal, initial_θ, n_samples, adaptor, n_adapts; progress=true)

    #--------------------------------------------------------
    # Reconstruct curves on latent points
    #--------------------------------------------------------

    calculatekernelmatrix!(K, dist2, rbf, [0.0; logs])

    K = (K+K')*0.5 + JITTER*I # numerical stability

    blockvector!(bvec, idx, b)

    diagonalmatrix!(A, idx, a)

    L    = Σdiag\I

    Λ    = K\I

    Σ = (Λ + A'*L*A)\I

    μ = Σ * (A' * L * (flux .- bvec) )

    figure(101)

    cla()

    clr = ["b","g","r","m","k","c"]

    for (i,l) in enumerate(unique(idx))

        aux = findall(idx .== l)

        plot(t[aux], flux[aux], clr[i]*".", alpha=0.5)

        plot(t, b[i] .+ a[i]*μ, clr[i]*".")

    end

    # #--------------------------------------------------------
    # # Reconstruct curves on test points
    # #--------------------------------------------------------
    #
    #
    #
    #     xtest = collect(LinRange(minimum(x), maximum(x), 200))
    #
    #     # dimensions: N × Ntest
    #     k = calculatekernelmatrix(x, xtest, rbf, [0.0;logs])
    #
    #     # Ntest × 1
    #     c = calculatekernelmatrix(xtest, xtest, rbf, [0.0;logs])
    #
    #     figure(102)
    #
    #     cla()
    #
    #     # predicted mean
    #
    #     μ = k' * ((A'*K*A + Σdiag) \ (flux .- bvec))
    #
    #     # predicted covariance
    #
    #     #C
    #
    #     for (i,l) in enumerate(unique(idx))
    #
    #         plot(xtest, a[i]*μ .+ b[i], clr[i]*".-")
    #
    #         aux = findall(idx .== l)
    #
    #         plot(x[aux], flux[aux], clr[i]*".", alpha=0.5)
    #
    #     end

    return a, b

end


text:
S = 35
N = 41

Theta = ones(S)
tauhat = rand(N*S, N*S)
tauFhat = rand(N*S, N)

wagehat = rand(N)
iter = 1

wagehat0 = wagehat
PrsjHat0 = rand(S, N*S)
PrFjHat0 = rand(S, N)

tol = 1e-3
maxit = 100
pfmax = 1
iterph = 1

cost = ones(S, N)

for s in 1:S
    for j in 1:N
            cost[s,j] = wagehat[j]^VAcoeff[s,j]
        for r in 1:S
            cost[s,j] = cost[s,j]*PrsjHat0[r,(j-1)*S+s]^gammas[r,(j-1)*S+s]
        end
    end
end

PrsijHat = zeros(N*S, N*S)

for i in 1:N
    for r in 1:S
        for j in 1:N
            for s in 1:S
                iirr = (i-1)*S+r
                jjss = (j-1)*S+s
                PrsijHat[iirr,jjss] = (tauhat[iirr,jjss]*cost[r,i])^(-Theta[r])
            end
        end
    end
end

PrsijHat

cost_Z = Pi .* PrsijHat
PrsjHat = [sum(cost_Z[i:S:(N-1)*S+i, j])^(-1/θ[i,ceil(Int,j/S)]) for i in 1:S, j in 1:N*S]

PrFijHat = zeros(N*S,N)
for i in 1:N
    for r in 1:S
        for j in 1:N
            iirr = (i-1)*S+r
            PrFijHat[iirr,j] = (tauFhat[iirr,j]*cost[r,i])^(-Theta[r])
        end
    end
end

PrFijHat

############### mine
cost_hat_w = [wagehat[ceil(Int,i/S)]^VA_coeff[i] for i in 1:N*S] # NS×1, wages
cost_hat_p = PrsjHat0.^γ
cost_hat_p = [prod(cost_hat_p[:,i]) for i in 1:N*S] # NS×1, price indices (prod is the same as sum just for multiplication)
cost_hat = cost_hat_w .* cost_hat_p # NS×1
cost_hat2 = reshape(cost_hat, S, N)

θ = ones(N*S)
cost_hat_Z = zeros(N*S, N*S)
θ = reshape(θ,S,N)

for i in 1:N
    for r in 1:S
        for j in 1:N
            for s in 1:S
                iirr = (i-1)*S+r
                jjss = (j-1)*S+s
                cost_hat_Z[iirr,jjss] = (cost_hat2[r,i]*tauhat[iirr,jjss])^(-θ[r,i])
            end
        end
    end
end

count(abs.(cost_hat_Z .- PrsijHat) .> 1e-5)

cost_Z = π_Z .* cost_hat_Z
P_hat_Z = [sum(cost_Z[i:S:(N-1)*S+i, j])^(-1/θ[i,ceil(Int,j/S)]) for i in 1:S, j in 1:N*S]

count(abs.(PrsjHat .- P_hat_Z) .> 1e-15)
count(isinf.(P_hat_Z))

cost_hat_F = zeros(N*S,N)
for i in 1:N
    for r in 1:S
        for j in 1:N
            iirr = (i-1)*S+r
            cost_hat_F[iirr,j] = (cost_hat2[r,i]*tauFhat[iirr,j])^(-θ[r,i])
        end
    end
end
cost_hat_F

count(abs.(PrFijHat .- cost_hat_F) .> 1e-10)

cost_F = π_F .* cost_hat_F

P_hat_F = [sum(cost_F[i:S:(N-1)*S+i, j])^(-1/θ[i,j]) for i in 1:S, j in 1:N] # S×N

θ = vec(θ)
θ = repeat(θ, 1, N*S)
π_hat_Z = (cost_hat_Z ./ repeat(P_hat_Z, N)) .^ (.-θ) # NS×NS

# -----------------------------------------------------------------------------------------------------------------------------
# wages

# Antras
# -----------------

CtyConsp = wagehat .* VAn .- Surplusn
FGS = π_prime_F .* repeat(alphas, N) .* repeat(CtyConsp', N*S)
FGS = [sum(FGS[i,:]) for i in 1:N*S]
IGS = π_prime_Z .* repeat(gammas, N)
Youtput = inv(I - IGS)*FGS

LHS = zeros(N)
RHS = zeros(N)
for j in 1:N
    for s in 1:N
        for i in 1:N
            iiss = (i-1)*S+s
            jjss = (j-1)*S+s
            for r in 1:S
                jjrr = (j-1)*S+r
                iirr = (i-1)*S+r
                LHS[j] += π_prime_Z[iiss,jjrr]*gammas[s,jjrr]*Youtput[jjrr]
                RHS[j] += π_prime_Z[jjss,iirr]*gammas[s,iirr]*Youtput[iirr]
            end
            LHS[j] += π_prime_F[iiss,j]*alphas[s,j]*CtyConsp[j]
            RHS[j] += π_prime_F[jjss,i]*alphas[s,i]*CtyConsp[i]
        end
    end
end

ExcessSurplus = RHS .- LHS .- Surplusn # N×1, excess trade balance

vfactor = 0.4
NES = ExcessSurplus ./ VAn
delta = sign.(NES) .* abs.(vfactor .* NES)
wagehat_new = wagehat .* (1 .+ delta ./ wagehat)

# MINE ---------------------

F_prime_ctry = wagehat .* VA_ctry .- TB_ctry # N×1, counterfactual country final goods consumption
# Goods market clearing from equation (35)
total_sales_F = π_prime_F .* repeat(α,N) .* repeat(transpose(F_prime_ctry),N*S) # NS×N
total_sales_F = [sum(total_sales_F[i,:]) for i in 1:N*S] # N×1

A_prime = π_prime_Z .* repeat(γ, N) # NS×NS, intermediate input coefficient matrix

Y_prime = inv(I - A_prime) * total_sales_F #  NS×1

count(abs.(Surplusn .- TB_ctry) .> 0.1)

count(abs.(VAn .- VA_ctry) .> 0.1)

count(abs.(Youtput .- Y_prime) .> 0.1)

Z_prime = π_prime_Z .* repeat(γ, N) .* repeat(transpose(Y_prime), N*S) # NS×NS
F_prime = π_prime_F .* repeat(α, N) .* repeat(transpose(F_prime_ctry), N*S) # NS×N

E_prime = [sum(Y_prime[i:i+S-1]) - sum(Z_prime[i:i+S-1,i:i+S-1]) - sum(F_prime[i:i+S-1,ceil(Int, i/S)]) for i in 1:S:N*S] # N×1

M_prime_Z = [sum(Z_prime[:,j:j+S-1]) - sum(Z_prime[j:j+S-1,j:j+S-1]) for j in 1:S:N*S]
M_prime_F = [sum(F_prime[:,ceil(Int, j/S)]) - sum(F_prime[j:j+S-1,ceil(Int, j/S)]) for j in 1:S:N*S]
M_prime = M_prime_Z .+ M_prime_F # N×1

ETB_ctry = E_prime .- M_prime .- TB_ctry # N×1

a = ExcessSurplus .- ETB_ctry

norm_ETB_ctry =  ExcessSurplus ./ VA_ctry # N×1, normalized excess trade balance
δ = sign.(norm_ETB_ctry) .* abs.(vfactor.*norm_ETB_ctry) # i.e. if norm_ETB_ctry > 0 => too much exports in model => wages should increase (sign fⁿ gives ± of surplus)
w_hat = wagehat .* (1.0 .+ δ ./ wagehat) # N×1, increase/decrease wages for countries with an excess surplus/deficit

F_prime_ctry = w_hat .* VA_ctry .- TB_ctry # N×1, counterfactual country final goods consumption
# Goods market clearing from equation (35)
total_sales_F = π_prime_F .* repeat(α,N) .* repeat(transpose(F_prime_ctry),N*S) # NS×N
total_sales_F = [sum(total_sales_F[i,:]) for i in 1:N*S] # NS×1
A_prime = π_prime_Z .* repeat(γ, N) # NS×NS, intermediate input coefficient matrix
Y_prime = inv(I - A_prime) * total_sales_F #  NS×1

ts = π_F .* repeat(α,N) .* repeat(transpose(F_ctry),N*S)
ts = [sum(ts[i,:]) for i in 1:N*S] # NS×1
A = π_Z .* repeat(γ, N)
Y2 = inv(I - A) * ts

LHS = zeros(N)
RHS = zeros(N)
for j in 1:N
    for s in 1:N
        for i in 1:N
            iiss = (i-1)*S+s
            jjss = (j-1)*S+s
            for r in 1:S
                jjrr = (j-1)*S+r
                iirr = (i-1)*S+r
                LHS[j] += π_prime_Z[iiss,jjrr]*γ[s,jjrr]*Y_prime[jjrr]
                RHS[j] += π_prime_Z[jjss,iirr]*γ[s,iirr]*Y_prime[iirr]
            end
            LHS[j] += π_prime_F[iiss,j]*α[s,j]*F_prime_ctry[j]
            RHS[j] += π_prime_F[jjss,i]*α[s,i]*F_prime_ctry[i]
        end
    end
end

ETB_ctry = RHS .- LHS .- TB_ctry # N×1, excess trade balance

# adjust wages to excess trade balance
norm_ETB_ctry =  ETB_ctry ./ VA_ctry # N×1, normalized excess trade balance
δ = sign.(norm_ETB_ctry) .* abs.(vfactor.*norm_ETB_ctry) # i.e. if norm_ETB_ctry > 0 => too much exports in model => wages should increase (sign fⁿ gives ± of surplus)
w_hat3 = w_hat .* (1.0 .+ δ ./ w_hat) # N×1, increase/decrease wages for countries with an excess surplus/deficit

w_hat = rand(N)
P0_hat_Z = rand(S, N*S)

cost_hat_w = [w_hat[ceil(Int,i/S)]^VA_coeff[i] for i in 1:N*S] # NS×1, wages
cost_hat_p = P0_hat_Z.^γ
cost_hat_p = [prod(cost_hat_p[:,i]) for i in 1:N*S] # NS×1, price indices (prod is the same as sum just for multiplication)
cost_hat = cost_hat_w .* cost_hat_p # NS×1
cost_hat = reshape(cost_hat, S, N)

cost_Z = π_Z .* cost_hat_Z # NS×NS, origin country-industry destination country-industry price index (inside of summation)

# sum over origin countries => price index of country-industry composite good in industry
P_hat_Z = [sum(cost_Z[i:S:(N-1)*S+i, j])^(-1/θ[i,ceil(Int,j/S)]) for i in 1:S, j in 1:N*S]

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

cost = ones(S,N)
VAcoeff = reshape(VA_coeff, S, N)
for s in 1:S
    for j in 1:N
        cost[s,j] = w_hat_prev[j]^VAcoeff[s,j]
        for r in 1:S
            cost[s,j] = cost[s,j]*P_hat_Z[r,(j-1)*S+s]^γ[r,(j-1)*S+s]
        end
    end
end

PrsijHat = zeros(N*S, N*S)
θ = reshape(θ, S, N)

for i in 1:N
    for r in 1:S
        for j in 1:N
            for s in 1:S
                iirr = (i-1)*S+r
                jjss = (j-1)*S+s
                PrsijHat[iirr,jjss] = (τ_hat_Z[iirr,jjss]*cost[r,i])^(-θ[r,i])
            end
        end
    end
end

PrsijHat

PrsjHat = zeros(S,N*S)

for r in 1:S
    for j in 1:N
        for s in 1:S
            jjss = (j-1)*S+s
            PrsjHat[r,jjss] = (sum(π_Z[r:S:(N-1)*S+r,jjss] .* PrsijHat[r:S:(N-1)*S+r,jjss]))^(-1/θ[r,j])
        end
    end
end
PrsjHat

PiHat = PrsijHat ./ (repmat(PrsjHat,N,1))

w_hat = ones(N)

F_prime_ctry = w_hat .* VA_ctry .- TB_ctry # N×1, counterfactual country final goods consumption

# Goods market clearing from equation (35)
total_sales_F = π_F .* repeat(α,N) .* repeat(transpose(F_prime_ctry),N*S) # NS×N
total_sales_F = [sum(total_sales_F[i,:]) for i in 1:N*S] # N×1

A_prime = π_Z .* repeat(γ, N) # NS×NS, intermediate input coefficient matrix

Y_prime = inv(I - A_prime) * total_sales_F #  NS×1
Y_prime = ifelse.(Y_prime .< 0.0, 0.0, Y_prime) # Antras and Chor (2018)

F_prime_ctry ≈ F_ctry

count(abs.(F_prime_ctry.-F_ctry).>0.7)

Z_prime = π_Z .* repeat(γ, N) .* repeat(transpose(Y), N*S) # NS×NS
F_prime = π_F .* repeat(α, N) .* repeat(transpose(F_ctry), N*S) # NS×N
count(abs.(Z_prime.-Z).>1)
count(abs.(F_prime.-F).>1)

max_stars_repo_path: old stuff/prices.jl

max_issues_repo_path: old stuff/prices.jl

max_forks_repo_path: old stuff/prices.jl


text:
"""
Cache for graphics pipelines.

!!! warning
    This cache is not thread-safe.
"""
struct GraphicsPipelineCache
    device::Device
    pipelines::Dictionary{GraphicsPipelineCreateInfo,Pipeline}
end

GraphicsPipelineCache(device) = GraphicsPipelineCache(device, Dictionary{GraphicsPipelineCreateInfo,Pipeline}())

"""
Retrieve or create graphics pipelines, depending on whether they were already cached.
Note that because of caching, you shouldn't expect duplicate pipelines if you provide `GraphicsPipelineCreateInfo`.
If duplicates are present, the behavior differs from `Vulkan.create_graphics_pipelines` in that
they will all be associated with the same pipeline.
"""
function Base.get!(cache::GraphicsPipelineCache, create_infos::AbstractVector{GraphicsPipelineCreateInfo}; allocator = C_NULL, pipeline_cache = C_NULL)
    info_uncached = filter(Base.Fix1(!haskey, cache.pipelines), create_infos)
    if !isempty(info_uncached)
        (pipelines, _) = unwrap(create_graphics_pipelines(cache.device, info_uncached; allocator, pipeline_cache))
        foreach(zip(info_uncached, pipelines)) do (info, pipeline)
            insert!(cache.pipelines, info, pipeline)
        end
    end
    map(Base.Fix1(getindex, cache.pipelines), create_infos)
end

struct PipelineLayoutCache
    device::Device
    layouts::Dictionary{PipelineLayoutCreateInfo,PipelineLayout}
end

PipelineLayoutCache(device) = PipelineLayoutCache(device, Dictionary{PipelineLayoutCreateInfo,PipelineLayout}())

@forward PipelineLayoutCache.layouts Base.getindex, Base.haskey, Base.insert!

function Base.get!(cache::PipelineLayoutCache, info::PipelineLayoutCreateInfo; allocator = C_NULL)
    if haskey(cache, info)
        cache[info]
    else
        layout = unwrap(create_pipeline_layout(cache.device, info; allocator))
        insert!(cache, info, layout)
        layout
    end
end