text:
function _power_flow_solution!(
    initial_guess::Vector{Float64},
    sys::PSY.System,
    inputs::SimulationInputs,
)
    res = PSY.solve_powerflow!(sys)
    if !res
        @error("PowerFlow failed to solve")
        return BUILD_FAILED
    end
    bus_size = length(PSY.get_bus_numbers(sys))
    @debug "Updating bus voltage magnitude and angle to match power flow result"
    for bus in PSY.get_components(PSY.Bus, sys)
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(inputs)[bus_n]
        initial_guess[bus_ix] = PSY.get_magnitude(bus) * cos(PSY.get_angle(bus))
        initial_guess[bus_ix + bus_size] = PSY.get_magnitude(bus) * sin(PSY.get_angle(bus))
        @debug "$(PSY.get_name(bus)) V_r = $(initial_guess[bus_ix]), V_i = $(initial_guess[bus_ix + bus_size])"
    end
    return BUILD_INCOMPLETE
end

function _initialize_static_injection!(inputs::SimulationInputs)
    @debug "Updating Source internal voltage magnitude and angle"
    static_injection_devices = get_static_injectors(inputs)
    if !isempty(static_injection_devices)
        try
            for s in static_injection_devices
                initialize_static_device!(s)
            end
        catch e
            bt = catch_backtrace()
            @error "Static Injection Failed to Initialize" exception = e, bt
            return BUILD_FAILED
        end
    end
    return BUILD_INCOMPLETE
end

function _initialize_dynamic_injection!(
    initial_guess::Vector{Float64},
    inputs::SimulationInputs,
    system::PSY.System,
)
    @debug "Updating Dynamic Injection Component Initial Guess"
    initial_inner_vars = zeros(get_inner_vars_count(inputs))
    try
        for dynamic_device in get_dynamic_injectors(inputs)
            static = PSY.get_component(
                dynamic_device.static_type,
                system,
                PSY.get_name(dynamic_device),
            )
            @debug "Initializing $(PSY.get_name(dynamic_device)) - $(typeof(dynamic_device.device))"
            n_states = PSY.get_n_states(dynamic_device)
            _inner_vars = @view initial_inner_vars[get_inner_vars_index(dynamic_device)]
            x0_device = initialize_dynamic_device!(dynamic_device, static, _inner_vars)
            @assert length(x0_device) == n_states
            ix_range = get_ix_range(dynamic_device)
            initial_guess[ix_range] = x0_device
        end
    catch e
        bt = catch_backtrace()
        @error "Dynamic Injection Failed to Initialize" exception = e, bt
        return BUILD_FAILED
    end
    return BUILD_INCOMPLETE
end

function _initialize_dynamic_branches!(
    initial_guess::Vector{Float64},
    inputs::SimulationInputs,
)
    try
        @debug "Initializing Dynamic Branches"
        for br in get_dynamic_branches(inputs)
            @debug "$(PSY.get_name(br)) -  $(typeof(br))"
            n_states = PSY.get_n_states(br)
            x0_branch = initialize_dynamic_device!(br)
            @assert length(x0_branch) == n_states
            ix_range = get_ix_range(br)
            initial_guess[ix_range] = x0_branch
        end
    catch e
        bt = catch_backtrace()
        @error "Dynamic Branches Failed to Initialize" exception = e, bt
        return BUILD_FAILED
    end
    return BUILD_INCOMPLETE
end

function check_valid_values(initial_guess::Vector{Float64}, inputs::SimulationInputs)
    invalid_initial_guess = String[]
    for i in get_bus_range(inputs)
        if initial_guess[i] > 1.3 || initial_guess[i] < -1.3
            push!(invalid_initial_guess, "Voltage entry $i")
        end
    end

    for device in get_dynamic_injectors(inputs)
        device_initial_guess = initial_guess[get_ix_range(device)]
        device_index = get_global_index(device)
        if haskey(device_index, :ω)
            dev_freq = initial_guess[device_index[:ω]]
            if dev_freq > 1.2 || dev_freq < 0.8
                push!(invalid_initial_guess, "$(PSY.get_name(device)) - :ω")
            end
        end
        all(isfinite, initial_guess) && continue
        invalid_set = findall(!isfinite, device_initial_guess)
        for state in get_global_index(device)
            if state.second ∈ i
                push!(invalid_initial_guess, "$device - $(p.first)")
            end
        end
    end

    if !isempty(invalid_initial_guess)
        @error("Invalid initial condition values $invalid_initial_guess")
        return BUILD_FAILED
    end
    return BUILD_IN_PROGRESS
end

# Default implementation for both models. This implementation is to future proof if there is
# a divergence between the required build methods
function _calculate_initial_guess!(sim::Simulation)
    inputs = get_simulation_inputs(sim)
    @assert sim.status == BUILD_INCOMPLETE
    while sim.status == BUILD_INCOMPLETE
        @debug "Start state intialization routine"
        TimerOutputs.@timeit BUILD_TIMER "Power Flow solution" begin
            sim.status = _power_flow_solution!(sim.x0_init, get_system(sim), inputs)
        end
        TimerOutputs.@timeit BUILD_TIMER "Initialize Static Injectors" begin
            sim.status = _initialize_static_injection!(inputs)
        end
        TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Injectors" begin
            sim.status =
                _initialize_dynamic_injection!(sim.x0_init, inputs, get_system(sim))
        end
        if has_dyn_lines(inputs)
            TimerOutputs.@timeit BUILD_TIMER "Initialize Dynamic Branches" begin
                sim.status = _initialize_dynamic_branches!(sim.x0_init, inputs)
            end
        else
            @debug "No Dynamic Branches in the system"
        end
        sim.status = check_valid_values(sim.x0_init, inputs)
    end
    return
end

function precalculate_initial_conditions!(sim::Simulation)
    _calculate_initial_guess!(sim)
    return sim.status != BUILD_FAILED
end

"""
Returns a Dictionary with the resulting initial conditions of the simulation
"""
function read_initial_conditions(sim::Simulation)
    system = get_system(sim)
    simulation_inputs = get_simulation_inputs(sim)
    bus_size = get_bus_count(simulation_inputs)
    V_R = Dict{Int, Float64}()
    V_I = Dict{Int, Float64}()
    Vm = Dict{Int, Float64}()
    θ = Dict{Int, Float64}()
    for bus in PSY.get_components(PSY.Bus, system)
        bus_n = PSY.get_number(bus)
        bus_ix = get_lookup(simulation_inputs)[bus_n]
        V_R[bus_n] = sim.x0_init[bus_ix]
        V_I[bus_n] = sim.x0_init[bus_ix + bus_size]
        Vm[bus_n] = sqrt(sim.x0_init[bus_ix]^2 + sim.x0_init[bus_ix + bus_size]^2)
        θ[bus_n] = atan(sim.x0_init[bus_ix + bus_size] / sim.x0_init[bus_ix])
    end
    results = Dict{String, Any}("V_R" => V_R, "V_I" => V_I, "Vm" => Vm, "θ" => θ)
    for device in get_dynamic_injectors(simulation_inputs)
        states = PSY.get_states(device)
        name = PSY.get_name(device)
        global_index = get_global_index(device)
        x0_device = Dict{Symbol, Float64}()
        for s in states
            x0_device[s] = sim.x0_init[global_index[s]]
        end
        results[name] = x0_device
    end
    dyn_branches = get_dynamic_branches(simulation_inputs)
    if !isempty(dyn_branches)
        for br in dyn_branches
            states = PSY.get_states(br)
            name = PSY.get_name(br)
            global_index = get_global_index(br)
            x0_br = Dict{Symbol, Float64}()
            for s in states
                x0_br[s] = sim.x0_init[global_index[s]]
            end
            printed_name = "Line " * name
            results[printed_name] = x0_br
        end
    end
    return results
end


text:
# Unit tests for utils.jl

using NumericExtensions
using Base.Test

@test eachrepeat(1:3, 2) == [1, 1, 2, 2, 3, 3]
@test eachrepeat(1:4, 1:4) == [1, 2, 2, 3, 3, 3, 4, 4, 4, 4]

rmat = [1 1 2 2 3 3;
        1 1 2 2 3 3;
        1 1 2 2 3 3;
        4 4 5 5 6 6;
        4 4 5 5 6 6;
        4 4 5 5 6 6]

@test eachrepeat([1 2 3; 4 5 6], (3, 2)) == rmat

x = [1, 1, 1, 2, 2, 3, 3, 3, 3, 1, 1, 2, 1]
sx, cx = sortindexes(x)
@test sx == sortperm(x) == [1, 2, 3, 10, 11, 13, 4, 5, 12, 6, 7, 8, 9]
@test cx == [6, 3, 4]

g = groupindexes(x)
@test size(g) == (3,)
@test g[1] == [1, 2, 3, 10, 11, 13]
@test g[2] == [4, 5, 12]
@test g[3] == [6, 7, 8, 9]

ext: jl


text:
using Cairo
using LSystems

c = CairoRGBSurface(964, 964);
cr = CairoContext(c);

save(cr);

set_source_rgb(cr, 1.0, 1.0, 1.0);
rectangle(cr, 0.0, 0.0, 964.0, 964.0);
fill(cr);

restore(cr);
save(cr);

x, y, a, d = 10., 10., 180, 15
set_line_width(cr, 1);
move_to(cr, x, y);

hilbert_start = "A"
hilbert_trans = Dict([
('A', "-BF+AFA+FB-"),
('B', "+AF-BFB-FA+")
])

function determine_new_position(x, y, a, d)
    x -= sind(a) * d
    y -= cosd(a) * d
    return x, y
end

level = 6
drawing_instructions = @task dol_system(hilbert_start, hilbert_trans, level)

for v in drawing_instructions
    if v == 'F'
        x, y = determine_new_position(x, y, a, d)
        line_to(cr, x, y)
    elseif v == '+'
        a -= 90
    elseif v == '-'
        a += 90
    elseif v == '['
        push!(s, [x, y, a])
    elseif v == ']'
        x, y, a = pop!(s)
        move_to(cr, x, y)
    end
end

stroke(cr);
write_to_png(c, "hilbert_$level.png");


text:
function paralax(weapon::Gun,tank::Tank)
    h = gun.hg*(target.position-gun.zerotage)/(target.position*gun.zerotage)
    v = gun.vg*(target.position-gun.zerotage)/(target.position*gun.zerotage)
    return h,v
end

function jump()
end

function drift()
end

function ballistic()
end

function fixedBias()

end


text:
@testset "Memory Utils Tests" begin
    @test ALifeBenchmark._determine_overlap(0, 1, 0, 1) == ALifeBenchmark.ExactMatch()
    @test ALifeBenchmark._determine_overlap(0, 10, 0, 10) == ALifeBenchmark.ExactMatch()
    @test ALifeBenchmark._determine_overlap(199, 20008, 199, 20008) == ALifeBenchmark.ExactMatch()

    @test ALifeBenchmark._determine_overlap(0, 1, 1, 1) == ALifeBenchmark.LeftNeighbor()
    @test ALifeBenchmark._determine_overlap(0, 10, 10, 20) == ALifeBenchmark.LeftNeighbor()
    @test ALifeBenchmark._determine_overlap(2245, 100, 2345, 20) == ALifeBenchmark.LeftNeighbor()

    @test ALifeBenchmark._determine_overlap(1, 1, 0, 1) == ALifeBenchmark.RightNeighbor()
    @test ALifeBenchmark._determine_overlap(10, 1, 0, 10) == ALifeBenchmark.RightNeighbor()
    @test ALifeBenchmark._determine_overlap(2000, 100, 1000, 1000) == ALifeBenchmark.RightNeighbor()

    @test ALifeBenchmark._determine_overlap(0, 1, 2, 1) == ALifeBenchmark.Disjunct()
    @test ALifeBenchmark._determine_overlap(0, 10, 100, 20) == ALifeBenchmark.Disjunct()
    @test ALifeBenchmark._determine_overlap(2983, 10276, 23, 50) == ALifeBenchmark.Disjunct()

    @test ALifeBenchmark._determine_overlap(1, 8, 0, 10) == ALifeBenchmark.Within()
    @test ALifeBenchmark._determine_overlap(5, 1, 0, 10) == ALifeBenchmark.Within()
    @test ALifeBenchmark._determine_overlap(500, 207, 213, 14320) == ALifeBenchmark.Within()

    @test ALifeBenchmark._determine_overlap(0, 1, 0, 10) == ALifeBenchmark.WithinFromStart()
    @test ALifeBenchmark._determine_overlap(0, 8, 0, 10) == ALifeBenchmark.WithinFromStart()
    @test ALifeBenchmark._determine_overlap(293, 208, 293, 1340) == ALifeBenchmark.WithinFromStart()

    @test ALifeBenchmark._determine_overlap(9, 1, 0, 10) == ALifeBenchmark.WithinToEnd()
    @test ALifeBenchmark._determine_overlap(2, 8, 0, 10) == ALifeBenchmark.WithinToEnd()
    @test ALifeBenchmark._determine_overlap(300, 200, 200, 300) == ALifeBenchmark.WithinToEnd()

    @test ALifeBenchmark._determine_overlap(0, 10, 1, 8) == ALifeBenchmark.CompleteOverlap()
    @test ALifeBenchmark._determine_overlap(0, 10, 5, 1) == ALifeBenchmark.CompleteOverlap()
    @test ALifeBenchmark._determine_overlap(213, 14320, 500, 207) == ALifeBenchmark.CompleteOverlap()

    @test ALifeBenchmark._determine_overlap(0, 10, 5, 5) == ALifeBenchmark.LeftOverlap()
    @test ALifeBenchmark._determine_overlap(0, 15, 10, 10) == ALifeBenchmark.LeftOverlap()
    @test ALifeBenchmark._determine_overlap(0, 15, 14, 10) == ALifeBenchmark.LeftOverlap()
    @test ALifeBenchmark._determine_overlap(325, 123, 425, 200) == ALifeBenchmark.LeftOverlap()

    @test ALifeBenchmark._determine_overlap(0, 10, 0, 5) == ALifeBenchmark.RightOverlap()
    @test ALifeBenchmark._determine_overlap(5, 10, 0, 10) == ALifeBenchmark.RightOverlap()
    @test ALifeBenchmark._determine_overlap(9, 10, 0, 10) == ALifeBenchmark.RightOverlap()
    @test ALifeBenchmark._determine_overlap(239, 1000, 0, 1000) == ALifeBenchmark.RightOverlap()
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
@symbol_func function d(cur_reactor::AbstractReactor)
  cur_d = K_DR(cur_reactor)

  cur_d *= R_CS(cur_reactor)

  cur_d += K_DO(cur_reactor)

  cur_d
