text:
using Interact, Distributed
import Plots, Printf, SearchBall
const pl = Plots
const pr = Printf
const sb = SearchBall

pl.gr()
addprocs(8)
@everywhere using SharedArrays

"""
robot lies in (0,0,0)
0 Rad -> orientation in positive x direction
"""
const v = 200.0
const ω = π/16.0
@everywhere function duration(x::Number, y::Number, start::Tuple{<:Number, <:Number, <:Number}=(0.0, 0.0, 0.0); velocity::Number=v, angular_velocity::Number=ω)
    (xi, yi, ai) = start
    a = mod(atan(y-yi, x-xi) - ai + π, 2π) - π
    sqrt((x-xi)^2 + (y-yi)^2)/velocity + abs(a)/angular_velocity
end

@everywhere const xs = -200:200
@everywhere const ys = -200:200
const maxxs = maximum(xs)
const maxys = maximum(ys)
const hmap1 = SharedArray(zeros(Float64, length(xs), length(ys)))
const hmap2 = SharedArray(zeros(Float64, length(xs), length(ys)))
const hmap3 = SharedArray(zeros(Float64, length(xs), length(ys)))

const init1 = (-100, 0, -π/2.0)
const init2 = (100, 0, π/2.0)

function calc_area!(hmap1, hmap2, hmap3, init1, init2, v=v, ω=ω)
    @inbounds @sync @distributed for ix in 1:length(xs)
        for iy in 1:length(ys)
            x = xs[ix]
            y = ys[iy]
            hmap1[iy, ix] = duration(x, y, init1, velocity=v, angular_velocity=ω)
            hmap2[iy, ix] = duration(x, y, init2, velocity=v, angular_velocity=ω)
            hmap3[iy, ix] = hmap1[iy, ix] < hmap2[iy, ix] ? 0.0 : 1.0
        end
    end
    # hmap1[3.5.<hmap1.<3.6] .= maximum(hmap1)
    # hmap2[3.5.<hmap2.<3.6] .= maximum(hmap2)
    return
end

pl.heatmap(hmap1, color=:ice_r, aspect_ratio=:equal, title=pr.@sprintf("v=%.1fm/s ω=%.0f°/s", v/100.0, rad2deg(ω)))
pl.heatmap(hmap2, color=:ice_r, aspect_ratio=:equal, title=pr.@sprintf("v=%.1fm/s ω=%.0f°/s", v/100.0, rad2deg(ω)))
pl.heatmap(hmap3, color=:ice, aspect_ratio=:equal, title=pr.@sprintf("v=%.1fm/s ω=%.0f°/s", v/100.0, rad2deg(ω)))

const mx = [maxxs, maxxs]
const my = [maxys, maxys]
const i11 = [init1[1], init1[1]]
const i12 = [init1[2], init1[2]]
const i21 = [init2[1], init2[1]]
const i22 = [init2[2], init2[2]]
ui = @manipulate for angle_red in -π:π/32:π, angle_blue in -π:π/32:π, vel in 1:400, ang in π/32:π/32:2π, show_field in [:voronoi, :red_player, :blue_player]
    calc_area!(hmap1, hmap2, hmap3, (init1[1], init1[2], angle_red), (init2[1], init2[2], angle_blue), vel, ang)
    if show_field == :voronoi
        myplot = pl.heatmap(hmap3, color=:ice, aspect_ratio=:equal, title=pr.@sprintf("v=%.1fm/s ω=%.0f°/s", vel/100.0, rad2deg(ang)), colorbar=false)
        # myplot = pl.plot(colorview(Gray, hmap3), aspect_ratio=:equal, title=pr.@sprintf("v=%.1fm/s ω=%.0f°/s", vel/100.0, rad2deg(ang)))
        dir1 = sb.rotate2d([30.0, 0.0], angle_red)
        pl.plot!(myplot, [0.0, dir1[1]] .+ mx .+ i11, [0.0, dir1[2]] .+ my .+ i12, lw=3, legend=false, color=:red)
        # pl.scatter!(myplot, init1[1], init1[2], lw=8, col=:white)
        dir2 = sb.rotate2d([30.0, 0.0], angle_blue)
        pl.plot!(myplot, [0.0, dir2[1]] .+ mx .+ i21, [0.0, dir2[2]] .+ my .+ i22, lw=3, legend=false, color=:blue)
        # pl.scatter!(myplot, init2[1], init2[2], lw=8, col=:black)
    elseif show_field == :red_player
        pl.heatmap(hmap1, color=:ice, aspect_ratio=:equal)
    elseif show_field == :blue_player
        pl.heatmap(hmap2, color=:ice, aspect_ratio=:equal)
    end
end


text:
module RocketSkipErrorOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: skip_error()" begin

    println("Testing: operator skip_error()")

    run_proxyshowcheck("SkipError", skip_error())

    run_testset([
        (
            source = from(1:5) |> skip_error(),
            values = @ts([ 1:5, c ])
        ),
        (
            source = completed() |> skip_error(),
            values = @ts(c)
        ),
        (
            source = faulted(1) |> skip_error(),
            values = @ts()
        ),
        (
            source = never() |> skip_error(),
            values = @ts()
        )
    ])

end

end

ext: jl


text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold:
    Household,
    get_benefit_units,
    get_head

using .STBIncomes

using .Results:
    BenefitUnitResult,
    HouseholdResult,
    HousingResult,
    IndividualResult,
    ITResult,
    LMTCanApplyFor,
    LMTIncomes,
    LMTResults,
    LocalTaxes,
    NIResult,

    aggregate_tax,
    aggregate!,
    has_any,
    init_benefit_unit_result,
    init_household_result,
    map_incomes,
    to_string,
    total
using .ExampleHelpers
#
#
#

@testset "Results record creation and summation" begin
    sph = get_example( single_parent_hh )
    sph.gross_rent = 200.0
    sph.water_and_sewerage = 0
    sph.mortgage_payment = 0

    hrs :: HouseholdResult = init_household_result( sph )
    head :: Person = get_head(sph)
    hp = head.pid
    headr :: IndividualResult = hrs.bus[1].pers[hp]
    headr.income[INCOME_TAX] = 100.0
    aggregate!( sph, hrs )
    @test hrs.income[INCOME_TAX] == 100.0
    @test hrs.bus[1].income[INCOME_TAX] == 100.0
    # check we don't ever double-count
    aggregate!( sph, hrs )
    @test hrs.income[INCOME_TAX] == 100.0
    @test hrs.bus[1].income[INCOME_TAX] == 100.0
    headr.income[WAGES] = 1_000.0
    aggregate!( sph, hrs )
    @test hrs.income[WAGES] == 1_000.0
    @test hrs.bus[1].income[WAGES] == 1_000.0
    @test has_any( hrs, WAGES )
    @test ! has_any( hrs, SELF_EMPLOYMENT_INCOME )

    @test total( hrs, WAGES ) == 1_000.0
    @test total( hrs, SELF_EMPLOYMENT_INCOME ) == 0.0
    println( inctostr( hrs.income ))

    others = 0.2868852459+3.1912568306+30.0+9.95

    @test hrs.bhc_net_income ≈ 900.0+others
    @test hrs.ahc_net_income ≈ 700.0+others

    # multiple benefit units example
    mbhh = get_example( mbu )
    bus = get_benefit_units(mbhh)
    mhead =  get_head(bus[3])
    mhp = mhead.pid

    mhead2 =  get_head(bus[2])
    mhp2 = mhead2.pid

    mhrs  = init_household_result( mbhh )
    mheadr = mhrs.bus[3].pers[mhp]
    mheadr.income[WAGES] = 500.0

    mheadr.ni.class_1_primary = 123.12
    mheadr.it.non_savings_tax = 321.21

    mheadr2 = mhrs.bus[2].pers[mhp2]
    mheadr2.ni.class_1_primary = 100
    mheadr2.it.non_savings_tax = 200

    t_wages = total( mhrs, WAGES )
    mheadr.income[LOCAL_TAXES] = 100
    aggregate!( mbhh, mhrs )

    @test mhrs.bhc_net_income ≈ t_wages
    @test mhrs.ahc_net_income ≈ t_wages - 100 - mbhh.gross_rent

    @test mhrs.income[LOCAL_TAXES] == 100
    @test mhrs.bus[1].income[LOCAL_TAXES] == 0
    @test mhrs.bus[3].income[LOCAL_TAXES] == 100

    @test total( mhrs, LOCAL_TAXES ) == 100
    @test has_any( mhrs, LOCAL_TAXES )
    @test has_any( mhrs, LOCAL_TAXES, CHILD_BENEFIT )
    @test has_any( mhrs, [LOCAL_TAXES, CHILD_BENEFIT] )
    @test ! has_any( mhrs, CHILD_BENEFIT,INCOME_SUPPORT )
    @test ! has_any( mhrs, [CHILD_BENEFIT,INCOME_SUPPORT] )

    it, ni = aggregate_tax( mhrs.bus[2] )
    @test ni.class_1_primary == 100
    @test it.non_savings_tax == 200
    @test it.taxable_income == 0

end


text:
using Test, YaoBlocks, LuxurySparse, BitBasis, YaoBlocks.ConstGate

@testset "test builtin gates" begin

    @testset "test $each" for each in [X, Y, Z, H, CNOT, SWAP, Toffoli]
        @test isunitary(each) == true
        @test isreflexive(each) == true
        @test ishermitian(each) == true
    end

    @testset "test $each nqubits" for each in [X, Y, Z, H, T, S, Tdag, Sdag, P0, P1, Pu, Pd]
        @test nqubits(each) == 1
    end

    @testset "test $each" for each in [S, T]
        @test isunitary(each)
        @test isreflexive(each) == false
        @test ishermitian(each) == false
    end

    @testset "test $each" for each in [P0, P1]
        @test isunitary(each) == false
        @test isreflexive(each) == false
        @test ishermitian(each) == true
    end

    @testset "test $each" for each in [Pu, Pd]
        @test isunitary(each) == false
        @test isreflexive(each) == false
        @test ishermitian(each) == false
    end

    @test nqubits(CNOT) == 2
    @test nqubits(CZ) == 2
    @test nqubits(SWAP) == 2
    @test nqubits(Toffoli) == 3

end

@testset "matrix" begin
    CNOT_R = PermMatrix([1, 2, 4, 3], ones(ComplexF64, 4))
    Toffoli_R = PermMatrix([1, 2, 3, 4, 5, 6, 8, 7], ones(ComplexF64, 8))

    for (each, MAT) in [
        (X, [0 1; 1 0]),
        (Y, [0 -im; im 0]),
        (Z, [1 0; 0 -1]),
        (H, (elem = 1 / sqrt(2); [elem elem; elem -elem])),
    ]
        @test mat(each) ≈ MAT

    end
    @test mat(CNOT) |> invorder == CNOT_R
    @test mat(Toffoli) |> invorder == Toffoli_R
    @test mat(T) * mat(T) ≈ mat(S)

    @test mat(T)' ≈ mat(T')
    @test mat(Tdag)' ≈ mat(Tdag')
    @test T' isa TdagGate
    @test Tdag' isa TGate

    @test mat(S)' ≈ mat(S')
    @test mat(Sdag)' ≈ mat(Sdag')
    @test S' isa SdagGate
    @test Sdag' isa SGate
end

@testset "test @const_gate" begin

    @testset "bind new type" begin
        @test @allocated(mat(X)) == 0
        @test @allocated(mat(ComplexF32, X)) > 0
        @const_gate X::ComplexF32
        @test @allocated(mat(ComplexF32, XGate())) == 0
    end

    @testset "define new" begin
        @const_gate TEST = rand(ComplexF64, 2, 2)

        # errors if given matrix is not a square matrix
        @test_throws DimensionMismatch @const_gate TEST::ComplexF32 = rand(2, 3)

        @const_gate TEST::ComplexF32
        @test @allocated(mat(ComplexF32, TEST)) == 0
    end

end

@testset "I gate" begin
    g = ConstGate.IGate{2}()
    @test mat(g) ≈ IMatrix{4,ComplexF64}()
    @test ishermitian(g)
    @test isunitary(g)
end

@testset "test adjoints" begin
    adjoint(Pu) == Pd
    adjoint(Pd) == Pu
end


text:
##############
photodatafolder = joinpath(splitdir(@__FILE__())[1],"../../data")
photodatatestfile = joinpath(photodatafolder,"naturallight.jld")
println("Loading photon sequence from file")
gc() ; photons = photonsequence(photodatatestfile)
for i=1:3
    gc() ; @time photons = photonsequence(photodatatestfile)
end
println("Loading light-induced current from file")
gc() ; lic = lightinducedcurrent(photodatatestfile)
for i=1:3
    gc() ; @time lic = lightinducedcurrent(photodatatestfile)
end
###############
testunique() = (for j=1:2000 unique(sort!(rand!(Array{Int}(150),1:30000))) end ; nothing)
testuniquesorted() = (for i=1:2000 GMHPhotoReceptor.uniquesorted(sort!(rand!(Array{Int}(150),1:30000))) end ; nothing)
testunique()
testuniquesorted()
gc()
println("Testing unique")
for i=1:3
  gc()
  @time testunique()
end
println("Testing uniquesorted")
for i=1:3
  gc()
  @time testuniquesorted()
end
###############
gc()
nvilli = 30000
nsteps = numvalues(photons)
receptorimpl = [:full,:row,:cell]
receptors = Dict{Symbol,Any}()

for args in [(:full,(UInt8,)),(:row,(UInt16,Int)),(:cell,(UInt16,))]
  receptors[args[1]] = photoreceptor(args[1],datavalues(photons),nvilli,args[2]...)
  println("Creating AbstractPhotoReceptor of type ",args[1])
  for j=1:3
    receptors[args[1]] = []
    gc()
    @time receptors[args[1]] = photoreceptor(args[1],datavalues(photons),nvilli,args[2]...)
  end
