text:
@testset "indexing" begin
    @on_device threadIdx().x
    @on_device blockDim().x
    @on_device blockIdx().x
    @on_device gridDim().x

    @on_device threadIdx().y
    @on_device blockDim().y
    @on_device blockIdx().y
    @on_device gridDim().y

    @on_device threadIdx().z
    @on_device blockDim().z
    @on_device blockIdx().z
    @on_device gridDim().z

    @on_device warpsize()
    @on_device laneid()
    @on_device active_mask()

    @testset "range metadata" begin
        foobar() = threadIdx().x
        ir = sprint(io->CUDA.code_llvm(io, foobar, Tuple{}; raw=true))

        @test occursin(r"call .+ @llvm.nvvm.read.ptx.sreg.tid.x.+ !range", ir)
    end
end

############################################################################################

@testset "assertion" begin
    function kernel(i)
        @cuassert i > 0
        @cuassert i > 0 "test"
        return
    end

    @cuda kernel(1)
end

############################################################################################

@testset "data movement and conversion" begin

if capability(device()) >= v"3.0"

@testset "shuffle idx" begin
    function kernel(d)
        i = threadIdx().x
        j = 32 - i + 1

        d[i] = shfl_sync(FULL_MASK, d[i], j)

        return
    end

    warpsize = CUDA.warpsize(device())

    @testset for T in [UInt8, UInt16, UInt32, UInt64, UInt128,
                       Int8, Int16, Int32, Int64, Int128,
                       Float16, Float32, Float64,
                       ComplexF32, ComplexF64, Bool]
        a = rand(T, warpsize)
        d_a = CuArray(a)
        @cuda threads=warpsize kernel(d_a)
        @test Array(d_a) == reverse(a)
    end
end

@testset "shuffle down" begin
    n = 14

    function kernel(d::CuDeviceArray, n)
        t = threadIdx().x
        if t <= n
            d[t] += shfl_down_sync(FULL_MASK, d[t], n÷2, 32)
        end
        return
    end

    @testset for T in [Int32, Float32]
        a = T[T(i) for i in 1:n]
        d_a = CuArray(a)

        threads = nextwarp(device(), n)
        @cuda threads=threads kernel(d_a, n)

        a[1:n÷2] += a[n÷2+1:end]
        @test a == Array(d_a)
    end
end

end

end

############################################################################################

@testset "clock and nanosleep" begin
@on_device clock(UInt32)
@on_device clock(UInt64)

if CUDA.release() >= v"10.0" && v"6.2" in CUDA.ptx_support()
    @on_device nanosleep(UInt32(16))
end
end

@testset "parallel synchronization and communication" begin

@on_device sync_threads()
@on_device sync_threads_count(Int32(1))
@on_device sync_threads_count(true)
@on_device sync_threads_and(Int32(1))
@on_device sync_threads_and(true)
@on_device sync_threads_or(Int32(1))
@on_device sync_threads_or(true)