end


text:
using Base.Meta

macro constants(array, stripprefix, expr)
    ret = Expr(:block)
    # Initialize the name lookup array
    push!(ret.args,:(const $array = Dict{UInt32,String}()))
    for e in expr.args
        if !isexpr(e,:const)
            continue
        end
        eq = e.args[1]
        @assert isexpr(eq,:(=))
        name = string(eq.args[1])
        name = replace(name,stripprefix,"",1)
        push!(ret.args,e)
        push!(ret.args,:($array[UInt32($(eq.args[1]))] = $name))
    end
    return esc(ret)
end
@constants PERF_EVENT_TYPE "" begin
    # Kernel defined events
    const PERF_RECORD_MMAP              = 1
    const PERF_RECORD_LOST              = 2
    const PERF_RECORD_COMM              = 3
    const PERF_RECORD_EXIT              = 4
    const PERF_RECORD_THROTTLE          = 5
    const PERF_RECORD_UNTHROTTLE        = 6
    const PERF_RECORD_FORK              = 7
    const PERF_RECORD_READ              = 8
    const PERF_RECORD_SAMPLE            = 9
    const PERF_RECORD_MMAP2             = 10
    const PERF_RECORD_AUX               = 11
    const PERF_RECORD_ITRACE_START      = 12
    const PERF_RECORD_LOST_SAMPLES      = 13
    const PERF_RECORD_SWITCH            = 14
    const PERF_RECORD_SWITCH_CPU_WIDE   = 15
    # perf (the user space tool) defined events
    const PERF_RECORD_HEADER_ATTR          = 64
	const PERF_RECORD_HEADER_EVENT_TYPE    = 65
	const PERF_RECORD_HEADER_TRACING_DATA  = 66
	const PERF_RECORD_HEADER_BUILD_ID      = 67
	const PERF_RECORD_FINISHED_ROUND       = 68
	const PERF_RECORD_ID_INDEX             = 69
	const PERF_RECORD_AUXTRACE_INF         = 70
	const PERF_RECORD_AUXTRACE             = 71
	const PERF_RECORD_AUXTRACE_ERRO        = 72
	const PERF_RECORD_THREAD_MAP           = 73
	const PERF_RECORD_CPU_MAP              = 74
	const PERF_RECORD_STAT_CONFIG          = 75
	const PERF_RECORD_STAT                 = 76
	const PERF_RECORD_STAT_ROUND           = 77
	const PERF_RECORD_EVENT_UPDAT          = 78
	const PERF_RECORD_TIME_CONV            = 79
end
@constants PERF_TYPE "" begin
    const PERF_TYPE_HARDWARE            = 0
    const PERF_TYPE_SOFTWARE            = 1
    const PERF_TYPE_TRACEPOINT          = 2
    const PERF_TYPE_HW_CACHE            = 3
    const PERF_TYPE_RAW                 = 4
    const PERF_TYPE_BREAKPOINT          = 5
end
@constants PERF_HW "PERF_COUNT_HW_" begin
    const PERF_COUNT_HW_CPU_CYCLES              = 0
    const PERF_COUNT_HW_INSTRUCTIONS            = 1
    const PERF_COUNT_HW_CACHE_REFERENCES        = 2
    const PERF_COUNT_HW_CACHE_MISSES            = 3
    const PERF_COUNT_HW_BRANCH_INSTRUCTIONS     = 4
    const PERF_COUNT_HW_BRANCH_MISSES           = 5
    const PERF_COUNT_HW_BUS_CYCLES              = 6
    const PERF_COUNT_HW_STALLED_CYCLES_FRONTEND = 7
    const PERF_COUNT_HW_STALLED_CYCLES_BACKEND  = 8
    const PERF_COUNT_HW_REF_CPU_CYCLES          = 9
end
@constants PERF_SAMPLE_TYPE "PERF_SAMPLE_" begin
    const PERF_SAMPLE_IP                  = UInt64(1) << 0
    const PERF_SAMPLE_TID                 = UInt64(1) << 1
    const PERF_SAMPLE_TIME                = UInt64(1) << 2
    const PERF_SAMPLE_ADDR                = UInt64(1) << 3
    const PERF_SAMPLE_READ                = UInt64(1) << 4
    const PERF_SAMPLE_CALLCHAIN           = UInt64(1) << 5
    const PERF_SAMPLE_ID                  = UInt64(1) << 6
    const PERF_SAMPLE_CPU                 = UInt64(1) << 7
    const PERF_SAMPLE_PERIOD              = UInt64(1) << 8
    const PERF_SAMPLE_STREAM_ID           = UInt64(1) << 9
    const PERF_SAMPLE_RAW                 = UInt64(1) << 10
    const PERF_SAMPLE_BRANCH_STACK        = UInt64(1) << 11
    const PERF_SAMPLE_REGS_USER           = UInt64(1) << 12
    const PERF_SAMPLE_STACK_USER          = UInt64(1) << 13
    const PERF_SAMPLE_WEIGHT              = UInt64(1) << 14
    const PERF_SAMPLE_DATA_SRC            = UInt64(1) << 15
    const PERF_SAMPLE_IDENTIFIER          = UInt64(1) << 16
    const PERF_SAMPLE_TRANSACTION         = UInt64(1) << 17
    const PERF_SAMPLE_REGS_INTR           = UInt64(1) << 18
end
@constants PERF_FORMAT "" begin
    const PERF_FORMAT_TOTAL_TIME_ENABLED  = UInt64(1) << 0
    const PERF_FORMAT_TOTAL_TIME_RUNNING  = UInt64(1) << 1
    const PERF_FORMAT_ID                  = UInt64(1) << 2
    const PERF_FORMAT_GROUP               = UInt64(1) << 3
end

const PERF_CONTEXT_HV			     = -32   % UInt64
const PERF_CONTEXT_KERNEL		     = -128  % UInt64
const PERF_CONTEXT_USER              = -512  % UInt64
const PERF_CONTEXT_GUEST             = -2048 % UInt64
const PERF_CONTEXT_GUEST_KERNEL      = -2176 % UInt64
const PERF_CONTEXT_GUEST_USER        = -2560 % UInt64

const PERF_RECORD_MISC_CPUMODE_MASK      = UInt16(7) << 0
const PERF_RECORD_MISC_CPUMODE_UNKNOWN   = UInt16(0) << 0
const PERF_RECORD_MISC_KERNEL            = UInt16(1) << 0
const PERF_RECORD_MISC_USER              = UInt16(2) << 0
const PERF_RECORD_MISC_HYPERVISOR        = UInt16(3) << 0
const PERF_RECORD_MISC_GUEST_KERNEL      = UInt16(4) << 0
const PERF_RECORD_MISC_GUEST_USER        = UInt16(5) << 0

const PERF_RECORD_MISC_PROC_MAP_PARSE_TIMEOUT = UInt16(1) << 12

const PERF_RECORD_MISC_MMAP_DATA              = UInt16(1) << 13
const PERF_RECORD_MISC_COMM_EXEC              = UInt16(1) << 13
const PERF_RECORD_MISC_SWITCH_OUT             = UInt16(1) << 13

const PERF_RECORD_MISC_EXACT_IP               = UInt16(1) << 14

const PERF_RECORD_MISC_EXT_RESERVED           = UInt16(1) << 15


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

(:)(a::Real, b::Real) = (:)(promote(a,b)...)

(:)(start::T, stop::T) where {T<:Real} = UnitRange{T}(start, stop)

(:)(start::T, stop::T) where {T} = (:)(start, oftype(stop-start, 1), stop)

# promote start and stop, leaving step alone
(:)(start::A, step, stop::C) where {A<:Real,C<:Real} =
    (:)(convert(promote_type(A,C),start), step, convert(promote_type(A,C),stop))

# AbstractFloat specializations
(:)(a::T, b::T) where {T<:AbstractFloat} = (:)(a, T(1), b)

(:)(a::T, b::AbstractFloat, c::T) where {T<:Real} = (:)(promote(a,b,c)...)
(:)(a::T, b::AbstractFloat, c::T) where {T<:AbstractFloat} = (:)(promote(a,b,c)...)
(:)(a::T, b::Real, c::T) where {T<:AbstractFloat} = (:)(promote(a,b,c)...)

(:)(start::T, step::T, stop::T) where {T<:AbstractFloat} =
    _colon(OrderStyle(T), ArithmeticStyle(T), start, step, stop)
(:)(start::T, step::T, stop::T) where {T<:Real} =
    _colon(OrderStyle(T), ArithmeticStyle(T), start, step, stop)
_colon(::Ordered, ::Any, start::T, step, stop::T) where {T} = StepRange(start, step, stop)
# for T<:Union{Float16,Float32,Float64} see twiceprecision.jl
_colon(::Ordered, ::ArithmeticRounds, start::T, step, stop::T) where {T} =
    StepRangeLen(start, step, floor(Integer, (stop-start)/step)+1)
_colon(::Any, ::Any, start::T, step, stop::T) where {T} =
    StepRangeLen(start, step, floor(Integer, (stop-start)/step)+1)

"""
    (:)(start, [step], stop)

Range operator. `a:b` constructs a range from `a` to `b` with a step size of 1 (a [`UnitRange`](@ref))
, and `a:s:b` is similar but uses a step size of `s` (a [`StepRange`](@ref)).

`:` is also used in indexing to select whole dimensions
 and for [`Symbol`](@ref) literals, as in e.g. `:hello`.
"""
(:)(start::T, step, stop::T) where {T} = _colon(start, step, stop)
(:)(start::T, step, stop::T) where {T<:Real} = _colon(start, step, stop)
# without the second method above, the first method above is ambiguous with
# (:)(start::A, step, stop::C) where {A<:Real,C<:Real}
function _colon(start::T, step, stop::T) where T
    T′ = typeof(start+zero(step))
    StepRange(convert(T′,start), step, convert(T′,stop))
end

"""
    range(start, stop, length)
    range(start, stop; length, step)
    range(start; length, stop, step)
    range(;start, length, stop, step)

Construct a specialized array with evenly spaced elements and optimized storage (an [`AbstractRange`](@ref)) from the arguments.
Mathematically a range is uniquely determined by any three of `start`, `step`, `stop` and `length`.
Valid invocations of range are:
* Call `range` with any three of `start`, `step`, `stop`, `length`.
* Call `range` with two of `start`, `stop`, `length`. In this case `step` will be assumed
  to be one. If both arguments are Integers, a [`UnitRange`](@ref) will be returned.
* Call `range` with one of `stop` or `length`. `start` and `step` will be assumed to be one.

See Extended Help for additional details on the returned type.

# Examples
```jldoctest
julia> range(1, length=100)
1:100

julia> range(1, stop=100)
1:100

julia> range(1, step=5, length=100)
1:5:496

julia> range(1, step=5, stop=100)
1:5:96

julia> range(1, 10, length=101)
1.0:0.09:10.0

julia> range(1, 100, step=5)
1:5:96

julia> range(stop=10, length=5)
6:10

julia> range(stop=10, step=1, length=5)
6:1:10

julia> range(start=1, step=1, stop=10)
1:1:10

julia> range(; length = 10)
Base.OneTo(10)

julia> range(; stop = 6)
Base.OneTo(6)

julia> range(; stop = 6.5)
1.0:1.0:6.0
```
If `length` is not specified and `stop - start` is not an integer multiple of `step`, a range that ends before `stop` will be produced.
```jldoctest
julia> range(1, 3.5, step=2)
1.0:2.0:3.0
```

Special care is taken to ensure intermediate values are computed rationally.
To avoid this induced overhead, see the [`LinRange`](@ref) constructor.

!!! compat "Julia 1.1"
    `stop` as a positional argument requires at least Julia 1.1.

!!! compat "Julia 1.7"
    The versions without keyword arguments and `start` as a keyword argument
    require at least Julia 1.7.

!!! compat "Julia 1.8"
    The versions with `stop` as a sole keyword argument,
    or `length` as a sole keyword argument require at least Julia 1.8.

# Extended Help

`range` will produce a `Base.OneTo` when the arguments are Integers and
* Only `length` is provided
* Only `stop` is provided

`range` will produce a `UnitRange` when the arguments are Integers and
* Only `start`  and `stop` are provided
* Only `length` and `stop` are provided

A `UnitRange` is not produced if `step` is provided even if specified as one.
"""
function range end

range(start; stop=nothing, length::Union{Integer,Nothing}=nothing, step=nothing) =
    _range(start, step, stop, length)
range(start, stop; length::Union{Integer,Nothing}=nothing, step=nothing) = _range(start, step, stop, length)
range(start, stop, length::Integer) = _range(start, nothing, stop, length)

range(;start=nothing, stop=nothing, length::Union{Integer, Nothing}=nothing, step=nothing) =
    _range(start, step, stop, length)

