text:
#= This dictionary translate numbers into the corresponding information:
    values [1,16] are state-info values [17,35] are event-info
=#
InfoDict = Dict(
    "1"=>"delay_forage_left",
    "2"=>"start_forage_left",
    "3"=>"in_left",
    "4"=>"out_left",
    "5"=>"reward_left_available",
    "6"=>"reward_left",
    "7"=>"reward_consumption_left",
    "8"=>"delay_forage_right",
    "9"=>"start_forage_right",
    "10"=>"in_right",
    "11"=>"out_right",
    "12"=>"reward_right_available",
    "13"=>"reward_right",
    "14"=>"reward_consumption_right",
    "15"=>"travel",
    "16"=>"travel_available",
    "17"=>"poke_4",
    "18"=>"poke_4_out",
    "19"=>"poke_2",
    "20"=>"poke_2_out",
    "21"=>"poke_3",
    "22"=>"poke_3_out",
    "23"=>"poke_6",
    "24" =>"poke_6_out",
    "25"=>"poke_9",
    "26"=>"poke_9_out",
    "27"=>"session_timer",
    "28"=>"forage_timer",
    "29"=>"travel_tone_increment",
    "30"=>"travel_out_of_poke",
    "31"=>"travel_resumed",
    "32"=>"travel_complete",
    "33"=>"task_disengagment",
    "34"=>"post_travel_delay_timer",
    "35"=>"reward_consumption_timer",
    "36"=>"rsync"
)

PortDict = Dict(
    "poke_2"=>"PokeLeft",
    "poke_4"=>"RewLeft",
    "poke_3"=>"PokeRight",
    "poke_6"=>"RewRight",
    "poke_9"=>"TravPoke"
)

"""
    process_raw_session(session_path::String; observe = false)

Given a string indicating the path to a rawdata text file produce a table
where each row is a poke in a specific port and columns containing timing
and state info about that poke. If observe is true opens intermediate and
final tables in html for review
"""
function process_raw_session(session_path::String; observe = false)
    !ispath(session_path) && error("\"$(session_path)\" is not a viable path")
    lines = readlines(session_path)
    eachlines = eachline(session_path)
    prov = table_raw_data(lines, eachlines)
    observe && open_html_table(prov)
    pokes = combine(groupby(prov, :Poke)) do dd
        adjustevents(dd)
    end
    sort!(pokes,:PokeIn)
    # pokes[!,:Duration] = pokes.PokeOut - pokes.PokeIn
    RichnessDict_keys = sort(collect(keys(countmap(pokes.IFT))))
    RichnessDict = Dict(x => y for (x,y) in zip(RichnessDict_keys,["rich", "medium", "poor", missing]))
    pokes[!,:State] = find_task_state(pokes)
    pokes[:,:Port] = [get(PortDict,x, x) for x in pokes.Port]
    idx = findall(ismatch.(r"^reward_consumption_",pokes.Port))
    for i in idx
        port = ismatch(r"left",pokes[i,:Port]) ? "RewLeft" : "RewRight"
        rewardedpoke = findprev(ismatch.(Regex(port),pokes.Port),i)
        pokes[rewardedpoke,:RewardConsumption] = pokes[i,:PokeIn]
    end
    transform!(pokes,[:Port,:TravelComplete] => ((p,t)->activeside(p,t)) => :ActivePort)
    # pokes[!,:Incorrect] = @. !ismatch(Regex(pokes.ActivePort),pokes.Port) && ismatch(r"^Poke", pokes.Port)
    incorrectpokes!(pokes)
    transform!(pokes, :State => count_patches => :Patch)
    transform!(groupby(pokes,:Patch), :T => determine_travel => :Travel)
    transform!(groupby(pokes,:Patch), :IFT => (x -> determine_richness(x,RichnessDict)) => :Richness)

    #occasionally pokeout have missing values in that case it search for the poke_out value accoriding to the poke_out number
    checkPokeOut  = ismissing.(pokes.PokeOut) .&& ismatch.(r"^Poke",pokes.Port)
    if any(checkPokeOut)
        for i in findall(checkPokeOut)
            pnum = pokes[i,:Poke]
            if ismissing(pokes[i,:Poke])
                println("poke deleted skipping poke out time correction")
                continue
            end
            pos = findfirst(prov.PokeOut_count .== pnum)
            pokes[i,:PokeOut] = prov[pos,:Time]
        end
    end

    pokes[!,:Duration] = pokes.PokeOut - pokes.PokeIn

    session_info = parseGlobalInfo(lines)
    for (x,y) in session_info
        pokes[!, x] .= y
    end
    pokes =  pokes[:,[:Richness, :Travel,:State,:Poke,:ActivePort,
        :Incorrect,:Port, :PokeIn, :PokeOut, :Duration,
        :Patch, :P,#:Bout,:AlternativePatch,
        :TravelOnset,:TravelComplete,
        :RewardAvailable, :RewardDelivery, :RewardConsumption,
        :Rsync_count, :Rsync_time, :StateIn, :StateOut, :T, :AFT, :IFT,
        :SubjectID, :Taskname, :Experimentname, :Startdate]]
    sort!(pokes,:PokeIn)
    observe && open_html_table(pokes)
    return pokes
end

"""
    parseGlobalInfo(lines::AbstractVector{String})

When feeded a raw data text file in a vector of string format, uses the first
five elements to extract session's global information. The output is a tuple
therefore the returning variable must respect the correct order:
Experiment, Task, TaskFileHash, SetUpID, SubjectID

"""
function parseGlobalInfo(lines::AbstractVector{String})
    # session_info_string = filter(l-> !isempty(l) && l[1] == 'I',lines)
    # tuple([replace(split(x,":")[2]," "=>"") for x in session_info_string]...)
    session_info_string = filter(l-> !isempty(l) && l[1] == 'I',lines)
    tuple([replace(split(x,":")[2]," "=>"") for x in session_info_string]...)
    session_info_string = replace.(session_info_string,"I "=> "")
    session_info_string = replace.(session_info_string," : "=> "|")
    session_info_string = replace.(session_info_string," "=> "")
    session_info = Dict()#(;)
    for x in session_info_string
        n,v = split(x,"|")
        session_info = merge(session_info,Dict(Symbol(n)=>v))
    end
    return session_info
end

"""
    adjustevents(df::AbstractDataFrame)
Using a GroupedDataFrame per poke, originating from
`table_raw_data(session_path::String)`extracts the info about identity
or timing of variable of interest and collect them in separate columns
"""
function adjustevents(df::AbstractDataFrame)
    if !ismissing(df[1,:Poke])
        port = findport(df.Event)
        DataFrame(
            Port = port,#findport(df.Event),
            PokeIn = event_time(Regex(port), df),#event_time(r"^poke_\d$", df),
            PokeOut = event_time(Regex(port*"_out"), df),#event_time(r"^poke_\d_out$",df),
            StateIn = event_time(r"(^in_((right)|(left))$)|(^travel$)", df),
            StateOut = event_time(r"^out_((right)|(left))$", df),
            RewardAvailable = event_time(r"^reward_((right)|(left))_available$",df),
            RewardDelivery = event_time(r"^reward_((right)|(left))$",df),
            RewardConsumption = event_time(r"^reward_consumption_((right)|(left))$",df),
            TravelOnset = event_time(r"^travel$",df),
            TravelComplete = event_time(r"^travel_complete$",df),
            Rsync_count = findrsynccount(df),
            Rsync_time = event_time(r"^rsync$",df),
            P = df[end,:P], #always take the latest Patch value
            T = df[end,:T], #always take the latest T value
            AFT = df[1,:AFT] == df[end,:AFT] ? df[1,:AFT] : missing, #if AFT changes report missing val
            IFT = df[1,:IFT] == df[end,:IFT] ? df[1,:IFT] : missing #if IFT changes report missing val
            )
    else
        pre = DataFrame(
            Port = df.Event,
            PokeIn = df.Time,
            PokeOut = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            StateIn = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            StateOut = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            RewardAvailable = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            RewardDelivery = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            RewardConsumption = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            TravelOnset = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            TravelComplete = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            Rsync_count = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            Rsync_time = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            P = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            T = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            AFT = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            IFT = Vector{Union{Float64,Missing}}(missing,nrow(df)),
            )
        pre[pre.Port .== "rsync", :Rsync_time] .= pre[pre.Port .== "rsync", :PokeIn]
        pre[pre.Port .== "rsync", :Rsync_count] .= df[df.Event .== "rsync", :Rsync_count]
        return pre
    end
end

"""
    findport(vec::AbstractVector{String})

Using Event column resulting from a by-poke subgroup (from
`table_raw_data(session_path::String)`) filters all element not matching
string containing exactly "poke_" followed by one digit (/d) and NOT concluding
there exactly (in case there is a poke_out without poke_in at the beginning of
the session). Then use the matching result to identify to which port the
subgroup information belong
"""
function findport(vec::AbstractVector{String})
    poke_id = match.(r"^poke_\d",vec)
    all(isnothing.(poke_id)) ? missing : filter(!isnothing,poke_id)[1].match
end

"""
    event_time(expression::Regex,df::AbstractDataFrame)
Use a regular expression to find a matching string in the Event column
in the by-poke subgroup (from `table_raw_data(session_path::String)`) and
return the corresponding value from the Time column
"""
function event_time(expression::Regex,df::AbstractDataFrame)
    pos = findstring(expression,df.Event)
    isnothing(pos) ? missing : df[pos, :Time]
end

"""
    findstring(expression::Regex,vec::AbstractVector{String})
Use a regular expression to find the index of the first element in a vector of
Strings that matches the expression
"""
function findstring(expression::Regex,vec::AbstractVector)
    # findfirst(.!isnothing.(match.(expression,vec)))
    findfirst(ismatch.(expression,vec))
end

ismatch(r::Regex,s) = !isnothing(match(r, s))

function findrsynccount(df::AbstractDataFrame)
    pos = findfirst(.!isnothing.(match.(r"^rsync$",df.Event)))
    isnothing(pos) ? missing : df[pos,:Rsync_count]
end