end

@inline function testmicrovilliaccess(receptor::GMHPhotoReceptor.AbstractPhotoReceptor)
    a = 0
    for t=1:numtimesteps(receptor)
        microvilli = GMHPhotoReceptor.microvilliwithphotons(receptor,t)
        for i=1:length(microvilli)
            a += microvilli[i]
        end
    end
    a
end

for i in [:full,:row,:cell]
  testmicrovilliaccess(receptors[i])
  println("Testing microvilliwithphotons access of type ",i)
  for j=1:3
    gc() ; @time testmicrovilliaccess(receptors[i])
  end
end

gc()
for i in [:full,:row,:cell]
  println("Memory taken by Microvilli object of type ",i,": ",GeneralizedMetropolisHastings.objectsizetostr(receptors[i]))
end

gc()
for T in [Float32,Float64]
    println("=================")
    println("Testing ",T)
    println("=================")
    ###Testing of bump functions
    timepoints = collect(T(1.0):T(1.0):T(30.0))
    shape = T[4.0,3.0,2.5]
    cutoff = T(0.025)
    bumpvals = zeros(T,30)
    bump!(bumpvals,timepoints,shape[1],shape[2],shape[3])
    bump(timepoints,shape[1],shape[2],shape[3])
    bump(shape[1],shape[2],shape[3],cutoff)
    println("Testing bump!{",T,"}")
    for j=1:3
        @time bumpvals = bump!(bumpvals,timepoints,shape[1],shape[2],shape[3])
    end
    println("Testing bump{",T,"} with given length")
    for j=1:3
        @time bumpvals = bump(timepoints,shape[1],shape[2],shape[3])
    end
    println("Testing bump{",T,"} with tolerance")
    for j=1:3
        @time bumpvals = bump(shape[1],shape[2],shape[3],cutoff)
    end
    bumplength = length(bumpvals)
    ###Testing of macrocurrent functions
    current = zeros(T,nsteps)
    macrocurrent!(current,datavalues(photons),bumpvals)
    macrocurrent(datavalues(photons),bumpvals)
    gc()
    println("Testing macrocurrent!{",T,"}")
    for j=1:3
        current = zeros(T,nsteps)
        @time macrocurrent!(current,datavalues(photons),bumpvals)
    end
    gc()
    println("Testing macrocurrent{",T,"}")
    for j=1:3
        @time current = macrocurrent(datavalues(photons),bumpvals)
    end

    latparas = T[2.7137,0.4568]
    refparas = T[log(150),log(1.7)]
    latency = zeros(T,sum(datavalues(photons)))
    refractory = zeros(T,sum(datavalues(photons)))
    GMHPhotoReceptor.randlognormal!(latency,latparas[1],latparas[2])
    GMHPhotoReceptor.randlognormal!(refractory,refparas[1],refparas[2])
    println("Generating latency and refractory values")
    for j=1:3
        gc() ; @time (GMHPhotoReceptor.randlognormal!(latency,latparas[1],latparas[2]) ; GMHPhotoReceptor.randlognormal!(refractory,refparas[1],refparas[2]))
    end

    for impl in [:cell]
        GMHPhotoReceptor._filterphotons!(zeros(Int,nsteps),zeros(Int,nvilli),receptors[impl],latency,refractory,bumplength)
        println("Testing _filterphotons!{$(T)} for $(impl)")
        for i=1:5
            gc()
            fphotons = zeros(Int,nsteps)
            rtimes = zeros(Int,nvilli)
            @time GMHPhotoReceptor._filterphotons!(fphotons,rtimes,receptors[impl],latency,refractory,bumplength)
        end
        filterphotons!(zeros(Int,nsteps),receptors[impl],latency,refractory,bumplength)
        println("Testing filterphotons!{$(T)} for $(impl)")
        for i=1:5
            gc()
            fphotons = zeros(Int,nsteps)
            @time filterphotons!(fphotons,receptors[impl],latency,refractory,bumplength)
        end
        filterphotons(receptors[impl],latency,refractory,bumplength)
        println("Testing filterphotons{$(T)} for $(impl)")
        for i=1:5
            gc() ; @time filterphotons(receptors[impl],latency,refractory,bumplength)
        end
        ######################################################################################
        GMHPhotoReceptor._lightinducedcurrent!(zeros(T,nsteps),zeros(Int,nsteps),zeros(Int,nvilli),receptors[impl],latency,refractory,bumpvals)
        println("Testing _lightinducedcurrent!{$(T)} for $(impl)")
        for i=1:5
            gc()
            current = zeros(T,nsteps)
            fphotons = zeros(Int,nsteps)
            rtimes = zeros(Int,nvilli)
            @time GMHPhotoReceptor._lightinducedcurrent!(current,fphotons,rtimes,receptors[impl],latency,refractory,bumpvals)
        end
        lightinducedcurrent!(zeros(T,nsteps),receptors[impl],latency,refractory,bumpvals)
        println("Testing lightinducedcurrent!{$(T)} for $(impl)")
        for i=1:5
            gc()
            current = zeros(T,nsteps)
            @time lightinducedcurrent!(current,receptors[impl],latency,refractory,bumpvals)
        end
        lightinducedcurrent(receptors[impl],latency,refractory,bumpvals)
        println("Testing lightinducedcurrrent{$(T)} for $(impl)")
        for i=1:5
            gc() ; @time lightinducedcurrent(receptors[impl],latency,refractory,bumpvals)
        end
    end

    licit = 10
    licvec = Array(T,numvalues(photons),licit)

    println("$licit iterations of lightinducedcurrent for :cell")
    for i=1:licit
        gc() ; @time (GMHPhotoReceptor.randlognormal!(latency,latparas[1],latparas[2]) ; GMHPhotoReceptor.randlognormal!(refractory,refparas[1],refparas[2]) ; licvec[:,i] = lightinducedcurrent(receptors[:cell],latency,refractory,bumpvals))
    end

    println("100 iterations of lightinducedcurrent! for :cell")
    gc()
    @time for i=1:100
        GMHPhotoReceptor.randlognormal!(latency,latparas[1],latparas[2])
        GMHPhotoReceptor.randlognormal!(refractory,refparas[1],refparas[2])
        lightinducedcurrent!(current,receptors[:cell],latency,refractory,bumpvals)
    end
end

current = data(:array,0.001f0(1:nsteps),zeros(Float32,nsteps))
p1 = policy(:photoreceptor)
p2 = policy(:photoreceptor,bump =:sample)

params1 = parameters(:photoreceptor,p1,(2.0,4.0),(0.2,0.8),(4.0,6.0),(1.5,2.5))
params2 = parameters(:photoreceptor,p2,(2.0,4.0),(0.2,0.8),(4.0,6.0),(1.5,2.5),(3.0,5.0),(log(3.0),0.3),(log(2.5),0.3))

m1 = model(:photoreceptor,params1,photons,current,[1.0],30000,p1)
m2 = model(:photoreceptor,params2,photons,current,[1.0],30000,p2)

println("Evaluating the model with latency and refractory parameters")
evaluate!(m1,[2.7,0.5,5.0,2.0])
for i=1:5
    gc() ; @time evaluate!(m1,[2.7,0.5,5.0,2.0])
end

println("Evaluating the model with latency, refractory and bump parameters")
evaluate!(m2,[2.7,0.5,5.0,2.0,4.0,3.0,2.5])
for i=1:5
    gc() ; @time evaluate!(m2,[2.7,0.5,5.0,2.0,4.0,3.0,2.5])
end

println("Evaluating the model with latency and refractory parameters 1000 times")
gc()
@time for i=1:1000
    evaluate!(m1,[2.7,0.5,5.0,2.0,4.0,3.0,2.5])
end

results1 = evaluate!(m1,[2.7,0.5,5.0,2.0],1)
println("Evaluating the model 1000 times with latency and refractory parameters, storing the results")
@time results1 = evaluate!(m1,[2.7,0.5,5.0,2.0],1000)

println("Evaluating the model 1000 times with latency, refractory and bump parameters, storing the results")
@time results2 = evaluate!(m2,[2.7,0.5,5.0,2.0,4.0,3.0,2.5],1000)


text:
@testset "statistics.jl" begin
    @testset "TalliedStats" begin
        global tallied = Variable{Int64}(data=:tally, history=false)

        @process test_process(value_durations) begin
            for (value, duration) in value_durations
                tallied.value = value
                work(duration)
            end
        end

        @simulation begin
            value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

            @schedule now test_process(value_durations)
            @schedule at 1.0 test_process(value_durations)
            @schedule in 3.0 test_process(value_durations)

            start_simulation()
        end

        @test tallied.stats.min == 1
        @test tallied.stats.max == 4
        @test tallied.stats.sum == 30.0
        @test tallied.stats.mean == 2.5
        @test tallied.stats.variance == 1.25
        @test tallied.stats.stddev == 1.118033988749895
        @test tallied.stats.sum_squares == 90.0
    end

    @testset "AccumulatedStats" begin
        global tallied = Variable{Int64}(0, history=false)

        @process test_process(value_durations) begin
            for (value, duration) in value_durations
                tallied.value = value
                work(duration)
            end
        end

        @simulation begin
            value_durations = [(1, 2.0), (2, 1.0), (3, 2.0), (4, 3.0)]

            @schedule now test_process(value_durations)
            @schedule at 1.0 test_process(value_durations)
            @schedule in 3.0 test_process(value_durations)

            start_simulation()
        end

        @test tallied.stats.min == 1
        @test tallied.stats.max == 4
        @test tallied.stats.sum == 32.0
        @test tallied.stats.mean == 2.909090909090909
        @test tallied.stats.variance == 1.1735537190082646
        @test tallied.stats.stddev == 1.083306844346635
        @test tallied.stats.sum_squares == 106.0
    end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using DelimitedFiles
using Random
Random.seed!(233)

function halo_exchange_neighbor_two(u,fill_value,m,n,tag,w)
    halo_exchange_neighbor_two_ = load_op_and_grad("./build/libHaloExchangeNeighborTwo","halo_exchange_neighbor_two")
    u,fill_value,m,n,tag,w = convert_to_tensor(Any[u,fill_value,m,n,tag,w], [Float64,Float64,Int64,Int64,Int64,Float64])
    halo_exchange_neighbor_two_(u,fill_value,m,n,tag,w)
end

mpi_init()
u = rand(5, 5)
fill_value = 10.0
m = 1
n = 1
tag = 1
w = 1.0
# TODO: specify your input parameters
uext = halo_exchange_neighbor_two(u,fill_value,m,n,tag,w)
sess = Session(); init(sess)
uval =  run(sess, uext)
writedlm(stdout, round.(uval, digits=3))
# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(halo_exchange_neighbor_two(x,fill_value,m,n,tag,w)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(5, 5))
v_ = rand(5, 5)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")
savefig("test.png")


text:
module Market
using DataFrames

export Demand, LogitDemand, NestedLogitDemand
export estimate

abstract type Demand end
abstract type LogitDemand2 <: Demand end

struct LogitDemand <: Demand
   foo::String
   bar::Int
end

struct NestedLogitDemand <: Demand
   foo::String
   bar::Int
end

function estimate(x::LogitDemand)
   foo = x.foo
   "foo: $(x.foo)"
end

function estimate(x::NestedLogitDemand)
   "foo: $(x.foo)"
end

end # module


text:
module CUTENSOR

using ..CuArrays
using ..CuArrays: libcutensor, @argout, @workspace

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const handles_lock = ReentrantLock()
const created_handles = Dict{Tuple{UInt,Int},Ref{cutensorHandle_t}}()
const active_handles = Vector{Union{Nothing,Ref{cutensorHandle_t}}}()

function handle()
    tid = Threads.threadid()
    if @inbounds active_handles[tid] === nothing
        ctx = context()
        key = (objectid(ctx), tid)
        lock(handles_lock) do
            active_handles[tid] = get!(created_handles, key) do
                handle = Ref{cutensorHandle_t}()
                cutensorInit(handle)
                handle
            end
        end
    end
    @inbounds active_handles[tid]
end

