text:
export
    AbstractDisturbance,
    DisturbanceSequence

abstract type AbstractDisturbance end
const DisturbanceSequence = Vector{Pair{Int,AbstractDisturbance}}

export
    StochasticProblem,
    stochastic_problem

"""
    StochasticProblem{P<:AbstractPC_TAPF}

Defines a stochastic version of PC_TAPF, wherein different disturbances can
cause unexpected problems in the factory.
"""
struct StochasticProblem{P<:AbstractPC_TAPF}
    prob::P
    disturbances::DisturbanceSequence
end
CRCBS.get_env(spc_tapf::StochasticProblem) = get_env(spc_tapf.prob)

function stochastic_problem(ptype::Type{P},solver,prob,clean_up_bot_ICS,disturbances) where {P<:AbstractPC_MAPF}
    env = get_env(deepcopy(prob))
    env_graph = get_graph(env)
    sched = get_schedule(env)
    cache = get_cache(env)
    for pred in clean_up_bot_ICS
        add_new_robot_to_schedule!(sched,pred,get_problem_spec(env))
    end
    robot_ICs = get_robot_ICs(sched)
    prob_spec = get_problem_spec(env)
    # populate_environment_dict_layers!(env.env_layers,sched,env_graph,prob_spec)
    new_prob = P(prob,construct_search_env(solver,sched,prob_spec,env_graph,cache))
    StochasticProblem(new_prob,disturbances)
end

export OilSpill
"""
    OilSpill

An obstruction that affects vertices `vtxs` and edges `edges`.
"""
struct OilSpill <: AbstractDisturbance
    vtxs::Set{Int}
    edges::Set{Edge}
end
function OilSpill(G::AbstractGraph,vtxs::Vector{Int})
    OilSpill(Set(vtxs),edge_cover(G,vtxs))
end

export Intruder
"""
    Intruder

An intruder that begins at location `start_vtx` and follows policy `policy`
"""
struct Intruder{P} <: AbstractDisturbance
    policy::P
    start_vtx::Int # start vtx
end

export DeadRobot
"""
    DeadRobot

Robot `id` is frozen.

Effect:
- Freeze robot
- Add "no-go" constraint to CBS/PIBT (how to do consistently? Perhaps place in
    SearchEnv and add directly to PCCBSEnv) OR temporarily remove vertex from
    graph
- Set robot state to NULL state? How to avoid having CBS complain about
    conflicts? Maybe set State to NULL State and place DeadRobotObject at the
    collection site?
- Dispatch CleanUpBot to collect frozen robot
- When CleanUpBot returns to "garage", regenerate frozen Robot's ROBOT_AT node
    and valid state.
"""
struct DeadRobot <: AbstractDisturbance
    id::Int
end

export DelayedRobot
"""
    DelayedRobot

Robot `id` is delayed by `dt` timesteps.
"""
struct DelayedRobot <: AbstractDisturbance
    id::Int
    dt::Int
end

export DroppedObject
"""
    DroppedObject

Object `id` dropped.
"""
struct DroppedObject <: AbstractDisturbance
    id::Int
end

export handle_disturbance!

function remove_vtxs(sched,remove_set::Set{A}) where {A<:AbstractID}
    remove_vtxs(sched,Set{Int}(get_vtx(sched,id) for id in remove_set))
end
function remove_vtxs!(sched,remove_set)
    new_sched = sched
    rem_nodes!(new_sched,remove_set)
    set_leaf_operation_vtxs!(new_sched)
    process_schedule!(new_sched)
    new_sched
end
function remove_vtxs(sched,remove_set)
    # Construct new graph
    # new_sched = sched
    # rem_nodes!(new_sched,remove_set)
    # set_leaf_operation_vtxs!(new_sched)
    # process_schedule!(new_sched)

    new_sched = OperatingSchedule()
    keep_vtxs = setdiff(Set{Int}(collect(vertices(sched))), remove_set)
    # add all non-deleted nodes to new project schedule
    for v in keep_vtxs
        # node_id = get_vtx_id(sched,v)
        # node = get_node_from_id(sched, node_id)
        # path_spec = get_path_spec(sched,v)
        # add_to_schedule!(new_sched,path_spec,node,node_id)
        add_to_schedule!(new_sched,get_node(sched,v))
    end
    # add all edges between nodes that still exist
    for e in edges(get_graph(sched))
        add_edge!(new_sched, get_vtx_id(sched, e.src), get_vtx_id(sched, e.dst))
    end
    set_leaf_operation_vtxs!(new_sched)
    process_schedule!(new_sched)
    # @assert sanity_check(new_sched)
    new_sched
end

"""
    get_delivery_task_vtxs(sched::OperatingSchedule,o::ObjectID)

Return all vertices that correspond to the delivery task
(`COLLECT` → `CARRY` → `DEPOSIT`) of object `o`.
"""
function get_delivery_task_vtxs(sched::OperatingSchedule,o::ObjectID)
    v = get_vtx(sched,o)
    vtxs = collect(capture_connected_nodes(sched,v,
        v->check_object_id(get_node_from_vtx(sched,v),o)))
    setdiff!(vtxs,v)
    return vtxs
end
"""
    isolate_delivery_task_vtxs(sched,o,vtxs=get_delivery_task_vtxs(sched,o))

Returns:
- incoming: a set of all incoming Action ScheduleNodes
- outgoing: a set of all outgoing Action ScheduleNodes
- op: the target Operation
"""
function isolate_delivery_task_vtxs(sched::OperatingSchedule,o::ObjectID,
        vtxs=get_delivery_task_vtxs(sched,o)
    )
    node_ids = Set{AbstractID}(get_vtx_id(sched,v) for v in vtxs)
    # incoming and outgoing robot nodes -- for single and team robot tasks
    in_edges = exclusive_edge_cover(sched,vtxs,:in)
    incoming = filter(node->matches_template(AbstractRobotAction,node),
        map(e->get_node(sched,e.src),collect(in_edges)))
    out_edges = exclusive_edge_cover(sched,vtxs,:out)
    outgoing = filter(node->matches_template(AbstractRobotAction,node),
        map(e->get_node(sched,e.dst),collect(out_edges)))
    # target operation node
    op = filter(node->matches_template(Operation,node),
        map(e->get_node(sched,e.dst),collect(out_edges)))[1]
    return incoming, outgoing, op
end

"""
    stitch_disjoint_node_sets!(sched,incoming,outgoing)

Finds and adds the appropriate edges between two sets of nodes. It is assumed
that `size(incoming) == size(outgoing)`, that each node in `incoming` has
exactly one feasible successor in `outgoing`, and that each node in `outgoing`
has exactly one feasible predecessor in `incoming`.
"""
function stitch_disjoint_node_sets!(sched,incoming,outgoing,env_state)
    @assert length(incoming) == length(outgoing)
    out_idxs = collect(1:length(outgoing))
    for node1 in incoming
        r = get_robot_id(node1)
        for i in 1:length(outgoing)
            idx = out_idxs[i]
            node2 = outgoing[idx]
            node2 = outgoing[i]
            if get_robot_id(node2) == r
                x = get_location_id(robot_positions(env_state)[r])
                add_edge!(sched,node1,node2)
                n1 = BOT_GO(r,get_initial_location_id(node1),x)
                n2 = BOT_GO(r,x,get_destination_location_id(node2))
                @log_info(-1,0,"connecting ",string(n1)," → ",string(n2))
                replace_in_schedule!(sched,n1,node1.id)
                replace_in_schedule!(sched,n2,node2.id)
                deleteat!(out_idxs,i)
                # deleteat!(outgoing,i)
                break
            end
        end
    end
    @assert length(out_idxs) == 0
    # @assert length(outgoing) == 0
    return sched
end

"""
    handle_disturbance!(solver,prob,env,d::DroppedObject,t,env_state=get_env_state(env,t))

Returns a new `SearchEnv` with a modified `OperatingSchedule`. The new schedule
replaces the previous delivery task (`OBJECT_AT(o,old_x)` → `COLLECT` → `CARRY`
→ `DEPOSIT`) with a new `CleanUpBot` delivery task (`OBJECT_AT(o,new_x)` →
`CUB_COLLECT` → `CUB_CARRY` → `CUB_DEPOSIT`).
It is assumed that the time `t` corresponds to a moment when the object referred
to by `d.id` is being `CARRIED`.
"""
function handle_disturbance!(solver,prob,env::SearchEnv,d::DroppedObject,t,
        env_state = get_env_state(env,t),
        )
    # Schedule surgery
    # - remove Delivery Task: OBJECT_AT(o,old_x) -> COLLECT -> CARRY -> DEPOSIT
    # - replace with equivalent CleanUpBot Delivery Task: OBJECT_AT(o,new_x) -> BOT_COLLECT{CleanUpBot} -> CARRY -> DEPOSIT
    sched = get_schedule(env)
    o = ObjectID(d.id)
    vtxs = get_delivery_task_vtxs(sched,o)
    # incoming and outgoing robot nodes -- for single and team robot tasks
    incoming, outgoing, op = isolate_delivery_task_vtxs(sched,o,vtxs)
    # Add GO->GO edges for affected robot(s), as if they were never assigned
    stitch_disjoint_node_sets!(sched,incoming,outgoing,env_state)
    # remove old nodes
    new_sched = remove_vtxs(sched,vtxs)
    # new_sched = remove_vtxs!(sched,vtxs)
    # add new CleanUpBot task nodes
    x = get_location_ids(object_position(env_state,o))
    replace_in_schedule!(new_sched,OBJECT_AT(o,x),o)
    set_t0!(new_sched,o,t)
    add_headless_delivery_task!(
        new_sched,get_problem_spec(env,:CleanUpBot),o,op.id,CleanUpBot;
        t0=t
    )
    process_schedule!(new_sched)
    construct_search_env(
        solver,
        new_sched,
        env,
        initialize_planning_cache(new_sched)
    )
end

function apply_disturbance!(env_graph::GridFactoryEnvironment,d::OilSpill)
    remove_edges!(env_graph,d.edges)
    return env_graph
end
function remove_disturbance!(env_graph::GridFactoryEnvironment,d::OilSpill)
    add_edges!(env_graph,d.edges)
    return env_graph
end

"""
    regenerate_path_specs!(solver,env)

Recompute all paths specs (to account for changes in the env_graphs that will
be propagated to the ProblemSpec's distance function as well.)
"""
function regenerate_path_specs!(solver,env)
    sched = get_schedule(env)
    for v in vertices(sched)
        node = get_node_from_vtx(sched,v)
        prob_spec = get_problem_spec(env,graph_key(node))
        set_path_spec!(sched,v,generate_path_spec(sched,prob_spec,node))
    end
    update_planning_cache!(solver,env)
    return env
end

function add_headless_cleanup_task!(
        sched::OperatingSchedule,
        spec::ProblemSpec,
        d::OilSpill,
        return_vtx::LocationID=LocationID(-1)
        ;
        t0=0,
        )

    robot_id = CleanUpBotID(-1)
    action_id = get_unique_action_id()
    add_to_sched!(sched, spec, CLEAN_UP(robot_id, d.vtxs), action_id)
    set_t0!(sched,action_id,t0)

    prev_action_id = action_id
    action_id = get_unique_action_id()
    add_to_sched!(sched, spec, CUB_GO(robot_id, d.vtxs[1], return_vtx), action_id)
    set_t0!(sched,action_id,t0)
    add_edge!(sched, prev_action_id, action_id)

    return
end

function remove_robot!(env::SearchEnv,id::BotID,t::Int)
    sched = get_schedule(env)
    G = get_graph(sched)
    # schedule
    to_remove = Set{AbstractID}(id)
    v = get_vtx(sched,id)
    for vp in edges(bfs_tree(G,v))
        node_id = get_vtx_id(sched,vp)
        if isa(node_id,ActionID)
            node = get_node_from_id(sched,node_id)
            if isa(node,BOT_GO)
                push!(to_remove,node_id)
            else
                new_node = replace_robot_id(node,-1)
                replace_in_schedule!(sched,get_problem_spec(env),new_node,node_id)
            end
        end
    end
    for node_id in to_remove
        delete_node!(sched,node_id)
    end
    # Verify that the robot is no longer in schedule
    for v in vertices(G)
        node_id = get_vtx_id(sched,v)
        if isa(node_id,ActionID)
            node = get_node_from_id(sched,node_id)
            @assert !(get_id(id) in get_valid_robot_ids(node)) "Robot id $(string(id)) should be wiped from schedule, but is present in $(string(node))"
        end
    end
    # Remap other robots' IDs? TODO refactor OperatinSchedule and RoutePlan so that I don't have to remap ids? i.e., have both uids and idxs?
    # - remove BOT_AT node
    # set assignment ids to -1 with replace_robot_id()
    # if in the middle of CARRY, another robot needs to get the object