function alternative_count_patches(vec::AbstractVector{AbstractString})
    res = Vector{Int64}(undef,0)
    count = 1
    foraging = true
    for x in vec
        if ismatch(r"^poke_[0-8]$",x) && !foraging
            count += 1
            foraging = true
        elseif ismatch(r"^poke_9$",x) && foraging
            foraging = false
        end
        push!(res,count)
    end
    return res
end

function find_task_state(df::AbstractDataFrame)
    state = "Forage"
    travel = false
    res = []
    # the order of updating travel and pushing is important to get the right state in each row
    for x in eachrow(df)
        if ismatch(r"^poke_9", x.Port) && (!travel) && !ismissing(x.TravelOnset)
            travel = true
            push!(res, travel ? "Travel" : "Forage")
        # elseif (travel) && (!ismissing(x.TravelComplete))
        elseif (travel) && ((ismatch(r"^start_forage_",x.Port)) || (!ismissing(x.TravelComplete)))
            push!(res, travel ? "Travel" : "Forage")
            travel = false
        else
            push!(res, travel ? "Travel" : "Forage")
        end
    end
    return res
end

function activeside(ev_vec::AbstractVector{AbstractString},trav_vec)
    first_idx = findstring(r"^start_forage",ev_vec)
    side = activeside(ev_vec[first_idx])
    ActivePort = repeat([side],length(ev_vec))
    # change_idx = findall(.!ismissing.(trav_vec))
    change_idx = findall(ismatch.(r"^start_forage_",ev_vec))
    for i in change_idx
        ActivePort[i:end] .= activeside(ev_vec[i])
        # side == "Left" ? (side = "Right") : (side = "Left")
        # ActivePort[i:end] .= side #activeside(ev_vec[i])
    end
    return ActivePort
end

activeside(x::String) = ismatch(r"start_forage_left",x) ? "Left" : "Right"

function incorrectpokes!(df::AbstractDataFrame)
    wrongside = @. !ismatch(Regex(df.ActivePort),df.Port) && ismatch(r"^Poke", df.Port)
    travelling = @. df.State == "Travel" && ismatch(r"(^Poke)|(^Rew)", df.Port)
    df[!,:Incorrect] = @. wrongside || travelling
end

function count_patches(state_vec::AbstractVector)
    c = 1
    res = [1]
    for i in 2:length(state_vec)
        (state_vec[i-1] == "Travel") && (state_vec[i] == "Forage") && (c += 1)
        push!(res,c)
    end
    return res
end

function determine_travel(vec)
    clean = filter(x -> !ismissing(x), vec)
    isempty(clean) ? missing : minimum(clean) < 2500 ? "Short" : "Long"
end

function determine_richness(vec, Rdict)
    clean = filter(x -> !ismissing(x), vec)
    isempty(clean) ? missing : get(Rdict,minimum(clean),"nope")
end


text:
######## StanSample example  ###########

using StanSample, Test

bernoullimodel = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
"

bernoullidata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
inittheta = Dict("theta" => 0.60)

sm = SampleModel("init_dict", bernoullimodel,
  seed=StanSample.RandomSeed(seed=-1));

rc = stan_sample(sm, data=bernoullidata, init=inittheta)