function __init__()
    resize!(active_handles, Threads.nthreads())
    fill!(active_handles, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        # we don't eagerly initialize handles, but do so lazily when requested
        active_handles[tid] = nothing
    end
end

end


text:
include("BCQ.jl")
include("BEAR.jl")
include("behavior_cloning.jl")
include("CRR.jl")
include("DiscreteBCQ.jl")
include("FisherBRC.jl")
include("PLAS.jl")
include("common.jl")


text:
#!/usr/bin/env julia

using Test

using NodalDG.jl

for file in ["test-utils.jl", "test-mesh1d.jl", "test-grid1d.jl"]
    println("Tests in ", file)
    @time include(file)
end


text:
function commander(c::Client, m::Message, ::Val{:discourse})
    @info "discourse_commander called"
    startswith(m.content, COMMAND_PREFIX * "discourse") || return

    regex = Regex(COMMAND_PREFIX * raw"discourse( .*)$")
    matches = match(regex, m.content)
    if matches === nothing
        help_commander(c, m, :discourse)
        return
    end
    arg = strip(matches.captures[1])
    if arg == "latest"
        endpoint = "posts.json"
        params = Dict{Symbol,String}()
        topic_fn = x -> x.latest_posts
        make_post_fn = x -> DiscoursePost(x.topic_id, x.topic_slug)
    elseif arg == "top"
        endpoint = "top/weekly.json"
        params = Dict{Symbol,String}()
        topic_fn = x -> x.topic_list.topics
        make_post_fn = x -> DiscoursePost(x.id, x.slug)
    else
        endpoint = "search.json"
        params = Dict(:q => arg)
        topic_fn = x -> x.topics
        make_post_fn = x -> DiscoursePost(x.id, x.slug)
    end
    discourse_search(c, m, endpoint, params, topic_fn, make_post_fn)
    return nothing
end

function help_commander(c::Client, m::Message, ::Val{:discourse})
    reply(c, m, """
        Perform a search on Julia Discourse.
        ```
        discourse <query>
        discourse top
        discourse latest
        ```
        """)
end

"""
    discourse_search(
        c::Client,
        m::Message,
        endpoint::AbstractString,
        params::Dict{Symbol, S},
        topic_fn::Function,
        make_post_fn::Function
    ) where {S <: AbstractString}

Perform a search on Julia Discourse with the provided `endpoint`
and query parameters `params`. Use `topic_fn` to index into the
JSON array and `make_post_fn` to construct a `DiscoursePost`
object for each hit.

Create a Discord message with the result and save the data.
"""
function discourse_search(
    c::Client,
    m::Message,
    endpoint::AbstractString,
    params::Dict{Symbol, S},
    topic_fn::Function,
    make_post_fn::Function
) where {S <: AbstractString}
    data = discourse_execute_query(endpoint, params, topic_fn, make_post_fn)
    new_message = reply(c, m, data)
    discourse_save(new_message.id, data)
end

"""
    discourse_execute_query(
        endpoint::AbstractString,
        params::Dict{Symbol, AbstractString},
        topic_fn::Function,
        make_post_fn::Function
    )

Execute a query and return `DiscourseData` object.
"""
function discourse_execute_query(
    endpoint::AbstractString,
    params::Dict{Symbol, S},
    topic_fn::Function,
    make_post_fn::Function
) where {S <: AbstractString}
    json = discourse_run_query(endpoint, params)
    return DiscourseData(topic_fn(json), make_post_fn)
end

"""
    discourse_run_query(endpoint::AbstractString, params::Dict{Symbol, AbstractString})

Run Discourse query and return result as a JSON3 object.
"""
function discourse_run_query(endpoint::AbstractString, params::Dict{Symbol, S}) where {
    S <: AbstractString
}
    sanitized_kv = Dict(HTTP.escapeuri(k) => HTTP.escapeuri(v) for (k,v) in params)
    params_str = join(["$k=$v" for (k,v) in sanitized_kv], "&")
    response = HTTP.get(
        "$DISCOURSE_DOMAIN/$endpoint?$params_str",
        headers = ["Accept" => "application/json"])
    return JSON3.read(response.body)
end

function discourse_save(message_id::UInt64, r::DiscourseData)
    @info "Saving discourse search data for message_id=$message_id"
    path = discourse_file_path(message_id)
    write(ensurepath!(path), JSON3.write(r))
end

function discourse_load(message_id::UInt64)
    @info "Loading discourse search data for message_id=$message_id"
    path = discourse_file_path(message_id)
    bytes = read(path)
    return JSON3.read(bytes, DiscourseData)
end

"Location of the data file"
discourse_file_path(message_id::UInt64) = "data/discourse/$message_id.json"

"""
Reply to user query with a new message containing the current post
from the Discourse search results.
"""
function Discord.reply(c::Client, m::Message, r::DiscourseData)
    @info "Creating result message for message_id=$(m.id)"

    new_message = @discord reply(c, m, message(r))
    @info "New message_id=$(new_message.id)"
    create(c, Reaction, new_message, "👈")
    sleep(0.01) # try to order these emojis properly with a delay
    create(c, Reaction, new_message, "👉")
    return new_message
end


text:
using StochasticDiffEq, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_lorenz
Random.seed!(100)
prob = prob_sde_lorenz

## Solve and plot
println("Plot 1")
sol1 =solve(prob,SRI();dt=1/2^(4),abstol=10,reltol=0)
#=
p1 = plot(sol1[:,1],sol1[:,2],sol1[:,3],title="Absolute Tolerance = 10",leg=false,
          top_margin=50px,right_margin=50px,xguide="X",yguide="Y",zguide="Z",bottom_margin=50px,
          guidefont=font(16),titlefont=font(18),tickfont=font(16))
=#

println("Plot 2")
sol2 =solve(prob,SRI();dt=1/2^(4),abstol=1,reltol=0)
#=
p2 = plot(sol2[:,1],sol2[:,2],sol2[:,3],title="Absolute Tolerance = 10",leg=false,
          top_margin=50px,right_margin=50px,xguide="X",yguide="Y",zguide="Z",bottom_margin=50px,
          guidefont=font(16),titlefont=font(18),tickfont=font(16))
=#

println("Plot 3")
sol3 =solve(prob,SRI();dt=1/2^(4),abstol=1/10,reltol=0)
#=
p3 = plot(sol3[:,1],sol3[:,2],sol3[:,3],title="Absolute Tolerance = 10",leg=false,
          top_margin=50px,right_margin=50px,xguide="X",yguide="Y",zguide="Z",bottom_margin=50px,
          guidefont=font(16),titlefont=font(18),tickfont=font(16))
=#

println("Plot 4")
sol4 =solve(prob,SRI();dt=1/2^(4),abstol=1/100,reltol=0)
#=
p4 = plot(sol4[:,1],sol4[:,2],sol4[:,3],title="Absolute Tolerance = 10",leg=false,
          top_margin=50px,right_margin=50px,xguide="X",yguide="Y",zguide="Z",bottom_margin=50px,
          guidefont=font(16),titlefont=font(18),tickfont=font(16))
=#

#plot(p1,p2,p3,p4,size=(1200,800),plot_title="Lorenz Attractors")
#gui()

true


text:
# this doesn't test `hijack` stuff, just `load` stuff. We put it here
# because we test Revise stuff
module HijackTestsLoad1
f() = 1

# HijackTestsLoad1 is returned by load() even if only a submodule uses ReTest
module Sub
using ReTest
@testset "HijackTestsLoad1" begin end
end
end

module HijackTestsLoad2
f() = 1
using ReTest
@testset "HijackTestsLoad2" begin end
end

# not using ReTest
module HijackTestsLoad3
f() = 1
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# # An unstable Bickley jet in Shallow Water model
#
# This example uses Oceananigans.jl's `ShallowWaterModel` to simulate
# the evolution of an unstable, geostrophically balanced, Bickley jet
# The example is periodic in ``x`` with flat bathymetry and
# uses the conservative formulation of the shallow water equations.
# The initial conditions superpose the Bickley jet with small-amplitude perturbations.
# See ["The nonlinear evolution of barotropically unstable jets," J. Phys. Oceanogr. (2003)](https://doi.org/10.1175/1520-0485(2003)033<2173:TNEOBU>2.0.CO;2)
# for more details on this problem.
#
# The mass transport ``(uh, vh)`` is the prognostic momentum variable
# in the conservative formulation of the shallow water equations,
# where ``(u, v)`` are the horizontal velocity components and ``h``
# is the layer height.
#
# ## Install dependencies
#
# First we make sure that we have all of the packages that are required to
# run the simulation.
#
# ```julia
# using Pkg
# pkg"add Oceananigans, NCDatasets, Plots, Printf, Polynomials"
# ```

using Oceananigans
using Oceananigans.Models: ShallowWaterModel

# ## Two-dimensional domain
#
# The shallow water model is a two-dimensional model and thus the number of vertical
# points `Nz` must be set to one.  Note that ``L_z`` is the mean depth of the fluid.

Lx, Ly, Lz = 2π, 20, 1
Nx, Ny = 128, 128

grid = RectilinearGrid(size = (Nx, Ny),
                       x = (0, Lx), y = (-Ly/2, Ly/2),
                       topology = (Periodic, Bounded, Flat))

# ## Physical parameters
#
# We choose non-dimensional parameters

const U = 1.0         # Maximum jet velocity

f = 1           # Coriolis parameter
g = 9.8         # Gravitational acceleration
Δη = f * U / g  # Maximum free-surface deformation as dictated by geostrophy

# ## Building a `ShallowWaterModel`
#
# We build a `ShallowWaterModel` with the `WENO5` advection scheme and
# 3rd-order Runge-Kutta time-stepping,

model = ShallowWaterModel(
                          timestepper = :RungeKutta3,
                          advection = WENO5(),
                          grid = grid,
                          gravitational_acceleration = g,
                          coriolis = FPlane(f=f))

# Use `architecture = GPU()` to run this problem on a GPU.

# ## Background state and perturbation
#
# The background velocity ``ū`` and free-surface ``η̄`` correspond to a
# geostrophically balanced Bickely jet with maximum speed of ``U`` and maximum
# free-surface deformation of ``Δη``,

h̄(x, y, z) = Lz - Δη * tanh(y)
ū(x, y, z) = U * sech(y)^2

# The total height of the fluid is ``h = L_z + \eta``. Linear stability theory predicts that
# for the parameters we consider here, the growth rate for the most unstable mode that fits
# our domain is approximately ``0.139``.

# The vorticity of the background state is

ω̄(x, y, z) = 2 * U * sech(y)^2 * tanh(y)

# The initial conditions include a small-amplitude perturbation that decays away from the
# center of the jet.

small_amplitude = 1e-4

 uⁱ(x, y, z) = ū(x, y, z) + small_amplitude * exp(-y^2) * randn()
uhⁱ(x, y, z) = uⁱ(x, y, z) * h̄(x, y, z)

# We first set a "clean" initial condition without noise for the purpose of discretely
# calculating the initial 'mean' vorticity,

ū̄h(x, y, z) = ū(x, y, z) * h̄(x, y, z)

set!(model, uh = ū̄h, h = h̄)

# We next compute the initial vorticity and perturbation vorticity,

uh, vh, h = model.solution

## Build velocities
u = uh / h
v = vh / h

## Build and compute mean vorticity discretely
ω = ComputedField(∂x(v) - ∂y(u))
compute!(ω)

## Copy mean vorticity to a new field
ωⁱ = Field(Face, Face, Nothing, model.architecture, model.grid)
ωⁱ .= ω

## Use this new field to compute the perturbation vorticity
ω′ = ComputedField(ω - ωⁱ)

# and finally set the "true" initial condition with noise,

set!(model, uh = uhⁱ)

# ## Running a `Simulation`
#
# We pick the time-step so that we make sure we resolve the surface gravity waves, which
# propagate with speed of the order ``\sqrt{g L_z}``. That is, with `Δt = 1e-2` we ensure
# that `` \sqrt{g L_z} Δt / Δx,  \sqrt{g L_z} Δt / Δy < 0.7``.

simulation = Simulation(model, Δt = 1e-2, stop_time = 150)

# ## Prepare output files
#
# Define a function to compute the norm of the perturbation on the cross channel velocity.
# We obtain the `norm` function from `LinearAlgebra`.

using LinearAlgebra: norm

perturbation_norm(args...) = norm(v)

# Build the `output_writer` for the two-dimensional fields to be output.
# Output every `t = 1.0`.

simulation.output_writers[:fields] = NetCDFOutputWriter(model, (; ω, ω′),
                                                        filepath = joinpath(@__DIR__, "shallow_water_Bickley_jet_fields.nc"),
                                                        schedule = TimeInterval(1),
                                                        mode = "c")

# Build the `output_writer` for the growth rate, which is a scalar field.
# Output every time step.

simulation.output_writers[:growth] = NetCDFOutputWriter(model, (; perturbation_norm),
                                                        filepath = joinpath(@__DIR__, "shallow_water_Bickley_jet_perturbation_norm.nc"),
                                                        schedule = IterationInterval(1),
                                                        dimensions = (; perturbation_norm = ()),
                                                        mode = "c")

# And finally run the simulation.

run!(simulation)

# ## Visualize the results

# Load required packages to read output and plot.

using NCDatasets, Plots, Printf
nothing # hide

# Define the coordinates for plotting.

x, y = xnodes(ω), ynodes(ω)
nothing # hide

# Define keyword arguments for plotting the contours.

kwargs = (
         xlabel = "x",
         ylabel = "y",
         aspect = 1,
           fill = true,
         levels = 20,
      linewidth = 0,
          color = :balance,
       colorbar = true,
           ylim = (-Ly/2, Ly/2),
           xlim = (0, Lx)
)
nothing # hide

# Read in the `output_writer` for the two-dimensional fields and then create an animation
# showing both the total and perturbation vorticities.

ds = NCDataset(simulation.output_writers[:fields].filepath, "r")

iterations = keys(ds["time"])

anim = @animate for (iter, t) in enumerate(ds["time"])
    ω = ds["ω"][:, :, 1, iter]
    ω′ = ds["ω′"][:, :, 1, iter]

    ω′_max = maximum(abs, ω′)

    plot_ω = contour(x, y, ω',
                     clim = (-1, 1),
                     title = @sprintf("Total vorticity, ω, at t = %.1f", t); kwargs...)

    plot_ω′ = contour(x, y, ω′',
                      clim = (-ω′_max, ω′_max),
                      title = @sprintf("Perturbation vorticity, ω - ω̄, at t = %.1f", t); kwargs...)

    plot(plot_ω, plot_ω′, layout = (1, 2), size = (800, 440))
end

close(ds)

mp4(anim, "shallow_water_Bickley_jet.mp4", fps=15)

# Read in the `output_writer` for the scalar field (the norm of ``v``-velocity).

ds2 = NCDataset(simulation.output_writers[:growth].filepath, "r")

iterations = keys(ds2["time"])

     t = ds2["time"][:]
norm_v = ds2["perturbation_norm"][:]

close(ds2)
nothing # hide

# We import the `fit` function from `Polynomials.jl` to compute the best-fit slope of the
# perturbation norm on a logarithmic plot. This slope corresponds to the growth rate.

using Polynomials: fit

I = 6000:7000

degree = 1
linear_fit_polynomial = fit(t[I], log.(norm_v[I]), degree, var = :t)

# We can get the coefficient of the ``n``-th power from the fitted polynomial by using `n`
# as an index, e.g.,

constant, slope = linear_fit_polynomial[0], linear_fit_polynomial[1]

# We then use the computed linear fit coefficients to construct the best fit and plot it
# together with the time-series for the perturbation norm for comparison.

best_fit = @. exp(constant + slope * t)

plot(t, norm_v,
        yaxis = :log,
        ylims = (1e-3, 30),
           lw = 4,
        label = "norm(v)",
       xlabel = "time",
       ylabel = "norm(v)",
        title = "growth of perturbation norm",
       legend = :bottomright)

plot!(t[I], 2 * best_fit[I], # factor 2 offsets fit from curve for better visualization
           lw = 4,
        label = "best fit")

# The slope of the best-fit curve on a logarithmic scale approximates the rate at which instability
# grows in the simulation. Let's see how this compares with the theoretical growth rate.

println("Numerical growth rate is approximated to be ", round(slope, digits=3), ",\n",
        "which is very close to the theoretical value of 0.139.")


text:
module ScoringTests
  using Base.Test
  include("../src/scoring.jl")

  # Offensive scoring
  @test score_off() == 0
  @test score_off(rush_yds=1) == 0
  @test score_off(rush_yds=10) == 1
  @test score_off(rush_td=1) == 6
  @test score_off(rec_yds=1) == 0
  @test score_off(rec_yds=10) == 1
  @test score_off(rec_td=1) == 6
  @test score_off(pass_yds=10) == 0
  @test score_off(pass_yds=25) == 1
  @test score_off(pass_td=1) == 4
  @test score_off(pass_int=1) == -2
  @test score_off(sack=1) == 0
  @test score_off(fum_l=1) == -2
  @test score_off(pr_td=1) == 6
  @test score_off(kr_td=1) == 6

  @test score_off(rush_yds=10, rush_td=1, rec_yds=10, rec_td=1, pass_yds=25, pass_td=1, pass_int=1, sack=1, fum_l=1, pr_td=1, kr_td=1) == (1 + 6 + 1 + 6 + 1 + 4 + -2 + 0 + -2 + 6 + 6)

  # Kicker scoring
  @test score_k() == 0
  @test score_k(xp_made=1) == 1
  @test score_k(xp_att=1) == -1
  @test score_k(fg_made=1) == 0
  @test score_k(fg_att=1) == -1
  @test score_k(f0019=1) == 3
  @test score_k(f2029=1) == 3
  @test score_k(f3039=1) == 3
  @test score_k(f4049=1) == 4
  @test score_k(f50=1) == 5

  @test score_k(xp_made=1, xp_att=2, fg_made=1, fg_att=2, f0019=1, f2029=1, f3039=1, f4049=1, f50=1) == (1 + -1 + -1 + 3 + 3 + 3 + 4 + 5)

  # Defense / Special Teams scoring
  @test score_def_st() == 0
  @test score_def_st(games=1) == 5
  @test score_def_st(games=2) == 10
  @test score_def_st(papg=1) == 0
  @test score_def_st(games=1, papg=0) == 5
  @test score_def_st(games=1, papg=6) == 4
  @test score_def_st(games=1, papg=13) == 3
  @test score_def_st(games=1, papg=17) == 1
  @test score_def_st(games=1, papg=27) == 0
  @test score_def_st(games=1, papg=34) == -1
  @test score_def_st(games=1, papg=45) == -3
  @test score_def_st(games=1, papg=46) == -5
  @test score_def_st(kr_td=1) == 6
  @test score_def_st(pr_td=1) == 6
  @test score_def_st(int_td=1) == 6
  @test score_def_st(fum_td=1) == 6
  @test score_def_st(sack=1) == 1
  @test score_def_st(sfty=1) == 2

  @test score_def_st(games=1, papg=1, kr_td=1, pr_td=1, int_td=1, fum_td=1, sack=1, sfty=1) == (4 + 6 + 6 + 6 + 6 + 1 + 2)

end

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
using GAMS
using Test

ws = GAMSWorkspace()

const TEST_GAMSLIB_MODELS = (
   "trnsport",
   "blend",
   "prodmix",
   "whouse",
   "jobt",
   "sroute",
   "diet",
   "aircraft",
   "prodsch",
   "pdi",
   "uimp",
   "magic",
   "ferts",
   "fertd",
   "mexss",
   "mexsd",
   "mexls",
   "weapons",
   "bid",
   "process",
   "chem",
   "ship",
   "linear",
   "least",
   "like",
   "chance",
   "sample",
   "pindyck",
   "vietman",
   "marco",
   "chenery",
   "pak",
   "himmel16",
   "robert",
   "rdata",
   "mine",
   "orani",
   "cube",
   "chakra",
   "andean",
   "copper",
   "otpop",
   "korpet",
   "sarf",
   "port",
   "prodschx",
   "bidsos",
   "wallmcp",
   "transmcp",
   # "nash"
)

# we need GAMS 32 for JuMP output in Convert
if GAMS.get_version(ws)[1] < 32
   return
end

# create temporary working directory
curdir = pwd()
tempdir = mktempdir(prefix = "gams_jl_")
cd(tempdir)

# prepare GAMS Convert option file to write JuMP model
open("convert.opt", "w") do io
   println(io, "jump jump.jl")
end

@testset "gamslib        " begin
   for model in TEST_GAMSLIB_MODELS
      @testset "$model" begin

         # get model and convert to scalar GAMS and JuMP model
         Base.run(`gamslib -q $model`)
         Base.run(`gams $model.gms lo=2 solver=convert`)
         Base.run(`gams $model.gms lo=2 solver=convert optfile=1`)

         # instruct GAMS model to export objective function
         open("gams.gms", "a") do io
            println(io, "File result_file / objval.txt /;")
            println(io, "result_file.nd = 7;")
            println(io, "put result_file;")
            println(io, "put m.objval;")
            println(io, "putclose;")
         end

         # solve GAMS model
         objval_gams = NaN;
         Base.run(`gams gams.gms lo=2`)
         open("objval.txt", "r") do io
            objval_gams = parse(Float64, read(io, String))
         end

         # solve JuMP model using GAMS
         include(joinpath(tempdir, "jump.jl"))
         set_optimizer(m, GAMS.Optimizer)
         set_optimizer_attribute(m, GAMS.WorkDir(), tempdir)
         set_optimizer_attribute(m, MOI.Silent(), true)
         optimize!(m)
         objval_jump = objective_value(m)

         @test objval_gams ≈ objval_jump rtol=1e-4
      end
   end
end

cd(curdir)


text:
module ActionsTestPackage

greet() = print("Hello World!")

function cov(b::Bool)
    if b
        @info "This branch is hit in tests."
    else
        @info "This branch isn't hit in tests."
    end
    b
end

end # module


text:
# the package is otherwise pretty lightweight, would be nice not to depend on Flux
export Conv1d
using Flux: glorot_uniform, Flux

struct Conv1d{F,A,V}
    σ::F
    weight::A
    bias::V
end

function Conv1d(width::Integer, cin_cout::Pair, σ = identity; init = glorot_uniform)
    weight = Flux.convfilter((width,), cin_cout, init=init)
    bias = similar(weight, cin_cout[2])
    fill!(bias, zero(eltype(bias)))
    return Conv1d(σ, weight, bias)
end

Flux.trainable(o::Conv1d) = (o.weight, o.bias)

(o::Conv1d)(x) = o.σ.(conv(x, o.weight) .+ o.bias)


text:
module AWS4AuthRequest

using ..Base64
using ..Dates
using MbedTLS: digest, MD_SHA256, MD_MD5
import ..Layer, ..request, ..Headers
using ..URIs
using ..Pairs: getkv, setkv, rmkv
import ..@debug, ..DEBUG_LEVEL

"""
    request(AWS4AuthLayer, ::URI, ::Request, body) -> HTTP.Response

Add a [AWS Signature Version 4](http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html)
`Authorization` header to a `Request`.

Credentials are read from environment variables `AWS_ACCESS_KEY_ID`,
`AWS_SECRET_ACCESS_KEY` and `AWS_SESSION_TOKEN`.
"""
abstract type AWS4AuthLayer{Next <: Layer} <: Layer end
export AWS4AuthLayer

function request(::Type{AWS4AuthLayer{Next}},
                 url::URI, req, body; kw...) where Next

    @static if VERSION > v"0.7.0-DEV.2915"
    if !haskey(kw, :aws_access_key_id) &&
       !haskey(ENV, "AWS_ACCESS_KEY_ID")
        kw = merge(dot_aws_credentials(), kw)
    end
    end

    sign_aws4!(req.method, url, req.headers, req.body; kw...)

    return request(Next, url, req, body; kw...)
end

function sign_aws4!(method::String,
                    url::URI,
                    headers::Headers,
                    body::Vector{UInt8};
                    body_sha256::Vector{UInt8}=digest(MD_SHA256, body),
                    body_md5::Vector{UInt8}=digest(MD_MD5, body),
                    t::DateTime=now(Dates.UTC),
                    aws_service::String=String(split(url.host, ".")[1]),
                    aws_region::String=String(split(url.host, ".")[2]),
                    aws_access_key_id::String=ENV["AWS_ACCESS_KEY_ID"],
                    aws_secret_access_key::String=ENV["AWS_SECRET_ACCESS_KEY"],
                    aws_session_token::String=get(ENV, "AWS_SESSION_TOKEN", ""),
                    kw...)

    # ISO8601 date/time strings for time of request...
    date = Dates.format(t, dateformat"yyyymmdd")
    datetime = Dates.format(t, dateformat"yyyymmddTHHMMSS\Z")

    # Authentication scope...
    scope = [date, aws_region, aws_service, "aws4_request"]

    # Signing key generated from today's scope string...
    signing_key = string("AWS4", aws_secret_access_key)
    for element in scope
        signing_key = digest(MD_SHA256, element, signing_key)
    end

    # Authentication scope string...
    scope = join(scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(body_sha256)

    # HTTP headers...
    rmkv(headers, "Authorization")
    setkv(headers, "x-amz-content-sha256",  content_hash)
    setkv(headers, "x-amz-date",  datetime)
    setkv(headers, "Content-MD5", base64encode(body_md5))
    if aws_session_token != ""
        setkv(headers, "x-amz-security-token", aws_session_token)
    end

    # Sort and lowercase() Headers to produce canonical form...
    canonical_headers = ["$(lowercase(k)):$(strip(v))" for (k,v) in headers]
    signed_headers = join(sort([lowercase(k) for (k,v) in headers]), ";")

    # Sort Query String...
    query = queryparams(url.query)
    query = Pair[k => query[k] for k in sort(collect(keys(query)))]

    # Create hash of canonical request...
    canonical_form = string(method, "\n",
                            aws_service == "s3" ? url.path
                                                : escapepath(url.path), "\n",
                            escapeuri(query), "\n",
                            join(sort(canonical_headers), "\n"), "\n\n",
                            signed_headers, "\n",
                            content_hash)
    @debug 3 "AWS4 canonical_form: $canonical_form"

    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    @debug 3 "AWS4 string_to_sign: $string_to_sign"
    @debug 3 "AWS4 signature: $signature"

    # Append Authorization header...
    setkv(headers, "Authorization", string(
        "AWS4-HMAC-SHA256 ",
        "Credential=$aws_access_key_id/$scope, ",
        "SignedHeaders=$signed_headers, ",
        "Signature=$signature"
    ))
end

@static if VERSION > v"0.7.0-DEV.2915"

using IniFile

credentials = NamedTuple()

"""
Load Credentials from [AWS CLI ~/.aws/credentials file]
(http://docs.aws.amazon.com/cli/latest/userguide/cli-config-files.html).
"""
function dot_aws_credentials()::NamedTuple

    global credentials
    if !isempty(credentials)
        return credentials
    end

    f = get(ENV, "AWS_CONFIG_FILE", joinpath(homedir(), ".aws", "credentials"))
    p = get(ENV, "AWS_DEFAULT_PROFILE", get(ENV, "AWS_PROFILE", "default"))

    if !isfile(f)
        return NamedTuple()
    end

    ini = read(Inifile(), f)

    credentials = (
        aws_access_key_id = String(get(ini, p, "aws_access_key_id")),
        aws_secret_access_key = String(get(ini, p, "aws_secret_access_key")))
end

end

end # module AWS4AuthRequest


text:
# Original implementations by @nantonel
# https://github.com/maxmouchet/HMMBase.jl/pull/6

# *log! methods use the samples log-likelihood instead of the likelihood.

function viterbi!(T1::AbstractMatrix, T2::AbstractMatrix, z::AbstractVector, a::AbstractVector, A::AbstractMatrix, L::AbstractMatrix)
    @argcheck size(T1, 1) == size(T2, 1) == size(L, 1) == size(z, 1)
    @argcheck size(T1, 2) == size(T2, 2) == size(L, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(L)
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)

    c = 0.0

    for i in OneTo(K)
        T1[1,i] = a[i] * L[1,i]
        c += T1[1,i]
    end

    for i in OneTo(K)
        T1[1,i] /= c
    end

    @inbounds for t in 2:T
        c = 0.0

        for j in OneTo(K)
            # TODO: If there is NaNs in T1 this may
            # stay to 0 (NaN > -Inf == false).
            # Hence it will crash when computing z[t].
            # Maybe we should check for NaNs beforehand ?
            amax = 0
            vmax = -Inf

            for i in OneTo(K)
                v = T1[t-1,i] * A[i,j]
                if v > vmax
                    amax = i
                    vmax = v
                end
            end

            T1[t,j] = vmax * L[t,j]
            T2[t,j] = amax
            c += T1[t,j]
        end

        for i in OneTo(K)
            T1[t,i] /= c
        end
    end

    z[T] = argmax(T1[T,:])
    for t in T-1:-1:1
        z[t] = T2[t+1,z[t+1]]
    end
end

function viterbilog!(T1::AbstractMatrix, T2::AbstractMatrix, z::AbstractVector, a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix)
    T, K = size(LL)
    (T == 0) && return

    fill!(T1, 0.0)
    fill!(T2, 0)

    al = log.(a)
    Al = log.(A)

    for i in OneTo(K)
        T1[1,i] = al[i] + LL[1,i]
    end

    @inbounds for t in 2:T
        for j in OneTo(K)
            amax = 0
            vmax = -Inf

            for i in OneTo(K)
                v = T1[t-1,i] + Al[i,j]
                if v > vmax
                    amax = i
                    vmax = v
                end
            end

            T1[t,j] = vmax + LL[t,j]
            T2[t,j] = amax
        end
    end

    z[T] = argmax(T1[T,:])
    for t in T-1:-1:1
        z[t] = T2[t+1,z[t+1]]
    end
end


text:
using MyExampleJl
using Test

@testset "MyExampleJl.jl" begin
    # Write your own tests here.
end


text:
#InspectDR: Numeric formatting
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
import NumericIO: UTF8_MINUS_SYMBOL, UTF8_INF_STRING
import NumericIO: UEXPONENT

#==Main types
===============================================================================#

#Convenience aliases
const NumericFormatting = NumericIO.IOFormattingReal
const ExponentFormatting = NumericIO.IOFormattingExp

abstract type RangeDisplayInfo end

#Hints on how to display numbers in a linear range:
struct LinearRangeDisplayInfo <: RangeDisplayInfo
	exp_min::Float64 #Base-10 exponent for max value
	exp_max::Float64 #Base-10 exponent for min value
	exp_step::Float64 #Base-10 exponent for step size
end
LinearRangeDisplayInfo() = LinearRangeDisplayInfo(0,0,0)

struct NoRangeDisplayInfo <: RangeDisplayInfo; end

#How to display tick labels:
struct TickLabelFormatting
	fmt::NumericFormatting
	splitexp::Bool #Only possible when !fmt.decfloating
end

mutable struct TickLabelStyle
	expdisplay::ExponentFormatting
	ndigits::NullOr{Int} #Force # of digits displayed
	decpos::NullOr{Int} #User-defined exponent for all displayed numbers
	commonexp::Bool #Use same exponent value for all displayed numbers
	showexp0::Bool
	splitexp::Bool #Only possible with Integer decpos or commonexp=true
end
TickLabelStyle(expdisplay::ExponentFormatting=UEXPONENT; ndigits=nothing,
		decpos=nothing, commonexp=true, showexp0=false, splitexp=true) =
	TickLabelStyle(expdisplay, ndigits, decpos, commonexp, showexp0, splitexp)

#==Helper functions
===============================================================================#
finiteornan(x::Float64) = isfinite(x) ? x : NaN

base10exp(v::Float64) = floor(log10(abs(v)))

#Get exponent string only:
function formatted_exp(fmt::NumericFormatting)
	s = IOBuffer()
	NumericIO.print_formatted_exp(s, fmt)
	d = s.data
	resize!(d, s.size)
	return String(d)
end

#==Other "constructors"
===============================================================================#

number_fmt(expdisplay::ExponentFormatting=UEXPONENT; ndigits=3,
		decpos=1, decfloating=true, eng=true, showexp0=false) =
	NumericIO.IOFormattingReal(expdisplay, ndigits=ndigits, decpos=decpos, decfloating=decfloating,
		showexp0=showexp0, eng=eng, minus=UTF8_MINUS_SYMBOL, inf=UTF8_INF_STRING,
	)

TickLabelFormatting(::NoRangeDisplayInfo) =
	TickLabelFormatting(number_fmt(ndigits=3, decfloating=true), false)

TickLabelFormatting(style::TickLabelStyle, rng::NoRangeDisplayInfo) =
	TickLabelFormatting(rng)

function TickLabelFormatting(style::TickLabelStyle, rng::LinearRangeDisplayInfo)
#TODO: Simplify/improve logic in this function.  Also: improve test coverage.
	decpos = style.decpos
	finite_expmin = finiteornan(rng.exp_min)
	finite_expmax = finiteornan(rng.exp_max)
	expmax = max(finite_expmin, finite_expmax)
	if !isfinite(expmax)
		expmax = rng.exp_step+1
	end

	decfloating = (nothing == decpos && !style.commonexp)
	splitexp = style.splitexp && !decfloating #Only possible when !decfloating

	ndigits = style.ndigits
	if nothing == ndigits
		ndigits = 3 #Default
		digitdelta1 = abs(expmax-rng.exp_step)
		digitdelta2 = 0
		if decpos != nothing
			digitdelta2 = decpos-rng.exp_step #In case decpos fixed much higher than data
		end
		digitdelta = max(digitdelta1, digitdelta2)
		if isfinite(digitdelta)
			ndigits = round(Int, digitdelta)+1
		end
	end

	if decfloating
		decpos = 0 #Formatting tool requires integer value
	elseif !decfloating
		decpos = expmax #Starting point for optimal decimal position.
		if !isfinite(decpos); decpos = 0; end #Algorithm requires finite decpos.
		decpos = round(Int, decpos) #Just in case
		decpos = floor(Int, decpos/3)*3

#=
		const MAXDIGITS = 5 #Before using exponential notation
		if decpos >= 0 && decpos <= MAXDIGITS && expmax <= MAXDIGITS
			decpos = 0
		elseif decpos < 0 && decpos >= MAXDIGITS && rng.exp_step >= MAXDIGITS
			decpos = 0
		end
=#
	end

	fmt = number_fmt(style.expdisplay, ndigits=ndigits,
		decpos=decpos, decfloating=decfloating, showexp0=style.showexp0
	)
	return TickLabelFormatting(fmt, splitexp)
end

#Last line


text:
# This file was generated, do not modify it. # hide
surrogate = Deterministic()
mach = machine(surrogate, Xs, ys; predict=yhat)

@from_network mach begin
    mutable struct OneHundredModels
        atom=atom
    end
end

one_hundred_models = OneHundredModels()


text:
#!/usr/bin/env julia

pnglist = filter!(r"\.png", readdir())

function expand(locci)
    list = []
    pos = split(locci, ':')
        chr = pos[1]
        region = split(pos[2], '-')
        n = parse(Int64, region[1])
        while n <= parse(Int64, region[end])
             out = "$chr" * ":" * "$n"
            push!(list, out)
            n +=  1
        end
    return list
end

function getSNPs(pnglist=pnglist)
	snps = []
	for img in pnglist
		img = replace(img, ".png", "")
		parts = split(img, '\_')
		chr = replace(parts[1], "chr", "")
		start = parse(Int, replace(parts[2], ",", ""))
		endp  = parse(Int, replace(parts[3], ",", ""))
		pos = "$chr" * ":" * "$start" * "-" * "$endp"
		vars = expand(pos)
		append!(snps, vars)
	end
	snps = unique(snps)
	for snp in snps
		println(snp)
	end
end

getSNPs()


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

using Random
isdefined(Main, :OffsetArrays) || @eval Main include("testhelpers/OffsetArrays.jl")
using .Main.OffsetArrays

# fold(l|r) & mapfold(l|r)
@test foldl(+, Int64[]) === Int64(0) # In reference to issues #7465/#20144 (PR #20160)
@test foldl(+, Int16[]) === Int16(0) # In reference to issues #21536
@test foldl(-, 1:5) == -13
@test foldl(-, 1:5; init=10) == -5

@test Base.mapfoldl(abs2, -, 2:5) == -46
@test Base.mapfoldl(abs2, -, 2:5; init=10) == -44

@test Base.mapfoldl(abs2, /, 2:5) ≈ 1/900
@test Base.mapfoldl(abs2, /, 2:5; init=10) ≈ 1/1440

@test Base.mapfoldl((x)-> x ⊻ true, &, [true false true false false]) == false
@test Base.mapfoldl((x)-> x ⊻ true, &, [true false true false false]; init=true) == false

@test Base.mapfoldl((x)-> x ⊻ true, |, [true false true false false]) == true
@test Base.mapfoldl((x)-> x ⊻ true, |, [true false true false false]; init=false) == true

@test foldr(+, Int64[]) === Int64(0) # In reference to issue #20144 (PR #20160)
@test foldr(+, Int16[]) === Int16(0) # In reference to issues #21536
@test foldr(-, 1:5) == 3
@test foldr(-, 1:5; init=10) == -7
@test foldr(+, [1]) == 1 # Issue #21493

@test Base.mapfoldr(abs2, -, 2:5) == -14
@test Base.mapfoldr(abs2, -, 2:5; init=10) == -4

@test foldr((x, y) -> ('⟨' * x * '|' * y * '⟩'), "λ 🐨.α") == "⟨λ|⟨ |⟨🐨|⟨.|α⟩⟩⟩⟩" # issue #31780
let x = rand(10)
    @test 0 == @allocated(sum(Iterators.reverse(x)))
    @test 0 == @allocated(foldr(-, x))
end

# reduce
@test reduce(+, Int64[]) === Int64(0) # In reference to issue #20144 (PR #20160)
@test reduce(+, Int16[]) === Int16(0) # In reference to issues #21536
@test reduce((x,y)->"($x+$y)", 9:11) == "((9+10)+11)"
@test reduce(max, [8 6 7 5 3 0 9]) == 9
@test reduce(+, 1:5; init=1000) == (1000 + 1 + 2 + 3 + 4 + 5)
@test reduce(+, 1) == 1

# mapreduce
@test mapreduce(-, +, [-10 -9 -3]) == ((10 + 9) + 3)
@test mapreduce((x)->x[1:3], (x,y)->"($x+$y)", ["abcd", "efgh", "01234"]) == "((abc+efg)+012)"

# mapreduce with multiple iterators
@test mapreduce(*, +, (i for i in 2:3), (i for i in 4:5)) == 23
@test mapreduce(*, +, (i for i in 2:3), (i for i in 4:5); init = 2) == 25
@test mapreduce(*, (x,y)->"($x+$y)", ["a", "b", "c"], ["d", "e", "f"]) == "((ad+be)+cf)"
@test mapreduce(*, (x,y)->"($x+$y)", ["a", "b", "c"], ["d", "e", "f"]; init = "gh") ==
    "(((gh+ad)+be)+cf)"

@test mapreduce(*, +, [2, 3], [4, 5]) == 23
@test mapreduce(*, +, [2, 3], [4, 5]; init = 2) == 25
@test mapreduce(*, +, [2, 3], [4, 5]; dims = 1) == [23]
@test mapreduce(*, +, [2, 3], [4, 5]; dims = 1, init = 2) == [25]
@test mapreduce(*, +, [2, 3], [4, 5]; dims = 2) == [8, 15]
@test mapreduce(*, +, [2, 3], [4, 5]; dims = 2, init = 2) == [10, 17]

@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]) == 110
@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]; init = 2) == 112
@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]; dims = 1) == [44 66]
@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]; dims = 1, init = 2) == [46 68]
@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]; dims = 2) == reshape([33, 77], :, 1)
@test mapreduce(*, +, [2 3; 4 5], [6 7; 8 9]; dims = 2, init = 2) == reshape([35, 79], :, 1)

