text:
randspin()                               = [1,-1][rand(1:2)] # Generate a random spin
spingrid(n::Int)                         = [randspin() for i in 1:n, j in 1:n] # Generate a random spin array
magnetization(a::Array{Int, 2})          = mean(a) |> abs # Get magnetizations of the grid
namefunc(f::Function)                    = "$f"[1:end-1]  # Get the simplified name of a function
nspins(a::Array{Int, 2}, i::Int, j::Int) = [a[x,y] for (x,y) in neighbors(a,i,j)] # Get surrounding spins

# Flip a spin
function flip!(grid::Array{Int, 2}, cluster::BitArray{2})
    for j in 1:size(grid, 2), i in 1:size(grid, 1)
        if cluster[i, j] == true grid[i, j] *= -1 end
    end
end

# Return the sum of spins on a cluster
function clusterspin(grid::Array{Int, 2}, cluster::BitArray{2})
    spin = 0
    for j in 1:size(grid, 2), i in 1:size(grid, 1)
        if cluster[i,j] == true spin += grid[i,j] end
    end
    return spin
end

# Return an array of neighbors coordinates (each one in a tuple of the kind (x,y))
function neighbors(grid::Array{Int, 2}, i::Integer, j::Int)
    n = Array((Int,Int), 0)

    if i > 1             push!(n, (i-1, j  )) end
    if j > 1             push!(n, (i,   j-1)) end
    if i < size(grid, 1) push!(n, (i+1, j  )) end
    if j < size(grid, 2) push!(n, (i,   j+1)) end

    return n
end

# Phase diagram (magnetization by temperature) using given algorithm
function diagram(func::Function;
                 size::Integer      = 10,    # Size of the grid
                 ensembles::Integer = 50,    # Number of ensembles
                 h::Float64         = 0.0,   # External field
                 mintemp::Float64   = 0.5,   # Starting temperature
                 step::Float64      = 0.2,   # Step of temperatures
                 maxtemp::Float64   = 6.0,   # Final temperature
                 iters::Integer     = 50000, # Number of the iterations
                 plot::Bool         = true,  # Plot flag
                 verbose::Bool      = true)  # Verbose flag

    name  = namefunc(func)
    temps = Float64[]
    mags  = Float64[]
    for t in mintemp:step:maxtemp
        m = mean([func(spingrid(size), h=h, temp=t, iters=iters, plot=false, verbose=false)[end] for i in 1:ensembles])
        if verbose println("(T=$t) Avg. magnetization after $name: $m") end

        push!(temps, t)
        push!(mags,  m)
    end

    if plot
        PyPlot.plot(temps, mags, "o", color="blue")
        PyPlot.plot(temps, mags, "-", color="blue")
        PyPlot.title("Phase Diagram ($name for size $size and h=$h)")
        PyPlot.xlabel("Temperature")
        PyPlot.ylabel("Mean Magnetization")
        PyPlot.ylim(0,1.1)
        PyPlot.savefig("diag_$(name)_H=$h.png")
        PyPlot.close()
    end

    return temps, mags
end


text:
ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

function timeit()
    ticks = 0
    try
        while true
            sleep(0.5)
            ticks += 1
            println(ticks)
        end
    catch
    end
end

@time timeit()
println("Done.")

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestCount

using Test
using LazyGroupBy

@testset begin
    @test count.(_ -> true, grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) ==
          Dict(false => 3, true => 7)
    @test count.(grouped(identity, [true, false, true, true, true])) ==
          Dict(false => 0, true => 4)
    @test_throws(
        Exception,
        count.(_ -> 0.5, grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])),
    )
end

end  # module


text:
"""
    test_ClustData(t::ClustData,ref::ClustData)

Test if the two structs ClustData are identical in key properties
"""
function test_ClustData(t::ClustData,ref::ClustData)
    @test t.region == ref.region
    @test t.K == ref.K
    @test t.T == ref.T
    @test all(t.weights .≈ ref.weights)
    @test all(t.delta_t .≈ ref.delta_t)
    @test all(t.k_ids .≈ ref.k_ids)
    for (k,v) in t.data
        @test all(t.data[k] .≈ ref.data[k])
        @test all(t.mean[k] .≈ ref.mean[k])
        @test all(t.sdv[k] .≈ ref.sdv[k])
    end
end

"""
    test_ClustResult

Tests if the two structs ClustResult are identical in key properties
"""
function test_ClustResult(t::AbstractClustResult,ref::AbstractClustResult)
    @test t.cost ≈ ref.cost
    test_ClustData(t.clust_data,ref.clust_data)
end


text:
function plotAffProp(Dict_AffProp, ; chosen=dists, ymin = 0.001, anotar = 0)
    L = 1
    figArray = Any[]
    Kolors = palette(:tab10, length(chosen))
    Xmax = maximum(collect(Dict_AffProp[d*"_TimeVec"][end] for d in chosen))
    # Xmax = maximum(collect(Dict_AffProp[d*"_memInf_TimeVec"][end] for d in dists))

    for k in 1:3
        fig = plot()
        # Legends box
        if k ==4
            fig = hline!([0.001], line = L, label = "0.001 %", legend = :outertop)
            fig = hline!([0.01], line = L, label = "0.01 %")
            fig = hline!([0.1], line = L, label = "0.1 %")
            fig = hline!([1], line = L, label = "1 %")
            fig= scatter!([0],[0], label = "Memory: 0", m = :circle, markersize = 8, color = :white)
            fig= scatter!([0],[0], label = "Memory: 1", m = :rect, markersize = 8, color = :white)
            fig= scatter!([0],[0], label = "Memory: Inf", m = :utriangle, markersize = 8, color = :white)
            else
                fig = hline!([0.001], line = L, label = "")
                fig = hline!([0.01], line = L, label = "")
                fig = hline!([0.1], line = L, label = "")
                fig = hline!([1], line = L, label = "")

            end

        mem = ["0" "1" "Inf"][k]
        for i in 1:length(chosen)
            phont = "Times New Roman"
            phontsize = 14
            d = chosen[i]

            if d == "GAPMVec"
                marquer = :diamond
                elseif (d != "GAPMVec")&&(!issubset([d], dists))
                        marquer = :star
                end

            kol = Kolors[i]
            d0 = d

            if issubset([d0], dists) # Check if its a distance
                d = d0*"_mem$(mem)"
                Gap = Dict_AffProp[d*"_Gap"][end]
                Time = Dict_AffProp[d*"_TimeVec"][end]
                # Markers
                if mem == "0"
                    marquer = :circle
                    elseif mem == "1"
                        marquer = :rect
                    elseif mem == "Inf"
                        marquer = :utriangle
                end

                if (Gap > 0.1)&&(Gap<1)
                    fig = scatter!([Time], [Gap], m = marquer, xlabel = "Time [s]", ylabel = "Gap [%]", yaxis = :log, label = "", ylim=[ymin,100], c = kol, markerstrokewidth=1, markersize = 8)
                    elseif (Gap < 0.1)&&(Gap>10^-6)
                        fig= scatter!([Time], [Gap], m = marquer, xlabel = "Time [s]", ylabel = "Gap [%]", yaxis = :log, label = "", ylim=[ymin,100], c = kol, markerstrokewidth=1, markersize = 8)
                        if (d0 == "jaccard")&&(anotar==1)
                        fig=annotate!((Time, 2*Gap, (d0, :top, phontsize, phont)))
                        elseif (d0 == "whammingvar")&&(anotar==1)
                            fig=annotate!((Time-17, 0.75*Gap, (d0, :left, phontsize, phont)))
                        elseif (anotar==1)
                            fig=annotate!((Time+2, Gap, (d0, :left, phontsize, phont)))
                        end
                    elseif Gap > 1
                        fig= scatter!([Time], [Gap], m = marquer, xlabel = "Time [s]", ylabel = "Gap [%]", yaxis = :log, label = "", ylim=[ymin,100], c = kol, markerstrokewidth=1, markersize = 8)
                            d0 == "kl_divergence" ? dlabel = "kl_div" : nothing
                            d0 == "gkl_divergence" ? dlabel = "gkl_div" : nothing
                            if Time > 20
                                if (d0 == "gkl_divergence")&&(mem=="0")&&(anotar==1)
                                        fig=annotate!((Time-2, 2*Gap, (dlabel, :bottom, phontsize, phont)))
                                    end
                                else
                                    if (d0 == "kl_divergence")&&(mem=="0")&&(anotar==1)
                                        fig=annotate!((Time+10, 1.1*Gap, (dlabel, :left, phontsize, phont)))
                                    end
                                end
                end

            end
            if !(issubset([d0], dists))&&(k==1)
                println()
                println("d: $d")
                Gap = Dict_AffProp[d*"_Gap"][end]
                Time = Dict_AffProp[d*"_TimeVec"][end]
                fig = scatter!([Time], [Gap], m = marquer, xlabel = "Time [s]", ylabel = "Gap [%]", yaxis = :log, label = "", ylim=[ymin,100],  c = kol, markerstrokewidth=1, markersize = 8)
                anotar == 1 ? fig=annotate!((Time+2, Gap, (d0 == "GAPMVec" ? d0 : d0[6:end], :left, phontsize, phont))) : nothing

            end
        end
        fig = xlims!(-5, 1.05*Xmax)
        fig = plot!(size=(1000*upscale, 600*upscale), bottom_margin = -upscale*6mm, top_margin = upscale*0mm, left_margin = -12.5mm, right_margin = -1mm, titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm, guidefontsize=fntsm)
            push!(figArray,fig)

            display(fig)
            # display(fig)
            savefig(fig, figsDir*"/AffProp_$mem.png")
            savefig(fig, figsDir*"/AffProp_$mem.svg")
    end

    # fig1 = plot(figArray..., layout= grid(3,1))#, widths=[0.3 ,0.3, 0.4]))
    end
    plotAffProp(Dict_AffProp, chosen = DictAffProp, ymin= 5*10^(-6))

    # plotAffProp(chosen = DictAffProp, ymin= 0.001)

y = rand(10)
    plot(y)
dists[1][1]=="b"
indexin(dists[1],dists)

filter( x -> x == DictAffProp[1], dists)
issubset(["GAPM"], dists)
maximum(collect(Dict_AffProp[d*"_TimeVec"][end] for d in DictAffProp))

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
println("What is your name?")
name = readline()
println("Hello $name") #string interpolation sytax

max_stars_repo_name: j-hayes/intro-to-computing-sci-eng

max_issues_repo_name: j-hayes/intro-to-computing-sci-eng

max_forks_repo_name: j-hayes/intro-to-computing-sci-eng


text:
#!/usr/bin/env julia --project=.

using PowerModelsSecurityConstrained
pmsc_path = joinpath(dirname(pathof(PowerModelsSecurityConstrained)), "..")

InFile1=joinpath(pmsc_path, "test/data/scenario_01/case.con")
InFile2=joinpath(pmsc_path, "test/data/scenario_01/case.inl")
InFile3=joinpath(pmsc_path, "test/data/scenario_01/case.raw")
InFile4=joinpath(pmsc_path, "test/data/scenario_01/case.rop")

TimeLimitInSeconds=600
ScoringMethod=2
NetworkModel="IEEE 14"

include("MyJulia1.jl")
MyJulia1(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds, ScoringMethod, NetworkModel)

include("MyJulia2.jl")
MyJulia2(InFile1, InFile2, InFile3, InFile4, TimeLimitInSeconds, ScoringMethod, NetworkModel)


text:
#########################Initial Condition Updating#########################################
# TODO: Consider when more than one UC model is used for the stages that the counts need
# to be scaled.
function calculate_ic_quantity(
    initial_condition_key::ICKey{TimeDurationON, T},
    ic::InitialCondition,
    var_value::Float64,
    cache::TimeStatusChange,
) where {T <: PSY.Component}
    name = device_name(ic)
    time_cache = cache_value(cache, name)

    current_counter = time_cache[:count]
    last_status = time_cache[:status]
    var_status = isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
    @debug last_status, var_status, abs(last_status - var_status)
    @assert abs(last_status - var_status) < ABSOLUTE_TOLERANCE

    return last_status >= 1.0 ? current_counter : 0.0
end

function calculate_ic_quantity(
    initial_condition_key::ICKey{TimeDurationOFF, T},
    ic::InitialCondition,
    var_value::Float64,
    cache::TimeStatusChange,
) where {T <: PSY.Component}
    name = device_name(ic)
    time_cache = cache_value(cache, name)

    current_counter = time_cache[:count]
    last_status = time_cache[:status]
    var_status = isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
    @debug last_status, var_status, abs(last_status - var_status)
    @assert abs(last_status - var_status) < ABSOLUTE_TOLERANCE

    return last_status >= 1.0 ? 0.0 : current_counter