_range(start::Nothing, step::Nothing, stop::Nothing, len::Nothing) = range_error(start, step, stop, len)
_range(start::Nothing, step::Nothing, stop::Nothing, len::Any    ) = range_length(len)
_range(start::Nothing, step::Nothing, stop::Any    , len::Nothing) = range_stop(stop)
_range(start::Nothing, step::Nothing, stop::Any    , len::Any    ) = range_stop_length(stop, len)
_range(start::Nothing, step::Any    , stop::Nothing, len::Nothing) = range_error(start, step, stop, len)
_range(start::Nothing, step::Any    , stop::Nothing, len::Any    ) = range_error(start, step, stop, len)
_range(start::Nothing, step::Any    , stop::Any    , len::Nothing) = range_error(start, step, stop, len)
_range(start::Nothing, step::Any    , stop::Any    , len::Any    ) = range_step_stop_length(step, stop, len)
_range(start::Any    , step::Nothing, stop::Nothing, len::Nothing) = range_error(start, step, stop, len)
_range(start::Any    , step::Nothing, stop::Nothing, len::Any    ) = range_start_length(start, len)
_range(start::Any    , step::Nothing, stop::Any    , len::Nothing) = range_start_stop(start, stop)
_range(start::Any    , step::Nothing, stop::Any    , len::Any    ) = range_start_stop_length(start, stop, len)
_range(start::Any    , step::Any    , stop::Nothing, len::Nothing) = range_error(start, step, stop, len)
_range(start::Any    , step::Any    , stop::Nothing, len::Any    ) = range_start_step_length(start, step, len)
_range(start::Any    , step::Any    , stop::Any    , len::Nothing) = range_start_step_stop(start, step, stop)
_range(start::Any    , step::Any    , stop::Any    , len::Any    ) = range_error(start, step, stop, len)

# Length as the only argument
range_length(len::Integer) = OneTo(len)

# Stop as the only argument
range_stop(stop) = range_start_stop(oneunit(stop), stop)
range_stop(stop::Integer) = range_length(stop)

# Stop and length as the only argument
range_stop_length(a::Real,          len::Integer) = UnitRange{typeof(a)}(oftype(a, a-len+1), a)
range_stop_length(a::AbstractFloat, len::Integer) = range_step_stop_length(oftype(a, 1), a, len)
range_stop_length(a,                len::Integer) = range_step_stop_length(oftype(a-a, 1), a, len)

range_step_stop_length(step, stop, length) = reverse(range_start_step_length(stop, -step, length))

range_start_length(a::Real,          len::Integer) = UnitRange{typeof(a)}(a, oftype(a, a+len-1))
range_start_length(a::AbstractFloat, len::Integer) = range_start_step_length(a, oftype(a, 1), len)
range_start_length(a,                len::Integer) = range_start_step_length(a, oftype(a-a, 1), len)

range_start_stop(start, stop) = start:stop

function range_start_step_length(a::AbstractFloat, step::AbstractFloat, len::Integer)
    range_start_step_length(promote(a, step)..., len)
end

function range_start_step_length(a::Real, step::AbstractFloat, len::Integer)
    range_start_step_length(float(a), step, len)
end

function range_start_step_length(a::AbstractFloat, step::Real, len::Integer)
    range_start_step_length(a, float(step), len)
end

function range_start_step_length(a::T, step::T, len::Integer) where {T <: AbstractFloat}
    _rangestyle(OrderStyle(T), ArithmeticStyle(T), a, step, len)
end

function range_start_step_length(a::T, step, len::Integer) where {T}
    _rangestyle(OrderStyle(T), ArithmeticStyle(T), a, step, len)
end

_rangestyle(::Ordered, ::ArithmeticWraps, a::T, step::S, len::Integer) where {T,S} =
    StepRange{typeof(a+zero(step)),S}(a, step, a+step*(len-1))
_rangestyle(::Any, ::Any, a::T, step::S, len::Integer) where {T,S} =
    StepRangeLen{typeof(a+zero(step)),T,S}(a, step, len)

range_start_step_stop(start, step, stop) = start:step:stop

function range_error(start, step, stop, length)
    hasstart  = start !== nothing
    hasstep   = step  !== nothing
    hasstop   = stop  !== nothing
    haslength = start !== nothing

    hint = if hasstart && hasstep && hasstop && haslength
        "Try specifying only three arguments"
    elseif !hasstop && !haslength
        "At least one of `length` or `stop` must be specified."
    elseif !hasstep && !haslength
        "At least one of `length` or `step` must be specified."
    elseif !hasstart && !hasstop
        "At least one of `start` or `stop` must be specified."
    else
        "Try specifying more arguments."
    end

    msg = """
    Cannot construct range from arguments:
    start = $start
    step = $step
    stop = $stop
    length = $length
    $hint
    """
    throw(ArgumentError(msg))
end

## 1-dimensional ranges ##

"""
    AbstractRange{T}

Supertype for ranges with elements of type `T`.
[`UnitRange`](@ref) and other types are subtypes of this.
"""
abstract type AbstractRange{T} <: AbstractArray{T,1} end

RangeStepStyle(::Type{<:AbstractRange}) = RangeStepIrregular()
RangeStepStyle(::Type{<:AbstractRange{<:Integer}}) = RangeStepRegular()

convert(::Type{T}, r::AbstractRange) where {T<:AbstractRange} = r isa T ? r : T(r)

## ordinal ranges

"""
    OrdinalRange{T, S} <: AbstractRange{T}

Supertype for ordinal ranges with elements of type `T` with
spacing(s) of type `S`. The steps should be always-exact
multiples of [`oneunit`](@ref), and `T` should be a "discrete"
type, which cannot have values smaller than `oneunit`. For example,
`Integer` or `Date` types would qualify, whereas `Float64` would not (since this
type can represent values smaller than `oneunit(Float64)`.
[`UnitRange`](@ref), [`StepRange`](@ref), and other types are subtypes of this.
"""
abstract type OrdinalRange{T,S} <: AbstractRange{T} end

"""
    AbstractUnitRange{T} <: OrdinalRange{T, T}

Supertype for ranges with a step size of [`oneunit(T)`](@ref) with elements of type `T`.
[`UnitRange`](@ref) and other types are subtypes of this.
"""
abstract type AbstractUnitRange{T} <: OrdinalRange{T,T} end

"""
    StepRange{T, S} <: OrdinalRange{T, S}

Ranges with elements of type `T` with spacing of type `S`. The step
between each element is constant, and the range is defined in terms
of a `start` and `stop` of type `T` and a `step` of type `S`. Neither
`T` nor `S` should be floating point types. The syntax `a:b:c` with `b > 1`
and `a`, `b`, and `c` all integers creates a `StepRange`.

# Examples
```jldoctest
julia> collect(StepRange(1, Int8(2), 10))
5-element Vector{Int64}:
 1
 3
 5
 7
 9

julia> typeof(StepRange(1, Int8(2), 10))
StepRange{Int64, Int8}

julia> typeof(1:3:6)
StepRange{Int64, Int64}
```
"""
struct StepRange{T,S} <: OrdinalRange{T,S}
    start::T
    step::S
    stop::T

    function StepRange{T,S}(start, step, stop) where {T,S}
        sta = convert(T, start)
        ste = convert(S, step)
        sto = convert(T, stop)
        new(sta, ste, steprange_last(sta,ste,sto))
    end
end

# to make StepRange constructor inlineable, so optimizer can see `step` value
function steprange_last(start::T, step, stop) where T
    if isa(start,AbstractFloat) || isa(step,AbstractFloat)
        throw(ArgumentError("StepRange should not be used with floating point"))
    end
    if isa(start,Integer) && !isinteger(start + step)
        throw(ArgumentError("StepRange{<:Integer} cannot have non-integer step"))
    end
    z = zero(step)
    step == z && throw(ArgumentError("step cannot be zero"))

    if stop == start
        last = stop
    else
        if (step > z) != (stop > start)
            last = steprange_last_empty(start, step, stop)
        else
            # Compute absolute value of difference between `start` and `stop`
            # (to simplify handling both signed and unsigned T and checking for signed overflow):
            absdiff, absstep = stop > start ? (stop - start, step) : (start - stop, -step)

            # Compute remainder as a nonnegative number:
            if T <: Signed && absdiff < zero(absdiff)
                # handle signed overflow with unsigned rem
                remain = convert(T, unsigned(absdiff) % absstep)
            else
                remain = absdiff % absstep
            end
            # Move `stop` closer to `start` if there is a remainder:
            last = stop > start ? stop - remain : stop + remain
        end
    end
    last
end

function steprange_last_empty(start::Integer, step, stop)
    # empty range has a special representation where stop = start-1
    # this is needed to avoid the wrap-around that can happen computing
    # start - step, which leads to a range that looks very large instead
    # of empty.
    if step > zero(step)
        last = start - oneunit(stop-start)
    else
        last = start + oneunit(stop-start)
    end
    last
end
# For types where x+oneunit(x) may not be well-defined
steprange_last_empty(start, step, stop) = start - step

StepRange{T}(start, step::S, stop) where {T,S} = StepRange{T,S}(start, step, stop)
StepRange(start::T, step::S, stop::T) where {T,S} = StepRange{T,S}(start, step, stop)

"""
    UnitRange{T<:Real}

A range parameterized by a `start` and `stop` of type `T`, filled
with elements spaced by `1` from `start` until `stop` is exceeded.
The syntax `a:b` with `a` and `b` both `Integer`s creates a `UnitRange`.

# Examples
```jldoctest
julia> collect(UnitRange(2.3, 5.2))
3-element Vector{Float64}:
 2.3
 3.3
 4.3

julia> typeof(1:10)
UnitRange{Int64}
```
"""
struct UnitRange{T<:Real} <: AbstractUnitRange{T}
    start::T
    stop::T
    UnitRange{T}(start, stop) where {T<:Real} = new(start, unitrange_last(start,stop))
end
UnitRange(start::T, stop::T) where {T<:Real} = UnitRange{T}(start, stop)

unitrange_last(::Bool, stop::Bool) = stop
unitrange_last(start::T, stop::T) where {T<:Integer} =
    ifelse(stop >= start, stop, convert(T,start-oneunit(stop-start)))
unitrange_last(start::T, stop::T) where {T} =
    ifelse(stop >= start, convert(T,start+floor(stop-start)),
                          convert(T,start-oneunit(stop-start)))

unitrange(x) = UnitRange(x)

if isdefined(Main, :Base)
    # Constant-fold-able indexing into tuples to functionally expose Base.tail and Base.front
    function getindex(@nospecialize(t::Tuple), r::AbstractUnitRange)
        @inline
        require_one_based_indexing(r)
        if length(r) <= 10
            return ntuple(i -> t[i + first(r) - 1], length(r))
        elseif first(r) == 1
            last(r) == length(t)   && return t
            last(r) == length(t)-1 && return front(t)
            last(r) == length(t)-2 && return front(front(t))
        elseif last(r) == length(t)
            first(r) == 2 && return tail(t)
            first(r) == 3 && return tail(tail(t))
        end
        return (eltype(t)[t[ri] for ri in r]...,)
    end
end

"""
    Base.OneTo(n)

Define an `AbstractUnitRange` that behaves like `1:n`, with the added
distinction that the lower limit is guaranteed (by the type system) to
be 1.
"""
struct OneTo{T<:Integer} <: AbstractUnitRange{T}
    stop::T
    function OneTo{T}(stop) where {T<:Integer}
        throwbool(r)  = (@noinline; throw(ArgumentError("invalid index: $r of type Bool")))
        T === Bool && throwbool(stop)
        return new(max(zero(T), stop))
    end

    function OneTo{T}(r::AbstractRange) where {T<:Integer}
        throwstart(r) = (@noinline; throw(ArgumentError("first element must be 1, got $(first(r))")))
        throwstep(r)  = (@noinline; throw(ArgumentError("step must be 1, got $(step(r))")))
        throwbool(r)  = (@noinline; throw(ArgumentError("invalid index: $r of type Bool")))
        first(r) == 1 || throwstart(r)
        step(r)  == 1 || throwstep(r)
        T === Bool && throwbool(r)
        return new(max(zero(T), last(r)))
    end
end
OneTo(stop::T) where {T<:Integer} = OneTo{T}(stop)
OneTo(r::AbstractRange{T}) where {T<:Integer} = OneTo{T}(r)
oneto(r) = OneTo(r)

## Step ranges parameterized by length

"""
    StepRangeLen(         ref::R, step::S, len, [offset=1]) where {  R,S}
    StepRangeLen{T,R,S}(  ref::R, step::S, len, [offset=1]) where {T,R,S}
    StepRangeLen{T,R,S,L}(ref::R, step::S, len, [offset=1]) where {T,R,S,L}

A range `r` where `r[i]` produces values of type `T` (in the first
form, `T` is deduced automatically), parameterized by a `ref`erence
value, a `step`, and the `len`gth. By default `ref` is the starting
value `r[1]`, but alternatively you can supply it as the value of
`r[offset]` for some other index `1 <= offset <= len`. In conjunction
with `TwicePrecision` this can be used to implement ranges that are
free of roundoff error.
"""
struct StepRangeLen{T,R,S,L<:Integer} <: AbstractRange{T}
    ref::R       # reference value (might be smallest-magnitude value in the range)
    step::S      # step value
    len::L       # length of the range
    offset::L    # the index of ref

    function StepRangeLen{T,R,S,L}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S,L}
        if T <: Integer && !isinteger(ref + step)
            throw(ArgumentError("StepRangeLen{<:Integer} cannot have non-integer step"))
        end
        len = convert(L, len)
        len >= zero(len) || throw(ArgumentError("length cannot be negative, got $len"))
        offset = convert(L, offset)
        L1 = oneunit(typeof(len))
        L1 <= offset <= max(L1, len) || throw(ArgumentError("StepRangeLen: offset must be in [1,$len], got $offset"))
        return new(ref, step, len, offset)
    end
end

StepRangeLen{T,R,S}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S} =
    StepRangeLen{T,R,S,promote_type(Int,typeof(len))}(ref, step, len, offset)
StepRangeLen(ref::R, step::S, len::Integer, offset::Integer = 1) where {R,S} =
    StepRangeLen{typeof(ref+zero(step)),R,S,promote_type(Int,typeof(len))}(ref, step, len, offset)
StepRangeLen{T}(ref::R, step::S, len::Integer, offset::Integer = 1) where {T,R,S} =
    StepRangeLen{T,R,S,promote_type(Int,typeof(len))}(ref, step, len, offset)

## range with computed step

"""
    LinRange{T,L}

A range with `len` linearly spaced elements between its `start` and `stop`.
The size of the spacing is controlled by `len`, which must
be an `Integer`.

# Examples
```jldoctest
julia> LinRange(1.5, 5.5, 9)
9-element LinRange{Float64, Int64}:
 1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5
```

Compared to using [`range`](@ref), directly constructing a `LinRange` should
have less overhead but won't try to correct for floating point errors:
```jldoctest
julia> collect(range(-0.1, 0.3, length=5))
5-element Vector{Float64}:
 -0.1
  0.0
  0.1
  0.2
  0.3

julia> collect(LinRange(-0.1, 0.3, 5))
5-element Vector{Float64}:
 -0.1
 -1.3877787807814457e-17
  0.09999999999999999
  0.19999999999999998
  0.3
```
"""
struct LinRange{T,L<:Integer} <: AbstractRange{T}
    start::T
    stop::T
    len::L
    lendiv::L

    function LinRange{T,L}(start::T, stop::T, len::L) where {T,L<:Integer}
        len >= 0 || throw(ArgumentError("range($start, stop=$stop, length=$len): negative length"))
        onelen = oneunit(typeof(len))
        if len == onelen
            start == stop || throw(ArgumentError("range($start, stop=$stop, length=$len): endpoints differ"))
            return new(start, stop, len, len)
        end
        lendiv = max(len - onelen, onelen)
        if T <: Integer && !iszero(mod(stop-start, lendiv))
            throw(ArgumentError("LinRange{<:Integer} cannot have non-integer step"))
        end
        return new(start, stop, len, lendiv)
    end
end

function LinRange{T,L}(start, stop, len::Integer) where {T,L}
    LinRange{T,L}(convert(T, start), convert(T, stop), convert(L, len))
end

function LinRange{T}(start, stop, len::Integer) where T
    LinRange{T,promote_type(Int,typeof(len))}(start, stop, len)
end

function LinRange(start, stop, len::Integer)
    T = typeof((stop-start)/len)
    LinRange{T}(start, stop, len)
end

range_start_stop_length(start, stop, len::Integer) =
    range_start_stop_length(promote(start, stop)..., len)
range_start_stop_length(start::T, stop::T, len::Integer) where {T} = LinRange{T}(start, stop, len)
range_start_stop_length(start::T, stop::T, len::Integer) where {T<:Integer} =
    _linspace(float(T), start, stop, len)
## for Float16, Float32, and Float64 we hit twiceprecision.jl to lift to higher precision StepRangeLen
# for all other types we fall back to a plain old LinRange
_linspace(::Type{T}, start::Integer, stop::Integer, len::Integer) where T = LinRange{T}(start, stop, len)

function show(io::IO, r::LinRange)
    print(io, "range(")
    show(io, first(r))
    print(io, ", stop=")
    show(io, last(r))
    print(io, ", length=")
    show(io, length(r))
    print(io, ')')
end

"""
`print_range(io, r)` prints out a nice looking range r in terms of its elements
as if it were `collect(r)`, dependent on the size of the
terminal, and taking into account whether compact numbers should be shown.
It figures out the width in characters of each element, and if they
end up too wide, it shows the first and last elements separated by a
horizontal ellipsis. Typical output will look like `1.0,2.0,3.0,…,4.0,5.0,6.0`.

`print_range(io, r, pre, sep, post, hdots)` uses optional
parameters `pre` and `post` characters for each printed row,
`sep` separator string between printed elements,
`hdots` string for the horizontal ellipsis.
"""
function print_range(io::IO, r::AbstractRange,
                     pre::AbstractString = " ",
                     sep::AbstractString = ",",
                     post::AbstractString = "",
                     hdots::AbstractString = ",\u2026,") # horiz ellipsis
    # This function borrows from print_matrix() in show.jl
    # and should be called by show and display
    sz = displaysize(io)
    if !haskey(io, :compact)
        io = IOContext(io, :compact => true)
    end
    screenheight, screenwidth = sz[1] - 4, sz[2]
    screenwidth -= length(pre) + length(post)
    postsp = ""
    sepsize = length(sep)
    m = 1 # treat the range as a one-row matrix
    n = length(r)
    # Figure out spacing alignments for r, but only need to examine the
    # left and right edge columns, as many as could conceivably fit on the
    # screen, with the middle columns summarized by horz, vert, or diag ellipsis
    maxpossiblecols = div(screenwidth, 1+sepsize) # assume each element is at least 1 char + 1 separator
    colsr = n <= maxpossiblecols ? (1:n) : [1:div(maxpossiblecols,2)+1; (n-div(maxpossiblecols,2)):n]
    rowmatrix = reshape(r[colsr], 1, length(colsr)) # treat the range as a one-row matrix for print_matrix_row
    nrow, idxlast = size(rowmatrix, 2), last(axes(rowmatrix, 2))
    A = alignment(io, rowmatrix, 1:m, 1:length(rowmatrix), screenwidth, screenwidth, sepsize, nrow) # how much space range takes
    if n <= length(A) # cols fit screen, so print out all elements
        print(io, pre) # put in pre chars
        print_matrix_row(io,rowmatrix,A,1,1:n,sep,idxlast) # the entire range
        print(io, post) # add the post characters
    else # cols don't fit so put horiz ellipsis in the middle
        # how many chars left after dividing width of screen in half
        # and accounting for the horiz ellipsis
        c = div(screenwidth-length(hdots)+1,2)+1 # chars remaining for each side of rowmatrix
        alignR = reverse(alignment(io, rowmatrix, 1:m, length(rowmatrix):-1:1, c, c, sepsize, nrow)) # which cols of rowmatrix to put on the right
        c = screenwidth - sum(map(sum,alignR)) - (length(alignR)-1)*sepsize - length(hdots)
        alignL = alignment(io, rowmatrix, 1:m, 1:length(rowmatrix), c, c, sepsize, nrow) # which cols of rowmatrix to put on the left
        print(io, pre)   # put in pre chars
        print_matrix_row(io, rowmatrix,alignL,1,1:length(alignL),sep,idxlast) # left part of range
        print(io, hdots) # horizontal ellipsis
        print_matrix_row(io, rowmatrix,alignR,1,length(rowmatrix)-length(alignR)+1:length(rowmatrix),sep,idxlast) # right part of range
        print(io, post)  # post chars
    end
end

## interface implementations

length(r::AbstractRange) = error("length implementation missing") # catch mistakes
size(r::AbstractRange) = (length(r),)

isempty(r::StepRange) =
    # steprange_last_empty(r.start, r.step, r.stop) == r.stop
    (r.start != r.stop) & ((r.step > zero(r.step)) != (r.stop > r.start))
isempty(r::AbstractUnitRange) = first(r) > last(r)
isempty(r::StepRangeLen) = length(r) == 0
isempty(r::LinRange) = length(r) == 0

"""
    step(r)

Get the step size of an [`AbstractRange`](@ref) object.

# Examples
```jldoctest
julia> step(1:10)
1

julia> step(1:2:10)
2

julia> step(2.5:0.3:10.9)
0.3

julia> step(range(2.5, stop=10.9, length=85))
0.1
```
"""
step(r::StepRange) = r.step
step(r::AbstractUnitRange{T}) where {T} = oneunit(T) - zero(T)
step(r::StepRangeLen) = r.step
step(r::StepRangeLen{T}) where {T<:AbstractFloat} = T(r.step)
step(r::LinRange) = (last(r)-first(r))/r.lendiv

# high-precision step
step_hp(r::StepRangeLen) = r.step
step_hp(r::AbstractRange) = step(r)

axes(r::AbstractRange) = (oneto(length(r)),)

# Needed to fold the `firstindex` call in SimdLoop.simd_index
firstindex(::UnitRange) = 1
firstindex(::StepRange) = 1
firstindex(::LinRange) = 1

# n.b. checked_length for these is defined iff checked_add and checked_sub are
# defined between the relevant types
function checked_length(r::OrdinalRange{T}) where T
    s = step(r)
    # s != 0, by construction, but avoids the division error later
    start = first(r)
    if s == zero(s) || isempty(r)
        return Integer(start - start + zero(s))
    end
    stop = last(r)
    if isless(s, zero(s))
        diff = checked_sub(start, stop)
        s = -s
    else
        diff = checked_sub(stop, start)
    end
    a = div(diff, s)
    return Integer(checked_add(a, oneunit(a)))
end

function checked_length(r::AbstractUnitRange{T}) where T
    # compiler optimization: remove dead cases from above
    if isempty(r)
        return Integer(first(r) - first(r))
    end
    a = checked_sub(last(r), first(r))
    return Integer(checked_add(a, oneunit(a)))
end

function length(r::OrdinalRange{T}) where T
    s = step(r)
    # s != 0, by construction, but avoids the division error later
    start = first(r)
    if s == zero(s) || isempty(r)
        return Integer(div(start-start, oneunit(s)))
    end
    stop = last(r)
    if isless(s, zero(s))
        diff = start - stop
        s = -s
    else
        diff = stop - start
    end
    a = div(diff, s)
    return Integer(a + oneunit(a))
end

function length(r::AbstractUnitRange{T}) where T
    @inline
    a = last(r) - first(r) # even when isempty, by construction (with overflow)
    return Integer(a + oneunit(a))
end

length(r::OneTo) = Integer(r.stop - zero(r.stop))
length(r::StepRangeLen) = r.len
length(r::LinRange) = r.len

let bigints = Union{Int, UInt, Int64, UInt64, Int128, UInt128}
    global length, checked_length
    # compile optimization for which promote_type(T, Int) == T
    length(r::OneTo{T}) where {T<:bigints} = r.stop
    # slightly more accurate length and checked_length in extreme cases
    # (near typemax) for types with known `unsigned` functions
    function length(r::OrdinalRange{T}) where T<:bigints
        s = step(r)
        s == zero(s) && return zero(T) # unreachable, by construction, but avoids the error case here later
        isempty(r) && return zero(T)
        diff = last(r) - first(r)
        # if |s| > 1, diff might have overflowed, but unsigned(diff)÷s should
        # therefore still be valid (if the result is representable at all)
        # n.b. !(s isa T)
        if s isa Unsigned || -1 <= s <= 1 || s == -s
            a = div(diff, s)
        elseif s < 0
            a = div(unsigned(-diff), -s) % typeof(diff)
        else
            a = div(unsigned(diff), s) % typeof(diff)
        end
        return Integer(a) + oneunit(a)
    end
    function checked_length(r::OrdinalRange{T}) where T<:bigints
        s = step(r)
        s == zero(s) && return zero(T) # unreachable, by construction, but avoids the error case here later
        isempty(r) && return zero(T)
        stop, start = last(r), first(r)
        # n.b. !(s isa T)
        if s > 1
            diff = stop - start
            a = convert(T, div(unsigned(diff), s))
        elseif s < -1
            diff = start - stop
            a = convert(T, div(unsigned(diff), -s))
        elseif s > 0
            a = div(checked_sub(stop, start), s)
        else
            a = div(checked_sub(start, stop), -s)
        end
        return checked_add(a, oneunit(a))
    end
end

# some special cases to favor default Int type
let smallints = (Int === Int64 ?
                Union{Int8, UInt8, Int16, UInt16, Int32, UInt32} :
                Union{Int8, UInt8, Int16, UInt16})
    global length, checked_length
    # n.b. !(step isa T)
    length(r::OrdinalRange{<:smallints}) = div(Int(last(r)) - Int(first(r)), step(r)) + 1
    length(r::AbstractUnitRange{<:smallints}) = Int(last(r)) - Int(first(r)) + 1
    length(r::OneTo{<:smallints}) = Int(r.stop)
    checked_length(r::OrdinalRange{<:smallints}) = length(r)
    checked_length(r::AbstractUnitRange{<:smallints}) = length(r)
    checked_length(r::OneTo{<:smallints}) = length(r)
end

first(r::OrdinalRange{T}) where {T} = convert(T, r.start)
first(r::OneTo{T}) where {T} = oneunit(T)
first(r::StepRangeLen) = unsafe_getindex(r, 1)
first(r::LinRange) = r.start

last(r::OrdinalRange{T}) where {T} = convert(T, r.stop)
last(r::StepRangeLen) = unsafe_getindex(r, length(r))
last(r::LinRange) = r.stop

minimum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : first(r)
maximum(r::AbstractUnitRange) = isempty(r) ? throw(ArgumentError("range must be non-empty")) : last(r)
minimum(r::AbstractRange)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : min(first(r), last(r))
maximum(r::AbstractRange)  = isempty(r) ? throw(ArgumentError("range must be non-empty")) : max(first(r), last(r))

"""
    argmin(r::AbstractRange)

Ranges can have multiple minimal elements. In that case
`argmin` will return a minimal index, but not necessarily the
first one.
"""
function argmin(r::AbstractRange)
    if isempty(r)
        throw(ArgumentError("range must be non-empty"))
    elseif step(r) > 0
        firstindex(r)
    else
        lastindex(r)
    end
end

"""
    argmax(r::AbstractRange)

Ranges can have multiple maximal elements. In that case
`argmax` will return a maximal index, but not necessarily the
first one.
"""
function argmax(r::AbstractRange)
    if isempty(r)
        throw(ArgumentError("range must be non-empty"))
    elseif step(r) > 0
        lastindex(r)
    else
        firstindex(r)
    end
end

extrema(r::AbstractRange) = (minimum(r), maximum(r))

# Ranges are immutable
copy(r::AbstractRange) = r

## iteration

function iterate(r::Union{StepRangeLen,LinRange}, i::Integer=zero(length(r)))
    @inline
    i += oneunit(i)
    length(r) < i && return nothing
    unsafe_getindex(r, i), i
end

iterate(r::OrdinalRange) = isempty(r) ? nothing : (first(r), first(r))

function iterate(r::OrdinalRange{T}, i) where {T}
    @inline
    i == last(r) && return nothing
    next = convert(T, i + step(r))
    (next, next)
end

## indexing

_in_unit_range(v::UnitRange, val, i::Integer) = i > 0 && val <= v.stop && val >= v.start

function getindex(v::UnitRange{T}, i::Integer) where T
    @inline
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    val = convert(T, v.start + (i - 1))
    @boundscheck _in_unit_range(v, val, i) || throw_boundserror(v, i)
    val
end

const OverflowSafe = Union{Bool,Int8,Int16,Int32,Int64,Int128,
                           UInt8,UInt16,UInt32,UInt64,UInt128}

function getindex(v::UnitRange{T}, i::Integer) where {T<:OverflowSafe}
    @inline
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    val = v.start + (i - 1)
    @boundscheck _in_unit_range(v, val, i) || throw_boundserror(v, i)
    val % T
end

function getindex(v::OneTo{T}, i::Integer) where T
    @inline
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    @boundscheck ((i > 0) & (i <= v.stop)) || throw_boundserror(v, i)
    convert(T, i)
end

function getindex(v::AbstractRange{T}, i::Integer) where T
    @inline
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    ret = convert(T, first(v) + (i - 1)*step_hp(v))
    ok = ifelse(step(v) > zero(step(v)),
                (ret <= last(v)) & (ret >= first(v)),
                (ret <= first(v)) & (ret >= last(v)))
    @boundscheck ((i > 0) & ok) || throw_boundserror(v, i)
    ret
end

function getindex(r::Union{StepRangeLen,LinRange}, i::Integer)
    @inline
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    @boundscheck checkbounds(r, i)
    unsafe_getindex(r, i)
end

# This is separate to make it useful even when running with --check-bounds=yes
function unsafe_getindex(r::StepRangeLen{T}, i::Integer) where T
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    u = i - r.offset
    T(r.ref + u*r.step)
end

function _getindex_hiprec(r::StepRangeLen, i::Integer)  # without rounding by T
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    u = i - r.offset
    r.ref + u*r.step
end

function unsafe_getindex(r::LinRange, i::Integer)
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    lerpi(i-1, r.lendiv, r.start, r.stop)
end

function lerpi(j::Integer, d::Integer, a::T, b::T) where T
    @inline
    t = j/d
    T((1-t)*a + t*b)
end

getindex(r::AbstractRange, ::Colon) = copy(r)

function getindex(r::AbstractUnitRange, s::AbstractUnitRange{T}) where {T<:Integer}
    @inline
    @boundscheck checkbounds(r, s)

    if T === Bool
        range(first(s) ? first(r) : last(r), length = Integer(last(s)))
    else
        f = first(r)
        st = oftype(f, f + first(s)-firstindex(r))
        return range(st, length=length(s))
    end
end

function getindex(r::OneTo{T}, s::OneTo) where T
    @inline
    @boundscheck checkbounds(r, s)
    OneTo(T(s.stop))
end

function getindex(r::AbstractUnitRange, s::StepRange{T}) where {T<:Integer}
    @inline
    @boundscheck checkbounds(r, s)

    if T === Bool
        range(first(s) ? first(r) : last(r), step=oneunit(eltype(r)), length = Integer(last(s)))
    else
        st = oftype(first(r), first(r) + s.start-firstindex(r))
        return range(st, step=step(s), length=length(s))
    end
end

function getindex(r::StepRange, s::AbstractRange{T}) where {T<:Integer}
    @inline
    @boundscheck checkbounds(r, s)

    if T === Bool
        if length(s) == 0
            return range(first(r), step=step(r), length=0)
        elseif length(s) == 1
            if first(s)
                return range(first(r), step=step(r), length=1)
            else
                return range(first(r), step=step(r), length=0)
            end
        else # length(s) == 2
            return range(last(r), step=step(r), length=1)
        end
    else
        st = oftype(r.start, r.start + (first(s)-1)*step(r))
        return range(st, step=step(r)*step(s), length=length(s))
    end
end

function getindex(r::StepRangeLen{T}, s::OrdinalRange{S}) where {T, S<:Integer}
    @inline
    @boundscheck checkbounds(r, s)

    len = length(s)
    sstep = step_hp(s)
    rstep = step_hp(r)
    L = typeof(len)
    if S === Bool
        rstep *= one(sstep)
        if len == 0
            return StepRangeLen{T}(first(r), rstep, zero(L), oneunit(L))
        elseif len == 1
            if first(s)
                return StepRangeLen{T}(first(r), rstep, oneunit(L), oneunit(L))
            else
                return StepRangeLen{T}(first(r), rstep, zero(L), oneunit(L))
            end
        else # len == 2
            return StepRangeLen{T}(last(r), rstep, oneunit(L), oneunit(L))
        end
    else
        # Find closest approach to offset by s
        ind = LinearIndices(s)
        offset = L(max(min(1 + round(L, (r.offset - first(s))/sstep), last(ind)), first(ind)))
        ref = _getindex_hiprec(r, first(s) + (offset-1)*sstep)
        return StepRangeLen{T}(ref, rstep*sstep, len, offset)
    end
end

function getindex(r::LinRange{T}, s::OrdinalRange{S}) where {T, S<:Integer}
    @inline
    @boundscheck checkbounds(r, s)

    len = length(s)
    L = typeof(len)
    if S === Bool
        if len == 0
            return LinRange{T}(first(r), first(r), len)
        elseif len == 1
            if first(s)
                return LinRange{T}(first(r), first(r), len)
            else
                return LinRange{T}(first(r), first(r), zero(L))
            end
        else # length(s) == 2
            return LinRange{T}(last(r), last(r), oneunit(L))
        end
    else
        vfirst = unsafe_getindex(r, first(s))
        vlast  = unsafe_getindex(r, last(s))
        return LinRange{T}(vfirst, vlast, len)
    end
end

show(io::IO, r::AbstractRange) = print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r)))
show(io::IO, r::UnitRange) = print(io, repr(first(r)), ':', repr(last(r)))
show(io::IO, r::OneTo) = print(io, "Base.OneTo(", r.stop, ")")
function show(io::IO, r::StepRangeLen)
    if step(r) != 0
        print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r)))
    else
        # ugly temporary printing, to avoid 0:0:0 etc.
        print(io, "StepRangeLen(", repr(first(r)), ", ", repr(step(r)), ", ", repr(length(r)), ")")
    end
end

function ==(r::T, s::T) where {T<:AbstractRange}
    isempty(r) && return isempty(s)
    _has_length_one(r) && return _has_length_one(s) & (first(r) == first(s))
    (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
end

function ==(r::OrdinalRange, s::OrdinalRange)
    isempty(r) && return isempty(s)
    _has_length_one(r) && return _has_length_one(s) & (first(r) == first(s))
    (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
end

==(r::AbstractUnitRange, s::AbstractUnitRange) =
    (isempty(r) & isempty(s)) | ((first(r) == first(s)) & (last(r) == last(s)))

==(r::OneTo, s::OneTo) = last(r) == last(s)

==(r::T, s::T) where {T<:Union{StepRangeLen,LinRange}} =
    (isempty(r) & isempty(s)) | ((first(r) == first(s)) & (length(r) == length(s)) & (last(r) == last(s)))

function ==(r::Union{StepRange{T},StepRangeLen{T,T}}, s::Union{StepRange{T},StepRangeLen{T,T}}) where {T}
    isempty(r) && return isempty(s)
    _has_length_one(r) && return _has_length_one(s) & (first(r) == first(s))
    (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))
end

_has_length_one(r::OrdinalRange) = first(r) == last(r)
_has_length_one(r::AbstractRange) = isone(length(r))

function ==(r::AbstractRange, s::AbstractRange)
    lr = length(r)
    if lr != length(s)
        return false
    elseif iszero(lr)
        return true
    end
    yr, ys = iterate(r), iterate(s)
    while yr !== nothing
        yr[1] == ys[1] || return false
        yr, ys = iterate(r, yr[2]), iterate(s, ys[2])
    end
    return true
end

intersect(r::OneTo, s::OneTo) = OneTo(min(r.stop,s.stop))
union(r::OneTo, s::OneTo) = OneTo(max(r.stop,s.stop))

intersect(r::AbstractUnitRange{<:Integer}, s::AbstractUnitRange{<:Integer}) = max(first(r),first(s)):min(last(r),last(s))

intersect(i::Integer, r::AbstractUnitRange{<:Integer}) = range(max(i, first(r)), length=in(i, r))

intersect(r::AbstractUnitRange{<:Integer}, i::Integer) = intersect(i, r)

function intersect(r::AbstractUnitRange{<:Integer}, s::StepRange{<:Integer})
    T = promote_type(eltype(r), eltype(s))
    if isempty(s)
        StepRange{T}(first(r), +step(s), first(r)-step(s))
    else
        sta, ste, sto = first_step_last_ascending(s)
        lo = first(r)
        hi = last(r)
        i0 = max(sta, lo + mod(sta - lo, ste))
        i1 = min(sto, hi - mod(hi - sta, ste))
        StepRange{T}(i0, ste, i1)
    end
end

function first_step_last_ascending(r::StepRange)
    if step(r) < zero(step(r))
        last(r), -step(r), first(r)
    else
        first(r), +step(r), last(r)
    end
end

function intersect(r::StepRange{<:Integer}, s::AbstractUnitRange{<:Integer})
    if step(r) < 0
        reverse(intersect(s, reverse(r)))
    else
        intersect(s, r)
    end
end

function intersect(r::StepRange, s::StepRange)
    T = promote_type(eltype(r), eltype(s))
    S = promote_type(typeof(step(r)), typeof(step(s)))
    if isempty(r) || isempty(s)
        return StepRange{T,S}(first(r), step(r), first(r)-step(r))
    end

    start1, step1, stop1 = first_step_last_ascending(r)
    start2, step2, stop2 = first_step_last_ascending(s)
    a = lcm(step1, step2)

    g, x, y = gcdx(step1, step2)

    if !iszero(rem(start1 - start2, g))
        # Unaligned, no overlap possible.
        if  step(r) < zero(step(r))
            return StepRange{T,S}(stop1, -a, stop1+a)
        else
            return StepRange{T,S}(start1, a, start1-a)
        end
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    n = min(stop1 - mod(stop1 - b, a), stop2 - mod(stop2 - b, a))
    step(r) < zero(step(r)) ? StepRange{T,S}(n, -a, m) : StepRange{T,S}(m, a, n)
end

function intersect(r1::AbstractRange, r2::AbstractRange, r3::AbstractRange, r::AbstractRange...)
    i = intersect(intersect(r1, r2), r3)
    for t in r
        i = intersect(i, t)
    end
    i
end

# _findin (the index of intersection)
function _findin(r::AbstractRange{<:Integer}, span::AbstractUnitRange{<:Integer})
    local ifirst
    local ilast
    fspan = first(span)
    lspan = last(span)
    fr = first(r)
    lr = last(r)
    sr = step(r)
    if sr > 0
        ifirst = fr >= fspan ? 1 : ceil(Integer,(fspan-fr)/sr)+1
        ilast = lr <= lspan ? length(r) : length(r) - ceil(Integer,(lr-lspan)/sr)
    elseif sr < 0
        ifirst = fr <= lspan ? 1 : ceil(Integer,(lspan-fr)/sr)+1
        ilast = lr >= fspan ? length(r) : length(r) - ceil(Integer,(lr-fspan)/sr)
    else
        ifirst = fr >= fspan ? 1 : length(r)+1
        ilast = fr <= lspan ? length(r) : 0
    end
    r isa AbstractUnitRange ? (ifirst:ilast) : (ifirst:1:ilast)
end

issubset(r::OneTo, s::OneTo) = r.stop <= s.stop

issubset(r::AbstractUnitRange{<:Integer}, s::AbstractUnitRange{<:Integer}) =
    isempty(r) || first(r) >= first(s) && last(r) <= last(s)

## linear operations on ranges ##

-(r::OrdinalRange) = range(-first(r), step=-step(r), length=length(r))
-(r::StepRangeLen{T,R,S,L}) where {T,R,S,L} =
    StepRangeLen{T,R,S,L}(-r.ref, -r.step, r.len, r.offset)
function -(r::LinRange)
    start = -r.start
    LinRange{typeof(start)}(start, -r.stop, length(r))
end

# promote eltype if at least one container wouldn't change, otherwise join container types.
el_same(::Type{T}, a::Type{<:AbstractArray{T,n}}, b::Type{<:AbstractArray{T,n}}) where {T,n}   = a
el_same(::Type{T}, a::Type{<:AbstractArray{T,n}}, b::Type{<:AbstractArray{S,n}}) where {T,S,n} = a
el_same(::Type{T}, a::Type{<:AbstractArray{S,n}}, b::Type{<:AbstractArray{T,n}}) where {T,S,n} = b
el_same(::Type, a, b) = promote_typejoin(a, b)

promote_rule(a::Type{UnitRange{T1}}, b::Type{UnitRange{T2}}) where {T1,T2} =
    el_same(promote_type(T1, T2), a, b)
UnitRange{T}(r::UnitRange{T}) where {T<:Real} = r
UnitRange{T}(r::UnitRange) where {T<:Real} = UnitRange{T}(r.start, r.stop)

promote_rule(a::Type{OneTo{T1}}, b::Type{OneTo{T2}}) where {T1,T2} =
    el_same(promote_type(T1, T2), a, b)
OneTo{T}(r::OneTo{T}) where {T<:Integer} = r
OneTo{T}(r::OneTo) where {T<:Integer} = OneTo{T}(r.stop)

promote_rule(a::Type{UnitRange{T1}}, ::Type{UR}) where {T1,UR<:AbstractUnitRange} =
    promote_rule(a, UnitRange{eltype(UR)})
UnitRange{T}(r::AbstractUnitRange) where {T<:Real} = UnitRange{T}(first(r), last(r))
UnitRange(r::AbstractUnitRange) = UnitRange(first(r), last(r))

AbstractUnitRange{T}(r::AbstractUnitRange{T}) where {T} = r
AbstractUnitRange{T}(r::UnitRange) where {T} = UnitRange{T}(r)
AbstractUnitRange{T}(r::OneTo) where {T} = OneTo{T}(r)

OrdinalRange{T1, T2}(r::StepRange) where {T1, T2<: Integer} = StepRange{T1, T2}(r)
OrdinalRange{T1, T2}(r::AbstractUnitRange{T1}) where {T1, T2<:Integer} = r
OrdinalRange{T1, T2}(r::UnitRange) where {T1, T2<:Integer} = UnitRange{T1}(r)
OrdinalRange{T1, T2}(r::OneTo) where {T1, T2<:Integer} = OneTo{T1}(r)

function promote_rule(::Type{StepRange{T1a,T1b}}, ::Type{StepRange{T2a,T2b}}) where {T1a,T1b,T2a,T2b}
    Tb = promote_type(T1b, T2b)
    # el_same only operates on array element type, so just promote second type parameter
    el_same(promote_type(T1a, T2a), StepRange{T1a,Tb}, StepRange{T2a,Tb})
end
StepRange{T1,T2}(r::StepRange{T1,T2}) where {T1,T2} = r

promote_rule(a::Type{StepRange{T1a,T1b}}, ::Type{UR}) where {T1a,T1b,UR<:AbstractUnitRange} =
    promote_rule(a, StepRange{eltype(UR), eltype(UR)})
StepRange{T1,T2}(r::AbstractRange) where {T1,T2} =
    StepRange{T1,T2}(convert(T1, first(r)), convert(T2, step(r)), convert(T1, last(r)))
StepRange(r::AbstractUnitRange{T}) where {T} =
    StepRange{T,T}(first(r), step(r), last(r))
(StepRange{T1,T2} where T1)(r::AbstractRange) where {T2} = StepRange{eltype(r),T2}(r)

function promote_rule(::Type{StepRangeLen{T1,R1,S1,L1}},::Type{StepRangeLen{T2,R2,S2,L2}}) where {T1,T2,R1,R2,S1,S2,L1,L2}
    R, S, L = promote_type(R1, R2), promote_type(S1, S2), promote_type(L1, L2)
    el_same(promote_type(T1, T2), StepRangeLen{T1,R,S,L}, StepRangeLen{T2,R,S,L})
end
StepRangeLen{T,R,S,L}(r::StepRangeLen{T,R,S,L}) where {T,R,S,L} = r
StepRangeLen{T,R,S,L}(r::StepRangeLen) where {T,R,S,L} =
    StepRangeLen{T,R,S,L}(convert(R, r.ref), convert(S, r.step), convert(L, r.len), convert(L, r.offset))
StepRangeLen{T}(r::StepRangeLen) where {T} =
    StepRangeLen(convert(T, r.ref), convert(T, r.step), r.len, r.offset)

promote_rule(a::Type{StepRangeLen{T,R,S,L}}, ::Type{OR}) where {T,R,S,L,OR<:AbstractRange} =
    promote_rule(a, StepRangeLen{eltype(OR), eltype(OR), eltype(OR), Int})
StepRangeLen{T,R,S,L}(r::AbstractRange) where {T,R,S,L} =
    StepRangeLen{T,R,S,L}(R(first(r)), S(step(r)), length(r))
StepRangeLen{T}(r::AbstractRange) where {T} =
    StepRangeLen(T(first(r)), T(step(r)), length(r))
StepRangeLen(r::AbstractRange) = StepRangeLen{eltype(r)}(r)

function promote_rule(a::Type{LinRange{T1,L1}}, b::Type{LinRange{T2,L2}}) where {T1,T2,L1,L2}
    L = promote_type(L1, L2)
    el_same(promote_type(T1, T2), LinRange{T1,L}, LinRange{T2,L})
end
LinRange{T,L}(r::LinRange{T,L}) where {T,L} = r
LinRange{T,L}(r::AbstractRange) where {T,L} = LinRange{T,L}(first(r), last(r), length(r))
LinRange{T}(r::AbstractRange) where {T} = LinRange{T,typeof(length(r))}(first(r), last(r), length(r))
LinRange(r::AbstractRange{T}) where {T} = LinRange{T}(r)

promote_rule(a::Type{LinRange{T,L}}, ::Type{OR}) where {T,L,OR<:OrdinalRange} =
    promote_rule(a, LinRange{eltype(OR),L})

promote_rule(::Type{LinRange{A,L}}, b::Type{StepRangeLen{T2,R2,S2,L2}}) where {A,L,T2,R2,S2,L2} =
    promote_rule(StepRangeLen{A,A,A,L}, b)

## concatenation ##

function vcat(rs::AbstractRange{T}...) where T
    n::Int = 0
    for ra in rs
        n += length(ra)
    end
    a = Vector{T}(undef, n)
    i = 1
    for ra in rs, x in ra
        @inbounds a[i] = x
        i += 1
    end
    return a
end

Array{T,1}(r::AbstractRange{T}) where {T} = vcat(r)
collect(r::AbstractRange) = vcat(r)

_reverse(r::OrdinalRange, ::Colon) = (:)(last(r), -step(r), first(r))
function _reverse(r::StepRangeLen, ::Colon)
    # If `r` is empty, `length(r) - r.offset + 1 will be nonpositive hence
    # invalid. As `reverse(r)` is also empty, any offset would work so we keep
    # `r.offset`
    offset = isempty(r) ? r.offset : length(r)-r.offset+1
    return typeof(r)(r.ref, -r.step, length(r), offset)
end
_reverse(r::LinRange{T}, ::Colon) where {T} = typeof(r)(r.stop, r.start, length(r))

## sorting ##

issorted(r::AbstractUnitRange) = true
issorted(r::AbstractRange) = length(r) <= 1 || step(r) >= zero(step(r))

sort(r::AbstractUnitRange) = r
sort!(r::AbstractUnitRange) = r

sort(r::AbstractRange) = issorted(r) ? r : reverse(r)

sortperm(r::AbstractUnitRange) = 1:length(r)
sortperm(r::AbstractRange) = issorted(r) ? (1:1:length(r)) : (length(r):-1:1)

function sum(r::AbstractRange{<:Real})
    l = length(r)
    # note that a little care is required to avoid overflow in l*(l-1)/2
    return l * first(r) + (iseven(l) ? (step(r) * (l-1)) * (l>>1)
                                     : (step(r) * l) * ((l-1)>>1))
end

function _in_range(x, r::AbstractRange)
    if !isfinite(x)
        return false
    elseif iszero(step(r))
        return !isempty(r) && first(r) == x
    else
        n = round(Integer, (x - first(r)) / step(r)) + 1
        return n >= 1 && n <= length(r) && r[n] == x
    end
end
in(x::Real, r::AbstractRange{<:Real}) = _in_range(x, r)
# This method needs to be defined separately since -(::T, ::T) can be implemented
# even if -(::T, ::Real) is not
in(x::T, r::AbstractRange{T}) where {T} = _in_range(x, r)

in(x::Integer, r::AbstractUnitRange{<:Integer}) = (first(r) <= x) & (x <= last(r))

in(x::Real, r::AbstractRange{T}) where {T<:Integer} =
    isinteger(x) && !isempty(r) && x >= minimum(r) && x <= maximum(r) &&
        (mod(convert(T,x),step(r))-mod(first(r),step(r)) == 0)
in(x::AbstractChar, r::AbstractRange{<:AbstractChar}) =
    !isempty(r) && x >= minimum(r) && x <= maximum(r) &&
        (mod(Int(x) - Int(first(r)), step(r)) == 0)

# Addition/subtraction of ranges

function _define_range_op(@nospecialize f)
    @eval begin
        function $f(r1::OrdinalRange, r2::OrdinalRange)
            r1l = length(r1)
            (r1l == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match: length of r1 is $r1l, length of r2 is $(length(r2))")))
            StepRangeLen($f(first(r1), first(r2)), $f(step(r1), step(r2)), r1l)
        end

        function $f(r1::LinRange{T}, r2::LinRange{T}) where T
            len = r1.len
            (len == r2.len ||
             throw(DimensionMismatch("argument dimensions must match: length of r1 is $len, length of r2 is $(r2.len)")))
            LinRange{T}(convert(T, $f(first(r1), first(r2))),
                        convert(T, $f(last(r1), last(r2))), len)
        end

        $f(r1::Union{StepRangeLen, OrdinalRange, LinRange},
           r2::Union{StepRangeLen, OrdinalRange, LinRange}) =
               $f(promote(r1, r2)...)
    end
end
_define_range_op(:+)
_define_range_op(:-)

function +(r1::StepRangeLen{T,S}, r2::StepRangeLen{T,S}) where {T,S}
    len = length(r1)
    (len == length(r2) ||
     throw(DimensionMismatch("argument dimensions must match: length of r1 is $len, length of r2 is $(length(r2))")))
    StepRangeLen(first(r1)+first(r2), step(r1)+step(r2), len)
end

-(r1::StepRangeLen, r2::StepRangeLen) = +(r1, -r2)

# Modular arithmetic on ranges

"""
    mod(x::Integer, r::AbstractUnitRange)

Find `y` in the range `r` such that ``x ≡ y (mod n)``, where `n = length(r)`,
i.e. `y = mod(x - first(r), n) + first(r)`.

See also [`mod1`](@ref).

# Examples
```jldoctest
julia> mod(0, Base.OneTo(3))  # mod1(0, 3)
3

julia> mod(3, 0:2)  # mod(3, 3)
0
```

!!! compat "Julia 1.3"
     This method requires at least Julia 1.3.
"""
mod(i::Integer, r::OneTo) = mod1(i, last(r))
mod(i::Integer, r::AbstractUnitRange{<:Integer}) = mod(i-first(r), length(r)) + first(r)

max_stars_repo_name: RelationalAI-oss/julia

max_issues_repo_name: RelationalAI-oss/julia

max_forks_repo_name: RelationalAI-oss/julia


text:
# shows some data analysis using DataFrames
using CSV, StatsModels, DataFrames, StatPlots
nerlove = CSV.read("nerlove.csv")
nerlove[:lnC] = log.(nerlove[:cost])
nerlove[:lnQ] = log.(nerlove[:output])
nerlove[:lnPL] = log.(nerlove[:labor])
nerlove[:lnPF] = log.(nerlove[:fuel])
nerlove[:lnPK] = log.(nerlove[:capital])
show(describe(nerlove))
@df nerlove scatter(:lnC, :lnQ)
gui()
f = @formula(lnC ~ 1 + lnQ + lnPL + lnPF + lnPK)
ols(f, nerlove)


text:
# This file was generated, do not modify it.

using MLJ, RDatasets, ScientificTypes,
      DataFrames, Statistics, StatsBase

smarket = dataset("ISLR", "Smarket")
@show size(smarket)
@show names(smarket)

describe(smarket, :mean, :std, :eltype)

y = smarket.Direction
X = select(smarket, Not(:Direction));

cm = X |> Matrix |> cor
round.(cm, sigdigits=1)

using PyPlot
figure(figsize=(8,6))
plot(X.Volume)
xlabel("Tick number", fontsize=14)
ylabel("Volume", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig("assets/ISL-volume.svg") # hide

yc = coerce(y, Multiclass)
unique(yc)

@load LogisticClassifier pkg=MLJLinearModels
X2 = select(X, Not([:Year, :Today]))
clf = machine(LogisticClassifier(), X2, y)

fit!(clf)
ŷ = predict(clf, X2)
cross_entropy(ŷ, y) |> mean

ŷ = predict_mode(clf, X2)
misclassification_rate(ŷ, y)

TN = down_down = sum(ŷ .== y .== "Down")
FN = down_up = sum(ŷ .!= y .== "Up")
FP = up_down = sum(ŷ .!= y .== "Down")
TP = up_up = sum(ŷ .== y .== "Up")

conf_mat = [down_down down_up; up_down up_up]

acc = (TN + TP) / length(y)
prec = TP /  (TP + FP)
rec  = TP / (TP + FN)
@show round(acc, sigdigits=3)
@show round(prec, sigdigits=3)
@show round(rec, sigdigits=3)

train = 1:findlast(X.Year .< 2005);
test = last(train)+1:length(y)

fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mcr = misclassification_rate(ŷ, y[test])
accuracy = 1 - mcr

X3 = select(X2, [:Lag1, :Lag2])
clf = machine(LogisticClassifier(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
mean(ŷ .== y[test])

Xnew = (Lag1 = [1.2, 1.5], Lag2 = [1.1, -0.8])
@show ŷ = predict(clf, Xnew)

mode.(ŷ)

@load BayesianLDA pkg=MultivariateStats

clf = machine(BayesianLDA(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

@load LDA pkg=MultivariateStats
using Distances

clf = machine(LDA(dist=CosineDist()), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

@load BayesianQDA pkg=ScikitLearn

clf = machine(BayesianQDA(), X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

acc = mean(ŷ .== y[test])

@load KNNClassifier pkg=NearestNeighbors

knnc = KNNClassifier(K=1)
clf = machine(knnc, X3, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
@show mean(ŷ .== y[test])

knnc.K = 3
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)
@show mean(ŷ .== y[test])

caravan  = dataset("ISLR", "Caravan")
@show size(caravan)

purchase = caravan.Purchase
vals     = unique(purchase)

nl1 = sum(purchase .== vals[1])
nl2 = sum(purchase .== vals[2])
println("#$(vals[1]) ", nl1)
println("#$(vals[2]) ", nl2)

y, X = unpack(caravan, ==(:Purchase), col->true)

std = machine(Standardizer(), X)
fit!(std)
Xs = transform(std, X)

var(Xs[:,1])

test = 1:1000
train = last(test)+1:nrows(Xs);

clf = machine(KNNClassifier(K=3), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

misclassification_rate(ŷ, y[test])

mean(y[test] .!= "No")

clf = machine(LogisticClassifier(), Xs, y)
fit!(clf, rows=train)
ŷ = predict_mode(clf, rows=test)

misclassification_rate(ŷ, y[test])

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
module Benchmarks

export @sync_gpu,
       Slab, Pencil,
       print_system_info,
       run_benchmarks,
       benchmarks_dataframe,
       benchmarks_pretty_table,
       gpu_speedups_suite,
       speedups_suite,
       speedups_dataframe

using Logging
using BenchmarkTools
using DataFrames
using PrettyTables
using CUDA

using BenchmarkTools: prettytime, prettymemory
using Oceananigans: OceananigansLogger
using Oceananigans.Architectures: CPU, GPU
using Oceananigans.Utils: oceananigans_versioninfo, versioninfo_with_gpu

abstract type AbstractDomainDecomposition end

struct Slab <: AbstractDomainDecomposition end
struct Pencil <: AbstractDomainDecomposition end

function __init__()
    Logging.global_logger(OceananigansLogger())
end

macro sync_gpu(expr)
    return CUDA.has_cuda() ? :($(esc(CUDA.@sync expr))) : :($(esc(expr)))
end

function print_system_info()
    println()
    println(oceananigans_versioninfo())
    println(versioninfo_with_gpu())
    println()
    return nothing
end

function run_benchmarks(benchmark_fun; kwargs...)
    keys = [p.first for p in kwargs]
    vals = [p.second for p in kwargs]

    cases = Iterators.product(vals...)
    n_cases = length(cases)

    tags = string.(keys)
    suite = BenchmarkGroup(tags)
    for (n, case) in enumerate(cases)
        @info "Benchmarking $n/$n_cases: $case..."
        suite[case] = benchmark_fun(case...)
    end
    return suite
end

function benchmarks_dataframe(suite)
    names = Tuple(Symbol(tag) for tag in suite.tags)
    df_names = (names..., :min, :median, :mean, :max, :memory, :allocs, :samples)
    empty_cols = Tuple([] for k in df_names)
    df = DataFrame(; NamedTuple{df_names}(empty_cols)...)

    for case in keys(suite)
        trial = suite[case]
        entry = NamedTuple{names}(case) |> pairs |> Dict{Any,Any}

        entry[:min] = minimum(trial.times) |> prettytime
        entry[:median] = median(trial.times) |> prettytime
        entry[:mean] = mean(trial.times) |> prettytime
        entry[:max] = maximum(trial.times) |> prettytime
        entry[:memory] = prettymemory(trial.memory)
        entry[:allocs] = trial.allocs
        entry[:samples] = length(trial)

        push!(df, entry)
    end

    return df
end

function benchmarks_pretty_table(df; title="")
    header = propertynames(df) .|> String
    pretty_table(df, header, nosubheader=true, title=title, title_alignment=:c,
                 title_autowrap = true, title_same_width_as_table = true)

    html_filename = replace(title, ' ' => '_') * ".html"
    @info "Writing $html_filename..."
    open(html_filename, "w") do io
        html_table = pretty_table(String, df, header, nosubheader=true,
                                  title=title, title_alignment=:c,
                                  backend=:html, tf=tf_html_simple)
        write(io, html_table)
    end

    return nothing
end

is_arch_type(e) = e == CPU || e == GPU
cpu_case(case) = Tuple(is_arch_type(e) ? CPU : e for e in case)
gpu_case(case) = Tuple(is_arch_type(e) ? GPU : e for e in case)

function gpu_speedups_suite(suite)
    tags = filter(e -> !occursin("arch", lowercase(e)), suite.tags)
    suite_speedup = BenchmarkGroup(tags)

    for case in keys(suite)
        case_cpu = cpu_case(case)
        case_gpu = gpu_case(case)
        case_speedup = filter(!is_arch_type, case_cpu)

        if case_speedup ∉ keys(suite_speedup)
            suite_speedup[case_speedup] = ratio(median(suite[case_gpu]), median(suite[case_cpu]))
        end
    end

    return suite_speedup
end

function speedups_suite(suite; base_case)
    suite_speedup = BenchmarkGroup(suite.tags)

    for case in keys(suite)
        suite_speedup[case] = ratio(median(suite[case]), median(suite[base_case]))
    end

    return suite_speedup
end

function speedups_dataframe(suite; slowdown=false, efficiency=nothing, base_case=nothing, key2rank=identity)
    names = Tuple(Symbol(tag) for tag in suite.tags)
    speed_type = slowdown ? :slowdown : :speedup
    df_names = isnothing(efficiency) ? (names..., speed_type, :memory, :allocs) : (names..., speed_type, :efficiency, :memory, :allocs)
    empty_cols = Tuple([] for k in df_names)
    df = DataFrame(; NamedTuple{df_names}(empty_cols)...)

    for case in keys(suite)
        trial_ratio = suite[case]
        entry = NamedTuple{names}(case) |> pairs |> Dict{Any,Any}

        entry[speed_type] = slowdown ? trial_ratio.time :  1/trial_ratio.time
        entry[:memory] = trial_ratio.memory
        entry[:allocs] = trial_ratio.allocs

        if efficiency == :strong
            R = key2rank(case)
            t₁ = suite[base_case].time
            tₙ = suite[case].time
            entry[:efficiency] = t₁ / (R * tₙ)
        end

        if efficiency == :weak
            R = key2rank(case)
            t₁ = suite[base_case].time
            tₙ = suite[case].time
            entry[:efficiency] = t₁ / tₙ
        end

        push!(df, entry)
    end

    return df
end

end # module


text:
# example showing how to do maximum likelihood estimation
# for data generated by CLRM with normality
# sample size
using Econometrics, Distributions
function main()
n = 1000
# random true parameters
theta = [1.0,-1.0,0.0,1.0]
b = theta[1:3]
sig = theta[4]
# generate random data
x = [ones(n,1) rand(n,2)]
e = sig*randn(n,1)
y = x*b + e
# do ML
obj = theta -> -mean(normal(theta, y, x))
thetahat, logL, junk = fminunc(obj, theta)
logL = -logL
# restricted ML
R = [0 1 0 0]
r = -1.5
thetahatr, logLr, junk = fmincon(obj, theta, R, r)
logLr = -logLr
# LR rest
LR = 2.0*n*(logL-logLr)
# results
println("the true parameters: ", theta)
println("the ML estimates: ", thetahat)
println("the restricted ML estimates: ", thetahatr)
println("the LR test: ", LR);
println("p-value of the LR test: ", 1.0-cdf(Chisq(1),LR))
end
main()

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
using Flux
using OceanParameterizations
using WindMixing
using OrdinaryDiffEq, DiffEqSensitivity
using JLD2
using FileIO

train_files = ["-1e-3", "-8e-4"]

PATH = pwd()
OUTPUT_PATH = joinpath(PATH, "training_output")

# FILE_PATH = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "testNN.jld2")

FILE_PATH_uw = joinpath(OUTPUT_PATH, "uw_NN_training_2sim_-8e-4_-1e-3_large.jld2")
FILE_PATH_vw = joinpath(OUTPUT_PATH, "vw_NN_training_2sim_-8e-4_-1e-3_large.jld2")
FILE_PATH_wT = joinpath(OUTPUT_PATH, "wT_NN_training_2sim_-8e-4_-1e-3_large.jld2")

FILE_PATH_NDE = joinpath(OUTPUT_PATH, "NDE_training_convective_adjustment_2sim_-8e-4_-1e-3_large.jld2")

# FILE_PATH_uw = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "uw_NN_training_2sim_-2.5e-4_-7.5e-4_large.jld2")
# FILE_PATH_vw = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "vw_NN_training_2sim_-2.5e-4_-7.5e-4_large.jld2")
# FILE_PATH_wT = joinpath("D:\\University Matters\\Massachusetts Institute of Technology\\CLiMA Project\\OceanParameterizations.jl\\training_output", "wT_NN_training_2sim_-2.5e-4_-7.5e-4_large.jld2")

𝒟train = data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=true)

N_inputs = 96
hidden_units = 400
N_outputs = 31
uw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))
vw_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))
wT_NN = Chain(Dense(N_inputs, hidden_units, relu), Dense(hidden_units, hidden_units, relu), Dense(hidden_units, N_outputs))

train_optimizers = [ADAM(0.01), Descent()]
train_epochs = [10,10]

write_metadata_NN_training(FILE_PATH_uw, train_files, train_epochs, train_optimizers, uw_NN, "uw")
uw_weights = train_NN(uw_NN, 𝒟train.uvT_scaled, 𝒟train.uw.scaled, train_optimizers, train_epochs, FILE_PATH_uw, "uw")

write_metadata_NN_training(FILE_PATH_vw, train_files, train_epochs, train_optimizers, vw_NN, "vw")
vw_weights = train_NN(vw_NN, 𝒟train.uvT_scaled, 𝒟train.vw.scaled, train_optimizers, train_epochs, FILE_PATH_vw, "vw")

write_metadata_NN_training(FILE_PATH_wT, train_files, train_epochs, train_optimizers, wT_NN, "wT")
wT_weights = train_NN(wT_NN, 𝒟train.uvT_scaled, 𝒟train.wT.scaled, train_optimizers, train_epochs, FILE_PATH_wT, "wT")

uw_NN = Flux.destructure(uw_NN)[2](uw_weights)
vw_NN = Flux.destructure(vw_NN)[2](vw_weights)
wT_NN = Flux.destructure(wT_NN)[2](wT_weights)

train_epochs = [1, 1, 1, 1, 1, 1, 1, 1, 10]
train_tranges = [1:5:50, 1:5:100, 1:10:200, 1:20:400, 1:20:500, 1:25:700, 1:25:900, 1:30:1000, 1:40:1153]
train_optimizers = [[ADAM(0.01)] for i in 1:length(train_epochs)]
timestepper = ROCK4()

function train(FILE_PATH, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper)
    write_metadata_NDE_training(FILE_PATH, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN)

    for i in 1:length(train_epochs)
        uw_NN, vw_NN, wT_NN = train_NDE_convective_adjustment(uw_NN, vw_NN, wT_NN, 𝒟train, train_tranges[i], timestepper, train_optimizers[i], train_epochs[i], FILE_PATH, i, 1)
    end

end

train(FILE_PATH_NDE, train_files, train_epochs, train_tranges, train_optimizers, uw_NN, vw_NN, wT_NN, 𝒟train, timestepper)


text:
module PandorasBoxPlayerClone

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/playerClone" PlayerClone(x::Integer, y::Integer, mode::String="Inventory", flag::String="")

const placements = Ahorn.PlacementDict(
    "Player Clone (Unholy, Pandora's Box)" => Ahorn.EntityPlacement(
        PlayerClone
    )
)

const modes = String[
    "Inventory",
    "Backpack",
    "NoBackpack",
    "MadelineAsBadeline"
]

Ahorn.editingOptions(entity::PlayerClone) = Dict{String, Any}(
    "mode" => modes
)

sprite = "characters/player/sitDown00.png"

function Ahorn.selection(entity::PlayerClone)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PlayerClone) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
function numdivisors(n)
	d = 0
	limit = n;

    i = 1
	while i < limit
		if n % i == 0
			limit = div(n,i) # converge to sqrt(n)
			if limit != i
				d += 1
			end
			d += 1
		end
		i += 1
	end
	return d
end

function main()
	N = 5

	n = 0
	i = 0
	while true
		i += 1
		n += i
		d = numdivisors(n)
		println("$d $n")
		if d > N
			println(n)
			return
		end
	end
end

main()

max_stars_repo_name: cmey/euler

max_issues_repo_name: cmey/euler

max_forks_repo_name: cmey/euler


text:
#clustered nearest neaighbor
function Kclearest_neighbor{T<:Real}(distmat::Matrix{T},n::Int,
                                    firstcity::Int = 0)
  scan = kmeans(distmat,n)
  asgn = scan.assignments
  bitClust = trues(n)
	numCities = size(distmat, 1)
  bitNode = trues(numCities)
  master = collect(1:numCities)
	# put first city on path
	path = Vector{Int}()
  sizehint!(path,numCities)
	if firstcity == 0
		firstcity = rand(1:numCities)
	end
	push!(path, firstcity)
	# cities to visit
  while in(true,bitNode)
    curCity = path[end]
    bitNode[curCity] = false
    Cid = asgn[curCity]
    if Cid == 0
      citiesToVisit = master[bitNode]
    elseif bitClust[Cid]
      bitPals = asgn.==Cid
      citiesToVisit = master[bitNode&bitPals]
      bitClust[Cid] = false
    else
      citiesToVisit = master[bitNode]
    end
	# nearest neighbor loop
  	while asgn[curCity] == Cid && !isempty(citiesToVisit)
      dists = distmat[citiesToVisit,curCity]
  		_, nextInd = findmin(dists)
  		nextCity = citiesToVisit[nextInd]
  		push!(path, nextCity)
      bitNode[nextCity] = false
  		deleteat!(citiesToVisit, nextInd)
      curCity = path[end]
  	end
  end

	return path
end


text:
# BraidingTensor:
# special (2,2) tensor that implements a standard braiding operation
#====================================================================#
"""
    struct BraidingTensor{S<:IndexSpace} <: AbstractTensorMap{S, 2, 2}

Specific subtype of [`AbstractTensorMap`](@ref) for representing the braiding tensor that
braids the first input over the second input; its inverse can be obtained as the adjoint.
"""
struct BraidingTensor{S<:IndexSpace, A} <: AbstractTensorMap{S, 2, 2}
    V1::S
    V2::S
    adjoint::Bool
    function BraidingTensor(V1::S, V2::S, adjoint::Bool = false, ::Type{A} = Matrix{ComplexF64}) where
                {S<:IndexSpace, A<:DenseMatrix}
        for a in sectors(V1)
            for b in sectors(V2)
                for c in (a ⊗ b)
                    Nsymbol(a, b, c) == Nsymbol(b, a, c) ||
                        throw(ArgumentError("Cannot define a braiding between $a and $b"))
                end
            end
        end
        return new{S,A}(V1, V2, adjoint)
        # partial construction: only construct rowr and colr when needed
    end
end

Base.adjoint(b::BraidingTensor{S,A}) where {S<:IndexSpace, A<:DenseMatrix} =
    BraidingTensor(b.V1, b.V2, !b.adjoint, A)

domain(b::BraidingTensor) = b.adjoint ? b.V2 ⊗ b.V1 : b.V1 ⊗ b.V2
codomain(b::BraidingTensor) = b.adjoint ? b.V1 ⊗ b.V2 : b.V2 ⊗ b.V1

storagetype(::Type{BraidingTensor{S,A}}) where {S<:IndexSpace, A<:DenseMatrix} = A

blocksectors(b::BraidingTensor) = blocksectors(b.V1 ⊗ b.V2)
hasblock(b::BraidingTensor, s::Sector) = s ∈ blocksectors(b)

function fusiontrees(b::BraidingTensor)
    codom = codomain(b)
    dom = domain(b)
    I = sectortype(b)
    F = fusiontreetype(I, 2)
    rowr = SectorDict{I, FusionTreeDict{F, UnitRange{Int}}}()
    colr = SectorDict{I, FusionTreeDict{F, UnitRange{Int}}}()
    for c in blocksectors(codom)
        rowrc = FusionTreeDict{F, UnitRange{Int}}()
        colrc = FusionTreeDict{F, UnitRange{Int}}()
        offset1 = 0
        for s1 in sectors(codom)
            for f1 in fusiontrees(s1, c, map(isdual, codom.spaces))
                r = (offset1 + 1):(offset1 + dim(codom, s1))
                push!(rowrc, f1 => r)
                offset1 = last(r)
            end
        end
        dim1 = offset1
        offset2 = 0
        for s2 in sectors(dom)
            for f2 in fusiontrees(s2, c, map(isdual, dom.spaces))
                r = (offset2 + 1):(offset2 + dim(dom, s2))
                push!(colrc, f2 => r)
                offset2 = last(r)
            end
        end
        dim2 = offset2
        push!(data, c=>f((dim1, dim2)))
        push!(rowr, c=>rowrc)
        push!(colr, c=>colrc)
    end
    return TensorKeyIterator(rowr, colr)
end

@inline function Base.getindex(b::BraidingTensor, f1::FusionTree{I,2}, f2::FusionTree{I,2}) where {I<:Sector}
    I == sectortype(b) || throw(SectorMismatch())
    c = f1.coupled
    V1, V2 = domain(b)
    @boundscheck begin
        c == f2.coupled || throw(SectorMismatch())
        ((f1.uncoupled[1] ∈ V2) && (f2.uncoupled[1] ∈ V1)) || throw(SectorMismatch())
        ((f1.uncoupled[2] ∈ V1) && (f2.uncoupled[2] ∈ V2)) || throw(SectorMismatch())
    end
    @inbounds begin
        d = (dims(V2 ⊗ V1, f1.uncoupled)..., dims(V1 ⊗ V2, f2.uncoupled)...)
        n1 = d[1]*d[2]
        n2 = d[3]*d[4]
        a1, a2 = f2.uncoupled
        data = fill!(storagetype(b)(undef, (n1, n2)), zero(eltype(b)))
        if f1.uncoupled == (a2, a1)
            braiddict = artin_braid(f2, 1; inv = b.adjoint)
            r = get(dict, f1, zero(valtype(braiddict)))
            data[1:(n1+1):end] .= r # set diagonal
        end
        return permutedims(sreshape(StridedView(data), d), (1,2,4,3))
    end
end

Base.copy(b::BraidingTensor) = copy!(similar(b), b)
function Base.copy!(t::TensorMap, b::BraidingTensor)
    space(t) == space(b) || throw(SectorMismatch())
    fill!(t, zero(eltype(t)))
    for (f1, f2) in fusiontrees(t)
        a1, a2 = f2.uncoupled
        c = f2.coupled
        f1.uncoupled == (a2, a1) || continue
        data = t[f1, f2]
        r = artin_braid(f2, 1; inv = b.adjoint)[f1]
        for i = 1:size(data, 1), j = 1:size(data, 2)
            data[i, j, j, i] = r
        end
    end
    return t
end
TensorMap(b::BraidingTensor) = copy(b)
Base.convert(::Type{TensorMap}, b::BraidingTensor) = copy(b)

function block(b::BraidingTensor, s::Sector)
    sectortype(b) == typeof(s) || throw(SectorMismatch())
    return block(copy(b), s) # stopgap awaiting actual implementation

    # V1, V2 = codomain(b)
    # n = blockdim(V1 ⊗ V2, s)
    # data = fill!(storagetype(b)(undef, (n, n)), eltype(b))
    # n == 0 && return data
    # for (f1, f2) in fusiontrees(b)
    #     f2.coupled == s || continue
    #     a1, a2 = f2.uncoupled
    #
end


text:
# build a single GAP binary which supports Julia 1.3 - 1.5
include("../common.jl")

platforms = configure(v"1.5", v"1.5.3")

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", compat="6.1.2"),
    Dependency("Readline_jll", compat="8.0.4"),
    Dependency("Zlib_jll"),
    BuildDependency(PackageSpec(name="libjulia_jll", version=v"1.5.3")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"7", init_block=init_block)

max_stars_repo_name: singularitti/Yggdrasil

max_issues_repo_name: singularitti/Yggdrasil

max_forks_repo_name: singularitti/Yggdrasil


text:
module MNIST

    using Pkg
    using MLDatasets
    using JLD2
    using ..Misc

    const FILEPATH = joinpath(Misc.DATASETS_FOLDER, "bitsmnist.jld2")

    function avgwhite(x)
        s = sum(x)
        cnt = count(x -> x > 0, x)
        s / cnt
    end

    function convert2binary(x)
        avg = avgwhite(x)
        rx = reshape(x, Misc.IMG_LENGTH, size(x, 3))
        [BitVector(rx[:, i] .> avg) for i in axes(x, 3)]
    end

    """
        download()

    Forces downloading the bits mnist dataset
    """
    function download()
        MLDatasets.MNIST.download(i_accept_the_terms_of_use=true)
        train_x, train_y = MLDatasets.MNIST.traindata()
        test_x, test_y = MLDatasets.MNIST.testdata()

        btrain_x, btest_x = convert2binary(train_x), convert2binary(test_x)
        set = Dict("train_x" => btrain_x,
                                    "train_y" => train_y,
                                    "test_x" => btest_x,
                                    "test_y" => test_y)
        save(FILEPATH, set)
        set
    end

    """
        mnist()

    Downloads and return the bits mnist dataset
    """
    function mnist()
        if isfile(FILEPATH)
            return load(FILEPATH)
        else
            return download()
        end
    end

    export mnist

end


text:
@with_kw struct SpikingSynapseParameter
    τpre::SNNFloat = 20ms
    τpost::SNNFloat = 20ms
    Wmax::SNNFloat = 0.01
    ΔApre::SNNFloat = 0.01 * Wmax
    ΔApost::SNNFloat = -ΔApre * τpre / τpost * 1.05
end

@with_kw mutable struct SpikingSynapse
    param::SpikingSynapseParameter = SpikingSynapseParameter()
    rowptr::Vector{SNNInt} # row pointer of sparse W
    colptr::Vector{SNNInt} # column pointer of sparse W
    I::Vector{SNNInt}      # postsynaptic index of W
    J::Vector{SNNInt}      # presynaptic index of W
    index::Vector{SNNInt}  # index mapping: W[index[i]] = Wt[i], Wt = sparse(dense(W)')
    W::Vector{SNNFloat}  # synaptic weight
    tpre::Vector{SNNFloat} = simzeros(W) # presynaptic spiking time
    tpost::Vector{SNNFloat} = simzeros(W) # postsynaptic spiking time
    Apre::Vector{SNNFloat} = simzeros(W) # presynaptic trace
    Apost::Vector{SNNFloat} = simzeros(W) # postsynaptic trace
    fireI::Vector{Bool} # postsynaptic firing
    fireJ::Vector{Bool} # presynaptic firing
    g::Vector{SNNFloat} # postsynaptic conductance
    records::Dict = Dict()
end

function SpikingSynapse(pre, post, sym; σ = 0.0, p = 0.0)
    w = σ * sprand(post.N, pre.N, p)
    rowptr, colptr, I, J, index, W = dsparse(w)
    fireI, fireJ = post.fire, pre.fire
    g = getfield(post, sym)
    SpikingSynapse(;@symdict(rowptr, colptr, I, J, index, W, fireI, fireJ, g)...)
end

function forward!(c::SpikingSynapse, param::SpikingSynapseParameter)
    @unpack colptr, I, W, fireJ, g = c
    @inbounds for j in 1:(length(colptr) - 1)
        if fireJ[j]
            for s in colptr[j]:(colptr[j+1] - 1)
                g[I[s]] += W[s]
            end
        end
    end
end

function plasticity!(c::SpikingSynapse, param::SpikingSynapseParameter, dt::SNNFloat, t::SNNFloat)
    @unpack rowptr, colptr, I, J, index, W, tpre, tpost, Apre, Apost, fireI, fireJ, g = c
    @unpack τpre, τpost, Wmax, ΔApre, ΔApost = param
    @inbounds for j in 1:(length(colptr) - 1)
        if fireJ[j]
            for s in colptr[j]:(colptr[j+1] - 1)
                Apre[s] *= exp32(- (t - tpre[s]) / τpre)
                Apost[s] *= exp32(- (t - tpost[s]) / τpost)
                Apre[s] += ΔApre
                tpre[s] = t
                W[s] = clamp(W[s] + Apost[s], 0f0, Wmax)
            end
        end
    end
    @inbounds for i in 1:(length(rowptr) - 1)
        if fireI[i]
            for st in rowptr[i]:(rowptr[i+1] - 1)
                s = index[st]
                Apre[s] *= exp32(- (t - tpre[s]) / τpre)
                Apost[s] *= exp32(- (t - tpost[s]) / τpost)
                Apost[s] += ΔApost
                tpost[s] = t
                W[s] = clamp(W[s] + Apre[s], 0f0, Wmax)
            end
        end
    end
end

max_stars_repo_name: russelljjarvis/juliaSNN

max_issues_repo_name: russelljjarvis/juliaSNN

max_forks_repo_name: russelljjarvis/juliaSNN


text:
function hello_world_and_no_arguments(;) # Attention! The function must accept keyword arguments, therefore semicolon. Here just a zero number of kwargs
    println("Hello World!")
    return (;) # Attention! If no values are to be returned, then return an empty named tuple.
end

fns = (;hello_world_and_no_arguments, )


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the Greenwich Mean Sideral Time (GMST).
#
# Remarks
# ==============================================================================
#
#   This files was copied from SatelliteToolbox.jl because the functionality is
#   required by SGP4 model. The functions here **must not** be exported to avoid
#   interferences.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    j2000_to_gmst(J2000_UT1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] given the instant
`J2000_UT1` in J2000.0 reference [UT1].

# Remarks

Based on algorithm in [2] (http://www.navipedia.net/index.php/CEP_to_ITRF),
accessed at 2015-12-01.

"""
function j2000_to_gmst(J2000_UT1::Number)
	# Julian centuries elapsed from the epoch J2000.0.
	T_UT1 = J2000_UT1 / 36525

	# Greenwich Mean Sideral Time at T_UT1 [s].
    θ_GMST = @evalpoly(
        T_UT1,
        + 67310.54841,
        + (876600.0*3600 + 8640184.812866),
        + 0.093104,
        - 6.2e-6
    )

    # Reduce to the interval [0, 86400]s.
    θ_GMST = mod(θ_GMST, 86400)

    # Convert to radian and return.
    return θ_GMST * π / 43200
end

"""
    jd_to_gmst(JD_UT1::Number)

Compute the Greenwich Mean Sideral Time (GMST) \\[rad] for the Julian Day
`JD_UT1` [UT1].

# Remarks

Based on algorithm in [1, pp. 188].

"""
jd_to_gmst(JD_UT1::Number) = j2000_to_gmst(JD_UT1 - JD_J2000)