# mapreduce() for 1- 2- and n-sized blocks (PR #19325)
@test mapreduce(-, +, [-10]) == 10
@test mapreduce(abs2, +, [-9, -3]) == 81 + 9
@test mapreduce(-, +, [-9, -3, -4, 8, -2]) == (9 + 3 + 4 - 8 + 2)
@test mapreduce(-, +, Vector(range(1.0, stop=10000.0, length=10000))) == -50005000.0
# empty mr
@test mapreduce(abs2, +, Float64[]) === 0.0
@test mapreduce(abs2, *, Float64[]) === 1.0
@test mapreduce(abs2, max, Float64[]) === 0.0
@test mapreduce(abs, max, Float64[]) === 0.0
@test_throws ArgumentError mapreduce(abs2, &, Float64[])
@test_throws ArgumentError mapreduce(abs2, |, Float64[])

# mapreduce() type stability
@test typeof(mapreduce(*, +, Int8[10])) ===
      typeof(mapreduce(*, +, Int8[10, 11])) ===
      typeof(mapreduce(*, +, Int8[10, 11, 12, 13]))
@test typeof(mapreduce(*, +, Float32[10.0])) ===
      typeof(mapreduce(*, +, Float32[10, 11])) ===
      typeof(mapreduce(*, +, Float32[10, 11, 12, 13]))