end


text:
module Test04

import Automa
import Automa.RegExp: @re_str
using Test

@testset "Test04" begin
    re = Automa.RegExp
    beg_a = re.cat('a', re"[ab]*")
    end_b = re.cat(re"[ab]*", 'b')
    beg_a_end_b = re.isec(beg_a, end_b)

    machine = Automa.compile(beg_a_end_b)

    for generator in (:table, :inline, :goto), checkbounds in (true, false), clean in (true, false)
        ctx = Automa.CodeGenContext(generator=generator, checkbounds=checkbounds, clean=clean)
        init_code = Automa.generate_init_code(ctx, machine)
        exec_code = Automa.generate_exec_code(ctx, machine)
        validate = @eval function (data)
            $(init_code)
            p_end = p_eof = lastindex(data)
            $(exec_code)
            return cs == 0
        end
        @test validate(b"") == false
        @test validate(b"a") == false
        @test validate(b"aab") == true
        @test validate(b"ab") == true
        @test validate(b"aba") == false
        @test validate(b"abab") == true
        @test validate(b"abb") == true
        @test validate(b"abbb") == true
        @test validate(b"b") == false
        @test validate(b"bab") == false
    end
end

end

ext: jl


text:
function settle{T<:Blotter}(b::Vector{T})
    ###  array of tickers
    ###  syms  = Symbol[s.value.asset.ticker for s in b]
    ###  iterate across this array with call to method below
end

function settle{T<:Blotter}(b::Vector{T}, ticker::Symbol)
    # index out by ticker

    bt = b[ticker]

    # aggregate on asset class and sort

    # stocks (Stock, LongStock, ShortStock)
    vs   = [Stock, LongStock, ShortStock]
    stck = bt[vs]

    # call options (LongCall, ShortCall),
    vc  = [LongCall, ShortCall]
    cll = bt[vc]

    # put options (LongPut, ShortPut)
    vp  = [LongPut, ShortPut]
    pt  = bt[vp]

    # futures (LongFuture, ShortFuture)
    vf  = [LongFuture, ShortFuture]
    ft  = bt[vf]

    # return settled blotter and un-settled blotter
    ### sb  = PnL(opened, closed, quantity, commission, pnl, financial asset)
    ### usb = Blotter(x[2].timestamp, AssetTransaction(x[2].value.quantity - 10, x[2].value.price, x[2].value.asset))

    return sb, usb
end
## function settle{T<:(Blotter.Stock, LongStock, ShortStock)}(b::Vector{T}
## end


text:
using RecSys
using Blobs
using MAT
#include("/home/tan/Work/julia/packages/Blobs/examples/matrix.jl")
using RecSys.MatrixBlobs

function split_sparse(S, chunkmax, metadir)
    isdir(metadir) || mkdir(metadir)
    spblobs = SparseMatBlobs(Float64, Int64, metadir)
    chunknum = 1
    count = 1
    colstart = 1
    nzval = S.nzval
    rowval = S.rowval
    for col in 1:size(S,2)
        npos = S.colptr[col+1]
        if (npos >= (count + chunkmax)) || (col == size(S,2))
            print("\tchunk $chunknum ... ")
            append!(spblobs, S[:, colstart:col])
            colstart = col+1
            count = npos
            chunknum += 1
            println("done")
        end
    end
    Blobs.save(spblobs)
    nothing
end

function splitall(inp::DlmFile, output_path::AbstractString, nsplits::Int)
    println("reading inputs...")
    ratings = RecSys.read_input(inp)

    users   = convert(Vector{Int64},   ratings[:,1]);
    items   = convert(Vector{Int64},   ratings[:,2]);
    ratings = convert(Vector{Float64}, ratings[:,3]);
    R = sparse(users, items, ratings);
    splitall(R, output_path, nsplits)
end

function splitall(R::SparseMatrixCSC, output_path::AbstractString, nsplits::Int)
    nratings = nnz(R)
    println("$nratings ratings in $(size(R)) sized sparse matrix")
    println("removing empty...")
    R, non_empty_items, non_empty_users = RecSys.filter_empty(R)
    nratings = nnz(R)
    println("$nratings ratings in $(size(R)) sized sparse matrix")

    nsplits_u = round(Int, nratings/nsplits)
    nsplits_i = round(Int, nratings/nsplits)

    println("splitting R itemwise at $nsplits_i items...")
    split_sparse(R, nsplits_i, joinpath(output_path, "R_itemwise"))
    RT = R'
    println("splitting RT userwise at $nsplits_u users...")
    split_sparse(RT, nsplits_u, joinpath(output_path, "RT_userwise"))
    nothing
end

function split_netflix(dataset_path = "/data/Work/datasets/netflix-full")
    println("reading dataset")
    f = matopen(joinpath(dataset_path, "netflixAll.mat"))
    a = read(f, "R_with_probe")
    println("randomizing items to remove skew")
    a = a[:, randperm(size(a,2))]
    splitall(a, joinpath(dataset_path, "splits"), 20)
end

function split_movielens(dataset_path = "/data/Work/datasets/movielens/ml-20m")
    ratings_file = DlmFile(joinpath(dataset_path, "ratings.csv"); dlm=',', header=true)
    splitall(ratings_file, joinpath(dataset_path, "splits"), 10)
end

function split_lastfm(dataset_path = "/data/Work/datasets/last_fm_music_recommendation/profiledata_06-May-2005")

    function read_artist_map(artist_map::FileSpec)
        t1 = time()
        RecSys.@logmsg("reading artist map")
        A = read_input(artist_map)
        valid = map(x->isa(x, Integer), A)
        valid = valid[:,1] & valid[:,2]
        Avalid = convert(Matrix{Int64}, A[valid, :])

        amap = Dict{Int64,Int64}()
        for idx in 1:size(Avalid,1)
            bad_id = Avalid[idx, 1]
            good_id = Avalid[idx, 2]
            amap[bad_id] = good_id
        end
        RecSys.@logmsg("read artist map in $(time()-t1) secs")
        amap
    end

    function read_trainingset(trainingset::FileSpec, amap::Dict{Int64,Int64})
        t1 = time()
        RecSys.@logmsg("reading trainingset")
        T = read_input(trainingset)
        for idx in 1:size(T,1)
            artist_id = T[idx,2]
            if artist_id in keys(amap)
                T[idx,2] = amap[artist_id]
            end
        end
        users   = convert(Vector{Int64},   T[:,1])
        artists = convert(Vector{Int64},   T[:,2])
        ratings = convert(Vector{Float64}, T[:,3])
        S = sparse(users, artists, ratings)
        RecSys.@logmsg("read trainingset in $(time()-t1) secs")
        S
    end

    amap = read_artist_map(DlmFile(joinpath(dataset_path, "artist_alias.txt")))
    T = read_trainingset(DlmFile(joinpath(dataset_path, "user_artist_data.txt")), amap)

    println("randomizing items to remove skew")
    T = T[:, randperm(size(T,2))]

    splitall(T, joinpath(dataset_path, "splits"), 20)
end

function load_splits(dataset_path = "/data/Work/datasets/last_fm_music_recommendation/profiledata_06-May-2005/splits")
    sp = SparseMatBlobs(joinpath(dataset_path, "R_itemwise"))
    for idx in 1:length(sp.splits)
        p = sp.splits[idx]
        r = p.first
        part, _r = Blobs.load(sp, first(r))
        RecSys.@logmsg("got part of size: $(size(part)), with r: $r, _r:$_r")
    end
    println("finished")
end

##
# an easy way to generate somewhat relevant test data is to take an existing dataset and replicate
# items and users based on existing data.
function Blobs.append!{Tv,Ti}(sp::SparseMatBlobs{Tv,Ti}, blob::Blob)
    S = blob.data.value
    m,n = size(S)
    if isempty(sp.splits)
        sp.sz = (m, n)
        idxrange = 1:n
    else
        (sp.sz[1] == m) || throw(BoundsError("SparseMatBlobs $(sp.sz)", (m,n)))
        old_n = sp.sz[2]
        idxrange = (old_n+1):(old_n+n)
        sp.sz = (m, old_n+n)
    end

    push!(sp.splits, idxrange => blob.id)
    RecSys.@logmsg("appending blob $(blob.id) of size: $(size(S)) for idxrange: $idxrange, sersz: $(blob.metadata.size)")
    blob
end

function generate_test_data(setname::AbstractString, generated_data_path, original_data_path, mul_factor)
    sp1 = SparseMatBlobs(joinpath(original_data_path, setname))
    metapath = joinpath(generated_data_path, setname)
    sp2 = SparseMatBlobs(Float64, Int64, metapath)
    isdir(metapath) || mkdir(metapath)
    for idx in 1:length(sp1.splits)
        p = sp1.splits[idx]
        r = p.first
        part, _r = Blobs.load(sp1, first(r))
        RecSys.@logmsg("got part of size: $(size(part)), with r: $r, _r:$_r")
        part_out = part
        for x in 1:(mul_factor-1)
            part_out = vcat(part_out, part)
        end
        blob = append!(sp2, part_out)
        for x in 1:(mul_factor-1)
            append!(sp2, blob)
            RecSys.@logmsg("generated part of size: $(size(part_out))")
        end
    end
    Blobs.save(sp2)
end

function generate_test_data(generated_data_path = "/data/Work/datasets/last_fm_music_recommendation/profiledata_06-May-2005/splits2",
        original_data_path = "/data/Work/datasets/last_fm_music_recommendation/profiledata_06-May-2005/splits",
        mul_factor=2)
    generate_test_data("R_itemwise", generated_data_path, original_data_path, mul_factor)
    generate_test_data("RT_userwise", generated_data_path, original_data_path, mul_factor)
    println("finished")
end


text:
######################################
# Example of kpik function
# Based on example_kpik.m on V. Simoncini's Website
######################################
workspace()
include("../src/LargeMatrixEquations.jl")
using LargeMatrixEquations

nh = 20
n = nh^2
T = diagm(ones(nh)*2)-diagm(ones(nh-1),-1)-diagm(ones(nh-1),1)
I = eye(nh,nh)
A = sparse(-(kron(T,I)+kron(I,T)))

srand(1234)
B = randn(n,2)
E = eye(n)#sparse(diagm(rand(n),0))

m=100
tol=1e-9
tolY=1e-12
Z,er2=kpik(A,B,tolY=tolY)

println(norm(A*Z*Z'*E+E*Z*Z'*A'+B*B'))

max_stars_repo_name: garrettthomaskth/LargeMatrixEquations

max_issues_repo_name: garrettthomaskth/LargeMatrixEquations

max_forks_repo_name: garrettthomaskth/LargeMatrixEquations


text:
@testset "core" begin
    @test uconvert(u"ppm", 1.0u"Hz"/1.0u"MHz") == 1.0u"ppm"
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
using JSON

include("network.jl")

"""
    works(; kwargs...)

Arguments:
- `doi::Union{Nothing, String}`: DOI of the work.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- Dict
"""
function works(; doi::Union{Nothing, String}=nothing,
               query::Union{Nothing, String}=nothing,
               filter::Union{Nothing, Dict{String, String}}=nothing,
               offset::Union{Nothing, Int}=nothing,
               limit::Union{Nothing, Int}=nothing,
               sample::Union{Nothing, Int}=nothing,
               sort::Union{Nothing, String}=nothing,
               order::Union{Nothing, String}=nothing,
               facet::Union{Nothing, Dict{String, Int}}=nothing,
               select::Union{Nothing, String}=nothing,
               cursor::Union{Nothing, String}=nothing)

    url = "https://api.crossref.org/works"

    if !isnothing(doi)
            url = string(url, "/", doi)
            r = request(url=url)
            return JSON.parse(String(r.body))
    end

    r = request(url=url,
                query=query,
                filter=filter,
                offset=offset,
                limit=limit,
                sample=sample,
                sort=sort,
                order=order,
                facet=facet,
                select=select,
                cursor=cursor)

    return JSON.parse(String(r.body))

end

