text:
export midpoint, intermediate_point, destination_point, intersection_point,
max_latitude, opposite_point, closest_point_to_pole

"""
    midpoint(pos₁::Point, pos₂::Point)

Return the half-way point `midpoint` [deg] on the great circle line between
the positions `pos₁` and `pos₂` [deg] on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function midpoint(pos₁::Point, pos₂::Point)
    Δpos = pos₂ - pos₁
    Bx = cosd(pos₂.ϕ) * cosd(Δpos.λ)
    By = cosd(pos₂.ϕ) * sind(Δpos.λ)
    ϕₘ = atand(sind(pos₁.ϕ) + sind(pos₂.ϕ), √((cosd(pos₁.ϕ) + Bx)^2 + By^2))
    λₘ = pos₁.λ + atand(By, cosd(pos₁.ϕ) + Bx)
    return normalize(Point(ϕₘ, λₘ))
end

"""
    midpoint(section::RouteSection)

Return the half-way point `midpoint` [deg] on the route section on a sphere.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
midpoint(section::RouteSection) = midpoint(section.pos₁, section.pos₂)

"""
    intermediate_point(pos₁::Point, pos₂::Point[, fraction::Float64 = 0.5])

Return the `intermediate_point` [deg] at any `fraction` along the great circle
path between two points with positions `pos₁` and `pos₂` [deg]. The fraction along the
great circle route is such that `fraction` = 0.0 is at `pos₁` [deg] and `fraction` 1.0
is at `pos₂` [deg].

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function intermediate_point(pos₁::Point, pos₂::Point, fraction::Float64 = 0.5)
    δ = angular_distance(pos₁, pos₂)
    a = sind((1.0 - fraction) * δ) / sind(δ)
    b = sind(fraction * δ) / sind(δ)
    x = a * cosd(pos₁.ϕ) * cosd(pos₁.λ) + b * cosd(pos₂.ϕ) * cosd(pos₂.λ)
    y = a * cosd(pos₁.ϕ) * sind(pos₁.λ) + b * cosd(pos₂.ϕ) * sind(pos₂.λ)
    z = a * sind(pos₁.ϕ) + b * sind(pos₂.ϕ)
    ϕᵢ = atand(z, √(x*x + y*y))
    λᵢ = atand(y, x)
    return Point(ϕᵢ, λᵢ)
end

"""
    intermediate_point(section::RouteSection [, fraction::Float64 = 0.5])

Return the `intermediate_point` (Point) at any `fraction` along the route
section. The fraction along the great circle route is such that `fraction` =
0.0 is at `pos₁` and `fraction` 1.0 is at `pos₂`.

Source: www.movable-type.co.uk/scripts/latlong.html
"""
intermediate_point(section::RouteSection, fraction::Float64 = 0.5) =
intermediate_point(section.pos₁, section.pos₂, fraction)

"""
    destination_point(start_pos::Point, distance::Float64, bearing::Float64[,
     radius::Float64=Rₑ_m])

Given a `start_pos` [deg], initial `bearing` (clockwise from North) [deg],
`distance` [m], and the earth `radius` [m] calculate the `destination_point`
[deg] travelling along a (shortest distance) great circle arc.
(The distance must use the same radius as the radius of the earth.)

Source: www.movable-type.co.uk/scripts/latlong.html
"""
function destination_point(start_pos::Point, distance::Float64,
    bearing::Float64, radius::Float64=Rₑ_m)
    δ = rad2deg(distance / radius)  # Angular distance
    ϕ₂ = asind(sind(start_pos.ϕ) * cosd(δ) + cosd(start_pos.ϕ) * sind(δ) *
    cosd(bearing))
    λ₂ = start_pos.λ + atand(sind(bearing) * sind(δ) * cosd(start_pos.ϕ),
    cosd(δ) - sind(start_pos.ϕ) * sind(ϕ₂))
    return Point(ϕ₂, λ₂)
end

"""
    intersection_point(pos₁::Point, pos₂::Point, bearing₁₃::Float64, bearing₂₃::Float64)

Return the intersection point `pos₃` [deg] of two great circle paths given two
start points [deg] `pos₁` and `pos₂` [deg] and two bearings [deg] from `pos₁` to
`pos₃` [deg], and from `pos₂` to `pos₃` [deg].

Under certain circumstances the results can be an ∞ or *ambiguous solution*.

Source: edwilliams.org/avform.htm
"""
function intersection_point(pos₁::Point, pos₂::Point, bearing₁₃::Float64,
    bearing₂₃::Float64)
    Δpos = pos₂ - pos₁
    δ₁₂ = 2asind(√(sind(Δpos.ϕ / 2)^2 + cosd(pos₁.ϕ) * cosd(pos₂.ϕ) * sind(Δpos.λ/2)^2))
    if sind(δ₁₂) * cosd(pos₁.ϕ) == 0.0 || sind(δ₁₂) * cosd(pos₂.ϕ) == 0.0
        return Point(NaN, NaN)
    end
    θₐ = acosd((sind(pos₂.ϕ) - sind(pos₁.ϕ) * cosd(δ₁₂)) / (sind(δ₁₂) * cosd(pos₁.ϕ)))
    θᵦ = acosd((sind(pos₁.ϕ) - sind(pos₂.ϕ) * cosd(δ₁₂)) / (sind(δ₁₂) * cosd(pos₂.ϕ)))
    local θ₁₂, θ₂₁
    if sind(Δpos.λ) > 0.0
        θ₁₂ = θₐ
        θ₂₁ = 360.0 - θᵦ
    else
        θ₁₂ = 360.0 - θₐ
        θ₂₁ = θᵦ
    end
    α₁ = (bearing₁₃ - θ₁₂ + 180.0)%(360.0) - 180.0
    α₂ = (θ₂₁ - bearing₂₃ + 180.0)%(360.0) - 180.0

    if sind(α₁) ≈ 0.0 && sind(α₂) ≈ 0.0
        return Point(Inf, Inf)  # infinity of intersections
    elseif sind(α₁) * sind(α₂) < 0.0
        # return Point_rad(Inf, Inf)  # intersection ambiguous
        return Point(NaN, NaN)
    else
        α₃ = acosd(-cosd(α₁) * cosd(α₂) + sind(α₁) * sind(α₂) * cosd(δ₁₂))
        δ₁₃ = atand(sind(δ₁₂) * sind(α₁) * sind(α₂), cosd(α₂) + cosd(α₁) * cosd(α₃))
        ϕ₃ = asind(sind(pos₁.ϕ) * cosd(δ₁₃) + cosd(pos₁.ϕ) * sind(δ₁₃) * cosd(bearing₁₃))
        Δλ₁₃ = atand(sind(bearing₁₃) * sind(δ₁₃) * cosd(pos₁.ϕ), cosd(δ₁₃)
        - sind(pos₁.ϕ) * sind(ϕ₃))
        λ₃ = pos₁.λ + Δλ₁₃
        return Point(ϕ₃, λ₃)
    end
end

function intersection_point(pos₁::Point, pos₂::Point, pos₃::Point, pos₄::Point)
    inter_pnt = intersection_point(pos₁, pos₃, bearing(pos₁, pos₂),
    bearing(pos₃, pos₄))
    if isinf(inter_pnt.ϕ)
        return inter_pnt
    elseif distance(pos₁, pos₂) + tolerance_m ≥ distance(pos₁, inter_pnt) &&
        distance(pos₃, pos₄) + tolerance_m ≥ distance(pos₃, inter_pnt)
        return inter_pnt
    else
        return Point(NaN, NaN)
    end
end

"""
    intersection_point(pos₁::Point, pos₂::Point, polygon::Vector{Point{Float64}})

Return the intersection point `pos₃` [deg] of a great circle section with a
polygon. The great circle section is defined by two points [deg] `pos₁` and
`pos₂` [deg].
"""
function intersection_point(pos₁::Point, pos₂::Point,
    polygon::Vector{Point{Float64}})
    min_distance = Inf
    inter_point = Point(NaN, NaN)
    current_point = Point(NaN, NaN)
    poly₁ = polygon[1]
    for poly₂ in polygon[2:end]
        current_point = intersection_point(pos₁, pos₂, poly₁, poly₂)
        if current_point.ϕ ∉ [Inf, NaN]
            current_distance = distance(pos₁, current_point)
            if min_distance > current_distance
                min_distance = current_distance
                inter_point = current_point
            end
        end
        poly₁ = poly₂
    end
    return inter_point
end

"""
    intersection_point(pos₁::Point, pos₂::Point, airspace::Airspace)

Return the intersection point `pos₃` [deg] of a great circle section with an
`Airspace`. The great circle section is defined by two points [deg] `pos₁` and
`pos₂` [deg].
"""
function intersection_point(pos₁::Point, pos₂::Point,
    airspace::Airspace)
    inter_point = Point(NaN, NaN)
    if isinside(pos₁, airspace.bounding_box) ||
        isinside(pos₂, airspace.bounding_box)
        inter_point = intersection_point(pos₁, pos₂, airspace.polygon)
    else
        inter_box_point = intersection_point(pos₁, pos₂, airspace.bounding_box)
        if !isnan(inter_box_point.ϕ)
            inter_point = intersection_point(pos₁, pos₂, airspace.polygon)
        end
    end
    return inter_point
end

#TODO Source for longitude
"""
closest_point_to_pole(starting_point::Point, bearing::Float64)

Given a `starting_point` [deg] and initial `bearing` [deg] calculate the
closest point [deg] to the next pole encountered.