@testset "llvm ir barrier int" begin
    function kernel_barrier_count(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_count(an_array_of_1[i]) == 3
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    b_in = Int32[1, 1, 1, 0, 0]  # 3 true
    d_b = CuArray(b_in)
    @cuda threads=5 kernel_barrier_count(d_b)
    b_out = Array(d_b)

    @test b_out == b_in .+ 1

    function kernel_barrier_and(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_and(an_array_of_1[i]) > 0
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    a_in = Int32[1, 1, 1, 1, 1]  # all true
    d_a = CuArray(a_in)
    @cuda threads=5 kernel_barrier_and(d_a)
    a_out = Array(d_a)

    @test a_out == a_in .+ 1

    function kernel_barrier_or(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_or(an_array_of_1[i]) > 0
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    c_in = Int32[1, 0, 0, 0, 0]  # 1 true
    d_c = CuArray(c_in)
    @cuda threads=5 kernel_barrier_or(d_c)
    c_out = Array(d_c)

    @test c_out == c_in .+ 1
end

@testset "llvm ir barrier bool" begin
    function kernel_barrier_count(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_count(an_array_of_1[i] > 0) == 3
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    b_in = Int32[1, 1, 1, 0, 0]  # 3 true
    d_b = CuArray(b_in)
    @cuda threads=5 kernel_barrier_count(d_b)
    b_out = Array(d_b)

    @test b_out == b_in .+ 1

    function kernel_barrier_and(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_and(an_array_of_1[i] > 0)
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    a_in = Int32[1, 1, 1, 1, 1]  # all true
    d_a = CuArray(a_in)
    @cuda threads=5 kernel_barrier_and(d_a)
    a_out = Array(d_a)

    @test a_out == a_in .+ 1

    function kernel_barrier_or(an_array_of_1)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        if sync_threads_or(an_array_of_1[i] > 0)
            an_array_of_1[i] += Int32(1)
        end
        return nothing
    end

    c_in = Int32[1, 0, 0, 0, 0]  # 1 true
    d_c = CuArray(c_in)
    @cuda threads=5 kernel_barrier_or(d_c)
    c_out = Array(d_c)

    @test c_out == c_in .+ 1
end

@on_device sync_warp()
@on_device sync_warp(0xffffffff)
@on_device threadfence_block()
@on_device threadfence()
@on_device threadfence_system()

@testset "voting" begin

@testset "any" begin
    d_a = CuArray([false])

    function kernel(a, i)
        vote = vote_any_sync(FULL_MASK, threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel(d_a, 1)
    @test Array(d_a)[]

    @cuda threads=2 kernel(d_a, 2)
    @test Array(d_a)[]

    @cuda threads=2 kernel(d_a, 3)
    @test !Array(d_a)[]
end

@testset "all" begin
    d_a = CuArray([false])

    function kernel(a, i)
        vote = vote_all_sync(FULL_MASK, threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel(d_a, 1)
    @test Array(d_a)[]

    @cuda threads=2 kernel(d_a, 2)
    @test !Array(d_a)[]

    @cuda threads=2 kernel(d_a, 3)
    @test !Array(d_a)[]
end

@testset "uni" begin
    d_a = CuArray([false])

    function kernel1(a, i)
        vote = vote_uni_sync(FULL_MASK, threadIdx().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel1(d_a, 1)
    @test Array(d_a)[]

    @cuda threads=2 kernel1(d_a, 2)
    @test !Array(d_a)[]

    function kernel2(a, i)
        vote = vote_uni_sync(FULL_MASK, blockDim().x >= i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    @cuda threads=2 kernel2(d_a, 3)
    @test Array(d_a)[]
end

@testset "ballot" begin
    d_a = CuArray(UInt32[0])

    function kernel(a, i)
        vote = vote_ballot_sync(FULL_MASK, threadIdx().x == i)
        if threadIdx().x == 1
            a[1] = vote
        end
        return
    end

    len = 4
    for i in 1:len
        @cuda threads=len kernel(d_a, i)
        @test Array(d_a) == [2^(i-1)]
    end
end

end

end

############################################################################################

@testset "libcudadevrt" begin
    kernel() = (CUDA.device_synchronize(); nothing)
    @cuda kernel()
end

############################################################################################


text:
using Base.Test

tests = Base.Test.DefaultTestSet("all")
Base.Test.push_testset(tests)

dir = dirname(@__FILE__)
for f in readdir(dir)
    m = match(r"^test_(.*)\.jl$", f)
    fp = joinpath(dir, f)
    if m === nothing || !isfile(fp)
        continue
    end

    name = m.captures[1]
    @printf("%-12s ... ", name)

    ts = Base.Test.DefaultTestSet(name)
    Base.Test.push_testset(ts)
    try
        evalfile(fp)
    catch err
        res = Base.Test.Error(:nontest_error, :(), err, catch_backtrace())
        Base.Test.record(ts, res)
    end
    Base.Test.pop_testset()
    Base.Test.finish(ts)

    # lazy: use get_test_counts to set ts.anynonpass
    Base.Test.get_test_counts(ts)
    ts.anynonpass || println("OK")
end

Base.Test.pop_testset()
Base.Test.finish(tests)

ext: jl


text:
using CodedComputing, MPI, MPIStragglers
using HDF5, LinearAlgebra, SparseArrays
using ArgParse
using MKLSparse # important for performance of sparse matrix multiplication on Intel-based systems

MPI.Init()
const comm = MPI.COMM_WORLD
const root = 0
const rank = MPI.Comm_rank(comm)
const isroot = MPI.Comm_rank(comm) == root
const data_tag = 0
const control_tag = 1

"""

Parse command line arguments
"""
function parse_commandline(isroot::Bool)

    # setup error handling so that only the root node prints error messages
    s = ArgParseSettings("Latency tracing kernel", autofix_names=true)
    function root_handler(settings::ArgParseSettings, err, err_code::Integer=1)
        println(Base.stderr, err.text)
        println(Base.stderr, usage_string(settings))
        MPI.Finalize()
        exit(err_code)
    end
    function worker_handler(settings::ArgParseSettings, err)
        MPI.Finalize()
        exit(0)
    end
    if isroot
        s.exc_handler = root_handler
    else
        s.exc_handler = worker_handler
        s.add_help = false
    end

    # command-line arguments common to all pca implementations
    @add_arg_table s begin
        "outputfile"
            help = "HFD5 file to write the output to"
            required = true
            arg_type = String
            range_tester = (x) -> !isfile(x) || HDF5.ishdf5(x)
        "--niterations"
            help = "Number of iterations to run"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--nbytes"
            help = "Number of bytes communicated in each direction per iteration"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 8 # to fit a Float64
        "--nrows"
            help = "Number of rows of the data matrix"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--ncols"
            help = "Number of columns of the data matrix"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--ncomponents"
            help = "Number of columns of the iterate matrix"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--density"
            help = "Density of the data matrix"
            required = true
            arg_type = Float64
            range_tester = (x) -> 0 < x <= 1
        "--nwait"
            help = "Number of workers to wait for in each iteration"
            required = true
            arg_type = Int
            range_tester = (x) -> x >= 1
        "--timeout"
            help = "Amount of time to sleep between iterations"
            required = true
            arg_type = Float64
            range_tester = (x) -> 0 <= x
    end

    # optionally add implementation-specific arguments
    if @isdefined update_argsettings!
        update_argsettings!(s)
    end

    # common parsing
    parsed_args = parse_args(s, as_symbols=true)

    # optional implementation-specific parsing
    if @isdefined update_parsed_args!
        update_parsed_args!(s, parsed_args)
    end

    return parsed_args
end

function worker_task!(V, W, X)
    mul!(W, X, V)
    mul!(V, X', W)
end

"""

Receive from the coordinator, perform a matrix-matrix multiplication, and respond with some bytes.
"""
function worker_main()
    parsed_args = parse_commandline(isroot)
    nworkers = MPI.Comm_size(comm) - 1
    parsed_args[:nworkers] = nworkers

    # communication setup
    recvbuf = zeros(UInt8, parsed_args[:nbytes])
    sendbuf = zeros(UInt8, parsed_args[:nbytes])
    tbuf = reinterpret(Float64, view(sendbuf, 1:8))

    # worker task setup
    X = sprand(parsed_args[:nrows], parsed_args[:ncols], parsed_args[:density])
    V = randn(parsed_args[:ncols], parsed_args[:ncomponents])
    W = zeros(parsed_args[:nrows], parsed_args[:ncomponents])

    # control channel, to tell the workers when to exit
    crreq = MPI.Irecv!(zeros(1), root, control_tag, comm)

    # main loop
    while true
        rreq = MPI.Irecv!(recvbuf, root, data_tag, comm)
        index, _ = MPI.Waitany!([crreq, rreq])
        if index == 1 # exit message on control channel
            break
        end
        tbuf[1] = @elapsed worker_task!(V, W, X)
        MPI.Isend(sendbuf, root, data_tag, comm)
    end
    return
end

function shutdown(pool::AsyncPool)
    for i in pool.ranks
        MPI.Isend(zeros(1), i, control_tag, comm)
    end
end

"""

Send some bytes to the workers, wait for their response, sleep for a bit, and repeat.
"""
function coordinator_main()
    parsed_args = parse_commandline(isroot)
    nworkers = MPI.Comm_size(comm) - 1
    parsed_args[:nworkers] = nworkers
    niterations::Int = parsed_args[:niterations]
    niterations > 0 || throw(DomainError(niterations, "The number of iterations must be non-negative"))
    nwait::Int = isnothing(parsed_args[:nwait]) ? nworkers : parsed_args[:nwait]
    0 < nwait <= nworkers || throw(DomainError(nwait, "nwait must be in [1, nworkers]"))
    timeout::Float64 = parsed_args[:timeout]

    # communication setup
    pool = AsyncPool(nworkers, epoch0=-1)
    recvbuf = zeros(UInt8, nworkers*parsed_args[:nbytes])
    sendbuf = zeros(UInt8, parsed_args[:nbytes])
    irecvbuf = similar(recvbuf)
    isendbuf = similar(sendbuf, nworkers*length(sendbuf))

    # views into recvbuf for each worker
    n = div(length(recvbuf), nworkers)
    recvbufs = [view(recvbuf, (i-1)*n+1:i*n) for i in 1:nworkers]
    tbufs = [reinterpret(Float64, view(buf, 1:8)) for buf in recvbufs]

    # store which workers responded in each iteration
    worker_repochs = zeros(Int, nworkers, niterations) # receive epoch for each worker and iteration
    worker_latency = zeros(nworkers, niterations) # latency of individual workers
    worker_compute_latency = zeros(nworkers, niterations) # compute latency recorded by the worker
    timestamps = zeros(UInt64, niterations)
    latency = zeros(niterations)

    # dummy iteration to force compilation
    timestamps[1] = time_ns()
    latency[1] = @elapsed begin
        repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm, nwait=nworkers, tag=data_tag)
    end
    worker_repochs[:, 1] .= repochs
    worker_latency[:, 1] .= pool.latency
    if timeout > 0
        sleep(timeout)
    end

    # main loop
    for i in 1:parsed_args[:niterations]
        timestamps[i] = time_ns()
        latency[i] = @elapsed begin
            repochs = asyncmap!(pool, sendbuf, recvbuf, isendbuf, irecvbuf, comm, nwait=nwait, tag=data_tag)
        end
        worker_repochs[:, i] .= repochs
        worker_latency[:, i] .= pool.latency
        for j in 1:nworkers
            worker_compute_latency[j, i] = tbufs[j][1]
        end
        if timeout > 0
            sleep(timeout)
        end
    end
    shutdown(pool)

    # write statistics to file
    h5open(parsed_args[:outputfile], "w") do fid
        for (key, val) in parsed_args
            if isnothing(val)
                continue
            end
            fid["parameters/$key"] = val
        end
        fid["worker_repochs"] = worker_repochs
        fid["worker_latency"] = worker_latency
        fid["worker_compute_latency"] = worker_compute_latency
        fid["timestamps"] = timestamps
        fid["latency"] = latency
    end
    return
end

if isroot
    coordinator_main()
else
    # try/catch to eliminate stacktrace
    try
        worker_main()
    catch e
        printstyled(stderr,"ERROR: ", bold=true, color=:red)
        printstyled(stderr,sprint(showerror,e), color=:light_red)
        println(stderr)
    end
end
MPI.Barrier(comm)


text:
function bisect(v, x, imin::Int, imax::Int)
    # https://docs.python.org/2/library/bisect.html
    # The returned insertion point i partitions the sorted vector v into two halves so that
    # all(val <= x for val in v[1:i]) for the left side and
    # all(val > x  for val in v[i:hi]) for the right side.
    if imin == imax
        x <= v[imin] ? imin : imin + 1
    else
        imid = div(imin+imax, 2)
        if x == v[imid]
            return imid
        elseif x > v[imid]
            bisect(v, x, imid + 1, imax)
        else
            bisect(v, x, imin, max(imid - 1, imin))
        end
    end
end
function bisect(v, x)
    bisect(v,x,1,length(v))
end

function genSamples(values, weights::Vector{Float64}, n::Int64)
    @assert length(values) == length(weights)
    if sum(weights) == 1.0
        probs = weights
    else
        probs = weights/sum(weights)
    end
    cdf = cumsum(probs)
    samples = []
    for i in 1:n
        u   = rand()
        i   = bisect(cdf, u)
        push!(samples, values[i])
    end
    return samples
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
"""
    abstract type Environment end

Supertype for all reinforcement learning (RL) environments. In
essence a subtype of an `Environment` is expected to behave
according to some underlying Markov Decision Process (MDP) that
only can only be interacted with in a stepwise fashion.

There are some notable conceptual differences between a subtype
of `Environment` and a subtype of `MDP`.

- While an environment is assumed to be governed by an MDP,
  this MDP is generally only partially known (emphasis on
  "known"; this doesn't imply it is only partially observable).
  The typical unknowns are the dynamics (transitions) and/or the
  reward function.

- Consequently an environment does not need to be able expose
  arbitrary state transitions. Instead, it can be stateful and
  only ever expose a sampled trajectory starting from initial
  states of its own choosing.

- Any subtype of an `MDP` can easily be converted to an
  `Environment` by keeping track of its own state. The other
  direction is not as simple.

Any subtype of `Environment` has to implement the following
functions: [`actions`](@ref), [`step!`](@ref), [`reset!`](@ref),
[`isterminal`](@ref).
"""
abstract type Environment end

"""
    step!(env::Environment, a) -> (s´, r)

Apply the given action `a` to the current state of the
environment `env` and return the next state `s´` as well as the
received award `r` caused by this transition.

Note that this function does not expect the user to provide a
state. This is because the environment is in general expected to
be stateful and thus knows its actual state. The state `s´`
observed from the environment is not required to have a causal
relationship to the next state. This does not imply that the
underlying MDP is therefor partially observed. For example the
next screen in an Atari game is not actually determined by
looking at the last 5 screens. Yet for Atari games it is still a
fair assumption that the last 5 screens correspond to a unique
game state.
"""
function step! end

"""
    reset!(env::Environment, [...]) -> s₀

Reset the given environment to an initial state and return it.
Note that individual environment may support additional
parameters.
"""
function reset! end

"""
    isterminal(env::Environment) -> Bool

Check is `env` is in a terminal state.
"""
isterminal(env::Environment) = false


text:
module GLAbstraction

using StaticArrays
using GeometryBasics
using ModernGL
using AbstractPlotting
using FixedPointNumbers
using ColorTypes
using FileIO
using ..GLMakie.GLFW
using Printf
using LinearAlgebra
using Observables
using ShaderAbstractions
using ShaderAbstractions: current_context, is_context_active, context_alive

import FileIO: load, save

import FixedPointNumbers: N0f8, N0f16, N0f8, Normed

import AbstractPlotting: update!

import Base: merge, resize!, similar, length, getindex, setindex!

include("AbstractGPUArray.jl")

#Methods which get overloaded by GLExtendedFunctions.jl:
import ModernGL.glShaderSource
import ModernGL.glGetAttachedShaders
import ModernGL.glGetActiveUniform
import ModernGL.glGetActiveAttrib
import ModernGL.glGetProgramiv
import ModernGL.glGetIntegerv
import ModernGL.glGenBuffers
import ModernGL.glGetProgramiv
import ModernGL.glGenVertexArrays
import ModernGL.glGenTextures
import ModernGL.glGenFramebuffers
import ModernGL.glGetTexLevelParameteriv
import ModernGL.glGenRenderbuffers
import ModernGL.glDeleteTextures
import ModernGL.glDeleteVertexArrays
import ModernGL.glDeleteBuffers
import ModernGL.glGetShaderiv
import ModernGL.glViewport
import ModernGL.glScissor

include("GLUtils.jl")
export @materialize #splats keywords from a dict into variables
export @materialize!  #splats keywords from a dict into variables and deletes them from the dict
export close_to_square
export AND, OR, isnotempty

include("shaderabstraction.jl")
include("GLTypes.jl")
export GLProgram                # Shader/program object
export Texture                  # Texture object, basically a 1/2/3D OpenGL data array
export TextureParameters
export TextureBuffer            # OpenGL texture buffer
export update!                  # updates a gpu array with a Julia array
export gpu_data                 # gets the data of a gpu array as a Julia Array

export RenderObject             # An object which holds all GPU handles and datastructes to ready for rendering by calling render(obj)
export prerender!               # adds a function to a RenderObject, which gets executed befor setting the OpenGL render state
export postrender!              # adds a function to a RenderObject, which gets executed after setting the OpenGL render states
export std_renderobject            # creates a renderobject with standard parameters
export instanced_renderobject    # simplification for creating a RenderObject which renders instances
export extract_renderable
export set_arg!
export GLVertexArray            # VertexArray wrapper object
export GLBuffer                 # OpenGL Buffer object wrapper
export indexbuffer              # Shortcut to create an OpenGL Buffer object for indexes (1D, cardinality of one and GL_ELEMENT_ARRAY_BUFFER set)
export opengl_compatible        # infers if a type is opengl compatible and returns stats like cardinality and eltype (will be deprecated)
export cardinality              # returns the cardinality of the elements of a buffer

export Style                    # Style Type, which is used to choose different visualization/editing styles via multiple dispatch
export mergedefault!            # merges a style dict via a given style
export TOrSignal, VecOrSignal, ArrayOrSignal, MatOrSignal, VolumeOrSignal, ArrayTypes, VectorTypes, MatTypes, VolumeTypes
export MouseButton, MOUSE_LEFT, MOUSE_MIDDLE, MOUSE_RIGHT

include("GLExtendedFunctions.jl")
export glTexImage # Julian wrapper for glTexImage1D, glTexImage2D, glTexImage3D
include("GLShader.jl")
export Shader                 #Shader Type
export readshader             #reads a shader
export glsl_variable_access # creates access string from julia variable for the use in glsl shaders
export createview #creates a view from a templated shader
export TemplateProgram # Creates a shader from a Mustache view and and a shader file, which uses mustache syntax to replace values.
export @comp_str #string macro for the different shader types.
export @frag_str # with them you can write frag""" ..... """, returning shader object
export @vert_str
export @geom_str
export AbstractLazyShader, LazyShader

include("GLUniforms.jl")
export gluniform                # wrapper of all the OpenGL gluniform functions, which call the correct gluniform function via multiple dispatch. Example: gluniform(location, x::Matrix4x4) = gluniformMatrix4fv(location, x)
export toglsltype_string        # infers a glsl type string from a julia type. Example: Matrix4x4 -> uniform mat4
# Also exports Macro generated GLSL alike aliases for Float32 Matrices and Vectors
# only difference to GLSL: first character is uppercase uppercase
export gl_convert

include("GLRender.jl")
export render  #renders arbitrary objects
export enabletransparency # can be pushed to an renderobject, enables transparency
export renderinstanced # renders objects instanced

include("GLInfo.jl")
export getUniformsInfo
export getProgramInfo
export getAttributesInfo

end # module


text:
"""
    ParseCSV(filename::String, header::Bool=true)

This function parses a CSV file containing input for the MCMC compuation.
The file should follow the conventions used for MrBayes. For example:

    Swedish_0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,?,0,0,?,0,0
    Welsh_N_0,0,0,0,0,0,0,?,0,0,0,0,?,?,0,0,?,0,0,0,1,?,?,0
    Sardinian_N_0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,?,0,0,0,0,0

 Set input for "header" to false if no header is present in the file.

Returns ntax, nchar as Integer values, gap as a String; df is returned as a DataFrame, and stores language names and data.

* `filename` : Name of CSV file to parse.

* `header` : Boolean value. "true" denotes there is a header to be skipped; input "false" if the file does not contain a header.
"""
function ParseCSV(filename::String, header::Bool=true)
    open(filename, "r") do file
        global content = readlines(file)
    end # do
    # remove the header if necessary
    header ||  popfirst!(content)
    df = create_csvdf(content)
    dimensions = size(df)
    ntax = dimensions[1]
    nchar = dimensions[2]
    symbols = get_alphabet(df)
    return ntax, nchar, gap, missing_representation, symbols, df
end # function ParseCSV

"""
    function create_csvdf(filecontnt::Array{String}, separator::AbstractString=",")::DataFrame

This function parses a CSV file and returns its content as a DataFrame.

NOTE: May later on be replaced by the respective DataFrames function.
"""
function create_csvdf(filecontent::Array{String}, separator::AbstractString=",")::DataFrame
    language = []
    data = []
    #df's created after editing the filecontent into separate language/data objects
    for line in filecontent
        splitline = split(line,separator)
        push!(language,splitline[1])
        push!(data,splitline[2:end])
    end #for
    df = DataFrame(Language = language,Data = data)
    return df
end # function create_csvdf


text:
using GooglePlaces

import ConfParser: ConfParse, parse_conf!, retrieve

conf = ConfParse("configs/config.ini")
parse_conf!(conf)

# get and store config parameters
appname = retrieve(conf, "app", "name")
appenv = retrieve(conf, "app", "env")
logtofile = retrieve(conf, "log", "log_to_file", Bool)
loglevel = retrieve(conf, "log", "level", Int32)
logpath = retrieve(conf, "log", "path")
key = retrieve(conf, "googleAPI", "key")

fields = ["name","rating","formatted_phone_number","address_component","opening_hours","geometry"]
placeid = "ChIJMQOS9VovmEcRTPG4ZBnGOi0"

GooglePlaces.initialise_logging(logtofile, logpath, loglevel, appenv)

df = GooglePlaces.placedetails(key=key, placeid=placeid, fields=fields)


text:
"""Abstract type of Devices that inject current/power """
abstract type StaticInjection <: Device end

function supports_services(::T) where {T <: Device}
    return true
end

function get_services(::Device)
    return Vector{Service}()
end

get_dynamic_injector(d::StaticInjection) = nothing


text:
using KLIEPInference
using JLD

scratch_dir = ARGS[1]

p = 0.05:0.05:0.95

for m in 25
    for sgn in 1
        for numChanges in [1, 3, 5]
            println("computing power for $(m) $(sgn) $(numChanges) ...")

            power = zeros(11, 4)
            for lbInd = 1:11
                file = jldopen("graphs/params_exp4_$(m)_$(sgn)_$(numChanges)_$(lbInd).jld", "r")
                γx = read(file, "γx")
                γy = read(file, "γy")
                θ = γx - γy
                close(file)

                coverage = zeros(length(p), 4)

                nrep = 0
                for rep in 1:1000
                    try
                        file = jldopen("$(scratch_dir)/res_$(m)_$(sgn)_$(numChanges)_$(lbInd)_$(rep).jld", "r")
                        boot1 = read(file, "boot1")
                        boot2 = read(file, "boot2")
                        close(file)

                        # coverage
                        q = boot_quantile(boot1, p)
                        coverage[:, 1] .+= maximum(abs.(boot1.θhat - θ)) .<= q

                        q = boot_quantile(boot2, p)
                        coverage[:, 2] .+=  maximum(abs.(boot2.θhat - θ)) .<= q

                        q, w = boot_quantile_studentized(boot1, p)
                        coverage[:, 3] .+= maximum(abs.((boot1.θhat - θ) ./ w)) .<= q

                        q, w = boot_quantile_studentized(boot2, p)
                        coverage[:, 4] .+= maximum(abs.((boot2.θhat - θ) ./ w)) .<= q

                        # power
                        q = boot_quantile(boot1, 0.95)
                        power[lbInd, 1] += maximum(abs.(boot1.θhat)) > q

                        q = boot_quantile(boot2, 0.95)
                        power[lbInd, 2] += maximum(abs.(boot2.θhat)) > q

                        q, w = boot_quantile_studentized(boot1, 0.95)
                        power[lbInd, 3] += maximum(abs.(boot1.θhat ./ w)) > q

                        q, w = boot_quantile_studentized(boot2, 0.95)
                        power[lbInd, 4] += maximum(abs.(boot2.θhat ./ w)) > q

                        nrep += 1
                    catch
                        @warn "could not open $(m) $(sgn) $(numChanges) $(lbInd) $(rep)"
                    end
                end

                coverage = broadcast(/, coverage, nrep)

                @save "res/coverage_$(m)_$(sgn)_$(numChanges)_$(lbInd).jld" coverage

                power[lbInd, :] ./= nrep
            end

            @save "res/power_$(m)_$(sgn)_$(numChanges).jld" power

            println("... done!")
        end
    end
end


text:
module FeatherLib

using Arrow, FlatBuffers, CategoricalArrays, Mmap

export featherread, featherwrite

import Dates

const FEATHER_VERSION = 2
# wesm/feather/cpp/src/common.h
const FEATHER_MAGIC_BYTES = Vector{UInt8}(codeunits("FEA1"))
const MIN_FILE_LENGTH = 12

include("metadata.jl")  # flatbuffer defintions
include("loadfile.jl")
include("read.jl")
include("write.jl")

end # module


text:
using Luxor

function visualize_state(file_name, state, num_pred, num_prey, step)
    Drawing(610, 610, file_name)
    background("white")
    origin()

    pred_color = "red"
    prey_color = "blue"
    agent_point_size = 8

    # List of grid index
    a = [2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 20, 21, 22, 23]

    # List of hex cell index
    b = [13, 12, 11, 14, 4, 3, 10, 15, 5, 1, 2, 9, 16, 6, 7, 8, 17, 18, 19, 20]

    # Map from grid index to hex cell index
    d = Dict([(i, j) for (i, j) in zip(a, b)])

    # Map from hex cell index to grid index
    rd = Dict(value => key for (key, value) in d)

    # Map from grid index to a tuple t, where
    #  t[1]: a list of predators at this grid index
    #  t[2]: a list of prey at this grid index
    agents_to_place = Dict([(i, ([], [])) for i=1:23])

    preds = state[1:num_pred]
    preys = state[(num_pred+1):(num_pred+num_prey)]

    for (pidx, p) in enumerate(preds)
        push!(agents_to_place[rd[p]][1], pidx)
    end

    for (pidx, p) in enumerate(preys)
        push!(agents_to_place[rd[p]][2], pidx)
    end

    radius = 50
    grid = GridHex(Point(-250, -125), radius, 100)
    for i in 1:23
        # Determine if this center point belong to the map
        hex_cell_idx = get(d, i, nothing)
        if hex_cell_idx !== nothing
            sethue("black")
        else
            # Don't draw the hex tile (i.e. color it white) if this is not in the map
            sethue("white")
        end

        # Get the center of the grid cell
        grid_point = nextgridpoint(grid)

        # Draw a hexagon at this center, rotate it 45 degree
        setline(3)
        ngon(grid_point, radius - 10, 6, pi/2, :stroke)

        if hex_cell_idx !== nothing
            # Get a list of agents to place at this cell
            (p, q) = agents_to_place[i]

            # Vector of predator colors
            pc = Array{Any}(undef, length(p))
            fill!(pc, pred_color)

            # Vector of prey colors
            qc = Array{Any}(undef, length(q))
            fill!(qc, prey_color)

            # Drawing multiple agents in the same grid cell
            agents = vcat(p, q)
            agent_colors = vcat(pc, qc)
            j = 1
            if !isempty(agents)
                # Place the agents at the vertices of a transparent polygon
                sethue("white")
                poly = ngon(grid_point, 10, length(agents), 0)
                prettypoly(
                    poly,
                    :fill,
                    () ->
                    begin
                        # Draw the agents as circles with their corresponding colors
                        sethue(agent_colors[j])
                        circle(O, agent_point_size, :fill)
                        # Label the index of the agents
                        sethue("black")
                        text(string(agents[j]), O, halign=:center, valign=:bottom)
                        j += 1
                    end,
                    close = true)
            end
        end
    end

    sethue("black")
    fontsize(20)
    text("Step $step", Point(-260, -200))
    finish()
end


text:
function main()
    # 点数
    N = 100000000

    # 刻み
    h = 0.00000004

    # 初期化
    y = zeros(N)

    # 初期値
    y[1] = 1

    # 数値計算
     for n = 1:N-1
        y[n+1] = (1-h)y[n]
    end
end
@time main()

max_stars_repo_name: terasakisatoshi/pythonCodes

max_issues_repo_name: terasakisatoshi/pythonCodes

max_forks_repo_name: terasakisatoshi/pythonCodes


text:
# load profiles

import JSON

struct LangProfile
    name::String
    n_words::Vector{Int}
    freq::Dict{String, Int}
end

function load_profiles(dir::String)
    profiles = Vector{LangProfile}()
    for name in readdir(dir)
        fpath = joinpath(dir, name)
        open(fpath) do reader
           data = JSON.parse(reader)
           push!(profiles, LangProfile(data["name"], data["n_words"], data["freq"]))
        end
    end
    return profiles
end


text:
list_test_str = split(
"1163751742
1381373672
2136511328
3694931569
7463417111
1319128137
1359912421
3125421639
1293138521
2311944581",
'\n'
)

list_str = let input = @__DIR__() * "/day15_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list)
    risk_map = parse.(Int, permutedims(reduce(hcat, collect.(list))))
    mem = Dict{NTuple{4, Int}, Tuple{Int, Matrix{Int}}}()
    h, w = size(risk_map)
    function risk!(mem, risk_map, n, i, m, j)
        (n, i, m, j) ∈ keys(mem) && return mem[(n, i, m, j)]
        if m == n + 1
            a = sum(risk_map[m, i:j])
            b = sum(risk_map[n, i+1:j]) + risk_map[m, j]
            risk_value = min(a, b)
            path = a ≤ b ? ( i < j ? [m i; m j] : [m j] ) : [n j; m j]
        else
            risks = Vector{Int}()
            paths = Vector{Matrix{Int}}()
            for k in i:j
                r1, p1 = risk!(mem, risk_map, n, i, m-1, k)
                r2, p2 = risk!(mem, risk_map, m-1, k, m, j)
                push!(risks, r1 + r2)
                push!(paths, [p1; p2])
            end
            k = argmin(risks)
            risk_value = risks[k]
            path = paths[k]
        end
        push!(mem, (n, i, m, j) => (risk_value, path))
        return risk_value, path
    end
    risk_value = risk!(mem, risk_map, 1, 1, h, w)
end

testit("Test Part 1", part1(list_test_str)[1], ==(40))
testit("Challenge Part 1", part1(list_str)[1], ==(447))

# Part 2

function part2(list, r)
    risk_map = parse.(Int, permutedims(reduce(hcat, collect.(list))))
    h, w = size(risk_map)
    risk_map = repeat(risk_map, r, r)
    for k in 2:r
        risk_map[1:h, (k-1)*w+1:k*w] .= mod1.(@view(risk_map[1:h, (k-2)*w+1:(k-1)*w]) .+ 1, 9)
    end
    for k in 2:r
        risk_map[(k-1)*h+1:k*h, :] .= mod1.(@view(risk_map[(k-2)*h+1:(k-1)*h, :]) .+ 1, 9)
    end
    h, w = size(risk_map)
    cum_risk = Matrix{Int}(undef, h, w)
    cum_risk[1, 1] = 0
    for j in 2:w
        cum_risk[1, j] = risk_map[1, j] + cum_risk[1, j-1]
    end
    for i in 2:h
        cum_risk[i, 1] = risk_map[i, 1] + cum_risk[i-1, 1]
    end
    for i in 2:h
        for j in 2:w
            cum_risk[i, j] = risk_map[i, j] + min(cum_risk[i-1, j], cum_risk[i, j-1])
        end
    end
    converged = false
    while !converged
        converged = true
        for j in 1:w
            for i in 1:h
                if j < w && cum_risk[i, j+1] < cum_risk[i, j] - risk_map[i, j]
                    converged = false
                    cum_risk[i, j] = cum_risk[i, j+1] + risk_map[i, j]
                end
                if j > 1 && cum_risk[i, j-1] < cum_risk[i, j] - risk_map[i, j]
                    converged = false
                    cum_risk[i, j] = cum_risk[i, j-1] + risk_map[i, j]
                end
            end
        end
        for i in 1:h
            for j in 1:w
                if i < h && cum_risk[i+1, j] < cum_risk[i, j] - risk_map[i, j]
                    converged = false
                    cum_risk[i, j] = cum_risk[i+1, j] + risk_map[i, j]
                end
                if i > 1 && cum_risk[i-1, j] < cum_risk[i, j] - risk_map[i, j]
                    converged = false
                    cum_risk[i, j] = cum_risk[i-1, j] + risk_map[i, j]
                end
            end
        end
    end
    return cum_risk[end]
end

function part2_nowiggling_rm(list, r)
    risk_map = parse.(Int, permutedims(reduce(hcat, collect.(list))))
    h, w = size(risk_map)
    risk_map = repeat(risk_map, r, r)
    for k in 2:r
        risk_map[1:h, (k-1)*w+1:k*w] .= mod1.(@view(risk_map[1:h, (k-2)*w+1:(k-1)*w]) .+ 1, 9)
    end
    for k in 2:r
        risk_map[(k-1)*h+1:k*h, :] .= mod1.(@view(risk_map[(k-2)*h+1:(k-1)*h, :]) .+ 1, 9)
    end
    mem = Dict{NTuple{4, Int}, Int}()
    h, w = size(risk_map)
    function risk!(mem, risk_map, n, i, m, j)
        (n, i, m, j) ∈ keys(mem) && return mem[(n, i, m, j)]
        if m == n + 1
            a = sum(@view(risk_map[m, i:j]))
            b = sum(@view(risk_map[n, i+1:j])) + risk_map[m, j]
            risk_value = min(a, b)
        else
            risk_value = typemax(Int)
            for k in i:j
                r1 = risk!(mem, risk_map, n, i, m-1, k)
                r2 = risk!(mem, risk_map, m-1, k, m, j)
                risk_value = min(risk_value, r1 + r2)
            end
        end
        push!(mem, (n, i, m, j) => risk_value)
        return risk_value
    end
    return risk!(mem, risk_map, 1, 1, h, w)
end

function part2_nowiggling(list, r)
    risk_map = parse.(Int, permutedims(reduce(hcat, collect.(list))))
    h, w = size(risk_map)
    risk_map = repeat(risk_map, r, r)
    for k in 2:r
        risk_map[1:h, (k-1)*w+1:k*w] .= mod1.(@view(risk_map[1:h, (k-2)*w+1:(k-1)*w]) .+ 1, 9)
    end
    for k in 2:r
        risk_map[(k-1)*h+1:k*h, :] .= mod1.(@view(risk_map[(k-2)*h+1:(k-1)*h, :]) .+ 1, 9)
    end
    risk_map[1, 1] = 0
    h, w = size(risk_map)
    for j in 3:w
        risk_map[1, j] += risk_map[1, j-1]
    end
    for i in 3:h
        risk_map[i, 1] += risk_map[i-1, 1]
    end
    for i in 2:h
        for j in 2:w
            risk_map[i, j] += min(risk_map[i-1, j], risk_map[i, j-1])
        end
    end
    return risk_map[end]
end

testit("Test Part 1", part2_nowiggling_rm(list_test_str, 1), ==(40))
testit("Challenge Part 1", part2_nowiggling_rm(list_str, 1), ==(447))
testit("Test Part 2", part2_nowiggling_rm(list_test_str, 5), ==(315))

testit("Test Part 1", part2_nowiggling(list_test_str, 1), ==(40))
testit("Challenge Part 1", part2_nowiggling(list_str, 1), ==(447))
testit("Test Part 2", part2_nowiggling(list_test_str, 5), ==(315))

testit("Test Part 1", part2(list_test_str, 1), ==(40))
testit("Challenge Part 1", part2(list_str, 1), ==(447))
testit("Test Part 2", part2(list_test_str, 5), ==(315))
testit("Challenge Part 2", part2(list_str, 5), ==(2825))

# got 2835 with two implementations of the nowiggling version (just to make sure).
# Had to write a second version allowing up/down and left/right paths, to get it right.


text:
using Test

@testset "All tests" begin
    include("types.jl")
    include("synthesis.jl")
    include("properties.jl")
    include("conversion.jl")
end


text:
# full analytic solution for
# global minimum variance portfolio
# returns the tip of the bullet

function globalmin(df::DataFrame)

  df     = copy(df)
  names  = colnames(df)
  onevec = ones(size(df, 2)-1)
  sigma  = cov(df)
  invsig = inv(sigma)
  port   = (invsig * onevec) / (onevec' * invsig * onevec)
  # port   = rowsums(DataArray(invsig)) / sum(invsig) # alternate method

  gminvar = DataFrame()

  for i in 2:length(names)
    colname  = names[i]
    within!(gminvar, :($colname = $port[$i-1]))
  end

  gminvar
end

function modelparams(df::AbstractDataFrame)
  df = copy(df)

  muhat  = mean(EachRow(df))
  #muhat  = float([mean(df[:,i]) for i in 1:size(df, 2)])
  varhat = float([var(df[:,i]) for i in 1:size(df, 2)])
  sdhat  = float([std(df[:,i]) for i in 1:size(df, 2)])
  covmat = cov(df)
  cormat = cor(df)

  return muhat, varhat, sdhat, covmat, cormat
end

function efficient(x)
  #code here
end

function tangency(x)
  #code here
end

function frontier()
  #code here
end


text:
using COSMO, Random, Test
rng = Random.MersenneTwister(12345)

include("./UnitTests/COSMOTestUtils.jl")

@testset "All Unit Tests" begin

  include("./UnitTests/simple.jl")
  include("./UnitTests/sets.jl")
  include("./UnitTests/constraints.jl")
  include("./UnitTests/model.jl")
  include("./UnitTests/qp-lasso.jl")
  include("./UnitTests/qp-box.jl")
  include("./UnitTests/socp-lasso.jl")
  include("./UnitTests/closestcorr.jl")
  include("./UnitTests/print.jl")
  include("./UnitTests/InfeasibilityTests/runTests.jl")
  include("./UnitTests/algebra.jl")
  include("./UnitTests/splitvector.jl")
  include("./UnitTests/interface.jl")
  include("./UnitTests/moi_wrapper.jl")
end
nothing


text:
int_rules_4_1_1_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 4.1.1.3*(g*tan)^p*(a+b*sin)^m =#
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~(p') / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => (1 / ~a) * Antiderivative(sec(~e + ~f * ~x) ^ 2 * (~g * tan(~e + ~f * ~x)) ^ ~p, ~x) - (1 / (~b * ~g)) * Antiderivative(sec(~e + ~f * ~x) * (~g * tan(~e + ~f * ~x)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && NeQ(~p, -1))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (1 / ~f) * Subst(Antiderivative(((~x) ^ ~p * (~a + ~x) ^ (~m - (~p + 1) / 2)) / (~a - ~x) ^ ((~p + 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ((~p + 1) / 2))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (~a) ^ ~p * Antiderivative(sin(~e + ~f * ~x) ^ ~p / (~a - ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(~m, ~p) && EqQ(~p, 2 * ~m)))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~a) ^ ~p * Antiderivative(ExpandIntegrand((sin(~e + ~f * ~x) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m - ~p / 2)) / (~a - ~b * sin(~e + ~f * ~x)) ^ (~p / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(~m, ~p / 2) && (LtQ(~p, 0) || GtQ(~m - ~p / 2, 0))))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => Antiderivative(ExpandIntegrand((~g * tan(~e + ~f * ~x)) ^ ~p, (~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~a) ^ (2 * ~m) * Antiderivative(ExpandIntegrand((~g * tan(~e + ~f * ~x)) ^ ~p * sec(~e + ~f * ~x) ^ -(~m), (~a * sec(~e + ~f * ~x) - ~b * tan(~e + ~f * ~x)) ^ -(~m), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && ILtQ(~m, 0))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ 2 * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (~b * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~a * ~f * (2 * ~m - 1) * cos(~e + ~f * ~x)) - (1 / ((~a) ^ 2 * (2 * ~m - 1))) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~a * ~m - ~b * (2 * ~m - 1) * sin(~e + ~f * ~x))) / cos(~e + ~f * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && LtQ(~m, 0)))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ 2 * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => -((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * ~m * cos(~e + ~f * ~x)) + (1 / (~b * ~m)) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~b * (~m + 1) + ~a * sin(~e + ~f * ~x))) / cos(~e + ~f * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && Not(LtQ(~m, 0))))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ 4 * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) - Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 - 2 * sin(~e + ~f * ~x) ^ 2)) / cos(~e + ~f * ~x) ^ 4, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ(~m - 1 / 2))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 2, ~x) => -((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~a * ~f * tan(~e + ~f * ~x)) + (1 / (~b) ^ 2) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) * (~b * ~m - ~a * (~m + 1) * sin(~e + ~f * ~x))) / sin(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m - 1 / 2) && LtQ(~m, -1)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m') / tan(~(e') + ~(f') * ~x) ^ 2, ~x) => -((~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (~f * tan(~e + ~f * ~x)) + (1 / ~a) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (~b * ~m - ~a * (~m + 1) * sin(~e + ~f * ~x))) / sin(~e + ~f * ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m - 1 / 2) && Not(LtQ(~m, -1))))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 4, ~x) => (-2 / (~a * ~b)) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 2) / sin(~e + ~f * ~x) ^ 3, ~x) + (1 / (~a) ^ 2) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 2) * (1 + sin(~e + ~f * ~x) ^ 2)) / sin(~e + ~f * ~x) ^ 4, ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m - 1 / 2) && LtQ(~m, -1)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 4, ~x) => Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) + Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 - 2 * sin(~e + ~f * ~x) ^ 2)) / sin(~e + ~f * ~x) ^ 4, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegerQ(~m - 1 / 2) && Not(LtQ(~m, -1))))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => ((sqrt(~a + ~b * sin(~e + ~f * ~x)) * sqrt(~a - ~b * sin(~e + ~f * ~x))) / (~b * ~f * cos(~e + ~f * ~x))) * Subst(Antiderivative(((~x) ^ ~p * (~a + ~x) ^ (~m - (~p + 1) / 2)) / (~a - ~x) ^ ((~p + 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && IntegerQ(~p / 2)))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => (((~g * tan(~e + ~f * ~x)) ^ (~p + 1) * (~a - ~b * sin(~e + ~f * ~x)) ^ ((~p + 1) / 2) * (~a + ~b * sin(~e + ~f * ~x)) ^ ((~p + 1) / 2)) / (~f * ~g * (~b * sin(~e + ~f * ~x)) ^ (~p + 1))) * Subst(Antiderivative(((~x) ^ ~p * (~a + ~x) ^ (~m - (~p + 1) / 2)) / (~a - ~x) ^ ((~p + 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && (EqQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(IntegerQ(~m)) && Not(IntegerQ(~p))))
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (1 / ~f) * Subst(Antiderivative(((~x) ^ ~p * (~a + ~x) ^ ~m) / ((~b) ^ 2 - (~x) ^ 2) ^ ((~p + 1) / 2), ~x), ~x, ~b * sin(~e + ~f * ~x)) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegerQ((~p + 1) / 2))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => Antiderivative(ExpandIntegrand((~g * tan(~e + ~f * ~x)) ^ ~p, (~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~p], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 2, ~x) => Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m * (1 - sin(~e + ~f * ~x) ^ 2)) / sin(~e + ~f * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 4, ~x) => ((-(cos(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (3 * ~a * ~f * sin(~e + ~f * ~x) ^ 3) - ((3 * (~a) ^ 2 + (~b) ^ 2 * (~m - 2)) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (3 * (~a) ^ 2 * ~b * ~f * (~m + 1) * sin(~e + ~f * ~x) ^ 2)) - (1 / (3 * (~a) ^ 2 * ~b * (~m + 1))) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1) / sin(~e + ~f * ~x) ^ 3) * Simp(((6 * (~a) ^ 2 - (~b) ^ 2 * (~m - 1) * (~m - 2)) + ~a * ~b * (~m + 1) * sin(~e + ~f * ~x)) - (3 * (~a) ^ 2 - (~b) ^ 2 * ~m * (~m - 2)) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (LtQ(~m, -1) && IntegerQ(2 * ~m)))
    #= Antiderivative(((~a)+(~b')*sin((~e')+(~f')*(~x)))^(~m)/tan((~e')+(~f')*(~x))^4,~x) :=  -cos(e+f*x)*(a+b*sin(e+f*x))^(m+1)/(3*a*f*sin(e+f*x)^3) - cos(e+f*x)*(a+b*sin(e+f*x))^(m+1)/(b*f*m*sin(e+f*x)^2) - 1/(3*a*b*m)*Antiderivative((a+b*sin(e+f*x))^m/sin(e+f*x)^3* Simp(6*a^2-b^2*m*(m-2)+a*b*(m+3)*sin(e+f*x)-(3*a^2-b^2*m*(m-1))* sin(e+f*x)^2,x),x) <-- FreeQ([a,b,e,f,m],x) && NeQ(a^2-b^2,0) &&  Not(LtQ(m,-1)) && IntegerQ(2*m) =#
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 4, ~x) => ((-(cos(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (3 * ~a * ~f * sin(~e + ~f * ~x) ^ 3) - (~b * (~m - 2) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (6 * (~a) ^ 2 * ~f * sin(~e + ~f * ~x) ^ 2)) - (1 / (6 * (~a) ^ 2)) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / sin(~e + ~f * ~x) ^ 2) * Simp(((8 * (~a) ^ 2 - (~b) ^ 2 * (~m - 1) * (~m - 2)) + ~a * ~b * ~m * sin(~e + ~f * ~x)) - (6 * (~a) ^ 2 - (~b) ^ 2 * ~m * (~m - 2)) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (Not(LtQ(~m, -1)) && IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m / tan(~(e') + ~(f') * ~x) ^ 6, ~x) => ((-(cos(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (5 * ~a * ~f * sin(~e + ~f * ~x) ^ 5) - (~b * (~m - 4) * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (20 * (~a) ^ 2 * ~f * sin(~e + ~f * ~x) ^ 4)) + (~a * cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / ((~b) ^ 2 * ~f * ~m * (~m - 1) * sin(~e + ~f * ~x) ^ 3) + (cos(~e + ~f * ~x) * (~a + ~b * sin(~e + ~f * ~x)) ^ (~m + 1)) / (~b * ~f * ~m * sin(~e + ~f * ~x) ^ 2) + (1 / (20 * (~a) ^ 2 * (~b) ^ 2 * ~m * (~m - 1))) * Antiderivative(((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / sin(~e + ~f * ~x) ^ 4) * Simp(((60 * (~a) ^ 4 - 44 * (~a) ^ 2 * (~b) ^ 2 * (~m - 1) * ~m) + (~b) ^ 4 * ~m * (~m - 1) * (~m - 3) * (~m - 4) + ~a * ~b * ~m * (20 * (~a) ^ 2 - (~b) ^ 2 * ~m * (~m - 1)) * sin(~e + ~f * ~x)) - ((40 * (~a) ^ 4 + (~b) ^ 4 * ~m * (~m - 1) * (~m - 2) * (~m - 4)) - 20 * (~a) ^ 2 * (~b) ^ 2 * (~m - 1) * (2 * ~m + 1)) * sin(~e + ~f * ~x) ^ 2, ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~m], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (NeQ(~m, 1) && IntegerQ(2 * ~m)))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~p / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => ((~a / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ ~p / sin(~e + ~f * ~x) ^ 2, ~x) - ((~b * ~g) / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ (~p - 1) / cos(~e + ~f * ~x), ~x)) - (((~a) ^ 2 * (~g) ^ 2) / ((~a) ^ 2 - (~b) ^ 2)) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ (~p - 2) / (~a + ~b * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(2 * ~p) && GtQ(~p, 1)))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~p / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => ((1 / ~a) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ ~p / cos(~e + ~f * ~x) ^ 2, ~x) - (~b / ((~a) ^ 2 * ~g)) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ (~p + 1) / cos(~e + ~f * ~x), ~x)) - (((~a) ^ 2 - (~b) ^ 2) / ((~a) ^ 2 * (~g) ^ 2)) * Antiderivative((~g * tan(~e + ~f * ~x)) ^ (~p + 2) / (~a + ~b * sin(~e + ~f * ~x)), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && (IntegersQ(2 * ~p) && LtQ(~p, -1)))
    @apply_utils Antiderivative(sqrt(~(g') * tan(~(e') + ~(f') * ~x)) / (~a + ~(b') * sin(~(e') + ~(f') * ~x)), ~x) => ((sqrt(cos(~e + ~f * ~x)) * sqrt(~g * tan(~e + ~f * ~x))) / sqrt(sin(~e + ~f * ~x))) * Antiderivative(sqrt(sin(~e + ~f * ~x)) / (sqrt(cos(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(1 / (sqrt(~g * tan(~(e') + ~(f') * ~x)) * (~a + ~(b') * sin(~(e') + ~(f') * ~x))), ~x) => (sqrt(sin(~e + ~f * ~x)) / (sqrt(cos(~e + ~f * ~x)) * sqrt(~g * tan(~e + ~f * ~x)))) * Antiderivative(sqrt(cos(~e + ~f * ~x)) / (sqrt(sin(~e + ~f * ~x)) * (~a + ~b * sin(~e + ~f * ~x))), ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g], ~x) && NeQ((~a) ^ 2 - (~b) ^ 2, 0)
    @apply_utils Antiderivative(tan(~(e') + ~(f') * ~x) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~m, ~x) => Antiderivative(ExpandIntegrand((sin(~e + ~f * ~x) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m) / (1 - sin(~e + ~f * ~x) ^ 2) ^ (~p / 2), ~x), ~x) <-- FreeQ([~a, ~b, ~e, ~f], ~x) && (NeQ((~a) ^ 2 - (~b) ^ 2, 0) && IntegersQ(~m, ~p / 2))
    @apply_utils Antiderivative((~(g') * tan(~(e') + ~(f') * ~x)) ^ ~(p') * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => Unintegrable((~g * tan(~e + ~f * ~x)) ^ ~p * (~a + ~b * sin(~e + ~f * ~x)) ^ ~m, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x)
    @apply_utils Antiderivative((~(g') * cot(~(e') + ~(f') * ~x)) ^ ~p * (~a + ~(b') * sin(~(e') + ~(f') * ~x)) ^ ~(m'), ~x) => (~g) ^ (2 * IntPart(~p)) * (~g * cot(~e + ~f * ~x)) ^ FracPart(~p) * (~g * tan(~e + ~f * ~x)) ^ FracPart(~p) * Antiderivative((~a + ~b * sin(~e + ~f * ~x)) ^ ~m / (~g * tan(~e + ~f * ~x)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~e, ~f, ~g, ~m, ~p], ~x) && Not(IntegerQ(~p))
end

max_stars_repo_path: src/rules/4 Trig functions/4.1 Sine/4.1.1.3 (g tan)^p (a+b sin)^m.jl


text:
# ---
# title: 485. Max Consecutive Ones
# id: problem485
# author: Indigo
# date: 2021-02-03
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/max-consecutive-ones/description/>
# hidden: true
# ---
#
# Given a binary array, find the maximum number of consecutive 1s in this array.
#
# **Example 1:**
#
#
#
#     Input: [1,1,0,1,1,1]
#     Output: 3
#     Explanation: The first two digits or the last three digits are consecutive 1s.
#         The maximum number of consecutive 1s is 3.
#
#
# **Note:**
#
#   * The input array will only contain `0` and `1`.
#   * The length of input array is a positive integer and will not exceed 10,000
#
#
## @lc code=start
using LeetCode

function find_max_consecutive_ones(nums::Vector{Int})
    cnt, tmp = 0, 0
    for num in nums
        num == 1 ? tmp += 1 : tmp = 0
        cnt = max(tmp, cnt)
    end
    return max(tmp, cnt)
end
## @lc code=end


text:
# # [Example: Treiber stack](@id ex-treiberstack)

using Reagents

# [Treiber stack](https://en.wikipedia.org/wiki/Treiber_stack) is a simple
# concurrent data structure where an immutable list holds the data items

struct TSNode{T}
    head::T
    tail::Union{TSNode{T},Nothing}
end

const TSList{T} = Union{TSNode{T},Nothing}

# ... which is referenced through an atomically updatable memory location (here,
# a [`Reagents.Ref`](@ref)):

struct TreiberStack{T,Ref<:Reagents.Ref{TSList{T}}}
    head::Ref
end

Base.eltype(::Type{<:TreiberStack{T}}) where {T} = T

# An empty stack can be constructed as a reference to the empty list; i.e.,
# `nothing`:

TreiberStack{T}() where {T} = TreiberStack(Reagents.Ref{TSList{T}}(nothing))

# To push an element to the stack, we can use [`Reagents.Update`](@ref)

pushing(stack::TreiberStack) =
    Reagents.Update((xs, x) -> (TSNode(x, xs), nothing), stack.head)

# Let's see how it works.

using Test

function test_pushing()
    #=
    When a stack is created, its head points to `nothing`:
    =#
    stack = TreiberStack{Int}()
    @test stack.head[] === nothing  # empty

    #=
    We can create a reagent for pushing a value to the stack:
    =#
    reagent = pushing(stack)

    #=
    However, note that reagent does nothing when it's created.
    =#
    @test stack.head[] === nothing  # empty

    #=
    The reagent must be executed ("react") to invoke its side-effect:
    =#
    reagent(111)
    @test stack.head[] === TSNode(111, nothing)
end

# Similarly, we can pop off an element from the stack, again by
# [`Reagents.Update`](@ref). To support empty stack, we return `nothing` when
# it's empty and return `Some(value)` when we find a `value`:

trypopping(stack::TreiberStack) =
    Reagents.Update(stack.head) do xs, _ignored
        if xs === nothing
            return (nothing, nothing)
        else
            return (xs.tail, Some(xs.head))
        end
    end

# Here's how it works.

function test_trypopping()
    #=
    Let's push 111 and then 222 using the `pushing` reagent:
    =#
    stack = TreiberStack{Int}()
    pushing(stack)(111)
    pushing(stack)(222)

    #=
    We can invoke the reagent `trypopping(stack)` by calling it. Since
    `trypopping` ignores the input (see the argument `_ignored` in `trypopping`
    definition above), we can pass an arbitrary value to the reagent, e.g.,
    `nothing`.
    =#
    @test trypopping(stack)(nothing) === Some(222)
    #=
    For convinience, `nothing` is the default argument when the reagent is
    called without an argument:
    =#
    @test trypopping(stack)() === Some(111)
    #=
    Now that all values are popped, invoking `trypopping` returns `nothing`:
    =#
    @test trypopping(stack)() === nothing
end

# It is simple to wrap these reagents into the `Base` API:

Base.push!(stack::TreiberStack, value) = pushing(stack)(convert(eltype(stack), value))
Base.pop!(stack::TreiberStack) = something(trypopping(stack)())

# Note that this version of `pop!(stack::TreiberStack)` throws when the `stack`
# is empty.  See [Blocking containers](@ref ex-blockingcontainers) for a
# *generic* derivation of a blocking version of `pop!` that waits for the value
# to be `push!`ed.
#
# For more usage examples, see
# [`/test/ReagentsTests/src/test_treiberstack.jl`](https://github.com/tkf/Reagents.jl/blob/master/test/ReagentsTests/src/test_treiberstack.jl).


text:
@testset "Automatic differentiation" begin
    @testset "Real -> Real" begin
        for T in (Float32, Float64, Int)
            x = rand(T)
            @test (@inferred gradient(x -> 2x^2 - 3x + 2, x))::T == 4x - 3
            @test (@inferred gradient(x -> 2x^2 - 3x + 2, x, :all))::Tuple{T, T} == (4x - 3, 2x^2 - 3x + 2)
            @test (@inferred gradient(x -> 3, x))::Int == zero(T)
            @test (@inferred gradient(x -> 3, x, :all))::Tuple{Int, Int} == (zero(T), 3)
        end
    end
    @testset "Real -> Tensor" begin
        for RetType in (Tensor{Tuple{2,3}}, Tensor{Tuple{3, @Symmetry{3,3}}})
            for T in (Float32, Float64)
                @eval begin
                    x = rand($T)
                    @test (@inferred gradient(x -> $RetType((ij...) -> prod(ij) * x^2), x))::$RetType{$T} == $RetType((ij...) -> prod(ij) * 2x)
                    @test (@inferred gradient(x -> $RetType((ij...) -> prod(ij) * x^2), x, :all))::Tuple{$RetType{$T}, $RetType{$T}} == ($RetType((ij...) -> prod(ij) * 2x), $RetType((ij...) -> prod(ij) * x^2))
                    @test (@inferred gradient(x -> $RetType((ij...) -> prod(ij)), x))::$RetType{Int} == zero($RetType)
                    @test (@inferred gradient(x -> $RetType((ij...) -> prod(ij)), x, :all))::Tuple{$RetType{Int}, $RetType{Int}} == (zero($RetType), $RetType((ij...) -> prod(ij)))
                end
            end
        end
    end
    @testset "Tensor -> Real" begin
        I₁ = x -> tr(x)
        I₂ = x -> (tr(x)^2 - tr(x^2)) / 2
        I₃ = x -> det(x)
        for RetType in (SecondOrderTensor, SymmetricSecondOrderTensor)
            for T in (Float32, Float64), dim in 1:4
                Random.seed!(1234)
                x = rand(RetType{dim, T})
                @test (@inferred gradient(I₁, x))::RetType{dim, T} ≈ one(x)
                @test (@inferred gradient(I₂, x))::RetType{dim, T} ≈ I₁(x)*one(x) - x'
                @test (@inferred gradient(I₃, x))::RetType{dim, T} ≈ det(x)*inv(x)'
                @test (@inferred gradient(I₁, x, :all))[1]::RetType{dim, T} ≈ one(x)
                @test (@inferred gradient(I₂, x, :all))[1]::RetType{dim, T} ≈ I₁(x)*one(x) - x'
                @test (@inferred gradient(I₃, x, :all))[1]::RetType{dim, T} ≈ det(x)*inv(x)'
                @test (@inferred gradient(I₁, x, :all))[2]::T ≈ I₁(x)
                @test (@inferred gradient(I₂, x, :all))[2]::T ≈ I₂(x)
                @test (@inferred gradient(I₃, x, :all))[2]::T ≈ I₃(x)
                @test (@inferred gradient(x -> 1, x))::RetType{dim, Int} == zero(x)
                @test (@inferred gradient(x -> 1, x, :all))::Tuple{RetType{dim, Int}, Int} == (zero(x), 1)
            end
        end
    end
    @testset "Tensor -> Tensor" begin
        for T in (Float32, Float64), dim in 1:4
            x = rand(SecondOrderTensor{dim, T})
            y = rand(SymmetricSecondOrderTensor{dim, T})
            @test (@inferred gradient(identity, x))::FourthOrderTensor{dim, T} == one(FourthOrderTensor{dim, T})
            @test (@inferred gradient(identity, y))::SymmetricFourthOrderTensor{dim, T} == one(SymmetricFourthOrderTensor{dim, T})
            @eval begin
                @test (@inferred gradient(x -> one(SecondOrderTensor{$dim, Int}), $x))::FourthOrderTensor{$dim, Int} == zero(FourthOrderTensor{$dim, Int})
                @test (@inferred gradient(x -> one(SecondOrderTensor{$dim, Int}), $y))::Tensor{Tuple{$dim,$dim,@Symmetry{$dim,$dim}}, Int} == zero(FourthOrderTensor{$dim, Int})
            end
        end
    end
end


text:
using Compat.Test
import Roots.newton, Roots.halley

@test abs(newton(sin, cos, 0.5) - 0.0) <= 100*eps(1.0)
@test newton(cos, x -> -sin(x), 1.0)  ≈ pi/2
@test newton(x -> x^2 - 2x - 1, x -> 2x - 2, 3.0)  ≈ 2.414213562373095
@test abs(newton(x -> exp(x) - cos(x), x -> exp(x) + sin(x), 3.0) - 0.0) <= 1e-14
@test halley(x -> x^2 - 2x - 1,x -> 2x - 2,x -> 2, 3.0)  ≈ 2.414213562373095
a = halley(x -> exp(x) - cos(x),
           x -> exp(x) + sin(x),
           x -> exp(x) + cos(x), 3.0)
@test abs(a - 0.0) <= 1e-14

## test with Complex input

@test real(newton(x ->  x^3 - 1, x ->  3x^2, 1+im)) ≈ 1.0
@test real(newton(x ->  x^3 - 1, x ->  3x^2, 1+10im)) ≈ (-1/2)


text:
"""
    simulate_detector(sim_config_filename)

AbstractString -> <detector simulation>

Simulate detector according to given simulation configuration.
What type <detector simulation> is depends on the given method of simulation.
"""
function simulate_detector(sim_config::PropDict)
    det_meta = propdict(sim_config.detector_metadata)
    env = Environment(sim_config)
    ps_simulator = PSSimulator(sim_config)

    simulate_detector(det_meta, env, sim_config.pss.cached_name, ps_simulator)
end

function simulate_detector(sim_config_filename::AbstractString)
    simulate_detector(propdict(sim_config_filename))
end

function simulate_detector(det_metadata::AbstractString, sim_config_filename::AbstractString)
    # merge detector info with the rest
    sim_config = load_config(det_metadata, sim_config_filename)
    simulate_detector(sim_config)
end

"""
    simulate_detector(det_meta, env, config_name, ::SSDSimulator)

PropDict, Environment, AbstractString -> SSD.Simulation

Look up cached SSD simulation .h5f file, and if does not exist,
    simulate the detector according to geometry as given in LEGEND metadata <det_meta>
    and environment settings given in <env> using SolidStateDetectors.
The simulation will be cached based on <config_name>.

"""
function simulate_detector(det_meta::PropDict, env::Environment, cached_name::AbstractString, ssd_sim::SSDSimulator)
    det_h5 = joinpath("cache", cached_name*"_ssd.h5f")
    if isfile(det_h5)
        @info("Reading SSD simulation from cached file $det_h5")
        simulation = SolidStateDetectors.ssd_read(det_h5, Simulation)
    else
        @info("Simulating $cached_name with SSD from scratch for given settings")
        ssd_conf = ssd_config(det_meta, env, ssd_sim)
        simulation = simulate_ssd(ssd_conf)

        if !ispath(dirname(det_h5)) mkpath(dirname(det_h5)) end
        SolidStateDetectors.ssd_write(det_h5, simulation)
        @info("-> Saved cached simulation to $det_h5")
    end

    simulation
end

"""
    simulate_detector(det_meta, env, config_name, ps_simulator)

PropDict, Environment, AbstractString, Siggen -> SigGenSetup

Construct a fieldgen/siggen configuration file
    according to geometry as given in LEGEND metadata <det_meta>
    and environment settings given in <env>.
Look up fieldgen generated electric field and weighting potential files
    based on given <config_name>, and if not found, call fieldgen.
"""
function simulate_detector(det_meta::PropDict, env::Environment, cached_name::AbstractString, ps_simulator::SiggenSimulator)
    @info "Constructing Fieldgen/Siggen configuration file"
    # returns the name of the resulting siggen config file
    # and the name of (already or to be) generated weighting potential file
    siggen_config_name, fieldgen_wp_name = siggen_config(det_meta, env, ps_simulator, cached_name)
    fieldgen_wp_name = joinpath("cache", fieldgen_wp_name)

    # if the WP file with such a name exists...
    if isfile(fieldgen_wp_name)
        #...do nothing, siggen will later read the files based on the generated conifg
        @info "Reading cached fields from $fieldgen_wp_name"
    else
        #...otherwise call fieldgen
        @info "_|~|_|~|_|~|_ Fieldgen simulation"
        fieldgen(siggen_config_name)
        @info "_|~|_|~|_|~|_ Fieldgen simulation complete"
    end

    # a SigGenSetup object
    SigGenSetup(siggen_config_name)
end

# ------------------------------------------------------------------- SSD

"""
    simulate_ssd(ssd_config)

Dict -> SSD.Simulation

Simulate detector based on given SSD config
"""
function simulate_ssd(ssd_config::Dict)

    @info "_|~|_|~|_|~|_ SSD simulation"

    simulation = Simulation(SolidStateDetector{Float32}(ssd_config))

    println("...electric potential")
    calculate_electric_potential!( simulation,
                               max_refinements = 4)

    println("...electric field")
    calculate_electric_field!(simulation, n_points_in_φ = 72)

    println("...drift field")
    calculate_drift_fields!(simulation)

    println("...weighting potential")
    for contact in simulation.detector.contacts
        calculate_weighting_potential!(simulation, contact.id, max_refinements = 4, n_points_in_φ = 2, verbose = false)
    end

    @info "_|~|_|~|_|~|_ SSD detector simulation complete"

    simulation
end

"""
    ssd_config(meta, env)

PropDict, Environment -> Dict

Construct SSD type config based on
    geometry as given in LEGEND metadata <meta>
    and environment settings given in <env>.
"""
function ssd_config(meta::PropDict, env::Environment, ssd_sim::SSDSimulator)
    ## some parameters that will be used multiple times later
    # or need a calculation
    cylinder_height = meta.geometry.height_in_mm
    cylinder_radius = meta.geometry.radius_in_mm
    borehole_height = cylinder_height - meta.geometry.well.gap_in_mm
    borehole_r_bottom = meta.geometry.well.radius_in_mm
    borehole_r_top = borehole_r_bottom + borehole_height * tan(meta.geometry.taper.top.inner.angle_in_deg / 180. * π)
    cone_r_top = cylinder_radius - meta.geometry.taper.top.outer.height_in_mm * tan(meta.geometry.taper.top.outer.angle_in_deg / 180. * π)
    cone_height = meta.geometry.taper.top.outer.height_in_mm

    ## grid
    # kind of initialize so that it exists out of the loop?
    grid = Dict()

    if ssd_sim.coord == "cylindrical"
        grid = Dict(
            "axes" => Dict(
                "r" => Dict("to" => cylinder_radius*1.2, "boundaries" => "inf"), # leave some margin
                "phi" => Dict("from" => 0, "to" => ssd_sim.comp == "3D" ? 360 : 0, "boundaries" => "periodic"),
                "z" => Dict(
                    "from" => -0.2*cylinder_height,
                    "to" => cylinder_height*1.2, # leave some margin
                    "boundaries" => Dict("left" => "inf", "right" => "inf")
                )
            )
        )
    else # cartesian
        grid = Dict(
            "axes" => Dict(
                "x" => Dict("from" => -cylinder_radius*1.2, "to" => cylinder_radius*1.2, "boundaries" => "inf"),
                "y" => Dict("from" => -cylinder_radius*1.2, "to" => cylinder_radius*1.2, "boundaries" => "inf"),
                "z" => Dict("from" => -0.2*cylinder_height, "to" => cylinder_height*1.2, "boundaries" => "inf")
            )
        )
    end

    grid = merge(grid, Dict("coordinates" => ssd_sim.coord))

    ## charge density model
    # kind of initialize so that it exists out of the loop?
    charge_density_model = Dict()

    if ssd_sim.coord == "cylindrical"
        charge_density_model = Dict(
            "r" => Dict("init" => 0, "gradient" => 0),
            "phi" => Dict("init" => 0.0, "gradient" => 0.0),
            "z" => Dict("init" => -1e7, "gradient" => 1e5)
        )
    else # cartesian
        charge_density_model = Dict(
            "x" => Dict("init" => 0.0, "gradient" => 0.000),
            "y" => Dict("init" => 0.0, "gradient" => 0.0),
            "z" => Dict("init" => 1e7, "gradient" => 5e4)
        )
    end

    # for now we only use linear charge density model for either cartesian or cylindrical coordinates
    # charge_density_model["name"] = "linear"
    charge_density_model = merge(charge_density_model, Dict("name" => "linear"))

    ## geometry
    dct = Dict(
        "name" => meta.det_name,
        "medium" => env.medium,
        "units" => Dict(
            "length" => "mm",
            "angle" => "deg",
            "potential" => "V",
            "temperature" => "K"
        ),
        "grid" => grid,
        "objects" => [
            # crystal
            Dict(
                "type" => "semiconductor",
                "material" => "HPGe",
                "bulk_type" => "p",
                "temperature" => env.crystal_t,
                "charge_density_model" => charge_density_model,
                "geometry" => Dict(
                    "type" => "difference",
                    "parts" => [
                        # crystal
                        Dict(
                            "type" => "difference",
                            "parts" => [
                                # cylinder
                                Dict(
                                    "name" => "Initial Cylinder",
                                    "type" => "tube",
                                    "r" => Dict(
                                        "from" => 0.0,
                                        "to" => cylinder_radius,
                                    ),
                                    "phi" => Dict("from" => 0.0, "to" => 360.0),
                                    "h" => cylinder_height,
                                    "translate" => Dict("z" => 0.0)
                                )#, # cylinder
                                # upper cone
                                Dict(
                                    "name" => "Upper cone",
                                    "type" => "cone",
                                    "r" => Dict(
                                        "bottom" => Dict("from" => cylinder_radius, "to" => cylinder_radius+1),
                                        "top" => Dict("from" => cone_r_top, "to" => cylinder_radius+1)
                                    ),
                                    "phi" => Dict("from" => 0.0, "to" => 360.0),
                                    "h" => meta.geometry.taper.top.outer.height_in_mm,
                                    "translate" => Dict("z" => cylinder_height - cone_height)
                                ) # upper cone
                            ] # crystal parts
                        ), # crystal
                        # borehole
                        Dict(
                            "name" => "borehole",
                            "type" => "cone",
                            "r" => Dict(
                                "bottom" => Dict(
                                    "from" => 0.0,
                                    "to" => borehole_r_bottom
                                ),
                                "top" => Dict(
                                    "from" => 0.0,
                                    "to" => borehole_r_top
                                )
                            ),
                            "phi" => Dict("from" => 0.0, "to" => 360.0),
                            "h" => borehole_height,
                            "translate" => Dict("z" => meta.geometry.well.gap_in_mm)
                        ), # borehole
                        # groove
                        Dict(
                            "name" => "ditch",
                            "type" => "tube",
                            "r" => Dict(
                                "from" => meta.geometry.groove.outer_radius_in_mm - meta.geometry.groove.width_in_mm,
                                "to" => meta.geometry.groove.outer_radius_in_mm
                            ),
                            "phi" => Dict("from" => 0.0, "to" => 360.0),
                            "h" => meta.geometry.groove.depth_in_mm
                        ) # groove
                    ] # parts
                ) # geometry
            ), # crystal
            # p+ contact
            Dict(
                "type" => "contact",
                "material" => "HPGe",
                "channel" => 1,
                "potential" => 0.0,
                "geometry" => Dict(
                    "type" => "tube",
                    "r" => Dict(
                        "from" => 0.0,
                        "to" => meta.geometry.contact.radius_in_mm
                    ),
                    "phi" => Dict("from" => 0.0, "to" => 360.0),
                    "h" => meta.geometry.contact.depth_in_mm
                )
            ), # p+ contact
            # n+ contact
            Dict(
                "type" => "contact",
                "material" => "HPGe",
                "channel" => 2,
                "potential" => env.op_voltage == 0 ? meta.characterization.manufacturer.op_voltage_in_V : env.op_voltage,
                "geometry" => Dict(
                    "type" => "union",
                    "parts" => [Dict() for i in 1:6] # to be filled later
                )
            ) # n+ contact
        ] # objects
    )

    ## Construct n+ surface

    # bottom lid
    dct["objects"][3]["geometry"]["parts"][1] = Dict(
        "name" => "bottom lid",
        "type" => "tube",
        "r" => Dict(
            "from" => meta.geometry.groove.outer_radius_in_mm,
            "to" => cylinder_radius # cylinder radius
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => 0
    )

    # side lining
    dct["objects"][3]["geometry"]["parts"][2] = Dict(
        "name" => "side lining",
        "type" => "tube",
        "r" => Dict(
            "from" => cylinder_radius,
            "to" => cylinder_radius
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => cylinder_height - cone_height
    )

    # upper cone
    dct["objects"][3]["geometry"]["parts"][3] = Dict(
        "name" => "Upper cone",
        "type"=> "cone",
        "r" => Dict(
            "bottom" => Dict("from" => cylinder_radius, "to" => cylinder_radius),
            "top" => Dict("from" => cone_r_top, "to" => cone_r_top)
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => cone_height,
        "translate" => Dict("z" => cylinder_height - cone_height)
    )

    # top lid
    dct["objects"][3]["geometry"]["parts"][4] = Dict(
        "name" => "top lid",
        "type" => "tube",
        "r" => Dict(
            "from" => borehole_r_top,
            "to" => cone_r_top
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => 0.0,
        "translate" => Dict("z" => cylinder_height)
    )

    # borehole lining
    dct["objects"][3]["geometry"]["parts"][5] = Dict(
        "name" => "borehole linning",
        "type" => "cone",
        "r" => Dict(
            "bottom" => Dict(
                "from" => borehole_r_bottom,
                "to" => borehole_r_bottom
            ),
            "top" => Dict(
                "from" => borehole_r_top,
                "to" => borehole_r_top
            )
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => borehole_height,
        "translate" => Dict("z" => meta.geometry.well.gap_in_mm)
    )

    # borehole bottom lid
    dct["objects"][3]["geometry"]["parts"][6] = Dict(
        "name" => "borehole bottom lid",
        "type" => "tube",
        "r" => Dict(
            "from" => 0.0,
            "to" => borehole_r_bottom
        ),
        "phi" => Dict("from" => 0.0, "to" => 360.0),
        "h" => 0.0,
        "translate" => Dict("z" => meta.geometry.well.gap_in_mm)
    )

    dct
end

"""
    ssd_config(det_meta)

AbstractString -> Dict

Construct SSD type config based on geometry as given in LEGEND metadata <det_meta>.
Set arbitrary envorinment variables and simulation settings (used only for geometry visualization)
"""
function ssd_config(det_meta::AbstractString)
    ssd_config(propdict(det_meta), Environment(), SSDSimulator())
end

"""
    siggen_config(meta, env, siggen_sim, config_name)

PropDict, Environment, SiggenSimulator, String -> String, String

Construct Siggen/Fieldgen config based on
    geometry as given in LEGEND metadata <meta>,
    environment settings given in <env>,
    and fieldgen settings contained in <siggen_sim>.

The resulting siggen config will be cached based on given <config_name>
The function returns the path to the cached siggen config file, and
    the relative path to the fieldgen WP file (to check if already exists later)
"""
function siggen_config(meta::PropDict, env::Environment, siggen_sim::SiggenSimulator, cached_name::AbstractString)
    # function pss_config(meta::PropDict, env::Environment, siggen_sim::SiggenSimulatorss, config_name::AbstractString)
    # construct siggen geometry based on LEGEND metadata JSON + extra environment settings
    println("...detector geometry")
    siggen_geom = meta2siggen(meta, env)

    # construct lines for future siggen config file
    detector_lines = Vector{String}()
    push!(detector_lines, "# detector related parameters")
    for (param, value) in siggen_geom
        push!(detector_lines, "$param   $value")
    end

    println("...fieldgen configuration")
    # create filenames for future/existing fieldgen output
    fieldgen_wp_name = joinpath("fieldgen", cached_name*"_fieldgen_WP.dat")
    fieldgen_names = Dict(
        "drift_name" => joinpath("..", siggen_sim.drift_vel),
        # "drift_name" => joinpath("..", "data", "drift_vel_tcorr.tab"),
        "field_name" => joinpath("fieldgen", cached_name*"_fieldgen_EV.dat"),
        "wp_name" => fieldgen_wp_name
    )

    # add corresponding lines to existing detector geometry lines
    push!(detector_lines, "")
    push!(detector_lines, "# fieldgen filenames")
    for (param, value) in fieldgen_names
        push!(detector_lines, "$param   $value")
    end

    # read lines from user input for fieldgen
    fieldgen_lines = readlines(open(siggen_sim.fieldgen_config))
    fieldgen_lines = replace.(fieldgen_lines, "\t" => "   ")

    # unite constructed detector geometry and fieldgen settings
    total_lines = vcat(detector_lines, [""], fieldgen_lines)

    # write a siggen/fieldgen config file
    sig_conf_name = joinpath("cache", cached_name*"_siggen_config.txt")
    writedlm(sig_conf_name, total_lines)
    println("...fieldgen/siggen config written to $sig_conf_name")

    # return resulting fieldgen/siggen config name
    # (currently kinda redundant but that's because we have no idea what's gonna happen later)
    sig_conf_name, fieldgen_wp_name
end

"""
    meta2siggen(meta, env)

PropDict, Environment -> Dict

Construct siggen type config in a form of Dict based on
    geometry as given in LEGEND metadata <meta>
    and environment settings given in <env>.
"""
function meta2siggen(meta::PropDict, env::Environment)
    Dict(
    # crystal
        # z length: crystal height
        "xtal_length" => meta.geometry.height_in_mm,
        # radius: crystal radius
        "xtal_radius" => meta.geometry.radius_in_mm,
    # contact
        # point contact length
        "pc_length" => meta.geometry.contact.depth_in_mm,
        # point contact radius
        "pc_radius" => meta.geometry.contact.radius_in_mm,
    # groove
        # wrap-around radius. Set to zero for ORTEC: groove outer radius
        "wrap_around_radius" => meta.geometry.groove.outer_radius_in_mm,
        # depth of ditch next to wrap-around for BEGes. Set to zero for ORTEC: groove height
        "ditch_depth" => meta.geometry.groove.depth_in_mm,
        # width of ditch next to wrap-around for BEGes. Set to zero for ORTEC: groove inner radius
        "ditch_thickness" => meta.geometry.groove.width_in_mm,
    # borehole
        # length of hole, for inverted-coax style
        "hole_length" => meta.geometry.height_in_mm - meta.geometry.well.gap_in_mm,
        # radius of hole, for inverted-coax style
        "hole_radius" => meta.geometry.well.radius_in_mm,
    # tapering

    # comment from David Radford:

    # The bottom taper is always at 45 degrees, but the top outer taper is not.
    # It's width/angle is defined using either the outer_taper_width parameter or the taper_angle parameter.
    # Likewise, the inner taper width/angle is defined using either the inner_taper_width parameter or the taper_angle parameter.
    # The inner_taper_length can be anything from zero to the hole_length
    # TODO: check if bottom taper angle is 45 deg, if not issue a warning for siggen case

        # size of 45-degree taper at bottom of ORTEC-type crystal (for r=z)
        "bottom_taper_length" => meta.geometry.taper.bottom.outer.height_in_mm,
        # z-length of outside taper for inverted-coax style
        "outer_taper_length" => meta.geometry.taper.top.outer.height_in_mm,
        # z-length of inside (hole) taper for inverted-coax style
        "inner_taper_length" => meta.geometry.taper.top.inner.height_in_mm,
        # taper angle in degrees, for inner or outer taper
        "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg,
        # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg,
        # depth of full-charge-collection boundary for Li contact (not currently used)
        "Li_thickness" => meta.geometry.dl_thickness_in_mm,

    # configuration for mjd_fieldgen (calculates electric fields & weighing potentials)
        # detector bias for fieldgen, in Volts
        "xtal_HV" => env.op_voltage == 0 ? meta.characterization.manufacturer.op_voltage_in_V : env.op_voltage,

    # configuration for signal calculation
        # crystal temperature in Kelvin
        "xtal_temp" => env.crystal_t
    )
end


text:
@testset "k-Means" begin
    seed!(42)
    N = 20
    nDim = 3
    M = 10
    X = rand(nDim, N)
    x = ColVecs(X)
    alg = KmeansAlg(M)

    @test alg.metric == SqEuclidean()
    Z = inducingpoints(alg, x)
    @test repr(alg) == "k-Means Selection of Inducing Points (m : $(M))"
    @test length(Z) == M
    Z = inducingpoints(alg, X; obsdim=2, weights=rand(N))
    @test length(Z) == M
    @test length(first(Z)) == nDim
end


text:
# Compares two methods for simulating a homogeneous Poisson point process,
# where the idea behind one method is faster than the other. In this code,
# the two methods are labelled A and B.
#
# Method A simulates all the Poisson ensembles randomly by first randomly
# generating all the Poisson random variables in one step. It then randomly
# positions all the points (across all ensembles) in one step. All the
# points are then are then separated accordingly into ensembles.
#
# Method B iterates through a for-loop, and for each iteration, it randomly
# generates a Poisson variable and positions the points for each ensemble.
#
# Method A uses more vectorization than Method B so it should be faster
# than Method B in general.
#
# Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts
# For more details, see the post:
# hpaulkeeler.com/testing-the-julia-language-with-point-process-simulations/

clearconsole();

using Random;
using Distributions; #for random simulations
using Plots; #for plotting

Random.seed!(1);

###START Parameters START###
numbSim=10^5; #number of simulations

#Point process parameters
lambda=5; #intensity (ie mean density) of Poisson point process

#Simulation window parameters
xMin=0;xMax=1;
yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta; #area of rectangle
massTotal=areaTotal*lambda;  #total measure/mass of the point process
###END Parameters END###

###START Simulation section START###
###START Method A: Generate *all* ensembles at once START###
@time begin
    xxCellA=Array{Array{Float64}}(undef,numbSim);
    yyCellA=Array{Array{Float64}}(undef,numbSim);

    numbPointsA=rand(Poisson(massTotal),numbSim); #Poisson number of points
    numbPointsCumA=cumsum(numbPointsA);
    numbPointsTotal=numbPointsCumA[end];

    #uniform x/y coordinates of Poisson points
    xxAll=xDelta.*(rand(numbPointsTotal)).+xMin;#x coordinates of Poisson points
    yyAll=yDelta.*(rand(numbPointsTotal)).+yMin;#y coordinates of Poisson points

    #create some indexing for the array reshaping step
    indexFirst=copy(numbPointsCumA);
    indexFirst[2:end]=indexFirst[1:end-1].+1;
    indexFirst[1]=1;
    indexSecond=copy(numbPointsCumA);
    #reshape point arrays into an array of arrays.
    for ii in 1:numbSim
        xxCellA[ii]=(xxAll[indexFirst[ii]:indexSecond[ii]]);
        yyCellA[ii]=(yyAll[indexFirst[ii]:indexSecond[ii]]);
    end

    #yyCellAllA=[(yyAll[indexFirst[ii]:indexSecond[ii]]) for ii in 1:numbSim]
end;
###END Method A: Generate *all* ensembles at once END###

###START Method B: Generate each ensemble separately START###
@time begin
    xxCellB=Array{Array{Float64}}(undef,numbSim);
    yyCellB=Array{Array{Float64}}(undef,numbSim);
    numbPointsB=zeros(numbSim);
    #loop through for all ensembles
    for ss=1:numbSim
        numbPointsTemp=rand(Poisson(massTotal));#Poisson number of points
        xxCellB[ss]=xDelta.*(rand(numbPointsTemp)).+xMin;#x coordinates of Poisson points
        yyCellB[ss]=yDelta.*(rand(numbPointsTemp)).+yMin;#y coordinates of Poisson points
        numbPointsB[ss]=numbPointsTemp;
    end
end;
###END Method B: Generate each ensemble separately END###
###END Simulation section   END###

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
#=
An core that can serve as an array input.
=#
@with_kw mutable struct InputHandler{T, F <: AbstractFifo}
    # Clock period of the input.
    clockperiod :: Int = 1

    # Data reserve to write out.
    data :: Vector{T} = T[]
    pointer :: Int = 1

    # Output fifo.
    output :: F = F()
end

function update!(c::InputHandler)
    # If there is data left to write, try to write it to the output.
    if checkbounds(Bool, c.data, c.pointer) && iswriteready(c.output)
        write(c.output, c.data[c.pointer])
        c.pointer += 1

        if mod(c.pointer, 500) == 0
            println("On input $(c.pointer) of $(length(c.data))")
        end
    end
    writeupdate!(c.output)
end

@with_kw mutable struct OutputHandler{T, F <: AbstractFifo}
    clockperiod :: Int = 1
    data :: Vector{T} = T[]
    stoplength :: Int = 0
    input :: F = F()
end

function update!(c::OutputHandler)
    # If new data is available, read from the input.
    if isreadready(c.input)
        push!(c.data, read(c.input))

        if (0 < c.stoplength) && (c.stoplength <= length(c.data))
            throw(StopSimulation())
        end
    end
    # Update input.
    readupdate!(c.input)
end

clockperiod(c::Union{InputHandler,OutputHandler}) = c.clockperiod


text:
#!/usr/bin/env julia

using ConcurrentCalls

function xz_compress_serial(files)
    compress(input, output) = run(pipeline(input, `xz`, output))
    for file in files
        compress(file, DevNull)
    end
end

function xz_compress_concurrent(files)
    compress(input, output) = run(pipeline(input, `xz`, output))
    @cc for file in files
        compress(file, DevNull)
    end
end

files = readlines(STDIN)
if isempty(files)
    error("no input files")
end

println("--- statistics ---")
filesizes = map(filesize, files)
println("number of files:     $(length(files))")
println("mean of filesizes:   $(@sprintf("%.1f", mean(filesizes)))")
println("median of filesizes: $(@sprintf("%.1f", median(filesizes)))")
println("std. of filesizes:   $(@sprintf("%.1f", std(filesizes)))")
println()

println("--- xz_compress_serial ---")
xz_compress_serial(files[1:1])
for _ in 1:3
    @time xz_compress_serial(files)
end
println()

println("--- xz_compress_concurrent ---")
xz_compress_concurrent(files[1:1])
for _ in 1:3
    @time xz_compress_concurrent(files)
end

#= 2017/08/29 on tostoint01
-bash-4.1$ find ~/anaconda3/ -type f -name "*.py" | head -3000 | julia -p4 examples/xz.jl
--- statistics ---
number of files:     3000
mean of filesizes:   10472.5
median of filesizes: 3743.5
std. of filesizes:   28214.8

--- xz_compress_serial ---
 25.834263 seconds (141.13 k allocations: 6.016 MiB)
 24.428520 seconds (141.13 k allocations: 6.016 MiB)
 24.663476 seconds (141.13 k allocations: 6.016 MiB)

--- xz_compress_concurrent ---
  6.535346 seconds (13.50 M allocations: 230.372 MiB, 1.61% gc time)
  5.665914 seconds (12.15 M allocations: 206.275 MiB, 1.12% gc time)
  5.747123 seconds (12.41 M allocations: 210.338 MiB, 1.28% gc time)
=#


text:
# NOTE: This script is used for generating data for the PPL comparisons.

import Pkg; Pkg.activate("../../../")

using Distributions
using JSON3
import Random

"""
Function to generate data for univariate Gaussian mixture model (GMM).

Args
====

- `n`: number of observations
- `mu`: Mixture component means
- `sig`: Standard deviations of each component
- `w`: Mixture component weights
- `seed`: Random number generator seed
"""
function generate_gmm_data(; n::Integer,
                           mu::Vector{Float64}, sig::Vector{Float64},
                           w::Vector{Float64},
                           seed=nothing) where T <: AbstractFloat
    if seed != nothing
        Random.seed!(seed)
    end

    # Number of mixture components
    num_components = length(mu)

    # Check that w is a proper probability vector
    @assert all(w .> 0) && sum(w) == 1

    # Generate random class labels
    rand_idx = wsample(1:num_components, w, n)

    # Generate random observations
    y = mu[rand_idx] + randn(n) .* sig[rand_idx]

    return Dict(:rand_idx => rand_idx,
                :mu => mu, :sig => sig, :w => w, :y => y)
end

function main(; seed=1, num_components=4,
              datasizes=[50, 100, 200, 400, 800, 1600],
              data_dir=joinpath(@__DIR__, "../data"))
  Random.seed!(seed)

  # simulate mu
  mu = let
      _mu = collect(1:num_components)
      centered_mu = _mu .- (num_components + 1) / 2
      jittered_mu = centered_mu .+ randn(num_components) * .1
      jittered_mu
  end

  # simulate sigma
  sig = rand(Uniform(.1, .2), num_components)

  # simulate weights
  w = let
      w_tmp = rand(Uniform(.8/num_components, 1.2/num_components),
                   num_components)
      w_tmp ./ sum(w_tmp)
  end

  # Create data sets and write to disk
  for n in datasizes
      # Generate GMM data.
      gmm_data = generate_gmm_data(n=n, mu=mu, sig=sig, w=w, seed=seed)

      # Write GMM data to json string.
      gmm_data_json = JSON3.write(gmm_data);

      # Write json string to disk.
      open("$(data_dir)/gmm-data-n$(n).json", "w") do io
        write(io, gmm_data_json)
      end
  end
end

# Generate data.
main()

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
# This should maybe be a supertype at some point?
const Transform = Union{Validator, DeclareMissings, Filter, Imputor}

"""
    Chain{T<:Tuple{Vararg{Transform}}} <: Function

Runs multiple `Validators`, `Filter` or `Imputor`s on the same data in the order they're
provided.

# Fields
* `transforms::Vector{Union{Validator, Filter, Imputor}}`
"""
struct Chain{T<:Tuple{Vararg{Transform}}} <: Function
    transforms::T
end

Chain(transforms::Vector{<:Transform}) = Chain(Tuple(transforms))

"""
    Chain(transforms::Union{Validator, Filter, Imputor}...) -> Chain

Creates a Chain using the transforms provided (ordering matters).
"""
Chain(transforms::Transform...) = Chain(tuple(transforms...))

"""
Compose new chains with the composition operator

# Example

```jldoctest
julia> using Impute: Impute, Interpolate, NOCB, LOCF

julia> M = [missing 2.0 missing missing 5.0; 1.1 2.2 missing 4.4 missing]
2×5 Matrix{Union{Missing, Float64}}:
  missing  2.0  missing   missing  5.0
 1.1       2.2  missing  4.4        missing

julia> C = Interpolate() ∘ NOCB() ∘ LOCF();

julia> C(M; dims=:rows)
2×5 Matrix{Union{Missing, Float64}}:
 2.0  2.0  3.0  4.0  5.0
 1.1  2.2  3.3  4.4  4.4
```
"""
Base.:(∘)(a::Transform, b::Transform) = Chain(a, b)
Base.:(∘)(C::Chain, b::Transform) = Chain(C.transforms..., b)

"""
    (C::Chain)(data; kwargs...)

Runnable the "callable" chain `C` on the supplied `data`.

# Arguments
* `data`: our data to impute

# Keyword Arguments
* `kwargs`: Keyword arguments that should be applied to each transform (ex `dims=:cols`)

# Returns
* our imputed data
"""
function (C::Chain)(data; kwargs...)
    # Since some operation like filtering can't consistently mutate the data we make a copy
    # and don't support a mutating form.
    X = trycopy(data)

    for t in C.transforms
        if isa(t, Validator)
            # Validators just return the input
            validate(X, t; kwargs...)
        elseif isa(t, DeclareMissings)
            # DeclareMissings isn't guaranteed to work in-place
            X = apply(X, t)
        elseif isa(t, Filter)
            # Filtering doesn't always work in-place
            X = apply(X, t; kwargs...)
        else
            # An in-place impute! method should always exist
            X = impute!(X, t; kwargs...)
        end
    end

    return X
end