"""
    members(; kwargs...)

Arguments:
- `member_id::Union{Nothing, Int}`: ID of a Crossref member.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `works::Bool`: If true, works returned as well.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- Dict
"""
function members(; member_id::Union{Nothing, Int}=nothing,
                 query::Union{Nothing, String}=nothing,
                 filter::Union{Nothing, Dict{String, String}}=nothing,
                 offset::Union{Nothing, Int}=nothing,
                 limit::Union{Nothing, Int}=nothing,
                 sample::Union{Nothing, Int}=nothing,
                 sort::Union{Nothing, String}=nothing,
                 order::Union{Nothing, String}=nothing,
                 facet::Union{Nothing, Dict{String, Int}}=nothing,
                 select::Union{Nothing, String}=nothing,
                 works::Bool=false,
                 cursor::Union{Nothing, String}=nothing)

        url="https://api.crossref.org/members"

        if !isnothing(member_id)
                url = string(url, "/", member_id)
                if works
                        url = string(url, "/works")
                else
                        r = request(url=url)
                        return JSON.parse(String(r.body))
                end

        end

        r = request(url=url,
                    query=query,
                    filter=filter,
                    offset=offset,
                    limit=limit,
                    sample=sample,
                    sort=sort,
                    order=order,
                    facet=facet,
                    select=select,
                    cursor=cursor)

        return JSON.parse(String(r.body))

end

"""
    funders(; kwargs...)

Arguments:
- `funder_id::Union{Nothing, String}`: ID of a funder.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `works::Bool`: If true, works returned as well.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- `Dict`
"""
function funders(; funder_id::Union{Nothing, String}=nothing,
                 query::Union{Nothing, String}=nothing,
                 filter::Union{Nothing, Dict{String, String}}=nothing,
                 offset::Union{Nothing, Int}=nothing,
                 limit::Union{Nothing, Int}=nothing,
                 sample::Union{Nothing, Int}=nothing,
                 sort::Union{Nothing, String}=nothing,
                 order::Union{Nothing, String}=nothing,
                 facet::Union{Nothing, Dict{String, Int}}=nothing,
                 select::Union{Nothing, String}=nothing,
                 works::Bool=false,
                 cursor::Union{Nothing, String}=nothing)

        url="https://api.crossref.org/funders"

        if !isnothing(funder_id)
                url = string(url, "/", funder_id)
                if works
                        url = string(url, "/works")
                else
                        r = request(url=url)
                        return JSON.parse(String(r.body))
                end

        end

        r =  request(url=url,
                     query=query,
                     filter=filter,
                     offset=offset,
                     limit=limit,
                     sample=sample,
                     sort=sort,
                     order=order,
                     facet=facet,
                     select=select,
                     cursor=cursor)

        return JSON.parse(String(r.body))

end

"""
    journals(; kwargs...)

Arguments:
- `issn::Union{Nothing, String}`: A issn of a journal.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `works::Bool`: If true, works returned as well.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- `Dict`
"""
function journals(; issn::Union{Nothing, String}=nothing,
                  query::Union{Nothing, String}=nothing,
                  filter::Union{Nothing, Dict{String, String}}=nothing,
                  offset::Union{Nothing, Int}=nothing,
                  limit::Union{Nothing, Int}=nothing,
                  sample::Union{Nothing, Int}=nothing,
                  sort::Union{Nothing, String}=nothing,
                  order::Union{Nothing, String}=nothing,
                  facet::Union{Nothing, Dict{String, Int}}=nothing,
                  select::Union{Nothing, String}=nothing,
                  works::Bool=false,
                  cursor::Union{Nothing, String}=nothing)

        url="https://api.crossref.org/journals"

        if !isnothing(issn)
                url = string(url, "/", issn)
                if works
                        url = string(url, "/works")
                else
                        r = request(url=url)
                        return JSON.parse(String(r.body))
                end

        end

        r = request(url=url,
                    query=query,
                    filter=filter,
                    offset=offset,
                    limit=limit,
                    sample=sample,
                    sort=sort,
                    order=order,
                    facet=facet,
                    select=select,
                    cursor=cursor)

        return JSON.parse(String(r.body))

end

"""
    types(; kwargs...)

Arguments:
- `type_id::Union{Nothing, String}`: ID of a metadata work type.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `works::Bool`: If true, works returned as well.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- `Dict`
"""
function types(; type_id::Union{Nothing, String}=nothing,
               query::Union{Nothing, String}=nothing,
               filter::Union{Nothing, Dict{String, String}}=nothing,
               offset::Union{Nothing, Int}=nothing,
               limit::Union{Nothing, Int}=nothing,
               sample::Union{Nothing, Int}=nothing,
               sort::Union{Nothing, String}=nothing,
               order::Union{Nothing, String}=nothing,
               facet::Union{Nothing, Dict{String, Int}}=nothing,
               select::Union{Nothing, String}=nothing,
               works::Bool=false,
               cursor::Union{Nothing, String}=nothing)

         url="https://api.crossref.org/types"

         if !isnothing(type_id)
                 url = string(url, "/", type_id)
                 if works
                         url = string(url, "/works")
                 else
                         r = request(url=url)
                         return JSON.parse(String(r.body))
                 end

         end

         r = request(url=url,
                     query=query,
                     filter=filter,
                     offset=offset,
                     limit=limit,
                     sample=sample,
                     sort=sort,
                     order=order,
                     facet=facet,
                     select=select,
                     cursor=cursor)

        return JSON.parse(String(r.body))

end

"""
    prefixes(; kwargs...)

Arguments:
- `owner_prefix::Union{Nothing, String}`: A DOI owner prefix.
- `query::Union{Nothing, String}`: A query string.
- `filter::Union{Nothing, Dict{String, String}}`: Filter options.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.
- `select::Union{Nothing, String}`: Select specific fields.
- `works::Bool`: If true, works returned as well.
- `cursor::Union{Nothing, String}`: Cursor character string to do deep paging.

Returns:
- `Dict`
"""
function prefixes(; owner_prefix::Union{Nothing, String}=nothing,
                  query::Union{Nothing, String}=nothing,
                  filter::Union{Nothing, Dict{String, String}}=nothing,
                  offset::Union{Nothing, Int}=nothing,
                  limit::Union{Nothing, Int}=nothing,
                  sample::Union{Nothing, Int}=nothing,
                  sort::Union{Nothing, String}=nothing,
                  order::Union{Nothing, String}=nothing,
                  facet::Union{Nothing, Dict{String, Int}}=nothing,
                  select::Union{Nothing, String}=nothing,
                  works::Bool=false,
                  cursor::Union{Nothing, String}=nothing)

         url="https://api.crossref.org/prefixes"

         if !isnothing(owner_prefix)
                 url = string(url, "/", owner_prefix)
                 if works
                         url = string(url, "/works")
                 else
                         r = request(url=url)
                         return JSON.parse(String(r.body))
                 end

         end

         r = request(url=url,
                     query=query,
                     filter=filter,
                     offset=offset,
                     limit=limit,
                     sample=sample,
                     sort=sort,
                     order=order,
                     facet=facet,
                     select=select,
                     cursor=cursor)

        return JSON.parse(String(r.body))

end

"""
    licenses(; kwargs...)

Arguments:
- `query::Union{Nothing, String}`: A query string.
- `offset::Union{Nothing, Int}`: Number of record to start at.
- `limit::Union{Nothing, Int}`: Number of results to return.
- `sample::Union{Nothing, Int}`: Number of random results to return.
- `sort::Union{Nothing, String}`: Field to sort on.
- `order::Union{Nothing, String}`: Sort order, one of ‘asc’ or ‘desc’.
- `facet::Union{Nothing, Dict{String, Int}}`: Facet options.

Returns:
- `Dict`
"""
function licenses(; query::Union{Nothing, String}=nothing,
                  offset::Union{Nothing, Int}=nothing,
                  limit::Union{Nothing, Int}=nothing,
                  sample::Union{Nothing, Int}=nothing,
                  sort::Union{Nothing, String}=nothing,
                  order::Union{Nothing, String}=nothing,
                  facet::Union{Nothing, Dict{String, Int}}=nothing)

         url="https://api.crossref.org/licenses"

         r = request(url=url,
                     query=query,
                     offset=offset,
                     limit=limit,
                     sample=sample,
                     sort=sort,
                     order=order,
                     facet=facet)

        return JSON.parse(String(r.body))

end

"""
    agency(doi)

Arguments:
- `doi::String`: DOI of the work.

Returns:
- Dict
"""
function agency(doi::String)

    url = string("https://api.crossref.org/works/", doi, "/agency")

    r = request(url=url)

    return JSON.parse(String(r.body))

end

max_stars_repo_name: naustica/Pitaya

max_issues_repo_name: naustica/Pitaya

max_forks_repo_name: naustica/Pitaya


text:
@testset "Assign bin labels" begin
    D = 3
    E = cembed([rand(30) for i = 1:D])
    npts = size(E.points, 2)

    @testset "ϵ is an Int" begin
        labels = assign_bin_labels(E, 10)
        @test size(labels, 1) == D
        @test size(labels, 2) == npts
    end

    @testset "ϵ is a Float64" begin
        labels = assign_bin_labels(E, 0.2)
        @test size(labels, 1) == D
        @test size(labels, 2) == npts
    end

    @testset "ϵ is a Vector{Float64}" begin
        labels = assign_bin_labels(E, [0.2, 0.3, 0.1])
        @test size(labels, 1) == D
        @test size(labels, 2) == npts
    end

    @testset "ϵ is a Vector{Int}" begin
        labels = assign_bin_labels(E, [3, 4, 2])
        @test size(labels, 1) == D
        @test size(labels, 2) == npts
    end
end

@testset "Marginal visitation frequency" begin
    o = rand(100, 3)
    x, y = o[:, 1], o[:, 2]
    E = cembed([x, y], [2, 2, 1], [1, 0, 0])
    n_bins = [4, 3, 4]
    # Which bins get visited by every point of the orbit?
    visited_bins_inds = assign_bin_labels(E, n_bins)
    along_which_axes = [1, 2]
    m = marginal_visitation_freq(along_which_axes, visited_bins_inds)

    @test sum(m) ≈ 1
end

# @testset "Simplex triangulation" begin
#     n_pts = 30
# 	@testset "Triangulation" begin
#         E_3D = cembed([rand(n_pts) for i = 1:3])
#         E_4D = cembed([rand(n_pts) for i = 1:4])
#
#         T_3D = triangulate(E_3D)
#         @test typeof(T_3D) <: Partitioning.AbstractTriangulation
#
#         T_4D = triangulate(E_4D)
#         @test typeof(T_4D) <: Partitioning.AbstractTriangulation
#     end
#
#     @testset "LinearlyInvariantTriangulation" begin
#         E_3D = invariantize(cembed([rand(n_pts) for i = 1:3]))
#         E_4D = invariantize(cembed([rand(n_pts) for i = 1:4]))
#
#         T_3D = triangulate(E_3D)
#         @test typeof(T_3D) == Partitioning.LinearlyInvariantTriangulation
#
#         T_4D = triangulate(E_4D)
#         @test typeof(T_4D) == Partitioning.LinearlyInvariantTriangulation
#
#     end
# end


text:
using FastTransforms
using Base.Test

println("Testing special functions")
n = 0:1000_000
λ = 0.123
@time FastTransforms.Cnλ(n,λ);
@time FastTransforms.Cnλ(n,λ);

x = linspace(0,20,81);
@test norm((FastTransforms.Λ(x)-FastTransforms.Λ(big(x)))./FastTransforms.Λ(big(x)),Inf) < 2eps()

x = 0:0.5:10_000
λ₁,λ₂ = 0.125,0.875
@test norm((FastTransforms.Λ(x,λ₁,λ₂)-FastTransforms.Λ(big(x),big(λ₁),big(λ₂)))./FastTransforms.Λ(big(x),big(λ₁),big(λ₂)),Inf) < 4eps()
λ₁,λ₂ = 1//3,2//3
@test norm((FastTransforms.Λ(x,Float64(λ₁),Float64(λ₂))-FastTransforms.Λ(big(x),big(λ₁),big(λ₂)))./FastTransforms.Λ(big(x),big(λ₁),big(λ₂)),Inf) < 4eps()

n = 0:1000
α = 0.125
β = 0.375

@time FastTransforms.Cnαβ(n,α,β);
@test norm(FastTransforms.Cnαβ(n,α,β)./FastTransforms.Cnαβ(n,big(α),big(β))-1,Inf) < 3eps()