if success(rc)
  samples = read_samples(sm)

  df = read_summary(sm)
  @test df[df.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.2

end

ext: jl


text:
######
# methods that converge quickly even when α is non-simple
# These usually take 1 additional function call compared to a similar ordered
# method
######

"""
    Roots.Order1B()
    Roots.King()

A superlinear (order `1.6...`) modification of the secant method for multiple roots.
Presented in A SECANT METHOD FOR MULTIPLE ROOTS, by RICHARD F. KING, BIT 17 (1977), 321-328

The basic idea is similar to Schroder's method: apply the secant method
to  `f/f'`. However, this uses `f' ~ fp = (fx - f(x-fx))/fx` (a Steffensen step). In
this implementation, `Order1B`, when `fx` is too big, a single secant step of `f`
is used.

The *asymptotic* error, `eᵢ = xᵢ - α`, is given by
`eᵢ₊₂ = 1/2⋅G''/G'⋅ eᵢ⋅eᵢ₊₁ + (1/6⋅G'''/G' - (1/2⋅G''/G'))^2⋅eᵢ⋅eᵢ₊₁⋅(eᵢ+eᵢ₊₁)`.

"""
struct King <: AbstractSecant end
struct Order1B <: AbstractSecant end

struct KingState{T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    fxn1::S
    fxn0::S
    G0::S
end

function init_state(M::Union{King, Order1B}, F, x₀, x₁::T, fx₀, fx₁::S) where {T,S}
    fₛ₀ = F(x₀ - fx₀*oneunit(T)/oneunit(S))
    G₀ = -fx₀^2 / (fₛ₀ - fx₀)
    KingState(x₁, x₀, fx₁, fx₀, G₀)
end
initial_fncalls(::Union{King, Order1B}) = 3

function update_state(method::Order1B, F, o::KingState{T,S}, options, l=NullTracks())  where {T, S}
    if do_guarded_step(Order1B(), o)
        state, flag = update_state(Order1(), F, o, options, l)

        x0,fx0 = state.xn0, state.fxn0 # clunky! Need to update G₀ after Order1() step
        fₛ = F(x0 - fx0*oneunit(T)/oneunit(S))
        incfn(l)
        G₀ = -fx0^2 / (fₛ - fx0)
        @set! state.G0 = G₀

        return (state, flag)
    else
        update_state(King(), F, o, options, l)
    end
end

function update_state(method::King, F,
                      o::KingState{T,S}, options, l=NullTracks())  where {T, S}

    x0, x1 = o.xn0, o.xn1
    fx0, fx1 = o.fxn0, o.fxn1
    G₀ = o.G0
    fₛ₁ = F(x1 - fx1*oneunit(T)/oneunit(S))
    incfn(l)
    G₁ = -fx1^2 / (fₛ₁ - fx1)

    m = (x1-x0) / (G₁-G₀) # approximate value of `m`, the multiplicity

    if abs(m) <= 1e-2 * oneunit(m)
        log_message(l, "Estimate for multiplicity has issues. ")
        return (o, true)
    end

    Δ = G₁ * (x1 - x0) / (G₁ - G₀)

    if isissue(Δ)
        log_message(l, "Increment `Δx` has issues. ")
        return o, true
    end

    x0, x1 = x1, x1 - Δ
    fx0, fx1 = fx1, F(x1)
    incfn(l)

    @set! o.xn0 = x0
    @set! o.xn1 = x1
    @set! o.fxn0 = fx0
    @set! o.fxn1 = fx1
    @set! o.G0 = G₁

    return o, false

end

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

### Order2B() Esser method
"""
    Roots.Order2B()
    Roots.Esser()

Esser's method. This is a quadratically convergent method that, like
Schroder's method, does not depend on the multiplicity of the
zero. Schroder's method has update step `x - r2/(r2-r1) * r1`, where `ri =
fⁱ⁻¹/fⁱ`. Esser approximates `f' ~ f[x-h, x+h], f'' ~
f[x-h,x,x+h]`, where `h = fx`, as with Steffensen's method, Requiring 3
function calls per step. The implementation `Order2B` uses a secant
step when `|fx|` is considered too large.

Esser, H. Computing (1975) 14: 367. DOI: [10.1007/BF02253547](https://doi.org/10.1007/BF02253547)
Eine stets quadratisch konvergente Modifikation des Steffensen-Verfahrens

Example
```
f(x) = cos(x) - x
g(x) = f(x)^2
x0 = pi/4
find_zero(f, x0, Order2(), verbose=true)        #  3 steps / 7 function calls
find_zero(f, x0, Roots.Order2B(), verbose=true) #  4 / 9
find_zero(g, x0, Order2(), verbose=true)        #  22 / 45
find_zero(g, x0, Roots.Order2B(), verbose=true) #  4 / 10
```
"""
struct Esser <: AbstractSecant end
struct Order2B <: AbstractSecant end

function update_state(method::Order2B, fs, o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}
     update_state_guarded(method, Secant(), Esser(), fs, o, options, l)
end

function update_state(method::Esser, F,
                      o::AbstractUnivariateZeroState{T,S}, options, l=NullTracks())  where {T, S}

    x1, fx1 = o.xn1, o.fxn1

    f0 = fx1

    f1  = F(x1 + f0 * oneunit(T) / oneunit(S))
    f_1 = F(x1 - f0 * oneunit(T) / oneunit(S))
    incfn(l, 2)

    # h = f0
    # r1 = f/f' ~ f/f[x+h,x-h]
    # r2 = f'/f'' ~ f[x+h, x-h]/f[x-h,x,x+h]
    r1 = f0 * 2*f0 / (f1 - f_1) * oneunit(T) / oneunit(S)
    r2 = (f1 - f_1)/(f1 - 2*f0 + f_1) * f0/2 * oneunit(T) / oneunit(S)

    k = r2/(r2-r1)  # ~ m

    if abs(k) <= 1e-2 * oneunit(k)
        log_message(l, "Estimate for multiplicity had issues. ")
        return o, true
    end

    delta = k * r1

    if isissue(delta)
        log_message(l, "Increment `Δx` has issues. ")
        return o, true
    end

    x0, x1 = x1, x1 - delta
    fx0, fx1 = fx1, F(x1)
    incfn(l)

    @set! o.xn0 = x0
    @set! o.xn1 = x1
    @set! o.fxn0 = fx0
    @set! o.fxn1 = fx1

    return o, false

end

"""
    Roots.Schroder()

Schröder's method, like Halley's method, utilizes `f`, `f'`, and
`f''`. Unlike Halley it is quadratically converging, but this is
independent of the multiplicity of the zero (cf. Schröder, E. "Über
unendlich viele Algorithmen zur Auflösung der Gleichungen."
Math. Ann. 2, 317-365, 1870;
[mathworld](http://mathworld.wolfram.com/SchroedersMethod.html)). Schröder's
method applies Newton's method to `f/f'`, a function with all
simple zeros.

Example
```
m = 2
f(x) = (cos(x)-x)^m
fp(x) = (-x + cos(x))*(-2*sin(x) - 2)
fpp(x) = 2*((x - cos(x))*cos(x) + (sin(x) + 1)^2)
find_zero((f, fp, fpp), pi/4, Roots.Halley())    # 14 steps
find_zero((f, fp, fpp), 1.0, Roots.Schroder())    # 3 steps
```

(Whereas, when `m=1`, Halley is 2 steps to Schröder's 3.)

If function evaluations are expensive one can pass in a function which
returns `(f, f/f',f'/f'')` as follows

```
find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.Schroder())
```

This can be advantageous if the derivatives are easily computed from
the value of `f`, but otherwise would be expensive to compute.

The error, `eᵢ = xᵢ - α`, is the same as `Newton` with `f` replaced by `f/f'`.

"""
struct Schroder <: AbstractHalleyLikeMethod
end
const Schroeder = Schroder # either spelling
const Schröder = Schroder

## Shroder-like methods
function init_state(M::Schroder, F::Callable_Function, x)
    x₀ = float(first(x))
    fx₀, (Δ, ΔΔ) = F(x₀)
    x₁ = x₀ - ΔΔ/(ΔΔ - Δ) * Δ # m*r1
    state = init_state(M, F, x₀, x₁, fx₀, fx₀)
end

function update_state(method::Schroder, F, o::AbstractUnivariateZeroState{T,S},
                      options::UnivariateZeroOptions, l=NullTracks()) where {T,S}
    xn = o.xn1
    fxn = o.fxn1
    r1, r2 = o.Δ, o.ΔΔ

    Δ =  r2 / (r2 - r1) * r1  # m * r1

    if isissue(Δ)
        log_message(l, "Issue with increment")
        return o, true
    end

    xn1::T = xn - Δ

    fxn1::S, (r1::T, r2::T) = F(xn1)
    incfn(l,3)

    @set! o.xn0 = xn
    @set! o.xn1 = xn1
    @set! o.fxn0 = fxn
    @set! o.fxn1 = fxn1
    @set! o.Δ = r1
    @set! o.ΔΔ = r2

    return o, false
end

## Thukral 3,4,5 (2) is Schroder
"""
    AbstractThukralBMethod

Abstract type for `ThukralXB` methods for `X` being `2`,`3`,`4`, or `5`.

These are a family of methods which are
* efficient (order `X`) for non-simple roots (e.g. `Thukral2B` is the `Schroder` method)
* taking `X+1` function calls per step
* require `X` derivatives. These can be passed as a tuple of functions, `(f, f', f'', …)`, *or* as
a function returning the ratios: `x -> (f(x), f(x)/f'(x), f'(x)/f''(x), …)`.

## Example

```julia
using ForwardDiff
Base.adjoint(f::Function)  = x  -> ForwardDiff.derivative(f, float(x))
f(x) = (exp(x) + x - 2)^6
x0 = 1/4
find_zero((f, f', f''), x0, Roots.Halley())               # 14 iterations; 45 function evaluations
find_zero((f, f', f''), big(x0), Roots.Thukral2B())       #  4 iterations; 15 function evaluations
find_zero((f, f', f'', f'''), big(x0), Roots.Thukral3B()) #  3 iterations; 16 function evaluations
```

## Refrence

*Introduction to a family of Thukral ``k``-order method for finding multiple zeros of nonlinear equations*,
R. Thukral, JOURNAL OF ADVANCES IN MATHEMATICS 13(3):7230-7237, DOI: [10.24297/jam.v13i3.6146](https://doi.org/10.24297/jam.v13i3.6146).
"""
abstract type AbstractThukralBMethod <: AbstractHalleyLikeMethod end
initial_fncalls(M::AbstractThukralBMethod) = fn_argout(M)

struct Thukral2B <: AbstractThukralBMethod end
fn_argout(::Thukral2B) = 3

struct Thukral3B <: AbstractThukralBMethod end
fn_argout(::Thukral3B) = 4

struct Thukral4B <: AbstractThukralBMethod end
fn_argout(::Thukral4B) = 5

struct Thukral5B <: AbstractThukralBMethod end
fn_argout(::Thukral5B) = 6

struct ThukralBState{N,T,S} <: AbstractUnivariateZeroState{T,S}
    xn1::T
    xn0::T
    Δs::NTuple{N, T}
    fxn1::S
    fxn0::S
end

function init_state(M::AbstractThukralBMethod, F::Callable_Function, x)
    x₁ = float(first(x))
    fx₁, Δs = F(x₁)
    state = init_state(M, F, nan(x₁), x₁, nan(fx₁), fx₁; Δs=Δs)
end

function init_state(M::AbstractThukralBMethod, F, x₀, x₁::T, fx₀, fx₁; Δs = nothing) where {T}
    ThukralBState(x₁, x₀, NTuple{fn_argout(M)-1,T}(Δs), fx₁, fx₀)
end

function update_state(M::AbstractThukralBMethod, F, o::AbstractUnivariateZeroState{T,S},
                      options::UnivariateZeroOptions, l=NullTracks()) where {T,S}

    x₀ = o.xn1

    Δ = compute_thukral_Δ(M, o)

    isissue(Δ) && return (o, true)
    x₁ = x₀ - Δ
    fx₁, Δs = F(x₁)
    incfn(l, fn_argout(M))

    @set! o.xn0 = x₀
    @set! o.fxn0 = o.fxn1
    @set! o.Δs =  NTuple{fn_argout(M)-1,T}(Δs)
    @set! o.xn1 = x₁
    @set! o.fxn1 = fx₁

    return (o, false)
end

function compute_thukral_Δ(M::Thukral2B, o)
    r₁, r₂ = o.Δs
    t₁, t₂ = 1/r₁, 1/r₂
    Δ = one(o.xn1)
    Δ /= (t₁ - t₂)
    Δ
end

function compute_thukral_Δ(M::Thukral3B, o)
    r₁, r₂, r₃ = o.Δs
    t₁, t₂, t₃ = 1/r₁, 1/r₂, 1/r₃
    Δ = (2t₁ - 2t₂)
    Δ /= (2t₁^2 - 3t₁*t₂ + t₂*t₃)
    Δ
end

function compute_thukral_Δ(M::Thukral4B, o)
    r₁, r₂, r₃, r₄ = o.Δs
    t₁, t₂, t₃, t₄ = 1/r₁, 1/r₂, 1/r₃, 1/r₄
    Δ = 6t₁^2 - 9t₁*t₂ +3t₂*t₃
    Δ /= 6t₁^3 - 12 * t₁^2*t₂ + 4t₁*t₂*t₃ - t₂*t₃*t₄ + 3*t₁*t₂^2
    Δ
end

function compute_thukral_Δ(M::Thukral5B, o)
    r₁, r₂, r₃, r₄, r₅ = o.Δs
    t₁, t₂, t₃, t₄, t₅ = 1/r₁, 1/r₂, 1/r₃, 1/r₄, 1/r₅
    Δ = 24*t₁^3 - 48t₁^2*t₂  + 16*t₁*t₂*t₃ - 4*t₂*t₃*t₄ + 12t₁*t₂^2
    Δ /= 24t₁^4  - 60t₁^3*t₂ + 20*t₁^2*t₂*t₃ - 5*t₁*t₂*t₃*t₄ + 30t₁^2*t₂^2 - 10*t₁*t₂^2*t₃ + t₂*t₃*t₄*t₅
    Δ
end


text:
const test = "3,4,3,1,2"

function simulate(fish::AbstractVector{<: Integer},
                  days::Integer;
                  fertility_cycle::Integer = 7,
                  spawn_age::Integer = 8
                 )::Integer
    ages = [count(==(age), fish) for age in 0:spawn_age]
    for day in 1:days
        spawning = ages[1]
        ages = [ages[2:end]; spawning]
        ages[fertility_cycle] += spawning
    end
    return sum(ages)
end

let v = parse.(Int, split(test, ','))
    @assert simulate(v, 80) == 5934
    @assert simulate(v, 256) == 26984457539
end

v = parse.(Int, split(readline("day06.in"), ','))
println("Q1: ", simulate(v, 80))
println("Q2: ", simulate(v, 256))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
include("run_scf_and_compare.jl")
include("testcases.jl")

function run_oxygen_pbe(T; kwargs...)
    # These values were computed using ABINIT at Ecut=13. They are not yet converged and
    # thus require the same discretization parameters to be obtained.
    ref_evals = [
        [ -1.214034187039492, -0.617266880915023, -0.370358278719988, -0.370358278712184,
          -0.352666293942291,  0.008939230950524,  0.008939230958503,  0.098061783836035],
        [ -1.177293790864051, -0.551559952680091, -0.318494822087423, -0.308762934378086,
          -0.308762934370473,  0.090371516203330,  0.090371516211105,  0.115096531941182],
    ]
    ref_etot = -29.505154342215235
    ref_magn = 1.99231275

    # Produce reference data and guess for this configuration
    O = ElementPsp(o2molecule.atnum, psp=load_psp("hgh/pbe/O-q6.hgh"))
    magnetic_moments = [O => [1., 1.]]
    model = model_PBE(Array{T}(o2molecule.lattice), [O => o2molecule.positions],
                      temperature=0.02, smearing=smearing=Smearing.Gaussian(),
                      magnetic_moments=magnetic_moments)
    basis = PlaneWaveBasis(model; Ecut=13, fft_size=[24, 24, 30], kgrid=[1, 1, 1])

    scfres = run_scf_and_compare(T, basis, ref_evals, ref_etot;
                                 ρ=guess_density(basis, magnetic_moments),
                                 test_etot=false,
                                 kwargs...)
    @test scfres.energies.total ≈ ref_etot atol=1e-4  # A little large a difference ...

    magnetization = sum(spin_density(scfres.ρ)) * basis.dvol
    @test magnetization ≈ ref_magn atol=1e-4
end

@testset "Oxygen PBE (Float64)" begin
    run_oxygen_pbe(Float64, test_tol=5e-5, scf_tol=1e-9, n_ignored=1)
end


text:
# # Gaussian process prior samples
#
# The kernels defined in this package can also be used to specify the
# covariance of a Gaussian process prior.
# A Gaussian process (GP) is defined by its mean function $m(\cdot)$ and its covariance function or kernel $k(\cdot, \cdot')$:
# ```math
#   f \sim \mathcal{GP}\big(m(\cdot), k(\cdot, \cdot')\big)
# ```
# In this notebook we show how the choice of kernel affects the samples from a GP (with zero mean).

## Load required packages
using KernelFunctions, LinearAlgebra
using Plots, Plots.PlotMeasures
default(; lw=1.0, legendfontsize=8.0)
using Random: seed!
seed!(42); # reproducibility

##
# ## Evaluation at finite set of points
#
# The function values $\mathbf{f} = \{f(x_n)\}_{n=1}^N$ of the GP at a finite number $N$ of points $X = \{x_n\}_{n=1}^N$ follow a multivariate normal distribution $\mathbf{f} \sim \mathcal{MVN}(\mathbf{m}, \mathrm{K})$ with mean vector $\mathbf{m}$ and covariance matrix $\mathrm{K}$, where
# ```math
# \begin{aligned}
#   \mathbf{m}_i &= m(x_i) \\
#   \mathrm{K}_{i,j} &= k(x_i, x_j)
# \end{aligned}
# ```
# with $1 \le i, j \le N$.
#
# We can visualize the infinite-dimensional GP by evaluating it on a fine grid to approximate the dense real line:
num_inputs = 101
xlim = (-5, 5)
X = range(xlim...; length=num_inputs);

# Given a kernel `k`, we can compute the kernel matrix as `K = kernelmatrix(k, X)`.

##
# ## Random samples
#
# To sample from the multivariate normal distribution $p(\mathbf{f}) = \mathcal{MVN}(0, \mathrm{K})$, we could make use of Distributions.jl and call `rand(MvNormal(K))`.
# Alternatively, we could use the [AbstractGPs.jl](https://github.com/JuliaGaussianProcesses/AbstractGPs.jl) package and construct a `GP` object which we evaluate at the points of interest and from which we can then sample: `rand(GP(k)(X))`.

# Here, we will explicitly construct samples using the Cholesky factorization $\mathrm{L} = \operatorname{cholesky}(\mathrm{K})$,
# with $\mathbf{f} = \mathrm{L} \mathbf{v}$, where $\mathbf{v} \sim \mathcal{N}(0, \mathbf{I})$ is a vector of standard-normal random variables.

# We will use the same randomness $\mathbf{v}$ to generate comparable samples across different kernels.
num_samples = 7
v = randn(num_inputs, num_samples);

# Mathematically, a kernel matrix is by definition positive semi-definite, but due to finite-precision inaccuracies, the computed kernel matrix might not be exactly positive definite. To avoid Cholesky errors, we add a small "nugget" term on the diagonal:
function mvn_sample(K)
    L = cholesky(K + 1e-6 * I)
    f = L.L * v
    return f
end;

##
# ## Visualization
# We now define a function that visualizes a kernel for us.

function visualize(k::Kernel)
    K = kernelmatrix(k, X)
    f = mvn_sample(K)

    p_kernel_2d = heatmap(
        X,
        X,
        K;
        yflip=true,
        colorbar=false,
        ylabel=string(nameof(typeof(k))),
        ylim=xlim,
        yticks=([xlim[1], 0, xlim[end]], ["\u22125", raw"$x'$", "5"]),
        vlim=(0, 1),
        title=raw"$k(x, x')$",
        aspect_ratio=:equal,
        left_margin=5mm,
    )

    p_kernel_cut = plot(
        X,
        k.(X, 0.0);
        title=string(raw"$k(x, x_\mathrm{ref})$"),
        label=raw"$x_\mathrm{ref}=0.0$",
        legend=:topleft,
        foreground_color_legend=nothing,
    )
    plot!(X, k.(X, 1.5); label=raw"$x_\mathrm{ref}=1.5$")

    p_samples = plot(X, f; c="blue", title=raw"$f(x)$", ylim=(-3, 3), label=nothing)

    return plot(
        p_kernel_2d,
        p_kernel_cut,
        p_samples;
        layout=(1, 3),
        xlabel=raw"$x$",
        xlim=xlim,
        xticks=collect(xlim),
    )
end;

##
# We can now visualize a kernel and show samples from
# a Gaussian process with a given kernel:

plot(visualize(SqExponentialKernel()); size=(800, 210), bottommargin=5mm, topmargin=5mm)

##
# ## Kernel comparison
# This also allows us to compare different kernels:

kernels = [
    Matern12Kernel(),
    Matern32Kernel(),
    Matern52Kernel(),
    SqExponentialKernel(),
    WhiteKernel(),
    ConstantKernel(),
    LinearKernel(),
    compose(PeriodicKernel(), ScaleTransform(0.2)),
    NeuralNetworkKernel(),
]
plot(
    [visualize(k) for k in kernels]...;
    layout=(length(kernels), 1),
    size=(800, 220 * length(kernels) + 100),
)


text:
function update!(mc::MarkovChain{T},
                 mcmc::MCMC{T},
                 Σ::Array{Float64, 2},
                 σ::Float64;
                 condition_on_network::Bool=false,
                 event_batches::Int64=1,
                 transmission_rate_cache::Union{Nothing, TransmissionRateCache}=nothing) where T <: EpidemicModel
  # Initialize
  new_tr_cache = transmission_rate_cache
  new_events = mc.events[end]
  new_events_array = new_events[_state_progressions[T][2:end]]
  new_params = mc.risk_parameters[end]
  new_network = mc.transmission_network[end]
  new_lposterior = mc.log_posterior[end]
  # Randomize event time augmentation
  event_indices = findall(new_events_array[:] .> -Inf)
  aug_order = sample(event_indices, length(event_indices), replace=false)
  @debug "$(length(aug_order)) events to augment"
  if event_batches < 0
    @error "Cannot have negative amount of event batches"
  end
  if event_batches > length(aug_order)
    @warn "More event batches than there are events to augment ($(event_batches) > $(length(aug_order))), setting to maximum ($(length(aug_order)))"
    event_batches = length(aug_order)
  end
  batch_size = fld(length(aug_order), event_batches)
  @debug "Performing data augmentation in batches of $batch_size events at a time"
  for i = 1:(1 + event_batches)
    if i == 1
      # Propose new risk parameters
      proposed_events = new_events
      proposed_events_array = new_events_array
      proposed_params = generate(RiskParameters{T}, new_params, Σ)
      @debug "Generating new TransmissionRateCache for TN-ILM parameter proposal"
      proposed_tr_cache = TransmissionRateCache(size(new_events_array, 1))
    else
      proposed_events_array = copy(new_events_array)
      proposed_events = Events{T}(proposed_events_array)
      for j = (batch_size*(i-2) + 1):min(batch_size*(i-1),length(aug_order))
        id, state_index = Tuple(CartesianIndices((new_events.individuals,
                                            length(_state_progressions[T][2:end])))[aug_order[j]])
        new_state = _state_progressions[T][state_index+1]
        time = new_events[new_state][id]

        # Conditioning on network means that only event times valid under current network will be proposed.
        # This is useful for models which may have additional contributions to the posterior based on network,
        # and require more modest proposals (e.g. phylodynamic models).
        if condition_on_network
          proposed_event = generate(Event,
                                    Event{T}(time, id, new_state),
                                    σ,
                                    mcmc.event_extents,
                                    mcmc.event_observations,
                                    proposed_events,
                                    new_network)
        else
          proposed_event = generate(Event,
                                    Event{T}(time, id, new_state),
                                    σ,
                                    mcmc.event_extents,
                                    mcmc.event_observations,
                                    proposed_events)
        end
        proposed_events_array = reshape([proposed_events_array[1:(aug_order[j]-1)]
                                         _time(proposed_event)
                                         proposed_events_array[(aug_order[j]+1):end]],
                                        size(new_events_array))
        proposed_events = Events{T}(proposed_events_array)
      end
      proposed_params = new_params
      proposed_tr_cache = new_tr_cache
    end
    proposed_lprior = logpriors(proposed_params, mcmc.risk_priors)
    # Based on the logprior and competiting MCMC iteration, this loglikelihood is required for acceptance
    # Calculating this in advance allows us to cut loglikelihood calculation short if it goes below threshold
    ll_acceptance_threshold = log(rand()) + new_lposterior - proposed_lprior
    if ll_acceptance_threshold < Inf
      proposed_llikelihood = loglikelihood(proposed_params,
                                           mcmc.risk_functions,
                                           proposed_events,
                                           mcmc.population,
                                           mcmc.starting_states,
                                           transmission_rates_output = false,
                                           transmissions_output = false,
                                           early_decision_value = ll_acceptance_threshold,
                                           transmission_rate_cache = proposed_tr_cache)[1]
      proposed_lposterior = proposed_lprior + proposed_llikelihood
    else
      proposed_llikelihood = -Inf
      proposed_lposterior = -Inf
    end
    if proposed_llikelihood >= ll_acceptance_threshold
      @debug "MCMC proposal accepted (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
      new_params = proposed_params
      new_tr_cache = proposed_tr_cache
      new_events = proposed_events
      new_events_array = proposed_events_array
      new_lposterior = proposed_lposterior
    else
      @debug "next!: MCMC proposal rejected (acceptance probability = $(round(min(1.0, exp(proposed_lposterior - new_lposterior)), digits=3)))"
    end
  end
  mc.iterations += 1
  tnr, tx = loglikelihood(new_params,
                          mcmc.risk_functions,
                          new_events,
                          mcmc.population,
                          mcmc.starting_states,
                          loglikelihood_output = false,
                          transmission_rate_cache = new_tr_cache)[[2; 3]]
  new_network = generate(TransmissionNetwork, tnr, mcmc, tx)
  push!(mc.events, new_events)
  push!(mc.transmission_network, new_network)
  push!(mc.risk_parameters, new_params)
  push!(mc.log_posterior, new_lposterior)
  return mc
end

function update!(mcmc::MCMC{T},
                 Σ::Array{Float64, 2},
                 σ::Float64) where T <: EpidemicModel
  @simd for mc in mcmc.markov_chains
    next!(mc, mcmc, Σ, σ)
  end
  return mcmc
end


text:
@testset "from glm" begin
    fdf = DataFrame(Carb=[0.1,0.3,0.5,0.6,0.7,0.9], OptDen=[0.086,0.269,0.446,0.538,0.626,0.782])
    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, req_stats="all")
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2
    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.ADJR2, 0.998808343507198)
    @test isapprox(lm1.AIC, -55.43694668654871) # using the SAS formula rather than the Julia-Statsmodel-GLM
    @test isapprox(lm1.stderrors,  [0.007833679251299831, 0.013534536322659505])
    @test isapprox(lm1.t_values, [0.6492114525712505, 64.74442074669666])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.ci_low, [-0.016664066127247305, 0.8387078171615459])
    @test isapprox(lm1.ci_up, [0.02683549469867456, 0.9138636114098853])
    @test isapprox(lm1.t_statistic, 2.7764451051977934)
    @test isapprox(lm1.VIF, [0.,  1.])
    @test isapprox(lm1.PRESS, 0.0010755278041106075)
    @test isapprox(lm1.Type1SS, [1.2576681666666665, 0.3135496333333334])
    @test isapprox(lm1.Type2SS, [3.152636815919868e-5, 0.31354963333333363])
