text:
export total_xray_fraction,
    radial_uv_fraction,
    disk_uv_fraction,
    spectral_band_fraction_frequency,
    spectral_band_fraction

"""
    spectral_band_fraction_frequency(bb, low, high)

Computes the amount of black body radiance that is emitted in
a particular frequency band.
"""
function spectral_band_fraction_frequency(bb::BlackBody, low, high)
    total_radiance = radiance(bb)
    band_radiance = spectral_band_radiance_frequency(bb, low, high)
    return band_radiance / total_radiance
end
"""
    spectral_band_fraction(bb, low, high)
Computes the amount of black body radiance that is emitted in
a particular energy band.
"""

function spectral_band_fraction(bb::BlackBody, low, high)
    total_radiance = radiance(bb)
    band_radiance = spectral_band_radiance(bb, low, high)
    return band_radiance / total_radiance
end

"""
    disk_uv_fraction(bh, r; uv_low_kev = 0.00387, uv_high_kev = 0.06)
Computes the fraction of luminosity in the UV band, for a ring at disk r.
The UV band is defined as uv_low_kev -- uv_high_kev.
"""
function disk_uv_fraction(bh, r; uv_low_kev = 0.00387, uv_high_kev = 0.06)
    if r <= bh.isco
        return 0.0
    end
    temperature = compute_disk_temperature(bh, r)
    bb = BlackBody(temperature)
    return spectral_band_fraction(bb, uv_low_kev, uv_high_kev)
end

"""
    warm_uv_fraction(warm, bh, r; uv_low_kev = 0.00387, uv_high_kev = 0.06)
Computes the fraction of luminosity in the UV band, for a ring at disk r.
The UV band is defined as uv_low_kev -- uv_high_kev.
"""
function compute_warm_uv_fraction(
    corona,
    warm,
    bh;
    uv_low_kev = 0.00387,
    uv_high_kev = 0.06,
    energy_low = 1e-5,
    energy_high = 1e3,
    n_energy = 5000,
)
    energy_range = 10 .^ range(log10(energy_low), log10(energy_high), length = n_energy)
    warm_spectral_luminosity = compute_warm_spectral_luminosity(corona, warm, energy_range)
    uv_lower_b = searchsortedfirst(energy_range, uv_low_kev)
    uv_upper_b = searchsortedfirst(energy_range, uv_high_kev)
    delta_energy_range = diff(energy_range)
    uv_lumin = sum(
        warm_spectral_luminosity[uv_lower_b:uv_upper_b] .*
        delta_energy_range[uv_lower_b:uv_upper_b],
    )
    total_lumin = sum(warm_spectral_luminosity[2:end] .* delta_energy_range)
    return uv_lumin / total_lumin
end

"""
    disk_xray_fraction(bh::BlackHole, r; xray_low_kev = 0.00387, xray_high_kev = 0.06)
Computes the fraction of luminosity in the UV band, for a ring at disk r.
The UV band is defined as uv_low_kev -- uv_high_kev.
"""
function disk_xray_fraction(bh::BlackHole; xray_low_kev = 0.06, xray_high_kev = 1)
    total_lumin = compute_bolometric_luminosity(bh)
    xray_lumin = disk_spectral_band_radiance(bh, xray_low_kev, xray_high_kev)
    xray_lumin / total_lumin
end

"""
    total_xray_fraction(corona, warm; xray_low_kev = 0.00387, xray_high_kev = 0.06)
Computes the fraction of luminosity that is in the X-ray fraction.
"""
function total_xray_fraction(
    corona,
    warm;
    energy_low = 1e-5,
    energy_high = 1e3,
    xray_low_kev = 0.06,
    xray_high_kev = 1,
    n_energy = 5000,
)
    energy_range = 10 .^ range(log10(energy_low), log10(energy_high), length = n_energy)
    cl = compute_corona_spectral_luminosity(corona, warm, energy_range)
    dl = Qsosed.compute_disk_spectral_luminosity(
        corona.bh,
        energy_range,
        r_min = warm.radius,
        r_max = gravity_radius(corona.bh),
    )
    wl = compute_warm_spectral_luminosity(corona, warm, energy_range)
    total_spectral_lumin = cl + dl + wl
    xray_lower_b = searchsortedfirst(energy_range, xray_low_kev)
    xray_upper_b = searchsortedfirst(energy_range, xray_high_kev)
    delta_energy_range = diff(energy_range)
    xray_lumin = sum(
        total_spectral_lumin[xray_lower_b:xray_upper_b] .*
        delta_energy_range[xray_lower_b:xray_upper_b],
    )
    total_lumin = compute_bolometric_luminosity(corona.bh)
    return xray_lumin / total_lumin
end
total_xray_fraction(
    model;
    energy_low = 1e-5,
    energy_high = 1e3,
    xray_low_kev = 0.06,
    xray_high_kev = 1,
    n_energy = 5000,
) = total_xray_fraction(
    model.corona,
    model.warm,
    energy_low = energy_low,
    energy_high = energy_high,
    xray_low_kev = xray_low_kev,
    xray_high_kev = xray_high_kev,
    n_energy = n_energy,
)

"""
    radial_uv_fraction(corona, warm; uv_low_kev = 0.00387, uv_high_kev = 0.06)
Fraction of power emitted in the UV band per radius.
"""
function radial_uv_fraction(
    corona,
    warm;
    r_min = "isco",
    r_max = "gravity_radius",
    uv_low_kev = 0.00387,
    uv_high_kev = 0.06,
    n_r = 5000,
)
    if r_min == "isco"
        r_min = corona.bh.isco
    end
    if r_max == "gravity_radius"
        r_max = gravity_radius(corona.bh)
    end
    r_range =
        10 .^ range(log10(r_min), log10(r_max), length = n_r)
    warm_uv_fraction = compute_warm_uv_fraction(
        corona,
        warm,
        corona.bh,
        uv_low_kev = uv_low_kev,
        uv_high_kev = uv_high_kev,
    )
    ret = Float64[]
    for r in r_range
        if r <= corona.radius
            push!(ret, 0)
        elseif corona.radius < r < warm.radius
            push!(ret, warm_uv_fraction)
        else
            push!(
                ret,
                disk_uv_fraction(
                    corona.bh,
                    r,
                    uv_low_kev = uv_low_kev,
                    uv_high_kev = uv_high_kev,
                ),
            )
        end
    end
    return r_range, ret
end
radial_uv_fraction(model; r_min = "isco", r_max="gravity_radius", n_r = 5000) = radial_uv_fraction(
    model.corona,
    model.warm,
    uv_low_kev = model.parameters.uv_min_energy,
    uv_high_kev = model.parameters.uv_max_energy,
    n_r = n_r,
    r_min=r_min,
    r_max=r_max
)


text:
export heisenberg, transverse_ising

"""
    heisenberg(nbit::Int; periodic::Bool=true)

1D heisenberg hamiltonian, for its ground state, refer `PRB 48, 6141`.
"""
function heisenberg(nbit::Int; periodic::Bool=true)
    map(1:(periodic ? nbit : nbit-1)) do i
        j=i%nbit+1
        repeat(nbit,X,(i,j)) + repeat(nbit, Y, (i,j)) + repeat(nbit, Z, (i,j))
    end |> sum
end

"""
    transverse_ising(nbit::Int; periodic::Bool=true)

1D transverse ising hamiltonian.
"""
function transverse_ising(nbit::Int; periodic::Bool=true)
    ising_term = map(1:(periodic ? nbit : nbit-1)) do i
        repeat(nbit,Z,(i,i%nbit+1))
    end |> sum
    ising_term + sum(map(i->put(nbit,i=>X), 1:nbit))
end

ext: jl

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
using SimpleSDMLayers

"""
    wombling(layer::T; convert_to::Type=Float64) where {T <: SimpleSDMLayer}

Performs a lattice wombling on a `SimpleSDMLayer`.
"""
function wombling(layer::T; convert_to::Type=Float64) where {T <: SimpleSDMLayer}
    try
        global nan = convert(convert_to, NaN)
    catch
        throw(ArgumentError("The type given as `convert_to` must have a `NaN` value."))
    end

    # Get the values for x and y
    y = collect(longitudes(layer))
    x = collect(latitudes(layer))

    # Get the grid
    z = convert(Matrix{Union{Nothing,convert_to}}, layer.grid)
    replace!(z, nothing => nan)
    return wombling(x, y, convert(Matrix{convert_to}, z))
end

# We want to make sure that the layers are returned without NaN values, and
# adding this method makes the code easier to write
Base.isnan(::Nothing) = false

function SimpleSDMLayers.SimpleSDMPredictor(W::T) where {T <: LatticeWomble}
    rate = SimpleSDMLayers.SimpleSDMPredictor(W.m, extrema(W.y)..., extrema(W.x)...)
    direction = SimpleSDMLayers.SimpleSDMPredictor(W.θ, extrema(W.y)..., extrema(W.x)...)
    rate.grid[findall(isnan, rate.grid)] .= nothing
    direction.grid[findall(isnan, direction.grid)] .= nothing
    return (rate, direction)
end

SimpleSDMLayers.SimpleSDMResponse(W::T) where {T <: LatticeWomble} = convert.(SimpleSDMResponse, SimpleSDMPredictor(W))


text:
function _bump(records::AbstractVector{R}, b::Int) where {T, R<:Record{T}}

    new_records = Vector{R}(undef, length(records))

    for (i, record) in enumerate(records)
        new_record  = Record{T}(record.chrom, record.first + b, record.last + b, record.value)
        new_records[i] = new_record
    end

    return new_records
end

_bump_forward(records::AbstractVector{<:Record}) = _bump(records, 1)
_bump_back(records::AbstractVector{<:Record}) = _bump(records, -1)

function _range(record::Record; right_open=true)

    pos_start = right_open ? record.first : record.first + 1
    pos_end = right_open ? record.last - 1 : record.last

    return pos_start : pos_end
end

function _range(records::AbstractVector{<:Record}; right_open=true)

    pos_start = _range(records[1], right_open=right_open)[1]
    pos_end = _range(records[end], right_open=right_open)[end]

    return  pos_start : pos_end
end

function compress(chroms::AbstractVector{<:AbstractString}, n::AbstractVector{Int}, values::AbstractVector{T}; right_open = true, bump_back=true) where {T<:Real}

    ranges = Vector{UnitRange{Int}}()
    compressed_values = Vector{T}()
    compressed_chroms = Vector{String}()

    range_start = 1
    push!(compressed_values, values[1])

    for (index, value) in enumerate(values)
        if value != compressed_values[end]
            push!(ranges, n[range_start] : n[index - 1] )
            push!(compressed_values, value)
            push!(compressed_chroms, chroms[index])
            range_start = index
        end

        if index == length(values)
            push!(ranges, n[range_start] : n[index] )
            push!(compressed_values, value)
            push!(compressed_chroms, chroms[index])
        end
    end

    if right_open
        for (index, value) in enumerate(ranges)
            ranges[index] = first(value) : last(value) + 1
        end
    else
        for (index, value) in enumerate(ranges)
            ranges[index] = first(value) -1 : last(value)
        end
    end

    len = length(ranges)

    new_records = Vector{Record{T}}(undef, len)

    for (index, chrom, range, value) in zip(1:len, compressed_chroms, ranges, compressed_values)
        new_records[index]  = Record{T}(chrom, first(range), last(range), value)
    end

    return bump_back ? _bump_back(new_records) : new_records

end

compress(chrom::AbstractString, n::AbstractVector{Int}, values::AbstractVector{<:Real}; right_open = true, bump_back=true) = compress(fill(chrom, length(n)), n, values, right_open = right_open, bump_back = bump_back)

function expand(records::AbstractVector{R}; right_open=true, bump_forward=true) where {T, R<:Record{T}}

    #TODO: ensure records are sorted with no overlap.

    if bump_forward
        records =  _bump_forward(records)
    end

    total_range =_range(records, right_open = right_open)

    values = Vector{T}(undef, length(total_range))
    chroms = Vector{String}(undef, length(total_range))

    for record in records
        values[indexin(_range(record, right_open = right_open), total_range)] .= record.value
        chroms[indexin(_range(record, right_open = right_open), total_range)] .= record.chrom
    end

    return collect(total_range), values, chroms
end