# mapreduce() type stability when f supports empty collections
@test typeof(mapreduce(abs, +, Int8[])) ===
      typeof(mapreduce(abs, +, Int8[10])) ===
      typeof(mapreduce(abs, +, Int8[10, 11])) ===
      typeof(mapreduce(abs, +, Int8[10, 11, 12, 13]))
@test typeof(mapreduce(abs, +, Float32[])) ===
      typeof(mapreduce(abs, +, Float32[10])) ===
      typeof(mapreduce(abs, +, Float32[10, 11])) ===
      typeof(mapreduce(abs, +, Float32[10, 11, 12, 13]))

# sum
@testset "sums promote to at least machine size" begin
    @testset for T in [Int8, Int16, Int32]
        @test sum(T[]) === Int(0)
    end
    @testset for T in [UInt8, UInt16, UInt32]
        @test sum(T[]) === UInt(0)
    end
    @testset for T in [Int, Int64, Int128, UInt, UInt64, UInt128,
                       Float16, Float32, Float64]
        @test sum(T[]) === T(0)
    end
    @test sum(BigInt[]) == big(0) && sum(BigInt[]) isa BigInt
end

@test sum(Bool[]) === sum(Bool[false]) === sum(Bool[false, false]) === 0
@test sum(Bool[true, false, true]) === 2