end

function calculate_ic_quantity(
    initial_condition_key::ICKey{DeviceStatus, T},
    ic::InitialCondition,
    var_value::Float64,
    cache::Union{Nothing, AbstractCache},
) where {T <: PSY.Component}
    return isapprox(var_value, 0.0, atol = ABSOLUTE_TOLERANCE) ? 0.0 : 1.0
end

function calculate_ic_quantity(
    initial_condition_key::ICKey{DevicePower, T},
    ic::InitialCondition,
    var_value::Float64,
    cache::Union{Nothing, AbstractCache},
) where {T <: PSY.ThermalGen}
    status_change_to_on =
        get_condition(ic) <= ABSOLUTE_TOLERANCE && var_value >= ABSOLUTE_TOLERANCE
    status_change_to_off =
        get_condition(ic) >= ABSOLUTE_TOLERANCE && var_value <= ABSOLUTE_TOLERANCE
    if status_change_to_on
        return ic.device.activepowerlimits.min
    end

    if status_change_to_off
        return 0.0
    end

    return var_value
end

function calculate_ic_quantity(
    initial_condition_key::ICKey{EnergyLevel, T},
    ic::InitialCondition,
    var_value::Float64,
    cache::Union{Nothing, AbstractCache},
) where {T <: PSY.Component}

    name = device_name(ic)
    energy_cache = cache_value(cache, name)
    if energy_cache != var_value
        return var_value
    end
    return energy_cache
end

"""
Status Init is always calculated based on the Power Output of the device
This is to make it easier to calculate when the previous model doesn't
contain binaries. For instance, looking back on an ED model to find the
IC of the UC model
"""
function status_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    _make_initial_conditions!(
        container,
        devices,
        ICKey(DeviceStatus, T),
        _make_initial_condition_active_power,
        _get_status_value,
    )

    return
end

function output_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    _make_initial_conditions!(
        container,
        devices,
        ICKey(DevicePower, T),
        _make_initial_condition_active_power,
        _get_active_power_output_value,
    )

    return
end

function duration_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    for key in (ICKey(TimeDurationON, T), ICKey(TimeDurationOFF, T))
        _make_initial_conditions!(
            container,
            devices,
            key,
            _make_initial_condition_active_power,
            _get_duration_value,
            TimeStatusChange,
        )
    end

    return
end

######################### Initialize Functions for Storage #################################
# TODO: This IC needs a cache for Simulation over long periods of tim
function storage_energy_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.Storage}
    key = ICKey(EnergyLevel, T)
    _make_initial_conditions!(
        container,
        devices,
        key,
        _make_initial_condition_energy,
        _get_energy_value,
    )

    return
end

######################### Initialize Functions for Hydro #################################
function status_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    _make_initial_conditions!(
        container,
        devices,
        ICKey(DeviceStatus, T),
        _make_initial_condition_active_power,
        _get_status_value,
        # Doesn't require Cache
    )
end

function output_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    _make_initial_conditions!(
        container,
        devices,
        ICKey(DevicePower, T),
        _make_initial_condition_active_power,
        _get_active_power_output_value,
        # Doesn't require Cache
    )

    return
end

function duration_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    for key in (ICKey(TimeDurationON, T), ICKey(TimeDurationOFF, T))
        _make_initial_conditions!(
            container,
            devices,
            key,
            _make_initial_condition_active_power,
            _get_duration_value,
            TimeStatusChange,
        )
    end

    return
end

function storage_energy_init(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.HydroGen}
    key = ICKey(EnergyLevel, T)
    _make_initial_conditions!(
        container,
        devices,
        key,
        _make_initial_condition_reservoir_energy,
        _get_reservoir_energy_value,
        StoredEnergy,
    )

    return
end

function _make_initial_conditions!(
    container::InitialConditions,
    devices::IS.FlattenIteratorWrapper{T},
    key::ICKey,
    make_ic_func::Function,
    get_val_func::Function,
    cache = nothing,
) where {T <: PSY.Device}
    length_devices = length(devices)

    if !has_initial_conditions(container, key)
        @debug "Setting $(key.ic_type) initial conditions for the status of all devices $(T) based on system data"
        ini_conds = Vector{InitialCondition}(undef, length_devices)
        set_initial_conditions!(container, key, ini_conds)
        for (ix, dev) in enumerate(devices)
            val = get_val_func(dev, key)
            ic = make_ic_func(container, dev, val, cache)
            ini_conds[ix] = ic
            IS.@record :simulation InitialConditionUpdateEvent(
                Dates.DateTime(1970, 1, 1),
                key,
                ic,
                val,
                0,
            )
        end
    else
        ini_conds = get_initial_conditions(container, key)
        ic_devices = Set((IS.get_uuid(ic.device) for ic in ini_conds))
        for dev in devices
            IS.get_uuid(dev) in ic_devices && continue
            @debug "Setting $(key.ic_type) initial conditions for the status device $(PSY.get_name(dev)) based on system data"
            val = get_val_func(dev, key)
            ic = make_ic_func(container, dev, val, cache)
            push!(ini_conds, ic)
            IS.@record :simulation InitialConditionUpdateEvent(
                Dates.DateTime(1970, 1, 1),
                key,
                ic,
                val,
                0,
            )
        end
    end

    @assert length(ini_conds) == length_devices
    return
end

function _make_initial_condition_active_power(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_active_power(T, container), value, cache)
end

function _make_initial_condition_energy(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_energy(T, container), value, cache)
end

function _make_initial_condition_reservoir_energy(
    container,
    device::T,
    value,
    cache = nothing,
) where {T <: PSY.Component}
    return InitialCondition(device, _get_ref_reservoir_energy(T, container), value, cache)
end

function _get_status_value(device, key)
    return PSY.get_status(device) ? 1.0 : 0.0
end

function _get_active_power_output_value(device, key)
    return PSY.get_activepower(device)
end

function _get_energy_value(device, key)
    return PSY.get_energy(device)
end

function _get_reservoir_energy_value(device, key)
    return PSY.get_initial_storage(device)
end

function _get_duration_value(dev, key)
    if key.ic_type == TimeDurationON
        value = PSY.get_status(dev) > 0 ? PSY.get_time_at_status(dev) : 0.0
    else
        @assert key.ic_type == TimeDurationOFF
        value = PSY.get_status(dev) <= 0 ? PSY.get_time_at_status(dev) : 0.0
    end

    return value
end

function _get_ref_active_power(
    ::Type{T},
    container::InitialConditions,
) where {T <: PSY.Component}
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, ACTIVE_POWER) :
           UpdateRef{T}(ACTIVE_POWER, "get_activepower")
end

function _get_ref_energy(::Type{T}, container::InitialConditions) where {T <: PSY.Component}
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, ENERGY) :
           UpdateRef{T}(ENERGY, "get_energy")
end

function _get_ref_reservoir_energy(
    ::Type{T},
    container::InitialConditions,
) where {T <: PSY.Component}
    return get_use_parameters(container) ? UpdateRef{JuMP.VariableRef}(T, ENERGY) :
           UpdateRef{T}(ENERGY, "get_storage_capacity")
end


text:
include("SchmidtAlg.jl")
# global isInit= false
using Base.Threads
function doesWinding( CumSumDataobj::spaceCumSumData )
    # global isInit
    temobj = CumSumDataobj
    rectsumObj = temobj.rectSum
    cfgdata = temobj.CfgInfoObj
    chnid = temobj.startID >1 ? 2 : 1
    if CumSumDataobj.isInit == false
        # println("??" , )
        SpaceSegs1= []
        WindAlarmSegs1=[]
        id1=CumSumDataobj.startID
        id2=CumSumDataobj.endID
        init!(SpaceSegs1 ,cfgdata.meterperSeg,UInt16(id1) ,UInt16(id2))
        init!(WindAlarmSegs1,SpaceSegs1 )
        CumSumDataobj.WindAlarmSegs  = WindAlarmSegs1
        CumSumDataobj.SpaceSegs = SpaceSegs1
        CumSumDataobj.isInit  = true
    end
     WindAlarmSegs  = CumSumDataobj.WindAlarmSegs
     SpaceSegs = CumSumDataobj.SpaceSegs
    # curPos = 2
    # for idx = curPos:curPos #test only

    l1 = length( WindAlarmSegs )
    @threads for idx = 1:l1
        # println( idx ," ",length( WindAlarmSegs ))
        windingObject = WindAlarmSegs[idx]
        if chnid ==2
            offset = cfgdata.chn2SegBegIdx-1
        else
            offset =0
        end
        windStartSegID = windingObject.SegID1 - offset
        windEndsegID = windingObject.SegID2 - offset
        SegNum = windEndsegID - windStartSegID

        # println( chnid ,":" ,windStartSegID ,"---" , windEndsegID,"; "," $(cfgdata.chn2SegBegIdx)  ",length( rectsumObj ))
        # println("-----------------------------------------------------------------------------")
        @views a = rectsumObj[ windStartSegID : windEndsegID  ]
        # try
            FeatureAnal(a ,windingObject ,temobj.dt)
        # catch(e)

            # println(e)
        # end
    end

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text: include("PETSc.jl")


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

JLLWrappers.@generate_main_file_header("XNNPACK")
JLLWrappers.@generate_main_file("XNNPACK", UUID("55708af3-2096-517b-81a3-3917441130f6"))
end  # module XNNPACK_jll


text:
using Flux: BilinearUpsample
using Test

@testset "BilinearUpsample" begin
  @test size(BilinearUpsample((2, 2))(rand(2, 2, 1, 1))) == (4, 4, 1, 1)
  @test size(BilinearUpsample((3, 3))(rand(2, 2, 1, 1))) == (6, 6, 1, 1)
  @test size(BilinearUpsample((2, 2))(rand(2, 2, 10, 10))) == (4, 4, 10, 10)
  @test size(BilinearUpsample((3, 3))(rand(2, 2, 10, 10))) == (6, 6, 10, 10)

  @test_throws BoundsError BilinearUpsample((2, 2))(rand(2, 2))
end


text:
# ---
# title: 217. Contains Duplicate
# id: problem217
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/contains-duplicate/description/>
# hidden: true
# ---
#
# Given an array of integers, find if the array contains any duplicates.
#
# Your function should return true if any value appears at least twice in the
# array, and it should return false if every element is distinct.
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,1]
#     Output: true
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,4]
#     Output: false
#
# **Example 3:**
#
#
#
#     Input: [1,1,1,3,3,4,3,2,4,2]
#     Output: true
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
convert_to_int_vector(line) = map(x -> parse(Int64, x), split(line, ","))

function read_data(lines)
    function read_interval(s)
        a, b = map(x -> parse(Int64, x), split(s, "-"))
        (a=a, b=b)
    end

    ranges = Dict{String, Set{NamedTuple{(:a, :b), Tuple{Int64, Int64}}}}()
    ix1, ix2 = findall(isempty, lines)
    for line in lines[1:ix1-1]
        label, rest = split(line, ": ")
        intervals = split(rest, " or ")
        ranges[label] = Set(map(read_interval, intervals))
    end

    own_ticket = convert_to_int_vector(lines[ix1+2])

    tickets = Vector{Int64}[]
    for line in lines[ix2+2:end]
        push!(tickets, convert_to_int_vector(line))
    end

    (ranges, own_ticket, tickets)
end

item_in_intervals(item, intervals) = any(a <= item <= b for (a, b) in intervals)

is_valid(x, ranges) = any(item_in_intervals(x, intervals) for intervals in values(ranges))

function collect_suitable_labels(ranges, column)
    suitable_labels = []
    for (label, intervals) in pairs(ranges)
        if all(item_in_intervals(x, intervals) for x in column)
            push!(suitable_labels, label)
        end
    end
    suitable_labels
end

function find_unique_assignments(labels)
    result = Array{String}(undef, length(labels))
    ix = findfirst(lst -> length(lst) == 1, labels)
    while ix !== nothing
        label = only(labels[ix])
        result[ix] = label
        labels = [filter(x -> x != label, lst) for lst in labels]
        ix = findfirst(lst -> length(lst) == 1, labels)
    end
    result
end

function departure_product(own_ticket, tickets, ranges)
    is_ticket_valid(ticket) = all(is_valid(x, ranges) for x in ticket)

    valid_tickets = permutedims(hcat(filter(is_ticket_valid, tickets)...))
    label_assignments = Vector{String}[]
    for column in eachcol(valid_tickets)
        labels = collect_suitable_labels(ranges, column)
        push!(label_assignments, labels)
    end

    unique_assignments = find_unique_assignments(label_assignments)
    prod(number for (label, number) in zip(unique_assignments, own_ticket) if startswith(label, "departure"))