end

@testset "from glm regresspredict" begin
    fdf = DataFrame(Carb=[0.1,0.3,0.5,0.6,0.7,0.9], OptDen=[0.086,0.269,0.446,0.538,0.626,0.782])
    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=[:default, :aic, :vif, :press, :t1ss, :t2ss])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2
    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.ADJR2, 0.998808343507198)
    @test isapprox(lm1.AIC, -55.43694668654871) # using the SAS formula rather than the Julia-Statsmodel-GLM
    @test isapprox(lm1.stderrors,  [0.007833679251299831, 0.013534536322659505])
    @test isapprox(lm1.t_values, [0.6492114525712505, 64.74442074669666])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])
    @test isapprox(lm1.ci_low, [-0.016664066127247305, 0.8387078171615459])
    @test isapprox(lm1.ci_up, [0.02683549469867456, 0.9138636114098853])
    @test isapprox(lm1.t_statistic, 2.7764451051977934)
    @test isapprox(lm1.VIF, [0.,  1.])
    @test isapprox(lm1.PRESS, 0.0010755278041106075)
    @test isapprox(lm1.Type1SS, [1.2576681666666665, 0.3135496333333334])
    @test isapprox(lm1.Type2SS, [3.152636815919868e-5, 0.31354963333333363])

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=["none"])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2

    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf, α=0.05, req_stats=["r2", "P_values"])
    target_coefs = [0.005085714285713629, 0.8762857142857156]
    @test isapprox(target_coefs, lm1.coefs)
    @test lm1.p == 2
    @test isapprox(lm1.R2, 0.9990466748057584)
    @test isapprox(lm1.p_values, [0.5515952883836446, 3.409192065429258e-7])

    df = DataFrame(y=[1., 3., 3., 2., 2., 1.], x1=[1., 2., 3., 1., 2., 3.], x2=[1., 1., 1., -1., -1., -1.])
    lm2 = regress(@formula(y ~ 1 + x1 + x2), df, req_stats=["r2"])
    @test [1.5, 0.25, 0.3333333333333333] == lm2.coefs
    @test 0.22916666666666663 == lm2.R2

    y = [3.547744106900422, 9.972950249405148, 16.471345464154027, 22.46768807351274, 20.369933318011807, 21.18590757820348, 29.962620198209024, 30.684400502954748, 29.28429078492597, 34.272759386588824, 33.05504692986838, 45.09273876302829, 45.28374262744938, 54.54563960566191, 46.86173948296966, 46.85926120310666, 67.54337216713414, 66.0400205145086, 64.77001443647681, 63.98759256558095, 15.016939388490687, 12.380885701920008, 10.221963402745288, 24.826987790646672, 22.231511892187548, 18.05125492502642, 21.809661866284717, 28.533306224702308, 24.70476800009685, 39.592181057942916, 37.425282624708906, 33.89372811236659, 44.07442335927311, 42.61943719318272, 52.32531447897055, 44.12096163880534, 49.19965880584543, 52.304114239221036, 59.79488104919937, 65.00419916894333, 10.978506745605673, 12.944637189081874, 16.96525561080181, 15.93489355798633, 33.559135307088305, 28.887571687420433, 28.296418324622593, 32.2405215537489, 31.466024917490223, 37.78855308849255, 48.09725215994402, 40.70190542438069, 43.2525436240948, 44.75159242105558, 52.553066965996074, 53.265851121437095, 60.092700204726015, 62.78014347092756, 69.07895282754228, 74.27890668517966, 10.57283140105129, 7.290600131361426, 16.10299402050687, 14.428954773831242, 22.418180226673464, 28.21022852582732, 29.60436622143203, 33.648588929669636, 37.451930576147994, 43.85548900583812, 41.44168340404242, 43.48815266671309, 54.72835160956764, 53.55177468229062, 50.62088969800169, 50.863408563713335, 47.40251347184729, 64.29413401802115, 64.20687412126479, 73.66653216085827, 15.820723594639162, 21.973463928234743, 23.804440715385162, 15.139822408433913, 30.015089890369985, 28.08454394421385, 33.041880065463566, 28.49429418531324, 38.33763660905526, 34.503176013521724, 48.748946870573235, 45.45351516824085, 53.522908096188765, 45.95216131836423, 63.13018379633756, 63.4236036208151, 65.28265579397677, 55.43500146374922, 76.50187470137375, 67.22421998359037, 19.289152114521762, 27.63557010588222, 17.700031078096686, 25.462368278660957, 22.94613580060685, 36.19731649621813, 35.22216995579936, 25.526434727578838, 45.882864925557726, 38.71433797181679, 50.617762276434554, 41.96951039650285, 52.265123328453214, 45.383991138243765, 62.7923270665056, 64.40670612696276, 74.89775274405821, 72.89056716347118, 66.37071343209973, 76.32913721410918, 29.923781174452596, 16.465832617450324, 30.530915733275403, 30.512411156491954, 28.04012351007911, 26.315140074869376, 37.18491928428231, 41.958551085353626, 48.370736387628895, 49.419917216561835, 48.67296268029715, 55.484477112881166, 51.120639229597025, 54.797092949987224, 61.92608065418044, 69.79495109420618, 64.43521939892794, 74.35280205157255, 77.22355341160723, 68.90654715174155, 10.87752158238755, 20.25006748725279, 31.71957876614495, 31.42898857400639, 32.505052996368256, 34.1225641368903, 35.66496173153403, 31.76371648643483, 42.97107454773545, 47.41827763710622]
    x2 = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
    x3 = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
    df = DataFrame(y=y, x2=x2, x3=x3)
    lrint = regress(@formula(y ~ x2 +x3), df, req_stats=["default", "pcorr1", "pcorr2", "scorr1", "scorr2", "vif"])