Source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function closest_point_to_pole(starting_point::Point, bearing::Float64)
    pnt = normalize(starting_point)
    brg = normalize(bearing, -180.0, 180.0)
    max_lat = max_latitude(pnt.ϕ, brg)
    max_lon = pnt.λ + atand(1.0, tand(brg)sind(pnt.ϕ))
    return normalize(Point(max_lat*(-90.0≤brg≤90.0 ? 1.0 : -1.0),
    max_lon-180.0*(brg<0.0)))
end

"""
max_latitude(latitude::Float64, bearing::Float64)

Using Clairaut's formula it is possible to calculate the maximum latitude [deg]
of a great circle path, given a `bearing` [deg] and `latitude` [deg] on the
great circle.

The minimum latitude is -max_latitude

Source: https://www.movable-type.co.uk/scripts/latlong.html
"""
function max_latitude(latitude::Float64, bearing::Float64)
    return acosd(abs(sind(bearing)cosd(latitude)))
end

"""
opposite_point(point::Point)

The point at the opposite site of the Sphere.
"""
function opposite_point(point::Point)
    normalize(Point(-point.ϕ, point.λ+180.0))
end


text:
using Printf

mutable struct VLQ
    quant::Vector{UInt8}
end

function VLQ(n::T) where T <: Integer
    quant = UInt8.(digits(n, 128))
    @inbounds for i in 2:length(quant) quant[i] |= 0x80 end
    VLQ(reverse(quant))
end

import Base.UInt64
function Base.UInt64(vlq::VLQ)
    quant = reverse(vlq.quant)
    n = shift!(quant)
    p = one(UInt64)
    for i in quant
        p *= 0x80
        n += p * ( i & 0x7f)
    end
    return n
end

const test = [0x00200000, 0x001fffff, 0x00000000, 0x0000007f,
              0x00000080, 0x00002000, 0x00003fff, 0x00004000,
              0x08000000, 0x0fffffff]

for i in test
    vlq = VLQ(i)
    j = UInt(vlq)
    @printf "0x%-8x => [%-25s] => 0x%x\n" i join(("0x" * hex(r, 2) for r in vlq.quant), ", ") j
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## MCMC functions for gaia.
##
##

### generate eps, min_neightbor and min_cluster
###
### peps , pminnei , pmincl are the parameter of the prior density for the DBSCAN parameters

function theta(p::GaiaClustering.abc)
    peps         = TruncatedNormal(p.epsmean, p.epsdisp, 0.1 , 1000)
    pminnei      = TruncatedNormal(p.min_nei, p.ncoredisp, 1 , 1000.)
    pmincl       = TruncatedNormal(p.min_cl, p.ncoredisp, 1 , 1000.)

    e = rand(peps)          ; pe = pdf(peps,e)
    n = trunc(Int,rand(pminnei))  ; pn = pdf(pminnei, n)
    c = trunc(Int,rand(pmincl))   ; pc = pdf(pmincl, c)

    ptotal = pe*pn*pc
    params = GaiaClustering.model(e, n, c)
    return(params, ptotal)
end

## iterate with random walk and yield the probability
##
function thetaiter(θi::GaiaClustering.model , p::GaiaClustering.abc)
    let
        peps         = TruncatedNormal(p.epsmean, p.epsdisp, 0.1 , 1000)
        pminnei      = TruncatedNormal(p.min_nei, p.ncoredisp, 1 , 1000.)
        pmincl       = TruncatedNormal(p.min_cl, p.ncoredisp, 1 , 1000.)

        ## Random Walk settings..
        eps_rw    = TruncatedNormal(θi.eps,0.25 , 0.1, 1000.)
        minnei_rw = TruncatedNormal(θi.min_nei , 4  , 1,  1000.)
        mincl_rw  = TruncatedNormal(θi.min_cl  , 4  , 1, 1000.)

        new_e = 0.
        new_mn = 0
        new_mcl = 0

        iternotfound =true
        while iternotfound
            new_e = rand(eps_rw)
            new_mn =  trunc(Int,rand(minnei_rw))
            new_mcl = trunc(Int,rand(mincl_rw))

            if new_e > 0 && new_mn > 0 && new_mcl > 0
                iternotfound = false
            end
        end

        pe = pdf(peps , new_e)
        pn = pdf(pminnei, new_mn)
        pc = pdf(pmincl, new_mcl)
        ptotal = pe * pn * pc

        params = GaiaClustering.model(new_e, new_mn, new_mcl)
        return(params, ptotal)
    end
end

#### function to adapt the minimum condition to evaluate the  solution
###
function check_qminqstar(df::GaiaClustering.Df, dfcart::GaiaClustering.Df,
        params::GaiaClustering.abc, minimumQ, minstars)
    let
        new_minq = minimumQ
        new_minstars = minstars
        notfound = true
        mingoodsolution = 50

        println("### Checking the minQ and minStars conditions...")
        println("### Minimum good solutions $mingoodsolution")
        while notfound
            goodsolutions = 0
            for i in 1:params.nburnout
                mi, probi = theta(params)
                qres , nstars = find_clusters(df, dfcart, mi)
                if qres > new_minq && nstars >= new_minstars
                    goodsolutions += 1
                end
                if goodsolutions > mingoodsolution
                    notfound = false
                    return(new_minq, new_minstars)
                end
            end
            if notfound
                new_minq *= 0.9
                new_minstars = trunc(Int, 0.9 * new_minstars)
                println("### MinQ not reached yet... testing with $new_minq")
            end

            if new_minstars == 0
                println("### Problem with minStars, return(0.5,5)")
                return(0.5,5)
            end
        end
        return(new_minq, new_minstars)
    end
end