@test sum(Int8(3)) === Int(3)
@test sum(3) === 3
@test sum(3.0) === 3.0

@test sum([Int8(3)]) === Int(3)
@test sum([3]) === 3
@test sum([3.0]) === 3.0

z = reshape(1:16, (2,2,2,2))
fz = float(z)
@test sum(z) === 136
@test sum(fz) === 136.0

@test_throws ArgumentError sum(sin, Int[])
@test sum(sin, 3) == sin(3.0)
@test sum(sin, [3]) == sin(3.0)
a = sum(sin, z)
@test a ≈ sum(sin, fz)
@test a ≈ sum(sin.(fz))

z = [-4, -3, 2, 5]
fz = float(z)
a = randn(32) # need >16 elements to trigger BLAS code path
b = complex.(randn(32), randn(32))

# check variants of summation for type-stability and other issues (#6069)
sum2(itr) = invoke(sum, Tuple{Any}, itr)
plus(x,y) = x + y
sum3(A) = reduce(plus, A)
sum4(itr) = invoke(reduce, Tuple{Function, Any}, plus, itr)
sum5(A) = reduce(plus, A; init=0)
sum6(itr) = invoke(Core.kwfunc(reduce), Tuple{NamedTuple{(:init,), Tuple{Int}}, typeof(reduce), Function, Any}, (init=0,), reduce, plus, itr)
sum7(A) = mapreduce(x->x, plus, A)
sum8(itr) = invoke(mapreduce, Tuple{Function, Function, Any}, x->x, plus, itr)
sum9(A) = mapreduce(x->x, plus, A; init=0)
sum10(itr) = invoke(Core.kwfunc(mapreduce), Tuple{NamedTuple{(:init,),Tuple{Int}}, typeof(mapreduce), Function, Function, Any}, (init=0,), mapreduce, x->x, plus, itr)
for f in (sum2, sum5, sum6, sum9, sum10)
    @test sum(z) == f(z)
    @test sum(Int[]) == f(Int[]) == 0
    @test sum(Int[7]) == f(Int[7]) == 7
    @test typeof(f(Int8[])) == typeof(f(Int8[1])) == typeof(f(Int8[1 7]))
end
for f in (sum3, sum4, sum7, sum8)
    @test sum(z) == f(z)
    @test_throws ArgumentError f(Int[])
    @test sum(Int[7]) == f(Int[7]) == 7
end
@test typeof(sum(Int8[])) == typeof(sum(Int8[1])) == typeof(sum(Int8[1 7]))

# check sum(abs, ...) for support of empty collections
@testset "sum(abs, [])" begin
    @test @inferred(sum(abs, Float64[])) === 0.0
    @test @inferred(sum(abs, Int[])) === 0
    @test @inferred(sum(abs, Set{Int}())) === 0
    @test_throws MethodError sum(abs, Any[])
end

# prod

@test prod(Int[]) === 1
@test prod(Int8[]) === Int(1)
@test prod(Float64[]) === 1.0

@test prod([3]) === 3
@test prod([Int8(3)]) === Int(3)
@test prod([UInt8(3)]) === UInt(3)
@test prod([3.0]) === 3.0

@test prod(z) === 120
@test prod(fz) === 120.0

@test prod(1:big(16)) == big(20922789888000)
@test prod(big(typemax(Int64)):big(typemax(Int64))+16) == parse(BigInt,"25300281663413827620486300433089141956148633919452440329174083959168114253708467653081909888307573358090001734956158476311046124934597861626299416732205795533726326734482449215730132757595422510465791525610410023802664753402501982524443370512346073948799084936298007821432734720004795146875180123558814648586972474376192000")

@test typeof(prod(Array(trues(10)))) == Bool

# check type-stability
prod2(itr) = invoke(prod, Tuple{Any}, itr)
@test prod(Int[]) === prod2(Int[]) === 1
@test prod(Int[7]) === prod2(Int[7]) === 7
@test typeof(prod(Int8[])) == typeof(prod(Int8[1])) == typeof(prod(Int8[1, 7])) == Int
@test typeof(prod2(Int8[])) == typeof(prod2(Int8[1])) == typeof(prod2(Int8[1 7])) == Int

# maximum & minimum & extrema

@test_throws ArgumentError maximum(Int[])
@test_throws ArgumentError minimum(Int[])

@test maximum(5) == 5
@test minimum(5) == 5
@test extrema(5) == (5, 5)
@test extrema(abs2, 5) == (25, 25)

let x = [4,3,5,2]
    @test maximum(x) == 5
    @test minimum(x) == 2
    @test extrema(x) == (2, 5)

    @test maximum(abs2, x) == 25
    @test minimum(abs2, x) == 4
    @test extrema(abs2, x) == (4, 25)
end

@test maximum([-0.,0.]) === 0.0
@test maximum([0.,-0.]) === 0.0
@test maximum([0.,-0.,0.]) === 0.0
@test minimum([-0.,0.]) === -0.0
@test minimum([0.,-0.]) === -0.0
@test minimum([0.,-0.,0.]) === -0.0

@testset "minimum/maximum checks all elements" begin
    for N in [2:20;150;300]
        for i in 1:N
            arr = fill(0., N)
            truth = rand()
            arr[i] = truth
            @test maximum(arr) == truth

            truth = -rand()
            arr[i] = truth
            @test minimum(arr) == truth

            arr[i] = NaN
            @test isnan(maximum(arr))
            @test isnan(minimum(arr))

            arr = zeros(N)
            @test minimum(arr) === 0.0
            @test maximum(arr) === 0.0

            arr[i] = -0.0
            @test minimum(arr) === -0.0
            @test maximum(arr) ===  0.0

            arr = -zeros(N)
            @test minimum(arr) === -0.0
            @test maximum(arr) === -0.0
            arr[i] = 0.0
            @test minimum(arr) === -0.0
            @test maximum(arr) === 0.0
        end
    end
end

@testset "maximum works on generic order #30320" begin
    for n in [1:20;1500]
        arr = randn(n)
        @test GenericOrder(maximum(arr)) === maximum(map(GenericOrder, arr))
        @test GenericOrder(minimum(arr)) === minimum(map(GenericOrder, arr))
        f = x -> x
        @test GenericOrder(maximum(f,arr)) === maximum(f,map(GenericOrder, arr))
        @test GenericOrder(minimum(f,arr)) === minimum(f,map(GenericOrder, arr))
    end
end

@testset "maximum no out of bounds access #30462" begin
    arr = fill(-Inf, 128,128)
    @test maximum(arr) == -Inf
    arr = fill(Inf, 128^2)
    @test minimum(arr) == Inf
    for center in [256, 1024, 4096, 128^2]
        for offset in -10:10
            len = center + offset
            x = randn()
            arr = fill(x, len)
            @test maximum(arr) === x
            @test minimum(arr) === x
        end
    end
end

@test isnan(maximum([NaN]))
@test isnan(minimum([NaN]))
@test isequal(extrema([NaN]), (NaN, NaN))

@test isnan(maximum([NaN, 2.]))
@test isnan(maximum([2., NaN]))
@test isnan(minimum([NaN, 2.]))
@test isnan(minimum([2., NaN]))
@test isequal(extrema([NaN, 2.]), (NaN,NaN))

@test isnan(maximum([NaN, 2., 3.]))
@test isnan(minimum([NaN, 2., 3.]))
@test isequal(extrema([NaN, 2., 3.]), (NaN,NaN))

@test isnan(maximum([4., 3., NaN, 5., 2.]))
@test isnan(minimum([4., 3., NaN, 5., 2.]))
@test isequal(extrema([4., 3., NaN, 5., 2.]), (NaN,NaN))

 # test long arrays
@test isnan(maximum([NaN; 1.:10000.]))
@test isnan(maximum([1.:10000.; NaN]))
@test isnan(minimum([NaN; 1.:10000.]))
@test isnan(minimum([1.:10000.; NaN]))
@test isequal(extrema([1.:10000.; NaN]), (NaN,NaN))
@test isequal(extrema([NaN; 1.:10000.]), (NaN,NaN))

@test maximum(abs2, 3:7) == 49
@test minimum(abs2, 3:7) == 9
@test extrema(abs2, 3:7) == (9, 49)

@test maximum(Int16[1]) === Int16(1)
@test maximum(Vector(Int16(1):Int16(100))) === Int16(100)
@test maximum(Int32[1,2]) === Int32(2)

A = circshift(reshape(1:24,2,3,4), (0,1,1))
@test extrema(A,dims=1) == reshape([(23,24),(19,20),(21,22),(5,6),(1,2),(3,4),(11,12),(7,8),(9,10),(17,18),(13,14),(15,16)],1,3,4)
@test extrema(A,dims=2) == reshape([(19,23),(20,24),(1,5),(2,6),(7,11),(8,12),(13,17),(14,18)],2,1,4)
@test extrema(A,dims=3) == reshape([(5,23),(6,24),(1,19),(2,20),(3,21),(4,22)],2,3,1)
@test extrema(A,dims=(1,2)) == reshape([(19,24),(1,6),(7,12),(13,18)],1,1,4)
@test extrema(A,dims=(1,3)) == reshape([(5,24),(1,20),(3,22)],1,3,1)
@test extrema(A,dims=(2,3)) == reshape([(1,23),(2,24)],2,1,1)
@test extrema(A,dims=(1,2,3)) == reshape([(1,24)],1,1,1)
@test size(extrema(A,dims=1)) == size(maximum(A,dims=1))
@test size(extrema(A,dims=(1,2))) == size(maximum(A,dims=(1,2)))
@test size(extrema(A,dims=(1,2,3))) == size(maximum(A,dims=(1,2,3)))
@test extrema(x->div(x, 2), A, dims=(2,3)) == reshape([(0,11),(1,12)],2,1,1)

# any & all

@test @inferred any([]) == false
@test @inferred any(Bool[]) == false
@test @inferred any([true]) == true
@test @inferred any([false, false]) == false
@test @inferred any([false, true]) == true
@test @inferred any([true, false]) == true
@test @inferred any([true, true]) == true
@test @inferred any([true, true, true]) == true
@test @inferred any([true, false, true]) == true
@test @inferred any([false, false, false]) == false

@test @inferred all([]) == true
@test @inferred all(Bool[]) == true
@test @inferred all([true]) == true
@test @inferred all([false, false]) == false
@test @inferred all([false, true]) == false
@test @inferred all([true, false]) == false
@test @inferred all([true, true]) == true
@test @inferred all([true, true, true]) == true
@test @inferred all([true, false, true]) == false
@test @inferred all([false, false, false]) == false

@test @inferred any(x->x>0, []) == false
@test @inferred any(x->x>0, Int[]) == false
@test @inferred any(x->x>0, [-3]) == false
@test @inferred any(x->x>0, [4]) == true
@test @inferred any(x->x>0, [-3, 4, 5]) == true

@test @inferred all(x->x>0, []) == true
@test @inferred all(x->x>0, Int[]) == true
@test @inferred all(x->x>0, [-3]) == false
@test @inferred all(x->x>0, [4]) == true
@test @inferred all(x->x>0, [-3, 4, 5]) == false

@test reduce((a, b) -> a .| b, fill(trues(5), 24))  == trues(5)
@test reduce((a, b) -> a .| b, fill(falses(5), 24)) == falses(5)
@test reduce((a, b) -> a .& b, fill(trues(5), 24))  == trues(5)
@test reduce((a, b) -> a .& b, fill(falses(5), 24)) == falses(5)

@test_throws TypeError any(x->0, [false])
@test_throws TypeError all(x->0, [false])

# short-circuiting any and all

let c = [0, 0], A = 1:1000
    any(x->(c[1]=x; x==10), A)
    all(x->(c[2]=x; x!=10), A)

    @test c == [10,10]
end

# 19151 - always short circuit
let c = Int[], d = Int[], A = 1:9
    all((push!(c, x); x < 5) for x in A)
    @test c == 1:5

    any((push!(d, x); x > 4) for x in A)
    @test d == 1:5
end

# any/all with non-boolean collections

let f(x) = x == 1 ? true : x == 2 ? false : 1
    @test any(Any[false,true,false])
    @test @inferred any(map(f, [2,1,2]))
    @test @inferred any([f(x) for x in [2,1,2]])

    @test all(Any[true,true,true])
    @test @inferred all(map(f, [1,1,1]))
    @test @inferred all([f(x) for x in [1,1,1]])

    @test_throws TypeError any([1,true])
    @test_throws TypeError all([true,1])
    @test_throws TypeError any(map(f,[3,1]))
    @test_throws TypeError all(map(f,[1,3]))
end

# any and all with functors

struct SomeFunctor end
(::SomeFunctor)(x) = true

@test @inferred any(SomeFunctor(), 1:10)
@test @inferred all(SomeFunctor(), 1:10)

# in

@test in(1, Int[]) == false
@test in(1, Int[1]) == true
@test in(1, Int[2]) == false
@test in(0, 1:3) == false
@test in(1, 1:3) == true
@test in(2, 1:3) == true

# occursin

@test occursin("fox", "quick fox") == true
@test occursin("lazy dog", "quick fox") == false

# count

@test count(x->x>0, Int[]) == count(Bool[]) == 0
@test count(x->x>0, -3:5) == count((-3:5) .> 0) == 5
@test count([true, true, false, true]) == count(BitVector([true, true, false, true])) == 3
@test_throws TypeError count(sqrt, [1])
@test_throws TypeError count([1])
let itr = (x for x in 1:10 if x < 7)
    @test count(iseven, itr) == 3
    @test_throws TypeError count(itr)
    @test_throws TypeError count(sqrt, itr)
end
@test count(iseven(x) for x in 1:10 if x < 7) == 3
@test count(iseven(x) for x in 1:10 if x < -7) == 0

@test count(!iszero, Int[]) == 0
@test count(!iszero, Int[0]) == 0
@test count(!iszero, Int[1]) == 1
@test count(!iszero, [1, 0, 2, 0, 3, 0, 4]) == 4

## cumsum, cummin, cummax

z = rand(10^6)
let es = sum(BigFloat.(z)), es2 = sum(BigFloat.(z[1:10^5]))
    @test (es - sum(z)) < es * 1e-13
    cs = cumsum(z)
    @test (es - cs[end]) < es * 1e-13
    @test (es2 - cs[10^5]) < es2 * 1e-13
end

@test sum(Vector(map(UInt8,0:255))) == 32640
@test sum(Vector(map(UInt8,254:255))) == 509

A = reshape(map(UInt8, 101:109), (3,3))
@test @inferred(sum(A)) == 945
@test @inferred(sum(view(A, 1:3, 1:3))) == 945

A = reshape(map(UInt8, 1:100), (10,10))
@test @inferred(sum(A)) == 5050
@test @inferred(sum(view(A, 1:10, 1:10))) == 5050

# issue #11618
@test sum([-0.0]) === -0.0
@test sum([-0.0, -0.0]) === -0.0
@test prod([-0.0, -0.0]) === 0.0

# containment
let A = Vector(1:10)
    @test A ∋ 5
    @test A ∌ 11
    @test any(y->y==6,A)
end

# issue #18695
test18695(r) = sum( t^2 for t in r )
@test @inferred(test18695([1.0,2.0,3.0,4.0])) == 30.0
@test_throws ArgumentError test18695(Any[])

# issue #21107
@test foldr(-,2:2) == 2

# test neutral element not picked incorrectly for &, |
@test @inferred(foldl(&, Int[1])) === 1
@test_throws ArgumentError foldl(&, Int[])

# prod on Chars
@test prod(Char[]) == ""
@test prod(Char['a']) == "a"
@test prod(Char['a','b']) == "ab"

@testset "optimized reduce(vcat/hcat, A) for arrays" begin
    for args in ([1:2], [[1, 2]], [1:2, 3:4], [[3, 4, 5], 1:2], [1:2, [3.5, 4.5]],
                 [[1 2], [3 4; 5 6]], [reshape([1, 2], 2, 1), 3:4])
        X = reduce(vcat, args)
        Y = vcat(args...)
        @test X == Y
        @test typeof(X) === typeof(Y)
    end
    for args in ([1:2], [[1, 2]], [1:2, 3:4], [[3, 4, 5], 1:3], [1:2, [3.5, 4.5]],
                 [[1 2; 3 4], [5 6; 7 8]], [1:2, [5 6; 7 8]], [[5 6; 7 8], [1, 2]])
        X = reduce(hcat, args)
        Y = hcat(args...)
        @test X == Y
        @test typeof(X) === typeof(Y)
    end
end

# offset axes
i = Base.Slice(-3:3)
x = [j^2 for j in i]
@test sum(x) == sum(x.parent) == 28
i = Base.Slice(0:0)
x = [j+7 for j in i]
@test sum(x) == 7

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
arr = [1, 2, 3, 4];

# length of array
println(length(arr));

# Dimensions in array

println(ndims(arr));
mat = [1 2; 3 4];
println(ndims(mat));

# Size of each Dimensions
println(size(mat));

# index range of each dimension

println(axes(mat));

# Array access methods

a1 = [1, 2, 3, 4];

# Push - append element to end of array
push!(a1, 5); # ! indicates that the object passed will be modified, no copy made
println(a1);

# Push First - append element to start of array
pushfirst!(a1, 0);
println(a1);

# Pop - remove last element from array
pop!(a1);
println(a1);

# Pop First - remove first element from array
popfirst!(a1);
println(a1);

# Insert at specified index
insert!(a1, 3, 6);
println(a1);

# Delete at specified index
deleteat!(a1, 3);
println(a1);

# Splice - like delete, but returns deleted element
splice!(a1, 3);
println(a1);

# Splice - can also be used as replacement
splice!(a1, 3, 6);
println(a1);

# Append - inserts at end of array
append!(a1, 5);
println(a1);

# Prepend - inserts at start of array
prepend!(a1, 0);
println(a1);

# Array operations

a = Int64[1, 2, 3, 4];
b = Int64[1, 2, 1, 1];

# Matrix operations
println(a + b);
println(a - b);
println(a / b);

c = [1 2 3 4];
println(a * c);

println(a == b);
println(a != b);

# Broadcasting

println(broadcast(*, a, b));
println(a .* b);
println(a .== b);

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using ApproxFun,Base.Test

import ApproxFun:interlace,Multiplication

d=Interval()
D=Derivative(d);
B=ldirichlet();
Bn=lneumann();
A=[B 0;
   0 B;
   D-I 2I;
   0I D+I];

f=Fun(x->[exp(x),cos(x)],d)

b=Any[0.,0.,f...]
f1,f2=vec(f)
u=A\b
u1=vec(u)[1];u2=vec(u)[2];