expand(chrom::AbstractString, firsts::AbstractVector{Int}, lasts::AbstractVector{Int}, values::AbstractVector{<:Real}; right_open=true, bump_forward=true) = expand(fill(chrom, length(firsts)), firsts, lasts, values, right_open=right_open, bump_forward=bump_forward)
expand(chroms::AbstractVector{<:AbstractString}, firsts::AbstractVector{Int}, lasts::AbstractVector{Int}, values::Vector{<:Real}; right_open=true, bump_forward=true) = expand(Record.(chroms, firsts, lasts, values), right_open=right_open, bump_forward=bump_forward)


text:
struct EnumFilterHierarchicalBatchQs <: MultiObjectiveBatchQs
    model::SVDD.OCClassifier
    inf_measure::SequentialPQs
    rep_measure::Function
    div_measure::Function
    k::Int

    function EnumFilterHierarchicalBatchQs(model::SVDD.OCClassifier, informativeness::SequentialPQs; representativeness::Symbol, diversity::Symbol,
        k::Int)::EnumFilterHierarchicalBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))

        representativeness_measure = get_rep_measure(representativeness)
        diversity_measure = get_div_measure(diversity)
        return new(model, informativeness, representativeness_measure, diversity_measure, k)
    end
end

"""
Idea: Find 4*k elements with largest representativeness,
    select 2*k elements with largest informativeness,
    select batch of size k with maximum diversity
"""
function select_batch(qs::EnumFilterHierarchicalBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end

    # representativeness
    rep_scores = qs.rep_measure(qs.model, x, labels, candidate_indices)
    most_representative_indices = candidate_indices[sortperm(rep_scores, rev=true)[1:min(4*qs.k, num_observations)]]

    #informativeness
    inf_scores = qs_score(qs.inf_measure, x, labels)[most_representative_indices]
    batch_candidate_indices = most_representative_indices[sortperm(inf_scores, rev=true)[1:min(2*qs.k, num_observations)]]

    best_batch = nothing
    best_score = -Inf
    for batch in subsets(batch_candidate_indices, qs.k)
        batch_diversity = qs.div_measure(qs.model, x, batch)
        if (batch_diversity > best_score)
            best_batch = batch
            best_score = batch_diversity
        end
    end
    return best_batch
end


text:
const def_exp_tim = 15

struct ObservationPlan{T1<:Integer, T2<:AbstractFloat}
    N_obs::T1
    N_nights::T1
    obs_per_night::T1
    time_per_obs::T2
    dead_time::T2
end

function ObservationPlan(;N_obs::T1=0, obs_per_night::T1=0, time_per_obs::T2=0.0,
                         dead_time::T2=0.0) where {T1<:Integer, T2<:AF}
    # assertions
    @assert N_obs > 0
    @assert obs_per_night > 0
    @assert time_per_obs >= def_exp_tim
    @assert dead_time >= 0.0

    # TODO only works if exposure times are 0 or factor of 15
    @assert iszero(time_per_obs % 15)
    @assert iszero(dead_time % 15)

    # get number of nights
    N_nights = Int(N_obs ÷ obs_per_night) + one(T1) * !iszero(N_obs % obs_per_night)
    return ObservationPlan(N_obs, N_nights, obs_per_night, time_per_obs, dead_time)
end

function get_nt(obs::ObservationPlan)
    nt_per_exposure = Int(obs.time_per_obs ÷ def_exp_tim) + 1 * !iszero(obs.time_per_obs % def_exp_tim)
    nt_per_deadtime = Int(obs.dead_time ÷ def_exp_tim) + 1 * !iszero(obs.dead_time % def_exp_tim)
    return nt_per_exposure, nt_per_deadtime
end

function get_skip_times(nt_per_night, nt_per_exposure, nt_per_deadtime)
    # get relevant numbers for iteration
    nt_per_iter = nt_per_exposure + nt_per_deadtime

    # idk fuck around
    iter_skips = BitArray{1}(undef, nt_per_iter)
    iter_skips[1:nt_per_exposure] .= false
    iter_skips[nt_per_exposure+1:end] .= true

    # allocate memory
    return repeat(iter_skips, nt_per_night ÷ nt_per_iter)
end

function simulate_observations(obs::ObservationPlan, spec::SpecParams;
                               N::Int=128, snr::T=Inf, new_res::T=NaN) where T<:AF
    # assertions
    @assert (isnan(new_res) | (new_res < 7e5))

    # get number of time steps needed per exposure
    nt_per_exposure, nt_per_deadtime = get_nt(obs)

    # get number of time steps needed per night
    nt_per_iter = nt_per_exposure + nt_per_deadtime
    nt_per_night = obs.obs_per_night * nt_per_iter

    # get BitArray of epochs to skip
    skip_times = get_skip_times(nt_per_night, nt_per_exposure, nt_per_deadtime)

    # set up disk params instance
    disk = DiskParams(N=N, Nt=nt_per_night)

    # synthesize the spectra at default resolution
    Nλ = length(spec.lambdas)
    prof = ones(Nλ)
    outspec = zeros(Nλ, disk.Nt)
    disk_sim(spec, disk, prof, outspec, skip_times=skip_times)

    # allocate memory for binned spectra
    flux_binned = zeros(Nλ, obs.N_obs)
    for i in 1:obs.N_obs
        inds = (i-1) * nt_per_iter + 1: i * nt_per_iter
        flux_binned[:, i] = sum(outspec[:, inds], dims=2) ./ sum(.!skip_times[inds])
    end

    # degrade the resolution
    if !isnan(new_res)
        flux_binned = convolve_gauss(spec.lambdas, flux_binned, new_res=new_res)
    end

    # add noise to specified snr per res element
    flux_binned = add_noise(flux_binned, snr)

    # calculate ccf and velocities
    v_grid, ccf1 = calc_ccf(spec.lambdas, flux_binned, spec, normalize=true)
    rvs, sigs = calc_rvs_from_ccf(v_grid, ccf1)
    return rvs, sigs
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
module PrivateMultiplicativeWeights

using
    Distributions: Laplace, wsample
using
    Printf,
    Hadamard,
    LinearAlgebra,
    Random,
    IterTools,
    Statistics,
    Distributed,

export
    mwem,
    MWParameters,
    Tabular,
    Histogram,
    HistogramQueries,
    SeriesRangeQueries,
    RangeQueries,
    RangeQuery,
    Parities,
    FactorParities,
    Interval,
    gosper,
    BinaryItr,

    maximum_error,
    kl_divergence_error,
    mean_squared_error,
    queriesMatrix,

    evaluate,
    get,

    kl_divergence,
    complete_way_marginals_indices,
    fourierCoefficients,
    calculateMarginal,
    Marginals,

    get_query_vector_from_interval,
    hadamard_basis_vector

import
    Base: eltype, length, iterate

include("interface.jl")
include("histogram.jl")
include("rangequeries.jl")
include("factors.jl")
include("gosper.jl")
include("parities.jl")
include("error.jl")
include("mwem.jl")

"""
    PrivateMultiplicativeWeights

A simple and practical algorithm for differentially private data analysis.
"""
PrivateMultiplicativeWeights

end # module

max_stars_repo_name: giladroyz/MWEM-project

max_issues_repo_name: giladroyz/MWEM-project

max_forks_repo_name: giladroyz/MWEM-project


text:
include("run.jl")
include(joinpath(SLIDE_DIR,"src","build","disagg.jl"))

io = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","disagg","partition.yml"));
shr = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","disagg","shr.yml"));
shr[:utd] = filter_with(shr[:utd], set; extrapolate = true)

cal = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","disagg","cal.yml"));
:va in propertynames(cal[:va0]) && (cal[:va0] = edit_with(unstack(copy(cal[:va0]), :va, :value), Replace.(Symbol.(set[:va]), missing, 0.0)))
set[:notrd] = setdiff(set[:g], unique(shr[:utd][:,:g]))

bdisagg_out = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","disagg","disagg_out.yml"));
bdisagg_int = read_to_check(joinpath(SLIDE_DIR,"data","readfiles","benchmark","disagg","disagg_int.yml"));
bdisagg = merge(bdisagg_out,bdisagg_int)

disagg_bench = Dict()

d = merge(copy(shr),copy(cal),Dict(
        :r => fill_with((r = set[:r],), 1.0),
        (:yr,:r,:g) => fill_with((yr = set[:yr], r = set[:r], g = set[:g]), 1.0)));

(d,set) = disagg!(d,set)

# function _disagg_hhadj(d::Dict)
#     dh = Dict(k => edit_with(copy(d[k]), Rename(:g,:s))
#         for k in [:c0,:ld0,:kd0,:yh0,:bopdef0,:ta0,:a0,:tm0,:ty0,:g0,:i0,:m0])
#     dh[:ys0] = copy(d[:ys0])

#     d[:hhadj] = dh[:c0] -
#         combine_over(dh[:ld0] + dh[:kd0] + dh[:yh0], :s) - dh[:bopdef0] -
#         combine_over(dh[:ta0]*dh[:a0] + dh[:tm0]*dh[:m0] + dh[:ty0]*combine_over(dh[:ys0],:g), :s) +
#         combine_over(dh[:g0] + dh[:i0], :s)

#     return d[:hhadj]
# end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
module Day13
export solvePart1, solvePart2

function getRawInput()
    open("day13input.txt") do f
        return read(f, String)
    end
end

function parseInput(s)
    lines = split(s, '\n')
    height = length(lines)
    width = max(map(length, lines)...)
    grid = Array{Any}(undef, height, width)
    for (row, line) in enumerate(lines)
        for (col, c) in enumerate(line)
            if c == '^'
                grid[row, col] = ('N', 'l', '|')
            elseif c == 'v'
                grid[row, col] = ('S', 'l', '|')
            elseif c == '<'
                grid[row, col] = ('W', 'l', '-')
            elseif c == '>'
                grid[row, col] = ('E', 'l', '-')
            else
                grid[row, col] = c
            end
        end
    end
    return grid
end

function getCartPositions(grid)
    carts = []
    for r in 1:size(grid, 1)
        for c in 1:size(grid, 2)
            ch = grid[r, c]
            if typeof(ch) == Tuple{Char, Char, Char}
                push!(carts, (r, c))
            end
        end
    end
    return carts
end

function getVectorForDirection(dir)
    if dir == 'N'
        return [-1, 0]
    elseif dir == 'W'
        return [0, -1]
    elseif dir == 'S'
        return [1, 0]
    elseif dir == 'E'
        return [0, 1]
    end
end

function getNextDirection(cur, turn)
    directions = ['N', 'E', 'S', 'W']
    ind = findfirst(dir -> dir == cur, directions)
    if turn == 'l'
        ind = ind - 1
        if ind < 1
            ind += 4
        end
        return directions[ind]
    elseif turn == 's'
        return cur
    elseif turn == 'r'
        ind = ind + 1
        if ind > 4
            ind -= 4
        end
        return directions[ind]
    end
    return ' '
end

function getNextTurn(cur)
    if cur == 'l'
        return 's'
    elseif cur == 's'
        return 'r'
    elseif cur == 'r'
        return 'l'
    end
    return ' '
end

function getNewCornerDirection(d, corner)
    if corner == '\\'
        if d == 'N'
            return 'W'
        elseif d == 'W'
            return 'N'
        elseif d == 'S'
            return 'E'
        elseif d == 'E'
            return 'S'
        end
    elseif corner == '/'
        if d == 'N'
            return 'E'
        elseif d == 'E'
            return 'N'
        elseif d == 'S'
            return 'W'
        elseif d == 'W'
            return 'S'
        end
    end
end

function moveCart(row, col, grid)
    if typeof(grid[row, col]) != Tuple{Char, Char, Char}
        # there's no cart here
        return false, nothing
    end
    dir, turn, under = grid[row, col]
    loc = [row, col]
    vec = getVectorForDirection(dir)
    oldLoc = loc
    loc = loc + vec

    # Restore the square we moved from
    grid[row, col] = under

    row, col = loc

    newSquare = grid[row, col]
    if newSquare == '|' || newSquare == '-'
        grid[row, col] = (dir, turn, newSquare)
    elseif newSquare == '+'
        dir = getNextDirection(dir, turn)
        turn = getNextTurn(turn)
        grid[row, col] = (dir, turn, newSquare)
    elseif newSquare == '/' || newSquare == '\\'
        dir = getNewCornerDirection(dir, newSquare)
        grid[row, col] = (dir, turn, newSquare)
    elseif typeof(newSquare) == Tuple{Char, Char, Char}
        grid[row, col] = grid[row, col][3]
        # We collided!
        x = col - 1
        y = row - 1
        println("collided at $x,$y")
        flush(stdout)
        return true, (x, y)
    end
    return false, nothing
end

function solvePart1()
    s = getRawInput()
    grid = parseInput(s)
    collided = false
    row, col = 1, 1
    xy = (-1, -1)
    while !collided
        cartPositions = getCartPositions(grid)
        for (row, col) in cartPositions
            collided, xy = moveCart(row, col, grid)
            if collided
                break
            end
        end
    end
    return xy
end

function solvePart2()
    s = getRawInput()
    grid = parseInput(s)
    row, col = 1, 1
    cartPositions = getCartPositions(grid)
    while length(cartPositions) > 1
        for (row, col) in cartPositions
            for (row, col) in cartPositions
                moveCart(row, col, grid)
            end
        end
        cartPositions = getCartPositions(grid)
    end
    row, col = cartPositions[1]
    x = col - 1
    y = row - 1
    println("Cart ends at $x,$y")
    return (x, y)
end

end


text:
### ============== ### ============== ###
##    Behavioural Rules Model          ##
##    Martin Zumaya Hernandez          ##
##    EXAMPLE SIMULATION SCRIPT        ##
### ============== ### ============== ###

### ============ INCLUDE PACKAGES ============ ###

@everywhere using CollectiveDynamics.BehaviouralRules

### ============ SYSTEM'S PARAMETERS ============ ###

@everywhere dt = 1.0 # time step
@everywhere ρ = 0.3 # initial density

@everywhere v0 = 1.0 # speed

@everywhere η = 0.15 # noise intensity

@everywhere θ = 40.0 # maximum turn
@everywhere δ = 0.05 # deviation from aligned velocity

### ============ METRIC BEHAVIORAL THRESHOLDS ============ ###

n = parse(Int64, ARGS[1]) # number of particles

o = parse(Float64, ARGS[2]) # size of the orientation zone relative to system size
a = parse(Float64, ARGS[3]) # size of the attraction zone relative to system size

T = parse(Int64, ARGS[4]) # 10 ^ T iterations

rep = parse(Int64, ARGS[5]) # ensemble index

init = ARGS[6] # random or aligned initial velocities

@eval @everywhere init_e = $init

@eval @everywhere N = $n
@eval @everywhere Δo = $o
@eval @everywhere Δa = $a

@everywhere L  = cbrt(N / ρ) # size of box

@everywhere zor = 0.1 # zone of repulsion
# @everywhere zor = 1.0 # zone of repulsion
@everywhere zoo = zor + Δo*L # zone of orientation
@everywhere zoa = zoo + Δa*L # zone of attraction

println("rep:\t", zor, "\torient:\t", zoo, "\tattr:\t", zoa )

### ============ SHARED ARRAY INITIALIZATION ============ ###

pos = SharedArray{Float64}(3N) # particles positions
vel = SharedArray{Float64}(3N) # array of particles' velocities

v_int = SharedArray{Float64}(3N) # total signal
v_r   = SharedArray{Float64}(3N) # repulsion metric interactions
v_o   = SharedArray{Float64}(3N) # orientation metric interactions
v_a   = SharedArray{Float64}(3N) # attraction metric interactions

Rij = SharedArray{Float64}(N,N)

output_path = ""

if init_e == "R" # random initial conditions

    output_path = set_output_data_structure("BEHAV_R_01_N_015", N, ARGS[2], ARGS[3])
    println(output_path)

    ### ============ RANDOM INITIAL CONDITIONS ============ ###
    for i in 1:length(pos)
        pos[i] = 2*rand()*L - L
        vel[i] = 2*rand() - 1
    end

elseif init_e == "A" # ordered state initial condition

    output_path = set_output_data_structure("BEHAV_R_01_N_015_AL", N, ARGS[2], ARGS[3])
    println(output_path)

    ### ============ RANDOM POSITIONS BUT ALIGNED VELOCITIES ============ ###

    vel_0 = [2*rand() - 1, 2*rand() - 1, 2*rand() - 1]

    for i in 1:3:length(pos)
        pos[i]   = 2*rand()*L - L
        pos[i+1] = 2*rand()*L - L
        pos[i+2] = 2*rand()*L - L
        vel[i]   = vel_0[1] + δ*rand() - δ
        vel[i+1] = vel_0[2] + δ*rand() - δ
        vel[i+2] = vel_0[3] + δ*rand() - δ
    end

end

### ============ VELOCITY NORMALIZATION ============ ###

for i in 1:3:length(vel)
    norm = sqrt(vel[i]^2 + vel[i+1]^2 + vel[i+2]^2)
    vel[i] /= norm
    vel[i+1] /= norm
    vel[i+2] /= norm
end

### ============ SET OUTPUT ============ ###

pos_file = open(joinpath(output_path,"pos_$(rep).dat"), "w+")
vel_file = open(joinpath(output_path,"vel_$(rep).dat"), "w+")

# write initial conditions
write(pos_file, pos)
write(vel_file, vel)
println("//////// ", 1)

### ============ TIME EVOLUTION ============ ###

times = [convert(Int, exp10(i)) for i in 0:T]

for i in 1:(length(times) - 1)

    for t in (times[i]+1):times[i+1]

        evolve_system(pos, vel, v_int, v_r, v_o, v_a, Rij, N, zor, zoo, zoa, η, θ, v0)

        if t % times[i] == 0 || t % div(times[i], exp10(1)) == 0
            println("//////// ", t)
            write(pos_file, pos)
            write(vel_file, vel)
        end
    end

end

close(pos_file)
close(vel_file)

rmprocs(workers())

println("Done all")

### ============== ### ============== ### ============== ###


text:
using MetaGen
include("example.jl")

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
get(ENV, "TRAVIS_OS_NAME", "")       == "linux" || exit()
get(ENV, "TRAVIS_JULIA_VERSION", "") == "1.0"   || exit()

using Coverage

cd(joinpath(dirname(@__FILE__), "..")) do
    Codecov.submit(Codecov.process_folder())
end


text:
# returns the list of subexpressions which a given tape depends on directly
function list_subexpressions(nd::Vector{NodeData})

    # TODO: code overlap with compute_gradient_sparsity
    indices = Set{Int}()
    for k in 1:length(nd)
        nod = nd[k]
        if nod.nodetype == SUBEXPRESSION
            push!(indices, nod.index)
        end
    end

    return sort(collect(indices))
end

# order the subexpressions which main_expressions depend on
# such that we can run forward mode in this order
function order_subexpressions(main_expressions::Vector{Vector{NodeData}},subexpressions::Vector{Vector{NodeData}})
    nsub = length(subexpressions)
    computed = falses(nsub)
    dependencies = Dict{Int,Vector{Int}}()
    to_visit = collect(nsub+1:nsub+length(main_expressions))
    depended_on_by = [Set{Int}() for i in 1:nsub]

    while !isempty(to_visit)
        idx = pop!(to_visit)
        if idx > nsub
            li = list_subexpressions(main_expressions[idx-nsub])
        else
            computed[idx] && continue
            li = list_subexpressions(subexpressions[idx])
            computed[idx] = true
        end
        dependencies[idx] = li
        for k in li
            if idx > nsub
                push!(depended_on_by[k], idx-nsub)
            else
                union!(depended_on_by[k], depended_on_by[idx])
            end
            push!(to_visit,k)
        end
    end

    # now order dependencies
    I = Int[]
    J = Int[]
    for (idx,li) in dependencies
        for k in li
            push!(I,idx)
            push!(J,k)
        end
    end
    N = nsub+length(main_expressions)
    sp = sparse(I,J,ones(length(I)),N,N)
    cmap = Array{Int}(undef, N)
    order = reverse(Coloring.reverse_topological_sort_by_dfs(sp.rowval,sp.colptr,N,cmap)[1])
    # remove the subexpressions which never appear anywhere
    # and the indices of the main expressions
    order_filtered = collect(filter(idx -> (idx <= nsub && computed[idx]), order))
    # also generate an individual order for each main expression
    individual_order = [Int[] for i in 1:length(main_expressions)]
    for o in order_filtered
        for i in depended_on_by[o]
            push!(individual_order[i],o)
        end
    end

    return order_filtered, individual_order
end

export list_subexpressions, order_subexpressions


text:
import Base: convert
import .Khepri
import .Khepri:
    Box,
    Circle,
    ClosedLine,
    Cuboid,
    Line,
    Loc,
    Path,
    Point,
    Polygon,
    Rectangle,
    RightCuboid,
    Surface,
    SurfaceCircle,
    SurfacePath,
    SurfacePolygon,
    SurfaceRectangle,
    SurfaceRegularPolygon,
    Sphere,
    Vec,

    box,
    circle,
    circular_path,
    closed_line,
    closed_polygonal_path,
    cuboid,
    cx,
    cy,
    cz,
    in_world,
    line,
    loc_from_o_vz,
    path_vertices,
    polygon,
    polygonal_path,
    rectangle,
    rectangular_path,
    right_cuboid,
    sphere,
    subtraction,
    surface,
    surface_circle,
    surface_path,
    surface_polygon,
    surface_rectangle,
    surface_regular_polygon,
    u0,
    vx,
    vxy,
    vxyz,
    vxz,
    vy,
    vz,
    world_cs,
    xy,
    xyz,
    xz,
    yz

include("khepri/kernel.jl")
include("khepri/polygon_2.jl")


text:
#!/usr/bin/env julia
using Printf
using Distributed
using TimerOutputs

function ζ_dist(N,a)
  s = @distributed (+) for i ∈ 1:N
    (1.0/i)^a
  end
  s
end

function ζ_simd(N,a)
  s = 0.0
  @simd for i ∈ 1:N
    s+=(1.0/i)^a
  end
  s
end

function ζ_gen_thr_range(mstart,mend,e)
    s = 0
    cls  = 8  # cache line size in units of Float64 storage size
              # so we can avoid false sharing
    Nthr = Threads.nthreads()
    psum = zeros(Float64,Nthr*cls)
    ms   = zeros(Float64,Nthr*cls)
    me   = zeros(Float64,Nthr*cls)

    # preparing each thread, start, end, and partial sums
    for i=1:Nthr
        ms[i*cls]   = (N/Nthr)*i
        me[i*cls]   = 1 + (N/Nthr)*(i-1)
        psum[i*cls] = 0.0
    end

    # Threaded inner loop, each thread has no dependence upon others
    Threads.@threads for t=1:Nthr
        inner_gen_cpu1!(psum,ms,me,cls,2)
    end

    # now sum up the partial sums and return s
    s = sum(psum)
    s
end

function ζ_2_sqr_thr_range(mstart,mend)
    s = 0
    cls  = 8  # cache line size in units of Float64 storage size
              # so we can avoid false sharing
    Nthr = Threads.nthreads()
    psum = zeros(Float64,Nthr*cls)
    ms   = zeros(Float64,Nthr*cls)
    me   = zeros(Float64,Nthr*cls)

    # preparing each thread, start, end, and partial sums
    for i=1:Nthr
        ms[i*cls]   = (N/Nthr)*i
        me[i*cls]   = 1 + (N/Nthr)*(i-1)
        psum[i*cls] = 0.0
    end

    # Threaded inner loop, each thread has no dependence upon others
    Threads.@threads for t=1:Nthr
        inner_2_sqr_cpu1!(psum,ms,me,cls)
    end

    # now sum up the partial sums and return s
    s = sum(psum)
    s
end

function inner_gen_cpu1!(psum,ms,me,cls,x)
    ndx     = Threads.threadid()
    @simd for i=ms[ndx*cls]:-1:me[ndx*cls]
        psum[ndx*cls] += (1.0/i)^x
    end
    nothing
end

function inner_2_sqr_cpu1!(psum,ms,me,cls)
    ndx     = Threads.threadid()
    @simd for i=ms[ndx*cls]:-1:me[ndx*cls]
        psum[ndx*cls] += (1.0/i)*(1.0/i)
    end
    nothing
end
to = TimerOutput()

N = 16000000000
#N = 16000000
#Nthr = 1024
st = 1
s =  zeros(Float64,32)
pd = zeros(Float64,32)
ndx = 1

to = TimerOutput()

s[1] = ζ_dist(N,2)
s[2] = ζ_simd(N,2)
s[3] = ζ_gen_thr_range(N,1,2)
s[4] = ζ_2_sqr_thr_range(N,1)

for i=1:4
  @printf("s[%i] = %f\n",i,s[i])
end

@timeit to "ζ_dist"       pd[ndx] = ζ_dist(N,2)
ndx+=1

@timeit to "ζ_simd"       pd[ndx] = ζ_simd(N,2)
ndx+=1

@timeit to "ζ_gen_thr_range"     pd[ndx] = ζ_gen_thr_range(N,1,2)
ndx+=1

@timeit to "ζ_2_sqr_thr_range"    pd[ndx] = ζ_2_sqr_thr_range(N,1)

for i=1:4
  @printf("pd[%i] = %f\n",i,pd[i])
end

show(to)
@printf("\n");

max_stars_repo_name: joelandman/rzf

max_issues_repo_name: joelandman/rzf

max_forks_repo_name: joelandman/rzf


text:
using Agents

mutable struct User <: AbstractAgent
    id::Int
    pos::Int
    activity
    concentration
    quality_perception
end

function draw_activity()
    # from beta distribution
end

function draw_concentration()
    # from Poisson distribution
end

function draw_quality_vec()
    # some sort of multivariate distribution?
end

# quality vector
# evaluation space
# User.quality_perception ~ Post.quality

function get_opinion_distribution()
end

function get_evaluation_satisfaction()

end

mutable struct Post
    quality
    creation_time
    evaluation
    score
    views
    relevance
end

function compute_relevance(model, post, gravity)
    numerator = sum(post.quality)
    denominator = (model.step - post.creation_time)^gravity
    return numerator / denominator
end

function quality_transformation(quality_vec)
    [1 / (1 + Base.MathConstants.e^(-0.5 * x)) for x in quality_vec]
end

function form_consensus_opinion(post, user)
    transformed_post_quality = quality_transformation(post.quality)
    transformed_quality_perception = quality_transformation(user.quality_perception)
    sum([x^y for (x, y) in zip(transformed_post_quality, transformed_quality_perception)]) / length(user.quality)
end

function form_dissent_opinion(post, user)
    q1 = quality_transformation(post.quality)
    q2 = quality_transformation(user.quality_perception)

    term = sqrt(sum([abs(i - j)^2 for (i, j) in zip(q1, q2)])) / length(q1)
    return 1 - term
end

q1 = [3, 4, 5]
q2 = [5, 1, 7]
term = sqrt(sum([abs(i - j)^2 for (i, j) in zip(q1, q2)])) / length(q1)

using Random

properties = Dict()

ABM(User, scheduler = random_activation, properties = properties)

# number of Users
# number of initial Posts
# new Posts per iteration

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


text:
using Radio, Winston

# generate 10,000 random QPSK symbols
symbols = pskmod( 10000, 4 )
# create some gaussian noise and add it to the symbols
noise  = wgn( length( symbols ), 10, "dBm", 1.0, true )
signal = symbols .+ noise

constellation = plot_constellation( signal )
setattr( constellation, title = "QPSK Modulation" )

display( constellation )
file( constellation, "QPSK.png" )


text:
"""
Bridges `CP.Element` to MILP constraints by using a unary encoding of the
index in the array.
"""
struct Element2MILPBridge{T} <: MOIBC.AbstractBridge
    vars::Vector{MOI.VariableIndex}
    vars_bin::Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    con_unary::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    con_choose_one::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
    con_value::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
end

function MOIBC.bridge_constraint(
    ::Type{Element2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Element{T},
) where {T}
    return MOIBC.bridge_constraint(
        Element2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Element2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Element{T},
) where {T}
    vars, vars_bin = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:length(s.values)])

    f_scalars = MOIU.scalarize(f)
    f_value = f_scalars[1]
    f_index = f_scalars[2]

    con_unary = MOI.add_constraint(
        model,
        sum(T(i) * vars[i] for i in 1:length(s.values)) - f_index,
        MOI.EqualTo(zero(T))
    )

    con_choose_one = MOI.add_constraint(
        model,
        sum(one(T) .* vars),
        MOI.EqualTo(one(T))
    )

    con_value = MOI.add_constraint(
        model,
        sum(s.values[i] * vars[i] for i in 1:length(s.values)) - f_value,
        MOI.EqualTo(zero(T))
    )

    return Element2MILPBridge(vars, vars_bin, con_unary, con_choose_one, con_value)