end

function main()
    lines = readlines("../resources/input_16.txt")
    ranges, own_ticket, nearby_tickets = read_data(lines)
    result = departure_product(own_ticket, nearby_tickets, ranges)
    println(result)
end

main()


text:
@testset "Parallel tests" begin

    model(x, θ) = θ[1]  + θ[2] * x[1]

    gmodel!(g, x, θ) = begin
        g[1] = 1.0
        g[2] = x[1]
    end

    n = 2

    np = 10

    p = 7

    Random.seed!(123456789)

    data, θSol, = RAFF.generate_noisy_data(model, n, np, p; std=0.0)

    # Remove outlier information
    data = data[:, 1:2]

    @testset "Updater" begin

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(0))

        bestθ = SharedArray{Float64, 1}(n)

        bestθ .= 0.0

        fut = @async RAFF.update_best(bqueue, bestθ)

        # Should update bestθ

        newbest1 = ones(Float64, n)

        put!(bqueue, newbest1)

        sleep(0.1)

        @test bestθ == newbest1

        # Should update bestθ again

        newbest2 = rand(Float64, n)

        put!(bqueue, newbest2)

        sleep(0.1)

        @test bestθ == (newbest1 + newbest2) / 2.0

        # Should not throw an error nor die with invalid element
        # The following test if failing in Julia Nightly. They do not
        # allow adding different types to specific-typed channels,
        # which is good. See how to fix this test.

        # put!(bqueue, 10)

        # @test !isready(bqueue)
        # @test !istaskdone(fut)

        # Should finish

        close(bqueue)

        sleep(0.1)

        @test istaskdone(fut)

    end

    @testset "Consumer" begin

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(0))

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        squeue = RemoteChannel(() -> Channel{RAFFOutput}(0))

        seedMS = MersenneTwister(1234)

        initguess = zeros(Float64, n)

        MAXMS = 1

        worker1 = @async @test begin
            RAFF.consume_tqueue(bqueue, tqueue, squeue,
                                model, gmodel!, data, n, p - 2, np,
                                MAXMS, seedMS, initguess, lmlovo)
            true
        end

        @test !istaskdone(worker1)

        # Should not do anything for an invalid interval
        put!(tqueue, p - 3:p)

        @test !isready(squeue)

        # Should work, since the problem is easy
        put!(tqueue, p - 2:p - 2)

        # Should return a vector with 1 solution pair
        rout = take!(squeue)

        @test !istaskdone(worker1)
        @test rout.p == p - 2
        @test rout.status == 1

        # Another test, with different p

        put!(tqueue, p:p)

        rout = take!(squeue)

        @test rout.p == p
        @test rout.status == 1
        @test rout.f ≈ 0.0 atol=1.0e-1
        @test rout.solution ≈ θSol atol=1.0e-1

        # Test with interval

        put!(tqueue, p - 1:p)

        rout = take!(squeue)

        @test rout.status == 1
        @test rout.p == p - 1

        rout = take!(squeue)

        @test rout.status == 1
        @test rout.p == p

        @test !isready(squeue)

        # Check if worker is alive when bqueue is closed

        close(bqueue)

        put!(tqueue, p:p)

        take!(squeue)

        @test !istaskdone(worker1)

        # Check if worker finishes when tqueue is closed

        close(tqueue)

        sleep(0.1)

        @test istaskdone(worker1)

        #  Worker should die if solution queue is closed

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        seedMS = MersenneTwister(1234)

        MAXMS = 1

        worker2 = @async @test begin
            RAFF.consume_tqueue(bqueue, tqueue, squeue,
                                model, gmodel!, data, n, p - 2, np,
                                MAXMS, seedMS, initguess, lmlovo)
            true
        end

        @test !istaskdone(worker2)

        close(squeue)

        put!(tqueue, p:p)

        sleep(0.1)

        @test !isready(tqueue)

        @test istaskdone(worker2)

    end

    @testset "Consumer Multistart" begin

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(4))

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        squeue = RemoteChannel(() -> Channel{RAFFOutput}(0))

        seedMS = MersenneTwister(1234)

        initguess = zeros(Float64, n)

        MAXMS = 1

        # Check if the worker dies after closing the task queue

        worker1 = @async @test begin
            RAFF.consume_tqueue(bqueue, tqueue, squeue,
                                model, gmodel!, data, n, p - 2, np,
                                MAXMS, seedMS, initguess, lmlovo)
            true
        end

        put!(tqueue, p:p)

        close(tqueue)

        rout1 = take!(squeue)

        sleep(0.1)

        @test istaskdone(worker1)

        # Save the objective function found and run multistart
        # strategy with the same initial random generator

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        MAXMS = 3

        seedMS = MersenneTwister(1234)

        worker = @async @test begin
            RAFF.consume_tqueue(bqueue, tqueue, squeue,
                                model, gmodel!, data, n, p - 2, np,
                                MAXMS, seedMS, initguess, lmlovo)
            true
        end

        put!(tqueue, p:p)

        rout2 = take!(squeue)

        close(tqueue)

        fetch(worker)

        close(bqueue)

        # Should find a better point

        @test rout1.f >= rout2.f
        @test rout1.p == rout2.p
        @test istaskdone(worker)

    end

    @testset "Worker checker" begin

        nworkers = 3

        # Test if all workers are dead

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(4))

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        squeue = RemoteChannel(() -> Channel{RAFFOutput}(0))

        futures = Vector{Future}(undef, nworkers)

        for i = 1:nworkers

            futures[i] = @spawn error()

        end

        RAFF.check_and_close(bqueue, tqueue, squeue, futures)

        @test !isopen(bqueue)
        @test !isopen(tqueue)
        @test !isopen(squeue)

        # Test if there is at least one live worker

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(4))

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        squeue = RemoteChannel(() -> Channel{RAFFOutput}(0))

        futures = Vector{Future}(undef, nworkers)

        for i = 1:nworkers - 1

            futures[i] = @spawn error()

        end

        futures[nworkers] = @spawn take!(tqueue)

        RAFF.check_and_close(bqueue, tqueue, squeue, futures)

        @test isopen(bqueue)
        @test isopen(tqueue)
        @test isopen(squeue)

        put!(tqueue, 1:1)

        RAFF.check_and_close(bqueue, tqueue, squeue, futures)

        @test !isopen(bqueue)
        @test !isopen(tqueue)
        @test !isopen(squeue)

        # Ensure that this checker does not closes the queue if the
        # workers have finished their job very fast

        bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(4))

        tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

        squeue = RemoteChannel(() -> Channel{RAFFOutput}(0))

        futures = Vector{Future}(undef, nworkers)

        for i = 1:nworkers

            futures[i] = @spawn(nothing)

        end

        # Simulates the case where all the tasks have already been
        # taken
        close(tqueue)

        RAFF.check_and_close(bqueue, tqueue, squeue, futures)

        @test isopen(bqueue)
        @test isopen(squeue)

    end

    @testset "PRAFF" begin

        model(x, θ) = θ[1] * exp(x[1] * θ[2])

        gmodel!(g, x, θ) = begin

            g[1] = exp(x[1] * θ[2])
            g[2] = x[1] * θ[1] * exp(x[1] * θ[2])

        end

        data = [-1.0   3.2974425414002564;
                -0.75  2.9099828292364025;
                -0.5    2.568050833375483;
                -0.25  2.2662969061336526;
                0.0                  2.0;
                0.25   1.764993805169191;
                0.5   1.5576015661428098;
                0.75  1.5745785575819442; #noise
                1.0   1.2130613194252668;
                1.25  1.0705228570379806;
                1.5   0.9447331054820294;
                1.75  0.8337240393570168;
                2.0   0.7357588823428847;
                2.25  0.6493049347166995;
                2.5   0.5730095937203802;
                2.75  0.5056791916094929;
                3.0  0.44626032029685964;
                3.25  0.5938233504083881; #noise
                3.5   0.3475478869008902;
                3.75 0.30670993368985694;
                4.0   0.5706705664732254; #noise
                ]

        answer = [2.0, -0.5]

        # Regular test

        rout = praff(model, data, 2)

        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.iter >= size(data)[1]
        @test rout.nf >= 1
        @test rout.nj >= 1

        rout = praff(model, gmodel!, data, 2)

        rgood = rout

        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.iter >= size(data)[1]
        @test rout.nf >= 1
        @test rout.nj >= 1

        # Multistart test

        rout = praff(model, data, 2; MAXMS=2)

        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18
        @test rout.f <= rgood.f
        @test rout.iter >= 1.1 * rgood.iter
        @test rout.nf >= 1.1 * rgood.nf
        @test rout.nj >= 1.1 * rgood.nj

    end

    @testset "Test parameters" begin

        data = [-1.0   3.2974425414002564;
                -0.75  2.9099828292364025;
                -0.5    2.568050833375483;
                -0.25  2.2662969061336526;
                 0.0                  2.0;
                 0.25   1.764993805169191;
                 0.5   1.5576015661428098;
                 0.75  1.5745785575819442; #noise
                 1.0   1.2130613194252668;
                 1.25  1.0705228570379806;
                 1.5   0.9447331054820294;
                 1.75  0.8337240393570168;
                 2.0   0.7357588823428847;
                 2.25  0.6493049347166995;
                 2.5   0.5730095937203802;
                 2.75  0.5056791916094929;
                 3.0  0.44626032029685964;
                 3.25  0.5938233504083881; #noise
                 3.5   0.3475478869008902;
                 3.75 0.30670993368985694;
                 4.0   0.5706705664732254; #noise
                ]

        answer = [2.0, -0.5]

        rout = praff(model, gmodel!, data, 2; noutliers=0)

        @test rout.p == 21

        rout = praff(model, data, 2; noutliers=5)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        rout = praff(model, data, 2; ftrusted=(21 - 5)/21)

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        rout = praff(model, data, 2; ftrusted=(18/21, 18/21))

        @test rout.f ≈ 0.0 atol=1.0e-5
        @test rout.solution ≈ answer atol=1.0e-5
        @test rout.p == 18

        @test praff(model, data, 2; ftrusted=(0.5, 1.1)) == RAFFOutput()
        @test praff(model, data, 2; ftrusted=-0.1) == RAFFOutput()

    end

end


text:
module MNIST

using GZip

function traindata(path::String)
    x = gzread(path,"train-images-idx3-ubyte.gz")[17:end]
    x = reshape(x/255, 28*28, 60000)
    x = Matrix{Float32}(x)
    y = gzread(path,"train-labels-idx1-ubyte.gz")[9:end]
    y = Vector{Int}(y)
    x, y
end

function testdata(path::String)
    data = gzread(path,"t10k-images-idx3-ubyte.gz")[17:end]
    x = reshape(data/255, 28*28, 10000)
    x = Matrix{Float32}(x)
    data = gzread(path,"t10k-labels-idx1-ubyte.gz")[9:end]
    y = Vector{Int}(data)
    x, y
end

function gzread(path::String, filename::String)
    mkpath(path)
    path = joinpath(path, filename)
    if !isfile(path)
        println("Downloading $filename...")
        download("http://yann.lecun.com/exdb/mnist/$filename", path)
    end
    gzopen(read, path)
end

end


text:
#push!(LOAD_PATH, "C:\\Users\\adamm\\Dropbox\\Development\\JuML\\src")
using JuML

importcsv("C:\\Users\\adamm_000\\Documents\\Julia\\airlinetrain.csv")

importcsv("C:\\Users\\adamm_000\\Documents\\Julia\\airlinetest.csv")

train_df = DataFrame("C:\\Users\\adamm_000\\Documents\\Julia\\airlinetrain") # note we are passing a path to a folder
test_df = DataFrame("C:\\Users\\adamm_000\\Documents\\Julia\\airlinetest")

factors = train_df.factors
covariates = train_df.covariates

distance = train_df["Distance"]
deptime = train_df["DepTime"]
dep_delayed_15min = train_df["dep_delayed_15min"]

summary(distance)
summary(deptime)
summary(dep_delayed_15min)

label = covariate(train_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)

deptime = factor(train_df["DepTime"], 1:2930)
distance = factor(train_df["Distance"], 11:4962)

factors = [train_df.factors; [deptime, distance]]

@time model = xgblogit(label, factors; η = 1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 1, maxdepth = 5, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);

pred = predict(model, test_df)
testlabel = covariate(test_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)
auc = getauc(pred, testlabel)
logloss = getlogloss(pred, testlabel)