end

@testset "weighted regression" begin
    tw = [
        2.3  7.4  0.058
        3.0  7.6  0.073
        2.9  8.2  0.114
        4.8  9.0  0.144
        1.3 10.4  0.151
        3.6 11.7  0.119
        2.3 11.7  0.119
        4.6 11.8  0.114
        3.0 12.4  0.073
        5.4 12.9  0.035
        6.4 14.0  0
    ] # data from https://blogs.sas.com/content/iml/2016/10/05/weighted-regression.html

    df = DataFrame(tw, [:y,:x,:w])
    f = @formula(y ~ x)
    lm = regress(f, df, weights="w", req_stats=["default", "t1ss", "t2ss"])

    @test isapprox([2.328237176867885, 0.08535712911515277], lm.coefs)
    @test isapprox(0.014954934572439349, lm.R2)
    @test isapprox(-0.10817569860600562, lm.ADJR2)
    @test isapprox([2.551864989438224, 0.24492357920520605], lm.stderrors)
    @test isapprox([0.3882424860021164, 0.7364546437428148], lm.p_values)
    @test isapprox([10.272024999999998, 0.02220919946016564], lm.Type1SS)
    @test isapprox([0.15221363313265734, 0.022209199460165682], lm.Type2SS)

    res = predict_in_sample(lm, df, req_stats="all")
    @test isapprox([2.9598799323200153, 2.976951358143046, 3.0281656356121376, 3.09645133890426, 3.2159513196654737, 3.326915587515172, 3.326915587515172, 3.335451300426688, 3.386665577895779, 3.4293441424533553], res.predicted)
    @test isapprox([0.2149107957203927, 0.24394048660618184, 0.27451117660671137, 0.22039739135761363, 0.15181541173049415, 0.19864023422787455, 0.19864023422787455, 0.20135117925989024, 0.18147639530421758, 0.1143166949587325], res.leverage)
    @test isapprox([1.957110958921382, 1.765206920504513, 1.4298044633303422, 1.244877562198597, 1.1810280676815839, 1.3571510747551465, 1.3571510747551465, 1.388160919947073, 1.7203164594367346, 2.412834568973736], res.stdi)
    @test isapprox([0.8231374655094712, 0.7816957716988308, 0.6635670917274183, 0.5290286945664605, 0.4287722865207785, 0.5524810400784785, 0.5524810400784785, 0.568305570562634, 0.6742266146509985, 0.7728194750528609], res.stdp)
    @test isapprox([1.5732681689017483, 1.3761754659974743, 1.0787484567368795, 0.9949760929369788, 1.0134771577471209, 1.1096794313749598, 1.1096794313749598, 1.1318340411046413, 1.4318958288930972, 2.151109196482762], res.stdr)
    @test isapprox([0.02407872915271305, 4.5252332002240634e-5, 0.0026705586843115184, 0.414368722166146, 0.31984274220916814, 0.0075059979227238794, 0.10614122824505004, 0.15735260690358094, 0.008083642818377243, 0.054162396113877506], res.cooksd)
    @test isapprox([1.0617215330136984, 1.1743576761339307, 1.4979771781033089, 1.8765089815923404, 2.2272006538876132, 2.0528920244723023, 2.0528920244723023, 2.024936304649691, 1.8318962164458708, 1.6472192372151455], res.lclp)
    @test isapprox([4.858038331626332, 4.779545040152161, 4.5583540931209665, 4.3163936962161795, 4.204701985443334, 4.600939150558042, 4.600939150558042, 4.645966296203685, 4.941434939345688, 5.2114690476915655], res.uclp)
    @test isapprox([-1.5532260320060822, -1.093623100031372, -0.2689693693610047, 0.225758532651414, 0.49249571179955565, 0.19731959703302637, 0.19731959703302637, 0.13434647869991867, -0.5803912914251206, -2.134662351163641], res.lcli)
    @test isapprox([7.472985896646113, 7.047525816317464, 6.325300640585279, 5.967144145157105, 5.939406927531392, 6.456511577997318, 6.456511577997318, 6.536556122153457, 7.353722447216679, 8.993350636070351], res.ucli)
    @test isapprox([-0.41943258330883165, 0.01674833073720578, -0.11880956567004286, 1.7121503453084896, -1.8904731152742318, 0.24609306504533465, -0.9254164387302033, 1.117256288156119, -0.2700375055877384, 0.9161114929771242], res.student)
    @test isapprox([-0.39672963851269843, 0.015666903524005342, -0.1112343500884851, 2.0120993810497283, -2.3774333761244986, 0.23107529120204373, -0.9160676175777195, 1.1376116334771753, -0.25375610179858576, 0.9057709902138592], res.rstudent)
    @test isapprox([-0.8405158658696843, 0.03048522166395697, -0.1766610752356951, 2.1851500267069777, -2.2588848537963426, 0.3407762956775175, -1.2814663667643953, 1.5833601286751433, -0.47239392447269685, 2.225011859577532], res.press)