@time FastTransforms.Anαβ(n,α,β);
@test norm(FastTransforms.Anαβ(n,α,β)./FastTransforms.Anαβ(n,big(α),big(β))-1,Inf) < 4eps()

N = 20
f(x) = exp(x)

x,w = FastTransforms.fejer1(N,0.,0.)
@test norm(dot(f(x),w)-2sinh(1)) ≤ 4eps()
x,w = FastTransforms.fejer2(N,0.,0.)
@test norm(dot(f(x),w)-2sinh(1)) ≤ 4eps()
x,w = FastTransforms.clenshawcurtis(N,0.,0.)
@test norm(dot(f(x),w)-2sinh(1)) ≤ 4eps()

#=
x = Fun(identity)
g = f(x)*(1-x)^.25*(1+x)^.35
val = sum(g)
=#

x,w = FastTransforms.fejer1(N,0.25,0.35)
@test norm(dot(f(x),w)-2.0351088204147243) ≤ 4eps()
x,w = FastTransforms.fejer2(N,0.25,0.35)
@test norm(dot(f(x),w)-2.0351088204147243) ≤ 4eps()
x,w = FastTransforms.clenshawcurtis(N,0.25,0.35)
@test norm(dot(f(x),w)-2.0351088204147243) ≤ 4eps()

println("Testing the Chebyshev–Jacobi transform")

Nr = 5
v = zeros(Nr)
Na,Nb = 5,5
V = zeros(Na,Nb)

for N in round(Int,logspace(1,3,3))
    println("")
    println("N = ",N)
    println("")
    for αi=1:Na, βi=1:Nb
        α,β = -0.45+(αi-1)/Na,-0.45+(βi-1)/Nb
        p1,p2 = plan_cjt(rand(N),α,β),plan_icjt(rand(N),α,β)
        for i=1:Nr
            c = rand(N)
            v[i] = log(norm(p2*(p1*c)-c,Inf)/(20N^(1+2max(α,β))*eps()))
        end
        mean(v) > 2 && println("This is the mean log(||Error||_∞/Estimate) that broke the test: ",mean(v)," and the standard deviation: ",std(v))
        V[αi,βi] = mean(v)
    end
    println("This is the mean log(||Error||_∞/Estimate): ",mean(V)," and the standard deviation: ",std(V))
end

println("Testing the special cases length(c) = 0,1,2:")

c = [1.5]
@test cjt(c,0.12,0.34) == c
@test icjt(c,0.12,0.34) == c

pop!(c)

@test cjt(c,0.12,0.34) == c
@test icjt(c,0.12,0.34) == c

c = [1.0;2.0]
@test norm(jjt(c,0.12,0.34,0.12,0.34)-c,Inf) ≤ 2eps()

# test same parameters

c=[0.,1.]
@test jjt(c,1.5,-0.5,1.5,-0.5) == c

println("Testing the special cases (α,β) = (±0.5,±0.5)")

# Chebyshev coefficients of exp
c_cheb = [1.2660658777520084,1.13031820798497,0.27149533953407656,0.044336849848663804,0.005474240442093732,0.0005429263119139438,4.497732295429515e-5,3.1984364624019905e-6,1.9921248066727958e-7,1.1036771725517344e-8,5.505896079673747e-10,2.4979566169849825e-11,1.03915223067857e-12,3.9912633564144015e-14,1.4237580108256572e-15,4.740926102561496e-17,1.4801800572082976e-18,4.3499194949441696e-20,1.2074289272797528e-21,3.175356737059445e-23,7.93367197163804e-25,1.8879484042289162e-26,4.2886738765925875e-28,9.318985281777577e-30,1.9406469749017399e-31,3.8798022492260015e-33,7.45850288739152e-35,1.3807477824110646e-36,2.4648623717710964e-38,4.2485421925059915e-40,7.079001176212895e-42,1.141486778254089e-43,1.783151037543284e-45,2.7011422638738024e-47,3.9714338657415316e-49,5.672351695612112e-51,7.876788130533199e-53,1.0642416282087358e-54,1.4000817885528905e-56,1.794689085390544e-58,2.243019482662972e-60,2.7349926016641967e-62,3.2554929833033215e-64,3.78495689402033e-66,4.300544446926769e-68,4.77780574839433e-70,5.192665133046628e-72,5.523442810988436e-74,5.717985746408847e-76]

# Jacobi coefficients of exp at (±0.5,±0.5)
c_00 = [1.2660658777520084,2.26063641596994,0.723987572090871,0.14187791951572415,0.020020079331085648,0.0022061767595233264,0.00019937999439479753,1.5268992715616027e-5,1.014420290055232e-6,5.9506859054175615e-8,3.124851418974203e-9,1.4852160704995775e-10,6.447143371527455e-12,2.5753240175931934e-13,9.526907336854026e-15,3.281725133829881e-16,1.0576496993856426e-17,3.20238467156949e-19,9.14299140846687e-21,2.4694550834939014e-22,6.3281713790321685e-24,1.542622093731937e-25,3.5857224408995815e-27,7.964666406806272e-29,1.6939043515668796e-30,3.4556190450066486e-32,6.773312266576111e-34,1.2775611026088504e-35,2.3221236341654026e-37,4.072731101981327e-39,6.901079409575277e-41,1.1310395206213411e-42,1.794875997727078e-44,2.760732681473006e-46,4.119631242417262e-48,5.969295899580336e-50,8.405881613439669e-52,1.1512859570474932e-53,1.5347890791678689e-55,1.9929132418968355e-57,2.522290550037395e-59,3.1134869971034086e-61,3.7506700259829663e-63,4.411970068315684e-65,5.0705898146983576e-67,5.696604001638848e-69,6.259278760298182e-71,6.729592140384465e-73,7.039950057933299e-75]
c_01 = [1.8312249817444934,1.4018135475190467,0.4211095858436539,0.07969774446521205,0.011002819207328322,0.0011944708771926122,0.00010677917238679192,8.110005868771403e-6,5.353106062476433e-7,3.1237733951100556e-8,1.6333110219427631e-9,7.735005972383578e-11,3.3473853404557083e-12,1.333595312027857e-13,4.92207038322879e-15,1.6920924742289326e-16,5.443658341283791e-18,1.6456374329092366e-19,4.691719175403533e-21,1.2655773772197326e-22,3.2393803393291905e-24,7.888321906112733e-26,1.8318188278743866e-27,4.065263937281931e-29,8.638847091039724e-31,1.7610248038105726e-32,3.449351261471638e-34,6.50183837007916e-36,1.1810743750838164e-37,2.070295858087742e-39,3.506179552173043e-41,5.743539156119835e-43,9.110325158556408e-45,1.4006615828866467e-46,2.0892357224279913e-48,3.026093616078655e-50,4.2597272086552854e-52,5.832159509538512e-54,7.772313549445791e-56,1.0089104305392274e-57,1.2765228632360985e-59,1.5752735945134057e-61,1.8971383534453687e-63,2.2310498815373167e-65,2.563461449357282e-67,2.8792582164274203e-69,3.1629293838222696e-71,3.399629156416382e-73,3.5199750289666496e-75]
c_10 = [0.7009067737595234,0.8588228684508934,0.30287798624721707,0.062180175050512106,0.009017260123757326,0.0010117058823307142,9.260082200800561e-5,7.1589868468446235e-6,4.791096838075886e-7,2.8269125103075062e-8,1.49154039703144e-9,7.117154732612197e-11,3.099758031071747e-12,1.2417287055653363e-13,4.604836953625235e-15,1.5896326596009483e-16,5.1328386525726346e-18,1.5567472386602532e-19,4.451272233063337e-21,1.2038777062741688e-22,3.088791039702978e-24,7.537899031206638e-26,1.753903613025195e-27,3.899402469524341e-29,8.300196424629072e-31,1.694594241196076e-32,3.3239610051044727e-34,6.273772656009344e-36,1.1410492590815863e-37,2.002435243893585e-39,3.3948998574022335e-41,5.566856050093576e-43,8.838434818714371e-45,1.3600710985863591e-46,2.0303955199892708e-48,2.943202283501681e-50,4.146154404784383e-52,5.68070006093642e-54,7.575577242232898e-56,9.84002811357608e-58,1.2457676868012966e-59,1.5382134025900027e-61,1.8535316725375976e-63,2.1809201867783675e-65,2.5071283653410757e-67,2.8173457852114277e-69,3.0963493764759127e-71,3.3299629839680836e-73,3.5199750289666496e-75]
c_11 = [1.13031820798497,0.7239875720908708,0.21281687927358628,0.040040158662171296,0.005515441898808317,0.0005981399831843923,5.34414745046561e-5,4.057681160220927e-6,2.677808657437903e-7,1.5624257094871016e-8,8.168688387747676e-10,3.8682860229164727e-11,1.6739606114355753e-12,6.668835135797819e-14,2.4612938503724105e-15,8.46119759508514e-17,2.722026970834067e-18,8.228692267620182e-20,2.3459823293192055e-21,6.328171379032168e-23,1.619753198418534e-24,3.944294684989539e-26,9.159366367827211e-28,2.0326852218802552e-29,4.31952380625831e-31,8.805305946548944e-33,1.7247074885219479e-34,3.2509730878315635e-36,5.905460097872925e-38,1.0351619114362916e-39,1.7531112569630786e-41,2.8718015963633246e-43,4.5552089244304594e-45,7.003373112109346e-47,1.044626782426559e-48,1.5130586904191403e-50,2.129879020537863e-52,2.916099250418951e-54,3.8861808216988285e-56,5.04458110007479e-58,6.382648344061989e-60,7.876407054564243e-62,9.485735646882059e-64,1.1155297592778338e-65,1.2817359024131178e-67,1.4396345613848593e-69,1.5814715428162693e-71,1.700107488097128e-73,1.7781317156635653e-75]

@test norm(cjt(c_00,-0.5,-0.5)-c_cheb,Inf) < eps()
@test norm(cjt(c_01,-0.5,0.5)-c_cheb,Inf) < eps()
@test norm(cjt(c_10,0.5,-0.5)-c_cheb,Inf) < eps()
@test norm(cjt(c_11,0.5,0.5)-c_cheb,Inf) < eps()

@test norm(icjt(c_cheb,-0.5,-0.5)-c_00,Inf) < eps()
@test norm(icjt(c_cheb,-0.5,0.5)-c_01,Inf) < eps()
@test norm(icjt(c_cheb,0.5,-0.5)-c_10,Inf) < eps()
@test norm(icjt(c_cheb,0.5,0.5)-c_11,Inf) < eps()

c = exp(-collect(1:1000)./30)

println("Testing increment/decrement operators for α,β ≤ -0.5")

α,β = -0.9,-0.6

@test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 10eps()
@test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 2000eps()
p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
@test norm(p2*(p1*c)-c,Inf) < 2000eps()

println("Testing increment/decrement operators for α ≤ -0.5, β > -0.5")

α,β = -0.9,3.1

@test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 10eps()
@test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 2000eps()
p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
@test norm(p2*(p1*c)-c,Inf) < 2000eps()

println("Testing increment/decrement operators for α > -0.5, β ≤ -0.5")

α,β = 2.4,-0.75

@test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 1200eps()
@test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 1e6eps()
p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
@test norm(p2*(p1*c)-c,Inf) < 1e6eps()

println("Testing increment/decrement operators for α > -0.5, β > -0.5")

α,β = 2.4,3.6

@test norm(FastTransforms.fromsquare!(FastTransforms.tosquare!(copy(c),α,β),α,β)-c,Inf) < 250eps()
@test norm(icjt(cjt(c,α,β),α,β)-c,Inf) < 2e5eps()
p1,p2 = plan_cjt(c,α,β),plan_icjt(c,α,β)
@test norm(p2*(p1*c)-c,Inf) < 2e5eps()

println("Testing for complex coefficients")

α,β = 0.12,0.34
c = complex(rand(100),rand(100))

@test cjt(c,α,β) == complex(cjt(real(c),α,β),cjt(imag(c),α,β))
@test icjt(c,α,β) == complex(icjt(real(c),α,β),icjt(imag(c),α,β))
@test jjt(c,α,β,α,β) == complex(jjt(real(c),α,β,α,β),jjt(imag(c),α,β,α,β))
@test norm(jjt(c,α,β,α,β)-c,Inf) < 200eps()

println("Testing for Vector{Float32}")

c64 = rand(100)
c32 = map(Float32,c64)

cL64 = cjt(c64,0.,0.)
cL32 = cjt(c32,0.f0,0.f0)

@test norm(cL32-cL64,Inf) < 20eps(Float32)