@time cv = cvxgblogit(label, factors, 5; aucmetric = true, loglossmetric = true, trainmetric = true, η = 0.3, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 2, maxdepth = 5, caching = true, usefloat64 = false, singlethread = true);
cv
model.pred[1:5]
@time pred = predict(model, test_df)

pred[1:5]
testlabel = covariate(test_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)

@time auc = getauc(pred, testlabel)

#importcsv("src\\Data\\agaricus_train.csv"; isnumeric = (colname, _) -> colname == "label")

#importcsv("src\\Data\\agaricus_test.csv"; isnumeric = (colname, _) -> colname == "label")

#importcsv("src\\Data\\agaricus.csv")

traindf = DataFrame("src\\Data\\agaricus_train"; preload = true)
testdf = DataFrame("src\\Data\\agaricus_test"; preload = true)

df = DataFrame("src\\Data\\agaricus"; preload = true)

#class = df["class"]
#summary(class)

#label = JuML.CachedCovariate(JuML.ParseFactorCovariate("label", class, level -> level == "p" ? 1.0 : 0.0))
#factors = filter((x -> JuML.getname(x) != "class"), df.factors)

#@time trees, pred = JuML.xgblogit(label, factors; η = 1, nrounds = 1, maxdepth = 1, minchildweight = 0.0, nthreads = 1);

@time model = JuML.xgblogit(traindf["label"], traindf.factors; η = 1, nrounds = 1, maxdepth = 1, minchildweight = 0.0, singlethread = true, slicelength = 5000);

R:
require(xgboost)
data(agaricus.train, package='xgboost')
train <- agaricus.train
data(agaricus.test, package='xgboost')
test <- agaricus.test
model <- xgboost(data = train$data, label = train$label, eta = 1, nrounds = 1, max_depth = 1, objective = "binary:logistic")
preds = predict(model, train$data)

airline:
suppressMessages({
library(data.table)
library(ROCR)
library(xgboost)
library(Matrix)
})
set.seed(123)
d_train <- fread("C:\\Users\\adamm_000\\Documents\\Julia\\airlinetrain.csv", showProgress=FALSE, stringsAsFactors=TRUE)
X_train <- sparse.model.matrix(dep_delayed_15min ~ .-1, data = d_train)
dxgb_train <- xgb.DMatrix(data = X_train, label = ifelse(d_train$dep_delayed_15min=='Y',1,0))
system.time(md <- xgb.train(data = dxgb_train, objective = "binary:logistic", nround = 1, max_depth = 2, eta = 1, tree_method='exact'))
preds = predict(md, dxgb_train)
xgb.plot.tree(model = md, feature_names = colnames(dxgb_train))

system.time(md <- xgb.cv(data = dxgb_train, objective = "binary:logistic", nround = 1, max_depth = 10, eta = 1, tree_method='exact', nfold = 5, metrics= list("auc")))

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
##### Surface models

abstract type SurfaceModel end

struct SurfaceFixedFlux{FT} <: SurfaceModel
    T::FT
    P::FT
    q_tot::FT
    shf::FT
    lhf::FT
    Tsurface::FT
    ρ_0_surf::FT
    α_0_surf::FT
    ρq_tot_flux::FT
    ρθ_liq_flux::FT
    ustar::FT
    windspeed_min::FT
    tke_tol::FT
    area::FT
    function SurfaceFixedFlux{FT}(
        param_set;
        T,
        P,
        q_tot,
        ustar,
        windspeed_min,
        tke_tol,
        area,
    ) where {FT}
        q_pt = PhasePartition(q_tot)
        ρ_0_surf = air_density(param_set, T, P, q_pt)
        α_0_surf = 1 / ρ_0_surf
        Tsurface = 299.1 * TD.exner_given_pressure(param_set, P)
        lhf = 5.2e-5 * ρ_0_surf * latent_heat_vapor(param_set, Tsurface)
        shf = 8.0e-3 * cp_m(param_set, q_pt) * ρ_0_surf
        ρ_tflux = shf / cp_m(param_set, q_pt)
        ρq_tot_flux = lhf / latent_heat_vapor(param_set, Tsurface)
        ρθ_liq_flux = ρ_tflux / TD.exner_given_pressure(param_set, P)
        return new{FT}(
            T,
            P,
            q_tot,
            shf,
            lhf,
            Tsurface,
            ρ_0_surf,
            α_0_surf,
            ρq_tot_flux,
            ρθ_liq_flux,
            ustar,
            windspeed_min,
            tke_tol,
            area,
        )
    end
end


text:
export
# Solve API
    solve,
    normalize,
    discretize,
    homogenize,
    ensemble,

# Algorithms
    A20,
    ASB07,
    BFFPSV18,
    BOX,
    GLGM06,
    INT,
    LGG09,
    ORBIT,
    QINT,
    TMJets,
    TMJets20,
    TMJets21a,
    TMJets21b,
    VREP,

# Approximation models
    FirstOrderddt,
    FirstOrderZonotope,
    CorrectionHull,
    FirstOrder,
    ForwardBackward,
    #Backward,
    Forward,
    NoBloating,
    StepIntersect,

# Flowpipes
    flowpipe,
    Flowpipe,
    ShiftedFlowpipe,
    MappedFlowpipe,
    HybridFlowpipe,
    MixedFlowpipe,
    MixedHybridFlowpipe,

# Reach-sets
    ReachSet,
    SparseReachSet,
    ShiftedReachSet,
    TaylorModelReachSet,
    TemplateReachSet,

# Getter functions
    set,
    tstart,
    tend,
    tspan,
    vars,
    sup_func, # TODO keep?
    setrep,
    rsetrep,
    reset_map,
    guard,
    source_invariant,
    target_invariant,
    # getter functions for Taylor model reach-sets
    domain, remainder, polynomial, get_order, expansion_point,
    numrsets,

# Concrete operations
    project,
    shift,
    complement,
    convexify,
    cluster,
    flatten,

# Lazy operations on flowpipes
    Projection,
    Shift,

# Hybrid types
    HACLD1,
    DiscreteTransition,

# Getter functions for hybrid systems
    jitter,
    switching,
    location,
    reset_map,
    guard,
    source_invariant,
    target_invariant,

# Algorithms for intersection operations
    FallbackIntersection,
    HRepIntersection,
    BoxIntersection,
    TemplateHullIntersection,

# Algorithms for disjointness operations
    FallbackDisjointness,
    NoEnclosure,
    BoxEnclosure,
    ZonotopeEnclosure,
    Dummy,

# Algorithm to check inclusion
    FallbackInclusion,

# Clustering methods
    NoClustering,
    LazyClustering,
    UnionClustering,
    BoxClustering,
    ZonotopeClustering


text:
# Calculates eps2 and eps4, the second and fourth-order artificial diffusion coefficients used in the JST method
function unstructured_JSTEps(mesh::Mesh, sln::SolutionState, fluid::Fluid, k2=0.5, k4=(1/32), c4=1)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)

    # Calc Pressure Gradient
    @views P = sln.cellPrimitives[:,1]
    P = reshape(P, nCells, :)
    gradP = greenGaussGrad(mesh, P, false)
    gradP = reshape(gradP, nCells, 3)

    sj = zeros(nCells) # 'Sensor' used to detect shock waves and apply second-order artificial diffusion to stabilize solution in their vicinity
    sjCount = zeros(nCells) # Store the number of sj's calculated for each cell, cell-center value will be the average of all of them
    @inbounds @fastmath for f in 1:nFaces-nBdryFaces
        # At each internal face, calculate sj, rj, eps2, eps4
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]
        d = mesh.cCenters[neighbourCell] .- mesh.cCenters[ownerCell]

        # 1. Find pressure at owner/neighbour cells
        oP = P[ownerCell]
        nP = P[neighbourCell]

        # 2. Calculate pressures at 'virtual' far-owner and far-neighbour cells using the pressure gradient (2nd-order)
        @views farOwnerP = nP - 2*dot(d, gradP[ownerCell, :])
        @views farNeighbourP = oP + 2*dot(d, gradP[neighbourCell, :])

        # 3. With the known and virtual values, can calculate sj at each cell center.
        sj[ownerCell] += (abs( nP - 2*oP + farOwnerP )/ max( abs(nP - oP) + abs(oP - farOwnerP), 0.0000000001))^2
        sjCount[ownerCell] += 1
        sj[neighbourCell] += (abs( oP - 2*nP + farNeighbourP )/ max( abs(farNeighbourP - nP) + abs(nP - oP), 0.0000000001))^2
        sjCount[neighbourCell] += 1
    end

    rj = zeros(nCells) # 'Spectral radius' -> maximum possible speed of wave propagation relative to mesh
    @inbounds @fastmath for c in 1:nCells
        @views rj[c] = mag(sln.cellPrimitives[c,3:5]) +  sqrt(fluid.gamma * fluid.R * sln.cellPrimitives[c,2]) # Velocity magnitude + speed of sound
        sj[c] /= sjCount[c] # Average the sj's computed at each face for each cell
    end

    # Values of rj, sj at faces is the maximum of their values at the two adjacent cell centers
    rjsjF = maxInterp(mesh, rj, sj) # column one is rj, column two is sj, both at face centers

    # Calculate eps2 and eps4
    eps2 = zeros(nFaces)
    eps4 = zeros(nFaces)
    for f in 1:nFaces-nBdryFaces
        eps2[f] = k2 * rjsjF[f,2] * rjsjF[f,1]
        eps4[f] = max(0, k4*rjsjF[f,1] - c4*eps2[f])
    end

    return eps2, eps4
end

#=
    Inputs: Expects that sln.cellState, sln.cellPrimitives and sln.cellFluxes are up-to-date
    Outputs: Updates sln.faceFluxes and sln.cellResiduals
    Returns: Updated sln.cellResiduals

    Applies classical JST method: central differencing + JST artificial diffusion. Each face treated as a 1D problem
    http://aero-comlab.stanford.edu/Papers/jst_2015_updated_07_03_2015.pdf -  see especially pg.5-6
=#
function unstructured_JSTFlux(mesh::Mesh, sln::SolutionState, boundaryConditions, fluid::Fluid)
    nCells, nFaces, nBoundaries, nBdryFaces = unstructuredMeshInfo(mesh)
    nVars = size(sln.cellState, 2)

    #### 1. Apply boundary conditions ####
    applyBoundaryConditions(mesh, sln, boundaryConditions, nBoundaries, fluid)

    #### 2. Centrally differenced fluxes ####
    linInterp_3D(mesh, sln.cellFluxes, sln.faceFluxes)

    #### 3. JST artificial Diffusion ####
    fDeltas = faceDeltas(mesh, sln)
    fDGrads = greenGaussGrad(mesh, fDeltas, false)
    eps2, eps4 = unstructured_JSTEps(mesh, sln, fluid)

    diffusionFlux = zeros(nVars)
    @inbounds @fastmath for f in 1:nFaces-nBdryFaces
        ownerCell = mesh.faces[f][1]
        neighbourCell = mesh.faces[f][2]
        d = mesh.cCenters[neighbourCell] .- mesh.cCenters[ownerCell]

        @views fD = fDeltas[f,:]
        @views farOwnerfD = fD .- dot(d, fDGrads[ownerCell,:,:])
        @views farNeighbourfD = fD .+ dot(d, fDGrads[ownerCell,:,:])

        diffusionFlux = eps2[f]*fD - eps4[f]*(farNeighbourfD - 2*fD + farOwnerfD)

        # Add diffusion flux in component form
        unitFA = normalize(mesh.fAVecs[f])
        for v in 1:nVars
            i1 = (v-1)*3+1
            i2 = i1+2
            sln.faceFluxes[f,i1:i2] .-= (diffusionFlux[v] .* unitFA)
        end
    end

    #### 4. Integrate fluxes at in/out of each cell (sln.faceFluxes) to get change in cell center values (sln.fluxResiduals) ####
    return integrateFluxes_unstructured3D(mesh, sln)
end

max_stars_repo_name: henrystoldt/fvCFD


text:
type NormalArm <: BanditArm
  mu::Float64
  sigma::Float64
end

function draw(arm::NormalArm)
  rand(Normal(mu, sigma))
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
include("./FieldFactory/rayclassgrp.jl")
include("./FieldFactory/ab_exts.jl")
include("./FieldFactory/rcf_brauer.jl")


text:
#--- input
#   nprow : [global] : num of rows in the process grid
#   npcol : [global] : num of cols in the process grid
#--- output
#   ctxt  : [global] : BLACS context
function sl_init(nprow::ScaInt, npcol::ScaInt)
    ctxt = zeros(ScaInt,1)
    ccall((:sl_init_, libscalapack), Nothing,
        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
        ctxt, Ref(nprow), Ref(npcol))
    return ctxt[1]
end

#--- input
#   n        : [global] : num of rows/cols
#   nb       : [global] : num of rows/cols in its block
#   iproc    : [global] : num of process rows and cols in the current process id
#   isrcproc : [global] : head grid address of its process grid
#   nprocs   : [global] : num of process grid
#--- output
#              [global] : num of process rows and cols in the current process id
function numroc(n::ScaInt, nb::ScaInt, iproc::ScaInt, isrcproc::ScaInt, nprocs::ScaInt)
    return ccall((:numroc_, libscalapack), ScaInt,
                (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                Ref(n), Ref(nb), Ref(iproc), Ref(isrcproc), Ref(nprocs))
end

#--- input
#   m/n         : [global] : num of rows/cols
#   mb/nb       : [global] : num of rows/cols in its block
#   irsrc/icsrc : [global] : head of row/col grid address of its process grid
#   ictxt       : [global] : BLACS context
#   lld         : [global] : local leading dimension
#--- output
#   desc        : [global] : array descriptor of the local matrix
function descinit(m::ScaInt, n::ScaInt, mb::ScaInt, nb::ScaInt, irsrc::ScaInt, icsrc::ScaInt, ictxt::ScaInt, lld::ScaInt)
    desc = zeros(ScaInt, 9)
    info = zeros(ScaInt, 1)
    ccall((:descinit_, libscalapack), Nothing,
            (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
            Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
            Ptr{ScaInt}, Ptr{ScaInt}),
            desc, Ref(m), Ref(n), Ref(mb),
            Ref(nb), Ref(irsrc), Ref(icsrc), Ref(ictxt),
            Ref(lld), info)
    if info[1] < 0
        error("input argument $(info[1]) has illegal value")
    end
    return desc
end

#--- input
#   A     : [global] : ( will be updated ! ) input matrix
#   ia/ja : [global] : first row/col index in the global matrix A
#   desca : [global] : descriptor of the local matrix A
#   α     : [global] : scalar value to substitute into the A
#--- output
#                      nothing
for (fname, elty) in ((:pselset_, :Float32),
                      (:pdelset_, :Float64),
                      (:pcelset_, :ComplexF32),
                      (:pzelset_, :ComplexF64))
    @eval begin
        function pXelset!(A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt}, α::$elty)
            ccall(($(string(fname)), libscalapack), Nothing,
                (Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{$elty}),
                A, Ref(ia), Ref(ja), desca, Ref(α))
        end
    end
end

#--- input
#   scope : [global] : BLACS scope in which alpha is returned
#   top   : [global] : topology to be used if broadcast is needed
#   A     : [local]  : local matrix
#   ia/ja : [global] : first row/col index in the global matrix A
#   desca : [global] : descriptor of the local matrix A
#--- output
#   α     : [local]  : scalar value which will be returned from the A
for (fname, elty) in ((:pselget_, :Float32),
                      (:pdelget_, :Float64),
                      (:pcelget_, :ComplexF32),
                      (:pzelget_, :ComplexF64))
    @eval begin
        function pXelget(scope::Char, top::Char, A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt})
            α = zeros($elty,1)
            ccall(($(string(fname)), libscalapack), Nothing,
                  (Ptr{Char}, Ptr{Char}, Ptr{$elty}, Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                  f_pchar(scope), f_pchar(top), α, A, Ref(ia), Ref(ja), desca)
            return α[1]
        end
    end
end

#
for (fname, elty) in ((:pslacpy_, :Float32),
                      (:pdlacpy_, :Float64),
                      (:pclacpy_, :ComplexF32),
                      (:pzlacpy_, :ComplexF64))
    @eval begin
        function pXlacpy!(uplo::Char, m::ScaInt, n::ScaInt,
                          A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt},
                          B::Matrix{$elty}, ib::ScaInt, jb::ScaInt, descb::Vector{ScaInt})
            ccall(($(string(fname)), libscalapack), Nothing,
                    (Ptr{Char}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                    f_pchar(uplo), Ref(m), Ref(n),
                    A, Ref(ia), Ref(ja), desca,
                    B, Ref(ib), Ref(jb), descb)
        end
    end
end

# input : num of rows/cols
#         local matrix A
#         head grid address of its process grid
#         array descriptor
#         ( will be updated ! ) matrix B
#         BLACS context
# output: Nothing
for (fname, elty) in ((:psgemr2d_, :Float32),
                      (:pdgemr2d_, :Float64),
                      (:pcgemr2d_, :ComplexF32),
                      (:pzgemr2d_, :ComplexF64))
    @eval begin
        function pXgemr2d!(m::ScaInt, n::ScaInt,
                           A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt},
                           B::Matrix{$elty}, ib::ScaInt, jb::ScaInt, descb::Vector{ScaInt},
                           ictxt::ScaInt)
            ccall(($(string(fname)), libscalapack), Nothing,
                 (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                  Ptr{ScaInt}, Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                  Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                  Ref(m), Ref(n),
                  A, Ref(ia), Ref(ja), desca,
                  B, Ref(ib), Ref(jb), descb,
                  Ref(ictxt))
        end
    end
end

# input : BLACS scope to operate; op(X) = X or X' or X*'
#         grid indices of the inputted global matrices
#         coefficients
#         local matrices, grid indices of the local matrices
#         array descriptors for the distributed matrices
# output: nothing
# === detail ===
# op(X) = X or X' or X*'
# sub(A) = A[ia:ia+m-1,ja:ja+k-1]
# sub(B) = B[ib:ib+k-1,jb:jb+n-1]
# sub(C) = C[ic:ic+m-1,jc:jc+n-1]
#        = α*op(sub(A))*op(sub(B))+β*sub(C)
# op(sub(A)) denotes A[ia:ia+m-1,ja:ja+k-1]   if transa = 'n',
#                    A[ia:ia+k-1,ja:ja+m-1]'  if transa = 't',
#                    A[ia:ia+k-1,ja:ja+m-1]*' if transa = 'c',
# op(sub(B)) denotes B[ib:ib+k-1,jb:jb+n-1]   if transb = 'n',
#                    B[ib:ib+n-1,jb:jb+k-1]'  if transb = 't',
#                    B[ib:ib+n-1,jb:jb+k-1]*' if transb = 'c',
for (fname, elty) in ((:psgemm_, :Float32),
                      (:pdgemm_, :Float64),
                      (:pcgemm_, :ComplexF32),
                      (:pzgemm_, :ComplexF64))
    @eval begin
        function pXgemm!(transa::Char, transb::Char,
                         m::ScaInt, n::ScaInt, k::ScaInt,
                         α::$elty,
                         A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt},
                         B::Matrix{$elty}, ib::ScaInt, jb::ScaInt, descb::Vector{ScaInt},
                         β::$elty,
                         C::Matrix{$elty}, ic::ScaInt, jc::ScaInt, descc::Vector{ScaInt})

            ccall(($(string(fname)), libscalapack), Nothing,
                (Ptr{Char}, Ptr{Char},
                 Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                 Ptr{$elty},
                 Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                 Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                 Ptr{$elty},
                 Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                 f_pchar(transa), f_pchar(transb),
                 Ref(m), Ref(n), Ref(k),
                 Ref(α),
                 A, Ref(ia), Ref(ja), desca,
                 B, Ref(ib), Ref(jb), descb,
                 Ref(β),
                 C, Ref(ic), Ref(jc), descc)
        end
    end
end

#
for (fname, elty) in ((:pslaset_, :Float32),
                      (:pclaset_, :Float64),
                      (:pdlaset_, :ComplexF32),
                      (:pzlaset_, :ComplexF64))
    @eval begin
        function pXlaset!(uplo::Char, m::ScaInt, n::ScaInt,
                          α::$elty, β::$elty,
                          A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt})
            ccall(($(string(fname)), libscalapack), Nothing,
                    (Ptr{Char}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{$elty},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}),
                    f_pchar(uplo), Ref(m), Ref(n),
                    Ref(α), Ref(β),
                    A, Ref(ia), Ref(ja), desca)
        end
    end
end

#
for (fname, elty, elty_s) in ((:psgebal_, :Float32, :Float32),
                              (:pdgebal_, :Float64, :Float64))
    @eval begin
        function pXgebal!(job::Char, n::ScaInt,
                          A::Matrix{$elty}, desca::Vector{ScaInt}, ilo::ScaInt, ihi::ScaInt,
                          scale::Vector{$elty_s})

            # inner variable
            info = zeros(ScaInt,1)
            ccall(($(string(fname)), libscalapack), Nothing,
                    (Ptr{Char}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty_s}, Ptr{ScaInt}),
                    f_pchar(job), Ref(n),
                    A, desca, Ref(ilo), Ref(ihi),
                    scale, info)

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

        end
    end
end

# ! assuming ! the A is a square matrix of n x n
# input : order of the local matrix A
#         lower/heigher range of the rows/cols in the global A
#         ( will be updated ! ) local matrix A
#         first row/col address of A
#         array descriptor for A
#         ( will be updated ! ) see 'the "τ" meaning'
#         ( will be updated ! ) see 'the "work" meaning'
#         see 'the "lwork" meaning'
# output: nothing
# === detail ===
# the "τ" meaning:
#     dim(τ) = numroc(ja+n-2,numblocks,mycol,csrc_a,npcolcols)
#     τ[ja+ilo-1:ja+ihi-2] = the scalar factors of the elementary reflectors
#     τ[ja:ja+ilo-2] = τ[ja+ihi-1:ja+n-2] = 0
# the "work" meaning:
#     dim(work) = lwork
#     if lwork == 0, work is ignored.
#     else,
#         if lwork != -1, its size is (at least) of length lwork.
#         else, its size is ( at least ) of length 1.
# the "lwork" meaning:
#     if lwork = -1, lwork is global
#     if lwork >= 0, lwork is local that value must be at least
#       lwork >= numblocks^2 + numblocks*max(ihip+1,ihlp+inlq)
#       where:
#           iroffa = (ia-1) % numblocks
#           icoffa = (ja-1) % numblocks
#           ioff = (ia+ilo-2) % numblocks
#           iarow = indxg2p(ia,numblocks,myrow,rsrc_a,npcolrows)
#                 = (rsrc_a+(ia-1)/numblocks) % npcolrows
#           ihip = numroc(ihi+iroffa,numblocks,myrow,iarow,npcolrows)
#           ilrow = indxg2p(ia+ilo-1,numblocks,myrow,rsrc_a,npcolrows)
#                 = (rsrc_a+(ia+ilo-2)/numblocks) % npcolrows
#           ihlp = numroc(ihi-ilo+ioff+1,numblocks,myrow,ilrow,npcolrows)
#           ilcol = indxg2p(ja+ilo-1,numblocks,mycol,csrc_a,npcolocls)
#                 = (csrc_a+(ja+ilo-2)/numblocks) % npcolcols
#           inlq = numroc(n-ilo+ioff+1,numblocks,mycol,ilcol,npcolcols)
for (fname, elty) in ((:psgehrd_, :Float32),
                      (:pdgehrd_, :Float64),
                      (:pcgehrd_, :ComplexF32),
                      (:pzgehrd_, :ComplexF64))
    @eval begin
        function pXgehrd!(n::ScaInt, ilo::ScaInt, ihi::ScaInt,
                          A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt},
                          τ::Vector{$elty})

            # inner variables
            info = zeros(ScaInt,1)
            work = zeros($elty,1); lwork = -1;
            # j = 1 for perform a workspace query
            # j = 2 for perform ccall
            for j = 1:2
                ccall(($(string(fname)), libscalapack), Nothing,
                        (Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}),
                        Ref(n), Ref(ilo), Ref(ihi),
                        A, Ref(ia), Ref(ja), desca,
                        τ, work, Ref(lwork), info)
                # allocate vector for j = 2
                if j == 1
                    lwork = convert(ScaInt, work[1])
                    work = zeros($elty, lwork)
                end
            end

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

        end
    end
end

#
for (fname, elty) in ((:psormhr_, :Float32),
                      (:pdormhr_, :Float64),
                      (:pcunmhr_, :ComplexF32),
                      (:pzunmhr_, :ComplexF64))
    @eval begin
        function pXYYmhr!(side::Char, trans::Char,
                          m::ScaInt, n::ScaInt, ilo::ScaInt, ihi::ScaInt,
                          A::Matrix{$elty}, ia::ScaInt, ja::ScaInt, desca::Vector{ScaInt},
                          τ::Vector{$elty},
                          C::Matrix{$elty}, ic::ScaInt, jc::ScaInt, descc::Vector{ScaInt})

            # inner variables
            info = zeros(ScaInt,1)
            work = zeros($elty,1); lwork = -1;
            # j = 1 for perform a workspace query
            # j = 2 for perform ccall
            for j = 1:2
                ccall(($(string(fname)), libscalapack), Nothing,
                    (Ptr{Char}, Ptr{Char},
                    Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}),
                    f_pchar(side), f_pchar(trans),
                    Ref(m), Ref(n), Ref(ilo), Ref(ihi),
                    A, Ref(ia), Ref(ja), desca,
                    τ,
                    C, Ref(ic), Ref(jc), descc,
                    work, Ref(lwork), info)
                # allocate vector for j = 2
                if j == 1
                    lwork = convert(ScaInt, work[1])
                    work = zeros($elty, lwork)
                end
            end

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

        end
    end
end

# ! assuming ! the A is an upper Hessenberg matrix of n x n
# input : booleans that requirement return form
#         order of the global Hessenberg matrix A
#         lower/heigher range of the rows/cols in the global A
#         ( will be updated ! ) global Hessenberg matrix A
#         array descriptor for A
#         index of the rows of Z that specified; Z[iloz:ihiz]
#         real/imag parts of the computed eigenvalues[ilo:ihi]
#         ( will be updated ! ) the Schur vector matrix Z; update part is [iloz:ihiz, ilo:ihi]
#         array descriptor for Z
#         ( will be updated ! ) see 'the "work" meaning'
#         see 'the "lwork" meaning'
#         temporary integer arrays and that order
# output: nothing
# === detail ===
# the "wantt" requires:
#   true ; the full Schur form
#   false; only eigenvalues
# the "wantz" requires:
#   true ; the matrix of Schur vectors Z
#   false; no Z
# the "work" meaning:
#     dim(work) = lwork
# the "lwork" meaning:
#     lwork >= 3*n + max(2*max(lld_z,lld_a)+2*locc(n),7*ceil(n/hbl)/lcm(npcolrows,npcolcols))
#     where:
#           lld_z = descz[9], lld_a = desca[9], csrc_a = desca[8]
#           locc(n) = numroc(n,numblocks,mycol=npcol,csrc_a,npcolcols)
#           hbl = numblocks ( = mblocks = nblocks )
for (fname, elty) in ((:pslaqr1_, :Float32),
                      (:pdlaqr1_, :Float64))
# for (fname, elty) in ((:pslahqr_, :Float32),
#                       (:pdlahqr_, :Float64))
        @eval begin
        function pXlahqr!(wantt::Bool, wantz::Bool,
                          n::ScaInt, ilo::ScaInt, ihi::ScaInt,
                          A::Matrix{$elty}, desca::Vector{ScaInt},
                          wr::Vector{$elty}, wi::Vector{$elty},
                          iloz::ScaInt, ihiz::ScaInt, Z::Matrix{$elty}, descz::Vector{ScaInt})

            # inner variables
            info = zeros(ScaInt,1)
            work = zeros($elty,1); lwork = convert(ScaInt, -1);
            iwork = zeros(ScaInt,1); ilwork = convert(ScaInt, -1);
            # j = 1 for perform a workspace query
            # j = 2 for perform ccall
            for j = 1:2
                ccall(($(string(fname)), libscalapack), Nothing,
                        (Ptr{Bool}, Ptr{Bool},
                        Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{$elty},
                        Ptr{ScaInt}, Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{ScaInt}),
                        Ref(wantt), Ref(wantz),
                        Ref(n), Ref(ilo), Ref(ihi),
                        A, desca, wr, wi,
                        Ref(iloz), Ref(ihiz), Z, descz,
                        work, Ref(lwork), iwork, Ref(ilwork),
                        info)
                # allocate vector for j = 2
                if j == 1
                    lwork = convert(ScaInt, work[1])
                    ilwork = convert(ScaInt, iwork[1])
                    # ilwork = convert(ScaInt, 1000)
                    work = zeros($elty, lwork)
                    iwork = zeros(ScaInt, ilwork)
                end
            end

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

        end
    end
end
for (fname, elty) in ((:pclahqr_, :ComplexF32),
                      (:pzlahqr_, :ComplexF64))
    @eval begin
        function pXlahqr!(wantt::Bool, wantz::Bool,
                            n::ScaInt, ilo::ScaInt, ihi::ScaInt,
                            A::Matrix{$elty}, desca::Vector{ScaInt},
                            w::Vector{$elty},
                            iloz::ScaInt, ihiz::ScaInt, Z::Matrix{$elty}, descz::Vector{ScaInt})

            # inner variables
            info = zeros(ScaInt,1)
            work = zeros($elty,1); lwork = convert(ScaInt, -1);
            iwork = zeros(ScaInt,1); ilwork = convert(ScaInt, -1);
            # j = 1 for perform a workspace query
            # j = 2 for perform ccall
            for j = 1:2
                ccall(($(string(fname)), libscalapack), Nothing,
                        (Ptr{Bool}, Ptr{Bool},
                        Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{ScaInt},
                        Ptr{$elty},
                        Ptr{ScaInt}, Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                        Ptr{$elty}, Ptr{ScaInt}, Ptr{ScaInt}, Ptr{ScaInt},
                        Ptr{ScaInt}),
                        Ref(wantt), Ref(wantz),
                        Ref(n), Ref(ilo), Ref(ihi),
                        A, desca, w,
                        Ref(iloz), Ref(ihiz), Z, descz,
                        work, Ref(lwork), iwork, Ref(ilwork),
                        info)
                # allocate vector for j = 2
                if j == 1
                    lwork = convert(ScaInt, work[1])
                    ilwork = convert(ScaInt, iwork[1])
                    work = zeros($elty, lwork)
                    iwork = zeros(ScaInt, ilwork)
                end
            end

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

        end
    end
end

#
for (fname, elty, elty_r) in ((:pctrevc_, :ComplexF32, :Float32),
                              (:pztrevc_, :ComplexF64, :Float64))
    @eval begin
        function pXtrevc!(side::Char, howmny::Char, select::Vector{Bool},
                           n::ScaInt, T::Matrix{$elty}, desct::Vector{ScaInt},
                           vl::Matrix{$elty}, descvl::Vector{ScaInt}, vr::Matrix{$elty}, descvr::Vector{ScaInt})

            # inner variables
            info = zeros(ScaInt, 1)
            mm = n
            work = zeros($elty, 2*n)
            rwork = zeros($elty_r, n)
            # output variable
            m = zeros(ScaInt, 1)
            ccall(($(string(fname)), libscalapack), Nothing,
                    (Ptr{Char}, Ptr{Char}, Ptr{Bool},
                    Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{ScaInt}, Ptr{$elty}, Ptr{ScaInt},
                    Ptr{ScaInt}, Ptr{ScaInt},
                    Ptr{$elty}, Ptr{$elty_r}, Ptr{ScaInt}),
                    f_pchar(side), f_pchar(howmny), select,
                    Ref(n), T, desct,
                    vl, descvl, vr, descvr,
                    Ref(mm), m,
                    work, rwork, info)

            if info[1] < 0
                error("input argument $(info[1]) has illegal value")
            end

            return m[1]

        end
    end
end


text:
using Genie.Configuration

export config

const config =  Settings(
                  output_length           = 100,
                  suppress_output         = false,
                  log_db                  = false,
                  log_queries             = false,
                  log_requests            = false,
                  log_responses           = false,
                  log_formatted           = true,
                  log_level               = :debug,
                  log_cache               = false,
                  log_views               = false,
                  assets_path             = "/",
                  cache_duration          = 0,
                  flax_compile_templates  = false,
                  websocket_server        = false,
                  session_auto_start      = false
                )

max_stars_repo_name: essenciary/genie-demo-hello-world

max_issues_repo_name: essenciary/genie-demo-hello-world

max_forks_repo_name: essenciary/genie-demo-hello-world


text:
function free_energy_nopool(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = 1)
    E_v_dense = energy(rbm.visible_dense, v_dense)
    E_v_conv = energy(rbm.visible_conv, v_conv)
    I = inputs_v_to_h(rbm, v_dense, v_conv)
    F = RBMs.free_energy(hidden(rbm), I; β)

    vsz_conv = vsizes(parts(rbm).conv, v_conv)
    @assert size(E_v_dense) == batch_size(rbm.visible_dense, v_dense)
    @assert size(E_v_conv) == (vsz_conv.input_size..., vsz_conv.batch_size...)
    @assert size(F) == (hsizes(rbm, I).output_size..., hsizes(rbm, I).batch_size...)

    # reduce over kernel dims
    E_vconv_reduced = reshape_maybe(sum(E_v_conv; dims=1:kernel_ndims(rbm)), vsz_conv.batch_size)
    F_reduced = reshape_maybe(sum(F; dims=1:kernel_ndims(rbm)), hsizes(rbm, I).batch_size)

    return E_v_dense .+ E_vconv_reduced .+ F_reduced
end

function free_energy_pooled(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = 1)
    E_v_dense = energy(rbm.visible_dense, v_dense)
    E_v_conv = energy(rbm.visible_conv, v_conv)
    I = inputs_v_to_h(rbm, v_dense, v_conv)
    F = RBMs.free_energies(hidden(rbm), I; β)

    vsz_conv = vsizes(parts(rbm).conv, v_conv)
    @assert size(E_v_dense) == batch_size(rbm.visible_dense, v_dense)
    @assert size(E_v_conv) == (vsz_conv.input_size..., vsz_conv.batch_size...)
    @assert size(F) == size(I) == (size(hidden(rbm))..., hsizes(rbm, I).output_size..., hsizes(rbm, I).batch_size...)
    G = -logsumexp(-β * F; dims=output_dims(rbm)) / β # pooling step <=> only one h active per group
    @assert size(G)[1:ndims(hidden(rbm))] == size(hidden(rbm))
    @assert all(size(G)[output_dims(rbm)] .== 1)

    # reduce over remaining dims
    E_vconv_reduced = reshape_maybe(sum(E_v_conv; dims=1:kernel_ndims(rbm)), vsz_conv.batch_size)
    F_reduced = reshape_maybe(sum(G; dims=1:ndims(hidden(rbm))), hsizes(rbm, I).batch_size)

    return E_v_dense .+ E_vconv_reduced .+ F_reduced
end

function sample_h_from_v_nopool(
    rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = true
)
    inputs = inputs_v_to_h(rbm, v_dense, v_conv)
    return transfer_sample(hidden(rbm), inputs; β)
end

function sample_h_from_v_pooled(rbm::DenseConvRBM, v_dense::AbstractArray, v_conv::AbstractArray; β::Real = 1)
    I = inputs_v_to_h(rbm, v_dense, v_conv)
    F = RBMs.free_energies(hidden(rbm), I; β)
    @assert size(F) == size(I) == (size(hidden(rbm))..., hsizes(rbm, I).output_size..., hsizes(rbm, I).batch_size...)
    selected_size = (size(hidden(rbm))..., hsizes(rbm, I).batch_size...)
    k = reshape(cartesian_sample_from_logits_gumbel(-β * F; dims=output_dims(rbm)), selected_size)

    selected_inputs = similar(I, selected_size)
    for n in CartesianIndices(hsizes(rbm, I).batch_size), μ in CartesianIndices(size(hidden(rbm)))
        @inbounds selected_inputs[μ,n] = I[k[μ,n]]
    end

    selected_h = RBMs.transfer_sample(hidden(rbm), selected_inputs; β)
    @assert size(selected_h) == size(selected_inputs) == selected_size

    h = fill!(similar(selected_h, size(F)), 0)
    for n in CartesianIndices(hsizes(rbm, I).batch_size), μ in CartesianIndices(size(hidden(rbm)))
        h[k[μ,n]] = selected_h[μ,n]
    end
    return h
end


text:
export ridgenorm, rdist

function rdist(;scale,freq)
  ((as,af),(bs,bf)) -> (as - bs)^2 / (scale^2) + (af - bf)^2 / (freq^2)
end

logdet_(x,n) = logdet(x)
logdet_(x::UniformScaling,n) = n*log(abs(x.λ))
factorize_(x) = factorize(x)
factorize_(x::UniformScaling) = x

function logpdf_mvt(v,μ,Σ,x)
  Σ_ = factorize_(Σ)
  d = length(x)
  C = (lgamma((v+d)/2)) - (lgamma(v/2)+(d/2)log(v*π)) - (0.5 * (logdet_(Σ_,d)))
  diff = abs.(μ.-x)
  # with poorly conditioned covariances, normdiff can be negative even when the
  # log determant is a valid, non-infinite number
  # NOTE: this somewhat akward notation is due to some missing method
  # definitions in the sparse matrix library (only right divide is defined)
  normdiff = max(0,(diff'*(Σ_'\diff))[1])
  result = C + -(v+d)/2*log(1+normdiff/v)
  @assert !isinf(result)
  result
end

mutable struct RidgeMultiNormalStats{T,C} <: Stats{T}
  μ::Vector{T}
  S::Vector{T}
  n::Float64
  x2_offset::Float64
  corr::C
end
Statistics.std(x::RidgeMultiNormalStats) = x.n > 0 ? x.S ./ x.n : Inf
Statistics.zero(x::RidgeMultiNormalStats{T},C::Coherence) where T =
  RidgeMultiNormalStats(zero(x.μ),zero(x.S),0.0,0.0,x.corr)

function findcorr(dims,dist,thresh)
  n = prod(dims)
  corr = spzeros(n,n)
  # TODO: consider using a ToeplitzMatrix
  # here
  for (i,ii) in enumerate(CartesianIndices(dims))
    for (j,jj) in enumerate(CartesianIndices(dims))
      if i <= j
        val = exp(-dist(ii.I,jj.I))
        if val > thresh
          corr[i,j] = val
        end
      end
    end
  end
  # @info "Sparsity of source-prior ($(size(corr,1)) × $(size(corr,2))) ridge: "*
  #   "$(nnz(corr)) ($(round(Int,100nnz(corr) / length(corr)))%)" maxlog=1

  Symmetric(corr)
end

function ridgenorm(prior::Coherence,n::Number,
                   x2_offset::Number=size(prior,2)*size(prior,3);
                   thresh=1e-3,scale=nothing,
                   freq=nothing)
  error("Outdated prior")
  ridgenorm(prior,n,rdist(scale=scale,freq=freq),x2_offset,thresh=thresh)
end

function ridgenorm(prior::Coherence,n::Number,dist::Function,
                   x2_offset::Number=size(prior,2)*size(prior,3);thresh=1e-3)
  error("Outdated prior")
  data = reshape(mean(prior,4),size(prior,1),:)
  μ = squeeze(mean(data,1),1)
  S = squeeze(sum(data.^2,1),1) .* n./size(data,1) .+ n.*thresh.^2
  corr = findcorr((size(prior,2),size(prior,3)),dist)
  RidgeMultiNormalStats{eltype(data)}(μ,S,n,x2_offset,corr)
end

# https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm
function update!(stats::RidgeMultiNormalStats{T},x::AbstractVector{T},
                 w=1.0) where T
  stats.n += w

  μ₀ = copy(stats.μ)
  stats.μ .+= (x .- stats.μ).*(w ./ stats.n)
  stats.S .+= (x .- stats.μ).*(x .- μ₀).*w

  stats
end

function downdate!(stats::RidgeMultiNormalStats{T},x::AbstractVector{T},
                   w=1.0) where T
  if w ≈ stats.n
    stats.μ .= 0
    stats.S .= 0
    stats.n = 0
  elseif w < stats.n
    μₙ = copy(stats.μ)
    stats.μ .= (stats.μ .- x.*(w ./ stats.n)) ./ (1 .- (w ./ stats.n))
    stats.S .-= (x .- stats.μ).*(x .- μₙ).*w

    stats.n -= w
  else
    error("Sum of weights would become negative. This happens when trying ",
          "to remove more samples then were added using `update!`.")
  end

  stats
end

function mult!(stats::RidgeMultiNormalStats,c)
  stats.n *= c

  stats
end

function Base.:(+)(a::RidgeMultiNormalStats{T},b::RidgeMultiNormalStats{T}) where T
  n = a.n + b.n
  RidgeMultiNormalStats(a.μ.*(a.n./n) .+ b.μ.*(b.n./n), a.S.+b.S,
                        n,a.x2_offset+b.x2_offset,a.corr)
end

function logpdf(stats::RidgeMultiNormalStats,x::AbstractVector)
  d = length(stats.μ)
  σ = sqrt.(stats.S ./ stats.n)
  nu = stats.n + stats.x2_offset - d + 1

  c = ((stats.n + 1)/(stats.n * nu))
  I,J,V = findnz(stats.corr.data)
  Λ = Symmetric(sparse(I,J,[v*σ[i]*σ[j] * c for (i,j,v) in zip(I,J,V)]))
  # Λ = ((Diagonal(σ) * stats.corr) * Diagonal(σ)) .* ((stats.n + 1)/(stats.n * nu))

  if nu < 1
    if length(x) != length(stats.μ)
      error("Dimension mistmatch between data ($(length(x))) and ",
            "the model ($(length(stats.μ))).")
    else
      error("Insufficient data for predictive distribution, ",
            "try collecting more data or increasing x2_offset")
    end
  end

  logpdf_mvt(nu,stats.μ,Λ,x)
end

struct ConstRidgePrior{T,C} <: Stats{T}
  S::T
  N::Float64
  x2_offset::Float64
  corr::C
end

function ridgenorm(prior::Number,N::Number,dims::Tuple,dist::Function,
                   x2_offset::Number=prod(dims);threshold)
  ConstRidgePrior(float(prior),Float64(N),Float64(x2_offset),
                  findcorr(dims,dist,threshold))
end

function ridgenorm(prior::Number,N::Number,dims::Tuple,
                   x2_offset::Number=prod(dims);scale,freq,threshold)
  dist = rdist(scale=scale,freq=freq)
  ridgenorm(float(prior),N,dims,dist,x2_offset,threshold=threshold)
end

function Base.zero(prior::ConstRidgePrior,d)
  RidgeMultiNormalStats(fill(zero(prior.S),d),fill(zero(prior.S),d),0.0,0.0,prior.corr)
end

function Base.:(+)(a::ConstRidgePrior{T},b::RidgeMultiNormalStats{T}) where T
  d = length(b.μ)
  n = b.n + a.N
  RidgeMultiNormalStats(b.μ.*(b.n/n),a.S.+b.S,n,a.x2_offset+b.x2_offset,
                        a.corr)
end

function logpdf(stats::ConstRidgePrior,x::AbstractVector)
  d = length(x)
  σ² = stats.S / stats.N
  Λ = stats.corr .* σ² .* ((stats.N + 1)/(stats.N * nu))
  nu = max(1,floor(Int,stats.N + stats.x2_offset - d + 1))

  logpdf_mvt(nu,0,Λ,x)
end

max_stars_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_issues_repo_name: JuliaTagBot/AuditoryBistabilityLE

max_forks_repo_name: JuliaTagBot/AuditoryBistabilityLE


text:
# some simpler (and faster) implementations for root finding
#
# Not exported
#
# These avoid the setup costs of the `find_zero` method, so should be faster
# though they will take similar number of function calls.
#
# `Roots.bisection(f, a, b)`  (Bisection).
# `Roots.secant_method(f, xs)` (Order1) secant method
# `Roots.dfree(f, xs)`  (Order0) more robust secant method
#

## Bisection
##
## Essentially from Jason Merrill https://gist.github.com/jwmerrill/9012954
## cf. http://squishythinking.com/2014/02/22/bisecting-floats/
## This also borrows a trick from https://discourse.julialang.org/t/simple-and-fast-bisection/14886
## where we keep x1 so that y1 is negative, and x2 so that y2 is positive
## this allows the use of signbit over y1*y2 < 0 which avoid < and a multiplication
## this has a small, but noticeable impact on performance.
"""
    bisection(f, a, b; [xatol, xrtol])

Performs bisection method to find a zero of a continuous
function.

It is assumed that (a,b) is a bracket, that is, the function has
different signs at a and b. The interval (a,b) is converted to floating point
and shrunk when a or b is infinite. The function f may be infinite for
the typical case. If f is not continuous, the algorithm may find
jumping points over the x axis, not just zeros.

If non-trivial tolerances are specified, the process will terminate
when the bracket (a,b) satisfies `isapprox(a, b, atol=xatol,
rtol=xrtol)`. For zero tolerances, the default, for Float64, Float32,
or Float16 values, the process will terminate at a value `x` with
`f(x)=0` or `f(x)*f(prevfloat(x)) < 0 ` or `f(x) * f(nextfloat(x)) <
0`. For other number types, the A42 method is used.

"""
function bisection(f, a::Number, b::Number; xatol=nothing, xrtol=nothing)

    x1, x2 = adjust_bracket(float.((a,b)))
    T = eltype(x1)

    atol = xatol == nothing ? zero(T) : abs(xatol)
    rtol = xrtol == nothing ? zero(one(T)) : abs(xrtol)
    CT = iszero(atol) && iszero(rtol) ?  Val(:exact) : Val(:inexact)

    x1, x2 = float(x1), float(x2)
    y1, y2 = f(x1), f(x2)

    _unitless(y1 * y2) >= 0  && error("the interval provided does not bracket a root")

    if isneg(y2)
        x1, x2, y1, y2 = x2, x1, y2, y1
    end

    xm = Roots._middle(x1, x2) # for possibly mixed sign x1, x2
    ym = f(xm)

    while true

        if has_converged(CT, x1, x2, xm, ym, atol, rtol)
            return xm
        end

        if isneg(ym)
            x1, y1 = xm, ym
        else
            x2, y2 = xm, ym
        end

        xm = Roots.__middle(x1,x2)
        ym = f(xm)

    end

end

# -0.0 not returned by __middle, so isneg true on [-Inf, 0.0)
@inline isneg(x::T) where {T <: AbstractFloat} = signbit(x)
@inline isneg(x) = _unitless(x) < 0

@inline function has_converged(::Val{:exact}, x1, x2, m, ym, atol, rtol)
    iszero(ym) && return true
    isnan(ym) && return true
    x1 != m && m != x2 && return false
    return true
end

@inline function has_converged(::Val{:inexact}, x1, x2, m, ym, atol, rtol)
    iszero(ym) && return true
    isnan(ym) && return true
    val = abs(x1 - x2) <= atol + max(abs(x1), abs(x2)) * rtol
    return val
end

"""
    secant_method(f, xs; [atol=0.0, rtol=8eps(), maxevals=1000])

Perform secant method to solve f(x) = 0.

The secant method is an iterative method with update step
given by b - fb/m where m is the slope of the secant line between
(a,fa) and (b,fb).

The inital values can be specified as a pair of 2, as in `(a,b)` or
`[a,b]`, or as a single value, in which case a value of `b` is chosen.

The algorithm returns m when `abs(fm) <= max(atol, abs(m) * rtol)`.
If this doesn't occur before `maxevals` steps or the algorithm
encounters an issue, a value of `NaN` is returned. If too many steps
are taken, the current value is checked to see if there is a sign
change for neighboring floating point values.

The `Order1` method for `find_zero` also implements the secant
method. This one will be faster, as there are fewer setup costs.

Examples:

```julia
Roots.secant_method(sin, (3,4))
Roots.secant_method(x -> x^5 -x - 1, 1.1)
```

Note:

This function will specialize on the function `f`, so that the inital
call can take more time than a call to the `Order1()` method, though
subsequent calls will be much faster.  Using `FunctionWrappers.jl` can
ensure that the initial call is also equally as fast as subsequent
ones.

"""
function secant_method(f, xs; atol=zero(float(real(first(xs)))), rtol=8eps(one(float(real(first(xs))))), maxevals=100)

    if length(xs) == 1 # secant needs x0, x1; only x0 given
        a = float(xs[1])

        h = eps(one(real(a)))^(1/3)
        da = h*oneunit(a) + abs(a)*h^2 # adjust for if eps(a) > h
        b = a + da

    else
        a, b = promote(float(xs[1]), float(xs[2]))
    end
    secant(f, a, b, atol, rtol, maxevals)
end

function secant(f, a::T, b::T, atol=zero(T), rtol=8eps(T), maxevals=100) where {T}
    nan = (0a)/(0a)
    cnt = 0

    fa, fb = f(a), f(b)
    fb == fa && return nan

    uatol = atol / oneunit(atol) * oneunit(real(a))
    adjustunit = oneunit(real(fb))/oneunit(real(b))

    while cnt < maxevals
        m = b - (b-a)*fb/(fb - fa)
        fm = f(m)

        iszero(fm) && return m
        isnan(fm) && return nan
        abs(fm) <= adjustunit * max(uatol, abs(m) * rtol) && return m
        if fm == fb
            sign(fm) * sign(f(nextfloat(m))) <= 0 && return m
            sign(fm) * sign(f(prevfloat(m))) <= 0 && return m
            return nan
        end

        a,b,fa,fb = b,m,fb,fm

        cnt += 1
    end

    return nan
end

"""
    muller(f, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)
    muller(f, xᵢ₋₂, xᵢ₋₁, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)

> *Muller’s method* generalizes the secant method, but uses quadratic
> interpolation among three points instead of linear interpolation between two.
> Solving for the zeros of the quadratic allows the method to find complex
> pairs of roots.
> Given three previous guesses for the root `xᵢ₋₂`, `xᵢ₋₁`, `xᵢ`, and the values
> of the polynomial `f` at those points, the next approximation `xᵢ₊₁` is produced.

Excerpt and the algorithm taken from

> W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery
> *Numerical Recipes in C*, Cambridge University Press (2002), p. 371

Convergence here is decided by `xᵢ₊₁ ≈ xᵢ` using the tolerances specified,
which both default to `eps(one(typeof(abs(xᵢ))))^4/5` in the appropriate units.
Each iteration performs three evaluations of `f`.
The first method picks two remaining points at random in relative proximity of `xᵢ`.

Note that the method may return complex result even for real intial values
as this depends on the function.

Examples:
```
muller(x->x^3-1, 0.5, 0.5im, -0.5) # → -0.500 + 0.866…im
muller(x->x^2+2, 0.0, 0.5, 1.0) # → ≈ 0.00 - 1.41…im
muller(x->(x-5)*x*(x+5), rand(3)...) # → ≈ 0.00
muller(x->x^3-1, 1.5, 1.0, 2.0) # → 2.0, Not converged
```

"""
muller(f, x₀::T; kwargs...) where T = muller(f, (rand(T, 2).*x₀)..., x₀; kwargs...)

muller(f, xᵢ₋₂, xᵢ₋₁, xᵢ; kwargs...) =
    muller(f, promote(xᵢ₋₂, xᵢ₋₁, xᵢ)...; kwargs...)

function muller(f, oldest::T, older::T, old::T;
    xatol=nothing, xrtol=nothing, maxevals=300) where {T}
    @assert old ≠ older ≠ oldest ≠ old # we want q to be non-degenerate
    xᵢ₋₂, xᵢ₋₁, xᵢ = oldest, older, old
    fxᵢ₋₂, fxᵢ₋₁ = f(xᵢ₋₂), f(xᵢ₋₁)

    RT = typeof(abs(oldest))
    atol = xatol !== nothing ? xatol : oneunit(RT) * (eps(one(RT)))^(4/5)
    rtol = xrtol !== nothing ? xrtol : eps(one(RT))^(4/5)

    for i in 1:maxevals÷3
        # one evaluation per iteration
        fxᵢ = f(xᵢ)
        x = muller_step(xᵢ₋₂, xᵢ₋₁, xᵢ, fxᵢ₋₂, fxᵢ₋₁, fxᵢ)

        if isnan(x)
            @warn "The algorithm might not have converged, stopping at i=$i:" abs(xᵢ - xᵢ₋₁)
            return xᵢ
        end

        # @debug "Iteration $i:" xᵢ₋₂ xᵢ₋₁ xᵢ x abs(x-xᵢ)
        xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, x
        fxᵢ₋₂, fxᵢ₋₁ = fxᵢ₋₁, fxᵢ
        #stopping criterion
        isapprox(xᵢ, xᵢ₋₁, atol=atol, rtol=rtol) && return xᵢ
    end
    @warn "The algorithm might not have converged, maxevals limit hit:" abs(xᵢ₋₁- xᵢ)
    return xᵢ
end

function muller_step(a,b,c,fa,fb,fc)
    a,b,c = promote(a,b,c)
    q = qq(a, b, c)
    q² = q^2
    q1 = q+one(q)

    A =      q*fc - q*q1*fb + q²*fa
    B = (q1+q)*fc - q1^2*fb + q²*fa
    C =     q1*fc

    den = let
        Δ = B^2 - 4A*C
        typeof(Δ) <: Real && Δ < 0 && throw(
            DomainError(Δ, "Discriminant is negative and the function most likely has complex roots. You might want to call muller with complex input."))
        Δ = √Δ
        d⁺ = B + Δ
        d⁻ = B - Δ
        abs(d⁺) > abs(d⁻) ? d⁺ : d⁻
    end
    return c - (c - b)*2C/den
end

@inline qq(a, b, c) = (c - b)/(b - a)

"""
    newton((f, f'), x0; xatol=nothing, xrtol=nothing, maxevals=100)
    newton(fΔf, x0; xatol=nothing, xrtol=nothing, maxevals=100)

Newton's method.

Function may be passed in as a tuple (f, f') *or* as function which returns (f,f/f').

Examples:
```
newton((sin, cos), 3.0)
newton(x -> (sin(x), sin(x)/cos(x)), 3.0, xatol=1e-10, xrtol=1e-10)
```

Note: unlike the call `newton(f, fp, x0)`--which dispatches to a method of `find_zero`, these
two interfaces will specialize on the function that is passed in. This means, these functions
will be faster for subsequent calls, but may be slower for an initial call.

Convergence here is decided by x_n ≈ x_{n-1} using the tolerances specified, which both default to
`eps(T)^4/5` in the appropriate units.

"""
struct TupleWrapper{F, Fp}
f::F
fp::Fp
end
(F::TupleWrapper)(x) = begin
    u, v = F.f(x), F.fp(x)
    return (u, u/v)
end

newton(f::Tuple, x0; kwargs...) = newton(TupleWrapper(f[1],f[2]), x0; kwargs...)
function newton(f, x0; xatol=nothing, xrtol=nothing, maxevals = 100)

    x = float(x0)
    T = typeof(x)
    atol = xatol != nothing ? xatol : oneunit(T) * (eps(one(T)))^(4/5)
    rtol = xrtol != nothing ? xrtol : eps(one(T))^(4/5)

    xo = Inf
    for i in 1:maxevals

        fx, Δx = f(x)
        iszero(fx) && return x

        x -= Δx

        if isapprox(x, xo, atol=atol, rtol=rtol)
            return x
        end

        xo = x
    end

    error("No convergence")
end

## This is basically Order0(), but with different, default, tolerances employed
## It takes more function calls, but works harder to find exact zeros
## where exact means either iszero(fx), adjacent floats have sign change, or
## abs(fxn) <= 8 eps(xn)
"""
    dfree(f, xs)

A more robust secant method implementation

Solve for `f(x) = 0` using an alogorithm from *Personal Calculator Has Key
to Solve Any Equation f(x) = 0*, the SOLVE button from the
[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).

This is also implemented as the `Order0` method for `find_zero`.

The inital values can be specified as a pair of two values, as in
`(a,b)` or `[a,b]`, or as a single value, in which case a value of `b`
is computed, possibly from `fb`.  The basic idea is to follow the
secant method to convergence unless:

* a bracket is found, in which case bisection is used;

* the secant method is not converging, in which case a few steps of a
  quadratic method are used to see if that improves matters.

Convergence occurs when `f(m) == 0`, there is a sign change between
`m` and an adjacent floating point value, or `f(m) <= 2^3*eps(m)`.

A value of `NaN` is returned if the algorithm takes too many steps
before identifying a zero.

# Examples

```julia
Roots.dfree(x -> x^5 - x - 1, 1.0)
```

"""
function dfree(f, xs)

    if length(xs) == 1
        a = float(xs[1])
        fa = f(a)

        h = eps(one(a))^(1/3)
        da = h*oneunit(a) + abs(a)*h^2 # adjust for if eps(a) > h
        b = float(a + da)
        fb = f(b)
    else
        a, b = promote(float(xs[1]), float(xs[2]))
        fa, fb = f(a), f(b)
    end

    nan = (0*a)/(0*a) # try to preserve type
    cnt, MAXCNT = 0, 5 * ceil(Int, -log(eps(one(a))))  # must be higher for BigFloat
    MAXQUAD = 3

    if abs(fa) > abs(fb)
        a,fa,b,fb=b,fb,a,fa
    end

    # we keep a, b, fa, fb, gamma, fgamma
    quad_ctr = 0
    while !iszero(fb)
        cnt += 1

        if sign(fa) * sign(fb) < 0
            return bisection(f, a, b)
        end

        # take a secant step
        gamma =  float(b - (b-a) * fb / (fb - fa))
        # modify if gamma is too small or too big
        if iszero(abs(gamma-b))
            gamma = b + 1/1000 * abs(b-a)  # too small
        elseif abs(gamma-b)  >= 100 * abs(b-a)
            gamma = b + sign(gamma-b) * 100 * abs(b-a)  ## too big
        end
        fgamma = f(gamma)

        # change sign
        if sign(fgamma) * sign(fb) < 0
            return bisection(f, gamma, b)
        end

        # decreasing
        if abs(fgamma) < abs(fb)
            a,fa, b,fb = b, fb, gamma, fgamma
            quad_ctr = 0
            cnt < MAXCNT && continue
        end

        gamma = float(quad_vertex(a,fa,b,fb,gamma,fgamma))
        fgamma = f(gamma)
        # decreasing now?
        if abs(fgamma) < abs(fb)
            a,fa, b,fb = b, fb, gamma, fgamma
            quad_ctr = 0
            cnt < MAXCNT && continue
        end

        quad_ctr += 1
        if (quad_ctr > MAXQUAD) || (cnt > MAXCNT) || iszero(gamma - b)  || isnan(gamma)
            bprev, bnext = prevfloat(b), nextfloat(b)
            fbprev, fbnext = f(bprev), f(bnext)
            sign(fb) * sign(fbprev) < 0 && return b
            sign(fb) * sign(fbnext) < 0 && return b
            for (u,fu) in ((b,fb), (bprev, fbprev), (bnext, fbnext))
                abs(fu)/oneunit(fu) <= 2^3*eps(u/oneunit(u)) && return u
            end
            return nan # Failed.
        end

        if abs(fgamma) < abs(fb)
            b,fb, a,fa = gamma, fgamma, b, fb
        else
            a, fa = gamma, fgamma
        end

    end
    b

end


text:
using BayesianDataFusion
using Base.Test

if nprocs() <= 1
  addprocs(2)
  @everywhere using BayesianDataFusion
end

#W = sprand(50, 10, 0.1)
X = sprand(50, 20, 0.1)

A = randn(50, 2);
B = randn(10, 2);

W = sparse(A * B')

rd = RelationData(W, class_cut = 0.5, feat1 = X)
assignToTest!(rd.relations[1], 50)

result = macau(rd, burnin = 10, psamples = 10, num_latent=5, latent_pids=workers())
@test result["latent_multi_threading"] == true

result0 = macau(rd, burnin = 10, psamples = 10, num_latent=5, latent_pids=Int[])
@test result0["latent_multi_threading"] == false


text:
using DECAES

# Silence console output
DECAES.redirect_to_devnull() do
    # Run code to be precompiled which is not included in test suite coverage
    main(["--help"])

    # Run test suite, excluding Matlab MWI Toolbox compatibility tests
    ENV["DECAES_RUN_MWI_TOOLBOX_TESTS"] = "0"
    include(joinpath(pkgdir(DECAES), "test", "runtests.jl"))
end


text:
module PenalizedRegression

using Compat
using Reexport, Compat.LinearAlgebra
@reexport using LearnBase, LossFunctions, PenaltyFunctions, LearningStrategies

import LearnBase: nobs
import LearningStrategies: setup!, update!, hook, finished, cleanup!, learn!

#-----------------------------------------------------------------------# Abstract types
abstract type AbstractMLModel end
function Base.show(io::IO, o::AbstractMLModel)
    println(io, typeof(o))
    println(io, "  > β0: ", o.β0)
    print(io,   "  > β : ", o.β')
end
_modelsetup!(o::AbstractMLModel, data) = (o.β = zeros(npredictors(data)))

abstract type MLAlgorithm <: LearningStrategy end
Base.show(io::IO, o::MLAlgorithm) = print(io, name(o, false, false))
function name(o, withmodule = false, withparams = true)
    s = string(typeof(o))
    if !withmodule
        s = replace(s, r"([a-zA-Z]*\.)" => "")  # remove text that ends in period
    end
    if !withparams
        s = replace(s, r"\{(.*)" => "")  # remove "{" to the end of the string
    end
    s
end

abstract type ClosedFormSolution <: MLAlgorithm end
finished(strat::ClosedFormSolution, model) = true

#-----------------------------------------------------------------------# includes
include("data.jl")
include("models/regerm.jl")
include("algorithms/gradientbased.jl")
# include("models/linreg.jl")
include("spec.jl")

end # module


text:
# Autogenerated wrapper script for Arb_jll for aarch64-apple-darwin
export libarb

using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Arb")
JLLWrappers.@declare_library_product(libarb, "@rpath/libarb-2.12.1.dylib")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libarb,
        "lib/libarb.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()