end

@testset "predictions statistics" begin
    t_carb = [0.1, 0.3, 0.5, 0.6, 0.7, 0.9]
    t_optden = [0.086, 0.269, 0.446, 0.538, 0.626, 0.782]
    t_leverage = [0.5918367346938774, 0.2816326530612245, 0.1673469387755102, 0.1836734693877552, 0.24897959183673485, 0.5265306122448984]
    t_predicted = [0.09271428571428536, 0.26797142857142836, 0.44322857142857136, 0.5308571428571429, 0.6184857142857143, 0.7937428571428573]
    t_residuals = [-0.006714285714285367, 0.0010285714285716563, 0.002771428571428647, 0.0071428571428571175, 0.007514285714285696, -0.011742857142857277]
    t_stdp = [0.0066535244611503905, 0.00458978457544483, 0.0035380151243903845, 0.003706585424647827, 0.004315515434962329, 0.006275706318490394]
    t_stdi = [0.01091189203370195, 0.009791124677432761, 0.009344386069745653, 0.009409504530540022, 0.009665592246181276, 0.010685714285717253]
    t_stdr = [0.00552545131594831, 0.00733033952495041, 0.007891923021648551, 0.007814168189246363, 0.0074950868260910365, 0.005951093194035971]
    t_student = [-1.2151560714878948, 0.1403170242074994, 0.3511727830880084, 0.9140905301586566, 1.0025615297914663, -1.9732268946192373]
    t_rstudent = [-1.3249515797718379, 0.12181828539462737, 0.3089239861926457, 0.8900235681358146, 1.003419757326401, -10.478921731163984]
    t_lcli = [0.06241801648886679, 0.24078690838638886, 0.4172843964641476, 0.5047321700609886, 0.5916497280049665, 0.7640745580187355]
    t_ucli = [0.12301055493970393, 0.2951559487564679, 0.4691727463929951, 0.5569821156532972, 0.6453217005664621, 0.8234111562669791]
    t_lclp = [0.07424114029181057, 0.2552281436530222, 0.43340546665434193, 0.520566011897882, 0.6065039225799076, 0.7763187030532258]
    t_uclp = [0.11118743113676015, 0.2807147134898345, 0.4530516762028008, 0.5411482738164038, 0.630467505991521, 0.8111670112324888]
    t_press = [-0.016449999999999146, 0.001431818181818499, 0.0033284313725491102, 0.00874999999999997, 0.010005434782608673, -0.02480172413793134]
    t_cooksd = [1.0705381016035724, 0.0038594654616162225, 0.012392684477580572, 0.09400066844914513, 0.16661116000566892, 2.1649894168822432]

    fdf = DataFrame([[0.1,0.3,0.5,0.6,0.7,0.9],[0.086,0.269,0.446,0.538,0.626,0.782]], [:Carb, :OptDen])
    lm1 = regress(@formula(OptDen ~ 1 + Carb), fdf)
    results = predict_in_sample(lm1, fdf, α=0.05, req_stats=["all"])
    @test isapprox(t_leverage, results.leverage)
    @test isapprox(t_predicted, results.predicted)
    @test isapprox(t_residuals, results.residuals)
    @test isapprox(t_stdp, results.stdp)
    @test isapprox(t_stdi, results.stdi)
    @test isapprox(t_stdr, results.stdr)
    @test isapprox(t_student, results.student)
    @test isapprox(t_rstudent, results.rstudent)
    @test isapprox(t_lcli, results.lcli)
    @test isapprox(t_ucli, results.ucli)
    @test isapprox(t_lclp, results.lclp)
    @test isapprox(t_uclp, results.uclp)
    @test isapprox(t_press, results.press)
    @test isapprox(t_cooksd, results.cooksd)

    results = predict_in_sample(lm1, fdf, α=0.05, req_stats=["none"])
    @test isapprox(t_predicted, results.predicted)
    @test_throws ArgumentError("column name :leverage not found in the data frame") t_leverage == results.leverage

    results = predict_out_of_sample(lm1, fdf)
    @test isapprox(t_predicted, results.predicted)
    @test_throws ArgumentError("column name :leverage not found in the data frame") t_leverage == results.leverage

end


text:
using Test
using Flux: onehotbatch, mse, crossentropy, logitcrossentropy,
            σ, binarycrossentropy, logitbinarycrossentropy, flatten,
            xlogx, xlogy

const ϵ = 1e-7

@testset "xlogx & xlogy" begin
  @test iszero(xlogx(0))
  @test isnan(xlogx(NaN))
  @test xlogx(2) ≈ 2.0 * log(2.0)
  @inferred xlogx(2)
  @inferred xlogx(0)

  @test iszero(xlogy(0, 1))
  @test isnan(xlogy(NaN, 1))
  @test isnan(xlogy(1, NaN))
  @test isnan(xlogy(NaN, NaN))
  @test xlogy(2, 3) ≈ 2.0 * log(3.0)
  @inferred xlogy(2, 3)
  @inferred xlogy(0, 1)
end