println("Testing for Matrix of coefficients")

c = rand(100,100)
@test maxabs(jjt(c,α,β,α,β)-c) < 10000eps()

println("Testing Gaunt coefficients")

include("gaunttest.jl")

println("Testing BigFloat support of FFT and DCT methods")

include("fftBigFloattest.jl")

println("Testing equivalence of CXN and ASY methods")

for k in round(Int,logspace(1,4,20))
    r = randn(k)./√(1:k) # Proven 𝒪(√(log N)) error for ASY method.
    @test_approx_eq leg2cheb(r) cjt(r,0.,0.)
end

@test_approx_eq leg2chebu([1.0,2,3,4,5])  [0.546875,0.5,0.5390625,1.25,1.3671875]

c = randn(1000)./√(1:1000);

@test_approx_eq leg2cheb(cheb2leg(c)) c
@test_approx_eq cheb2leg(leg2cheb(c)) c

@test norm(jac2jac(c,0.,√2/2,-1/4,√2/2)-jjt(c,0.,√2/2,-1/4,√2/2),Inf) < 10length(c)*eps()

@test norm(ultra2ultra(ultra2ultra(c,.5,.75),.75,.5)-c,Inf) < 10length(c)*eps()


text:
LOG_GROUP_COST_FUNCTIONS = :CostFunctionsConstructor
LOG_GROUP_OPTIMZATION_CONTAINER = :OptimizationContainer
LOG_GROUP_TYPE_REGISTRATIONS = :TypeRegistrations
LOG_GROUP_FEEDFORWARDS_CONSTRUCTION = :FeedforwardConstructor
LOG_GROUP_BRANCH_CONSTRUCTIONS = :BranchConstructor
LOG_GROUP_SERVICE_CONSTUCTORS = :ServicesConstructor
LOG_GROUP_MODELS_VALIDATION = :ModelValidation
LOG_GROUP_OPTIMIZATION_CONTAINER = :OptimizationContainer
LOG_GROUP_BUILD_INITIAL_CONDITIONS = :InitialConditionsBuild
LOG_GROUP_NETWORK_CONSTRUCTION = :NetworkConstructor
LOG_GROUP_MODEL_STORE = :ModelStore
LOG_GROUP_RESULTS = :Results
LOG_GROUP_SIMULATION_STORE = :SimulationStore


text:
module CorpusLoaders
using DataDeps
using WordTokenizers
using MultiResolutionIterators
using InternedStrings
using Glob
using StringEncodings
using CSV

export Document, TaggedWord, SenseAnnotatedWord, PosTaggedWord
export title, sensekey, word, named_entity, part_of_speech
export load

export WikiCorpus, SemCor, Senseval3, CoNLL, IMDB, Twitter, StanfordSentimentTreebank, WikiGold, CoNLL2000

function __init__()
    include(joinpath(@__DIR__, "WikiCorpus_DataDeps.jl"))
    include(joinpath(@__DIR__, "SemCor_DataDeps.jl"))
    include(joinpath(@__DIR__, "SemEval2007Task7_DataDeps.jl"))
    include(joinpath(@__DIR__, "Senseval3_DataDeps.jl"))
    include(joinpath(@__DIR__, "CoNLL_DataDeps.jl"))
    include(joinpath(@__DIR__, "IMDB_DataDeps.jl"))
    include(joinpath(@__DIR__, "Twitter_DataDeps.jl"))
    include(joinpath(@__DIR__, "StanfordSentimentTreebank_DataDeps.jl"))
    include(joinpath(@__DIR__, "WikiGold_DataDeps.jl"))
    include(joinpath(@__DIR__, "CoNLL2000_DataDeps.jl"))
end

include("types.jl")
include("apply_subparsers.jl")
include("WikiCorpus.jl")
include("SemCor.jl")
include("SemEval2007Task7.jl")
include("Senseval3.jl")
include("CoNLL.jl")
include("IMDB.jl")
include("Twitter.jl")
include("StanfordSentimentTreebank.jl")
include("WikiGold.jl")
include("CoNLL2000.jl")

end


text:
#
# This is the queue_mmc SimJulia example found on
# https://github.com/BenLauwens/SimJulia.jl/blob/master/examples/queue_mmc.ipynb
#
# only modified for benchmarking
#
#
# This is the queue_mmc SimJulia example found on
# https://github.com/BenLauwens/SimJulia.jl/blob/master/examples/queue_mmc.ipynb
#
# only modified for benchmarking
#
cd(@__DIR__)
using Distributions, SimJulia, ResumableFunctions, Random
using BenchmarkTools, Plots

const _bench = [false]

Random.seed!(8710) # set random number seed for reproducibility
num_customers = 10 # total number of customers generated
num_servers = 2 # number of servers
mu = 1.0 / 2 # service rate
lam = 0.9 # arrival rate
arrival_dist = Exponential(1 / lam) # interarrival time distriubtion
service_dist = Exponential(1 / mu); # service time distribution

@resumable function customer(env::Environment, server::Resource, id::Integer, time_arr::Float64, dist_serve::Distribution)
    @yield timeout(env, time_arr) # customer arrives
    _bench[end] || println("Customer $id arrived: ", now(env))
    @yield request(server) # customer starts service
    _bench[end] || println("Customer $id entered service: ", now(env))
    @yield timeout(env, rand(dist_serve)) # server is busy
    @yield release(server) # customer exits service
    _bench[end] || println("Customer $id exited service: ", now(env))
end

function run_model(arrival_dist, service_dist, num_customers, num_servers)
    sim = Simulation() # initialize simulation environment
    server = Resource(sim, num_servers) # initialize servers
    arrival_time = 0.0
    for i = 1:num_customers # initialize customers
        arrival_time += rand(arrival_dist)
        @process customer(sim, server, i, arrival_time, service_dist)
    end
    run(sim) # run simulation
end

function run_model(num_customers)
    result = @benchmark run_model($arrival_dist, $service_dist, $num_customers, $num_servers)
    return mean(result).time*1e-9
end

_bench[end] = true
N = [10,1000,2000,4000]

times = run_model.(N)
plot(N, times, xlabel="Customers", ylabel="Time (seconds)", leg=false, grid=false)
savefig("SimJulia_Queue.png")


text:
@testset "parser API" begin
    @testset "String and buffer input" begin
        # String
        @test parse(Expr, "x+y\nz") == (Expr(:toplevel, :(x+y), :z), [], 6)
        @test parse(Expr, "x+y\nz", rule=:statement) == (:(x+y), [], 4)
        @test parse(Expr, "x+y\nz", rule=:atom) == (:x, [], 2)
        @test parse(Expr, "x+y\nz", 5, rule=:atom) == (:z, [], 6)

        # Vector{UInt8}
        @test parse(Expr, Vector{UInt8}("x+y"), rule=:statement) == (:(x+y), [], 4)
        @test parse(Expr, Vector{UInt8}("x+y"), 3, rule=:statement) == (:y, [], 4)
        # Ptr{UInt8}, len
        code = "x+y"
        GC.@preserve code begin
            stream = ParseStream(pointer(code), 3)
            parse(stream, rule=:statement)
            @test JuliaSyntax.build_tree(Expr, stream) == :(x+y)
            @test JuliaSyntax.last_byte(stream) == 3
        end

        # SubString
        @test parse(Expr, SubString("x+y"), rule=:statement) == (:(x+y), [], 4)
        @test parse(Expr, SubString("x+y"), 1, rule=:atom) == (:x, [], 2)
        @test parse(Expr, SubString("x+y"), 3, rule=:atom) == (:y, [], 4)
        @test parse(Expr, SubString("x+y",3,3), 1, rule=:atom) == (:y, [], 2)
        @test parse(Expr, SubString("α+x"), rule=:statement) == (:(α+x), [], 5)
    end

    @testset "IO input" begin
        # IOBuffer
        io = IOBuffer("x+y")
        @test parse(Expr, io, rule=:statement) == (:(x+y), [])
        @test position(io) == 3
        io = IOBuffer("x+y")
        seek(io, 2)
        @test parse(Expr, io, rule=:atom) == (:y, [])
        @test position(io) == 3
        # A GenericIOBuffer, not actually IOBuffer
        io = IOBuffer(SubString("x+y"))
        @test parse(Expr, io, rule=:statement) == (:(x+y), [])
        @test position(io) == 3
        # Another type of GenericIOBuffer
        io = IOBuffer(codeunits("x+y"))
        @test parse(Expr, io, rule=:statement) == (:(x+y), [])
        @test position(io) == 3
        # IOStream
        mktemp() do path, io
            write(io, "x+y")
            close(io)

            open(path, "r") do io
                @test parse(Expr, io, rule=:statement) == (:(x+y), [])
                @test position(io) == 3
            end
        end
    end

    @testset "parseall" begin
        @test parseall(Expr, " x ") == Expr(:toplevel, :x)
        @test parseall(Expr, " x ", rule=:statement) == :x
        @test parseall(Expr, " x ", rule=:atom) == :x
        # TODO: Fix this situation with trivia here; the brackets are trivia, but
        # must be parsed to discover the atom inside. But in GreenTree we only
        # place trivia as siblings of the leaf node with identifier `x`, not as
        # children.
        @test_broken parseall(Expr, "(x)", rule=:atom) == :x

        @test_throws JuliaSyntax.ParseError parseall(Expr, "x+y", rule=:atom)
        @test_throws JuliaSyntax.ParseError parseall(Expr, "x+y\nz", rule=:statement)
    end
end


text:
function sieve(x::Int)
    sieve = fill(true, x)
    for i in 2:isqrt(x)
        if sieve[i]
            j = 2i
            while j <= x
                sieve[j] = false
                j += i
            end
        end
    end
    primes = findall(sieve)
    popfirst!(primes)
    return primes
end

max_stars_repo_name: ellygaytor/Projects

max_issues_repo_name: ellygaytor/Projects

max_forks_repo_name: ellygaytor/Projects


text:
# @inline function loopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{0}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{1}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] += Cₘₙ
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, ::StaticInt{1}, β, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{0}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{1}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] += Cₘₙ
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, ::StaticInt{-1}, β, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
#     nothing
# end
# @inline function loopmul!(C, A, B, α, ::StaticInt{0}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = α * Cₘₙ
#     end
#     nothing
# end

# @inline function loopmul!(C, A, B, α, ::StaticInt{1}, M, K, N)
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] += α * Cₘₙ
#     end
#     nothing
# end
@inline function loopmul!(C, A, B, α, β, M, K, N)
    @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
        Cₘₙ = zero(eltype(C))
        for k ∈ CloseOpen(K)
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n] = α * Cₘₙ + β * C[m,n]
    end
    nothing
end

# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{1}, ::StaticInt{0},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] = Cₘₙ
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{1}, ::StaticInt{1},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] += Cₘₙ
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{1}, β,
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{-1}, ::StaticInt{0},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ -= Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] = Cₘₙ
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{-1}, ::StaticInt{1},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ -= Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] += Cₘₙ
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     ::StaticInt{-1}, β,
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ -= Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     α, ::StaticInt{0},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] = α * Cₘₙ
#     end
# end
# @inline function packamul!(
#     C, Ãₚ, A, B,
#     α, ::StaticInt{1},
#     M, K, N
# )
#     @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Aₘₖ = A[m,k]
#             Cₘₙ += Aₘₖ * B[k,n]
#             Ãₚ[m,k] = Aₘₖ
#         end
#         C[m,n] += α * Cₘₙ
#     end
# end
@inline function packamul!(
    C, Ãₚ, A, B,
    α, β, M, K, N
)
    @avx for n ∈ CloseOpen(N), m ∈ CloseOpen(M)
        Cₘₙ = zero(eltype(C))
        for k ∈ CloseOpen(K)
            Aₘₖ = A[m,k]
            Cₘₙ += Aₘₖ * B[k,n]
            Ãₚ[m,k] = Aₘₖ
        end
        C[m,n] = α * Cₘₙ + β * C[m,n]
   end
end
@inline function alloc_a_pack(A, M, ::Type{T}) where {T}
    buffer = first_cache_buffer(T)
    bufferptr = if Sys.WORD_SIZE == 32
        buffer
    else
        align(pointer(buffer))
    end
    Apack = default_zerobased_stridedpointer(bufferptr, (One(), align(M, T)))
    Apack, buffer