## ABC MCMC (following Weyan et al. 2013)
## first simple scheme for testing
##
function abc_mcmc_dbscan(df::GaiaClustering.Df, dfcart::GaiaClustering.Df, params::GaiaClustering.abc)
    let
        println("## ABC/MCMC for DBSCAN...")
        mci = GaiaClustering.mc(zeros(Float64,0),zeros(Int32,0),zeros(Int32,0) , zeros(Float64,0), zeros(Int32,0))

        initial = true
        th = []
        mi = GaiaClustering.model(0.0,0,0)
        micurrent = GaiaClustering.model(0.0,0,0)
        probi = 0.
        probcurrent = 0.

        minimumQ = params.minQ
        minstars = params.minstars
        niter = params.niter
        nburn = params.nburnout

        println("### Minimum Q : $minimumQ")
        println("### Minimum nstars : $minstars")
        minimumQ , minstars = check_qminqstar(df, dfcart, params, minimumQ, minstars)
        #params.minQ = minimuQ
        #params.minstars = minstars
        println("### Minimum Q : $minimumQ")
        println("### Minimum nstars : $minstars")
        println("### Lower conditions lead to lower quality solutions...")
        println("###")

        while initial
            mi, probi = theta(params)
            qres , nstars = find_clusters(df, dfcart, mi)
            if qres > minimumQ && nstars >= minstars
                println("### init done ...")
                initial = false
                push!(mci.eps, mi.eps)
                push!(mci.mne, mi.min_nei)
                push!(mci.mcl, mi.min_cl)
                push!(mci.qc,  qres)
                push!(mci.qn,  nstars)
            end
        end

        ministats(100, df, dfcart, mi , params)

        nchain = 1
        loopAgain = true
        burndone = false

        while loopAgain
            micurrent, probcurrent = thetaiter(mi , params)
            qres , nstars = find_clusters(df, dfcart, micurrent)

            if qres > minimumQ && nstars >= minstars
            ### Metropolis-Hasting
                α = probcurrent / probi
                if α > rand()
                    mi = micurrent
                    probi = probcurrent
                    nchain += 1
                    if (nchain%500 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burn-in done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                else
                    nchain += 1
                    if (nchain%500 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burn-in done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                end
            end
        end
        println("## ABC/MCMC done")
        println("##")
        return(mci)
    end
end

## mini stats over Qc and Qn
function ministats(niter::Int, df::GaiaClustering.Df, dfcart::GaiaClustering.Df,mi::GaiaClustering.model, params::GaiaClustering.abc)
    println("### mini stats...")
    qcmini = []
    qnmini = []
    for i in 1:niter
        mitest, probtest = thetaiter(mi , params)
        qtest , ntest = find_clusters(df, dfcart, mitest)
        push!(qcmini,qtest)
        push!(qnmini,ntest)
    end
    println("### Qc : ",mean(qcmini))
    println("### Qn : ",mean(qnmini))

    return(0)
end

##############################################
#### MCMC FULL for dbscan
######
function getDfcartnorm(dfcart::GaiaClustering.Df, mc::modelfull)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [mc.w3d, mc.wvel, mc.whrd]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false, false)
    return(dfcartnorm)
end

## optimize dbscan parameters with ABC/MCMC method
##
function abc_mcmc_dbscan_full2(dfcart::GaiaClustering.Df, params::GaiaClustering.meta)
    let
        Random.seed!()
        println("## ABC/MCMC for DBSCAN FULL (parameters+weighting)...")
        println("## ABC/MCMC v2")

        mci = mcfull(zeros(Float64,0),zeros(Int32,0),zeros(Int32,0) , zeros(Float64,0), zeros(Int32,0),
        zeros(Float64,0), zeros(Float64,0), zeros(Float64,0))

        initial = true
        FLAG= -1            ## everything normal
        th = []
        mi = modelfull(0.0,0,0, 0., 0., 0.)
        micurrent = modelfull(0.0,0,0, 0., 0., 0.)
        probi = 0.
        probcurrent = 0.

        minimumQ= params.minQc
        minstars= params.minQn
        maxstars= params.maxQn
        niter= params.nchain             ## number of chains
        nburn= params.nburnout

        maxiter= params.maxiter          ## twice, for init and normal iteration

        println("### Burn-in : $nburn")
        println("### Chains  : $niter")
        println("### Minimum Qc : $minimumQ")
        println("### Minimum Qn : $minstars")
        println("### Maximum Qn : $maxstars")
        println("### Maximum iterations: $maxiter")
        println("###")

        minimumQ , minstars = check_qminqstar_full2(dfcart, params)
        @printf("### Minimum Qc: %3.3f \n",minimumQ)
        println("### Minimum Qn: $minstars")
        if minstars < params.forcedminstars
            minstars = params.forcedminstars
            println("### Minimum Qn forced to $minstars")
        end

        iter= 0
        while initial
            mi, probi = theta_full(params)
            dfcartnorm = getDfcartnorm(dfcart , mi)
            qres , nstars = find_clusters2(dfcartnorm, dfcart, mi, params)
            if qres > minimumQ && nstars >= minstars && nstars <= maxstars
                println("### init done ...")
                initial = false
                push!(mci.eps, mi.eps)
                push!(mci.mne, mi.min_nei)
                push!(mci.mcl, mi.min_cl)
                push!(mci.w3d, mi.w3d)
                push!(mci.wvel, mi.wvel)
                push!(mci.whrd, mi.whrd)
                push!(mci.qc,  qres)
                push!(mci.qn,  nstars)
            end
            iter += 1
            if iter > maxiter
                println("### Maximum iteration reached, no solution for init...")
                FLAG= 2      ## init not performed completely
                return(mci, iter,FLAG)
            end
        end

        ministats_full(100, dfcart, mi , params)

        nchain= 1
        loopAgain= true
        burndone= false
        tstart= now()

        iter= 0
        while loopAgain
            micurrent, probcurrent = thetaiter_full(mi , params)
            dfcartnorm = getDfcartnorm(dfcart , micurrent)
            qres , nstars = find_clusters2(dfcartnorm, dfcart, micurrent, params)

            if qres > minimumQ && nstars >= minstars && nstars <= maxstars
            ### Metropolis-Hasting
                α = probcurrent / probi
                if α > rand()
                    mi = micurrent
                    probi = probcurrent
                    nchain += 1
                    if (nchain%1000 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burn-in done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter && burndone loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.w3d, mi.w3d)
                        push!(mci.wvel, mi.wvel)
                        push!(mci.whrd, mi.whrd)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                else
                    nchain += 1
                    if (nchain%1000 == 0) println("### chain:",nchain) end
                    if nchain > nburn && !burndone println("### burn-in done...") ; nchain = 0 ; burndone = true end
                    if nchain > niter && burndone loopAgain = false end
                    if burndone
                        push!(mci.eps, mi.eps)
                        push!(mci.mne, mi.min_nei)
                        push!(mci.mcl, mi.min_cl)
                        push!(mci.w3d, mi.w3d)
                        push!(mci.wvel, mi.wvel)
                        push!(mci.whrd, mi.whrd)
                        push!(mci.qc , qres)
                        push!(mci.qn,  nstars)
                    end
                end
            end

            iter += 1
            if (iter%10000 == 0)
                titer= now()
                duration= Dates.value(titer-tstart) / (1000*3600)
                meanTime= duration / nchain
                eta= @sprintf("%3.3f",meanTime * (niter-nchain))
                println("### iteration: $iter (ETA:$eta h)")
            end
            if iter > maxiter
                println("### Maximum iteration reached, current solution returned...")
                if burndone   FLAG= 3 end    ## nchain not reached before maxiter but burnin done...
                if !burndone  FLAG= 4 end    ## burn-in not performed
                return(mci, iter, FLAG)
            end
        end
        println("## ABC/MCMC (full) done")
        println("##")
        return(mci,  iter, FLAG)
    end
end

## mini stats over Qc and Qn
function ministats_full(niter::Int, dfcart::GaiaClustering.Df, mi::modelfull, params::meta)
    println("### mini stats...")
    qcmini = []
    qnmini = []
    for i in 1:niter
        mitest, probtest = thetaiter_full(mi , params)
        dfcartnorm = getDfcartnorm(dfcart , mitest)

        qtest , ntest = find_clusters2(dfcartnorm, dfcart, mitest, params)
        push!(qcmini,qtest)
        push!(qnmini,ntest)
    end
    @printf("### Qc : %3.3f \n",mean(qcmini))
    @printf("### Qn : %3.3f \n",mean(qnmini))

    return(0)
end

function theta_full(p::meta)
    peps         = truncated(Normal(p.epsmean, p.epsdisp), 0.1 , 1000)
    pminnei      = truncated(Normal(p.min_nei, p.ncoredisp), 1 , 1000.)
    pmincl       = truncated(Normal(p.min_cl, p.ncoredisp), 1 , 1000.)
    pw3d         = truncated(Normal(p.w3dmean , p.w3ddisp), 0.1, 1000)
    pwvel        = truncated(Normal(p.wvelmean , p.wveldisp), 0.1, 1000)
    pwhrd        = truncated(Normal(p.whrdmean , p.whrddisp), 0.1, 1000)

    e = rand(peps)          ; pe = pdf(peps,e)
    n = trunc(Int,rand(pminnei))  ; pn = pdf(pminnei, n)
    c = trunc(Int,rand(pmincl))   ; pc = pdf(pmincl, c)
    w3 = rand(pw3d)          ; pw3 = pdf(pw3d,w3)
    wv = rand(pwvel)         ; pwv = pdf(pwvel,wv)
    wh = rand(pwhrd)         ; pwh = pdf(pwhrd,wh)

    ptotal = pe*pn*pc*pw3*pwv*pwh
    params = modelfull(e, n, c , w3, wv, wh)
    return(params, ptotal)
end

## iterate with random walk and yield the probability
##
function thetaiter_full(θi::modelfull , p::meta)
    let
        peps         = truncated(Normal(p.epsmean, p.epsdisp), 0.1, 1000)
        pminnei      = truncated(Normal(p.min_nei, p.ncoredisp), 1 , 1000.)
        pmincl       = truncated(Normal(p.min_cl, p.ncoredisp), 1 , 1000.)
        pw3d         = truncated(Normal(p.w3dmean , p.w3ddisp), 0.1, 1000)
        pwvel        = truncated(Normal(p.wvelmean , p.wveldisp), 0.1, 1000)
        pwhrd        = truncated(Normal(p.whrdmean , p.whrddisp), 0.1, 1000)

        eps_rw    = truncated(Normal(θi.eps,p.eps_rw) , 0.1, 1000.)
        minnei_rw = truncated(Normal(θi.min_nei , p.minnei_rw) , 1,  1000.)
        mincl_rw  = truncated(Normal(θi.min_cl  , p.mincl_rw)  , 1, 1000.)
        w3d_rw    = truncated(Normal(θi.w3d,  p.w3d_rw), 0.1, 1000)
        wvel_rw   = truncated(Normal(θi.wvel, p.wvel_rw), 0.1, 1000)
        whrd_rw   = truncated(Normal(θi.whrd, p.whrd_rw), 0.1, 1000)

        new_e = 0.
        new_mn = 0
        new_mcl = 0
        new_w3d = 0.
        new_wvel = 0.
        new_whrd = 0.

        iternotfound = true
        while iternotfound
            new_e    = rand(eps_rw)
            new_mn   = trunc(Int,rand(minnei_rw))
            new_mcl  = trunc(Int,rand(mincl_rw))
            new_w3d  = rand(w3d_rw)
            new_wvel = rand(wvel_rw)
            new_whrd = rand(whrd_rw)

            ### adding the condition on minclusterand > minnei..
            if new_e > 0 && new_mn > 0 && new_mcl > 0 && new_mn <= new_mcl
                iternotfound = false
            end
        end

        pe = pdf(peps , new_e)
        pn = pdf(pminnei, new_mn)
        pc = pdf(pmincl, new_mcl)
        pw3 = pdf(pw3d, new_w3d)
        pwv = pdf(pwvel, new_wvel)
        pwh = pdf(pwhrd, new_whrd)

        ptotal = pe * pn * pc * pw3 * pwv * pwh
        params = modelfull(new_e, new_mn, new_mcl , new_w3d, new_wvel, new_whrd)
        return(params, ptotal)
    end
end

## check maximum iterations
function check_qminqstar_full2(dfcart::GaiaClustering.Df, params::GaiaClustering.meta)
    let
        new_minq = params.minQc
        new_minstars = params.minQn
        notfound = true

        ## standard sofar 50
        mingoodsolution = params.mingoodsolution
        niter = params.niterqminq
        maxiter= niter*30           ## cycle numbers * niter

        println("#### Checking the minQc and minQn conditions...")
        println("#### Minimum good solutions: $mingoodsolution")
        println("#### Number of iterations: $niter, maxiter: $maxiter")

        totaliter= 0
        while notfound
            goodsolutions = 0
            for i in 1:niter
                # println("-- $i")
                mi, probi = theta_full(params)
                # println(mi)
                # println(probi)
                dfcartnorm = getDfcartnorm(dfcart , mi)
                qres , nstars = find_clusters2(dfcartnorm, dfcart, mi, params)
                if qres > new_minq && nstars >= new_minstars
                    goodsolutions += 1
                end
                if goodsolutions > mingoodsolution
                    notfound = false
                    return(new_minq, new_minstars)
                end
            end
            totaliter += niter
            if notfound
                new_minq *= 0.95
                new_minstars = trunc(Int, 0.95 * new_minstars)
                new_minqstr= @sprintf("%3.3f", new_minq)
                println("#### MinQ not reached yet... testing with $new_minqstr")
            end

            if new_minstars == 0 || totaliter > maxiter
                println("#### Maximum iteration reached..,")
                if new_minstars < 5
                    new_minstars= 5
                end
                if new_minq < 0.5
                    new_minq= 0.5
                end
                println("#### Maximum iterations reached.. minimums set to $new_minq, $new_minstars")
                return(new_minq, new_minstars)
            end
        end
        return(new_minq, new_minstars)
    end
end

function extraction_mcmc(mc, votname)
    epsm = median(mc.eps)
    epsd = std(mc.eps)
    mneim = median(mc.mne)
    mneid = std(mc.mne)
    mclm = median(mc.mcl)
    mcld = std(mc.mcl)
    qcm = median(mc.qc)
    qnm = median(mc.qn)
    qcd = std(mc.qc)
    qnd = std(mc.qn)
    w3dm = median(mc.w3d)
    w3dd = std(mc.w3d)
    wvelm = median(mc.wvel)
    wveld = std(mc.wvel)
    whrdm = median(mc.whrd)
    whrdd = std(mc.whrd)

    @printf("## DBSCAN/MCMC stats: \n")
    @printf("### ϵ : %3.3f +/- %3.3f \n", epsm, epsd)
    @printf("### min_nei  : %3.1f +/- %3.3f \n", mneim, mneid)
    @printf("### min_clus  : %3.1f +/- %3.3f \n", mclm, mcld)
    @printf("### W3d  : %3.3f +/- %3.3f \n", w3dm, w3dd)
    @printf("### Wvel  : %3.3f +/- %3.3f \n" , wvelm, wveld)
    @printf("### Whrd  : %3.3f +/- %3.3f \n", whrdm, whrdd)
    @printf("### Qn  : %3.3f +/- %3.3f \n",qnm, qnd)
    @printf("### Qc  : %3.3f +/- %3.3f \n",qcm, qcd)
    @printf("##")

    res = DataFrame(votname=votname, epsm = epsm, epsd=epsd, mneim=mneim,mneid=mneid,mclm=mclm,mcld=mcld,
            qcm=qcm,qcd=qcd, qnm=qnm,qnd=qnd,
            w3dm=w3dm,w3dd=w3dd,wvelm=wvelm,wveld=wveld,whrdm=whrdm,whrdd=whrdd)
    return(res)
end

### create DataFrame for the chaine
###
function create_DFchain(mc, votname, cycle)
    df= DataFrame(votname=votname, cycle=cycle, eps=mc.eps , mne=mc.mne, mcl=mc.mcl, w3d=mc.w3d, wvel=mc.wvel, whrd=mc.whrd, qc=mc.qc, qn=mc.qn)
    return(df)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Modflow6_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Modflow6")
JLLWrappers.@generate_main_file("Modflow6", UUID("81230ed9-c869-5493-bf63-2c1575eaf49b"))
end  # module Modflow6_jll


text:
# This file contains the API that users interact with.  Everything in here should
# depend only on `brew()` calls.

"""
`brew(cmd::Cmd; no_stderr=false, no_stdout=false, verbose=false, force=false, quiet=false)`

Run command `cmd` using the configured brew binary, optionally suppressing
stdout and stderr, and providing flags such as `--verbose` to the brew binary.
"""
function brew(cmd::Cmd; no_stderr=false, no_stdout=false, verbose::Bool=false, force::Bool=false, quiet::Bool=false)
    cmd = add_flags(`$brew_exe $cmd`, Dict(`--verbose` => verbose, `--force` => force, `--quiet` => quiet))

    if no_stderr
        cmd = pipeline(cmd, stderr=devnull)
    end
    if no_stdout
        cmd = pipeline(cmd, stdout=devnull)
    end
    return run(cmd)
end

"""
`brewchomp(cmd::Cmd; no_stderr=false, no_stdout=false, verbose=false, force=false, quiet=false))`

Run command `cmd` using the configured brew binary, optionally suppressing
stdout and stderr, and providing flags such as `--verbose` to the brew binary.

This function uses `readchomp()`, as opposed to `brew()` which uses `run()`
"""
function brewchomp(cmd::Cmd; no_stderr=false, no_stdout=false, verbose::Bool=false, force::Bool=false, quiet::Bool=false)
    cmd = add_flags(`$brew_exe $cmd`, Dict(`--verbose` => verbose, `--force` => force, `--quiet` => quiet))

    if no_stderr
        cmd = pipeline(cmd, stderr=devnull)
    end
    if no_stdout
        cmd = pipeline(cmd, stdout=devnull)
    end
    return readchomp(cmd)
end

"""
`update(;verbose::Bool=false)`

Runs `brew update` to update Homebrew itself and all taps.  Then runs `upgrade()`
to upgrade all formulae that have fallen out of date.
"""
function update(;verbose::Bool=false)
    # Just run `brew update`
    brew(`update`; force=true, verbose=verbose)

    # Ensure we're on the right tag
    update_tag(;verbose=verbose)

    # Finally, upgrade outdated packages.
    upgrade(;verbose=verbose)
end

"""
`prefix()`

Returns `brew_prefix`, the location where all Homebrew files are stored.
"""
function prefix()
    return brew_prefix
end

"""
`prefix(pkg::Union{AbstractString,BrewPkg})`

Returns the prefix for a particular package's latest installed version.
"""
function prefix(pkg::StringOrPkg) end

function prefix(name::AbstractString)
    path, tap_path = formula_tap(name)
    return joinpath(brew_prefix, "Cellar", path, info(name).version)
end

function prefix(pkg::BrewPkg)
    prefix(pkg.name)
end

"""
`list()`

Returns a list of all installed packages as a `Vector{BrewPkg}`
"""
function list()
    # Get list of fully-qualified installed packages
    names = brewchomp(`list --full-name`)

    # If we've got something, then return info on all those names
    if !isempty(names)
        return info(split(names,"\n"))
    end
    return BrewPkg[]
end

"""
`outdated()`

Returns a list of all installed packages that are out of date as a `Vector{BrewPkg}`
"""
function outdated()
    json_str = brewchomp(`outdated --json=v1`)
    if isempty(json_str)
        return BrewPkg[]
    end
    brew_outdated = JSON.parse(json_str)

    outdated_pkgs = BrewPkg[info(pkg["name"]) for pkg in brew_outdated]
    return outdated_pkgs
end

"""
`refresh!(;verbose=false)`

Forcibly remove all packages and add them again.  This should only be used to
fix a broken installation, normal operation should never need to use this.
"""
function refresh!(;verbose=false)
    pkg_list = list()
    for pkg in pkg_list
        rm(pkg,verbose=verbose)
    end
    for pkg in pkg_list
        add(pkg,verbose=verbose)
    end
end

"""
`upgrade(;verbose::Bool=false)`

Iterate over all packages returned from `outdated()`, removing the old version
and adding a new one.  Note that we do not simply call `brew upgrade` here, as
we have special logic inside of `add()` to install from our tap before trying to
install from mainline Homebrew.
"""
function upgrade(;verbose::Bool=false)
    # We have to manually upgrade each package, as `brew upgrade` will pull from mxcl/master
    for pkg in outdated()
        rm(pkg; verbose=verbose)
        add(pkg; verbose=verbose)
    end
end

const json_cache = Dict{String,Dict{AbstractString,Any}}()
"""
`json(names::Vector{AbstractString})`

For each package name in `names`, return the full JSON object for `name`, the
result of `brew info --json=v1 \$name`, stored in a dictionary keyed by the names
passed into this function. If `brew info` fails, throws an error. If `brew info`
returns an empty object "[]", that object is represented by an empty dictionary.

Note that running `brew info --json=v1` is somewhat expensive, so we cache the
results in a global dictionary, and batching larger requests with this function
similarly increases performance.
"""
function json(names::Vector{T}) where {T <: AbstractString}
    # This is the dictionary of responses we'll return
    objs = Dict{String,Dict{AbstractString,Any}}()

    # Build list of names we have to ask for, eschewing asking for caching when we can
    ask_names = String[]
    for name in names
        if haskey(json_cache, name)
            objs[name] = json_cache[name]
        else
            push!(ask_names, name)
        end
    end

    # Now ask for all these names if we have any
    if !isempty(ask_names)
        # Tap the necessary tap for each name we're asking about, if there is one
        for name in ask_names
            path, tap_path = formula_tap(name)

            if !isempty(tap_path)
                tap(tap_path)
            end
        end
        try
            jdata = JSON.parse(brewchomp(Cmd(String["info", "--json=v1", ask_names...])))
            for idx in 1:length(jdata)
                json_cache[ask_names[idx]] = jdata[idx]
                objs[ask_names[idx]] = jdata[idx]
            end
        catch
            throw(ArgumentError("`brew info` failed for $(ask_names)!"))
        end
    end

    # Return these hard-won objects
    return objs
end

"""
`json(pkg::Union{AbstractString,BrewPkg})`

Return the full JSON object for `pkg`, the result of `brew info --json=v1 \$pkg`.
If `brew info` fails, throws an error.  If `brew info` returns an empty object,
(e.g. "[]"), this returns an empty Dict.

Note that running `brew info --json=v1` is somewhat expensive, so we cache the
results in a global dictionary, and batching larger requests with the vectorized
`json()` function similarly increases performance.
"""
function json(pkg::StringOrPkg) end

function json(name::AbstractString)
    return json([name])[name]
end

function json(pkg::BrewPkg)
    return json([fullname(pkg)])[fullname(pkg)]
end

"""
`info(names::Vector{String})`

For each name in `names`, returns information about that particular package name
as a BrewPkg.  This is our batched `String` -> `BrewPkg` converter.
"""
function info(names::Vector{T}) where {T <: AbstractString}
    # Get the JSON representations of all of these packages
    objs = json(names)

    infos = BrewPkg[]
    for name in names
        obj = objs[name]

        # First, get name and version
        obj_name = obj["name"]
        version = obj["versions"]["stable"]

        # Manually append the revision to the end of the version, as brew is wont to do
        if obj["revision"] > 0
            version *= "_$(obj["revision"])"
        end

        tap = dirname(obj["full_name"])
        if isempty(tap)
            tap = "Homebrew/core"
        end

        # Push that BrewPkg onto our infos object
        push!(infos, BrewPkg(obj_name, tap, version))
    end

    # Return the list of infos
    return infos
end

"""
`info(name::AbstractString)`

Returns information about a particular package name as a BrewPkg.  This is our
basic `String` -> `BrewPkg` converter.
"""
function info(name::AbstractString)
    return info([name])[1]
end

"""
`direct_deps(pkg::Union{AbstractString,BrewPkg}; build_deps::Bool=false)`

Return a list of all direct dependencies of `pkg` as a `Vector{BrewPkg}`
If `build_deps` is `true` include formula depependencies marked as `:build`.
"""
function direct_deps(pkg::StringOrPkg; build_deps::Bool=false) end

function direct_deps(name::AbstractString; build_deps::Bool=false)
    obj = json(name)

    # Iterate over all dependencies, removing optional and build dependencies
    dependencies = String[dep for dep in obj["dependencies"]]
    dependencies = filter(x -> !(x in obj["optional_dependencies"]), dependencies)

    if !build_deps
        dependencies = filter(x -> !(x in obj["build_dependencies"]), dependencies)
    end
    return info(dependencies)
end

function direct_deps(pkg::BrewPkg; build_deps::Bool=false)
    return direct_deps(fullname(pkg); build_deps=build_deps)
end

"""
`deps_tree(pkg::Union{AbstractString,BrewPkg}; build_deps::Bool=false)`

Return a dictionary mapping every dependency (both direct and indirect) of `pkg`
to a `Vector{BrewPkg}` of all of its dependencies.  Used in `deps_sorted()`.

If `build_deps` is `true` include formula depependencies marked as `:build`.
"""
function deps_tree(pkg::StringOrPkg; build_deps::Bool=false) end

function deps_tree(name::AbstractString; build_deps::Bool=false)
    # First, get all the knowledge we need about dependencies
    deptree = Dict{String,Vector{BrewPkg}}()

    pending_deps = direct_deps(name; build_deps=build_deps)
    completed_deps = Set(String[])
    while !isempty(pending_deps)
        # Temporarily move pending_deps over to curr_pending_deps
        curr_pending_deps = pending_deps
        pending_deps = BrewPkg[]

        # Iterate over these currently pending deps, adding new ones to pending_deps
        for pkg in curr_pending_deps
            deptree[fullname(pkg)] = direct_deps(pkg; build_deps=build_deps)
            push!(completed_deps, fullname(pkg))
            for dpkg in deptree[fullname(pkg)]
                if !(fullname(dpkg) in completed_deps)
                    push!(pending_deps, dpkg)
                end
            end
        end
    end

    return deptree
end

function deps_tree(pkg::BrewPkg; build_deps::Bool=false)
    return deps_tree(fullname(pkg); build_deps=build_deps)
end

"""
`insert_after_dependencies(tree::Dict, sorted_deps::Vector{BrewPkg}, name::AbstractString)`

Given a mapping from names to dependencies in `tree`, and a list of sorted
dependencies in `sorted_deps`, insert a new dependency `name` into `sorted_deps`
after all dependencies of `name`.  If a dependency of `name` is not already in
`sorted_deps`, then recursively add that dependency as well.
"""
function insert_after_dependencies(tree::Dict, sorted_deps::Vector{BrewPkg}, name::AbstractString)
    # First off, are we already in sorted_deps?  If so, back out!
    self_idx = findfirst(x -> (fullname(x) == name), sorted_deps)
    if self_idx != nothing
        return self_idx
    end

    # This is the index at which we will insert ourselves
    insert_idx = 1
    # Iterate over all dependencies
    for dpkg in tree[name]
        # Is this dependency already in the sorted_deps?
        idx = findfirst(x -> (fullname(x) == fullname(dpkg)), sorted_deps)

        # If the dependency is not already in this list, then recurse into it!
        if self_idx == nothing
            idx = insert_after_dependencies(tree, sorted_deps, fullname(dpkg))
        end

        # Otherwise, update insert_idx
        insert_idx = max(insert_idx, idx + 1)
    end

    # Finally, insert ourselves
    insert!(sorted_deps, insert_idx, info(name))
    return insert_idx
end

"""
`deps_sorted(pkg::Union{AbstractString,BrewPkg}; build_deps::Bool)`

Return a sorted `Vector{BrewPkg}` of all dependencies (direct and indirect) such
that each entry in the list appears after all of its own dependencies.

If `build_deps` is `true` include formula depependencies marked as `:build`.
"""
function deps_sorted(pkg::StringOrPkg; build_deps::Bool=false) end

function deps_sorted(name::AbstractString; build_deps::Bool=false)
    tree = deps_tree(name; build_deps=build_deps)
    sorted_deps = BrewPkg[]

    # For each package in the tree, insert it only after all of its dependencies
    # Just for aesthetic purposes, sort the keys by the number of dependencies
    # they have first, so that packages with few deps end up on top
    for name in sort(collect(keys(tree)), by=k-> length(tree[k]))
        insert_after_dependencies(tree, sorted_deps, name)
    end

    return sorted_deps
end

function deps_sorted(pkg::BrewPkg; build_deps::Bool=false)
    return deps_sorted(fullname(pkg); build_deps=build_deps)
end

"""
`add(pkg::Union{AbstractString,BrewPkg}; verbose::Bool=false, keep_translations::Bool=false)`

Install package `pkg` and all dependencies, using bottles only, unlinking any
previous versions if necessary, and linking the new ones in place. Will attempt
to install non-relocatable bottles from `Homebrew/core` by translating formulae
and forcing `cellar :any` into the formulae.

Automatically deletes all translated formulae before adding formulae and after,
unless `keep_translations` is set to `true`.
"""
function add(pkg::StringOrPkg; verbose::Bool=false, keep_translations::Bool=false) end

function add(name::AbstractString; verbose::Bool=false, keep_translations::Bool=false)
    # We do this so that things are as unambiguous and fresh as possible.
    # It's not a bad situation because translation is very fast.
    if !keep_translations
        delete_all_translated_formulae()
    end

    # Begin by translating all dependencies of `name`, including :build deps.
    # We need to translate :build deps so that we don't run into the diamond of death
    sorted_deps = AbstractString[]
    build_deps_only = AbstractString[]
    runtime_deps = [x.name for x in deps_sorted(name)]
    if verbose
        println("runtime_deps: $runtime_deps")
    end
    for dep in deps_sorted(name; build_deps=true)
        translated_dep = translate_formula(dep; verbose=verbose)

        # Collect the translated runtime_deps into sorted_deps
        if dep.name in runtime_deps
            push!(sorted_deps, translated_dep)
        else
            push!(build_deps_only, translated_dep)
        end
    end
    if verbose
        println("sorted_deps: $sorted_deps")
        println("build_deps only: $build_deps_only")
    end

    # Translate the actual formula we're interested in
    name = translate_formula(name; verbose=verbose)

    # Push `name` onto the end of `sorted_deps` so we have one list of things to install
    push!(sorted_deps, name)

    # Ensure that we are able to install relocatable bottles for every package
    non_relocatable = filter(x -> !has_relocatable_bottle(x), sorted_deps)

    if !isempty(non_relocatable)
        @warn("""
The following packages do not have relocatable bottles, installation may fail!
Please report these packages to https://github.com/JuliaLang/Homebrew.jl:
  $(join(non_relocatable, "\n  "))""")
    end

    # Get list of outdated packages and remove them all
    for dep in filter(pkg -> pkg.name in sorted_deps, Homebrew.outdated())
        unlink(dep; verbose=verbose)
    end

    # Install this package and all dependencies, in dependency order
    for dep in sorted_deps
        install_and_link(dep; verbose=verbose)
    end

    # Cleanup translated formula once we're done, unless asked not to
    if !keep_translations
        delete_all_translated_formulae()
    end
end

function add(pkg::BrewPkg; verbose::Bool=false, keep_translations::Bool=false)
    add(fullname(pkg); verbose=verbose, keep_translations=keep_translations)
end

"""
`install_and_link(pkg::Union{AbstractString,BrewPkg}; verbose=false)`

Installs, and links package `pkg`.  Used by `add()`.  Don't call manually
unless you really know what you're doing, as this doesn't deal with
dependencies, and so can trigger compilation when you don't want it to.
"""
function install_and_link(pkg::StringOrPkg; verbose::Bool=false) end

function install_and_link(name::AbstractString; verbose::Bool=false)
    # If we're already linked, don't sweat it.
    if linked(name)
        return
    end

    # Install dependency and link it
    brew(`install $name`; verbose=verbose)
    link(name; verbose=verbose)
end

function install_and_link(pkg::BrewPkg; verbose::Bool=false)
    return install_and_link(fullname(pkg); verbose=verbose)
end

"""
`postinstall(pkg::Union{AbstractString,BrewPkg}; verbose=false)`

Runs `brew postinstall` against package `pkg`, useful for debugging complicated
formulae when a bottle doesn't install right and you want to re-run postinstall.
"""
function postinstall(pkg::StringOrPkg; verbose::Bool=false) end

function postinstall(name::AbstractString; verbose::Bool=false)
    brew(`postinstall $name`, verbose=verbose)
end

function postinstall(pkg::BrewPkg; verbose::Bool=false)
    postinstall(fullname(pkg), verbose=verbose)
end

"""
`link(pkg::Union{AbstractString,BrewPkg}; verbose=false, force=true)`

Link package `name` into the global namespace, uses `--force` if `force == true`
"""
function link(name::StringOrPkg; verbose::Bool=false, force::Bool=true) end

function link(name::AbstractString; verbose::Bool=false, force::Bool=true)
    brew(`link $name`, no_stdout=true, verbose=verbose, force=force)
end

function link(pkg::BrewPkg; verbose::Bool=false, force::Bool=true)
    return link(fullname(pkg); force=force, verbose=verbose)
end

"""
`unlink(pkg::Union{AbstractString,BrewPkg}; verbose::Bool=false, quiet::Bool=true)`

Unlink package `pkg` from the global namespace, uses `--quiet` if `quiet == true`
"""
function unlink(name::StringOrPkg; verbose::Bool=false) end

function unlink(name::AbstractString; verbose::Bool=false)
    brew(`unlink $name`; verbose=verbose)
end

function unlink(pkg::BrewPkg; verbose::Bool=false)
    return unlink(fullname(pkg); verbose=verbose)
end

"""
`rm(pkg::Union{AbstractString,BrewPkg}; verbose::Bool=false, force::Bool=true)`

Remove package `pkg`, use `--force` if `force` == `true`
"""
function rm(pkg::StringOrPkg; verbose::Bool=false, force::Bool=false) end

function rm(pkg::AbstractString; verbose::Bool=false, force::Bool=true)
    brew(`rm $pkg`; verbose=verbose, force=force)
end

function rm(pkg::BrewPkg; verbose::Bool=false, force::Bool=true)
    return rm(fullname(pkg); verbose=verbose, force=force)
end

"""
`rm(pkgs::Vector{String or BrewPkg}; verbose::Bool=false, force::Bool=true)`

Remove packages `pkgs`, use `--force` if `force` == `true`
"""
function rm(pkgs::Vector{T}; verbose::Bool=false, force::Bool=false) where {T <: StringOrPkg}
    for pkg in pkgs
        try
            rm(pkg; verbose=verbose, force=force)
        catch
        end
    end
end

"""
`cleanup()`

Cleans up old installed versions of formulae, as well as purging all downloaded bottles
"""
function cleanup()
    brew(`cleanup -s`)
end

"""
`tap(tap_name::AbstractString; full::Bool=true, verbose::Bool=false)`

Runs `brew tap \$tap_name` if the tap does not already exist.  If `full` is `true`
adds the flag `--full` to clone a full tap instead of Homebrew's default shallow.

If `git` is not available, manually tap it with curl and tar.  This tap will be
unshallowed at the next `brew update` when `git` is available.
"""
function tap(tap_name::AbstractString; full::Bool=true, verbose::Bool=false)
    if !tap_exists(tap_name)
        # If we have no git available, then do it oldschool style
        if !git_installed()
            user = dirname(tap_name)
            repo = "homebrew-$(basename(tap_name))"
            tarball_url = "https://github.com/$user/$repo/tarball/master"
            tap_path = joinpath(brew_prefix,"Library","Taps", user, repo)

            if verbose
                @info("Manually tapping $tap_name...")
            end
            try
                mkpath(tap_path)
                run(pipeline(`curl -\# -L $tarball_url`, `tar xz -m --strip 1 -C $tap_path`))
            catch
                @warn("Could not download/extract $tarball_url into $(tap_path)!")
                rethrow()
            end
        else
            if full
                brew(`tap --full $tap_name`; verbose=verbose)
            else
                brew(`tap $tap_name`; verbose=verbose)
            end
        end
    end
end

const OSX_VERSION = [""]
function osx_version_string()
    global OSX_VERSION
    if isempty(OSX_VERSION[1])
        OSX_VERSION[1] = join(split(readchomp(`sw_vers -productVersion`), ".")[1:2],".")
    end
    return Dict(
        "10.4" => "tiger",
        "10.5" => "leopard",
        "10.6" => "snow_leopard",
        "10.7" => "lion",
        "10.8" => "mountain_lion",
        "10.9" => "mavericks",
        "10.10" => "yosemite",
        "10.11" => "el_capitan",
        "10.12" => "sierra",
        "10.13" => "high_sierra",
        "10.14" => "mojave",
        "10.15" => "catalina"
    )[OSX_VERSION[1]]
end

"""
`has_bottle(name::AbstractString)`

Checks if a given formula has a bottle at all
"""
function has_bottle(name::AbstractString)
    return haskey(json(name)["bottle"], "stable") &&
           haskey(json(name)["bottle"]["stable"]["files"], osx_version_string())
end

"""
`has_relocatable_bottle(name::AbstractString)`

Checks to see if a given formula has a bottle that can be installed anywhere
"""
function has_relocatable_bottle(name::AbstractString)
    if !has_bottle(name)
        return false
    end

    return json(name)["bottle"]["stable"]["cellar"] in [":any", ":any_skip_relocation"]
end

function versioninfo(;verbose=false)
    InteractiveUtils.versioninfo(stdout; verbose=verbose)
    run(Cmd(`git log -1 '--pretty=format:Homebrew git revision %h; last commit %cr'`; dir=joinpath(prefix())))
    run(Cmd(`git log -1 '--pretty=format:homebrew/core git revision %h; last commit %cr'`; dir=joinpath(prefix(), "Library", "Taps", "homebrew", "homebrew-core")))
    run(Cmd(`git log -1 '--pretty=format:staticfloat/juliadeps revision %h; last commit %cr'`; dir=joinpath(prefix(), "Library", "Taps", "staticfloat", "homebrew-juliadeps")))

    installed_pkgs = list()
    println("\n$(length(installed_pkgs)) total packages installed:")
    println(join([x.name for x in installed_pkgs], ", "))
end


text:
#__precompile__(false)
"""
Main module for `LikelihoodProfiler.jl`.

Four functions are exported from this module for public use:

- [`get_endpoint`](@ref). Computes lower or upper endpoints of confidence interval.
- [`get_interval`](@ref). Computes confidence interval.
- [`profile`](@ref). Generates the profile function based on `loss_func`
- [`update_profile_points!`](@ref). Updates confidence interval with likelihood profile points.

"""
module LikelihoodProfiler

using NLopt, ForwardDiff
using Calculus
using LinearAlgebra
using RecipesBase
import PlotUtils.adapted_grid
using ProgressMeter

# include
include("structures.jl")
include("get_endpoint.jl")
include("get_interval.jl")
include("cico_one_pass.jl")
include("method_lin_extrapol.jl")
include("method_quadr_extrapol.jl")
include("profile.jl")
include("plot_interval.jl")
include("show.jl")
include("get_optimal.jl")

# export
export get_right_endpoint,
    get_endpoint,
    scaling,
    unscaling,
    ProfilePoint,
    EndPoint,
    ParamIntervalInput,
    ParamInterval,
    get_interval,
    profile,
    update_profile_points!,
    get_optimal
end #module


text:
module GeometricAlgebra

import Base: ==, *, /, \, +, -, ^
import InteractiveUtils: subtypes
import Combinatorics: parity, powerset, permutations
import OffsetArrays: OffsetArray # only ever used by basis(::OffsetSignature)
using SparseArrays
using StaticArrays
using TaylorSeries: Taylor1

include("unit_blades.jl")

export AbstractMultivector, HomogeneousMultivector
export Blade, Multivector, MixedMultivector
include("types.jl")

export signature, dimension
export OffsetSignature
include("signatures.jl")

export ∧, ⋅, ⨼, ⨽, ∗
export ~, reversion, involute
export grade, grades, isscalar, scalar
export unit_pseudoscalar, dual
include("algebra.jl")

include("special.jl")

export basis, @basis, @basisall
include("convenience.jl")

include("show.jl")

export best_type

end # module


text:
using Mimi

include("multi-region-model.jl")
using .MyModel
model = construct_MyModel()

run(model)

# show results
getdataframe(model, :emissions, :E_Global)


text:
module PerfMachine

include("sbatch.jl")
include("runner.jl")
include("analyze.jl")

end


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

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

module QueensProblems

export ModuleQueensProblems
export L319284, Queens

"""

For some background see: [Backtracking with profiles](https://wp.me/paipV7-E).

* Queens, L319284
"""
const ModuleQueensProblems = ""

function solve!(profile, level, size, start, cols, diag4, diag1)

    if level > 0
        for i ∈ start:size
            save = cols & (1 << i) + diag1 & (1 << (i + level)) +
                   diag4 & (1 << (32 + i - level))

            if save == 0

                cols = xor(cols, 1 << i)
                diag1 = xor(diag1, 1 << (i + level))
                diag4 = xor(diag4, 1 << (32 + i - level))

                solve!(profile, level - 1, size, 0, cols, diag4, diag1)

                cols = xor(cols, 1 << i)
                diag1 = xor(diag1, 1 << (i + level))
                diag4 = xor(diag4, 1 << (32 + i - level))

                profile[level+1] += 1
            end
        end
    else
        for i ∈ 0:size
            save = cols & (1 << i) + diag1 & (1 << i) + diag4 & (1 << (32 + i))
            save == 0 && (profile[1] += 1)
        end
    end
end

function search(n::Int)
    profile = zeros(Int, n + 1)
    cols = diag4 = diag1 = Int(0)
    solve!(profile, n - 1, n - 1, 0, cols, diag4, diag1)
    return profile
end

"""

Returns the profile of the backtrack tree for the n queens problem (see `A319284`).
"""
function Queens(n::Int)
    n == 0 && return [1]
    profile = search(n)
    profile[n+1] = 1  # add the root
    [profile[n-i+1] for i = 0:n]
end

"""

Returns the profile of the backtrack tree for the n queens problem.
"""
L319284(n) = Queens(n)

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

using Test

function test()
    @testset "QueensProblem" begin
        levels = [1, 10, 72, 364, 1400, 3916, 7552, 9632, 7828, 4040, 724]
        @test all(L319284(10) .== levels)
    end
end

function demo()

    up_to = 10
    for n ∈ 0:up_to
        print("elapsed: ")
        @time profile = Queens(n)
        println("size:      ", n)
        println("profile:   ", profile)
        println("nodes:     ", sum(profile))
        println("solutions: ", profile[n+1])
        println()
    end
end

"""

L319284(15) :: 18.731622 seconds (4 allocations: 480 bytes)
"""
function perf()
    @time L319284(15)
end

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

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# write some simple stuff to the map
function gen_spin0!(m, powerlaw_exponent=2)
    for i in axes(m, 1)
        for j in axes(m, 2)
            m[i,j] = ( ((j-1) * size(m, 1)) + i )^powerlaw_exponent
        end
    end
end

@testset "spin-0 SHTs" begin
    shape, wcs = fullsky_geometry(10.0 * Pixell.degree)
    m = Enmap(zeros(shape), wcs)
    gen_spin0!(m)

    alms = map2alm(m; lmax=18)
    data = readdlm("data/simple_analytic_sht.txt")
    ref_alms = data[:,1] + 1im .* data[:,2]
    @test (alms.alm ≈ ref_alms)

    @test length(map2alm(m).alm) == 190

    alms = map2alm(m[6:end-2, 5:end-3]; lmax=18)
    data = readdlm("data/simple_analytic_sht_sliced.txt")
    ref_alms = data[:,1] + 1im .* data[:,2]
    @test (alms.alm ≈ ref_alms)

    box = [10   -10;           # RA
           -5     5] * Pixell.degree  # DEC
    shape, wcs = geometry(CarClenshawCurtis, box, 1.0 * Pixell.degree)
    m = Enmap(zeros(shape), wcs)
    gen_spin0!(m, 2.5)
    alms = map2alm(m; lmax=100)
    data = readdlm("data/simple_box_analytic_sht.txt")
    ref_alms = data[:,1] + 1im .* data[:,2]
    @test (alms.alm ≈ ref_alms)
end

##
function gen_spin2!(m)
    for i in axes(m, 1)
        for j in axes(m, 2)
            m[i,j,1] = i^2 * j
            m[i,j,2] = i * j
        end
    end
end

@testset "spin-2 SHTs" begin
    shape, wcs = fullsky_geometry(10.0 * Pixell.degree; dims=(2,))
    m_QU = Enmap(zeros(shape), wcs)
    gen_spin2!(m_QU)
    alm_e, alm_b = map2alm(m_QU; lmax=3 * 36)
    data = readdlm("data/simple_pol_analytic_sht.txt")
    ref_alm_e = data[:,1] + 1im .* data[:,2]
    ref_alm_b = data[:,3] + 1im .* data[:,4]
    @test (ref_alm_e ≈ alm_e.alm)
    @test (ref_alm_b ≈ alm_b.alm)

    alm_e, alm_b = map2alm((m_QU[:,:,1], m_QU[:,:,2]); lmax=3 * 36)
    @test (ref_alm_e ≈ alm_e.alm)
    @test (ref_alm_b ≈ alm_b.alm)

    shape, wcs = fullsky_geometry(10.0 * Pixell.degree; dims=(3,))
    m_IQU = Enmap(zeros(shape), wcs)
    gen_spin0!(@view(m_IQU[:,:,1]))
    gen_spin2!(@view(m_IQU[:,:,2:3]))
    alm_t, alm_e, alm_b = map2alm(m_IQU; lmax=3 * 36)

    data = readdlm("data/simple_analytic_sht_fullalm.txt")
    ref_alm_t = data[:,1] + 1im .* data[:,2]
    @test (ref_alm_t ≈ alm_t.alm)
    @test (ref_alm_e ≈ alm_e.alm)
    @test (ref_alm_b ≈ alm_b.alm)

    alm_t, alm_e, alm_b = map2alm((m_IQU[:,:,1], m_IQU[:,:,2], m_IQU[:,:,3]); lmax=3 * 36)
end
##

@testset "alm2map" begin
    alm_r = readdlm("data/TEB_alms_real.dat")
    alm_i = readdlm("data/TEB_alms_imag.dat")
    lmax = 540
    alms = [Alm(lmax,lmax), Alm(lmax,lmax), Alm(lmax,lmax)]
    for i in 1:3
        alms[i].alm .+= alm_r[i,:] .+ 1im .* alm_i[i,:]
    end

    test_IQU = read_map("data/IQU.fits")
    I, Q, U = alm2map(alms, size(test_IQU), test_IQU.wcs)

    pyI = read_map("data/pyI.fits")[:,:,1]
    pyQU = read_map("data/pyQU.fits")

    @test maximum(abs.(I .- pyI)) < 1e-9
    @test maximum(abs.(Q .- pyQU[:,:,1])) < 1e-9
    @test maximum(abs.(U .- pyQU[:,:,2])) < 1e-9
end

@testset "map2alm CMB" begin
    test_IQU = read_map("data/IQU.fits")
    aT, aE, aB = map2alm(test_IQU; lmax=540)
    data = readdlm("data/test_cls_IQU.txt")
    @test maximum(abs.(alm2cl(aT) .- data[:,1])) < 1e-10
    @test maximum(abs.(alm2cl(aT, aE) .- data[:,2])) < 1e-10
    @test maximum(abs.(alm2cl(aE) .- data[:,3])) < 1e-10
    @test maximum(abs.(alm2cl(aB) .- data[:,4])) < 1e-10
end


text:
using Test

## Ratios
@testset "ef implementation for a single parameter" begin
    using GEEBRA
    using Random

    ## Ratio data
    struct ratio_data
        y::Vector
        x::Vector
    end

    ## Ratio nobs
    function ratio_nobs(data::ratio_data)
        nx = length(data.x)
        ny = length(data.y)
        if (nx != ny)
            error("length of x is not equal to the length of y")
        end
        nx
    end

    ## Ratio contribution to estimating function
    function ratio_ef(theta::Vector,
                      data::ratio_data,
                      i::Int64)
        p = length(theta)
        (data.y[i] .- theta * data.x[i])[1:p]
    end

    ## Set the ratio template
    ratio_template = estimating_function_template(ratio_nobs, ratio_ef)
    @inferred estimating_function_template(ratio_nobs, ratio_ef)

    ## Generate some data
    Random.seed!(123);
    my_data = ratio_data(randn(10), rand(10));

    ## Get M-estimator for the ratio
    result_m = fit(ratio_template, my_data, [0.1], estimation_method = "M")
    @inferred fit(ratio_template, my_data, [0.1], estimation_method = "M")
    ## Get reduced-bias estimator for the ratio
    result_br = fit(ratio_template, my_data, [0.1], estimation_method = "RBM")
    @inferred fit(ratio_template, my_data, [0.1], estimation_method = "RBM")
    ## Gere reduced-bias estimator for the ration using explicit RBM-estimation
    result_br1 = fit(ratio_template, my_data, [0.1], estimation_method = "RBM", br_method = "explicit_trace")

    ## Quantities for estimators
    sx = sum(my_data.x)
    sxx = sum(my_data.x .* my_data.x)
    sy = sum(my_data.y)
    sxy = sum(my_data.x .* my_data.y)

    @test isapprox(sy/sx, coef(result_m)[1])
    @test isapprox((sy + sxy/sx)/(sx + sxx/sx), coef(result_br)[1])
    @test isapprox(sy/sx * (1 - sxx / sx^2) + sxy/sx^2, coef(result_br1)[1])
    @test_throws ErrorException  result_br1 = fit(ratio_template, my_data, [0.1], estimation_method = "RBM", br_method = "magic_br_method")
end

## Instrumental variables
@testset "ef implementation for multiple parameters" begin
    using GEEBRA
    using Random
    using Distributions
    using NLsolve

    ## IV data
    struct iv_data
        y::Vector
        t::Vector
        w::Vector
    end

    ## IV nobs
    function iv_nobs(data::iv_data)
        nw = length(data.w)
        ny = length(data.y)
        nt = length(data.t)
        if (nw != ny)
            error("length of w is not equal to the length of y")
        elseif (nw != nt)
            error("length of w is not equal to the length of t")
        elseif (nw != nt)
            error("length of w is not equal to the length of t")
        end
        nw
    end

    ## Contributions to the estimating functions
    function iv_ef(theta::Vector,
                   data::iv_data,
                   i::Int64)
        [theta[1] - data.t[i], (data.y[i] - theta[2] * data.w[i]) * (theta[1] - data.t[i])]
    end

    ## Simulating IV data
    function simulate_iv(nobs::Int,
                         theta::Vector)
        alpha = theta[1]
        beta = theta[2]
        gamma = theta[3]
        delta = theta[4]
        mux = theta[5]
        sigmax = theta[6]
        sigmae = theta[7]
        sigmau = theta[8]
        sigmat = theta[9]
        e1 = rand(Normal(0, sigmae), nobs)
        e2 = rand(Normal(0, sigmau), nobs)
        e3 = rand(Normal(0, sigmat), nobs)
        x = rand(Normal(mux, sigmax), nobs)
        w = x + e2
        y = alpha .+ beta * x + e1
        t = gamma .+ delta * x + e3
        iv_data(y, t, w)
    end

    ## Set up IV GEEBRA template
    iv_template = estimating_function_template(iv_nobs, iv_ef)
    @inferred estimating_function_template(iv_nobs, iv_ef)

    ## Simulate data
    true_theta = [2.0, 2.0, 1.0, 3.0, 0.0, 1.0, 2.0, 1.0, 1.0]
    true_parameter = true_theta[[3, 2]]

    Random.seed!(123)
    my_data = simulate_iv(100, true_theta)

    o1_ml = fit(iv_template, my_data, true_parameter, estimation_method = "M")
    @inferred fit(iv_template, my_data, true_parameter, estimation_method = "M")
    o1_br = fit(iv_template, my_data, true_parameter, estimation_method = "RBM")
    @inferred fit(iv_template, my_data, true_parameter, estimation_method = "RBM")
    o1_br1 = fit(iv_template, my_data, true_parameter, estimation_method = "RBM", br_method = "explicit_trace")
    @inferred fit(iv_template, my_data, true_parameter, estimation_method = "RBM", br_method = "explicit_trace")

    ef_br = get_estimating_function(my_data, iv_template, true)
    @inferred get_estimating_function(my_data, iv_template, true)
    ef_ml = get_estimating_function(my_data, iv_template, false)
    @inferred get_estimating_function(my_data, iv_template, false)
    o2_ml = nlsolve(ef_ml, [0.1, 0.2])
    o2_br = nlsolve(ef_br, [0.1, 0.2])

    qs = GEEBRA.ef_quantities(coef(o1_ml), my_data, iv_template, true)
    @test isapprox(coef(o1_ml) + qs[2] * qs[1], coef(o1_br1))

    @test isapprox(coef(o1_ml), o2_ml.zero)
    @test isapprox(coef(o1_br), o2_br.zero)

    ## Estimating function at the estimates
    @test isapprox(estimating_function(o2_br.zero, my_data, iv_template, true),
                   zeros(Float64, 2, 1),
                   atol = 1e-10)
end

@testset "obj implementation for multiple parameters" begin
    using GEEBRA
    using Random
    using Distributions
    using Optim
    using LinearAlgebra

    ## Logistic regression data
    struct logistic_data
        y::Vector
        x::Array{Float64}
        m::Vector
    end

    ## Logistic regression nobs
    function logistic_nobs(data::logistic_data)
        nx = size(data.x)[1]
        ny = length(data.y)
        nm = length(data.m)
        if (nx != ny)
            error("number of rows in of x is not equal to the length of y")
        elseif (nx != nm)
            error("number of rows in of x is not equal to the length of m")
        elseif (ny != nm)
            error("length of y is not equal to the length of m")
        end
        nx
    end

    function logistic_loglik(theta::Vector,
                             data::logistic_data,
                             i::Int64)
        eta = sum(data.x[i, :] .* theta)
        mu = exp.(eta)./(1 .+ exp.(eta))
        ll = data.y[i] .* log.(mu) + (data.m[i] - data.y[i]) .* log.(1 .- mu)
        ll
    end

    Random.seed!(123);
    n = 100;
    m = 1;
    x = Array{Float64}(undef, n, 2);
    x[:, 1] .= 1.0;
    x[:, 2] .= rand(n);
    true_betas = [0.5, -1];
    y = rand.(Binomial.(m, cdf.(Logistic(), x * true_betas)));
    my_data = logistic_data(y, x, fill(m, n));

    logistic_template = objective_function_template(logistic_nobs, logistic_loglik)
    @inferred objective_function_template(logistic_nobs, logistic_loglik)

    o1_ml = optimize(b -> -objective_function(b, my_data, logistic_template, false),
                     true_betas, LBFGS())
    o1_br = optimize(b -> -objective_function(b, my_data, logistic_template, true),
                     true_betas, LBFGS())

    o2_ml = fit(logistic_template, my_data, true_betas, estimation_method = "M")
    o2_br = fit(logistic_template, my_data, true_betas, estimation_method = "RBM")

    o3_ml = optimize(b -> -objective_function(b, my_data, logistic_template, false),
                     true_betas, Optim.Options(iterations = 2))
    o3_br = optimize(b -> -objective_function(b, my_data, logistic_template, true),
                     true_betas, Optim.Options(iterations = 2))

    o4_ml = fit(logistic_template, my_data, true_betas, estimation_method = "M",
                               optim_method = NelderMead(),
                               optim_options = Optim.Options(iterations = 2))
    o4_br = fit(logistic_template, my_data, true_betas, estimation_method = "RBM",
                               optim_method = NelderMead(),
                               optim_options = Optim.Options(iterations = 2))

    @test isapprox(Optim.minimizer(o1_ml), Optim.minimizer(o2_ml.results))
    @test isapprox(Optim.minimizer(o1_br), Optim.minimizer(o2_br.results))
    @test isapprox(Optim.minimizer(o3_ml), Optim.minimizer(o4_ml.results))
    @test isapprox(Optim.minimizer(o3_br), Optim.minimizer(o4_br.results))

    @test isapprox(sqrt.(diag(vcov(o2_br))), stderror(o2_br))

end

@testset "agreement between obj and ef implementations" begin
    using GEEBRA
    using Random
    using Distributions
    using Optim
    using NLsolve

    ## Logisti regression data
    struct logistic_data
        y::Vector
        x::Array{Float64}
        m::Vector
    end

    ## Logistic regression nobs
    function logistic_nobs(data::logistic_data)
        nx = size(data.x)[1]
        ny = length(data.y)
        nm = length(data.m)
        if (nx != ny)
            error("number of rows in of x is not equal to the length of y")
        elseif (nx != nm)
            error("number of rows in of x is not equal to the length of m")
        elseif (ny != nm)
            error("length of y is not equal to the length of m")
        end
        nx
    end

    function logistic_loglik(theta::Vector,
                             data::logistic_data,
                             i::Int64)
        eta = sum(data.x[i, :] .* theta)
        mu = exp.(eta)./(1 .+ exp.(eta))
        ll = data.y[i] .* log.(mu) + (data.m[i] - data.y[i]) .* log.(1 .- mu)
        ll
    end

    function logistic_ef(theta::Vector,
                         data::logistic_data,
                         i::Int64)
        eta = sum(data.x[i, :] .* theta)
        mu = exp.(eta)./(1 .+ exp.(eta))
        data.x[i, :] * (data.y[i] - data.m[i] * mu)
    end

    Random.seed!(123);
    n = 100;
    m = 1;
    p = 5;
    x = Array{Float64}(undef, n, p);
    x[:, 1] .= 1.0;
    for j in 2:p
        x[:, j] .= rand(n);
    end
    true_betas = randn(p) * sqrt(p);
    y = rand.(Binomial.(m, cdf.(Logistic(), x * true_betas)));
    my_data = logistic_data(y, x, fill(m, n));

    logistic_obj_template = objective_function_template(logistic_nobs, logistic_loglik)
    logistic_ef_template = estimating_function_template(logistic_nobs, logistic_ef)

    o1_ml = fit(logistic_obj_template, my_data, true_betas, estimation_method = "M")
    e1_ml = fit(logistic_ef_template, my_data, true_betas, estimation_method = "M")
    o1_br = fit(logistic_obj_template, my_data, true_betas, estimation_method = "RBM")
    e1_br = fit(logistic_ef_template, my_data, true_betas, estimation_method = "RBM")
    o1_br1 = fit(logistic_obj_template, my_data, true_betas, estimation_method = "RBM", br_method = "explicit_trace")
    o1_br2 = fit(logistic_obj_template, my_data, coef(o1_ml), estimation_method = "RBM", br_method = "explicit_trace")
    e1_br1 = fit(logistic_ef_template, my_data, true_betas, estimation_method = "RBM", br_method = "explicit_trace")
    e1_br2 = fit(logistic_ef_template, my_data, coef(o1_ml), estimation_method = "RBM", br_method = "explicit_trace")

    @test isapprox(coef(o1_ml), coef(e1_ml), atol = 1e-05)
    @test isapprox(coef(o1_br), coef(e1_br), atol = 1e-05)
    @test isapprox(coef(o1_br1), coef(e1_br1), atol = 1e-05)
    @test isapprox(coef(o1_br2), coef(e1_br2), atol = 1e-05)
    @test isapprox(coef(o1_br1), coef(e1_br2), atol = 1e-05)

    @test isapprox(aic(o1_ml),
                   -2 * (objective_function(coef(o1_ml), my_data, logistic_obj_template, false) - p))

    @test isapprox(aic(o1_br),
                   -2 * (objective_function(coef(o1_br), my_data, logistic_obj_template) - p))

    quants_ml = GEEBRA.obj_quantities(coef(o1_ml), my_data, logistic_obj_template, true)
    quants_br = GEEBRA.obj_quantities(coef(o1_br), my_data, logistic_obj_template, true)

    @test isapprox(tic(o1_ml),
                   -2 * (objective_function(coef(o1_ml), my_data, logistic_obj_template) + 2 * quants_ml[1]))

    @test isapprox(tic(o1_br),
                   -2 * (objective_function(coef(o1_br), my_data, logistic_obj_template) + 2 * quants_br[1]))

    @test isapprox(vcov(o1_ml), vcov(e1_ml))
    @test isapprox(vcov(o1_br), vcov(e1_br))
    @test isapprox(vcov(o1_br1), vcov(e1_br1))
    @test isapprox(vcov(o1_br2), vcov(e1_br2))

    @test isapprox(coeftable(o1_ml).cols, coeftable(e1_ml).cols)
    @test isapprox(coeftable(o1_br).cols, coeftable(e1_br).cols)
    @test isapprox(coeftable(o1_br1).cols, coeftable(e1_br1).cols)
    @test isapprox(coeftable(o1_br2).cols, coeftable(e1_br2).cols)

end

# using Revise
# using Pkg
# Pkg.activate("/Users/yiannis/Repositories/GEEBRA.jl")


text:
function gen_player(model, number)

    if model == "uniform"
        # generate n uniformly spread over a given range
    elseif model == "random"
        # generate randomly according to given distribution
    end

end