@test norm(u1'-u1+2u2-f1)<10eps()
@test norm(u2'+u2-f2)<10eps()

Ai=Operator(A)
u=Ai\b
u1=vec(u)[1];u2=vec(u)[2];

@test norm(u1'-u1+2u2-f1)<10eps()
@test norm(u2'+u2-f2)<10eps()

A=[B 0;
   Bn 0;
   0 B;
   D^2-I 2I;
   0 D+I];

b=Any[0.,0.,0.,f...]

u=A\b
u1=vec(u)[1];u2=vec(u)[2];

@test norm(differentiate(u1,2)-u1+2u2-f1)<2eps()
@test norm(u2'+u2-f2)<2eps()

Ai=Operator(A)
u=Ai\b
u1=vec(u)[1];u2=vec(u)[2];

@test norm(u1''-u1+2u2-f1)<2eps()
@test norm(u2'+u2-f2)<2eps()

## Matrix exponential

n=4
d=fill(Interval(0.,1.),n)
B=Evaluation(d,0.)
D=Derivative(d)
A=rand(n,n)
L=[B;D-A]
u=L\eye(n)
@test norm(evaluate(u,1.)-expm(A))<eps(1000.)

n=4
d=fill(Interval(0.,1.),n)
B=Evaluation(d,0.)
D=Derivative(d)
A=rand(n,n)
L=[B;D-A]
u=L\eye(2)
@test norm(evaluate(u,1.)-expm(A)[:,1:2])<eps(1000.)

## Multiplication

d = Interval()
t=Fun(identity,d)
f = devec([t^2, sin(t)])
@test norm(((Derivative(space(f))*f)-Fun(t->[2t,cos(t)])).coefficients)<100eps()
@test norm((([1 2;3 4]*f)-Fun(t->[t^2+2sin(t),3t^2+4sin(t)])).coefficients)<100eps()

## Multiplication operator

Γ=Circle()∪Circle(0.5)
G=Fun(z->in(z,Γ[2])?[1 -z^(-1); 0 1]:
                   [z 0; 0 z^(-1)],Γ);
G1=demat(mat(G)[:,1])
M=Multiplication(G,space(G1))
u=M*G1
@test norm(u(exp(.1im))-[exp(.2im),0])<100eps()
@test norm(u(.5exp(.1im))-[1,0])<100eps()

# Vector operations
@test_approx_eq (Fun(x->[1., 2.]) + [2, 2])(0.) [3., 4.]

G=Fun(z->[-1 -3; -3 -1]/z +
         [ 2  2;  1 -3] +
         [ 2 -1;  1  2]*z,Circle())

@test G[1,1](exp(0.1im)) == G(exp(0.1im))[1,1]

F̃ = (G-I)[:,1]
F=Fun((G-I)[:,1])

@test_approx_eq F(exp(0.1im)) [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]
@test_approx_eq Fun(F̃,space(F))(exp(0.1im)) [-exp(-0.1im)+1+2exp(0.1im);-3exp(-0.1im)+1+1exp(0.1im)]

@test coefficients(F̃,space(F)) == F.coefficients
@test Fun(F̃,space(F)) == F

## Check conversion

f=Fun(t->[cos(t) 0;sin(t) 1],[-π,π])
g=Fun(f,Space(PeriodicInterval(-π,π)))
@test_approx_eq g(.1) f(.1)

## Operator * Matrix

D=Derivative()

u=D*[Fun(exp) Fun(cos)]
@test_approx_eq u(0.1) [exp(0.1) -sin(0.1)]


text:
"""
Types deriving from `AbstractSegment{T}` must have a `first` and `last` function each returning a value of type `T`,
and `first(i) <= last(i)` must always be true.
"""
abstract type AbstractSegment{T} end

Base.first(i::AbstractSegment{T}) where T = i.first
Base.last(i::AbstractSegment{T}) where T = i.last

function basic_isless(u::AbstractSegment, v::AbstractSegment)
    return first(u) < first(v) || (first(u) == first(v) && last(u) < last(v))
end

function Base.isless(u::AbstractSegment, v::AbstractSegment)
    return basic_isless(u, v)
end

"""
A basic Segment.
"""
struct Segment{T} <: AbstractSegment{T}
    first::T
    last::T
end

Base.convert(::Type{Segment{T}}, range::AbstractRange{T}) where T =
    Segment(first(range), last(range))
Segment(range::AbstractRange{T}) where T = convert(Segment{T}, range)

abstract type AbstractSegmentTree{T} end

"""
A basic Segment tree
"""
struct SegmentTree{T} <: AbstractSegmentTree{T}
    dataarray::Vector{Int}
    segments::Vector{Segment{T}}
end

"""
An Segment with some associated data.
"""
struct SegmentValue{K,V} <: AbstractSegment{K}
    first::K
    last::K
    value::V
end

SegmentValue(range::AbstractRange{K}, value::V) where {K,V} =
    SegmentValue(first(range), last(range), value)

valtype(::Type{SegmentValue{K,V}}) where {K,V} = V

value(i::SegmentValue{K,V}) where {K,V} = i.value

Base.print(io::IO, x::Segment) = print(io, "\n($(first(x)),$(last(x)))")
function Base.show(io::IO, x::Segment)
    show(io, typeof(x))
    print(io, x)
end

Base.print(io::IO, x::SegmentValue) = print(io, "\n($(first(x)),$(last(x))) => $(value(x))")
function Base.show(io::IO, x::SegmentValue)
    show(io, typeof(x))
    print(io, x)
end

# Each of these types is indexes by K, V, B, where
#   K : Segment type. Segments are represented as (K, K) tuples.
#   V : Segment value type. This an Segment type that may have associated
#        data. It must have `first`, `last`, and `isless` methods.
#   B : Integer giving the B-tree order.

abstract type Node{K,V,B} end

mutable struct InternalNode{K,V,B} <: Node{K,V,B}
    # Internal nodes are keyed by the minimum Segment in the right subtree.  We
    # need internal node keys to be Segments themselves, since ordering by
    # start value alone only works if start values are unique. We don't
    # force that restriction, so we must break ties in order to split nodes that
    # are full of Segments with the same start value.
    keys::Slice{Segment{K},B}

    # Maximum end ofd this node
    maxend::K

    # Maximum child subtree end-points
    maxends::Slice{K,B}

    children::Slice{Node{K,V,B},B}

    parent::Union{Nothing,InternalNode{K,V,B}}

    # Sibling/cousin pointers.
    left::Union{Nothing,InternalNode{K,V,B}}
    right::Union{Nothing,InternalNode{K,V,B}}

    function InternalNode{K,V,B}() where {K,V,B}
        return new{K,V,B}(
            Slice{Segment{K},B}(), zero(K), Slice{K,B}(), Slice{Node{K,V,B},B}(),
            nothing,
            nothing,
            nothing)
    end
end

function minstart(t::InternalNode)
    return t.keys[1].first
end

mutable struct LeafNode{K,V,B} <: Node{K,V,B}
    entries::Vector{V}
    keys::Vector{Segment{K}}
    count::UInt32 # number of stored entries

    # maximum Segment end-point in this leaf node.
    maxend::K

    parent::Union{Nothing,InternalNode{K,V,B}}

    # Sibling/cousin pointers.
    left::Union{Nothing,LeafNode{K,V,B}}
    right::Union{Nothing,LeafNode{K,V,B}}

    function LeafNode{K,V,B}() where {K,V,B}
        return new{K,V,B}(
            Array{V}(undef, B), Array{Segment{K}}(undef, B), 0,
            zero(K),
            nothing,
            nothing,
            nothing)
    end
end

function minstart(t::LeafNode)
    return first(t.keys[1])
end

function Base.resize!(leaf::LeafNode{K,V,B}, n) where {K,V,B}
    @assert 0 <= n <= B
    leaf.count = n
end

function Base.splice!(leaf::LeafNode{K,V,B}, i::Integer) where {K,V,B}
    slice_splice!(leaf.keys, leaf.count, i)
    x, leaf.count = slice_splice!(leaf.entries, leaf.count, i)
    return x
end

function Base.insert!(leaf::LeafNode{K,V,B}, i::Integer, value::V) where {K,V,B}
    slice_insert!(leaf.keys, leaf.count, i, Segment{K}(first(value), last(value)))
    slice_insert!(leaf.entries, leaf.count, i, value)
    leaf.count += 1
end

function Base.push!(leaf::LeafNode{K,V,B}, value::V) where {K,V,B}
    leaf.count += 1
    leaf.entries[leaf.count] = value
    leaf.keys[leaf.count] = Segment{K}(first(value), last(value))
end

function Base.pop!(leaf::LeafNode{K,V,B}) where {K,V,B}
    x = leaf.entries[leaf.count]
    leaf.count -= 1
    return x
end

function Base.iterate(leaf::LeafNode, i::Int=1)
    return i ≤ length(leaf) ? (leaf.entries[i], i + 1) : nothing
end

mutable struct SegmentBTree{K,V,B}
    root::Node{K,V,B}
    n::Int # Number of entries

    function SegmentBTree{K,V,B}() where {K,V,B}
        return new{K,V,B}(LeafNode{K,V,B}(), 0)
    end

    # Construct an Segment tree from a sorted array of Segments.
    #
    # This is generally much more efficient than constructing the same tree by
    # inserting Segments one by one.
    #
    # Args:
    #   entries: Segment entry values in sorted order.
    #
    function SegmentBTree{K,V,B}(entries::AbstractVector{V}) where {K,V,B}
        if !issorted(entries, lt=basic_isless)
            error("Segments must be sorted to construct an SegmentTree")
        end

        # Here's the plan: figure out how many leaf nodes we need, allocate
        # them all up front. Copy in the keys and values, then work up towards
        # the root.

        n = length(entries)

        if n == 0
            return new{K,V,B}(LeafNode{K,V,B}(), 0)
        end

        numleaves = cld(n, B - 2)
        leaves = [LeafNode{K,V,B}() for _ in 1:numleaves]

        maxends = Vector{K}(undef, numleaves)
        minkeys = Vector{Segment{K}}(undef, numleaves)

        # divy up the keys and values among the leaves
        keys_per_leaf = cld(n, numleaves)
        for i in 1:numleaves
            u = (i - 1) * keys_per_leaf + 1
            v = min(n, i * keys_per_leaf)
            minkeys[i] = Segment{K}(first(entries[u]), last(entries[u]))
            maxends[i] = last(entries[u])
            for j in u:v
                push!(leaves[i], entries[j])
                maxends[i] = max(maxends[i], last(entries[j]))
            end
            leaves[i].maxend = maxends[i]
        end

        children = leaves
        while length(children) > 1
            # sibling pointers
            for i in 1:length(children) - 1
                children[i].right = children[i + 1]
                children[i + 1].left = children[i]
            end

            # make parents
            numparents = cld(length(children), B - 2)
            parents = [InternalNode{K,V,B}() for _ in 1:numparents]

            # divy up children among parents
            children_per_parent = cld(length(children), numparents)
            for i in 1:numparents
                u = (i - 1) * keys_per_leaf + 1
                v = min(length(children), i * keys_per_leaf)
                maxend = maxends[u]
                for j in u:v
                    push!(parents[i].children, children[j])
                    push!(parents[i].maxends, maxends[j])
                    children[j].parent = parents[i]
                    maxend = max(maxend, maxends[j])
                    if j > u
                        push!(parents[i].keys, minkeys[j])
                    end
                end
                minkeys[i] = minkeys[u]
                parents[i].maxend = maxend
                maxends[i] = maxend
            end

            children = parents
        end
        @assert length(children) == 1

        return new{K,V,B}(children[1], n)
    end
end

# Default B-tree order
const SegmentTree{K,V} = SegmentBTree{K,V,64}

# Show

function Base.show(io::IO, it::SegmentTree)
    show(io, typeof(it))
    n = length(it)
    if length(it) > 6
        # Hacky random access ...
        for (i, x) in enumerate(it)
            i < 4 && print(x)
        end
        print("\n\u22EE") # Vertical ellipsis
        for (i, x) in enumerate(it)
            i > (n - 3) && print(x)
        end
    else
        for x in it
            print(x)
end
    end
end

# Length
# ------

function Base.length(t::SegmentBTree)
    return t.n
end

function depth(t::SegmentBTree)
    k = 1
    node = t.root
    while !isa(node, LeafNode)
        k += 1
        node = node.children[1]
    end
    return k
end

function Base.length(t::InternalNode)
    return length(t.children)
end

function Base.length(t::LeafNode)
    return t.count
end

function Base.isempty(t::SegmentBTree)
    return t.n == 0
end

function Base.isempty(t::InternalNode)
    return isempty(t.children)
end

function Base.isempty(t::LeafNode)
    return t.count == 0
end


text:
"""
Calculate the Legendre polynomials. TODO reference
"""
function legendre_polynomials(  j::Int,
                                ε::AbstractMatrix,
                                ε⁻¹::AbstractMatrix,
                                mx::Int,
                                nx::Int,
                                G::Geometry)

    @unpack coslat_NH, sinlat_NH = G

    alp = zeros(mx+1,nx)

    x = sinlat_NH[j]
    y = coslat_NH[j]

    # Start recursion with n = 1 (m = l) diagonal
    alp[1,1] = sqrt(0.5)
    for m in 2:mx+1
        alp[m,1] = sqrt(0.5*(2m - 1)/(m-1))*y*alp[m-1,1]
    end

    # Continue with other elements
    for m in 1:mx+1
        alp[m,2] = (x*alp[m,1])*ε⁻¹[m,2]
    end

    for n in 3:nx
        for m in 1:mx+1
            alp[m,n] = (x*alp[m,n-1] - ε[m,n-1]*alp[m,n-2])*ε⁻¹[m,n]
        end
    end

    # Zero polynomials with absolute values smaller than 10^(-30)
    # small = 1e-30
    # for n in 1:nx
    #     for m in 1:mx+1
    #         if abs(alp[m,n]) <= small
    #             alp[m,n] = 0.0
    #         end
    #     end
    # end

    # pick off the required polynomials
    return alp[1:mx,1:nx]
end

"""
Computes the inverse Legendre transform.
"""
function legendre_inverse(  input::Array{Complex{NF},2},
                            G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    @unpack leg_weight, nsh2, leg_poly = G.spectral
    @unpack trunc, mx, nx = G.spectral
    @unpack nlat, nlat_half = G.geometry

    # Initialize output array
    output = zeros(Complex{NF}, mx, nlat)

    # Loop over Northern Hemisphere, computing odd and even decomposition of incoming field
    for j in 1:nlat_half
        j1 = nlat + 1 - j

        # Initialise arrays TODO preallocate them in a separate struct
        even = zeros(Complex{NF}, mx)
        odd  = zeros(Complex{NF}, mx)

        # Compute even decomposition
        for n in 1:2:nx
            for m in 1:nsh2[n]
                even[m] = even[m] + input[m,n]*leg_poly[m,n,j]
            end
        end

        # Compute odd decomposition
        for n in 2:2:nx
            for m in 1:nsh2[n]
                odd[m] = odd[m] + input[m,n]*leg_poly[m,n,j]
            end
        end

        # Compute Southern Hemisphere
        output[:,j1] = even + odd

        # Compute Northern Hemisphere
        output[:,j]  = even - odd
    end
    return output
end

"""
Computes the Legendre transform
"""
function legendre(  input::Array{Complex{NF},2},
                    G::GeoSpectral{NF}) where {NF<:AbstractFloat}

    @unpack leg_weight, nsh2, leg_poly = G.spectral
    @unpack trunc, mx, nx = G.spectral
    @unpack nlat, nlat_half = G.geometry

    # Initialise output array
    output = zeros(Complex{NF}, mx, nx)

    even = zeros(Complex{NF}, mx, nlat_half)
    odd  = zeros(Complex{NF}, mx, nlat_half)

    # Loop over Northern Hemisphere, computing odd and even decomposition of
    # incoming field. The Legendre weights (leg_weight) are applied here
    for j in 1:nlat_half
        # Corresponding Southern Hemisphere latitude
        j1 = nlat + 1 - j

        even[:,j] = (input[:,j1] + input[:,j])*leg_weight[j]
        odd[:,j]  = (input[:,j1] - input[:,j])*leg_weight[j]
    end

    # The parity of an associated Legendre polynomial is the same
    # as the parity of n' - m'. n', m' are the actual total wavenumber and zonal
    # wavenumber, n and m are the indices used for SPEEDY's spectral packing.
    # m' = m - 1 and n' = m + n - 2, therefore n' - m' = n - 1

    # Loop over coefficients corresponding to even associated Legendre polynomials
    for n in 1:2:trunc+1
        for m in 1:nsh2[n]
            output[m,n] = dot(leg_poly[m,n,:], even[m,:])
        end
    end

    # Loop over coefficients corresponding to odd associated Legendre polynomials
    for n in 2:2:trunc+1
        for m in 1:nsh2[n]
            output[m,n] = dot(leg_poly[m,n,:], odd[m,:])
        end
    end

    return output
end


text:
# Quadratic decision function #
# ========================== #
struct QuadraticDecisionObjectiveBridge{T} <: MOIB.Objective.AbstractBridge
    decision_function::QuadraticDecisionFunction{T}
end

function MOIB.Objective.bridge_objective(::Type{QuadraticDecisionObjectiveBridge{T}}, model::MOI.ModelLike,
                                         f::QuadraticDecisionFunction{T}) where T
    # All decisions have been mapped to the variable part terms
    # at this point.
    F = MOI.ScalarQuadraticFunction{T}
    # Set the bridged objective
    MOI.set(model, MOI.ObjectiveFunction{F}(),
            MOI.ScalarQuadraticFunction(
                f.variable_part.affine_terms,
                f.variable_part.quadratic_terms,
                zero(T)))
    # Save decision function to allow modifications
    return QuadraticDecisionObjectiveBridge{T}(f)
end

function MOIB.Objective.supports_objective_function(
    ::Type{<:QuadraticDecisionObjectiveBridge}, ::Type{<:QuadraticDecisionFunction})
    return true
end
MOIB.added_constrained_variable_types(::Type{<:QuadraticDecisionObjectiveBridge}) = Tuple{DataType}[]
function MOIB.added_constraint_types(::Type{<:QuadraticDecisionObjectiveBridge})
    return Tuple{DataType, DataType}[]
end
function MOIB.Objective.concrete_bridge_type(::Type{<:QuadraticDecisionObjectiveBridge{T}},
                                             ::Type{QuadraticDecisionFunction{T}}) where T
    return QuadraticDecisionObjectiveBridge{T}
end
function MOIB.set_objective_function_type(::Type{QuadraticDecisionObjectiveBridge{T}}) where T
    return MOI.ScalarQuadraticFunction{T}
end

function MOI.get(::QuadraticDecisionObjectiveBridge, ::MOI.NumberOfVariables)
    return 0
end

function MOI.get(::QuadraticDecisionObjectiveBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

function MOI.delete(::MOI.ModelLike, ::QuadraticDecisionObjectiveBridge)
    # Nothing to delete
    return nothing
end

function MOI.set(::MOI.ModelLike, ::MOI.ObjectiveSense,
                 ::QuadraticDecisionObjectiveBridge, ::MOI.OptimizationSense)
    # Nothing to handle if sense changes
    return nothing
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOIB.ObjectiveFunctionValue{QuadraticDecisionFunction{T}},
                 bridge::QuadraticDecisionObjectiveBridge{T}) where T
    f = bridge.decision_function
    G = MOI.ScalarQuadraticFunction{T}
    obj_val = MOI.get(model, MOIB.ObjectiveFunctionValue{G}(attr.result_index))
    # Calculate and add constant
    constant = f.variable_part.constant +
        f.decision_part.constant
    return obj_val + f.variable_part.constant
end

function MOI.get(model::MOI.ModelLike,
                 attr::MOI.ObjectiveFunction{QuadraticDecisionFunction{T}},
                 bridge::QuadraticDecisionObjectiveBridge{T}) where T
    f = bridge.decision_function
    # Remove mapped variables to properly unbridge
    from_decision(v) = begin
        result = any(t -> t.variable_index == v,
                     f.decision_part.affine_terms)
        result |= any(t -> t.variable_index_1 == v || t.variable_index_2 == v,
                      f.decision_part.quadratic_terms)
        result |= any(t -> t.variable_index_1 == v,
                      f.cross_terms.affine_terms)
    end
    g = QuadraticDecisionFunction(
        MOIU.filter_variables(v -> !from_decision(v), f.variable_part),
        copy(f.decision_part),
        copy(f.cross_terms))
    return g
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::MOI.ScalarConstantChange) where T
    f = bridge.decision_function
    f = f.linear_quadratic_terms
    # Modify constant of variable part
    lq.variable_part.constant = change.new_constant
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::MOI.ScalarCoefficientChange) where T
    f = bridge.decision_function
    # Update variable part
    modify_coefficient!(f.variable_part.affine_terms, change.variable, change.new_coefficient)
    # Modify decision part if variable originated from a decision
    i = something(findfirst(t -> t.variable_index == change.variable,
                            f.decision_part.affine_terms), 0)
    if !iszero(i)
        modify_coefficient!(f.decision_part.affine_terms, change.variable, change.new_coefficient)
    end
    return nothing
end

function MOI.modify(model::MOI.ModelLike, bridge::QuadraticDecisionObjectiveBridge{T}, change::DecisionCoefficientChange) where T
    f = bridge.decision_function
    variable_terms = f.variable_part.affine_terms
    decision_terms = f.decision_part.affine_terms
    # Check if decision has been mapped
    i = something(findfirst(t -> t.variable_index == change.decision,
                            variable_terms), 0)
    if !iszero(i)
        # Update mapped variable through ScalarCoefficientChange
        F = QuadraticDecisionFunction{T}
        MOI.modify(model, MOI.ObjectiveFunction{F}(),
                   MOI.ScalarCoefficientChange(change.decision, change.new_coefficient))
    else
        i = something(findfirst(t -> t.variable_index == change.decision,
                                decision_terms), 0)
        # Query the fixed decision value
        unbridged = MOIB.unbridged_variable_function(model, change.decision)::MOI.SingleVariable
        decision_value = MOI.get(model, MOI.VariablePrimal(), unbridged.variable)
        # Update decision part of objective constant
        coefficient = iszero(i) ? zero(T) : decision_terms[i].coefficient
        lq.decision_part.constant +=
            (change.new_coefficient - coefficient) * decision_value
    end
    # Update coefficient in decision part
    modify_coefficient!(decision_terms, change.decision, change.new_coefficient)
    return nothing
end