@testset "losses" begin
  # First, regression-style y's
  y = [1, 1, 0, 0]
  ŷ = [.9, .1, .1, .9]

  @testset "mse" begin
    @test mse(ŷ, y) ≈ (.1^2 + .9^2)/2
  end

  @testset "mae" begin
    @test Flux.mae(ŷ, y) ≈ 1/2
  end

  @testset "huber_loss" begin
    @test Flux.huber_loss(ŷ, y) ≈ 0.20500000000000002
  end

  y = [123.0,456.0,789.0]
  ŷ = [345.0,332.0,789.0]
  @testset "msle" begin
    @test Flux.msle(ŷ, y) ≈ 0.38813985859136585
  end

  # Now onehot y's
  y = onehotbatch([1, 1, 0, 0], 0:1)
  ŷ = [.1 .9; .9 .1; .9 .1; .1 .9]'
  v = log(.1 / .9)
  logŷ = [v 0.0; 0.0 v; 0.0 v; v 0.0]'
  lossvalue = 1.203972804325936

  @testset "crossentropy" begin
    @test crossentropy([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ crossentropy([0.1,0.9], [0.1,0.9])
    @test crossentropy(ŷ, y) ≈ lossvalue
  end

  @testset "logitcrossentropy" begin
    @test logitcrossentropy(logŷ, y) ≈ lossvalue
  end

  @testset "weighted_crossentropy" begin
    @test crossentropy(ŷ, y, weight = ones(2)) ≈ lossvalue
    @test crossentropy(ŷ, y, weight = [.5, .5]) ≈ lossvalue/2
    @test crossentropy(ŷ, y, weight = [2, .5]) ≈ 1.5049660054074199
  end

  @testset "weighted_logitcrossentropy" begin
    @test logitcrossentropy(logŷ, y, weight = ones(2)) ≈ lossvalue
    @test logitcrossentropy(logŷ, y, weight = [.5, .5]) ≈ lossvalue/2
    @test logitcrossentropy(logŷ, y, weight = [2, .5]) ≈ 1.5049660054074199
  end

  logŷ, y = randn(3), rand(3)
  @testset "binarycrossentropy" begin
    @test binarycrossentropy.(σ.(logŷ), y; ϵ=0) ≈ -y.*log.(σ.(logŷ)) - (1 .- y).*log.(1 .- σ.(logŷ))
    @test binarycrossentropy.(σ.(logŷ), y) ≈ -y.*log.(σ.(logŷ) .+ eps.(σ.(logŷ))) - (1 .- y).*log.(1 .- σ.(logŷ) .+ eps.(σ.(logŷ)))
  end

  @testset "logitbinarycrossentropy" begin
    @test logitbinarycrossentropy.(logŷ, y) ≈ binarycrossentropy.(σ.(logŷ), y; ϵ=0)
  end

  y = [1 2 3]
  ŷ = [4.0 5.0 6.0]
  @testset "kldivergence" begin
    @test Flux.kldivergence([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ Flux.kldivergence([0.1,0.9], [0.1,0.9])
    @test Flux.kldivergence(ŷ, y) ≈ -1.7661057888493457
    @test Flux.kldivergence(y, y) ≈ 0
  end

  y = [1 2 3 4]
  ŷ = [5.0 6.0 7.0 8.0]
  @testset "hinge" begin
    @test Flux.hinge(ŷ, y) ≈ 0
    @test Flux.hinge(y, 0.5 .* y) ≈ 0.125
  end

  @testset "squared_hinge" begin
    @test Flux.squared_hinge(ŷ, y) ≈ 0
    @test Flux.squared_hinge(y, 0.5 .* y) ≈ 0.0625
  end

  y = [0.1 0.2 0.3]
  ŷ = [0.4 0.5 0.6]
  @testset "poisson" begin
    @test Flux.poisson(ŷ, y) ≈ 0.6278353988097339
    @test Flux.poisson(y, y) ≈ 0.5044459776946685
  end

  y = [1.0 0.5 0.3 2.4]
  ŷ = [0 1.4 0.5 1.2]
  @testset "dice_coeff_loss" begin
    @test Flux.dice_coeff_loss(ŷ, y) ≈ 0.2799999999999999
    @test Flux.dice_coeff_loss(y, y) ≈ 0.0
  end

  @testset "tversky_loss" begin
    @test Flux.tversky_loss(ŷ, y) ≈ -0.06772009029345383
    @test Flux.tversky_loss(ŷ, y, β = 0.8) ≈ -0.09490740740740744
    @test Flux.tversky_loss(y, y) ≈ -0.5576923076923075
  end

  @testset "no spurious promotions" begin
    for T in (Float32, Float64)
      y = rand(T, 2)
      ŷ = rand(T, 2)
      for f in (mse, crossentropy, logitcrossentropy, Flux.kldivergence, Flux.hinge, Flux.poisson,
              Flux.mae, Flux.huber_loss, Flux.msle, Flux.squared_hinge, Flux.dice_coeff_loss, Flux.tversky_loss)
        fwd, back = Flux.pullback(f, ŷ, y)
        @test fwd isa T
        @test eltype(back(one(T))[1]) == T
      end
    end
  end
end

@testset "helpers" begin
  @testset "flatten" begin
    x = randn(Float32, 10, 10, 3, 2)
    @test size(flatten(x)) == (300, 2)
  end
end


text:
module Formatters

using Niancat.Gameface

abstract type Formatter end

format(::Formatter, ::Response) = error("Implement format for Response in Formatter subtypes")
format(::Formatter, ::GameNotification) = error("Implement format for GameNotification in Formatter subtypes")

struct SlackFormatter <: Formatter end

# The default implementations just stringify the types, which isn't great
# but perhaps better than failing with a MethodError or error.
format(::SlackFormatter, response::Response) = string(response)
format(::SlackFormatter, notification::GameNotification) = string(notification)

export Formatter, SlackFormatter
export format

end


text:
abstract type AbstractOrder end

struct FirstOrder <: AbstractOrder end

struct SecondOrder <: AbstractOrder end

struct FirstOrderState{Tx, Tf}
    x::Tx
    f::Tf
    ∇f::Vector{Tf}
end

state(::FirstOrder, f, g, h, x) = FirstOrderState(x, f(x), g(x))

struct SecondOrderState{Tx, Tf}
    x::Tx
    f::Tf
    ∇f::Vector{Tf}
    ∇²f::Matrix{Tf}
end

state(::SecondOrder, f, g, h, x) = SecondOrderState(x, f(x), g(x), h(x))


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(AcuteML.render2file)),NamedTuple{(:id, :age, :field, :GPA, :courses),Tuple{Int64,Int64,String,Float64,Array{String,1}}},typeof(render2file),String,Bool})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Array{Union{Expr, Symbol},1},Array{Any,1},Array{Union{Expr, Symbol, Type},1},Array{Union{Expr, Symbol},1},Array{Union{Missing, String},1},Array{Union{Missing, Function, Symbol},1},Array{Union{Missing, Type},1},String,Type,Array{Union{Missing, Function, Symbol},0},Bool,Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Expr})
    Base.precompile(Tuple{typeof(AcuteML.aml_create),Expr,Array{Union{Expr, Symbol},1},Array{Any,1},Array{Union{Expr, Symbol, Type},1},Array{Union{Expr, Symbol},1},Array{Union{Missing, String},1},Array{Union{Missing, Function, Symbol},1},Array{Union{Missing, Type},1},String,Type,Array{Union{Missing, Function, Symbol},0},Bool,Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Array{Union{Nothing, Expr},1},Symbol})
    Base.precompile(Tuple{typeof(AcuteML.aml_parse),Expr})
    Base.precompile(Tuple{typeof(AcuteML.get_struct_xmlcreator),Expr,Array{Any,1},Array{Union{Expr, Symbol},1},Nothing,Expr,Array{Expr,1},Array{Expr,1},Nothing})
    Base.precompile(Tuple{typeof(AcuteML.multiString),Array{String,1}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Array{Float64,1},Type{AbsNormal}})
    Base.precompile(Tuple{typeof(addnode!),Node,String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Float64,1}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{Int64,1}},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(findcontent),Type{Array{String,1}},String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(get),Node,String,Nothing})
    Base.precompile(Tuple{typeof(updatecontent!),Array{Any,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{String,1},String,Node,Type{AbsNormal}})
    Base.precompile(Tuple{typeof(updatecontent!),Array{String,1},String,Node,Type{AbsText}})
    Base.precompile(Tuple{typeof(updatecontent!),Dict{String,String},String,Node,Type{AbsNormal}})
end

max_stars_repo_name: aminya/AML

max_forks_repo_name: aminya/AML


text:
@testset "456.132-pattern.jl" begin
    @test find132pattern([3, 1, 4, 2]) == true
    @test find132pattern([1, 2, 3, 4]) == false
    @test find132pattern([-1, 3, 2, 0]) == true
end


text:
module ConvexHull
using LinearAlgebra
using StaticArrays

import Polyhedra # for convex hull
import GLPK

module LibTriangle
	using Triangle
end
include("Projectors.jl")
using .Projectors
# using MiniQhull

# tools««1
# findextrema««2
"""
    findextrema(itr; lt=isless)

Like `findmin`, except that
 - it returns both extrema, as a `NamedTuple`;
 - it accepts an `lt` parameter, like `sort`.
"""
function findextrema(itr; lt=isless)
  p = pairs(itr); y = iterate(p)
  if y == nothing
    throw(ArgumentError("collection must be non-empty"))
  end
  (mi, m), s = y; (Mi, M) = (mi, m); i = mi
  while true
    y = iterate(p, s)
    y == nothing && break
    (ai, a), s = y
    if lt(a, m) m = a; mi = ai; end
    if lt(M, a) M = a; Mi = ai; end
  end
  return (min=(m, mi), max=(M, Mi))
end
# determinants««2
@inline det2(x,y) = x[1]*y[2]-x[2]*y[1]
@inline det2(x,y,z) = det2(y-x, z-x)
# triangulations««2
"""
    triangulate_face(points; direction, map)

Returns a triangulation of the face (assumed convex; points in any order)
Optional keyword arguments:
 - `direction` is a normal vector (used for projecting to 2d).
 - `map` is a labeling of points (default is identity map).
 - `convex` is a Val(true) or Val(false).

The triangulation is returned as a vector of StaticVector{3,Int},
containing the labels of three points of each triangle.
"""
function triangulate_face(points::AbstractVector{<:StaticVector{3}};
		direction::AbstractVector = face_normal(points),
		map::AbstractVector{<:Integer} = [1:length(points)...],
		)
	axis = main_axis(direction)
	N = length(points)
	# this common case deserves a shortcut:
	N == 3 && return [(map[1], map[2], map[3])]

	points2d = Matrix{Float64}(undef, N, 2)
	for (i, p) in pairs(points)
		points2d[i,:] .= project2d(axis, p)
	end
	tri = LibTriangle.basic_triangulation(points2d, map)
	return [(t[1], t[2], t[3]) for t in tri]
end
# Polyhedra interface««1
@inline polyhedra_lib(T::Type{<:Real}) =
	Polyhedra.DefaultLibrary{T}(GLPK.Optimizer)

# fixing an oversight in Polyhedra.jl: it has multiplications but no
# divisions
@inline Base.:(/)(h::Polyhedra.HyperPlane, α::Real) =
	Polyhedra.HyperPlane(h.a / α, h.β / α)

# HRepElement is the supertype of HalfSpace and HyperPlane
@inline direction(h::Polyhedra.HRepElement) = h.a
@inline function normalize(h::Polyhedra.HRepElement)
	n = norm(direction(h))
	(n ≠ 0) ? (h / n) : h
end
@inline (h::Polyhedra.HRepElement)(p) = dot(h.a, p) - h.β
@inline ∈(p, h::Polyhedra.HyperPlane) = iszero(h(p))
@inline Base.convert(T::Type{<:Polyhedra.HRepElement},
		h::Polyhedra.HRepElement) = T(h.a, h.β)

# 2d convex hull ««1

"""
    convex_hull_list(points)

Returns the convex hull of the points, as a list of indexes (in direct
order, starting at a reproducible index in the list of points).
"""
function convex_hull_list(points)
  # Uses Graham scan
  # In practice this is faster than using `Polyhedra.jl`.
#   println("points=$points, length=$(length(points))")
  i0 = findextrema(points;
    lt=(p,q)->(p[2]<q[2])|| (p[2]==q[2] && p[1]>q[1])).min[2]
  @inline detp2(i,j,k) = det2(points[[i,j,k]]...)
	# 1024 is an arbitrary value for detecting “aligned” points (i.e. up to
	# representation errors), which should be fast for both Float and Fixed
	# division
  @inline function are_aligned(i,j,k)
    v1 = points[j]-points[i]
    v2 = points[k]-points[j]
    d = det2(v1, v2)
    c = dot(v1, v2)
    return abs(d) < abs(c)/1024
  end
  scan = sort(filter(!isequal(i0), eachindex(points)),
    lt=(p,q)->detp2(i0,p,q) > 0)
#   println("i0=$i0, scan=$scan")
  stack = [i0, scan[1]]
  for h in scan[2:end]
#     println("scanning: $stack + $h")
    v1 = points[stack[end]] - points[stack[end-1]]
    v2 = points[h] - points[stack[end]]
    s = det2(v1, v2)
    c = dot(v1, v2)
    if abs(s) < abs(c)/1024 && c < 0 # points are aligned and backwards
			# here we know that we can insert at (end)
			# look for an insertion point i:
			i = length(stack)
			while i > 2
# 				println(" try to insert at $i")
				v1 = points[stack[i]] - points[stack[i-1]]
				v2 = points[h] - points[stack[i]]
				s = det2(v1, v2)
				c = dot(v1, v2)
				if s < -1e-3*abs(c)
# 					println(" break at $i")
					break
				end
				i -= 1
			end
# 			println("  inserting at $i")
			insert!(stack, i, h)
			continue
# 			println("  now stack=$stack")
    end
    while detp2(last(stack), h, stack[end-1]) < 0
      pop!(stack)
    end
    push!(stack, h)
  end
  return stack
end

"""
    convex_hull([vector of 2d points])

Returns the convex hull (as a vector of 2d points, ordered in direct
order).
"""
function convex_hull(points::AbstractVector{<:StaticVector{2}})
	upoints = unique(sort(points))
	return upoints[convex_hull_list(upoints)]
end

# """
#     convex_hull(x::Geometry{3}...)
#
# Returns the convex hull of the union of all the given solids, as a
# pair `(points, faces)`. `faces` is a list of triangles.
# """
# @inline convex_hull(x::AbstractGeometry{3}) =
# 	convex_hull(vcat([vertices(y) for y in x]...))

# this is the version using MiniQhull: #««
# convex_hull(points::AbstractVector{<:AnyVec(2)}) =
#		# Delaunay triangulation:
#		let T = MiniQhull.delaunay([p[i] for i in 1:2, p in points]),
#				N = length(points),
#				M = zeros(Bool, N, N) # ⚠ memory O(N²)
#		# mark all edges:
#		for (a,b,c) in eachcol(T)
#			b1 = points[b] - points[a]
#			c1 = points[c] - points[a]
#			d = b1[1]*c1[2] - b1[2]*c1[1] # determinant === orientation of triangle
#			if d < 0
#				M[a,b] = M[b,c] = M[c,a] = true
#			else
#				M[b,a] = M[c,b] = M[a,c] = true
#			end
#		end
#		# list of remaining edges (retrograde oriented)
#		L= sort([(i,j) for i in 1:N, j in 1:N if M[i,j] && ! M[j,i]], by=v->v[1])
#		next(i) = L[searchsorted(L, i, by=y->y[1])][1][2]
#		R = zeros(Int, length(L))
#		R[1:2] .= L[1] # initialize with first edge
#		for i in 3:length(R)
#			R[i] = next(R[i-1])
#		end
#		# returns in retrograde ordering (OpenSCAD convention):
#		points[R]
# end#»»
# this is the version using Polyhedra:««
# """
#     convex_hull([vector of 2d points])
#
# Returns the convex hull (as a vector of 2d points).
# """
# function convex_hull(points::Union{AbstractVector{<:Vec{2}},
# 	AbstractMatrix{<:Number}})
# # M is a matrix with the points as *columns*, hence the transpose
# # below:
# 	PH = Polyhedra
# 	poly = vpoly(points)
# 	PH.removevredundancy!(poly)
# 	return Vec{2}.(collect(PH.points(poly)))
# end
# »»

# 3d convex hull««1
"""
    convex_hull(vector of 3d points)

Returns the convex hull of these points, as a pair `(points, faces)`.
All the faces are triangles.
"""
function convex_hull(p::AbstractVector{<:StaticVector{3,T}}) where{T}
	M = reduce(hcat, p)
	PH = Polyhedra
	poly = PH.polyhedron(PH.vrep(transpose(M)), polyhedra_lib(T))
	R = PH.removevredundancy!(poly)
	V = SVector{3,T}.(collect(PH.points(poly)))

	triangles = SVector{3,Int}[]
	for i in PH.eachindex(PH.halfspaces(poly)) # index of halfspace
		h = PH.get(poly, i)
		pts = PH.incidentpointindices(poly, i) # vector of indices of points
		vlist = [SVector{3,T}(PH.get(poly, j)) for j in pts]
		for t in triangulate_face( vlist;
				direction = h.a,
				map = [j.value for j in pts])
			(a,b,c) = (V[j] for j in t)
			k = det([b-a c-a h.a])
			push!(triangles, (k > 0) ? t : (t[1], t[3], t[2]))
		end
	end
	return (points=V, faces=triangles)
end

#»»1
export convex_hull
end # module


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

## getindex
const rep = 20
const reps = rep^2

function sparse_getindex_perf()
    ## create some sparse matrices (need to be square):
    seed = 1
    srand(seed)

    small = 10^3
    med = 10^4
    large = 10^5
    huge = 10^6
    # # 1 entry per line
    # ss = []
    # push!(ss, sprand(small, small, 1e-3))
    # #push!(ss, sprand(med, med, 1e-4))
    # push!(ss, sprand(large, large, 1e-5))

    # 10 entries per line
    ts = []
    push!(ts, sprand(small, small, 1e-2))
    #push!(ts, sprand(med, med, 1e-3))
    push!(ts, sprand(large, large, 1e-4))

    # 100 entries per line
    us = []
    push!(us, sprand(small, small, 1e-1))
    #push!(us, sprand(med, med, 1e-2))
    push!(us, sprand(large, large, 1e-3))
    #push!(us, sprand(huge, huge, 1e-4))

    # # 1000 entries per line
    # vs = []
    # push!(vs, sprand(small, small, 1.0))
    # #push!(vs, sprand(med, med, 1e-1))
    # push!(vs, sprand(large, large, 1e-2))
    # #push!(vs, sprand(huge, huge, 1e-3))

    ## using uint32 (works up to 10^9)
    uus = []
    for u in us
        push!(uus, SparseMatrixCSC(u.m, u.n, map(UInt32,u.colptr), map(UInt32,u.rowval), u.nzval))
    end

    # test performance with matrices in us, uus and ts for
    # - integer indexing
    # - range indexing
    # - ordered array indexing
    # - disordered array indexing
    # - 1d indexing with integers and booleans

    # setup:
    # - indices
    intinds = nothing
    logicalinds = nothing # needs to be generated for a specific matrix size.
    rangeinds = 121:237
    orderedinds = collect(rangeinds)
    disorderedinds = orderedinds[randperm(length(orderedinds))]

    inds = [(intinds, "integers"), (logicalinds, "logical array"), (rangeinds, "a range"),
            (orderedinds, "a ordered array"), (disorderedinds, "a disordered array")]

    # - matrix sizes
    sizes = [(1, "small", "Small sparse matrix"), (2, "medium", "Medium sparse matrix")]

    # - matrix types
    mattyp = [(ts, "10 entries/column"), (us, "100 entries/column"), (uus, "100 entries/column uint32")]

    # - functions
    funs = [(integer_indexing, 1, "indexing"), (one_arg_indexing, 1, "1d indexing"),
            (row_indexing, 2, "indexing rows"), (col_indexing, 2, "indexing columns"),
            (row_col_indexing, 3, "indexing rows & columns")]

    # performance tests:
    counters = [1,1] # for small and medium matrix
    # integer indexing
    for (sz,s1,s2) in sizes # size of the matrix
        for (mt, ms) in mattyp  # type of the matrix
            m = mt[sz]
            c = counters[sz]
            @timeit integer_indexing(m) "sparse_getindex_$s1$c" "$s2 with $ms, indexing with integers"
            counters[sz] += 1
        end
    end

    # range & array indexing
    for (sz,s1,s2) in sizes # size of the matrix
        for (mt, ms) in mattyp  # type of the matrix
            m = mt[sz]
            for (fun, nargs, funstr) in funs[3:5] # indexing test function
                for (ind,indstr) in inds[3:5] # type of indices
                    c = counters[sz]
                    if indstr=="logical array"
                        # make a logical array of the right size
                        ind = sprand(Bool, size(m,1)..., 1e-5)
                    end
                    if nargs==2
                        @timeit fun(m, ind)  "sparse_getindex_$s1$c" "Sparse matrix with $ms, $funstr with $indstr"
                    elseif nargs==3
                        @timeit fun(m, ind, ind) "sparse_getindex_$s1$c" "Sparse matrix with $ms, $funstr with $indstr"
                    else
                        error("Something is amiss here.")
                    end
                    counters[sz] += 1
                end
            end
        end
    end

    # linear indexing
    for (sz,s1,s2) in sizes # size of the matrix
        for (mt, ms) in mattyp  # type of the matrix
            m = mt[sz]
            for (ind,indstr) in inds[2:5] # type of indices
                if indstr=="logical array"
                    if sz==2
                        continue # logical indexing with medium size sparse matrix takes too long
                    end
                    # make a logical array of the right size
                    ind = sprand(Bool, size(m)..., 1e-5)
                    c = counters[sz]
                    @timeit one_arg_indexing(m, ind) "sparse_getindex_$s1$c" "$s2 with $ms, linear indexing with $indstr"
                    counters[sz] += 1
                else
                    c = counters[sz]
                    @timeit one_arg_indexing(m, ind) "sparse_getindex_$s1$c" "$s2 with $ms, linear indexing with $indstr"
                    counters[sz] += 1
                end
            end
        end
    end

end

function integer_indexing(A)
    # index with two random integers
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for i in rand(rI, reps)
        for j in rand(rJ, rep)
            tmp += A[i,j]
        end
    end
    tmp
end

function row_indexing(A, rowinds)
    # index rows with rowinds and columns with a random integer
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for j in rand(rJ, reps)
        tmp += sum(A[rowinds,j])
    end
    tmp
end

function col_indexing(A, colinds)
    # index rows with a random integer and columns with colinds
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    tmp = zero(eltype(A))
    for i in rand(rI, div(reps,10) )
        tmp += sum(A[i,colinds])
    end
    tmp
end

function row_col_indexing(A, rowinds, colinds)
    # index rows with rowinds and columns with colinds
    # we need:
    (maximum(rowinds)+rep < size(A,1) && maximum(colinds)+rep < size(A, 2)) || error("bad rowinds or colinds")
    nI, nJ = size(A)
    rI = 1:nI
    rJ = 1:nJ
    for i in 1:10
        for j in 1:10
            tmp2 = A[rowinds.+i, colinds.+j]
        end
    end
end

function one_arg_indexing(A, lininds)
    # This is for 1d-indexing and indexing with one array of logicals.
    # Both return a nx1 sparse matrix.
    tmp = zero(eltype(A))
    if isa(eltype(A), Bool)
        tmp = sum(A[lininds])
    else
        for i in 1:rep
            tmp += sum(A[lininds])
        end
    end
    tmp
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia