text:
using StochasticDiffEq

## Set general parameters ##

## Define the parameters for the experiments:

# Path parameters:
path_params = (
    v̄ = 0.25,   # the average velocity of the path
    γ = 10.0,   # the convergence rate with which the momentary velocity relaxes back to the `v̄`
    aₐ= 0.25,   # the scale of acceleration in the angle component
    aᵥ= 0.1     # the scale of acceleration in the velocity component
)

# Grid parameters:
grid_params = (
    r = 13/45*0.1,      # aspect ratio of the grid
    xscale=0.1,
    yscale=0.1*85/90
)

# Plotting domain:
domain = ((0,0),(grid_params.xscale,grid_params.yscale))    # set the box domain containing the hex grid

# grid cell centers are located on a hexagonal grid (here 6 centers on a circle around the 1st)
gridcell_centers = [
    [cos(2π*4/6)*grid_params.r+0.5*grid_params.xscale,sin(2π*4/6)*grid_params.r+0.5*grid_params.yscale],
    [                          0.5*grid_params.xscale,                          0.5*grid_params.yscale],
    [cos(2π*1/6)*grid_params.r+0.5*grid_params.xscale,sin(2π*1/6)*grid_params.r+0.5*grid_params.yscale],
]

λ = 50.0                                #[Hz] rate at which population spikes are emitted
λ_background = 10.0                     #[Hz] rate at which population spikes are emitted
t_jitter = 5.0e-3
path_trange = (0.0, 0.2)              # duration of generated path
r = 1.5*grid_params.r                   # distance of path start-point from center
v_opt = 2r/(path_trange[2]-path_trange[1])  # velocity of the generated path (set such that path is symmetric around middle)
trials = 500                            # number of trials over which to average the response probability

## Define helper functions

"""
generate_stochastic_path(trange, params, u₀; dt=1e-3)

Generate one paths in the time interval `trange` with parameters specified by `param` and an initial condition `u₀` according to an sODE.
`params` must contain:
`v̄`    the average velocity of the path
`γ`    the convergence rate with which the momentary velocity relaxes back to the `v̄`
`aₐ`   the scale of acceleration in the angle component
`aᵥ`   the scale of acceleration in the velocity component

Return a JuliaDiffEq `solution` object.
"""
function generate_stochastic_path(trange, params, u₀, args...; kwargs...)
    """
    Drift term of the movement.
    """
    function f(du, u, p, t)
        du[1] = cos(2π*u[3])*u[4]
        du[2] = sin(2π*u[3])*u[4]
        du[3] = 0.0    # the angle doesn't change deterministically
        du[4] = p.γ*(p.v̄-u[4])# model smooth changes in velocity
    end

    """
    Diffusion term of the movement.
    """
    function g(du, u,p,t)
        du[1] = 0.0     # x positions don't change randomly
        du[2] = 0.0     # y positions don't change randomly
        du[3] = p.aₐ      # model fluctuations in angle
        du[4] = p.aᵥ      # model fluctuations in velocity
    end

    prob = SDEProblem(f, g, u₀, trange, params)
    sol = solve(prob, SRIW1(), args...; kwargs...)

    return (t, args...; idxs=1:2, kwargs...) -> trange[1] ≤ t ≤ trange[2] ? sol(t, args...; idxs, kwargs...) : missing
end

"""
generate_straight_path(trange, α, v, x₀; kwargs...)

Generate one straight-line path in the time interval `trange` with angle `α` (c.c.w from positive x-axis, in radians), velocity `v` and starting position `x₀`.
Returns a function `t -> [x(t),y(t)]`.
"""
generate_straight_path(trange, α, v, x₀) = ((t)-> trange[1] ≤ t ≤ trange[2] ? x₀ .+ [cos(α),sin(α)] .* (v*(t-trange[1])) : missing)

"""
generate_u₀(params,domain)

draw an initial condition from the bounded `domain` according to the specified `params`.
"""
generate_u₀(params,domain) = [rand()*(domain[2][1]-domain[1][1])+domain[1][1], rand()*(domain[2][2]-domain[1][2])+domain[1][2], rand(), randn()*params.aᵥ/√(2*params.γ)+params.v̄]

"""generate_path_spikes(path, populations, λ; group_size=ones(Int,length(rfs)), dt=0.0, background_rate=0.0, t_jitter=0.0)

Generates outputs resembling the activity of multiple homogeneous neuron `populations`.
`populations` must be an iterable of objects with fields:
- `rf` for the population's receptive field
- `neurons` for the population's neurons.

Emits volleys of spikes at a constant rate `λ` and encodes the receptive field responses
`populations[i].rf(input(t))` to the input at those times in the magnitudes of the volleys.
The spikes are jittered by a random time-lag uniformly drawn from `[0,t_jitter]`.
"""
function generate_path_spikes(trange, path, populations, λ; dt=0.0, background_rate=0.0, t_jitter=0.0)
    T = (trange[end]-trange[1])

    all_spikes = Event{Float64}[]
    for population ∈ populations
        N = rand(Poisson(λ*T))
        master_spikes = rand(N) .* T .+ trange[1]

        p_accept = population.rf.(path.(master_spikes))

        for neuron in population.neurons
            background_spikes = rand(rand(Poisson(background_rate*T))) .* T .+ trange[1]
            t_accepted = master_spikes[rand(length(master_spikes)) .≤ p_accept]
            spikes = sort!([t_accepted .+ t_jitter .* rand(length(t_accepted)); background_spikes])
            dels = Int[]
            for (k,t) ∈ enumerate(spikes[2:end])
                if t ≤ spikes[k] + dt
                    push!(dels, k+1)
                end
            end
            deleteat!(spikes, dels)
            append!(all_spikes, Event.(:input_spikes, 0.0, spikes, Ref(neuron)))
        end
    end

    return all_spikes
end

"""run_path!(net, trange, path, populations, λ; extra_events=[], kwargs...)

Simulates one traversal of the `path` by the given `net` over a time-interval `trange`.
Spikes are generated for each population in `populations` in volleys at fixed rate `λ`.
"""
function run_path!(net, trange, path, populations, λ; extra_events=Event{Float64}[], kwargs...)
    DPC.reset!(net)

    # convert path into spike trains and convert the spike trains into Events
    s = generate_path_spikes(trange, path, populations, λ; kwargs...)
    events = sort!([s; extra_events])

    # simulate the neuron's response to the input
    logger=simulate!(net, events)

    return s,logger
end


text:
module SymbolicRegression

# Types
export Population,
    PopMember,
    HallOfFame,
    Options,
    Node,

    #Functions:
    EquationSearch,
    SRCycle,
    calculateParetoFrontier,
    countNodes,
    copyNode,
    printTree,
    stringTree,
    evalTreeArray,
    differentiableEvalTreeArray,
    node_to_symbolic,
    symbolic_to_node,
    custom_simplify,
    simplifyWithSymbolicUtils,
    combineOperators,
    genRandomTree,

    #Operators:
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    log1p_abs,
    acosh_abs,
    sqrt_abs,
    neg,
    greater,
    relu,
    logical_or,
    logical_and,

    # special operators
    gamma,
    erf,
    erfc,
    atanh_clip

using Distributed
import JSON3
using Printf: @printf, @sprintf
using Pkg
using Random: seed!, shuffle!
using FromFile
using Reexport
@reexport using LossFunctions

@from "Core.jl" import CONST_TYPE, MAX_DEGREE, BATCH_DIM, FEATURE_DIM, RecordType, Dataset, Node, copyNode, Options, plus, sub, mult, square, cube, pow, div, log_abs, log2_abs, log10_abs, log1p_abs, sqrt_abs, acosh_abs, neg, greater, greater, relu, logical_or, logical_and, gamma, erf, erfc, atanh_clip, SRConcurrency, SRSerial, SRThreaded, SRDistributed
@from "Utils.jl" import debug, debug_inline, is_anonymous_function, recursive_merge, next_worker, @sr_spawner
@from "EquationUtils.jl" import countNodes, printTree, stringTree
@from "EvaluateEquation.jl" import evalTreeArray, differentiableEvalTreeArray
@from "CheckConstraints.jl" import check_constraints
@from "MutationFunctions.jl" import genRandomTree
@from "LossFunctions.jl" import EvalLoss, Loss, scoreFunc
@from "PopMember.jl" import PopMember, copyPopMember
@from "Population.jl" import Population, bestSubPop, record_population
@from "HallOfFame.jl" import HallOfFame, calculateParetoFrontier, string_dominating_pareto_curve
@from "SingleIteration.jl" import SRCycle, OptimizeAndSimplifyPopulation
@from "InterfaceSymbolicUtils.jl" import node_to_symbolic, symbolic_to_node
@from "CustomSymbolicUtilsSimplification.jl" import custom_simplify
@from "SimplifyEquation.jl" import simplifyWithSymbolicUtils, combineOperators, simplifyTree
@from "ProgressBars.jl" import ProgressBar, set_multiline_postfix
@from "Recorder.jl" import @recorder, find_iteration_from_record

include("Configure.jl")
include("Deprecates.jl")

"""
    EquationSearch(X, y[; kws...])

Perform a distributed equation search for functions `f_i` which
describe the mapping `f_i(X[:, j]) ≈ y[i, j]`. Options are
configured using SymbolicRegression.Options(...),
which should be passed as a keyword argument to options.
One can turn off parallelism with `numprocs=0`,
which is useful for debugging and profiling.

# Arguments
- `X::AbstractMatrix{T}`:  The input dataset to predict `y` from.
    The first dimension is features, the second dimension is rows.
- `y::Union{AbstractMatrix{T}, AbstractVector{T}}`: The values to predict. The first dimension
    is the output feature to predict with each equation, and the
    second dimension is rows.
- `niterations::Int=10`: The number of iterations to perform the search.
    More iterations will improve the results.
- `weights::Union{AbstractMatrix{T}, AbstractVector{T}, Nothing}=nothing`: Optionally
    weight the loss for each `y` by this value (same shape as `y`).
- `varMap::Union{Array{String, 1}, Nothing}=nothing`: The names
    of each feature in `X`, which will be used during printing of equations.
- `options::Options=Options()`: The options for the search, such as
    which operators to use, evolution hyperparameters, etc.
- `numprocs::Union{Int, Nothing}=nothing`:  The number of processes to use,
    if you want `EquationSearch` to set this up automatically. By default
    this will be `4`, but can be any number (you should pick a number <=
    the number of cores available).
- `procs::Union{Array{Int, 1}, Nothing}=nothing`: If you have set up
    a distributed run manually with `procs = addprocs()` and `@everywhere`,
    pass the `procs` to this keyword argument.
- `runtests::Bool=true`: Whether to run (quick) tests before starting the
    search, to see if there will be any problems during the equation search
    related to the host environment.

# Returns
- `hallOfFame::HallOfFame`: The best equations seen during the search.
    hallOfFame.members gives an array of `PopMember` objects, which
    have their tree (equation) stored in `.tree`. Their score (loss)
    is given in `.score`. The array of `PopMember` objects
    is enumerated by size from `1` to `options.maxsize`.
"""
function EquationSearch(X::AbstractMatrix{T}, y::AbstractMatrix{T};
        niterations::Int=10,
        weights::Union{AbstractMatrix{T}, AbstractVector{T}, Nothing}=nothing,
        varMap::Union{Array{String, 1}, Nothing}=nothing,
        options::Options=Options(),
        numprocs::Union{Int, Nothing}=nothing,
        procs::Union{Array{Int, 1}, Nothing}=nothing,
        multithreading::Bool=false,
        runtests::Bool=true
       ) where {T<:Real}

    nout = size(y, FEATURE_DIM)
    if weights != nothing
        weights = reshape(weights, size(y))
    end
    datasets = [Dataset(X, y[j, :],
                    weights=(weights == nothing ? weights : weights[j, :]),
                    varMap=varMap)
                for j=1:nout]

    return EquationSearch(datasets;
        niterations=niterations, options=options,
        numprocs=numprocs, procs=procs, multithreading=multithreading,
        runtests=runtests)
end

function EquationSearch(X::AbstractMatrix{T1}, y::AbstractMatrix{T2}; kw...) where {T1<:Real,T2<:Real}
    U = promote_type(T1, T2)
    EquationSearch(convert(AbstractMatrix{U}, X), convert(AbstractMatrix{U}, y); kw...)
end

function EquationSearch(X::AbstractMatrix{T1}, y::AbstractVector{T2}; kw...) where {T1<:Real,T2<:Real}
    EquationSearch(X, reshape(y, (1, size(y, 1))); kw...)
end

function EquationSearch(datasets::Array{Dataset{T}, 1};
        niterations::Int=10,
        options::Options=Options(),
        numprocs::Union{Int, Nothing}=nothing,
        procs::Union{Array{Int, 1}, Nothing}=nothing,
        multithreading::Bool=false,
        runtests::Bool=true
       ) where {T<:Real}

    noprocs = (procs == nothing && numprocs == 0)
    someprocs = !noprocs

    concurrency = if multithreading
        @assert procs == nothing && numprocs in [0, nothing]
        SRThreaded()
    elseif someprocs
        SRDistributed()
    else #noprocs, multithreading=false
        SRSerial()
    end

    return _EquationSearch(concurrency, datasets;
        niterations=niterations, options=options,
        numprocs=numprocs, procs=procs,
        runtests=runtests)
end

function _EquationSearch(::ConcurrencyType, datasets::Array{Dataset{T}, 1};
        niterations::Int=10,
        options::Options=Options(),
        numprocs::Union{Int, Nothing}=nothing,
        procs::Union{Array{Int, 1}, Nothing}=nothing,
        runtests::Bool=true,
       ) where {T<:Real,ConcurrencyType<:SRConcurrency}

    example_dataset = datasets[1]
    nout = size(datasets, 1)

    if runtests
        testOptionConfiguration(T, options)
        # Testing the first output variable is the same:
        testDatasetConfiguration(example_dataset, options)
    end

    if example_dataset.weighted
        avgys = [sum(dataset.y .* dataset.weights) / sum(dataset.weights)
                for dataset in datasets]
        baselineMSEs = [Loss(dataset.y, ones(T, dataset.n) .* avgy, dataset.weights, options)
                       for dataset in datasets, avgy in avgys]
    else
        avgys = [sum(dataset.y)/dataset.n for dataset in datasets]
        baselineMSEs = [Loss(dataset.y, ones(T, dataset.n) .* avgy, options)
                       for (dataset, avgy) in zip(datasets, avgys)]
    end

    if options.seed !== nothing
        seed!(options.seed)
    end
    # Start a population on every process
    #    Store the population, hall of fame
    allPopsType = if ConcurrencyType == SRSerial
        Tuple{Population,HallOfFame,RecordType}
    elseif ConcurrencyType == SRDistributed
        Future
    else
        Task
    end

    allPops = [allPopsType[] for j=1:nout]
    init_pops = [allPopsType[] for j=1:nout]
    # Set up a channel to send finished populations back to head node
    if ConcurrencyType in [SRDistributed, SRThreaded]
        if ConcurrencyType == SRDistributed
            channels = [[RemoteChannel(1) for i=1:options.npopulations] for j=1:nout]
        else
            channels = [[Channel(1) for i=1:options.npopulations] for j=1:nout]
        end
        tasks = [Task[] for j=1:nout]
    end

    # These initial populations are discarded:
    bestSubPops = [[Population(datasets[j], baselineMSEs[j], npop=1, options=options, nfeatures=datasets[j].nfeatures)
                    for i=1:options.npopulations]
                    for j=1:nout]
    hallOfFame = [HallOfFame(options) for j=1:nout]
    actualMaxsize = options.maxsize + MAX_DEGREE
    frequencyComplexities = [ones(T, actualMaxsize) for i=1:nout]
    curmaxsizes = [3 for j=1:nout]
    record = RecordType("options"=>"$(options)")

    if options.warmupMaxsizeBy == 0f0
        curmaxsizes = [options.maxsize for j=1:nout]
    end

    we_created_procs = false
    ##########################################################################
    ### Distributed code:
    ##########################################################################
    if ConcurrencyType == SRDistributed
        if numprocs == nothing && procs == nothing
            numprocs = 4
            procs = addprocs(4)
            we_created_procs = true
        elseif numprocs == nothing
            numprocs = length(procs)
        elseif procs == nothing
            procs = addprocs(numprocs)
            we_created_procs = true
        end

        if we_created_procs
            project_path = splitdir(Pkg.project().path)[1]
            activate_env_on_workers(procs, project_path, options)
            import_module_on_workers(procs, @__FILE__, options)
        end
        move_functions_to_workers(procs, options, example_dataset)
        if runtests
            test_module_on_workers(procs, options)
        end

        if runtests
            test_entire_pipeline(procs, example_dataset, options)
        end
    end
    # Get the next worker process to give a job:
    worker_assignment = Dict{Tuple{Int,Int}, Int}()

    for j=1:nout
        for i=1:options.npopulations
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end
            new_pop = @sr_spawner ConcurrencyType worker_idx (
                Population(datasets[j], baselineMSEs[j], npop=options.npop,
                           nlength=3, options=options, nfeatures=datasets[j].nfeatures),
                HallOfFame(options),
                RecordType()
            )
            push!(init_pops[j], new_pop)
        end
    end
    # 2. Start the cycle on every process:
    for j=1:nout
        dataset = datasets[j]
        baselineMSE = baselineMSEs[j]
        frequencyComplexity = frequencyComplexities[j]
        curmaxsize = curmaxsizes[j]
        for i=1:options.npopulations
            @recorder record["out$(j)_pop$(i)"] = RecordType()
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end

            # TODO - why is this needed??
            # Multi-threaded doesn't like to fetch within a new task:
            updated_pop = @sr_spawner ConcurrencyType worker_idx let
                in_pop = if ConcurrencyType in [SRDistributed, SRThreaded]
                    fetch(init_pops[j][i])[1]
                else
                    init_pops[j][i][1]
                end

                cur_record = RecordType()
                @recorder cur_record["out$(j)_pop$(i)"] = RecordType("iteration0"=>record_population(in_pop, options))
                tmp_pop, tmp_best_seen = SRCycle(dataset, baselineMSE, in_pop,
                                                 options.ncyclesperiteration, curmaxsize,
                                                 copy(frequencyComplexity)/sum(frequencyComplexity),
                                                 verbosity=options.verbosity, options=options,
                                                 record=cur_record)
                tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize, cur_record)
                if options.batching
                    for i_member=1:(options.maxsize + MAX_DEGREE)
                        tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                    end
                end
                (tmp_pop, tmp_best_seen, cur_record)
            end
            push!(allPops[j], updated_pop)
        end
    end

    debug(options.verbosity > 0 || options.progress, "Started!")
    total_cycles = options.npopulations * niterations
    cycles_remaining = [total_cycles for j=1:nout]
    sum_cycle_remaining = sum(cycles_remaining)
    if options.progress && nout == 1
        #TODO: need to iterate this on the max cycles remaining!
        progress_bar = ProgressBar(1:sum_cycle_remaining;
                                   width=options.terminal_width)
        cur_cycle = nothing
        cur_state = nothing
    end

    last_print_time = time()
    num_equations = 0.0
    print_every_n_seconds = 5
    equation_speed = Float32[]

    if ConcurrencyType in [SRDistributed, SRThreaded]
        for j=1:nout
            for i=1:options.npopulations
                # Start listening for each population to finish:
                t = @async put!(channels[j][i], fetch(allPops[j][i]))
                push!(tasks[j], t)
            end
        end
    end

    # Randomly order which order to check populations:
    # This is done so that we do work on all nout equally.
    all_idx = [(j, i) for j=1:nout for i=1:options.npopulations]
    shuffle!(all_idx)
    kappa = 0
    head_node_time = Dict("occupied"=>0, "start"=>time())
    while sum(cycles_remaining) > 0
        kappa += 1
        if kappa > options.npopulations * nout
            kappa = 1
        end
        # nout, npopulations:
        j, i = all_idx[kappa]

        # Check if error on population:
        if ConcurrencyType in [SRDistributed, SRThreaded]
            if istaskfailed(tasks[j][i])
                fetch(tasks[j][i])
                error("Task failed for population")
            end
        end
        # Non-blocking check if a population is ready:
        population_ready = ConcurrencyType in [SRDistributed, SRThreaded] ? isready(channels[j][i]) : true
        # Don't start more if this output has finished its cycles:
        # TODO - this might skip extra cycles?
        population_ready &= (cycles_remaining[j] > 0)
        if population_ready
            head_node_start_work = time()
            # Take the fetch operation from the channel since its ready
            (cur_pop, best_seen, cur_record) = ConcurrencyType in [SRDistributed, SRThreaded] ? take!(channels[j][i]) : allPops[j][i]
            cur_pop::Population
            best_seen::HallOfFame
            cur_record::RecordType
            bestSubPops[j][i] = bestSubPop(cur_pop, topn=options.topn)
            @recorder record = recursive_merge(record, cur_record)

            dataset = datasets[j]
            baselineMSE = baselineMSEs[j]
            curmaxsize = curmaxsizes[j]

            #Try normal copy...
            bestPops = Population([member for pop in bestSubPops[j] for member in pop.members])

            for (i_member, member) in enumerate(Iterators.flatten((cur_pop.members, best_seen.members[best_seen.exists])))
                part_of_cur_pop = i_member <= length(cur_pop.members)
                size = countNodes(member.tree)
                if part_of_cur_pop
                    frequencyComplexities[j][size] += 1
                end
                actualMaxsize = options.maxsize + MAX_DEGREE
                if size < actualMaxsize && all([member.score < hallOfFame[j].members[size2].score for size2=1:size])
                    hallOfFame[j].members[size] = copyPopMember(member)
                    hallOfFame[j].exists[size] = true
                end
            end

            # Dominating pareto curve - must be better than all simpler equations
            dominating = calculateParetoFrontier(dataset, hallOfFame[j], options)
            hofFile = options.hofFile
            if nout > 1
                hofFile = hofFile * ".out$j"
            end
            open(hofFile, "w") do io
                println(io,"Complexity|MSE|Equation")
                for member in dominating
                    println(io, "$(countNodes(member.tree))|$(member.score)|$(stringTree(member.tree, options, varMap=dataset.varMap))")
                end
            end
            cp(hofFile, hofFile*".bkup", force=true)

            # Try normal copy otherwise.
            if options.migration
                for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplaced))
                    to_copy = rand(1:size(bestPops.members, 1))
                    cur_pop.members[k] = PopMember(
                        copyNode(bestPops.members[to_copy].tree),
                        bestPops.members[to_copy].score)
                end
            end

            if options.hofMigration && size(dominating, 1) > 0
                for k in rand(1:options.npop, round(Int, options.npop*options.fractionReplacedHof))
                    # Copy in case one gets used twice
                    to_copy = rand(1:size(dominating, 1))
                    cur_pop.members[k] = PopMember(
                       copyNode(dominating[to_copy].tree), dominating[to_copy].score
                    )
                end
            end

            cycles_remaining[j] -= 1
            if cycles_remaining[j] == 0
                break
            end
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end
            @recorder begin
                key = "out$(j)_pop$(i)"
                iteration = find_iteration_from_record(key, record) + 1
            end

            allPops[j][i] = @sr_spawner ConcurrencyType worker_idx let
                cur_record = RecordType()
                @recorder cur_record[key] = RecordType("iteration$(iteration)"=>record_population(cur_pop, options))
                tmp_pop, tmp_best_seen = SRCycle(
                    dataset, baselineMSE, cur_pop, options.ncyclesperiteration,
                    curmaxsize, copy(frequencyComplexities[j])/sum(frequencyComplexities[j]),
                    verbosity=options.verbosity, options=options, record=cur_record)
                tmp_pop = OptimizeAndSimplifyPopulation(dataset, baselineMSE, tmp_pop, options, curmaxsize, cur_record)
                if options.batching
                    for i_member=1:(options.maxsize + MAX_DEGREE)
                        tmp_best_seen.members[i_member].score = scoreFunc(dataset, baselineMSE, tmp_best_seen.members[i_member].tree, options)
                    end
                end
                (tmp_pop, tmp_best_seen, cur_record)
            end
            if ConcurrencyType in [SRDistributed, SRThreaded]
                tasks[j][i] = @async put!(channels[j][i], fetch(allPops[j][i]))
            end

            cycles_elapsed = total_cycles - cycles_remaining[j]
            if options.warmupMaxsizeBy > 0
                fraction_elapsed = 1f0 * cycles_elapsed / total_cycles
                if fraction_elapsed > options.warmupMaxsizeBy
                    curmaxsizes[j] = options.maxsize
                else
                    curmaxsizes[j] = 3 + floor(Int, (options.maxsize - 3) * fraction_elapsed / options.warmupMaxsizeBy)
                end
            end
            num_equations += options.ncyclesperiteration * options.npop / 10.0

            if options.progress && nout == 1
                # set_postfix(iter, Equations=)
                equation_strings = string_dominating_pareto_curve(hallOfFame[j], baselineMSE,
                                                                  datasets[j], options,
                                                                  avgys[j])
                load_string = @sprintf("Head worker occupation: %.1f", 100*head_node_time["occupied"]/(time() - head_node_time["start"])) * "%\n"
                equation_strings = load_string * equation_strings
                set_multiline_postfix(progress_bar, equation_strings)
                if cur_cycle == nothing
                    (cur_cycle, cur_state) = iterate(progress_bar)
                else
                    (cur_cycle, cur_state) = iterate(progress_bar, cur_state)
                end
                sum_cycle_remaining = sum(cycles_remaining)
            end
            head_node_end_work = time()
            head_node_time["occupied"] += (head_node_end_work - head_node_start_work)
        end
        sleep(1e-6)

        ################################################################
        ## Printing code
        elapsed = time() - last_print_time
        #Update if time has passed, and some new equations generated.
        if elapsed > print_every_n_seconds && num_equations > 0.0
            # Dominating pareto curve - must be better than all simpler equations
            current_speed = num_equations/elapsed
            average_over_m_measurements = 10 #for print_every...=5, this gives 50 second running average
            push!(equation_speed, current_speed)
            if length(equation_speed) > average_over_m_measurements
                deleteat!(equation_speed, 1)
            end
            if (options.verbosity > 0) || (options.progress && nout > 1)
                @printf("\n")
                average_speed = sum(equation_speed)/length(equation_speed)
                @printf("Cycles per second: %.3e\n", round(average_speed, sigdigits=3))
                @printf("Head worker occupation: %.1f%%\n", 100 * head_node_time["occupied"]/(time() - head_node_time["start"]))
                cycles_elapsed = total_cycles * nout - sum(cycles_remaining)
                @printf("Progress: %d / %d total iterations (%.3f%%)\n",
                        cycles_elapsed, total_cycles * nout,
                        100.0*cycles_elapsed/total_cycles/nout)

                @printf("==============================\n")
                for j=1:nout
                    if nout > 1
                        @printf("Best equations for output %d\n", j)
                    end
                    equation_strings = string_dominating_pareto_curve(hallOfFame[j], baselineMSEs[j],
                                                                      datasets[j], options,
                                                                      avgys[j])
                    print(equation_strings)
                    @printf("==============================\n")
                end
            end
            last_print_time = time()
            num_equations = 0.0
        end
        ################################################################
    end
    if we_created_procs
        rmprocs(procs)
    end
    ##########################################################################
    ### Distributed code^
    ##########################################################################

    @recorder begin
        open(options.recorder_file, "w") do io
            JSON3.write(io, record, allow_inf = true)
        end
    end

    if nout == 1
        return hallOfFame[1]
    else
        return hallOfFame
    end
end

end #module SR

ext: jl


text:
# These tests are to be skipped in OSX builds, see ./interactive.jl for more info.

@testset "TravisCI" begin
    write(stdin.buffer, "\n")
    p = interactive(TravisCI)
    @test p.src == joinpath(DEFAULTS_DIR, "travis.yml")
    write(stdin.buffer, "$test_file\n")
    p = interactive(TravisCI)
    @test p.src == test_file
    write(stdin.buffer, "none\n")
    p = interactive(TravisCI)
    @test p.src === nothing
    write(stdin.buffer, "$fake_path\n")
    @test_throws ArgumentError interactive(TravisCI)
    println()
end

@testset "AppVeyor" begin
    write(stdin.buffer, "\n")
    p = interactive(AppVeyor)
    @test p.src == joinpath(DEFAULTS_DIR, "appveyor.yml")
    write(stdin.buffer, "$test_file\n")
    p = interactive(AppVeyor)
    @test p.src == test_file
    write(stdin.buffer, "none\n")
    p = interactive(AppVeyor)
    @test p.src === nothing
    write(stdin.buffer, "$fake_path\n")
    @test_throws ArgumentError interactive(AppVeyor)
    println()
end

@testset "GitLabCI" begin
    write(stdin.buffer, "\n\n")
    p = interactive(GitLabCI)
    @test p.src == joinpath(DEFAULTS_DIR, "gitlab-ci.yml")
    @test p.view == Dict("GITLABCOVERAGE" => true)
    write(stdin.buffer, "$test_file\nno\n")
    p = interactive(GitLabCI)
    @test p.src == test_file
    @test p.view == Dict("GITLABCOVERAGE" => false)
    write(stdin.buffer, "none\n\n")
    p = interactive(GitLabCI)
    @test p.src === nothing
    write(stdin.buffer, "$fake_path\n\n")
    @test_throws ArgumentError interactive(GitLabCI)
    println()
end

@testset "Codecov" begin
    write(stdin.buffer, "\n")
    p = interactive(Codecov)
    @test p.src === nothing
    write(stdin.buffer, "$test_file\n")
    p = interactive(Codecov)
    @test p.src == test_file
    write(stdin.buffer, "none\n")
    p = interactive(Codecov)
    @test p.src === nothing
    write(stdin.buffer, "$fake_path\n")
    @test_throws ArgumentError interactive(Codecov)
    println()
end

@testset "Coveralls" begin
    write(stdin.buffer, "\n")
    p = interactive(Coveralls)
    @test p.src === nothing
    write(stdin.buffer, "$test_file\n")
    p = interactive(Coveralls)
    @test p.src == test_file
    write(stdin.buffer, "none\n")
    p = interactive(Coveralls)
    @test p.src === nothing
    write(stdin.buffer, "$fake_path\n")
    @test_throws ArgumentError interactive(Coveralls)
    println()
end

@testset "GitHubPages" begin
    write(stdin.buffer, "\n")
    p = interactive(GitHubPages)
    @test isempty(p.assets)
    write(stdin.buffer, "$test_file\n")
    p = interactive(GitHubPages)
    @test p.assets == [test_file]
    write(stdin.buffer, "$fake_path\n")
    @test_throws ArgumentError interactive(GitHubPages)
    println()
end


text:
using Test
using Plots

opts = TrajectoryOptimization.SolverOptions()
opts.verbose = true
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-6
opts.square_root = true
opts.outer_loop_update_type = :default
opts.constraint_tolerance_second_order_dual_update = sqrt(opts.constraint_tolerance)
opts.use_second_order_dual_update = true
opts.penalty_max = 1e8
opts.iterations_outerloop = 25
opts.live_plotting = false
model,obj = TrajectoryOptimization.Dynamics.dubinscar_parallelpark

solver = Solver(model, obj, integration=:rk4, N=100, opts=opts)
U0 = rand(solver.model.m,solver.N-1)

results, stats = TrajectoryOptimization.solve(solver,U0)

@show stats["iterations"]
@show stats["outer loop iterations"]
@show stats["c_max"][end]
@show max_violation(results) < opts.constraint_tolerance

plot(stats["c_max"],yscale=:log10,ylabel="c_max",label="first")

plot(to_array(results.U)',label="")
plot(to_array(results.X)',label="")


text:
function polycoefs(n::Int64)
    pc = typeof(n)[]
    if n < 0
        return pc
    end
    sgn = one(n)
    for k in n:-1:0
        push!(pc, sgn*binomial(n, k))
        sgn = -sgn
    end
    return pc
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using PeriodicTable, Unitful
import Unitful: u, g, cm, K, J, mol
using Test
using Base64

@test eltype(elements) == Element
@test length(elements) == 119 == length(collect(elements))

# test element lookup
O = elements[8]
F = elements[9]
@test O === elements["oxygen"] == elements[:O]
@test haskey(elements, 8) && haskey(elements, "oxygen") && haskey(elements, :O)
@test !haskey(elements, -8) && !haskey(elements, "ooxygen") && !haskey(elements, :Oops)
@test F === get(elements, 9, O) === get(elements, "oops", F) === get(elements, :F, O)
@test elements[8:9] == [O, F]
@test O.name == "Oxygen"
@test O.symbol == "O"
@test nfields(O) == 23

# cpk colors
@test O.cpk_hex == "#ff0d0d"
@test F.cpk_hex == "#90e050"

# Unitful units
H = elements[1]
@test unit(H.density) === g/cm^3
@test unit(H.boil) === K
@test unit(H.melt) === K
@test unit(H.molar_heat) === J/(mol*K)
@test unit(H.atomic_mass) === u

# iteration protocol
if VERSION < v"0.7-"
    @test start(elements) == 1
    @test next(elements, 1) == (elements[:H], 2)
    @test done(elements, length(elements)+1)
else
    @test iterate(elements) == (elements[:H], 2)
    @test iterate(elements, 4) == (elements[:Be], 5)
    @test iterate(elements, length(elements)+1) === nothing
end

# 2-argument show functions
@test repr(elements) == "Elements(…119 elements…)"
if VERSION < v"0.7-"
    @test repr([O, F]) == "PeriodicTable.Element[Element(Oxygen), Element(Fluorine)]"
else
    @test repr([O, F]) == "Element[Element(Oxygen), Element(Fluorine)]"
end

@test stringmime("text/plain", elements) == "Elements(…119 elements…):\nH                                                  He \nLi Be                               B  C  N  O  F  Ne \nNa Mg                               Al Si P  S  Cl Ar \nK  Ca Sc Ti V  Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr \nRb Sr Y  Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I  Xe \nCs Ba    Hf Ta W  Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn \nFr Ra    Rf Db Sg Bh Hs Mt Ds Rg Cn Nh Fl Mc Lv Ts Og \nUue                                                   \n      La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu    \n      Ac Th Pa U  Np Pu Am Cm Bk Cf Es Fm Md No Lr    \n"
@test stringmime("text/plain", O) == "Oxygen (O), number 8:\n        category: diatomic nonmetal\n     atomic mass: 15.999 u\n         density: 1.429 g/cm³\n   melting point: 54.36 K\n   boiling point: 90.188 K\n           phase: Gas\n          shells: [2, 6]\ne⁻-configuration: 1s² 2s² 2p⁴\n         summary: Oxygen is a chemical element with symbol O and atomic number 8. It is a member of the chalcogen group on the periodic table and is a highly reactive nonmetal and oxidizing agent that readily forms compounds (notably oxides) with most elements. By mass, oxygen is the third-most abundant element in the universe, after hydrogen and helium.\n   discovered by: Carl Wilhelm Scheele\n        named by: Antoine Lavoisier\n          source: https://en.wikipedia.org/wiki/Oxygen\n  spectral image: https://en.wikipedia.org/wiki/File:Oxygen_spectre.jpg\n"
@test stringmime("text/html", O) == "\nOxygen (O), number 8:\n
\n
category	diatomic nonmetal
\n
atomic mass	15.999 u
\n
density	1.429 g/cm³
\n
melting point	54.36 K
\n
boiling point	90.188 K
\n
phase	Gas
\n
shells	[2, 6]
\n
electron configuration	1s² 2s² 2p⁴
\n
summary	Oxygen is a chemical element with symbol O and atomic number 8. It is a member of the chalcogen group on the periodic table and is a highly reactive nonmetal and oxidizing agent that readily forms compounds (notably oxides) with most elements. By mass, oxygen is the third-most abundant element in the universe, after hydrogen and helium.
\n
discovered by	Carl Wilhelm Scheele
\n
named by	Antoine Lavoisier
\n
source	https://en.wikipedia.org/wiki/Oxygen
\n
\n\n"

@test_throws ErrorException O.name = "Issue21"
@test O.name == "Oxygen"

@test isless(elements[28], elements[29])
@test !isless(elements[88], elements[88])
@test !isless(elements[92], elements[90])
@test isequal(elements[38], elements[38])
@test !isequal(elements[38], elements[39])

@test elements[28] < elements[29]
@test ! (elements[88] < elements[88])
@test ! (elements[92] < elements[90])
@test elements[92] > elements[91]
@test !(elements[92] > elements[92])
@test !(elements[92] > elements[93])
@test elements[90] <= elements[91]
@test elements[91] <= elements[91]
@test !(elements[92] <= elements[91])
@test elements[38] == elements[38]
@test elements[38] ≠ elements[39]

# Ensure that the hashcode works in Dict{}
elmdict = Dict{Element,Int}( elements[z] => z for z in eachindex(elements))
@test length(elmdict) == 119
for z in eachindex(elements)
    @test haskey(elmdict, elements[z])
    @test elmdict[elements[z]] == z
end


text:
module BridgeSPDE
using Bridge
import Bridge: outer, b, a

using GaussianDistributions, Trajectories, DynamicIterators
using LinearAlgebra, FileIO, Random, Statistics, SparseArrays, Colors

using DynamicIterators: dub
import DynamicIterators: evolve, dyniterate

export gridderiv, gridlaplacian, OrnsteinUhlenbeck, b, riccatti
export downsample, boundary
export Euler
export sparsity

include("sparse.jl")
include("laplace.jl")
include("ornsteinuhlenbeck.jl")

include("sidrs.jl")
include("gaussian.jl")
end # module


text:
# code for Video 4.1:

julia> fname = "iris.csv"
julia> using DataFrames
julia> data = readtable(fname, separator = ',')

julia> describe(data)
sepal_length
Min      4.3
1st Qu.  5.1
Median   5.8
Mean     5.843333333333334
3rd Qu.  6.4
Max      7.9
NAs      0
NA%      0.0%

sepal_width
Min      2.0
1st Qu.  2.8
Median   3.0
Mean     3.054
3rd Qu.  3.3
Max      4.4
NAs      0
NA%      0.0%

petal_length
Min      1.0
1st Qu.  1.6
Median   4.35
Mean     3.758666666666666
3rd Qu.  5.1
Max      6.9
NAs      0
NA%      0.0%

petal_width
Min      0.1
1st Qu.  0.3
Median   1.3
Mean     1.1986666666666665
3rd Qu.  1.8
Max      2.5
NAs      0
NA%      0.0%

species
Length  150
Type    UTF8String
NAs     0
NA%     0.0%
Unique  3

julia> summarystats(data[:sepal_length])
Summary Stats:
Mean:         5.843333
Minimum:      4.300000
1st Quartile: 5.100000
Median:       5.800000
3rd Quartile: 6.400000
Maximum:      7.900000

julia> minimum(data[:sepal_length]) # 4.3
julia> maximum(data[:sepal_length]) # 7.9
julia> extrema(data[:sepal_length]) # (4.3,7.9)

julia> mean(data[:sepal_length]) # 5.843333333333334
julia> length(data[:sepal_length]) # 150
julia> median(data[:sepal_length]) # 5.8
julia> std(data[:sepal_length]) # 0.828066127977863
julia> var(data[:sepal_length]) # 0.6856935123042507 = std * std
julia> variation(data[:sepal_length]) # 0.14171125977944032
julia> mode(data[:sepal_length]) # 5.0
julia> sem(data[:sepal_length]) # 0.0676113162275986
julia> trimmean(data[:sepal_length], 0.1) # 5.8096296296296295
julia> iqr(data[:sepal_length]) # 1.3000000000000007

julia> quantile(data[:sepal_length])
5-element Array{Float64,1}:
 4.3
 5.1
 5.8
 6.4
 7.9

julia> hist(data[:sepal_length])
(4.0:0.5:8.0,[5,27,27,30,31,18,6,6])

julia> countmap(data[:species])
Dict{Union{DataArrays.NAtype,UTF8String},Int64} with 3 entries:
  "I. virginica"  => 50
  "I. versicolor" => 50
  "I. setosa"     => 50

julia> proportionmap(data[:species])
Dict{Union{DataArrays.NAtype,UTF8String},Float64} with 3 entries:
  "I. virginica"  => 0.3333333333333333
  "I. versicolor" => 0.3333333333333333
  "I. setosa"     => 0.3333333333333333

julia> skewness(data[:sepal_length]) # 0.31175305850229673
julia> kurtosis(data[:sepal_length]) # -0.5735679489249748
julia> round(skewness(data[:sepal_length]),3) # 0.312

julia> cor(data[:sepal_width], data[:sepal_length]) # -0.10936924995064931
julia> cor(data[:sepal_length], data[:petal_length])# 0.8717541573048717

julia> cov(data[:sepal_width], data[:sepal_length]) # -0.03926845637583893
julia> cov(data[:sepal_length], data[:petal_length]) # 1.2736823266219242

max_stars_repo_path: Section 4/4_1.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 4/4_1.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 4/4_1.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using FinancialToolbox, DualNumbers, FinancialFFT, FinancialMonteCarlo

A = 600.0;
N = 16;

method = CarrMadanMethod(A, N);
S0 = 100.0;
K = 100.0;
r = dual(0.02, 1.0);
# r = 0.02;
T = 1.0;
d = 0.01;
# sigma = dual(0.2, 1.0);
sigma = 0.2;
lam = 15.0;
mu1 = 0.03;
sigma1 = 0.02;
zero_rate = ZeroRate(r);

Model = MertonProcess(sigma, lam, mu1, sigma1, Underlying(S0, d));
#Model = BlackScholesProcess(sigma, Underlying(S0, d));

EUData = EuropeanOption(T, K)
Nsim = 10000;
Nstep = 30;
mc = MonteCarloConfiguration(Nsim, Nstep);

@time @show pricer(Model, zero_rate, method, EUData);
@time @show pricer(Model, zero_rate, mc, EUData);
typeof(Model) <: BlackScholesProcess ? @time(@show(blsprice(S0, K, r, T, sigma, d))) : nothing;


text:
function _create_dictionary_for_join_general(f, v, ::Val{T}) where T
    flag = false
    n = length(v)
    sz = max(1 + ((5 * n) >> 2), 16)
    sz = 1 << (8 * sizeof(sz) - leading_zeros(sz - 1))
    @assert 4 * sz >= 5 * n
    gslots = zeros(T, sz)
    szm1 = sz - 1
    maxprobe = 0
    @inbounds for i in 1:length(v)
        _fv = f(v[i])
        slotix = hash(_fv) & szm1 + 1
        probe = 0
        while true
            g_row = gslots[slotix]
            if g_row == 0
                gslots[slotix] = i
                break
            elseif isequal(_fv, f(v[g_row]))
                flag = true
                return gslots, maxprobe, sz, flag, 1
            end
            slotix = slotix & szm1 + 1
            probe += 1
            @assert probe < sz
        end
        if probe > maxprobe
            maxprobe = probe
        end
    end
    gslots, maxprobe, sz, flag, 1
end

function _create_dictionary_for_join_int(f, v, minval, rangelen, ::Val{T}) where T
    flag = false
    offset = 1 - minval
    n = length(v)
    sz = rangelen + 1
    gslots = zeros(T, sz)
    @inbounds for i in 1:length(v)
        _fv = f(v[i])
        if ismissing(_fv)
            slotix = sz
        else
            slotix = _fv + offset
        end
        if gslots[slotix] == 0
            gslots[slotix] = i
        else
            flag = true
            break
        end
    end
    gslots, minval, sz, flag, 2
end

function _create_dictionary_for_join(f, v, fl, vl, ::Val{T}) where T
    if f == _date_value∘identity && DataAPI.refpool(v) !== nothing
        minval = hp_minimum(DataAPI.refarray(v))
        maxval = hp_maximum(DataAPI.refarray(v))
        rangelen = maxval - minval + 1
        _create_dictionary_for_join_int(identity, DataAPI.refarray(v), minval, rangelen, Val(T))
    elseif nonmissingtype(return_type(f, v)) <: AbstractVector{<:Union{Missing, INTEGERS}} && nonmissingtype(return_type(fl, vl)) <: AbstractVector{<:Union{Missing, INTEGERS}}
        minval = hp_minimum(f, v)
        # if minval is missing all values are missing
        if ismissing(minval)
            return _create_dictionary_for_join_general(f, v, Val(T))
        end
        maxval::Integer = hp_maximum(f, v)
        rnglen = BigInt(maxval) - BigInt(minval) + 1
        o1 = false
        if rnglen < typemax(Int)
            o1 = true
            rangelen = Int(rnglen)
        end
        if o1 && maxval < typemax(Int) && (rangelen < 2.0*length(v))
            _create_dictionary_for_join_int(f, v, Int(minval), Int(rangelen), Val(T))
        else
            _create_dictionary_for_join_general(f, v, Val(T))
        end
    else
        _create_dictionary_for_join_general(f, v, Val(T))
    end
end

function _query_dictionary_for_join_general(f, r_v, fv, gslots, maxprobe, sz)
    szm1 = sz - 1
    slotix = hash(fv) & szm1 + 1
    probe = 0
    @inbounds while true
        g_row = gslots[slotix]
        g_row == 0 && break
        if isequal(fv, f(r_v[g_row]))
            return g_row
        end
        slotix = slotix & szm1 + 1
        probe += 1
        probe > maxprobe && break
    end
    0
end
function _query_dictionary_for_join_int(f, r_v, fv, gslots, minval, sz)
    offset = 1 - minval
    if ismissing(fv)
        slotix = sz
    else
        slotix = fv + offset
        !(slotix in 1:sz-1) && return 0
    end

    if slotix in 1:sz
        rowid = gslots[slotix]
        if rowid == 0
            return 0
        else
            return rowid
        end
    end
    0
end

function _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, x_l, x_r, sz, type; threads = true)
    if _fr == _date_value∘identity && DataAPI.refpool(x_r) !== nothing
        invp = DataAPI.invrefpool(x_r)
        @_threadsfor threads for i in 1:length(x_l)
            revmap_paval_ref = get(invp, DataAPI.unwrap(_fl(x_l[i])), missing)
            loc = _query_dictionary_for_join_int(_fr, x_r,revmap_paval_ref, dict, maxprob, sz)
            if loc == 0
                ranges[i] = 1:0
            else
                ranges[i] = loc:loc
            end
        end
    elseif return_type(_fr, x_r) <: AbstractVector{<:Union{Missing, INTEGERS}} && type == 2
        @_threadsfor threads for i in 1:length(x_l)
            loc = _query_dictionary_for_join_int(_fr, x_r, _fl(x_l[i]), dict, maxprob, sz)
            if loc == 0
                ranges[i] = 1:0
            else
                ranges[i] = loc:loc
            end
        end
    elseif type == 1
        @_threadsfor threads for i in 1:length(x_l)
            loc = _query_dictionary_for_join_general(_fr, x_r, _fl(x_l[i]), dict, maxprob, sz)
            if loc == 0
                ranges[i] = 1:0
            else
                ranges[i] = loc:loc
            end
        end
    end
end

function _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, makeunique, ::Val{T}; threads = true) where T
    oncols_left = onleft
    oncols_right = onright
    right_cols = setdiff(1:length(index(dsr)), oncols_right)
    if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
        throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
    end

    cols = Any[]
    for j in 1:length(oncols_left)
        if mapformats[1]
            fl = getformat(dsl, oncols_left[j])
        else
            fl = identity
        end
        if mapformats[2]
            fr = getformat(dsr, oncols_right[j])
        else
            fr = identity
        end
        push!(cols, Cat2Vec(_columns(dsl)[oncols_left[j]], _columns(dsr)[oncols_right[j]], fl, fr))
    end
    newds = Dataset(cols, :auto, copycols = false)
    a = _gather_groups(newds, :, nrow(newds)< typemax(Int32) ? Val(Int32) : Val(Int64), stable = false, mapformats = false, threads = threads)

    reps = _find_counts_for_join(view(a[1], nrow(dsl)+1:length(a[1])), a[3])
    gslots, minval, sz = _create_dictionary_for_join_int(identity, view(a[1], nrow(dsl)+1:length(a[1])), reps, 1, a[3], Val(T))

    ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
    where = Vector{T}(undef, length(reps)+1)
    cumsum!(view(where, 2:length(where)), reps)
    where[1] = 0
    _find_range_for_join!(ranges, view(a[1], 1:nrow(dsl)), gslots, reps, where, 1, sz, threads = threads)
    ranges, a, gslots, minval, reps, sz, right_cols
end

function _join_left_dict(dsl, dsr, ranges, onleft, onright, right_cols, ::Val{T}; makeunique = makeunique, mapformats = mapformats, check = check, threads = true, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id ) where T
    _fl = _date_value∘identity
    _fr = _date_value∘identity
    if mapformats[1]
        _fl = _date_value∘getformat(dsl, onleft[1])
    end
    if mapformats[2]
        _fr = _date_value∘getformat(dsr, onright[1])
    end
    dict, maxprob, sz, fallback, type = _create_dictionary_for_join(_fr, _columns(dsr)[onright[1]], _fl, _columns(dsl)[onleft[1]], Val(T))
    # key is not unique, fall back to sort
    if fallback
        return false, Dataset()
    end

    _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, _columns(dsl)[onleft[1]], _columns(dsr)[onright[1]], sz, type, threads = threads)

    new_ends = map(x -> max(1, length(x)), ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if multiple_match
        multiple_match_col = _create_multiple_match_col_left(ranges, total_length)
    end

    res = []
    for j in 1:length(index(dsl))
        addmissing = false
        _res = allocatecol(_columns(dsl)[j], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            # fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_oncols_left_table_left!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length, missing, threads = threads)
        else
            _fill_oncols_left_table_left!(_res, _columns(dsl)[j], ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(res, _res)

    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, DataAPI.refarray(_columns(dsr)[right_cols[j]]), ranges, new_ends, total_length, fill_val, threads = threads)
        else
            _fill_right_cols_table_left!(_res, _columns(dsr)[right_cols[j]], ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(_columns(newds), _res)

        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing; threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, 1:nrow(dsr), ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    true, newds

end

function _join_left!_dict(dsl, dsr, ranges, onleft, onright, right_cols, ::Val{T}; makeunique = makeunique, mapformats = mapformats, check = check, threads = true, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    _fl = _date_value∘identity
    _fr = _date_value∘identity
    if mapformats[1]
        _fl = _date_value∘getformat(dsl, onleft[1])
    end
    if mapformats[2]
        _fr = _date_value∘getformat(dsr, onright[1])
    end
    dict, maxprob, sz, fallback, type = _create_dictionary_for_join(_fr, _columns(dsr)[onright[1]], _fl, _columns(dsl)[onleft[1]], Val(T))
    # key is not unique, fall back to sort
    if fallback
        return false, Dataset()
    end

    _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, _columns(dsl)[onleft[1]], _columns(dsr)[onright[1]], sz, type, threads = threads)
    if !all(x->length(x) <= 1, ranges)
        throw(ArgumentError("`leftjoin!` can only be used when each observation in left data set matches at most one observation from right data set"))
    end
    new_ends = map(x -> max(1, length(x)), ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if multiple_match
        multiple_match_col = _create_multiple_match_col_left(ranges, total_length)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, DataAPI.refarray(_columns(dsr)[right_cols[j]]), ranges, new_ends, total_length, fill_val, threads = threads)
        else
            _fill_right_cols_table_left!(_res, _columns(dsr)[right_cols[j]], ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(_columns(dsl), _res)
        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(dsl), new_var_name)
        setformat!(dsl, index(dsl)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(dsl, ncol(dsl)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(dsl, ncol(dsl)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, 1:nrow(dsr), ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(dsl, ncol(dsl)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    _modified(_attributes(dsl))
    true, dsl
end

function _join_inner_dict(dsl, dsr, ranges, onleft, onright, right_cols, ::Val{T}; makeunique = makeunique, mapformats = mapformats, check = check, threads = true, multiple_match=false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :obs_id) where T
    _fl = _date_value∘identity
    _fr = _date_value∘identity
    if mapformats[1]
        _fl = _date_value∘getformat(dsl, onleft[1])
    end
    if mapformats[2]
        _fr = _date_value∘getformat(dsr, onright[1])
    end

    dict, maxprob, sz, fallback, type = _create_dictionary_for_join(_fr, _columns(dsr)[onright[1]], _fl, _columns(dsl)[onleft[1]], Val(T))
    if fallback
        return false, Dataset()
    end
    _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, _columns(dsl)[onleft[1]], _columns(dsr)[onright[1]], sz, type, threads = threads)

    new_ends = map(length, ranges)
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end]

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if multiple_match
        multiple_match_col = _create_multiple_match_col_inner(ranges, nothing, total_length)
    end

    res = []
    for j in 1:length(index(dsl))
        _res = allocatecol(_columns(dsl)[j], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            _fill_oncols_left_table_inner!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length, threads = threads)
        else
            _fill_oncols_left_table_inner!(_res, _columns(dsl)[j], ranges, new_ends, total_length, threads = threads)
        end
        push!(res, _res)
    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length, addmissing = false)
        if DataAPI.refpool(_res) !== nothing
            _fill_right_cols_table_inner!(_res.refs, DataAPI.refarray(_columns(dsr)[right_cols[j]]), ranges, new_ends, total_length, threads = threads)
        else
            _fill_right_cols_table_inner!(_res, _columns(dsr)[right_cols[j]], ranges, new_ends, total_length, threads = threads)
        end
        push!(_columns(newds), _res)

        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_inner!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_inner!(obs_id_right, 1:nrow(dsr), ranges, new_ends, total_length, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    true, newds

end

function _join_outer_dict(dsl, dsr, ranges, onleft, onright, oncols_left, oncols_right, right_cols, ::Val{T}; makeunique = makeunique, mapformats = mapformats, check = check, threads = true, source::Bool = false, source_col_name = :source, multiple_match = false, multiple_match_name = :multiple, obs_id = [false, false], obs_id_name = :__OBSID__) where T
    _fl = _date_value∘identity
    _fr = _date_value∘identity
    if mapformats[1]
        _fl = _date_value∘getformat(dsl, onleft[1])
    end
    if mapformats[2]
        _fr = _date_value∘getformat(dsr, onright[1])
    end

    dict, maxprob, sz, fallback, type = _create_dictionary_for_join(_fr, _columns(dsr)[onright[1]], _fl, _columns(dsl)[onleft[1]], Val(T))
    if fallback
        return false, Dataset()
    end
    _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, _columns(dsl)[onleft[1]], _columns(dsr)[onright[1]], sz, type, threads = threads)
    new_ends = map(x -> max(1, length(x)), ranges)
    notinleft = _find_right_not_in_left(ranges, nrow(dsr), 1:nrow(dsr))
    cumsum!(new_ends, new_ends)
    total_length = new_ends[end] + length(notinleft)

    if check
        @assert total_length < 10*nrow(dsl) "the output data set will be very large ($(total_length)×$(ncol(dsl)+length(right_cols))) compared to the left data set size ($(nrow(dsl))×$(ncol(dsl))), make sure that the `on` keyword is selected properly, alternatively, pass `check = false` to ignore this error."
    end

    if source
        source_col = _create_source_for_outer(ranges, notinleft, total_length, new_ends)
    end
    if multiple_match
        multiple_match_col = _create_multiple_match_col_outer(ranges, notinleft, total_length, new_ends)
    end

    res = []
    for j in 1:length(index(dsl))
        _res = allocatecol(_columns(dsl)[j], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_oncols_left_table_left!(_res.refs, DataAPI.refarray(_columns(dsl)[j]), ranges, new_ends, total_length, fill_val, threads = threads)
        else
            _fill_oncols_left_table_left!(_res, _columns(dsl)[j], ranges, new_ends, total_length, missing, threads = threads)
        end
        push!(res, _res)
    end
    for j in 1:length(oncols_left)
        _fill_oncols_left_table_left_outer!(res[oncols_left[j]], _columns(dsr)[oncols_right[j]], notinleft, new_ends, total_length)
    end
    if dsl isa SubDataset
        newds = Dataset(res, copy(index(dsl)), copycols = false)
    else
        newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)), copycols = false)
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_left!(_res.refs, DataAPI.refarray(_columns(dsr)[right_cols[j]]), ranges, new_ends, total_length, fill_val, threads = threads)
            _fill_oncols_left_table_left_outer!(_res.refs, DataAPI.refarray(_columns(dsr)[right_cols[j]]), notinleft, new_ends, total_length)
        else
            _fill_right_cols_table_left!(_res, _columns(dsr)[right_cols[j]], ranges, new_ends, total_length, missing, threads = threads)
            _fill_oncols_left_table_left_outer!(_res, _columns(dsr)[right_cols[j]], notinleft, new_ends, total_length)
        end
        push!(_columns(newds), _res)
        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end
    if source
        insertcols!(newds, ncol(newds)+1, source_col_name => source_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if multiple_match
        insertcols!(newds, ncol(newds)+1, multiple_match_name => multiple_match_col, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[1]
        # Note that the name convention obs_id_name1 and name2 are used in other places
        obs_id_name1 = Symbol(obs_id_name, "_left")
        obs_id_left = allocatecol(nrow(dsl) < typemax(Int32) ? Int32 : Int64, total_length)
        _fill_oncols_left_table_left!(obs_id_left, 1:nrow(dsl), ranges, new_ends, total_length, missing, threads = threads)
        insertcols!(newds, ncol(newds)+1, obs_id_name1 => obs_id_left, unsupported_copy_cols = false, makeunique = makeunique)
    end
    if obs_id[2]
        obs_id_name2 = Symbol(obs_id_name, "_right")
        obs_id_right = allocatecol(T, total_length)
        _fill_right_cols_table_left!(obs_id_right, 1:nrow(dsr), ranges, new_ends, total_length, missing, threads = threads)
        _fill_oncols_left_table_left_outer!(obs_id_right, 1:nrow(dsr), notinleft, new_ends, total_length)
        insertcols!(newds, ncol(newds)+1, obs_id_name2 => obs_id_right, unsupported_copy_cols = false, makeunique = makeunique)
    end
    true, newds

end

function _update!_dict(dsl, dsr, ranges, onleft, onright, right_cols, ::Val{T}; allowmissing = true, mode = :all, mapformats = [true, true], stable = false, alg = HeapSort, threads = threads, op = nothing) where T
    _fl = _date_value∘identity
    _fr = _date_value∘identity
    if mapformats[1]
        _fl = _date_value∘getformat(dsl, onleft[1])
    end
    if mapformats[2]
        _fr = _date_value∘getformat(dsr, onright[1])
    end
    dict, maxprob, sz, fallback, type = _create_dictionary_for_join(_fr, _columns(dsr)[onright[1]], _fl, _columns(dsl)[onleft[1]], Val(T))
    # key is not unique, fall back to sort
    if fallback
        return false, Dataset()
    end

    _fill_ranges_for_dict_join!(ranges, dict, maxprob, _fl, _fr, _columns(dsl)[onleft[1]], _columns(dsr)[onright[1]], sz, type, threads = threads)

    if mode == :all
        f_mode = x->true
    elseif mode == :missing || mode == :missings
        f_mode = x->ismissing(x)
    else
        f_mode = x->mode(x)::Bool
    end

    for j in 1:length(right_cols)
        if haskey(index(dsl).lookup, _names(dsr)[right_cols[j]])
            left_cols_idx = index(dsl)[_names(dsr)[right_cols[j]]]
            TL = nonmissingtype(eltype(_columns(dsl)[left_cols_idx]))
            TR = nonmissingtype(eltype(_columns(dsr)[right_cols[j]]))
            if promote_type(TR, TL) <: TL
                _update_left_with_right!(_columns(dsl)[left_cols_idx], _columns(dsr)[right_cols[j]], ranges, allowmissing, f_mode, threads = threads, op = op)
            end
        end
    end
    _modified(_attributes(dsl))
    true, dsl
end

# a new idea for joining without sorting
function _in_hash(dsl::AbstractDataset, dsr::AbstractDataset, ::Val{T}; onleft, onright, mapformats = [true, true], threads = true) where T
    isempty(dsl) && return Bool[]
    oncols_left = onleft
    oncols_right = onright

    # use Set when there is only one column in `on`
    cols = Any[]
    for j in 1:length(oncols_left)
        if mapformats[1]
            fl = getformat(dsl, oncols_left[j])
        else
            fl = identity
        end
        if mapformats[2]
            fr = getformat(dsr, oncols_right[j])
        else
            fr = identity
        end
        push!(cols, Cat2Vec(_columns(dsl)[oncols_left[j]], _columns(dsr)[oncols_right[j]], fl, fr))
    end
    newds = Dataset(cols, :auto, copycols = false)
    a = _gather_groups(newds, :, nrow(newds)< typemax(Int32) ? Val(Int32) : Val(Int64), stable = false, mapformats = false, threads = threads)
    res = _in_use_Set_int(view(a[1], 1:nrow(dsl)), view(a[1], nrow(dsl)+1:length(a[1])), 1, a[3]; threads = threads)
end

function _create_Set_for_join_int(f, v, minval, rangelen)
    flag = false
    offset = 1 - minval
    n = length(v)
    sz = rangelen + 1
    gslots = falses(sz)
    @inbounds for i in 1:length(v)
        _fv = f(v[i])
        if ismissing(_fv)
            slotix = sz
        else
            slotix = _fv + offset
        end
        if !gslots[slotix]
            gslots[slotix] = true
        end
    end
    gslots, minval, sz
end

function _query_Set_for_join_int(f, fv, gslots, minval, sz)
    offset = 1 - minval
    if ismissing(fv)
        slotix = sz
    else
        slotix = fv + offset
        !(slotix in 1:sz-1) && return 0
    end

    if slotix in 1:sz
        rowid = gslots[slotix]
        return rowid
    end
    false
end

function _in_use_Set_int_barrier!(res, ldata, gslots, minval, sz; threads = true)
    @_threadsfor threads for i in 1:length(res)
        res[i] = _query_Set_for_join_int(identity, ldata[i], gslots, minval, sz)
    end
end

function _in_use_Set_int(ldata, rdata, minval, rangelen; threads = true)
    gslots, minval, sz  =  _create_Set_for_join_int(identity, rdata, minval, rangelen)
    res = Vector{Bool}(undef, length(ldata))
    _in_use_Set_int_barrier!(res, ldata, gslots, minval, sz; threads = threads)
    res
end

# f is a function which should be applied on each element of v
# v is a vector of Int with minimum minval and range length of rangelen
# reps gives how many times a specific integer will appear in v
# no missing in v
function _create_dictionary_for_join_int(f, v, reps, minval, rangelen, ::Val{T}) where T
    offset = 1 - minval
    n = length(v)
    where = Vector{T}(undef, length(reps)+1)
    cumsum!(view(where, 2:length(where)), reps)
    where[1] = 0
    gslots = zeros(T, where[end])
    sz = rangelen
    @inbounds for i in 1:length(v)
        _fv = f(v[i])
        slotix = _fv + offset
        gslots[where[slotix]+1] = i
        where[slotix] += 1
    end
    gslots, minval, sz
end
# there is no missing in v
# here can be defined as:
# where = Vector{T}(undef, length(reps)+1)
# cumsum!(view(where, 2:length(where)), reps)
# where[1] = 0
function _query_dictionary_for_join_int(f, v, gslots,reps, where, minval, sz)
    offset = 1 - minval
    slotix = v + offset
    !(slotix in 1:sz) && return 1:0
    if slotix in 1:sz
        rowid = reps[slotix]
        rowid == 0 && return 1:0
        return where[slotix]+1:where[slotix+1]
    end
end

function _in_use_Set(ldata, rdata, _fl, _fr; threads = true)

    ss = Set(Base.Generator(_fr, rdata));
    res = Vector{Bool}(undef, length(ldata))
    @_threadsfor threads for i in 1:length(res)
        res[i] = _fl(ldata[i]) in ss
    end
    res
end


text:
@testset "clustercoeffs" begin
    A = load_matrix_network("clique-10")
    v = ones(Int64,10)
    @testset "MatrixNetwork" begin
        cc = clustercoeffs(MatrixNetwork(A))
        @test v == cc
    end
    @testset "SparseMatrixCSC" begin
        cc = clustercoeffs(A)
        @test v == cc
    end
    @testset "CSR" begin
        cc = clustercoeffs(sparse_to_csr(A)...)
        @test v == cc
    end
    @testset "triplet" begin
        cc = clustercoeffs(findnz(A)[1], findnz(A)[2])
        @test v == cc
    end
    @testset "error throwing" begin
        # not undirected
        bad_mat = spdiagm(1=>[1;1]) #spdiagm(([1; 1],), 1, 3, 3)

        @test_throws ErrorException clustercoeffs(MatrixNetwork(bad_mat))
        @test_throws ErrorException clustercoeffs(bad_mat)
        # negative weights
        @test_throws ErrorException clustercoeffs(MatrixNetwork(-sparse(1.0I,4,4)))
        @test_throws ErrorException clustercoeffs(-sparse(1.0I,4,4))
    end
    cc_mn = clustercoeffs(MatrixNetwork(A), false, false)
    cc_csc = clustercoeffs(A, false, false)
    cc_csr = clustercoeffs(sparse_to_csr(A)..., false, false)
    cc_tri = clustercoeffs(findnz(A)[1], findnz(A)[2], false, false)
    @test cc_mn == cc_csc
    @test cc_csr == cc_csc
    @test cc_tri == cc_csc
end


text:
using DynamicHMCExamples: rel_path
using Literate
using Documenter
using Random
Random.seed!(UInt32[0x57a97f0d, 0x1a38664c, 0x0dddb228, 0x7dbba96f])

# generate using Literate

DOCROOT = rel_path("../docs")
DOCSOURCE = joinpath(DOCROOT, "src")
EXAMPLES = ["independent_bernoulli", "linear_regression"]

for example in EXAMPLES
    Literate.markdown(rel_path("example_$(example).jl"), DOCSOURCE)
end

# render & deploy using Documenter

makedocs(root = DOCROOT,
         modules = Module[],
         format = :html,
         sitename = "DynamicHMCExamples.jl",
         pages = vcat(Any["index.md"],
                      Any["example_$(example).md" for example in EXAMPLES]))

deploydocs(root = DOCROOT,
           repo = "github.com/tpapp/DynamicHMCExamples.jl.git",
           target = "build",
           osname = "linux",
           julia = "1.0",
           deps = nothing,
           make = nothing)


text:
module MultiwayNumberPartitioning

using JuMP, LinearAlgebra

export partition, partition_min_largest!, partition_max_smallest!, partition_min_range!

"""
    partition([objective! = partition_min_range!], S::AbstractVector, k::Int; optimizer, strict::Bool=true) -> Vector{Int}

Find a partition of `S` into `k` subsets such the sums of the subsets are "as near as possible".

Arguments:
* (optional): an objective, such as [`partition_min_largest!`](@ref), [`partition_max_smallest!`](@ref), or [`partition_min_range!`](@ref). Defaults to `partition_min_range!`.
* `S` a collection of values to partition,
* `k` a number of subsets to partition into

Keyword arguments:
* `optimizer` (required): a JuMP-compatible mixed-integer optimizer, such as GLPK, HiGHS, Cbc, Gurobi, etc.
* `strict` (optional): if `true`, an error will be thrown if the model was not solved optimally. If `false`, only a warning will be issued instead. Defaults to `true`.

Returns a vector of indices `v` such that `v[i] == j` if `S[i]` is assigned to subset `j`.
"""
function partition(S, k; optimizer, strict::Bool=true)
    return partition(partition_min_range!, S, k; optimizer, strict)
end

function partition(objective!, S::AbstractVector, k::Int; optimizer, strict::Bool=true)
    model = Model(optimizer)
    element = populate_model!(model, S, k)
    objective!(model)
    JuMP.optimize!(model)

    # Check that we've solved it optimally
    if !(termination_status(model) == MOI.OPTIMAL &&
         primal_status(model) == MOI.FEASIBLE_POINT)
        msg = "Problem was not solved optimally or solution was not feasible. Termination status: $(termination_status(model)). Primal status: $(primal_status(model))."
        if strict
            error(mg)
        else
            @warn msg
        end
    end

    # Now we've got a sort of one-hot encoded partitioning:
    onehot_partition = value.(element)

    # We wish to convert it to a vector of indexes so that the `i`th element of the
    # vector corresponds to the subset `j` that element `i` belongs to:
    partition = dropdims(mapslices(row -> only(findall(>(0.5), row)), onehot_partition;
                                   dims=2); dims=2)

    return partition
end

# Populates an empty `model` with variables and objectives
function populate_model!(model, S, k)
    # does element i belong to set j
    @variable(model, element[i=1:length(S), j=1:k], binary = true)

    # Element `i` belongs to exactly one set j
    @constraint(model, [i = 1:length(S)], sum(element[i, :]) == 1)

    # The sum of the `j`th subset can be expressed as a dot product between the `j`th column of `element` and `S`.
    # Why? Because the `j`th column of `element` is a boolean vector expressing which elements are in set `j`.
    # If element `i` is in set `j`, then `element[i,j]==1` and we add `element[i,j]*S[i] = S[i]` to our total.
    # If element `i` is *not* in set `j`, then `element[i,j]=0` and we add nothing to our total.
    # Thus, the total is the sum of subset `j`.
    @expression(model, subset_sum[j=1:k], dot(element[:, j], S))

    # We can remove some non-uniqueness by ordering our subset sums.
    # This is also convenient for getting the smallest and largest elements.
    @constraint(model, [j = 1:(k - 1)], subset_sum[j] <= subset_sum[j + 1])

    return element
end

"""
    partition_min_largest!(model)

Adds an objective function to the model to achieve the following:

* Given a partition of `S` into `k` subsets, consider the sum of each subset
* Minimize the *largest* such sum over all subsets

"""
partition_min_largest!(model) = @objective(model, Min, last(model[:subset_sum]))

"""
    partition_max_smallest!(model)

Adds an objective function to the model to achieve the following:

* Given a partition of `S` into `k` subsets, consider the sum of each subset
* Maximize the *smallest* such sum over all subsets

"""
partition_max_smallest!(model) = @objective(model, Max, first(model[:subset_sum]))

"""
    partition_min_range!(model)

Adds an objective function to the model to achieve the following:

* Given a partition of `S` into `k` subsets, consider the sum of each subset
* Minimize the difference between the largest and smallest such sum over all subsets

"""
function partition_min_range!(model)
    @objective(model, Min, last(model[:subset_sum]) - first(model[:subset_sum]))
end

end # module


text:
# Make a dataset with blob glandsizes per fram
# BlobID, first frame mean (x, y), all frames (size)

function extractglandsdataframe(A::AbstractArray{T}) where T <: Integer
    minimum_blob_size = 35 #sweat.minimum_blob_size
    # Create a Dict with new cluster ids -> Dict(old cluster ids -> count)
    dicts = [Dict{T, NTuple{3, Int}}() for _ in axes(A, 3)]
    for frame in axes(A, 3)
        dict = dicts[frame]
        for x in axes(A,1), y in axes(A, 2)
            a = A[x, y, frame]
            if a != 0
                count, sumx, sumy = get!(dict, a) do
                    (0,0,0)
                end
                dict[a] = (count+1, sumx+x, sumy+y)
            end
        end
    end
    seen = Dict{T, Any}()
    for (frame, dict) in enumerate(dicts)
        for (id, stuff) in dict
            (count, sumx, sumy) = stuff
            if count >= minimum_blob_size
                if !haskey(seen, id)
                    seen[id] = ((sumx/count, sumy/count),[])
                end
                a = seen[id]
                push!(a[2], (frame, count))
            end
        end
    end
    seen
end

# Remove glands that don't exist on last frame
# Area needs to grow after one half of the lifetime
function extractglands(A)
    a = extractglandsdataframe(A)
    glands2 = NamedTuple{(:id, :y, :x, :frame, :area)}[]
    id = 0

    foreach(values(a)) do (xy, t)
        lastframe, lastarea = t[end]
        if length(t) >= 2 && lastframe == size(A, 3)
            changes = count(i -> t[i][2] != t[i+1][2], 1:length(t)-1)
            if 4changes >= length(t) || lastarea >= 100
                id += 1
                append!(glands2, [NamedTuple{(:id, :y, :x, :frame, :area)}((id, xy..., t1...)) for t1 in t])
            end
        end
    end
    glands2
end

function split(f, A::TA) where TA
    d = Vector{TA}()
    index = Dict()
    for a in A
        fa = f(a)
        i = get!(index, fa) do
            push!(d, similar(A, 0))
            length(d)
        end
        push!(d[i], a)
    end
    d
end
# Remove glands that are too close to larger glands
function removetooclose(A, R)
    As = split(x -> x.frame, A)
    B = similar(A, 0)
    N = 0
    R2 = R^2
    d2((x1, y1), (x2, y2)) = (x1-x2)^2 + (y1-y2)^2
    foreach(As) do A
        xy = [(x.x, x.y) for x in A]
        for i in eachindex(xy)
            dists = [d2(x, xy[i]) for x in xy]
            competitors = [(x.area, i) for (i, x) in enumerate(A) if dists[i] < R2]
            if length(competitors) >= 2
                _, j = maximum(competitors)
                if j == i
                    push!(B, A[i])
                else
                    N += 1
                end
            else
                push!(B, A[i])
            end
        end
    end
    if N != 0
        println("Removed $N points being too close")
    end
    B, N
end


text:
using QuadEig
using Test, Random, LinearAlgebra, SparseArrays

function check(A₀, A₁, A₂; atol = sqrt(eps(real(eltype(A₀)))))
    A, B, V = linearize(A₀, A₁, A₂)
    Ad, Bd, Vd = deflate(A₀, A₁, A₂; atol)
    λ, φ = eigen(Matrix(A), Matrix(B), sortby = abs)
    λd, φd = eigen(Matrix(Ad), Matrix(Bd), sortby = sortby = abs)
    v = normalize_vecs(V * φ)
    vd = normalize_vecs(Vd * φd)
    tests = Bool[]
    for (id, ld) in enumerate(λd)
        (isnan(ld) || !(atol <= abs(ld) <= inv(atol))) && continue
        test = false
        for (i, l) in enumerate(λ)
            if ≈(l, ld) && ≈(v[:, i], vd[:, id])
                test = true
                break
            end
        end
        push!(tests, test)
    end
    return all(tests)
end

function normalize_vecs(φ)
    for col in eachcol(φ)
        normalize!(col)
        col ./= sum(col) # phase fixing
    end
    return φ
end

@testset "QuadEig.jl" begin
    Random.seed!(1234)

    A₀ = rand(6,6); A₁ = rand(6, 6); A₂ = rand(6, 6);
    A₀[:, 1:3] .= A₀[:, 4:6]; A₂[:, 2:3] .= A₂[:, 4:5];
    @test check(A₀, A₁, A₂)

    A₀ = rand(ComplexF64, 6, 6); A₁ = rand(6, 6); A₂ = rand(6, 6);
    A₀[:, 1:3] .= A₀[:, 4:6]; A₂[:, 2:3] .= A₂[:, 4:5];
    @test check(A₀, A₁, A₂)

    A₀ = sprand(ComplexF64, 100, 100, 0.1); A₁ = sprand(ComplexF64, 100, 100, 0.1); A₂ = sprand(ComplexF64, 100, 100, 0.1);
    A₀[:, 1:3] .= A₀[:, 4:6]; A₂[:, 2:3] .= A₂[:, 4:5];
    @test check(A₀, A₁, A₂)
end


text:
#unique!(push!(LOAD_PATH, "~/Documents/UGent/PhD/Code_cMPS/julia1.0/CMPSKit.jl/src/"))
unique!(push!(LOAD_PATH, joinpath(pwd(), "src")))

using Revise
using CMPSKit
using KrylovKit
using OptimKit
using LinearAlgebra
using JLD2
using TensorOperations
using Plots

D = 4
k = 1.
μ = 2.3
c = 1.5

x = randn(D,D)
f = qr(x)
diagR = sign.(real(diag(f.R)))
diagR[diagR.==0] .= 1
diagRm = diagm(diagR)
V = f.Q * diagRm

V = Constant(V)
S1 = Constant(diagm(rand(D)/D))
S2 = Constant(diagm(rand(D)/D))
Q = Constant(Matrix(Symmetric(randn(D,D))))
#Q = Constant(randn(D,D))
R1 = V*S1*inv(V)
R2 = V*S2*inv(V)
Ss = (S1,S2)
Rs = (R1,R2)

#Put them in cMPS form
Ψ = InfiniteCMPS(Q, Rs)

h = k * (∂ψ[1]'*∂ψ[1] + ∂ψ[2]'*∂ψ[2]) - μ * (ψ[1]'*ψ[1] + ψ[2]'*ψ[2]) + c * ((ψ[1]')^2*ψ[1]^2 + (ψ[2]')^2*ψ[2]^2)
H = ∫(h, (-Inf,+Inf))

alg1 = LBFGS(; verbosity = 2, maxiter = 1000000, gradtol = 1e-3);
linalg = GMRES(krylovdim = 80; tol = 1e-5)

Ψ, ρL, ρR, E, e, normgrad, numfg, history = groundstate7(H, Ψ, V, Ss; optalg = alg1, linalg = linalg)
# αs,fs, dfs1, dfs2 = groundstate7(H, Ψ, V, Ss; optalg = alg1, linalg = linalg)
# display(plot(αs,[dfs1,dfs2]))
# gui()


text:
using Revise
module QFactors
using Revise

using  CSV, DataFrames, Serialization, Dates, Statistics,
  Distributed, Finometrics, Measures, StatsBase,
  Formatting, CodecLz4, CodecZlib,  Distributions, LinearAlgebra

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

import Base.eltype, Base.getindex, Base.iterate, Base.length

#necessary due to exporting conflict
describe(df::DataFrame) = DataFrames.describe(df)

############################CONSTANTS

const OVERLEAF_PATH = "C:\\Users\\Clinton\\Dropbox\\Apps\\Overleaf\\QFactors"
const TABLE_PATH = "$OVERLEAF_PATH\\sub\\tab"

const OUT_PATH = pwd() * "\\output"
const DATA_PATH = pwd() * "\\data"

const CRSP_TYPE = :monthly
const CRSP_SUFFIX = CRSP_TYPE == :monthly ? 'M' : 'D'
const CRSP_NAME = "CRSP-$CRSP_SUFFIX"
const CRSP_DATE_FORMAT = dateformat"yyyymmdd"

const COMP_NAME = "COMP"
const COMP_A_NAME = "COMP-A"
const COMP_Q_NAME = "COMP-Q"
const COMP_DATE_FORMAT = dateformat"yyyymmdd"

const CCM_NAME = "CCM"
const CCM_DATE_FORMAT = dateformat"yyyymmdd"

const LAST_DATE_STRING = "20191231"
const LAST_DATE = Date(LAST_DATE_STRING, CCM_DATE_FORMAT)

const UNIV_NAME = "$(CRSP_NAME)_$(COMP_NAME)_UNIV"

const FACT_NAME = "$(CRSP_NAME)_$(COMP_NAME)_FACT"

const TEST_OUTPUT = true

const TIME_2_STALE = Month(6)
const MIN_DATE = Date(1972,1,1)
const MAX_DATE = Date(2018,12,31)

const VWCRSP_NAME = "VWCRSP-D"
const RFCRSP_NAME = "RFCRSP-M"
const IDXCRSP_PATH = "$DATA_PATH\\crspidx"
const IDXCRSP_DATE_FORMAT = CRSP_DATE_FORMAT
const IDX_DATE_COL = :caldt

const FF_PORT_NAME = "FF-25-SIZE-BM-D"
const FF_PATH = "$DATA_PATH\\french"
const FF_PORT_RENAME =   Dict(:smalllobm=>:me1bm1, :smallhibm=>:me1bm5,
  :biglobm=>:me5bm1, :bighibm=>:me5bm5)
const FF_PORT_VALS = [Symbol("me$(i)bm$(j)") for i ∈ 1:5 for j ∈ 1:5]
const FF_DATE_FORMAT = CRSP_DATE_FORMAT

#############################Multi-threading

#available lock if needed- use lock(spinlock) unlock(spinlock) etc
const spinlock = Threads.SpinLock()

const PARALLEL = Ref{Bool}(true)

#NOTE: This macro allows for the parallelization routines to be shut down
#by flicking the above PARALLEL switch
macro mpar(expr)
  if PARALLEL[]
    quote
      :($(Threads.@threads($expr)))
    end
  else
    quote
      :($($expr))
    end
  end
end

###########################IO Compression/Decompression routines
IN_CSV_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_CSV_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

IN_JLS_STREAM(p::String) = LZ4DecompressorStream(open(p))
IN_JLS_STREAM(F::Function, p::String) = open(F, LZ4DecompressorStream, p)

OUT_JLS_STREAM(p::String) = LZ4CompressorStream(open(p, "w"))
OUT_JLS_STREAM(F::Function, p::String) = open(F, LZ4CompressorStream, p, "w")

##############################EXPORTS

############################FILES
include("QUtilities.jl")
include("FormFactors\\PreprocessCRSP.jl")
include("FormFactors\\PreprocessCOMP.jl")
include("FormFactors\\MergeCOMPCCM.jl")
include("FormFactors\\MergeCRSPCOMP.jl")
include("FormFactors\\PortfolioConstruction.jl")
include("FormFactors\\FormFactors.jl")
include("FormFactors\\MergeFACTOthers.jl")
include("AnalyzeFactors.jl")
include("AnalyzeQFactors.jl")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#=

This file contains implementations of batched linear algebra code

=#

"Given 3D CuArrays A and B, fill in X[:,:,k] = A[:,:,k] \\ B[:,:,k]"
function batch_solve!(X::CuArray{FT,3}, A::CuArray{FT,3}, B::CuArray{FT,3}) where {FT}

    # Temporary factorization matrix
    temp = similar(A)

    # LU-factorize A
    pivot, info   = CUBLAS.getrf_strided_batched!(A, true);
    synchronize()

    # Invert LU factorization of A
    CUBLAS.getri_strided_batched!(A, temp, pivot);

    # X = inv(A) * B
    NNlib.batched_mul!(X, temp, B)
    synchronize()

end

"Given 3D Julia Arrays A and B, fill in X[:,:,k] = A[:,:,k] \\ B[:,:,k]"
function batch_solve!(X::AbstractArray{FT,3}, A::AbstractArray{FT,3}, B::AbstractArray{FT,3}) where {FT}
    for i = 1:size(A, 3)
        @views ldiv!(X[:,:,i], qr!(A[:,:,i]), B[:,:,i])
    end
end

"Given 3D CuArray A, fill in X[:,:,k] = A[:,:,k] \\ I"
function batch_inv!(X::CuArray{FT,3}, A::CuArray{FT,3}) where {FT}

    # LU-factorize A
    pivot, info   = CUBLAS.getrf_strided_batched!(A, true);synchronize()
    # Invert LU factorization of A
    CUBLAS.getri_strided_batched!(A, X, pivot); synchronize()
end

"Given 3D Julia Array A, fill in X[:,:,k] = A[:,:,k] \\ I"
function batch_inv!(X::AbstractArray{FT,3}, A::AbstractArray{FT,3}) where {FT}
    for i = 1:size(A, 3)
        @views X[:,:,i] = A[:,:,i]\I;
    end
end

"Batched matrix multiply (overwrite NNlib definition)"
function batched_mul(A::CuArray{FT,3}, B::CuArray{FT,3}) where {FT}
    CUBLAS.gemm_strided_batched('N', 'N', A, B)
end

"Define batched matrix multiply for GPU and Duals"
function batched_mul(A::CuArray{ForwardDiff.Dual{T,V,N},3}, B::CuArray{ForwardDiff.Dual{T,V,N},3}) where {T,V,N}
    # Extract values:
    Av = ForwardDiff.value.(A)
    Bv = ForwardDiff.value.(B)
    # Use strided batch for A*B (defined as gemm_strided_batched):
    Cv = Av ⊠ Bv
    # Compute derivatives ∂(AB)/∂x = A * ∂B/∂x + ∂A/∂x * B;
    dABdx = [Av ⊠ ForwardDiff.partials.(B,i) + ForwardDiff.partials.(A,i) ⊠ Bv for i=1:N];
    dABdx = ForwardDiff.Partials.(tuple.(dABdx...));
    return eltype(A).(Cv,dABdx);
end

"Overload of batch_inv! for Dual numbers"
function batch_inv!(X::CuArray{ForwardDiff.Dual{T,V,N},3}, A::CuArray{ForwardDiff.Dual{T,V,N},3}) where {T,V,N}
    #@show typeof(ForwardDiff.value.(A))
    #@show T,V,N
    Atemp = ForwardDiff.value.(A)
    invA  = 0 * Atemp;

    # Set invA=A⁻¹
    batch_inv!(invA,Atemp)

    # Find sparsity (brute force)
    doIt = zeros(Bool,N)
    K    = [ForwardDiff.partials.(A,i) for i=1:N]
    doIt = [~all(iszero.(K[i])) for i=1:N]
    #@show doIt
    #dummy = 0*similar(K[1])

    # Compute derivatives ∂A⁻¹/∂x = -A⁻¹ * ∂A/∂x * A⁻¹; using NNlib batched matrix multiply
    @timeit "InvDerivs" dAdx = [doIt[i] ? -invA ⊠ K[i] ⊠ invA : K[i] for i=1:N];
    # Pack into tuples again
    dAdx = ForwardDiff.Partials.(tuple.(dAdx...));
    X .= eltype(X).(invA,dAdx);
end


text:
# # Julia is fast
#
# Very often, benchmarks are used to compare languages.  These benchmarks can lead to long discussions, first as to exactly what is being benchmarked and secondly what explains the differences.  These simple questions can sometimes get more complicated than you at first might imagine.
#
# The purpose of this notebook is for you to see a simple benchmark for yourself.  One can read the notebook and see what happened on the author's Macbook Pro with a 4-core Intel Core I7, or run the notebook yourself.
#
# (This material began life as a wonderful lecture by Steven Johnson at MIT: https://github.com/stevengj/18S096/blob/master/lectures/lecture1/Boxes-and-registers.ipynb.)

#-

# # Outline of this notebook
#
# - Define the sum function
# - Implementations & benchmarking of sum in...
#     - C (hand-written)
#     - C (hand-written with -ffast-math)
#     - python (built-in)
#     - python (numpy)
#     - python (hand-written)
#     - Julia (built-in)
#     - Julia (hand-written)
#     - Julia (hand-written with SIMD)
# - Summary of benchmarks

#-

# # `sum`: An easy enough function to understand

#-

# Consider the  **sum** function `sum(a)`, which computes
# $$
# \mathrm{sum}(a) = \sum_{i=1}^n a_i,
# $$
# where $n$ is the length of `a`.

a = rand(10^7) # 1D vector of random numbers, uniform on [0,1)

#-

sum(a)

# The expected result is 0.5 * 10^7, since the mean of each entry is 0.5

#-

# # Benchmarking a few ways in a few languages

@time sum(a)

#-

@time sum(a)

#-

@time sum(a)

# The `@time` macro can yield noisy results, so it's not our best choice for benchmarking!
#
# Luckily, Julia has a `BenchmarkTools.jl` package to make benchmarking easy and accurate:

## using Pkg
## Pkg.add("BenchmarkTools")

#-

using BenchmarkTools

# #  1. The C language
#
# C is often considered the gold standard: difficult on the human, nice for the machine. Getting within a factor of 2 of C is often satisfying. Nonetheless, even within C, there are many kinds of optimizations possible that a naive C writer may or may not get the advantage of.
#
# The current author does not speak C, so he does not read the cell below, but is happy to know that you can put C code in a Julia session, compile it, and run it. Note that the `"""` wrap a multi-line string.

using Libdl
C_code = """
#include <stddef.h>
double c_sum(size_t n, double *X) {
    double s = 0.0;
    for (size_t i = 0; i < n; ++i) {
        s += X[i];
    }
    return s;
}
"""

const Clib = tempname()   # make a temporary file

## compile to a shared library by piping C_code to gcc
## (works only if you have gcc installed):

open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(Clib * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

## define a Julia function that calls the C function:
c_sum(X::Array{Float64}) = ccall(("c_sum", Clib), Float64, (Csize_t, Ptr{Float64}), length(X), X)

#-

c_sum(a)

#-

c_sum(a) ≈ sum(a) # type \approx and then <TAB> to get the ≈ symbolb

#-

c_sum(a) - sum(a)

#-

≈  # alias for the `isapprox` function

#-

?isapprox

# We can now benchmark the C code directly from Julia:

c_bench = @benchmark c_sum($a)

#-

println("C: Fastest time was $(minimum(c_bench.times) / 1e6) msec")

#-

d = Dict()  # a "dictionary", i.e. an associative array
d["C"] = minimum(c_bench.times) / 1e6  # in milliseconds
d

#-

using Plots
gr()

#-

using Statistics # bring in statistical support for standard deviations
t = c_bench.times / 1e6 # times in milliseconds
m, σ = minimum(t), std(t)

histogram(t, bins=500,
    xlim=(m - 0.01, m + σ),
    xlabel="milliseconds", ylabel="count", label="")
## # 2. C with -ffast-math
#
## If we allow C to re-arrange the floating point operations, then it'll vectorize with SIMD (single instruction, multiple data) instructions.
const Clib_fastmath = tempname()   # make a temporary file

## The same as above but with a -ffast-math flag added
open(`gcc -fPIC -O3 -msse3 -xc -shared -ffast-math -o $(Clib_fastmath * "." * Libdl.dlext) -`, "w") do f
    print(f, C_code)
end

## define a Julia function that calls the C function:
c_sum_fastmath(X::Array{Float64}) = ccall(("c_sum", Clib_fastmath), Float64, (Csize_t, Ptr{Float64}), length(X), X)

#-

c_fastmath_bench = @benchmark $c_sum_fastmath($a)

#-

d["C -ffast-math"] = minimum(c_fastmath_bench.times) / 1e6  # in milliseconds

# # 3. Python's built in `sum`

#-

# The `PyCall` package provides a Julia interface to Python:

## using Pkg; Pkg.add("PyCall")
using PyCall

#-

## get the Python built-in "sum" function:
pysum = pybuiltin("sum")

#-

pysum(a)

#-

pysum(a) ≈ sum(a)

#-

py_list_bench = @benchmark $pysum($a)

#-

d["Python built-in"] = minimum(py_list_bench.times) / 1e6
d

# # 4. Python: `numpy`
#
# ## Takes advantage of hardware "SIMD", but only works when it works.
#
# `numpy` is an optimized C library, callable from Python.
# It may be installed within Julia as follows:

## using Pkg; Pkg.add("Conda")
using Conda

#-

## Conda.add("numpy")

#-

numpy_sum = pyimport("numpy")["sum"]

py_numpy_bench = @benchmark $numpy_sum($a)

#-

numpy_sum(a)

#-

numpy_sum(a) ≈ sum(a)

#-

d["Python numpy"] = minimum(py_numpy_bench.times) / 1e6
d

# # 5. Python, hand-written

py"""
def py_sum(A):
    s = 0.0
    for a in A:
        s += a
    return s
"""

sum_py = py"py_sum"

#-

py_hand = @benchmark $sum_py($a)

#-

sum_py(a)

#-

sum_py(a) ≈ sum(a)

#-

d["Python hand-written"] = minimum(py_hand.times) / 1e6
d

# # 6. Julia (built-in)
#
# ## Written directly in Julia, not in C!

@which sum(a)

#-

j_bench = @benchmark sum($a)

#-

d["Julia built-in"] = minimum(j_bench.times) / 1e6
d

# # 7. Julia (hand-written)

function mysum(A)
    s = 0.0 # s = zero(eltype(a))
    for a in A
        s += a
    end
    s
end

#-

j_bench_hand = @benchmark mysum($a)

#-

d["Julia hand-written"] = minimum(j_bench_hand.times) / 1e6
d

# # 8. Julia (hand-written w. simd)

function mysum_simd(A)
    s = 0.0 # s = zero(eltype(A))
    @simd for a in A
        s += a
    end
    s
end

#-

j_bench_hand_simd = @benchmark mysum_simd($a)

#-

mysum_simd(a)

#-

d["Julia hand-written simd"] = minimum(j_bench_hand_simd.times) / 1e6
d

# # Summary

for (key, value) in sort(collect(d), by=last)
    println(rpad(key, 25, "."), lpad(round(value; digits=1), 6, "."))
end

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# ---
# title: 945. Minimum Increment to Make Array Unique
# id: problem945
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/minimum-increment-to-make-array-unique/description/>
# hidden: true
# ---
#
# Given an array of integers A, a _move_ consists of choosing any `A[i]`, and
# incrementing it by `1`.
#
# Return the least number of moves to make every value in `A` unique.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,2]
#     Output: 1
#     Explanation: After 1 move, the array could be [1, 2, 3].
#
#
# **Example 2:**
#
#
#
#     Input: [3,2,1,2,1,7]
#     Output: 6
#     Explanation: After 6 moves, the array could be [3, 4, 1, 2, 5, 7].
#     It can be shown with 5 or less moves that it is impossible for the array to have all unique values.
#
#
#
#
# **Note:**
#
#   1. `0 <= A.length <= 40000`
#   2. `0 <= A[i] < 40000`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
module PanGraph

ENV["PYTHON"] = ""

using GZip
using Rematch
using Random: seed!

include("graph.jl")
using .Graphs

import .Graphs.Utility: read_fasta, write_fasta

include("simulate.jl")
using .Simulation

export Graphs, Simulation

# ------------------------------------------------------------------------
# errors

struct Error <: Exception
    msg::AbstractString
end
Base.showerror(io::IO, e::Error) = print(io, "PanGraph Error: ", e.msg)

function panic(msg...)
    print(stderr, string(msg...))
    exit(2)
end

function open(func, path, args...)
    endswith(path, ".gz") && return GZip.open(func, path, args...)
    return Base.open(func, path, args...)
end

function load(path, cmd)
   if path === nothing
       unmarshal(stdin)
   elseif length(path) == 1
       path = path[1]
       !isfile(path) && error("file '$(path)' not found")
       open(unmarshal, path)
   else
       usage(cmd)
       return 2
   end
end

# ------------------------------------------------------------------------
# subcommands and arguments

include("args.jl")
using .Commands

# ---------------------------
# All subcommands go here

include("build.jl")
include("generate.jl")
include("help.jl")
include("polish.jl")
include("marginalize.jl")
include("export.jl")

Dispatch = Command(
    "pangraph",
    "pangraph <command> [arguments]",
    "a tool for aligning large sets of closely related genomes in the presence of horizontal gene transfer",
    "passed directly to the chosen command",
    [
     Build,
     Generate,
	 Help,
     Polish,
     Marginalize,
     Export,
    ],
)

function main(args)
    if length(args) == 0
        usage(Dispatch)
        return 2
    end

    return run(Dispatch, parse(Dispatch, args))
end

function julia_main()::Cint
    try
        main(ARGS)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end

    return 0
end

if abspath(PROGRAM_FILE) == @__FILE__
    main(ARGS)
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using SafeTestsets

@safetestset "SDDEProblem, solve" begin include("test_prob_sol.jl") end
@safetestset "Analyticless Convergence Tests" begin include("analyticless_convergence_tests.jl") end
@safetestset "Event handling" begin include("events.jl") end
@safetestset "Non-Diagonal Sparse Noise" begin include("nondiagonal_sparse_noise.jl") end


text:
using Test, SimpleDrawing, Plots

@test non_colinear_check(im, 2im, 2im+1)
z = 1+im
@test !non_colinear_check(z,-z,2z)

@test isinf(find_center(im,2im,3im))
@test find_center(im,-1+0im,1+im) == 0.5 - 0.5im

y = [1;5;2;-1;3]
S = Spline(y,:open)
@test S(2)==5
S = Spline(y,:closed)
@test S(7)==5