end
@inline function packaloopmul!(
    C::AbstractStridedPointer{T},
    A::AbstractStridedPointer,
    B::AbstractStridedPointer,
    α, β, M, K, N
) where {T}
    Ãₚ, buffer = alloc_a_pack(A, M, T)
    GC.@preserve buffer begin
        Mᵣ, Nᵣ = matmul_params()
        packamul!(C, Ãₚ, A, B, α, β, M, K, Nᵣ)
        loopmul!(gesp(C, (Zero(), Nᵣ)), Ãₚ, gesp(B, (Zero(), Nᵣ)), α, β, M, K, N - Nᵣ)
    end
    nothing
end

# @inline function inlineloopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{0}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, ::StaticInt{1}, ::StaticInt{1}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] += Cₘₙ
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, ::StaticInt{1}, β, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{0}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, ::StaticInt{1}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] += Cₘₙ
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, ::StaticInt{-1}, β, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ -= A[m,k] * B[k,n]
#         end
#         C[m,n] = Cₘₙ + β * C[m,n]
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, α, ::StaticInt{0}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] = α * Cₘₙ
#     end
#     C
# end
# @inline function inlineloopmul!(C, A, B, α, ::StaticInt{1}, M, K, N)
#     @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
#         Cₘₙ = zero(eltype(C))
#         for k ∈ CloseOpen(K)
#             Cₘₙ += A[m,k] * B[k,n]
#         end
#         C[m,n] += α * Cₘₙ
#     end
#     C
# end
@inline function inlineloopmul!(C, A, B, α, β, M, K, N)
    @avx inline=true for m ∈ CloseOpen(M), n ∈ CloseOpen(N)
        Cₘₙ = zero(eltype(C))
        for k ∈ CloseOpen(K)
            Cₘₙ += A[m,k] * B[k,n]
        end
        C[m,n]  = α * Cₘₙ + β * C[m,n]
    end
    C
end


text:
# 04_experiment_demo.jl

using Plots     # You should import Plots to use Experiments
using MAB

ex1 = Sutton2017.fig_2_2( agents = [ (GradientBandit,(0.1)) ] )
f = run( ex1, T = 100, rounds = 100 )
plot( f[1], reuse = false, show=true )
plot( f[2], reuse = false,  show=true  )


text:
export TR_SecA_ls
"""
A one dimensionnal trust region algorithm using the improved secant as the
approximation for the second derivative in the second order Taylor expansion.
For more documentation referer to TR_generic_ls.
"""
function TR_SecA_ls(h :: LineModel, stop_ls :: LS_Stopping; kwargs...)
    ls_meta = LS_Function_Meta(dir = :SecA)
    (state, stop_ls.meta.optimal) = TR_generic_ls(h, stop_ls; f_meta = ls_meta, kwargs...)
    return (state, stop_ls.meta.optimal)
end


text:
_init_logJ(data) = zeros(eltype(data), 1, size(data, 2))

sqnorm(x) = sum(abs2, x)
l2_reg(ps) = sum(sqnorm, ps)


text:
const UNIXEPOCH_TS = Dates.value(DateTime(1970))  # millisceonds
const UNIXEPOCH_DT = Dates.value(Date(1970))

value(x) = x

abstract type ArrowTime <: Dates.AbstractTime end

ArrowTime(t::Dates.TimeType) = convert(ArrowTime, t)

"""
    Timestamp{P<:Dates.TimePeriod} <: ArrowTime

Timestamp in which time is stored in units `P` as `Int64` for Arrow formatted data.
"""
struct Timestamp{P<:TimePeriod} <: ArrowTime
    value::Int64
end

Timestamp(t::P) where P<:TimePeriod = Timestamp{P}(Dates.value(t))
Timestamp{P}(t::DateTime) where P<:TimePeriod = convert(Timestamp{P}, t)
Timestamp(t::DateTime) = convert(Timestamp, t)

value(t::Timestamp) = t.value
unitvalue(t::Timestamp{P}) where P = P(value(t))

scale(::Type{D}, t::Timestamp{P}) where {D,P} = convert(D, unitvalue(t))

function Base.convert(::Type{DateTime}, t::Timestamp{P}) where P
    DateTime(Dates.UTM(UNIXEPOCH_TS + Dates.value(scale(Millisecond, t))))
end
Base.convert(::Type{TimeType}, t::Timestamp) = convert(DateTime, t)

function Base.convert(::Type{Timestamp{P}}, t::DateTime) where P
    Timestamp(convert(P, Millisecond(Dates.value(t) - UNIXEPOCH_TS)))
end
Base.convert(::Type{Timestamp}, t::DateTime) = convert(Timestamp{Millisecond}, t)
Base.convert(::Type{ArrowTime}, t::DateTime) = convert(Timestamp, t)

Base.show(io::IO, t::Timestamp) = show(io, convert(DateTime, t))

"""
    TimeOfDay{P<:Dates.TimePeriod,T<:Union{Int32,Int64}} <: ArrowTime

An arrow formatted object for representing the time of day.
Underlying data is `Int32` for seconds and milliseconds, `Int64` for microsecond and nanosecond.
"""
struct TimeOfDay{P<:TimePeriod,T<:Union{Int32,Int64}} <: ArrowTime
    value::T
end

function TimeOfDay{P,T}(t::P) where {P<:TimePeriod,T<:Union{Int32,Int64}}
    TimeOfDay{P,T}(Dates.value(convert(P, t)))
end
TimeOfDay{P}(t::P) where P<:Union{Second,Millisecond} = TimeOfDay{P,Int32}(t)
TimeOfDay{P}(t::P) where P<:Union{Microsecond,Nanosecond} = TimeOfDay{P,Int64}(t)
TimeOfDay(t::P) where P<:TimePeriod = TimeOfDay{P}(t)
TimeOfDay(t::Time) = convert(TimeOfDay, t)

value(t::TimeOfDay) = t.value
unitvalue(t::TimeOfDay{P}) where P = P(value(t))

scale(::Type{D}, t::TimeOfDay{P}) where {D,P} = convert(D, unitvalue(t))

function Base.convert(::Type{Time}, t::TimeOfDay{P}) where P
    Time(Nanosecond(scale(Nanosecond, t)))
end
Base.convert(::Type{TimeType}, t::TimeOfDay) = convert(Time, t)

Base.convert(::Type{TimeOfDay{P,T}}, t::Time) where {P,T} = TimeOfDay{P,T}(convert(P, t.instant))
Base.convert(::Type{TimeOfDay{P}}, t::Time) where P = TimeOfDay{P}(convert(P, t.instant))
Base.convert(::Type{TimeOfDay}, t::Time) = convert(TimeOfDay{Nanosecond}, t)
Base.convert(::Type{ArrowTime}, t::Time) = convert(TimeOfDay, t)

Base.show(io::IO, t::TimeOfDay) = show(io, convert(Time, t))

"""
    Datestamp <: ArrowTime

Stores a date as an `Int32` for Arrow formatted data.
"""
struct Datestamp <: ArrowTime
    value::Int32
end

Datestamp(t::Date) = convert(Datestamp, t)

value(t::Datestamp) = t.value

Base.convert(::Type{Date}, t::Datestamp) = Date(Dates.UTD(UNIXEPOCH_DT + value(t)))
Base.convert(::TimeType, t::Datestamp) = convert(Date, t)

Base.convert(::Type{Datestamp}, t::Date) = Datestamp(Dates.value(t) - UNIXEPOCH_DT)
Base.convert(::Type{ArrowTime}, t::Date) = convert(Datestamp, t)

Base.show(io::IO, t::Datestamp) = show(io, convert(Date, t))

#======================================================================================================
    some basic utilities for dates and datetimes...
======================================================================================================#
for symb ∈ [:(Base.:(>)), :(Base.:(≥)), :(Base.:(<)), :(Base.:(≤)), :(Base.:(==))]
    @eval $symb(t1::T, t2::T) where {T<:ArrowTime} = $symb(value(t1), value(t2))
end


text:
function mergelist(a, b)
    out = Vector{Int}()
    while !isempty(a) && !isempty(b)
        if a[1] < b[1]
            push!(out, popfirst!(a))
        else
            push!(out, popfirst!(b))
        end
    end
    append!(out, a)
    append!(out, b)
    out
end

function strand(a)
    i, s = 1, [popfirst!(a)]
    while i < length(a) + 1
        if a[i] > s[end]
            append!(s, splice!(a, i))
        else
            i += 1
        end
    end
    s
end

strandsort(a) = (out = strand(a); while !isempty(a) out = mergelist(out, strand(a)) end; out)

println(strandsort([1, 6, 3, 2, 1, 7, 5, 3]))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## Build menu tree structure for explorer

function tree_view_values(model::Model)
    all_subcomps = []
    for comp_def in compdefs(model)
        subcomp = tree_view_values(model, nameof(comp_def), comp_def)
        push!(all_subcomps, subcomp)
    end

    # Return sorted list so that the UI list of items will be in lexicographic order
    return sort(all_subcomps, by = x -> lowercase(x["name"]))
end

function tree_view_values(model::Model, comp_name::Symbol, comp_def::AbstractComponentDef)
    sub_comp_item = _tree_view_node(comp_name)
    for subcomp in compdefs(comp_def)
        push!(sub_comp_item["children"], tree_view_values(model, nameof(subcomp), subcomp));
    end
    return sub_comp_item
end

function _tree_view_node(comp_name::Symbol)
    return Dict("name" => "$comp_name", "children" => Dict[])
end

# Create the list of variables and parameters
function menu_item_list(model::Model)
    var_menuitems = []
    par_menuitems = []

    for comp_def in compdefs(model)
        all_subcomp_values = menu_item_list(model, nameof(comp_def), comp_def)
        append!(var_menuitems, all_subcomp_values["vars"])
        append!(par_menuitems, all_subcomp_values["pars"])
    end

    # Return sorted list so that the UI list of items will be in lexicographic order
    return Dict("vars" => sort(var_menuitems, by = x -> lowercase(x["name"])),"pars" => sort(par_menuitems, by = x -> lowercase(x["name"])))
end

# Create the list of variables and parameters
function menu_item_list(m::Model, comp_name::Symbol, comp_def::AbstractComponentDef)
    var_menu_items = map(var_name -> _menu_item(m, Symbol(comp_name), var_name), variable_names(comp_def));
    par_menu_items = map(par_name -> _menu_item(m, Symbol(comp_name), par_name), parameter_names(comp_def));

    # Return sorted list so that the UI list of items will be in lexicographic order
    return Dict("vars" => sort(var_menu_items, by = x -> lowercase(x["name"])),"pars" => sort(par_menu_items, by = x -> lowercase(x["name"])))
end

function menu_item_list(sim_inst::SimulationInstance)
    all_menuitems = []
    for datum_key in sim_inst.sim_def.savelist
        menu_item = _menu_item(sim_inst, datum_key)
        if menu_item !== nothing
            push!(all_menuitems, menu_item)
        end
    end

    # Return sorted list so that the UI list of items will be in lexicographic order
    return sort(all_menuitems, by = x -> lowercase(x["name"]))
end