end

function MOI.supports_constraint(
    ::Type{Element2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Element{T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{Element2MILPBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{Element2MILPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::Element2MILPBridge{T}, ::MOI.NumberOfVariables) where {T}
    return length(b.vars)
end

function MOI.get(
    b::Element2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_bin)
end

function MOI.get(
    ::Element2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return 3
end

function MOI.get(
    b::Element2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return copy(b.vars)
end

function MOI.get(
    b::Element2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return copy(b.vars_bin)
end

function MOI.get(
    b::Element2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return [b.con_unary, b.con_choose_one, b.con_value]
end


text:
module JungleHelperUITextTrigger

using ..Ahorn, Maple

@mapdef Trigger "JungleHelper/UITextTrigger" UITextTrigger(x::Integer, y::Integer, width::Integer=8, height::Integer=8, Dialog::String="",Flag::String="",TextX::Number=0.0,TextY::Number=0.0,FadeIn::Number=1.0,FadeOut::Number=1.0)

const placements = Ahorn.PlacementDict(
    "UI Text Display (Jungle Helper)" => Ahorn.EntityPlacement(
       UITextTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
module JuMPLPUtils
using SparseArrays
using MathProgBase
using JuMP
include("lp_course_utils.jl")
export getslack
export getoptranges
export getoptrange

end # module

max_stars_repo_name: EdvinAblad/JuMPLPUtils

max_issues_repo_name: EdvinAblad/JuMPLPUtils

max_forks_repo_name: EdvinAblad/JuMPLPUtils


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.
# This file includes parts from Combinatorics.jl in modified form.

# Version of: UTC 2019-05-26 10:30:37
# 89dc11ae-7f90-11e9-1984-5170e93e3915

# Do not edit this file, it is generated from the modules and will be overwritten!
# Edit the modules in the src directory and build this file with BuildSequences.jl!

tstdir = realpath(joinpath(dirname(@__FILE__)))
srcdir = joinpath(dirname(tstdir), "src")
tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)
module runtests
using Test, SeqTests, IntegerSequences, Nemo, IterTools, Combinatorics
# *** Abundant.jl *********
@testset "Abundant" begin
@test isAbundant(100800) == true
@test isAbundant(2402400) == true
@test isAbundant(49008960) == true
if is_oeis_installed()
SeqTest([L005101], 'L')
end
end
# *** AndreNumbers.jl *********
@testset "André" begin
@test isa(André(2, 10), Nemo.fmpz)
@test André(2, 10) == ZZ(-50521)
@test André(2, 50) == ZZ(-6053285248188621896314383785111649088103498225146815121)
@test V178963(30)  == ZZ(2716778010767155313771539)
@test V178964(40)  == ZZ(11289082167259099068433198467575829)
if is_oeis_installed()
V = [V000111, V178963, V178964, V181936, V250283]
for v in V SeqTest(v, 'V') end
end
V = [1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521]
generator = C000111()
for n in 1:10
v = take!(generator)
@test V[n]  == v
end
v = take!(generator)
close(generator)
@test isa(v, fmpz)
end
# *** BellNumbers.jl *********
@testset "Bell" begin
seq = fmpz[1, 1, 2, 5, 15, 52]
a = [BellTrans(6, k, seq) for k in 0:5]
b = [0, 52, 205, 210, 85, 15]
@test all(a .== b)
seq = [1, 1, 1, 1, 1, 1]
a = BellTriangle(6, seq)
b = fmpz[1, 0, 1, 0, 1, 1, 0, 1, 3, 1, 0, 1, 7, 6, 1, 0, 1, 15, 25]
@test all(a[1:5] .== b[1:5])
a = fmpz[1, 1, 2, 5, 15, 52, 203]
b = [BellNumber(n) for n in 0:6]
@test all(a .== b)
if is_oeis_installed()
SeqTest(V000110, 'V')
SeqTest([L000110], 'L')
T = [T264428, T137452, T132062, T265606]
SeqTest(T, 'T')
P = [T104556, T001497, T039683, T203412, T004747, T051141,
T119274, T000369, T051142, T137513]
SeqTest(P, 'P')
end
end
# *** BernoulliNumbers.jl *********
@testset "BernoulliNum" begin
@test BernoulliInt(2, 10) == ZZ(79360)
@test BernoulliInt(3, 30) == ZZ(-7708110416280010548302670)
@test BernoulliInt(4, 40) == ZZ(-44494882577309421077208834962882560)
@test isa(BernoulliInt(3, 30), Nemo.fmpz)
@test isa(BernoulliIntList(2, 20)[end], Nemo.fmpz)
@test isa(Bernoulli(0), Nemo.fmpq)
@test isa(V027641(10), Nemo.fmpz)
@test isa(V195441(10), Nemo.fmpz)
for m in 1:20
@test BernoulliInt(m, 200) == BernoulliIntList(m, 200 + 1)[end]
end
for n in 0:20
@test denominator(Rational(Bernoulli(2 * n))) == ClausenNumber(n)
end
t = ZZ(8622529719094842064796322984685715031642180319435676189471082876882585178647210)
@test V195441(10000) == t
l = fmpz[0, 1, 2, 3, 8, 25, 96, 427, 2176]
@test all(L065619(9) .== l)
if is_oeis_installed()
V = [V195441, V281586, V281588, V027641]
for v in V SeqTest(v, 'V') end
end
end
# *** BinaryInteger.jl *********
@testset "BinaryInteger" begin
@test V295513(0) == -1
@test V295513(1) == -1
@test V295513(2) == 0
@test V295513(3) == 2
if is_oeis_installed()
V = [V001855, V003314, V033156, V054248, V061168, V097383]
for v in V SeqTest(v , 'V', 1) end
V = [V083652, V123753]
for v in V SeqTest(v , 'V', 0) end
end
end
# *** BinaryQF.jl *********
Data = Dict{Int, Array{fmpz}}(
002476 => [7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97,103],
008784 => [1,  2,  5, 10, 13, 17, 25, 26, 29, 34, 37, 41],
031363 => [1,  4,  5,  9, 11, 16, 19, 20, 25, 29, 31, 36],
034017 => [1,  3,  7, 13, 19, 21, 31, 37, 39, 43, 49, 57],
035251 => [1,  2,  4,  7,  8,  9, 14, 16, 17, 18, 23, 25],
038872 => [5, 11, 19, 29, 31, 41, 59, 61, 71, 79, 89,101],
038873 => [2,  7, 17, 23, 31, 41, 47, 71, 73, 79, 89, 97],
042965 => [1,  3,  4,  5,  7,  8,  9, 11, 12, 13, 15, 16],
057126 => [1,  2,  7, 14, 17, 23, 31, 34, 41, 46, 47, 49],
057127 => [1,  2,  3,  6,  9, 11, 17, 18, 19, 22, 27, 33],
068228 => [13,37, 61, 73, 97,109,157,181,193,229,241,277],
084916 => [1,  4,  6,  9, 13, 16, 22, 24, 25, 33, 36, 37],
089270 => [1,  5, 11, 19, 29, 31, 41, 55, 59, 61, 71, 79],
141158 => [5, 11, 19, 29, 31, 41, 59, 61, 71, 79, 89,101],
242660 => [1,  4,  7,  9, 10, 13, 16, 18, 19, 22, 25, 27],
243655 => [1,  6, 13, 22, 33, 37, 46, 61, 69, 73, 78, 94],
244779 => [1,  2,  4,  7,  8, 11, 14, 16, 22, 23, 28, 29],
244780 => [1,  3,  5,  9, 11, 15, 23, 25, 27, 31, 33, 37],
244819 => [1,  3,  4,  7, 12, 13, 19, 21, 28, 31, 37, 39],
243168 => [1,  4,  9, 13, 16, 24, 25, 33, 36, 37, 49, 52],
244291 => [1,  4, 13, 24, 33, 37, 52, 61, 69, 73, 88, 97],
007522 => [7, 23, 31, 47, 71, 79,103,127,151,167,191,199]
)
L = [L002476, L008784, L031363, L034017, L035251, L038872,
L038873, L042965, L057126, L057127, L068228, L084916, L089270,
L141158, L242660, L243655, L244779, L244780, L244819, L243168,
L244291, L007522]
SeqNum(seq) =  parse(Int, string(seq)[end-5:end])
@testset "BinaryQF" begin
for seq in L
S = seq(300)
anum = SeqNum(seq)
data = Data[anum]
@test all(S[1:11] .== data[1:11])
end
end
# *** CantorMachines.jl *********
@testset "CantorMachines" begin
len = 64
println(); println("CantorMachine")
x, y, state = 0, 0, false
for n in 0:len
p = CantorPairing(x, y)
@test n == p
x, y, state = CantorMachine(x, y, state)
end
println(); println("CantorBoustrophedonicMachine")
x, y = 0, 0
for n in 0:len
p = CantorBoustrophedonicPairing(x, y)
println("$n -> ($x, $y) -> $p")
@test n == p
x, y = CantorBoustrophedonicMachine(x, y)
end
println(); println("RosenbergStrongBoustrophedonicMachine")
x, y, state = 0, 0, 0
for n in 0:len
p = RosenbergStrongBoustrophedonicPairing(x, y)
@test n == p
x, y, state = RosenbergStrongBoustrophedonicMachine(x, y, state)
end
println()
end
# *** CarmichaelNumbers.jl *********
@testset "Carmichael" begin
@test ! isCarmichael(560)
@test isCarmichael(561)
@test ! isCarmichael(563)
@test isweakCarmichael(561)
@test !isweakCarmichael(563)
@test isweakCarmichael(625)
if is_oeis_installed()
L = [L002997, L225498]
SeqTest(L, 'L')
end
end
# *** ClausenNumbers.jl *********
@testset "Clausen" begin
C = ClausenNumberList(800)
@test C[125] == 30
@test C[781] == 32695402455500348373810
@test C[794] == 6
@test isa(C[781], Nemo.fmpz)
@test isa(ClausenNumber(10), Nemo.fmpz)
end
if is_oeis_installed()
SeqTest([L002445], 'L')
SeqTest(V002445, 'V', 0)
SeqTest(V027642, 'V', 0)
end
# *** CombinationsIterator.jl *********
C = collect(Combinations([2, 3, 5]))
@testset "Combinations" begin
@test C == [[2], [3], [5], [2, 3], [2, 5], [3, 5], [2, 3, 5]]
end
# *** Compositions.jl *********
# *** Counts.jl *********
indicators = [isPositive, isEven, isSquare, isPrime]
indicatorNames = ["isPositive", "isEven", "isSquare", "isPrime"]
len = 14
@testset "Counts" begin
@test Nth(96, isPrime) == 503
@test Nth(97, isPrime) == 509
@test Nth(98, isPrime) == 521
@test Count(96, isPrime) == 24
@test Count(97, isPrime) == 25
@test Count(98, isPrime) == 25
@test List(24, isPrime)[end] == 89
@test List(25, isPrime)[end] == 97
for (i, isA) in enumerate(indicators)
for n in 1:len
@test isA(n) == (Nth(Count(n, isA), isA) == n)
@test     n  ==  Count(Nth(n, isA), isA)
end
end
a = [A257993(n) for n in 1:10]
b = [1, 2, 1, 2, 1, 3, 1, 2, 1, 2]
@test all(a .== b)
if is_oeis_installed()
L = [
L000961, L002808, L005117, L013928, L246547, L246655
]
SeqTest(L, 'L')
end
end
# *** CyclotomicBinaryForms.jl *********
@testset "CyclotomicForms" begin
if is_oeis_installed()
L = [L206942, L206864, L299498, L299928, L299929, L299930,
L299214, L296095, L293654]
SeqTest(L, 'L')
end
end
# *** DedekindEta.jl *********
@testset "DedekindEta" begin
@test DedekindEtaPowers(0, 1) == fmpz[]
@test isa(DedekindEtaPowers(30, 1)[10], fmpz)
@test isa(DedekindEtaPowers(30, -1)[10], fmpz)
@test RamanujanTau(20) == -7109760
P0 = [1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42, 56, 77, 101]
P1 = PartitionNumberList(14)
@test all(P0 .== P1)
if is_oeis_installed()
L = [L010815, L002107, L010816, L000727, L000728, L000729, L000730,
L000731, L010817, L010819, L000735, L010820, L010821, L010822,
L000739, L010823, L010824, L010825, L010826, L010827, L010828,
L010829, L000594, L010830, L010831, L010832, L010833, L010834,
L010835, L010836, L010837, L010840, L010838, L010839, L010841,
L000041, L000712, L000716, L023003, L023004, L023005, L023006,
L023007, L023008, L023009, L023010, L005758, L023011, L023012,
L023013, L023014, L023015, L023016, L023017, L023018, L023019,
L023020, L023021, L006922, L082556, L082557, L082558, L082559]
SeqTest(L, 'L')
end
end
# *** DelehamDelta.jl *********
Data = Dict{Int, Array{fmpz}}(
084938 => [1, 0, 1,  0,  1,  1,   0,   2,   2,  1],
060693 => [1, 1, 1,  2,  3,  1,   5,  10,   6,  1],
106566 => [1, 0, 1,  0,  1,  1,   0,   2,   2,  1],
094665 => [1, 0, 1,  0,  1,  3,   0,   4,  15, 15],
090238 => [1, 0, 1,  0,  2,  1,   0,   6,   4,  1],
225478 => [1, 3, 4, 21, 40, 16, 231, 524, 336, 64],
055883 => [1, 1, 1,  2,  4,  2,   5,  15,  15,  5],
184962 => [1, 0, 1,  0,  1,  1,   0,   3,   3,  1],
088969 => [1, 0, 1,  0,  1,  3,   0,   5,  20, 21],
090981 => [1, 1, 0,  1,  1,  0,   1,   4,   1,  0],
011117 => [1, 1, 0,  1,  1,  0,   1,   2,   3,  0]
)
@testset "Deléham" begin
n = 7
B = [bell(n) * binomial(n, j) for j in 0:n]
R = Row(T055883(n+1), n)
@test all(B .== R)
a = fmpz[0, 5040, 2208, 828, 272, 70, 12, 1]
b = Row(T090238(8), 7)
@test all(a .== b)
SeqNum(seq) =  parse(Int, string(seq)[end-5:end])
Seq = [T084938, T060693, T106566, T094665, T090238, T225478, T055883,
T184962, T088969]
for seq in Seq
S = seq(20)
anum = SeqNum(seq)
data = Data[anum]
@test all(S[1:9] .== data[1:9])
end
end
# *** EulerTransforms.jl *********
@testset "EulerTransform" begin
if is_oeis_installed()
L = [L107895, L190905, L061256, L275585, L006171, L290351, L052847]
SeqTest(L, 'L')
V = [V107895, V190905, V061256, V275585, V006171]
for v in V SeqTest(v, 'V', 0) end
end
end
# *** Fibonacci.jl *********
@testset "Fibonacci" begin
@test isa(V000045(30), Nemo.fmpz)
end
# *** FigurativeNumbers.jl *********
@testset "Figurative" begin
@test V002411(1000) == 500500000
if is_oeis_installed()
V = [V014107, V067998, V001477, V000217, V000290, V000326,
V000384, V000566, V000567, V001106, V001107, V000292,
V000330, V002411, V002412, V002413, V007584, V007585]
for v in V  SeqTest(v, 'V') end
end
end
# *** GaussFactorials.jl *********
@testset "GaussFactorial" begin
if is_oeis_installed()
V = [V001783, V124441, V124442, V066570]
for v in V SeqTest(v, 'V', 1) end
W = [V055634, V232980, V232981, V232982]
for w in W SeqTest(w, 'V', 0) end
end
end
# *** GeneralizedBinomial.jl *********
@testset "Binomial" begin
for n in 0:10, k in 0:n
@test Binomial(n, k) == div(fac(n), (fac(n - k) * fac(k)))
@test Binomial(n, k) == Pascal(n, k)
end
end
# *** HighlyAbundant.jl *********
# *** Hyper1F1.jl *********
@testset "Hyper1F1" begin
@test V000255(10) == ZZ(16019531)
@test V000262(10) == ZZ(58941091)
@test V001339(10) == ZZ(98641011)
@test V007060(10) == ZZ(871804170613555200)
@test V033815(10) == ZZ(1465957162768492800)
@test V099022(10) == ZZ(3984884716852972800)
@test V251568(10) == ZZ(123320412181)
if is_oeis_installed()
V = [V000255, V000262, V001339, V007060, V033815,
V099022, V251568]
for v in V SeqTest(v, 'V') end
end
end
# *** JacobiTheta.jl *********
@testset "JacobiTheta" begin
@test JacobiTheta3Powers(0, 1) == fmpz[]
@test isa(JacobiTheta3Powers(30, 1)[end], fmpz)
@test isa(JacobiTheta4Powers(30, 1)[end], fmpz)
@test L035016(999 + 1)[end]  == ZZ(-16565884160)
@test L035016(1000 + 1)[end] == ZZ(18365675328)
if is_oeis_installed()
L = [L000122, L002448, L004018, L104794, L005875, L213384, L000118,
L035016, L008452, L096727, L000132, L000141, L008451, L000143,
L000144, L008453, L000145, L276285, L276286, L276287, L004402,
L004406, L004407, L015128, L004403, L001934, L004404, L004405,
L004408, L004409, L004410, L004411, L004412, L004413, L004414,
L004420, L004421, L004415, L004416, L004417, L004418, L004419,
L004422, L004423, L004424, L004425]
SeqTest(L, 'L')
end
end
# *** Kolakoski.jl *********
@testset "Kolakoski" begin
K = KolakoskiList(100)
@test K[1]  == 1
@test K[33] == 2
@test K[72] == 2
generator = C000002()
for n in [1, 33, 72]
k = take!(generator)
@test K[n] == k
end
close(generator)
end
# *** NarayanaCows.jl *********
@testset "Narayana" begin
@test IterTools.nth(NarayanasCows(12), 12) == 5
end
# *** NumberTheory.jl *********
Data = Dict{Int,Array{fmpz}}(034386 => [1, 1, 2, 6, 6, 30, 30, 210, 210, 210],
061142 => [1, 1, 2, 2, 4, 2, 4, 2, 8, 4],
002110 => [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870],
000005 => [0, 1, 2, 2, 3, 2, 4, 2, 4, 3],
000010 => [0, 1, 1, 2, 2, 4, 2, 6, 4, 6],
000203 => [0, 1, 3, 4, 7, 6, 12, 8, 15, 13],
001222 => [0, 0, 1, 1, 2, 1, 2, 1, 3, 2],
001221 => [0, 0, 1, 1, 1, 1, 2, 1, 1, 1],
008683 => [0, 1, -1, -1, 0, -1, 1, -1, 0, 0],
181830 => [0, 1, 0, 0, 0, 1, 0, 2, 2, 2],
034444 => [1, 1, 2, 2, 2, 2, 4, 2, 2, 2])
@testset "NumTheory" begin
@test τ(7560) == 64
@test τ(46080) == 66
@test τ(25920) == 70
@test σ(7560) == 28800
@test σ(46080) == 159666
@test σ(25919) == 25920
@test σ(25920) == 92202
@test ω(7560) == 4
@test ω(46080) == 3
@test ω(25919) == 1
@test ω(25920) == 3
@test Ω(7560) == 8
@test Ω(46080) == 13
@test Ω(25919) == 1
@test Ω(25920) == 11
@test Radical(58564) == 22
@test Radical(58565) == 58565
FB(n::Int) = (r = 1; for k in 1:n ⊥(n, k) && (r = mods(r * k, n)) end; r)
FA(n::Int) = mods(∏([j for j in 1:n if ⊥(j, n)]), n)
for n in 1:20
@test FA(n) == FB(n)
end
if is_oeis_installed()
V = [V061142, V000005, V000010, V000203, V001222, V001221, V008683]
W = [V034386, V002110, V181830, V034444]
for v in V SeqTest(v, 'V', 1) end
L = [L003277]
SeqTest(L, 'L')
end
end
composita = [false, false, false, false, true, false, true, false]
@testset "Queries" begin
for n in 0:7
@test isComposite(n) == composita[n + 1]
end
end
# *** OEISUtils.jl *********
if oeis_isinstalled()
@info("OEIS data is installed as: " * oeis_path())
end
# *** OrthoPolynomials.jl *********
@testset "OrthoPoly" begin
@test isa(OrthoPoly(10, n -> 1, n -> n + 1)[end], fmpz)
@test isa(InvOrthoPoly(10, n -> 1, n -> n + 1)[end], fmpz)
@test RowSums(T217537(8)) == L217924(8)
if is_oeis_installed()
T = [ T066325, T049310, T137338, T104562, T037027, T049218, T159834,
T137286, T053120, T053117, T053121, T216916, T217537, T064189,
T202327, T111062, T099174, T111593, T064189 ]
SeqTest(T, 'T')
L = [L217924, L005425, L000085, L001464, L003723, L108624, L006229]
SeqTest(L, 'L')
end
end
# *** Partitions.jl *********
V080577(7)
println()
# *** PrimesIterator.jl *********
@testset "Primes" begin
@test PrimePi(2^14) == 1900
@test PrimePi(10^6) == 78498
end
# *** PrimeSwingFactorial.jl *********
@testset "PrimeSwingF" begin
for n in 0:999
S = PSfactorial(n)
B = Base.factorial(BigInt(n))
@test all(S .== B)
end
end
# *** Products.jl *********
@testset "FallingFact" begin
@test FallingFactorial(100, 100) == factorial(BigInt(100))
@test (100 ↓ 100) == factorial(BigInt(100))
@test FallingFactorial(333, 333) == factorial(BigInt(333))
@test FallingFactorial(111, 0) == 1
end
@testset "RisingFact" begin
@test RisingFactorial(11, 11) == 14079294028800
@test (11 ↑ 11) == 14079294028800
@test RisingFactorial(33, 33) == 31344295059422473624824839739793024055460338073600000000
@test RisingFactorial(111, 0) == 1
end
@testset "MultiFact" begin
a = [MultiFactorial(2, 1)(n) for n in 0:6]
b = [1, 1, 3, 15, 105, 945, 10395]
@test all(a .== b)
end
if is_oeis_installed()
V = [V000142, V000165, V007696, V001813, V047053, V001147, V008545,
V081125, V000407, V032031, V007559, V008544]
@testset "Products" begin
for v in V SeqTest(v, 'V') end
end
end
# *** QueensProblems.jl *********
@testset "QueensProblem" begin
levels = [1, 10, 72, 364, 1400, 3916, 7552, 9632, 7828, 4040, 724]
@test all(L319284(10) .== levels)
end
# *** RationalTrees.jl *********
@testset "EuclidTree" begin
S = [numerator(sum(r for r in EuclidTree(n))) for n in 1:9]
L = Nemo.fmpz[1, 5, 11, 23, 47, 95, 191, 383, 767] # A052940
@test all(S[1:9] .== L[1:9])
end
# *** RecordSearch.jl *********
# *** RiordanSquares.jl *********
@testset "RiordanSqr" begin
a = T039599(5)
b = [1, 1, 1, 2, 3, 1, 5, 9, 5, 1, 14, 28, 20, 7, 1]
@test all(a .== b)
a = T116392(5)
b = [1, 1, 1, 3, 4, 1, 7, 13, 7, 1, 19, 42, 32, 10, 1]
@test all(a .== b)
a = T172094(5)
b = [1, 1, 1, 3, 4, 1, 11, 17, 7, 1, 45, 76, 40, 10, 1]
@test all(a .== b)
a = T321620(5)
b = [1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 1]
@test all(a .== b)
a = T321621(5)
b = [1, 1, 1, 2, 3, 1, 4, 8, 5, 1, 9, 21, 18, 7, 1]
@test all(a .== b)
a = T321623(5)
b = T321623(5)
@test all(a .== b)
a = T321624(5)
b = [1, 1, 1, 3, 4, 1, 4, 10, 7, 1, 7, 24, 26, 10, 1]
@test all(a .== b)
a = T322942(5)
b = [1, 1, 1, 1, 2, 1, 3, 5, 3, 1, 5, 12, 10, 4, 1]
@test all(a .== b)
end
# *** SelfConvolutive.jl *********
@testset "SelfConvRec" begin
if is_oeis_installed()
L = [L005411, L005412, L006012, L006318, L047891, L062980, L082298,
L082301, L082302, L105523, L107716, L111529, L111530, L111531,
L111532, L111533, L146559, L167872]
SeqTest(L, 'L')
end
end
# *** SeqUtils.jl *********
# *** SeriesExpansion.jl *********
@testset "SeriesCoefficients" begin
a = L000032(14)
b = [1, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521]
@test all(a .== b)
a = L000073(15)
b = [1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927, 1705, 3136]
@test all(a .== b)
a = L000108(13)
b = [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012]
@test all(a .== b)
a = L000257(12)
b = [1, 1, 3, 12, 56, 288, 1584, 9152, 54912, 339456, 2149888, 13891584]
@test all(a .== b)
a = L000957(13)
b = [1, 1, 0, 1, 2, 6, 18, 57, 186, 622, 2120, 7338, 25724]
@test all(a .== b)
a = L001003(12)
b = [1, 1, 3, 11, 45, 197, 903, 4279, 20793, 103049, 518859, 2646723]
@test all(a .== b)
a = L001006(14)
b = [1, 1, 2, 4, 9, 21, 51, 127, 323, 835, 2188, 5798, 15511, 41835]
@test all(a .== b)
a = L001045(13)
b = [1, 1, 1, 3, 5, 11, 21, 43, 85, 171,341, 683, 1365]
@test all(a .== b)
a = L002426(12)
b = [1, 1, 3, 7, 19, 51, 141, 393, 1107, 3139, 8953, 25653]
@test all(a .== b)
a = L005043(12)
b = [1, 1, 0, 1, 1, 3, 6, 15, 36, 91, 232, 603]
@test all(a .== b)
xL006318(n) = Coefficients(G006318, n)
a = xL006318(12)
b = [1, 2, 6, 22, 90, 394, 1806, 8558, 41586, 206098, 1037718, 5293446]
@test all(a .== b)
a = L068875(10)
b = [1, 2, 4, 10, 28, 84, 264, 858, 2860, 9724]
@test all(a .== b)
end
# *** StirlingLahNumbers.jl *********
@testset "Stirling1" begin
for (n, S1) in enumerate(I132393(10))
SC = [stirlings1(n-1, k) for k in 0:n-1]
@test S1 == SC
end
end
@testset "Stirling2" begin
for (n, S2) in enumerate(I048993(10))
SC = [stirlings2(n-1, k) for k in 0:n-1]
@test S2 == SC
end
end
# *** SwingFactorial.jl *********
@testset "SwingFactorial" begin
for n in 0:999
S = Sfactorial(n)
B = Base.factorial(BigInt(n))
@test S == B
end
end
# *** Triangles.jl *********
BinTri(n) = [binomial(k, j) for k in 0:n - 1 for j in 0:k]
X(n, k) = ((k > n || k < 0) && return 0;
(n == 0 && k == 0) && return 1;
4 * X(n - 1, k - 1) + (4 * n - 1) * X(n - 1, k); )
T(n) = [X(j, k) for j in 0:n for k in 0:j]
@testset "Triangles" begin
t = T(6)
a = [3465, 8784, 7136, 2304, 256]
b = Row(t, 4)
@test all(a .== b)
t = BinTri(8)
a = [1, 2, 4, 8, 16, 32, 64, 128]
b = RowSums(t)
@test all(a .== b)
a = [1, 0, 0, 0, 0, 0, 0, 0]
b = RowSums(t, true)  # alternating sum
@test all(a .== b)
end
# *** UlamNumbers.jl *********
@testset "UlamNumbers" begin
@test UlamList(6)[end] == 8
end
# *** ZumkellerNumbers.jl *********
@testset "Zumkeller" begin
@test isZumkeller(17000) == true
@test isZumkeller(27472) == true
@test isZumkeller(29062) == false
@test isZumkeller(43464) == true
end
end # module


text:
################################################################################
# Copyright (C) Marcus Greiff 2020
#
# @file CodeOcean.jl
# @author Marcus Greiff
# @date January 2021
#
# This file simply one of the relevant Julia examples for the LCSS-publication,
# it runs all of the implemented examples pertaining to the FSF attitude
# control, and saves the plots to the ../results directory.
################################################################################

# This is needed to make GR work in CodeOcean
ENV["GKSwstype"] = "100";
# Load and activate the AerialVehicleControl.jl package
using Pkg; Pkg.activate(".");

################################################################################
############################# Run the tests ####################################
################################################################################
println("Running tests...")
include("test/runtests.jl")

################################################################################
################### Run the FSF controller Feedback examples ###################
################################################################################
fileNames = ["continuous_SO3", "robust_SO3", "continuous_SU2", "discontinuous_SU2", "robust_SU2"]
applyDist = [false, true, false, false, true];

for i = 1:5
    # Simulation with a sought controller type
    global controllerType  = i;

    # Plotting settings
    global showDisturbance = applyDist[i];
    global showPlot        = [1,2,3];
    global savePlot        = true;
    global namePlot        = fileNames[i];

    # Simulation setting
    global useDisturbance = showDisturbance;

    # Run example
    include("examples/example_FSF_attitude.jl")
end


text:
let
    using UnsteadyFlowSolvers

    # Set up a single vortex
    start_vc = 0.05 # Core radius
    start_vort = -1.
    re = 1e3
    n_particles = 100
    radius = 1
    nu = start_vort / re
    field = TwoDVFlowField(re)
    for i = 1 : n_particles
        theta = i * 2 * pi / n_particles
        vort = TwoDVVort(radius * sin(theta), radius * cos(theta),
            start_vort/n_particles, start_vc, 0., 0., 0.)
        push!(field.tev, vort)
    end
    # Solver settings
    dt = 0.015
    t0 = 0.
    tmax = 1.5
    t = t0
    nsteps = Int(ceil((tmax-t)/dt)) + 1

    for i = 1:nsteps
        t += dt
        UnsteadyFlowSolvers.wakeroll(field, dt)
        ave = mapreduce(x->x.vc, +, field.tev) / (n_particles)
        stddev = sqrt(mapreduce(x->(x.vc-ave)^2, +, field.tev) / (n_particles-1))
        println("Average core readius: ", ave, " std deviation = ", stddev)
    end
end


text:
#definition of warlpiris for ease.

#tests using the "modified warlpiri" numbers
W = Valid{4,0}
WP = Vnum{4,0}

const WP0000 = WP(0b0000)
const WP0001 = WP(0b0001)
const WP0010 = WP(0b0010)
const WP0011 = WP(0b0011)
const WP0100 = WP(0b0100)
const WP0101 = WP(0b0101)
const WP0110 = WP(0b0110)
const WP0111 = WP(0b0111)
const WP1000 = WP(0b1000)
const WP1001 = WP(0b1001)
const WP1010 = WP(0b1010)
const WP1011 = WP(0b1011)
const WP1100 = WP(0b1100)
const WP1101 = WP(0b1101)
const WP1110 = WP(0b1110)
const WP1111 = WP(0b1111)

#create a matrix that will be useful for testing purposes.
const wmatrix = [WP(a) → WP(b) for a in 0b0000:0b1111, b  in 0b0000:0b1111]

#=
  structure of wmatrix:

  [ 0b0000 → 0b0001 0b0000 → 0b0010 0b0000 → 0b0011 ... 0b0000 → 0b1111
    0b0001 → 0b0001 0b0001 → 0b0010 0b0001 → 0b0011 ... 0b0001 → 0b1111
    0b0010 → 0b0001 0b0010 → 0b0010 0b0010 → 0b0011 ... 0b0010 → 0b1111
    ...
    0b1111 → 0b0001 0b1111 → 0b0010 0b1111 → 0b0011 ... 0b1111 → 0b1111]
=#

#re-alias these numbers to make them even easier to read.
const W0_0 = WP0000 → WP0000
const W0_1 = WP0000 → WP0001
const W0_2 = WP0000 → WP0010
const W0_3 = WP0000 → WP0011
const W0_4 = WP0000 → WP0100
const W0_5 = WP0000 → WP0101
const W0_6 = WP0000 → WP0110
const W0_7 = WP0000 → WP0111
const W0_8 = WP0000 → WP1000
const W1_1 = WP0001 → WP0001
const W1_2 = WP0001 → WP0010
const W1_3 = WP0001 → WP0011
const W1_4 = WP0001 → WP0100
const W1_5 = WP0001 → WP0101
const W1_6 = WP0001 → WP0110
const W1_7 = WP0001 → WP0111
const W1_8 = WP0001 → WP1000
const W2_2 = WP0010 → WP0010
const W2_3 = WP0010 → WP0011
const W2_4 = WP0010 → WP0100
const W2_5 = WP0010 → WP0101
const W2_6 = WP0010 → WP0110
const W2_7 = WP0010 → WP0111
const W2_8 = WP0010 → WP1000
const W3_3 = WP0011 → WP0011
const W3_4 = WP0011 → WP0100
const W3_5 = WP0011 → WP0101
const W3_6 = WP0011 → WP0110
const W3_7 = WP0011 → WP0111
const W3_8 = WP0011 → WP1000
const W4_4 = WP0100 → WP0100
const W4_5 = WP0100 → WP0101
const W4_6 = WP0100 → WP0110
const W4_7 = WP0100 → WP0111
const W4_8 = WP0100 → WP1000
const W5_5 = WP0101 → WP0101
const W5_6 = WP0101 → WP0110
const W5_7 = WP0101 → WP0111
const W5_8 = WP0101 → WP1000
const W6_6 = WP0110 → WP0110
const W6_7 = WP0110 → WP0111
const W6_8 = WP0110 → WP1000
const W7_7 = WP0111 → WP0111
const W7_8 = WP0111 → WP1000
const W8_8 = WP1000 → WP1000
const WALL = W(ℝp)
const WNUL = W(∅)

const pos_warlpiri = [W0_0, W0_1, W0_2, W0_3, W0_4, W0_5, W0_6, W0_7, W0_8, W1_1, W1_2, W1_3, W1_4, W1_5, W1_6, W1_7, W1_8, W2_2, W2_3, W2_4, W2_5, W2_6, W2_7, W2_8, W3_3, W3_4, W3_5, W3_6, W3_7, W3_8, W4_4, W4_5, W4_6, W4_7, W4_8, W5_5, W5_6, W5_7, W5_8, W6_6, W6_7, W6_8, W7_7, W7_8, W8_8]

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
module AWS

using Compat: Compat, @something
using Base64
using Dates
using HTTP
using MbedTLS
using Mocking
using OrderedCollections: LittleDict, OrderedDict
using Retry
using Sockets
using URIs
using UUIDs: UUIDs
using XMLDict

export @service
export _merge
export AbstractAWSConfig, AWSConfig, AWSExceptions, AWSServices, Request
export ec2_instance_metadata, ec2_instance_region
export generate_service_url, global_aws_config, set_user_agent
export sign!, sign_aws2!, sign_aws4!
export JSONService, RestJSONService, RestXMLService, QueryService

const DEFAULT_REGION = "us-east-1"

include(joinpath("utilities", "utilities.jl"))

include("AWSExceptions.jl")
include("AWSCredentials.jl")
include("AWSConfig.jl")
include("AWSMetadata.jl")

include(joinpath("utilities", "request.jl"))
include(joinpath("utilities", "sign.jl"))

using ..AWSExceptions
using ..AWSExceptions: AWSException

const user_agent = Ref("AWS.jl/1.0.0")
const aws_config = Ref{AbstractAWSConfig}()

"""
    global_aws_config()

Retrieve the global AWS configuration.
If one is not set, create one with default configuration options.

# Keywords
- `kwargs...`: AWSConfig kwargs to be passed along if the global configuration is not already set

# Returns
- `AWSConfig`: The global AWS configuration
"""
function global_aws_config(; kwargs...)
    if !isassigned(aws_config) || !isempty(kwargs)
        aws_config[] = AWSConfig(; kwargs...)
    end

    return aws_config[]
end

"""
    global_aws_config(config::AbstractAWSConfig)

Set the global AWSConfig.

# Arguments
- `config::AWSConfig`: The AWSConfig to set in the global state

# Returns
- `AWSConfig`: Global AWSConfig
"""
function global_aws_config(config::AbstractAWSConfig)
    return aws_config[] = config
end

"""
    set_user_agent(new_user_agent::String)

Set the global user agent when making HTTP requests.

# Arguments
- `new_user_agent::String`: User agent to set when making HTTP requests

# Return
- `String`: The global user agent
"""
set_user_agent(new_user_agent::String) = return user_agent[] = new_user_agent

"""
    macro service(module_name::Symbol)

Include a high-level service wrapper based off of the module_name parameter.

When calling the macro you cannot match the predefined constant for the lowl level API.
The low level API constants are named in all lowercase, and spaces replaced with underscores.

Examples:
```julia
using AWS.AWSServices: secrets_manager
using AWS: @service

# This matches the constant and will error!
@service secrets_manager
> ERROR: cannot assign a value to variable AWSServices.secrets_manager from module Main

# This does NOT match the filename structure and will error!
@service secretsmanager
> ERROR: could not open file /.julia/dev/AWS.jl/src/services/secretsmanager.jl

# All of the examples below are valid!
@service Secrets_Manager
@service SECRETS_MANAGER
@service sECRETS_MANAGER
```

# Arguments
- `module_name::Symbol`: Name of the service to include high-level API wrappers in your namespace

# Return
- `Expression`: Base.include() call to introduce the high-level service API wrapper functions in your namespace
"""
macro service(module_name::Symbol)
    service_name = joinpath(@__DIR__, "services", lowercase(string(module_name)) * ".jl")

    return Expr(:toplevel,
    :(module($(esc(module_name)))
        Base.include($(esc(module_name)), $(esc(service_name)))
     end))
end

struct RestXMLService
    name::String
    api_version::String
end

struct QueryService
    name::String
    api_version::String
end

struct JSONService
    name::String
    api_version::String

    json_version::String
    target::String
end

struct RestJSONService
    name::String
    api_version::String

    service_specific_headers::LittleDict{String, String}
end

RestJSONService(name::String, api_version::String) = RestJSONService(name, api_version, LittleDict{String, String}())

# Needs to be included after the definition of struct otherwise it cannot find them
include("AWSServices.jl")

function generate_service_url(aws::AbstractAWSConfig, service::String, resource::String)
    SERVICE_HOST = "amazonaws.com"
    reg = region(aws)

    regionless_services = ("iam", "route53")

    if service in regionless_services || (service == "sdb" && reg == DEFAULT_REGION)
        reg = ""
    end

    return string("https://", service, ".", isempty(reg) ? "" : "$reg.", SERVICE_HOST, resource)
end

"""
    (service::RestXMLService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestXML request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::RestXMLService)(
    request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return_headers = _pop!(args, "return_headers", false)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        request_method=request_method,
        content=_pop!(args, "body", ""),
    )

    if request.service == "s3"
        request_uri = _clean_s3_uri(request_uri)
    end
    request.resource = _generate_rest_resource(request_uri, args)
    query_str = HTTP.escapeuri(args)

    if !isempty(query_str)
        if occursin('?', request.resource)
            request.resource *= "&$query_str"
        else
            request.resource *= "?$query_str"
        end
    end

    request.url = generate_service_url(aws_config, request.service, request.resource)

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::QueryService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a Query request to AWS.

# Arguments
- `operation::String`:
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::QueryService)(
    operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    POST_RESOURCE = "/"
    return_headers = _pop!(args, "return_headers", false)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        resource=POST_RESOURCE,
        request_method="POST",
        url=generate_service_url(aws_config, service.name, POST_RESOURCE),
    )

    request.headers["Content-Type"] = "application/x-www-form-urlencoded; charset=utf-8"

    args["Action"] = operation
    args["Version"] = service.api_version
    request.content = HTTP.escapeuri(_flatten_query(service.name, args))

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::JSONService)(
        operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a JSON request to AWS.

# Arguments
- `operation::String`: Name of the operation to perform
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::JSONService)(
    operation::String, args::AbstractDict{String, <:Any}=Dict{String, Any}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    POST_RESOURCE = "/"
    return_headers = _pop!(args, "return_headers", false)

    request = Request(;
        _extract_common_kw_args(service,args)...,
        resource=POST_RESOURCE,
        request_method="POST",
        content=json(args),
        url=generate_service_url(aws_config, service.name, POST_RESOURCE),
    )

    request.headers["Content-Type"] = "application/x-amz-json-$(service.json_version)"
    request.headers["X-Amz-Target"] = "$(service.target).$(operation)"

    return submit_request(aws_config, request; return_headers=return_headers)
end

"""
    (service::RestJSONService)(
        request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
        aws::AbstractAWSConfig=aws_config
    )

Perform a RestJSON request to AWS.

# Arguments
- `request_method::String`: RESTful request type, e.g.: `GET`, `HEAD`, `PUT`, etc.
- `request_uri::String`: AWS URI for the endpoint
- `args::AbstractDict{String, <:Any}`: Additional arguments to be included in the request

# Keywords
- `aws::AbstractAWSConfig`: AWSConfig containing credentials and other information for fulfilling the request, default value is the global configuration

# Returns
- `Tuple or Dict`: If `return_headers` is passed in through `args` a Tuple containing the Headers and Response will be returned, otherwise just a Dict
"""
function (service::RestJSONService)(
    request_method::String, request_uri::String, args::AbstractDict{String, <:Any}=Dict{String, String}();
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return_headers = _pop!(args, "return_headers", false)

    request = Request(;
        _extract_common_kw_args(service, args)...,
        request_method=request_method,
        resource=_generate_rest_resource(request_uri, args),
    )

    request.url = generate_service_url(aws_config, request.service, request.resource)

    if !isempty(service.service_specific_headers)
        merge!(request.headers, service.service_specific_headers)
    end

    request.headers["Content-Type"] = "application/json"
    request.content = json(args)

    return submit_request(aws_config, request; return_headers=return_headers)
end

end  # module AWS


text:
module CountMut

const data_path = "/home/guo/DeepSomaticSNV/data/CosmicCompleteExport.tsv"
#const data_path = "/home/guo/DeepSomaticSNV/data/1K.tsv"

### output format tsv
# gene_name mutation_id mutation_cds mutation_aa count

type MutStat
    gene_name::String
    id::String
    cds::String
    AA::String
    count::Int64
    samples::Array{String,1}
end

function check_mutationid(file_path=data_path)
    d = Dict{String,String}() # mutation_id => cds,AA
    open(file_path) do file
        i = 0
        while !eof(file)
            i += 1
            i % 1000_000 == 0 && info("processing $i lines")
            line = readline(file)
            row  = convert(Array{String,1}, split(chomp(line),"\t"))
            id = row[17] # mutation ID
            cds= row[18]
            AA = row[19]
            if !in(id, keys(d))
                d[id] = string(cds,AA)
            else
                @assert d[id] == string(cds,AA)
            end
        end
    end
end

function mutation_dict(file_path=data_path)
    d = Dict{String,String}() # mutation_id => mutation_id,cds,AA
    open(file_path) do file
        i = 0
        while !eof(file)
            i += 1
            i % 1000_000 == 0 && info("processing $i lines")
            line = readline(file)
            row  = convert(Array{String,1}, split(chomp(line),"\t"))
            id = row[17] # mutation ID
            cds= row[18]
            AA = row[19]
            if !in(id, keys(d))
                d[id] = string(id, ",", cds, ",", AA)
            else
                @assert d[id] == string(id, ",", cds, ",", AA)
            end
        end
    end
    d
end

function read_data(target_gene; file_path=data_path)
    mut_sample = Dict{String, Array{String,1}}() # id => sample_id
    i = 0
    open(file_path) do file
        header = convert(Array{String,1}, split(chomp(readline(file)), "\t"))
        while !eof(file)
            i += 1
            i % 1000_000 == 0 && info("processing $i lines")
            line = readline(file)
            row  = convert(Array{String,1}, split(chomp(line), "\t"))
            gene_name = row[1]
            if gene_name != target_gene
                continue
            end
            id = row[17] # mutation ID
            sample_id = row[6]
            if !in(id, keys(mut_sample))
                mut_sample[id] = [sample_id]
            else
                push!(mut_sample[id], sample_id)
            end
        end
    end
    mut_count = Dict{String,Int64}()
    for id in keys(mut_sample)
        mut_count[id] = length(mut_sample[id])
    end
    sort(collect(mut_count), by=x->x[2], rev=true)[1:20]
end

function main()
    #check_mutationid()
    mut_dict = mutation_dict(data_path)

    for gene in ["KRAS","NRAS","HRAS","BRAF","PIK3CA","EGFR","ERBB2","KIT","PDGFRA"]
        id_count = read_data(gene)
        open(string(gene, ".csv"), "w") do file
            header = "mutation_ID,mutation_CDS,mutation_AA,Count\n"
            write(file, header)
            for (id,count) in id_count
                line = string(mut_dict[id], ",", count, "\n")
                write(file, line)
            end
        end
    end

end

#=
# gene-name -- mut id --> sample id
function read_data(;file_path=data_path)
    gene_mutid = Dict{String, Array{String,1}}()
    target_genes = Set{String}(["KRAS","NRAS","HRAS","BRAF","PIK3CA","EGFR","ERBB2","KIT","PDGFRA"])
    i = 0
    open(file_path) do file
        header = convert(Array{String,1}, split(chomp(readline(file)), "\t"))
        #@show header
        while !eof(file)
            i += 1
            i % 1000_000 == 0 && info("processing $i lines")
            line = readline(file)
            row  = convert(Array{String,1}, split(chomp(line), "\t"))
            gene_name = row[1]
            if !in(gene_name, target_genes)
                continue
            end
            id = row[17] # mutation ID
            cds = row[18]
            AA = row[19]
            sample_id = row[6]
            @show sample_id

            break
        end
    end
    gene_mutid
end
=#

end


text:
module DifferentialGeometries

include("implementation.jl")

export AbstractCoordinateSystem, CoordinateSystem
export AbstractScalar, Scalar
export Covariant, Contravariant
export ToCartesian, FromCartesian
export CovariantVector, ContravariantVector
export ∇, ∂, grad, div, ∇o, curl, ∇x, J, dV, dl, gᵢⱼ, gⁱʲ, J⃗
export invert

end

max_stars_repo_name: jwscook/DifferentialGeometries


text:
# tests default choice model

@generator DCMGen begin # prefix DCM (for Default Choice Model) to avoid type name clashes
    start() = reps(a,1+1,2+2) # non-literal arguments to allow us to pass range to godelnumber during testing
    a() = choose(Int,-2-2,1+1)
    a() = choose(Float64,-2.0-2.0,1.0+1.0)
    a() = choose(Bool)
end

@generator DCMRepsGen begin
    start() = reps(a,1+1,2+2) # non-literal arguments to allow us to pass range to godelnumber during testing
    a() = 'a'
end

@generator DCMChooseBoolGen begin
    start() = choose(Bool) #
end

@generator DCMChooseIntGen begin
    start() = choose(Int,-2-2,1+1) # non-literal arguments to allow us to pass range to godelnumber during testing
end

@generator DCMChooseFloat64Gen begin
    start() = choose(Float64,-2.0-2.0,1.0+1.0) # non-literal arguments to allow us to pass range to godelnumber during testing
end

@generator DCMRuleGen begin
    start() = x()
    x() = 'a'
    x() = 'b'
    x() = 'c'
    x() = 'd'
end

# internally will use range of choice points, so convenient for testing choice model when more than one choice point
@generator DCMChooseStringGen begin
    start() = choose(String,"a(b|c)d+ef?")
end

@testset "default choice model" begin

	@testset "constructor" begin

	    gn = DCMGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
	    @test typeof(cm) == DataGenerators.SimpleChoiceModel

	end

	@testset "set/get parameters and ranges" begin

		gn = DCMGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)

		@testset "paramranges" begin
		    ranges = paramranges(cm)
		    @test typeof(ranges) <: Vector{Tuple{Float64,Float64}}
		    @test length(ranges) == 0
		end

		@testset "getparams" begin
		    params = getparams(cm)
		    @test typeof(params) <: Vector{Float64}
		    @test length(params) == 0
		end

		@testset "setparams" begin
		    newparams = [(paramrange[1]+paramrange[2])/2 for paramrange in paramranges(cm)]
		    setparams!(cm, newparams)
		    @test length(getparams(cm)) == 0 # can't check equality of params since some adjustment can be made by the cm
		end

	end

	@testset "rule choice point" begin

		gn = DCMRuleGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		# type ChoiceContext
		# 	derivationstate::DerivationState
		# 	cptype::Symbol
		# 	cpid::Uint64
		# 	datatype::DataType
		# 	lowerbound::Real
		# 	upperbound::Real
		# end
		cc = DataGenerators.ChoiceContext(DummyDerivationState(), :rule, cpids[1], Int, 1, 4)

		@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		    gnum, trace = DataGenerators.godelnumber(cm, cc)
		    @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		    @mtest_values_are [1,2,3,4] gnum
		end

	end

	@testset "sequence choice point" begin

		gn = DCMRepsGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "small finite range" begin

			cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 0, 2)

			@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
			    gnum, trace = DataGenerators.godelnumber(cm, cc)
			    @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
			    @test 0 <= gnum <= 2 # default choice model restricts sequence lengths to a maximum of 3 more than minimum
			    @mtest_values_are [0,1,2] gnum
			end

		end

		@testset "large finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 11, 16)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 11 <= gnum <= 13 # default choice model restricts sequence lengths to a maximum of 3 more than minimum
		        @mtest_values_are [11,12,13] gnum
		    end

		end

		@testset "semi-finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :sequence, cpids[1], Int, 1, typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 1 <= gnum <= 3 # default choice model restricts sequence lengths to a maximum of 3 more than minimum
		        @mtest_values_are [1,2,3] gnum
		    end

		end

	end

	@testset "Bool value choice point" begin

		gn = DCMChooseBoolGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Bool, false, true)

		@mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		    gnum, trace = DataGenerators.godelnumber(cm, cc)
		    @test convert(Bool,gnum) != nothing  # raises exception if value can't be converted
		    @mtest_values_are [false,true] gnum
		end

	end

	@testset "Int value choice point" begin

		gn = DCMChooseIntGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "small finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, -1, 2)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test -1 <= gnum <= 2
		        @mtest_values_are [-1,0,1,2] gnum
		    end

		end

		@testset "large finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, 11, 16)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 11 <= gnum <= 16
		        @mtest_values_include [11,13,16] gnum # just a selection of possible values including end points
		    end

		end

		@testset "semi-finite range (upper)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, 128, typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test 128 <= gnum <= typemax(Int)
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (lower)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, typemin(Int), 128)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test typemin(Int) <= gnum <= 128
		        @mtest_values_vary gnum
		    end

		end

		@testset "infinite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Int, typemin(Int), typemax(Int))

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Int,gnum) != nothing  # raises exception if value can't be converted
		        @test typemin(Int) <= gnum <= typemax(Int)
		        @mtest_values_vary gnum
		    end

		end

	end

	@testset "Float64 value choice point" begin

		gn = DCMChooseFloat64Gen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)
		cpi = choicepointinfo(gn)
		cpids = collect(keys(cpi))

		@testset "finite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -42.2, -8.7)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test -42.2 <= gnum <= -8.7
		        @mtest_that_sometimes round(gnum) != gnum
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (upper)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, 450001.6, Inf)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test 450001.6 <= gnum
		        @mtest_values_vary gnum
		    end

		end

		@testset "semi-finite range (lower)" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -Inf, 450001.6)

		    @mtestset "valid Godel numbers returned" begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @test gnum <= 450001.6
		        @mtest_values_vary gnum
		    end

		end

		@testset "infinite range" begin

		    cc = DataGenerators.ChoiceContext(DummyDerivationState(), :value, cpids[1], Float64, -Inf, Inf)

		    @mtestset "valid Godel numbers returned" reps=Main.REPS alpha=Main.ALPHA begin
		        gnum, trace = DataGenerators.godelnumber(cm, cc)
		        @test convert(Float64,gnum) != nothing  # raises exception if value can't be converted
		        @mtest_values_vary gnum
		    end

		end

	end

	@testset "generate for model with multiple choice points" begin

		gn = DCMChooseStringGen()
		setchoicemodel!(gn, SimpleChoiceModel())
	    cm = choicemodel(gn)

		@mtestset "full range of values generated using choice model" reps=Main.REPS alpha=Main.ALPHA begin
		    td = choose(gn, choicemodel=cm)
		    @test ismatch(r"^a(b|c)d+ef?$", td)
		    @mtest_values_include [1,2,3] count(x->x=='d', td)
		    @mtest_values_are ['b','c'] td[2]
		    @mtest_values_are ['e','f'] td[end]
		end

	end

end


text:
"""
    NeuralNetworkKernel()

Kernel of a Gaussian process obtained as the limit of a Bayesian neural network with a
single hidden layer as the number of units goes to infinity.

# Definition

Consider the single-layer Bayesian neural network
``f \\colon \\mathbb{R}^d \\to \\mathbb{R}`` with ``h`` hidden units defined by
```math
f(x; b, v, u) = b + \\sqrt{\\frac{\\pi}{2}} \\sum_{i=1}^{h} v_i \\mathrm{erf}\\big(u_i^\\top x\\big),
```
where ``\\mathrm{erf}`` is the error function, and with prior distributions
```math
\\begin{aligned}
b &\\sim \\mathcal{N}(0, \\sigma_b^2),\\\\
v &\\sim \\mathcal{N}(0, \\sigma_v^2 \\mathrm{I}_{h}/h),\\\\
u_i &\\sim \\mathcal{N}(0, \\mathrm{I}_{d}/2) \\qquad (i = 1,\\ldots,h).
\\end{aligned}
```
As ``h \\to \\infty``, the neural network converges to the Gaussian process
```math
g(\\cdot) \\sim \\mathcal{GP}\\big(0, \\sigma_b^2 + \\sigma_v^2 k(\\cdot, \\cdot)\\big),
```
where the neural network kernel ``k`` is given by
```math
k(x, x') = \\arcsin\\left(\\frac{x^\\top x'}{\\sqrt{\\big(1 + \\|x\\|^2_2\\big) \\big(1 + \\|x'\\|_2^2\\big)}}\\right)
```
for inputs ``x, x' \\in \\mathbb{R}^d``.[^CW]

[^CW]: C. K. I. Williams (1998). Computation with infinite neural networks.
"""
struct NeuralNetworkKernel <: Kernel end

function (κ::NeuralNetworkKernel)(x, y)
    return asin(dot(x, y) / sqrt((1 + sum(abs2, x)) * (1 + sum(abs2, y))))
end

function kernelmatrix(::NeuralNetworkKernel, x::ColVecs, y::ColVecs)
    validate_inputs(x, y)
    X_2 = sum(x.X .* x.X; dims=1)
    Y_2 = sum(y.X .* y.X; dims=1)
    XY = x.X' * y.X
    return asin.(XY ./ sqrt.((X_2 .+ 1)' * (Y_2 .+ 1)))
end

function kernelmatrix(::NeuralNetworkKernel, x::ColVecs)
    X_2_1 = sum(x.X .* x.X; dims=1) .+ 1
    XX = x.X' * x.X
    return asin.(XX ./ sqrt.(X_2_1' * X_2_1))
end

function kernelmatrix_diag(::NeuralNetworkKernel, x::ColVecs)
    x_2 = vec(sum(x.X .* x.X; dims=1))
    return asin.(x_2 ./ (x_2 .+ 1))
end

function kernelmatrix_diag(::NeuralNetworkKernel, x::ColVecs, y::ColVecs)
    validate_inputs(x, y)
    x_2 = vec(sum(x.X .* x.X; dims=1) .+ 1)
    y_2 = vec(sum(y.X .* y.X; dims=1) .+ 1)
    xy = vec(sum(x.X' .* y.X'; dims=2))
    return asin.(xy ./ sqrt.(x_2 .* y_2))
end

function kernelmatrix(::NeuralNetworkKernel, x::RowVecs, y::RowVecs)
    validate_inputs(x, y)
    X_2 = sum(x.X .* x.X; dims=2)
    Y_2 = sum(y.X .* y.X; dims=2)
    XY = x.X * y.X'
    return asin.(XY ./ sqrt.((X_2 .+ 1) * (Y_2 .+ 1)'))
end

function kernelmatrix(::NeuralNetworkKernel, x::RowVecs)
    X_2_1 = sum(x.X .* x.X; dims=2) .+ 1
    XX = x.X * x.X'
    return asin.(XX ./ sqrt.(X_2_1 * X_2_1'))
end

function kernelmatrix_diag(::NeuralNetworkKernel, x::RowVecs)
    x_2 = vec(sum(x.X .* x.X; dims=2))
    return asin.(x_2 ./ (x_2 .+ 1))
end

function kernelmatrix_diag(::NeuralNetworkKernel, x::RowVecs, y::RowVecs)
    validate_inputs(x, y)
    x_2 = vec(sum(x.X .* x.X; dims=2) .+ 1)
    y_2 = vec(sum(y.X .* y.X; dims=2) .+ 1)
    xy = vec(sum(x.X .* y.X; dims=2))
    return asin.(xy ./ sqrt.(x_2 .* y_2))
end

Base.show(io::IO, ::NeuralNetworkKernel) = print(io, "Neural Network Kernel")