function _menu_item(m::Model, comp_name::Symbol, item_name::Symbol)
    dims = dim_names(m, comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    if length(dims) == 0
        paths = _get_all_paths(m)
        comp_path = paths[comp_name];
        value = m[comp_path, item_name]
        name = "$comp_name : $item_name = $value"
    elseif length(dims) > 2
        @warn("$comp_name.$item_name has > 2 indexed dimensions, not yet implemented in explorer")
        name = "$comp_name : $item_name (CANNOT DISPLAY)"
    else
        name = "$comp_name : $item_name"          # the name is needed for the list label
    end

    menu_item = Dict("name" => name, "comp_name" => comp_name, "item_name" => item_name)
    return menu_item
end

function _menu_item(sim_inst::SimulationInstance, datum_key::Tuple{Symbol, Symbol})
    (comp_name, item_name) = datum_key
    dims = dim_names(sim_inst.models[1], comp_name, item_name)
    if length(dims) > 2
        # Drop references to singleton dimensions
        dims = tuple([dim for dim in dims if dim_count(m, dim) != 1]...)
    end

    if length(dims) > 2
        @warn("$comp_name.$item_name has >2 graphing dims, not yet implemented in explorer")
        name = "$comp_name : $item_name (CANNOT DISPLAY)"
    else
        name = "$comp_name : $item_name"          # the name is needed for the list label
    end

    menu_item = Dict("name" => "$item_name", "comp_name" => comp_name, "item_name" => item_name)
    return menu_item
end


text:
using Diban

path = "tmp.parquet"
#v1 = read_column(path, 1)
col_num=1

par = ParFile(path)
filemetadata = Diban.metadata(path)

T = TYPES[filemetadata.schema[col_num+1]._type+1]
# TODO detect if missing is necessary
res = Vector{Union{Missing, T}}(missing, nrows(par))
write_cursor = 1
#for row_group in filemetadata.row_groups
row_group =  filemetadata.row_groups[1]
pgs = pages(par, row_group.columns[col_num])
# the first page is always the dictionary page
dictionary_page = pgs[1]

# TODO different logic for different stuff
dictionary_of_values = T.(values(par, dictionary_page)[1])

data_page = pgs[2]

@which values(par, data_page)

#@time vals, repetition, decode = values(par, data_page)

################################################################################
# begin: expansion of functions `values`
################################################################################

page = data_page
ctype = coltype(page.colchunk)
rawbytes = bytes(page)
io = IOBuffer(rawbytes)
encs = page_encodings(page)
num_values = page_num_values(page)
typ = page.hdr._type

(typ == PageType.DATA_PAGE) || (typ == PageType.DATA_PAGE_V2)

@which read_levels_and_values(io, encs, ctype, num_values, par, page)
# below line is my code
vals, repetition, dd = read_levels_and_values(io, encs, ctype, num_values, par, page)

	################################################################################
	# begin: expansion of functions `read_levels_and_values`
	################################################################################
cname = colname(page.colchunk)
enc, defn_enc, rep_enc = encs

position(io)
#defn_levels = read(io, 8)

#@debug("before reading repn levels bytesavailable in page: $(bytesavailable(io))")
# read repetition levels. skipped if all columns are at 1st level
# TODO ZJ why is it checking for a "." in the name
repn_levels = ('.' in cname) ? read_levels(io, max_repetition_level(par.schema, cname), rep_enc, num_values) : Int[]

#@debug("before reading defn levels bytesavailable in page: $(bytesavailable(io))")
# read definition levels. skipped if column is required
defn_levels = isrequired(par.schema, cname) ? Int[] : read_levels(io, max_definition_level(par.schema, cname), defn_enc, num_values)

num_values
bw = 1

@which read_hybrid(io, num_values, bw)

read_hybrid(io, num_values, bw)

runhdr = 2
@which read_bitpacked_run(io, runhdr, bits, byt, typ, subarr)

# TODO I don't know how to deal with anything else; so below
@assert all(<=(1), defn_levels)

#@debug("before reading values bytesavailable in page: $(bytesavailable(io))")
# read values
# if there are missing values in the data then
# where defn_levels's elements == 1 are present and only
# sum(defn_levels) values can be read.
# because defn_levels == 0 are where the missing vlaues are
nmissing = sum(==(0), defn_levels)
#vals = read_values(io, enc, ctype, num_values - nmissing)

		################################################################################
		# begin: expansion of functions `read_values`
		################################################################################
@which read_values(io, enc, ctype, num_values - nmissing)

enc == Encoding.PLAIN
enc == Encoding.PLAIN_DICTIONARY || enc == Encoding.RLE_DICTIONARY

#vals = read_rle_dict(io, num_values)
@which read_rle_dict(io, num_values)

# begin: read_rle_dict

bits = read(io, UInt8)
Int(bits)
@debug("reading rle dictionary bits:$bits")
#arr = read_hybrid(io, num_values, bits; read_len=false)
@debug("read $(length(arr)) dictionary values")
#arr
@which read_hybrid(io, num_values, bits; read_len=false)
# end: read_rle_dict

# begin: read_hybrid
read_len = false
position(io)

#abc =  read(io, 2)
#Int.(abc)

using LittleEndianBase128

#LittleEndianBase128.decode(abc)[1] |> Int
#LittleEndianBase128.decodesigned(abc)[1] |> Int

isbitpack

len = read_len ? read_fixed(io, Int32) : Int32(0)
@debug("reading hybrid data", len, num_values, bits)
arrpos = 1
#while arrpos <= num_values

byt = bit2bytewidth(bits)
typ = byt2itype(byt)
@which _read_varint(io, Int)
Int(MSB)

read(io, UInt8)
read(io, UInt8)

using LittleEndianBase128
LittleEndianBase128.decode([0x02, 0x00])

Int(0x7f)
0x02

res = Int(0x02)

0x00 & MSB

runhdr = _read_varint(io, Int)
isbitpack = ((runhdr & 0x1) == 0x1)
runhdr >>= 1
nitems = isbitpack ? min(runhdr*8, num_values - arrpos + 1) : runhdr
#@debug("nitems=$nitems, isbitpack:$isbitpack, runhdr:$runhdr, remaininglen: $(num_values - arrpos + 1)")
#@debug("creating sub array for $nitems items at $arrpos, total length:$(length(arr))")
#subarr = pointer_to_array(pointer(arr, arrpos), nitems)
subarr = unsafe_wrap(Array, pointer(arr, arrpos), nitems, own=false)
isbitpack
position(io)
if isbitpack
	read_bitpacked_run(io, runhdr, bits, byt, typ, subarr)
else # rle
	read_rle_run(io, runhdr, bits, byt, typ, subarr)
end

runhdr
subarr

position(io)
read_rle_run(io, runhdr, bits, byt, typ, subarr)

arrpos += nitems
#end
arr

# end

		################################################################################
		# begin: expansion of functions `_read_varint`
		################################################################################
res = zero(T)
n = 0
byte = UInt8(MSB)
while (byte & MSB) != 0
	global res, n
    byte = read(io, UInt8)
    res |= (convert(T, byte & MASK7) << (7*n))
    n += 1
end
# in case of overflow, consider it as missing field and return default value
if (n-1) > sizeof(T)
    #@debug("overflow reading $T. returning 0")
    return zero(T)
end
res
		################################################################################
		# end: expansion of functions `_read_varint`
		################################################################################

while arrpos <= count
	runhdr = _read_varint(io, Int)
	isbitpack = ((runhdr & 0x1) == 0x1)
	runhdr >>= 1
	nitems = isbitpack ? min(runhdr*8, count - arrpos + 1) : runhdr
	#@debug("nitems=$nitems, isbitpack:$isbitpack, runhdr:$runhdr, remaininglen: $(count - arrpos + 1)")
	#@debug("creating sub array for $nitems items at $arrpos, total length:$(length(arr))")
	#subarr = pointer_to_array(pointer(arr, arrpos), nitems)
	subarr = unsafe_wrap(Array, pointer(arr, arrpos), nitems, own=false)

	if isbitpack
		read_bitpacked_run(io, runhdr, bits, byt, typ, subarr)
	else # rle
		read_rle_run(io, runhdr, bits, byt, typ, subarr)
	end
	arrpos += nitems
end
arr

		################################################################################
		# end: expansion of functions `read_values`
		################################################################################

	################################################################################
	# end: expansion of functions `read_levels_and_values`
	################################################################################if (typ == PageType.DATA_PAGE) || (typ == PageType.DATA_PAGE_V2)

# elseif typ == PageType.DICTIONARY_PAGE
# 	(read_plain_values(io, num_values, ctype),)
# else
# 	()
# end
################################################################################
# begin: expansion of functions `values`
################################################################################

cname = "Cl.thickness"
@which max_definition_level(par.schema, cname)
@which max_definition_level(par.schema, cname)

parentname("Cl.thickness")

@which isrequired(par.schema, "Cl.thickness")

schname = "Cl.thickness"
schelem = elem(par.schema, schname)

using Thrift
Thrift.isfilled(schelem, :repetition_type)

@which parentname(schname)

join(parentname(), '.')

parentname(split(schname, '.'))

(schelem.repetition_type == PAR2.FieldRepetitionType.REQUIRED)

# everything after the first data datapages
for data_page in Base.Iterators.drop(pages, 1)
	values, repetition, decode = values(par, data_page)
	l = sum(repetition)

	# if all repetition values are 1 then it's not used
	repetition_not_used = all(==(1), repetition)

	# data_page can be either
	# * dictionary-encoded in which we should look into the dictionary
	# * plained-encoded in which case just return the values
	page_encoding = Diban.page_encoding(data_page)

	if page_encoding == Encoding.PLAIN_DICTIONARY
		if repetition_not_used
			res[write_cursor:write_cursor+l-1] .= dictionary_of_values[values.+1]
		else
			for (offset, (repetition, value))  in enumerate(zip(repetition, values))
				if repetition != 0
					res[write_cursor+offset-1] = dictionary_of_values[value.+1]
				end
			end
		end
	elseif page_encoding == Encoding.PLAIN
		if repetition_not_used
			res[write_cursor:write_cursor+l-1] .= values
		else
			for (offset, (repetition, value))  in enumerate(zip(repetition, values))
				if repetition != 0
					res[write_cursor+offset-1] = value
				end
			end
		end
	else
		error("page encoding not supported yet")
	end

	write_cursor += l
end
#end
return res


text: struct _UnimplementedField <: Field end


text:
export absolute_genebodycoverage

"""
    absolute_genebodycoverage(path_bam::String, path_bed12::String;
							output_prefix::String = "", bin_size::Int=100)

Calculates absolute gene body coverage.

# Arguments
---------
- `path_bam::String`: Path to a BAM file.
- `path_bed12::String`: Path to a reference gene model file with BED12 format.
- `output_prefix::String`: Prefix for output files. If this keyword is not specified or set to "" (defalut), no output files are saved.
- `bin_size::Int`: Bin size for transcripts (Defalut: 100).
"""
function absolute_genebodycoverage(path_bam::String, path_bed12::String;
									output_prefix::String = "", bin_size::Int=100)
	if output_prefix != ""
		path_out = output_prefix * ".absoluteGeneBodyCoverage.txt"
		path_out_plot = output_prefix * ".absoluteGeneBodyCoverage.pdf"
		println(@sprintf "Start calculating absolute gene body coverage...\n- bam: %s\n- bed12: %s\n- output: %s\n          %s\n- bin_size: %d\n" path_bam path_bed12 path_out path_out_plot bin_size)
	end

	# File check
	if !isfile(path_bam)
		error(@sprintf "No such file: %s\n" path_bam)
	end
	if !isfile(path_bed12)
		error(@sprintf "No such file: %s\n" path_bed12)
	end
	if output_prefix != ""
		if !uv_access_writable(dirname(output_prefix))
			error(@sprintf "Output files are not writable to: %s\n" (dirname(output_prefix) != "" ? dirname(output_prefix) : "."))
		end
	end

	# Load transcripts information (BED12)
	transcripts = load_transcript(path_bed12)
	println(@sprintf "transcripts: %d" length(transcripts))

	# Calculate L_max (Maximum length of transcript)
	L_max = maximum(sum(BED.blocksizes(t)) for t in transcripts)
	println(@sprintf("L_max: %d", L_max))

	# Calculate K
	K = Int(ceil(L_max/bin_size))
	println(@sprintf("K: %d", K))

	# Calculate the number of transcripts in each bin
	N_transcript = zeros(Int, K)

	# Generate a L_max-length array for read coverage
	coverage = zeros(L_max)

	# Generate a K-length array for absolute coverage
	abcov = zeros(K)

	# Open BAM file
	open(BAM.Reader, path_bam, index=path_bam*".bai") do bam_reader

		# For each transcript
		println("Calculate absolute read coverage for each transcript...")

		for t in transcripts
			# Skip the transcript whose region does not exist in BAM
			if ! exists_region_bam(bam_reader, BED.chrom(t), BED.chromstart(t), BED.chromend(t))
				continue
			end

			# Get blockSizes (= exon lengths) and blockStarts (=start position of exons relative to chromStart of t)
			blockSizes = BED.blocksizes(t)
			L = sum(blockSizes)
			blockStarts = BED.blockstarts(t)

			# Count the number of transcripts whose 5'-end in each bin
			N_transcript[Int(ceil(sum(blockSizes)/bin_size))] += 1

			# Strand
			BED.strand(t)

			# Initialize coverage
			for p in 1:length(coverage)
				coverage[p] = 0
			end

			# Get the read coverage on each exon and save it to `coverage`
			# Note: Read coverage is saved in `coverage` in a 'Left justified' manner,
			#       and left is 5'-end.
			# If the strand of a transcript is '-' STRAND_NEG, read coverage was reversed.
			p_start = 1
			p_end = 0
			exon_start = BED.chromstart(t) # Returns 1-based coordinate(though 0-based in BED12 format)
			exon_end = 0
			for i in 1:length(blockSizes) # For each exon
				# Calculate start and end positions of an exon on the genome
				exon_start = blockStarts[i] - 1 + exon_start  # Returns 1-based coordinate(though 0-based in BED12 format)
				exon_end = exon_start + blockSizes[i] - 1

				# Calculate start and end positions of an exon on mRNA
				p_start = p_end + 1
				p_end = p_start + blockSizes[i] - 1

				# Get read coverage and save it
				if BED.strand(t) == STRAND_POS
					pos_cov = p_start:p_end
				elseif BED.strand(t) == STRAND_NEG
					pos_cov = (L-p_start+1):-1:(L-p_end+1)
				else # STRAND_NA, STRAND_BOTH
					continue
				end

				coverage[pos_cov] = readcoverage_bam_base(bam_reader, BED.chrom(t), exon_start, exon_end)
			end

			# Calculate count per bin
			# (Distance from 3'-end is binned.)
			# (k-th bin: [(k-1)*b+1, k*b])
			for k in 1:K # for each bin
				if k*bin_size > L # Case when k-th bin exceeds boundary of transcript
					if sum(coverage[1:(L-((k-1)*bin_size+1)+1)]) > 0
						abcov[k] += 1
						break
					end
				else
					if sum(coverage[(L-k*bin_size+1):(L-((k-1)*bin_size+1)+1)]) > 0
						abcov[k] += 1
					end
				end
			end
		end
	end

	# Calculate number of transcripts in each bin
	for k in (K-1):-1:1
		N_transcript[k] += N_transcript[k+1]
	end

	# Divide count by the number of transcripts for each bin
	for k in 1:K
		abcov[k] = abcov[k] / N_transcript[k]
	end

	if output_prefix != ""
		# Write output (TSV) (Sample_ID, Bin, Coverage)
		println("Write results...")
		bin_numbers = collect(1:K)
		bin_starts = (bin_numbers .- 1) .* bin_size .+ 1
		bin_ends = bin_numbers .* bin_size
		open(path_out, "w") do io
			writedlm(io, [bin_numbers bin_starts bin_ends abcov N_transcript], '\t')
		end

		# Save plot
		plot_absolute_coverage(abcov, out_path=path_out_plot)

		# Message
		println(@sprintf "Finished! Check output files:\n- %s\n- %s" path_out path_out_plot)
	end

    return(abcov)
end


text:
using FillArrays

"""
    annular(alg, cube, angles; fwhm, ann_size=4, init=0, nframes=4, delta_rot=1)
    annular(algs::AbstractVector, cube, angles; fwhm, ann_size=4, init=0, nframes=4, delta_rot=1)

"""
function annular(alg::ADIAlgorithm, cube, angles; kwargs...)
    n, y, x = size(cube)
    n_annuli = round(Int, (y / 2 - init) / ann_size)
    annular(Fill(alg, n_annuli), cube, angles; kwargs...)
end

function annular(algs::AbstractVector, cube, angles; fwhm, ann_size=4, init=0, min_frames=2, delta_rot=1)
    n, y, x = size(cube)
    n_annuli = round(Int, (y / 2 - init) / ann_size)
    # TODO
    # if n_segments === :auto
    #     ld = 2tan(45)
    #     n_segments = [2, 3, ]
    # end
    @info "# Annuli: $n_annuli, FWHM: $fwhm"
    Threads.@threads for ann_idx in eachindex(algs)
        pa_thr, inner_rad, ann_center = define_annuli(angles, ann_idx, n_annuli, fwhm, init, ann_size, delta_rot)
        mask = get_annulus_segments(view(cube, 1, :, :), inner_rad, ann_size)
        indices = findall(mask)
        matrix = cube[:, indices]
        out = similar(matrix)
        res = alg(matrix, angles)
    end
end

# # Optimized median subtraction for a given annulus
function medsub_annular(cube::AbstractArray{T,3}, angles::AbstractVector; fwhm, ann_size=4, radius_int=0, nframes=4, delta_sep=(0.1, 1), delta_rot=1, kwargs...) where T
    angles_ = normalize_par_angles(angles)
    med = median(cube, dims=1)
    S = cube .- med
    n, y, x = size(cube)

    n_annuli = round(Int, (y / 2 - radius_int) / ann_size)

    # _median_subt_ann_adi
    ann = 0:n_annuli - 1
    pa_thresh, inner_r, _ = define_annuli(angles, ann, n_annuli, fwhm, radius_int, ann_size, delta_rot)
    mask = get_annulus_segments(@view(S[1, :, :]), first(inner_r), ann_size)
    indices = findall(mask)
    matrix = S[:, indices]
    out = similar(matrix)

    for frame in 1:n
        if !iszero(pa_thresh)
            inds_left = find_indices_adi.((angles,), frame, pa_thresh, nframes)
            @show size(pa_thresh)
            @show size(inds_left)
            matrix_disc = @view matrix[inds_left, :, :]
        else
            matrix_disc = matrix
        end
        ref_psf_opt = median(matrix_disc, dims=1)
        curr_frame = @view matrix[frame, :, :]
        @. out[frame, :, :] = curr_frame - ref_psf_opt
    end
    # return out, yy, xx, pa_thresh
    # end _median_subt_ann_adi

    for  ann in 1:n_annuli
        out[:, yy[ann], xx[ann]] .= out[ann, :, :]
    end

end

function find_indices_adi(angles, frame, pa_thresh, nframes)
    n = length(angles)

    index_prev = findfirst(1:frame - 1) do i
        abs(angles[frame] - angles[i]) < pa_thresh
    end
    isnothing(index_prev) && (index_prev = frame - 1)

    index_foll = findfirst(frame:n) do k
        abs(angles[k] - angles[frame]) > pa_thresh
    end
    isnothing(index_foll) && (index_foll = n)

    # nframesi s not None (median annular)
    window = nframes ÷ 2

    ind1 = max(index_prev - window, 1)
    ind2 = index_prev - 1
    ind3 = index_foll
    ind4 = min(index_foll + window, n) - 1
    indices = vcat(ind1:ind2, ind3:ind4)
end

"""
Function that defines the annuli geometry using the input parameters.
Returns the parallactic angle threshold, the inner radius and the annulus
center for each annulus.
"""
function define_annuli(angles, ann, N, fwhm, radius_int, width, delta_rot)
    inner_rad = @. ifelse(ann == N - 1, radius_int + (ann * width - 1), radius_int + ann * width)
    ann_center = @. inner_rad + (width / 2)
    thresh = @. rad2deg(2 * atan(delta_rot * fwhm / (2 * ann_center)))
    mid_range = abs(maximum(angles) - minimum(angles)) / 2
    @. thresh[thresh ≥ 0.9 * mid_range] = 0.9 * mid_range
    return thresh, inner_rad, ann_center
end


text:
#= -------------------------------------------------------------------------
# @file runtests.jl
#
# @date 08/07/19 08:25:17
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
# Testing
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
using Test

# NOTE(elsuizo:2019-08-07): desde aqui vamos a llamar un archivo de test por cada
# archivo que haya en el src
@testset "Robotics" begin
   include("testing_transformations.jl")
end

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot


text:
using DelimitedFiles, ProgressMeter, PyPlot

include("./util.jl")
include("./parsers.jl")

function scenarioCoverage(topology,
                          iterations,
                          cutoffs; shape=.8,
                          scale=.0001,
                          weibull=true,
                          plot=true,
                          dirname="./data/raw/scenario_coverage")
    y_vals = []
    weibull_probs = []
    labels = []
    optimal_vals = []

    # COMUTE SCENARIOS
    links, capacity, link_probs, nodes = readTopology(topology)
    for j in 1:iterations
        if weibull
            probs = weibullProbs(length(link_probs), shape=shape, scale=scale)
        else
            probs = link_probs
        end
        push!(weibull_probs, probs)
        scenarios, probabilities = subScenariosRecursion(probs, sum(probs)/length(probs)^2)
        push!(optimal_vals, sum(probabilities))
    end

    # COMPUTE COVERAGE
    progress = ProgressMeter.Progress(length(cutoffs)*iterations, .1, "Computing Scenario Coverage...", 50)
    for i in 1:length(cutoffs)
        vals = []
        for j in 1:iterations
            scenarios, probabilities = subScenariosRecursion(weibull_probs[j], cutoffs[i])
            push!(vals, sum(probabilities))
            ProgressMeter.next!(progress, showvalues = [(:cutoff, cutoffs[i]), (:iteration, "$(j)/$(iterations)"), (:coverage, sum(probabilities))])
        end
        push!(y_vals, sum(vals)/length(vals))
        push!(labels, cutoffs[i])
    end
    push!(labels, "optimal")
    push!(y_vals, sum(optimal_vals)/length(optimal_vals))

    # LOG RESULTS
    dir = nextRun(dirname)
    writedlm("$(dir)/weibull_probs", weibull_probs)
    writedlm("$(dir)/coverage", y_vals)
    writedlm("$(dir)/cutoffs", labels)
    writedlm("$(dir)/params", [["iterations"], [iterations]])

    # PLOT
    if plot
        PyPlot.clf()
        nbars = length(labels)
        barWidth = 1/(nbars)
        for i in 1:length(y_vals)
            PyPlot.bar((barWidth .* i), y_vals[i], alpha=0.8, width=barWidth)
        end
        PyPlot.ylabel("Total Coverage", fontweight="bold")
        PyPlot.legend(labels, loc="lower right")
        PyPlot.savefig("$(dir)/plot.png")
        PyPlot.show()
    end
end

# scenarioCoverage("./Data/IBM", 100, [.001,.0001,.00001,.000001,.0000001]; shape=.8, scale=.0001, weibull=true)


text:
module RestrictedTriangulationsTests

using Gridap.Arrays
using Gridap.ReferenceFEs
using Gridap.Geometry
using Test

domain = (-1,1,-1,1)
partition = (10,10)
oldmodel = CartesianDiscreteModel(domain,partition)
oldtrian = get_triangulation(oldmodel)

cell_to_oldcell = collect(1:34)
trian = RestrictedTriangulation(oldtrian,cell_to_oldcell)
test_triangulation(trian)

cell_to_oldcell = [2,9,7]
trian_portion = RestrictedTriangulation(oldtrian,cell_to_oldcell)
trian_portion_portion = RestrictedTriangulation(trian_portion,[3,1])
test_triangulation(trian_portion_portion)
@test get_cell_to_bgcell(trian_portion_portion)==[7,2]

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

oldcell_to_coods = get_cell_coordinates(oldtrian)
oldcell_to_mask = lazy_map(is_in,oldcell_to_coods)
trian = RestrictedTriangulation(oldtrian,oldcell_to_mask)
trian = Triangulation(oldmodel,oldcell_to_mask)
trian = Triangulation(oldtrian,oldcell_to_mask)
trian = Triangulation(oldmodel,tags="interior")
trian = Triangulation(oldmodel,get_face_labeling(oldmodel),tags="interior")

#using Gridap.Visualization
#writevtk(oldtrian,"oldtrian")
#writevtk(btrian,"btrian",cellfields=["normal"=>nb],celldata=["oldcell"=>get_cell_to_bgcell(btrian)])
#writevtk(strian,"strian",cellfields=["normal"=>ns],
#  celldata=["oldcell_left"=>get_cell_to_bgcell(strian).plus,"oldcell_right"=>get_cell_to_bgcell(strian).minus])

end # module


text:
using Documenter
using MeasureTheory

DocMeta.setdocmeta!(MeasureBase, :DocTestSetup, :(using MeasureBase); recursive=true)
DocMeta.setdocmeta!(MeasureTheory, :DocTestSetup, :(using MeasureTheory); recursive=true)

pages = [
    "Home" => "index.md",
    "Tutorials" => [
        "Adding a new measure" => "adding.md",
        "Affine transformations" => "affine.md",
    ],
    "API Reference" => [
        "MeasureBase" => "api_measurebase.md",
        "MeasureTheory" => "api_measuretheory.md",
        "Index" => "api_index.md",
    ],
]

makedocs(;
    modules=[MeasureBase, MeasureTheory],
    authors="Chad Scherrer <chad.scherrer@gmail.com> and contributors",
    repo="https://github.com/cscherrer/MeasureTheory.jl/blob/{commit}{path}#L{line}",
    sitename="MeasureTheory.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cscherrer.github.io/MeasureTheory.jl",
        assets=String[],
    ),
    pages=pages
)

deploydocs(;
    repo="github.com/cscherrer/MeasureTheory.jl",
)