text:
#==============================================================================#
# IotData.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/iot-data-2015-05-28.normal.json
#==============================================================================#

__precompile__()

module IotData

using AWSCore

"""
    using AWSSDK.IotData.delete_thing_shadow
    delete_thing_shadow([::AWSConfig], arguments::Dict)
    delete_thing_shadow([::AWSConfig]; thingName=)

    using AWSCore.Services.iot_data
    iot_data([::AWSConfig], "DELETE", "/things/{thingName}/shadow", arguments::Dict)
    iot_data([::AWSConfig], "DELETE", "/things/{thingName}/shadow", thingName=)

# DeleteThingShadow Operation

Deletes the thing shadow for the specified thing.

For more information, see [DeleteThingShadow](http://docs.aws.amazon.com/iot/latest/developerguide/API_DeleteThingShadow.html) in the *AWS IoT Developer Guide*.

# Arguments

## `thingName = ::String` -- *Required*
The name of the thing.

# Returns

`DeleteThingShadowResponse`

# Exceptions

`ResourceNotFoundException`, `InvalidRequestException`, `ThrottlingException`, `UnauthorizedException`, `ServiceUnavailableException`, `InternalFailureException`, `MethodNotAllowedException` or `UnsupportedDocumentEncodingException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-data-2015-05-28/DeleteThingShadow)
"""
@inline delete_thing_shadow(aws::AWSConfig=default_aws_config(); args...) = delete_thing_shadow(aws, args)

@inline delete_thing_shadow(aws::AWSConfig, args) = AWSCore.Services.iot_data(aws, "DELETE", "/things/{thingName}/shadow", args)

@inline delete_thing_shadow(args) = delete_thing_shadow(default_aws_config(), args)

"""
    using AWSSDK.IotData.get_thing_shadow
    get_thing_shadow([::AWSConfig], arguments::Dict)
    get_thing_shadow([::AWSConfig]; thingName=)

    using AWSCore.Services.iot_data
    iot_data([::AWSConfig], "GET", "/things/{thingName}/shadow", arguments::Dict)
    iot_data([::AWSConfig], "GET", "/things/{thingName}/shadow", thingName=)

# GetThingShadow Operation

Gets the thing shadow for the specified thing.

For more information, see [GetThingShadow](http://docs.aws.amazon.com/iot/latest/developerguide/API_GetThingShadow.html) in the *AWS IoT Developer Guide*.

# Arguments

## `thingName = ::String` -- *Required*
The name of the thing.

# Returns

`GetThingShadowResponse`

# Exceptions

`InvalidRequestException`, `ResourceNotFoundException`, `ThrottlingException`, `UnauthorizedException`, `ServiceUnavailableException`, `InternalFailureException`, `MethodNotAllowedException` or `UnsupportedDocumentEncodingException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-data-2015-05-28/GetThingShadow)
"""
@inline get_thing_shadow(aws::AWSConfig=default_aws_config(); args...) = get_thing_shadow(aws, args)

@inline get_thing_shadow(aws::AWSConfig, args) = AWSCore.Services.iot_data(aws, "GET", "/things/{thingName}/shadow", args)

@inline get_thing_shadow(args) = get_thing_shadow(default_aws_config(), args)

"""
    using AWSSDK.IotData.publish
    publish([::AWSConfig], arguments::Dict)
    publish([::AWSConfig]; topic=, <keyword arguments>)

    using AWSCore.Services.iot_data
    iot_data([::AWSConfig], "POST", "/topics/{topic}", arguments::Dict)
    iot_data([::AWSConfig], "POST", "/topics/{topic}", topic=, <keyword arguments>)

# Publish Operation

Publishes state information.

For more information, see [HTTP Protocol](http://docs.aws.amazon.com/iot/latest/developerguide/protocols.html#http) in the *AWS IoT Developer Guide*.

# Arguments

## `topic = ::String` -- *Required*
The name of the MQTT topic.

## `qos = ::Int`
The Quality of Service (QoS) level.

## `payload = blob`
The state information, in JSON format.

# Exceptions

`InternalFailureException`, `InvalidRequestException`, `UnauthorizedException` or `MethodNotAllowedException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-data-2015-05-28/Publish)
"""
@inline publish(aws::AWSConfig=default_aws_config(); args...) = publish(aws, args)

@inline publish(aws::AWSConfig, args) = AWSCore.Services.iot_data(aws, "POST", "/topics/{topic}", args)

@inline publish(args) = publish(default_aws_config(), args)

"""
    using AWSSDK.IotData.update_thing_shadow
    update_thing_shadow([::AWSConfig], arguments::Dict)
    update_thing_shadow([::AWSConfig]; thingName=, payload=)

    using AWSCore.Services.iot_data
    iot_data([::AWSConfig], "POST", "/things/{thingName}/shadow", arguments::Dict)
    iot_data([::AWSConfig], "POST", "/things/{thingName}/shadow", thingName=, payload=)

# UpdateThingShadow Operation

Updates the thing shadow for the specified thing.

For more information, see [UpdateThingShadow](http://docs.aws.amazon.com/iot/latest/developerguide/API_UpdateThingShadow.html) in the *AWS IoT Developer Guide*.

# Arguments

## `thingName = ::String` -- *Required*
The name of the thing.

## `payload = blob` -- *Required*
The state information, in JSON format.

# Returns

`UpdateThingShadowResponse`

# Exceptions

`ConflictException`, `RequestEntityTooLargeException`, `InvalidRequestException`, `ThrottlingException`, `UnauthorizedException`, `ServiceUnavailableException`, `InternalFailureException`, `MethodNotAllowedException` or `UnsupportedDocumentEncodingException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/iot-data-2015-05-28/UpdateThingShadow)
"""
@inline update_thing_shadow(aws::AWSConfig=default_aws_config(); args...) = update_thing_shadow(aws, args)

@inline update_thing_shadow(aws::AWSConfig, args) = AWSCore.Services.iot_data(aws, "POST", "/things/{thingName}/shadow", args)

@inline update_thing_shadow(args) = update_thing_shadow(default_aws_config(), args)

end # module IotData

#==============================================================================#
# End of file
#==============================================================================#


text:
using Gadfly, Distributions

set_default_plot_size(6inch, 3inch)

vals = rand(MultivariateNormal([0.0, 0.0], [1.0 0.5; 0.5 1.0]), 10000)
Gadfly.with_theme(:dark) do
    plot(x=vals[1,:], y=vals[2,:], Geom.hexbin)
end

ext: jl


text:
using Printf
function set_initQbase(xmax, ymax, restart_file, init_rho, init_u, init_v, init_p, init_N2, init_N, dt,
                        specific_heat_ratio, out_file_front, out_ext, out_dir, restart_num, R, nch, nval)
    Qbase=[]
    cellxmax = xmax - 1
    cellymax = ymax - 1

    init_step = 0
    init_time = 0.0

    restart_check = 0
    try Qbase = setup_restart_value(cellxmax, cellymax, out_dir, restart_file, nch, nval)
        println("Restart "*restart_file)
        restart_check = 2
        init_step = restart_num
        init_time = init_step*dt
    catch
        restart_check = 1
    end

    if restart_check == 1
        Qbase = setup_init_value(cellxmax, cellymax, init_rho, init_u, init_v, init_p, init_N2, init_N, nval)
        println("Start Initial condition")
        restart_num = 0
        Rhat = set_gasconst_hat(Qbase, cellxmax, cellymax, nval, nch, R)
        output_result(0, Qbase, cellxmax, cellymax, specific_heat_ratio, out_file_front, out_ext, out_dir, Rhat, nval)
    end

    return Qbase, cellxmax, cellymax, restart_num, init_step, init_time
end

function setup_init_value(cellxmax, cellymax, init_rho, init_u, init_v, init_p, init_N2, init_N, nval)
    Qbase=zeros(cellxmax, cellymax, nval)
    for i in 1:cellxmax
        for j in 1:cellymax
            Qbase[i,j,1] = init_rho
            Qbase[i,j,2] = init_u
            Qbase[i,j,3] = init_v
            Qbase[i,j,4] = init_p
            Qbase[i,j,5] = init_N2
            Qbase[i,j,6] = init_N
        end
    end
    return Qbase
end

function setup_restart_value(cellxmax, cellymax, out_dir, restart_file, nch, nval)
    Qbase = zeros(cellxmax, cellymax, nval)

    skipnum = 1
    fff = []
    open("result/"*restart_file, "r") do f
        fff = read(f, String)
    end
    fff = split(fff,"\n", keepempty = false)

    for i in 1+skipnum:length(fff)
        fff[i] = replace(fff[i]," \r" => "")
    end

    ite = 1
    npre = nval - nch
    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            temp = split(fff[ite+skipnum]," ")
            for l in 1:nval
                Qbase[i,j,l] = parse(Float64,temp[l])
            end
            for l in 1:nch # mass frac
                Qbase[i,j,npre+l] = Qbase[i,j,npre+l]
            end
            ite = ite+1
        end
    end
    return Qbase
end

function check_divrege(Qbase, cellxmax, cellymax, Rhat, fwrite)
    ite = 0
    for i in 1:cellxmax
        for j in 1:cellymax
            T = Qbase[i,j,4]/(Qbase[i,j,1]*Rhat[i,j])
            if T < 0 || isequal(T, NaN) == true
                open( fwrite, "a" ) do f
                    ai = @sprintf("%4.0f", i)
                    aj = @sprintf("%4.0f", j)

                    write(f, "\n")
                    write(f, " thermal diverge ")
                    write(f, "\n")
                    write(f, " i = "*ai)
                    write(f, "\n")
                    write(f, " j = "*aj)
                    write(f, "\n")
                end
                ite = 1
            end
            if 1+1e-3 < Qbase[i,j,5] || Qbase[i,j,5] < -1e-3
                open( fwrite, "a" ) do f
                    ai = @sprintf("%4.0f", i)
                    aj = @sprintf("%4.0f", j)

                    write(f, "\n")
                    write(f, " chemical diverge ")
                    write(f, "\n")
                    write(f, " i = "*ai)
                    write(f, "\n")
                    write(f, " j = "*aj)
                    write(f, "\n")
                end
                ite = 2
            end
        end
    end

    if ite == 1
        println("\n")
        println(" T < 0 or T = NaN ")
        println(" thermal diverge ")
        println("\n")
        throw(UndefVarError(:x))
    elseif ite == 2
        println("\n")
        println(" N > 1 or N < 0 ")
        println(" chemical diverge ")
        println("\n")
        throw(UndefVarError(:x))
    end
end

function check_small(Qbase, cellxmax, cellymax, nch)
    for i in 1:cellxmax
        for j in 1:cellymax
            for ns in 1:nch
                if Qbase[i,j,4+ns] < 1e-20
                    Qbase[i,j,4+ns] = 0.0
                end
            end
        end
    end
end

function GE(matrix, nval)
    # nval * nval 行列の逆行列を求める。

    for l in 1:nval
        a = 1.0/matrix[l,l]
        matrix[l,l] = 1.0

        for m in 1:nval
            matrix[l,m] = matrix[l,m]*a
        end

        for m in 1:nval
            if m == l
            else
                b = matrix[m,l]
                matrix[m,l] = 0.0
                for n in 1:nval
                    matrix[m,n] = matrix[m,n] - b*matrix[l,n]
                end
            end
        end
    end
    return matrix
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
# Deprecated during 0.7.*

@deprecate principal_gen principal_generator

@deprecate principal_gen_fac_elem principal_generator_fac_elem

@deprecate right_principal_gen right_principal_generator

@deprecate left_principal_gen left_principal_generator

@deprecate principal_gen_eichler principal_generator_eichler

# Deprecated during 0.10.*

@deprecate subgroup sub

@deprecate upper_bound(x::arb, y::Type{fmpz}) upper_bound(y, x)

@deprecate abs_upper_bound(x::arb, y::Type{fmpz}) abs_upper_bound(y, x)

@deprecate isequivalent isisometric

@deprecate isequivalent_with_isometry isisometric_with_isometry

# Things that moved to Nemo

# > 0.18.1
if isdefined(Nemo, :simplest_between)
  simplest_inside(x::arb) = simplest_rational_inside(x)
else
  function _fmpq_simplest_between(l_num::fmpz, l_den::fmpz,
                                  r_num::fmpz, r_den::fmpz)
     n = fmpz()
     d = fmpz()

     ccall((:_fmpq_simplest_between, libflint), Nothing,
           (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{fmpz}),
           n, d, l_num, l_den, r_num, r_den)

     return n//d
  end

  @doc Markdown.doc"""
        simplest_between(l::fmpq, r::fmpq)

  > Return the simplest fraction in the closed interval `[l, r]`. A canonical >
  > fraction `a_1/b_1` is defined to be simpler than `a_2/b_2` iff `b_1 < b_2` or
  > `b_1 = b_2` and `a_1 < a_2`.
  """
  function simplest_between(l::fmpq, r::fmpq)
     z = fmpq()
     ccall((:fmpq_simplest_between, libflint), Nothing,
           (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, l, r)
     return z
  end

  function simplest_inside(x::arb)
    a = fmpz()
    b = fmpz()
    e = fmpz()

    ccall((:arb_get_interval_fmpz_2exp, libarb), Nothing, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}, Ref{arb}), a, b, e, x)
    @assert fits(Int, e)
    _e = Int(e)
    if e >= 0
      return a << _e
    end
    _e = -_e
    d = one(fmpz) << _e
    return _fmpq_simplest_between(a, d, b, d)
  end
end


text:
using RandomTestMatrices
using Base.Test

d = rand(1:100)

let A = randsym(d)
    @Test.test isreal(A)
    @Test.test size(A) == (d,d)
    @Test.test issym(A)
end

let A = randantisym(d)
    i = rand(1:d)
    j = rand(1:d)
    @Test.test isreal(A)
    @Test.test size(A) == (d,d)
    @Test.test A[i,j] == -A[j,i]
end

let A = randposdef(d)
    @Test.test isreal(A)
    @Test.test size(A) == (d,d)
    @Test.test isposdef(A)
    @Test.test issym(A)
end

let A = randstable(d)
    @Test.test size(A) == (d,d)
    @Test.test isreal(A)
    @Test.test all(real(eig(A)[1]) .< 0.0)
end

let A = randunitary(d)
    @Test.test size(A) == (d,d)
    @Test.test !isreal(A)
    @Test.test norm(A*A' - eye(d)) < d*eps()
end

let A = randorth(d)
    @Test.test size(A) == (d,d)
    @Test.test isreal(A)
    @Test.test norm(A*A' - eye(d)) < d*eps()
end

let A = randnormal(d)
    @Test.test size(A) == (d,d)
    @Test.test isreal(A)
    @Test.test norm(A*A' - A'*A) < d*eps(100.)
end


text:
module PendulumWithBar1

using ModiaLang

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

import Modia3D
using  Modia3D.ModiaInterface

Bar = Model(
    m  = 0.1,
    Lx = 0.1,
    Ly = Par(value=:(0.2*Lx)),
    Lz = Par(value=:(0.2*Lx)),
    vmat1 = VisualMaterial(color="DeepSkyBlue2", transparency=0.5),
    vmat2 = VisualMaterial(color="Red"),
    frame0 = Object3D(feature=Solid(shape=Beam(axis=1, length=:Lx, width=:Ly, thickness=:Lz),
                                    massProperties=MassProperties(mass=:m),
                                    visualMaterial=:(vmat1))),
	frame1 = Object3D(parent=:frame0,
                      translation=:[-Lx/2, 0.0, 0.0],
                      feature=Visual(shape=Cylinder(axis=3, diameter=:(Ly/2), length=:(1.2*Lz)),
                                     visualMaterial=:(vmat2)))
)

Pendulum = Model(
    m = 1.0,
    Lx = 0.1,
    world = Object3D(feature=Scene() ),
    worldFrame = Object3D(parent=:world,
                          feature=Visual(shape=CoordinateSystem(length=:(Lx/2)))),
    bar = Bar | Map(m=:m, Lx=:Lx),
    rev = RevoluteWithFlange(obj1=:world, obj2=:(bar.frame1))
)

PendulumWithBar = Model(
    pendulum = buildModia3D(Pendulum | Map(Lx=1.0, m=2.0, rev=Map(phi=Var(init=1.0)))),

    damper = Damper | Map(d=0.5),
    fixed = Fixed,
    connect = :[(damper.flange_b, pendulum.rev.flange),
                (damper.flange_a, fixed.flange)]
)

pendulumWithBar = @instantiateModel(PendulumWithBar, unitless=true)

import DifferentialEquations
algorithm = DifferentialEquations.Tsit5()
stopTime = 10.0
requiredFinalStates = [-1.578178763749515, 0.06153191687388868]
simulate!(pendulumWithBar, algorithm, stopTime=stopTime, log=true, logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plot(pendulumWithBar, ["pendulum.rev.flange.phi", "pendulum.rev.variables[1]"], figure=1)

end


text:
# Autogenerated wrapper script for protoc_jll for x86_64-w64-mingw32-cxx11
export libprotoc, protoc

JLLWrappers.@generate_wrapper_header("protoc")
JLLWrappers.@declare_library_product(libprotoc, "libprotoc-30.dll")
JLLWrappers.@declare_executable_product(protoc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libprotoc,
        "bin\\libprotoc-30.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        protoc,
        "bin\\protoc.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Base.Iterators: product

abstract type Partitioning{N} <: AbstractArray{Int, N} end

"""
Store the distribution of the array indices over the different partitions.
This class forces a continuous, ordered distribution without overlap

    ContinuousPartitioning(partition_sizes...)

Construct a distribution using the number of elements per partition in each direction, e.g:

```
p = ContinuousPartitioning([2,5,3], [2,3])
```

will construct a distribution containing 6 partitions of 2, 5 and 3 rows and 2 and 3 columns.

"""
struct ContinuousPartitioning{N} <: Partitioning{N}
    ranks::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
    index_starts::NTuple{N,Vector{Int}}
    index_ends::NTuple{N,Vector{Int}}

    function ContinuousPartitioning(partition_sizes::Vararg{Any,N}) where {N}
        index_starts = Vector{Int}.(undef,length.(partition_sizes))
        index_ends = Vector{Int}.(undef,length.(partition_sizes))
        for (idxstart,idxend,nb_elems_dist) in zip(index_starts,index_ends,partition_sizes)
            currentstart = 1
            currentend = 0
            for i in eachindex(idxstart)
                currentend += nb_elems_dist[i]
                idxstart[i] = currentstart
                idxend[i] = currentend
                currentstart += nb_elems_dist[i]
            end
        end
        ranks = LinearIndices(length.(partition_sizes))
        return new{N}(ranks, index_starts, index_ends)
    end
end

"""
    CyclicPartitioning for ScaLAPACK.jl

```
p = CyclicPartitioning([2,5,3], [2,3], blocksize=(1, 1))
```

will construct a cyclic distribution containing 6 partitions of 2, 5 and 3 rows and 2 and 3 columns.

N is dimension of partition

"""
struct CyclicPartitioning{N} <: Partitioning{N}
    # rank partition mapping
    ranks::LinearIndices{N,NTuple{N,Base.OneTo{Int}}}
    partitions::Array{NTuple{N, CyclicRange}, N}

    function CyclicPartitioning(;array_sizes::NTuple{N, Integer}, n_procs::NTuple{N, Integer}, blocksizes::NTuple{N, Integer}) where {N}
        # process grid array
        ranks = LinearIndices(n_procs)
        # indices array
        partitions = Array{NTuple{N, CyclicRange}, N}(undef, n_procs...)

        for process_indices in product(Base.OneTo.(n_procs)...)
            cyclic_ranges = CyclicRange[]
            for (dim, array_size, n_proc, blocksize) in zip(1:N, array_sizes, n_procs, blocksizes)
                # Index of grid of process in dim dimension
                cyclic_range = CyclicRange(start=Int(1), stop=Int(array_size), blocksize=blocksize, n_split=Int(n_proc), split_index=Int(process_indices[dim]))
                push!(cyclic_ranges, cyclic_range)
            end
            partitions[process_indices...] = Tuple(cyclic_ranges)
        end
        # if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        # show(stdout, "text/plain", map(x->length.(x), partitions))
        # end
        return new{N}(ranks, partitions)
    end
end

Base.IndexStyle(::Type{T}) where T <: Partitioning = IndexCartesian()
# number of grid process each dimension
Base.size(p::T) where T <: Partitioning = size(p.ranks)

@inline function Base.getindex(p::ContinuousPartitioning{N}, I::Vararg{Int, N}) where N
    return UnitRange.(getindex.(p.index_starts,I), getindex.(p.index_ends,I))
end

@inline function Base.getindex(p::CyclicPartitioning{N}, I::Vararg{Int, N}) where N
    return p.partitions[I...]
end

function partition_sizes(p::ContinuousPartitioning)
    result = (p.index_ends .- p.index_starts)
    for v in result
        v .+= 1
    end
    return result
end

# Cyclic partion_sizes implemented only 1d and 2d
function partition_sizes(p::CyclicPartitioning{1})
    grids = map(x->length(p.partitions[x][1]), p.ranks[:, 1])
    return (grids, )
end
function partition_sizes(p::CyclicPartitioning{2})
    n_row = map(x->length(p.partitions[x][1]), p.ranks[:, 1])
    n_col = map(x->length(p.partitions[x][2]), p.ranks[1, :])
    result = (n_row, n_col)
    return result
end

"""
  (private method)

Get the rank and local 0-based index
"""
function local_index(p::ContinuousPartitioning, I::NTuple{N,Int}) where {N}
    proc_indices = searchsortedfirst.(p.index_ends, I)
    lininds = LinearIndices(Base.Slice.(p[proc_indices...]))
    return (p.ranks[proc_indices...]-1, lininds[I...] - first(lininds))
end

# 1-dimentional
function local_index(p::CyclicPartitioning{1}, I::NTuple{1,Int})
    cyclic_range = p.partitions[1][1]
    # map global index to local array
    # global_index_mapped = product(cyclic_range...)
    rank_indices = local_index(cyclic_range, I[1])[1]
    local_indices = local_index(cyclic_range, I[1])[2]
    # to 0-based index (rank, indices of local matrix)
    return p.ranks[rank_indices]-1, local_indices-1
end

# 2-dimentional
function local_index(p::CyclicPartitioning{2}, I::NTuple{2,Int})
    # return local indices as 1 dimensional
    rank, (i, j) = local_index_xd(p, I)
    n_local_row, n_local_col = length.(p.partitions[rank+1])
    localindex = i + j*n_local_row
    return rank, localindex
end

function local_index_xd(p::CyclicPartitioning, I::NTuple{N,Int}) where {N}
    # return local indices as Tuple
    # all partitions have same cyclic information
    cyclic_range = p.partitions[1]
    # map global index to local array
    # global_index_mapped = product(cyclic_range...)
    rank_indices = map(x->x[1], local_index.(cyclic_range, I))
    local_indices = map(x->x[2], local_index.(cyclic_range, I))
    # to 0-based index (rank, indices of local matrix)
    return p.ranks[rank_indices...]-1, local_indices .- 1
end

function global_size(p::CyclicPartitioning{1})
    n_index = sum(map(x->length(p.partitions[x]), p.ranks))
    return n_index
end

function global_size(p::CyclicPartitioning{2})
    n_row = sum(map(x->length(p.partitions[x][1]), p.ranks[:, 1]))
    n_col = sum(map(x->length(p.partitions[x][2]), p.ranks[1, :]))
    return (n_row, n_col)
end

blocksizes(p::CyclicPartitioning) = map(x->x.blocksize, p.partitions[1])


text:
"""
    cl5(x::Float64)::Float64

Returns the value of the Clausen function ``\\operatorname{Cl}_5(x)``
for a real angle ``x`` of type `Float64`.  This function is defined as

```math
\\operatorname{Cl}_5(x) = \\Re[\\operatorname{Li}_5(e^{ix})] = \\sum_{k=1}^\\infty \\frac{\\cos(kx)}{k^5}
```

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> cl5(1.0)
0.5228208076420943
```
"""
function cl5(x::Float64)::Float64
    zeta5 = 1.0369277551433699
    pi28 = pi*pi/8.0

    x = range_reduce_odd(x)

    if x == zero(x)
        return zeta5
    end

    if x < 0.5*pi
        P = (1.0369277551433699e+00, -6.1354800479984468e-01,
             9.4076401395712763e-02, -9.4056155866704436e-04)
        Q = (1.0000000000000000e+00, -1.2073698633244778e-02,
             1.3703409625482991e-05, -1.9701280330628469e-09,
             2.1944550184416500e-11)
        y = x*x
        y2 = y*y
        p = P[1] + y * P[2] + y2 * (P[3] + y * P[4])
        q = Q[1] + y * Q[2] + y2 * (Q[3] + y * Q[4] + y2 * Q[5])
        p/q - 1/24*y2*log(x)
    else
        P = (-4.5930112735784898e-01, 4.3720705508867954e-01,
             -7.5895226486465095e-02, 5.2244176912488065e-03,
             -1.5677716622013956e-04, 1.6641624171748576e-06)
        Q = ( 1.0000000000000000e+00, -1.2211486825401188e-01,
              3.8940070749313620e-03, -2.2674805547074318e-05,
             -7.4383354448335299e-08, -3.4131758392216437e-10)
        y = pi - x
        z = y*y - pi28
        z2 = z*z
        z4 = z2*z2
        p = P[1] + z * P[2] + z2 * (P[3] + z * P[4]) + z4 * (P[5] + z * P[6])
        q = Q[1] + z * Q[2] + z2 * (Q[3] + z * Q[4]) + z4 * (Q[5] + z * Q[6])
        p/q
    end
end


text:
# Use to provide proof that everything works
include("src/main.jl")

BUSES, LINES, GENERATORS, rPTDF = load_case_data(datafile="basecase")
# BUSES, LINES, GENERATORS, rPTDF = load_case_data(datafile="33buscase_pu")

η_v = .05 # 1 - Confidence for voltage limit cc
η_g = .05 # 1 - Confidence for generation limit cc
ξ = .0005# 1 - Confidence for variance estimator

# SND Quantiles for global use
snd = Normal(0,1)
z_v = quantile(snd, 1-η_v)
z_g = quantile(snd, 1-η_g)
z_g2 = quantile(snd, 1-(η_g/2))

# Set True Load Standard Deviation
ld = 0.2

ERROR_VARIANCES_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE =
    init_stochasticity(variance_opt="implicit")

# To test feasability
ld_objective, ld_bus_results, ld_line_results = Solve_LinDist(BUSES, LINES, GENERATORS)
info("Results of standard lindist. Objective: $ld_objective")
println(ld_bus_results)
println(ld_line_results)

# To test SOC
# socobjective, socbus_results, socline_results = Solve_Dist_SOC(BUSES, LINES, GENERATORS)
# info("Results of soc dist. Objective: $fobjective")
# println(fbus_results)
# println(fline_results)

# linobjective, linbus_results, linline_results = Solve_Linearized_Dist_SOC(BUSES, LINES, GENERATORS,
#                         socline_results[:f_P], socline_results[:f_Q], socline_results[:a_squared])
# info("Results of soc dist. Objective: $fobjective")
# println(fbus_results)
# println(fline_results)

# println("SOC Objective = $(socobjective)")
# println("LIN Objective = $(linobjective)")
# println("LinDist Objective = $(ld_objective)")
# compare_voltages(socbus_results, linbus_results, ld_bus_results, label1="soc", label2="linearized", label3="lindist")
# compare_voltages(socbus_results, linbus_results, label1="soc", label2="linearized")
# compare_voltages(ld_bus_results, linbus_results, label1="lindist", label2="linearized")
# compare_voltages(fbus_results, b_ref, label1="soc_dist", label2="ref" )
# compare_voltages(fbus_results, ld_bus_results, label1="soc_dist", label2="lin_dist" )

ccobjective, ccbus_results, ccline_results = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
info("Results of standard cc. Objective: $ccobjective")
println(ccbus_results)
println(ccline_results)

objective, bus_results, line_results = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_WC)
info("Results of dd cc. Objective: $objective")
println(bus_results)
println(line_results)

sample_size = 750
# Create global set of true observations to make tests comparable
OBSERVATIONS = create_data(ERROR_VARIANCES_TRUE, sample_size)
test_results(ccbus_results, ccline_results, ccobjective; samples = sample_size, message="Standard CC")
test_results(bus_results, line_results, objective; samples = sample_size, message="DD CC")


text:
using StaticArrays, LinearAlgebra

const SVec3 = SVector{3}

export compute_y, compute_y!, cYlm_from_cart, cYlm_from_cart!, index_y

"""
	sizeP(maxDegree)

Return the size of the set of Associated Legendre Polynomials ``P_l^m(x)`` of
degree less than or equal to the given maximum degree
"""
sizeP(maxDegree) = div((maxDegree + 1) * (maxDegree + 2), 2)

"""
	sizeY(maxDegree)

Return the size of the set of real spherical harmonics ``Y_{l,m}(θ,φ)`` of
degree less than or equal to the given maximum degree
"""
sizeY(maxDegree) = (maxDegree + 1) * (maxDegree + 1)

"""
	index_p(l,m)

Return the index into a flat array of Associated Legendre Polynomials ``P_l^m``
for the given indices ``(l,m)``.
``P_l^m`` are stored in l-major order i.e. [P(0,0), [P(1,0), P(1,1), P(2,0), ...]
"""
index_p(l,m) = m + div(l*(l+1), 2) + 1

"""
	index_y(l,m)

Return the index into a flat array of real spherical harmonics ``Y_{l,m}``
for the given indices ``(l,m)``.
``Y_{l,m}`` are stored in l-major order i.e.
[Y(0,0), [Y(1,-1), Y(1,0), Y(1,1), Y(2,-2), ...]
"""
index_y(l,m) = m + l + (l*l) + 1

"""
TODO: documentation
"""
struct ALPCoefficients
	A::Array{Float64}
	B::Array{Float64}
end

ALPCoefficients(maxDegree::Int) =
	ALPCoefficients( Array{Float64}(undef, sizeP(maxDegree)),
						  Array{Float64}(undef, sizeP(maxDegree)) )

"""
	compute_coefficients(L)

Precompute coefficients ``a_l^m`` and ``b_l^m`` for all l <= L, m <= l
"""
function compute_coefficients(L::Int)
	coeff = ALPCoefficients(L)
	for l in 2:L
		ls = l*l
		lm1s = (l-1) * (l-1)
		for m in 0:(l-2)
			ms = m * m
			coeff.A[index_p(l, m)] = sqrt((4 * ls - 1.0) / (ls - ms))
			coeff.B[index_p(l, m)] = -sqrt((lm1s - ms) / (4 * lm1s - 1.0))
		end
	end
	return coeff
end

"""
	compute_coefficients(L)

Create an array large enough to store an entire set of Associated Legendre
Polynomials ``P_l^m(x)`` of maximum degree L.
"""
allocate_p(L::Int) = Array{Float64}(undef, sizeP(L))

"""
	compute_p(L, x, coeff, P)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)``
using the given coefficients, and store in the array P.
"""
function compute_p!(L::Int, x::Float64, coeff::ALPCoefficients,
					     P::Array{Float64,1})
	@assert length(coeff.A) >= sizeP(L)
	@assert length(coeff.B) >= sizeP(L)
	@assert length(P) >= sizeP(L)

	sintheta = sqrt(1.0 - x * x)
	temp = 0.39894228040143267794 # = sqrt(0.5/M_PI)
	P[index_p(0, 0)] = temp

	if (L > 0)
		SQRT3 = 1.7320508075688772935
		P[index_p(1, 0)] = x * SQRT3 * temp
		SQRT3DIV2 = -1.2247448713915890491
		temp = SQRT3DIV2 * sintheta * temp
		P[index_p(1, 1)] = temp

		for l in 2:L
			for m in 0:(l-2)
				P[index_p(l, m)] = coeff.A[index_p(l, m)] *(x * P[index_p(l - 1, m)]
						     + coeff.B[index_p(l, m)] * P[index_p(l - 2, m)])
			end
			P[index_p(l, l - 1)] = x * sqrt(2 * (l - 1) + 3) * temp
			temp = -sqrt(1.0 + 0.5 / l) * sintheta * temp
			P[index_p(l, l)] = temp
		end
	end
	return P
end

"""
	compute_p(L, x)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)`` where
``0 ≤ l ≤ L`` and ``0 ≤ m ≤ l``. Assumes ``|x| ≤ 1``.
"""
function compute_p(L::Int, x::Float64)
	P = Array{Float64}(undef, sizeP(L))
	coeff = compute_coefficients(L)
	compute_p!(L, x, coeff, P)
	return P
end

"""
	compute_y(L, P, φ, Y)

Compute an entire set of real spherical harmonics ``Y_{l,m}(θ,φ)``
using the given Associated Legendre Polynomials ``P_l^m(cos θ)``
and store in array Y
"""
function compute_y!(L::Int, x::Float64, phi::Float64,
						  P::Array{Float64,1},
						  Y::Array{ComplexF64,1})
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)

	SQRT2 = 1.41421356237309504880
	for l in 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * 0.5 * SQRT2
	end

	# ------ real spherical harmonics code ----------------
	# NR2 5.5.4-5.5.5
	# c1 = 1.0; c2 = cos(phi)
	# s1 = 0.0; s2 = -sin(phi)
	# tc = 2.0 * c2
	# for m in 1:L
	# 	s = tc * s1 - s2
	# 	c = tc * c1 - c2
	# 	s2 = s1
	# 	s1 = s
	# 	c2 = c1
	# 	c1 = c
	# 	for l in m:L
	# 		Y[index_y(l, -m)] = P[index_p(l, m)] * s
	# 		Y[index_y(l, m)] = P[index_p(l, m)] * c
	# 	end
	# end
	# ------------------------------------------------------------

   sig = 1
	for m in 1:L
		sig *= -1
		ep = exp(im*m*phi) / sqrt(2)
		em = sig * conj(ep)
		for l in m:L
			p = P[index_p(l,m)]
			# Y[index_y(l, -m)] = (-1)^m * p * exp(-im*m*phi) / sqrt(2)
			# Y[index_y(l,  m)] = p * exp( im*m*phi) / sqrt(2)
			Y[index_y(l, -m)] = em * p
			Y[index_y(l,  m)] = ep * p
		end
	end

	return Y
end

"""
	compute_y(L, x, φ)

Compute an entire set of real spherical harmonics ``Y_{l,m}(θ,φ)`` for
``x = cos θ`` where ``0 ≤ l ≤ L`` and ``-l ≤ m ≤ l``.
"""
function compute_y(L::Int, x::Float64, phi::Float64)
	P = Array{Float64}(undef, sizeP(L))
	coeff = compute_coefficients(L)
	compute_p!(L, x, coeff, P)
	Y = Array{ComplexF64}(undef, sizeY(L))
	compute_y!(L, x, phi, P, Y)
	return Y
end

# ------------------------------------------------------------------------
#                  NEW CARTESIAN IMPLEMENTATION
# ------------------------------------------------------------------------

function compute_rxz(R::SVec3{T}) where {T}
   r = norm(R)
   z = R[3] / r
   x = R[2] / sqrt(1 - z^2) / r
   s = sign(R[1])
   return r, x, z, s
end

function cYlm_from_cart!(Y, L, r, x, z, s, P)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(z) <= 1.0

	SQRT2_DIV_2 = 0.5 * 1.41421356237309504880
	INVSQRT2 = 1 / 1.41421356237309504880

	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * SQRT2_DIV_2
	end

   sig = 1
   ep = INVSQRT2
   ep_fact = s * sqrt(1-x^2) + im * x
	for m in 1:L
		sig *= -1
		ep *= ep_fact        # ep =   exp(i *   m  * φ)
		em = sig * conj(ep)  # ep = ± exp(i * (-m) * φ)
		for l in m:L
			p = P[index_p(l,m)]
			Y[index_y(l, -m)] = em * p   # (-1)^m * p * exp(-im*m*phi) / sqrt(2)
			Y[index_y(l,  m)] = ep * p   # p * exp( im*m*phi) / sqrt(2)
		end
	end

	return Y
end

"""
	cYlm_from_xz(L, x, z)

Compute an entire set of real spherical harmonics ``Y_{l,m}(θ, φ)`` for
``x = cos θ, z = sin φ`` where ``0 ≤ l ≤ L`` and ``-l ≤ m ≤ l``.
"""
function cYlm_from_cart(L::Integer, R::SVec3{T}) where {T}
   r, x, z, s = compute_rxz(R)
	P = Vector{T}(undef, sizeP(L))
	coeff = compute_coefficients(L)
	compute_p!(L, z, coeff, P)
	Y = Vector{ComplexF64}(undef, sizeY(L))
	cYlm_from_cart!(Y, L, r, x, z, s, P)
	return Y
end

# ------- Experimental

function rYlm_from_cart!(Y, L, r, x, z, s, P)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(z) <= 1.0

	SQRT2_DIV_2 = 0.5 * 1.41421356237309504880
	INVSQRT2 = 1 / 1.41421356237309504880

	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * SQRT2_DIV_2
	end

   sig = 1
   cosφ = x
   sinφ = s * sqrt(1-x^2)
   cosmφ = 1.0
   sinmφ = 0.0
	for m in 1:L
		sig *= -1
      # recursive formula for trigs
      cosmφ, sinmφ = (cosmφ * cosφ - sinmφ * sinφ), (cosmφ * sinφ + sinmφ *  cosφ)
		for l in m:L
			p = P[index_p(l,m)]
			Y[index_y(l, -m)] = sig * p * cosmφ
			Y[index_y(l,  m)] =       p * sinmφ
		end
	end

	return Y
end

	# ------ real spherical harmonics code ----------------
	# NR2 5.5.4-5.5.5
	# c1 = 1.0; c2 = cos(phi)
	# s1 = 0.0; s2 = -sin(phi)
	# tc = 2.0 * c2
	# for m in 1:L
	# 	s = tc * s1 - s2
	# 	c = tc * c1 - c2
	# 	s2 = s1
	# 	s1 = s
	# 	c2 = c1
	# 	c1 = c
	# 	for l in m:L
	# 		Y[index_y(l, -m)] = P[index_p(l, m)] * s
	# 		Y[index_y(l, m)] = P[index_p(l, m)] * c
	# 	end
	# end
	# ------------------------------------------------------------


text:
module SpringCollab2020GroupedTriggerSpikes

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/GroupedTriggerSpikesUp" GroupedTriggerSpikesUp(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, type::String="default", behindMoveBlocks::Bool=false)
@mapdef Entity "SpringCollab2020/GroupedTriggerSpikesDown" GroupedTriggerSpikesDown(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, type::String="default", behindMoveBlocks::Bool=false)
@mapdef Entity "SpringCollab2020/GroupedTriggerSpikesLeft" GroupedTriggerSpikesLeft(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, type::String="default", behindMoveBlocks::Bool=false)
@mapdef Entity "SpringCollab2020/GroupedTriggerSpikesRight" GroupedTriggerSpikesRight(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, type::String="default", behindMoveBlocks::Bool=false)

const placements = Ahorn.PlacementDict()

groupedTriggerSpikes = Dict{String, Type}(
    "up" => GroupedTriggerSpikesUp,
    "down" => GroupedTriggerSpikesDown,
    "left" => GroupedTriggerSpikesLeft,
    "right" => GroupedTriggerSpikesRight
)

groupedTriggerSpikesUnion = Union{GroupedTriggerSpikesUp, GroupedTriggerSpikesDown, GroupedTriggerSpikesLeft, GroupedTriggerSpikesRight}

for variant in Maple.spike_types
    if variant != "tentacles"
        for (dir, entity) in groupedTriggerSpikes
            key = "Grouped Trigger Spikes ($(uppercasefirst(dir)), $(uppercasefirst(variant))) (Spring Collab 2020)"
            placements[key] = Ahorn.EntityPlacement(
                entity,
                "rectangle",
                Dict{String, Any}(
                    "type" => variant
                )
            )
        end
    end
end

Ahorn.editingOptions(entity::groupedTriggerSpikesUnion) = Dict{String, Any}(
    "type" => String[variant for variant in Maple.spike_types if variant != "tentacles"]
)

directions = Dict{String, String}(
    "SpringCollab2020/GroupedTriggerSpikesUp" => "up",
    "SpringCollab2020/GroupedTriggerSpikesDown" => "down",
    "SpringCollab2020/GroupedTriggerSpikesLeft" => "left",
    "SpringCollab2020/GroupedTriggerSpikesRight" => "right",
)

offsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (4, -4),
    "down" => (4, 4),
    "left" => (-4, 4),
    "right" => (4, 4),
)

groupedTriggerSpikesOffsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (0, 5),
    "down" => (0, -4),
    "left" => (5, 0),
    "right" => (-4, 0),
)

rotations = Dict{String, Number}(
    "up" => 0,
    "right" => pi / 2,
    "down" => pi,
    "left" => pi * 3 / 2
)

resizeDirections = Dict{String, Tuple{Bool, Bool}}(
    "up" => (true, false),
    "down" => (true, false),
    "left" => (false, true),
    "right" => (false, true),
)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::groupedTriggerSpikesUnion)
    direction = get(directions, entity.name, "up")
    theta = rotations[direction] - pi / 2

    width = Int(get(entity.data, "width", 0))
    height = Int(get(entity.data, "height", 0))

    x, y = Ahorn.position(entity)
    cx, cy = x + floor(Int, width / 2) - 8 * (direction == "left"), y + floor(Int, height / 2) - 8 * (direction == "up")

    Ahorn.drawArrow(ctx, cx, cy, cx + cos(theta) * 24, cy + sin(theta) * 24, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        x, y = Ahorn.position(entity)

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        direction = get(directions, entity.name, "up")

        ox, oy = offsets[direction]

        return Ahorn.Rectangle(x + ox - 4, y + oy - 4, width, height)
    end
end

Ahorn.minimumSize(entity::groupedTriggerSpikesUnion) = 8, 8

function Ahorn.resizable(entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        direction = get(directions, entity.name, "up")

        return resizeDirections[direction]
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::groupedTriggerSpikesUnion)
    if haskey(directions, entity.name)
        variant = get(entity.data, "type", "default")
        direction = get(directions, entity.name, "up")
        groupedTriggerSpikesOffset = groupedTriggerSpikesOffsets[direction]

        width = get(entity.data, "width", 8)
        height = get(entity.data, "height", 8)

        for ox in 0:8:width - 8, oy in 0:8:height - 8
            drawX, drawY = (ox, oy) .+ offsets[direction] .+ groupedTriggerSpikesOffset
            Ahorn.drawSprite(ctx, "danger/spikes/$(variant)_$(direction)00", drawX, drawY)
        end
    end
end

end


text:
using ExtremalOptimization
using Test, Random

@testset "ExtremalOptimization.jl" begin
    sol = optimize(x -> (x - 1)^2, i -> randn(), 50)
    @test sol.x ≈ 1 rtol = 1e-4
    sol = optimize(
        x -> (x[1] - 1)^2 + (x[2] - x[1]^2)^2,
        i -> randn(2),
        50,
        reps_per_particle = 200,
    )
    map(sol.x) do x
        @test x ≈ 1 rtol = 1e-4
    end
    sol = optimize(x -> ifelse(1.2 < x, Inf, (x - 1)^2), i -> rand() * 0.2 + 0.9, 50)
    @test sol.x ≈ 1 rtol = 1e-4
    @test_throws ErrorException optimize(x -> Inf, i -> nothing, 50)
end


text:
# Autogenerated wrapper script for GSL_jll for armv7l-linux-gnueabihf
export gsl_histogram, gsl_randist, libgsl, libgslcblas

JLLWrappers.@generate_wrapper_header("GSL")
JLLWrappers.@declare_executable_product(gsl_histogram)
JLLWrappers.@declare_executable_product(gsl_randist)
JLLWrappers.@declare_library_product(libgsl, "libgsl.so.25")
JLLWrappers.@declare_library_product(libgslcblas, "libgslcblas.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        gsl_histogram,
        "bin/gsl-histogram",
    )

    JLLWrappers.@init_executable_product(
        gsl_randist,
        "bin/gsl-randist",
    )

    JLLWrappers.@init_library_product(
        libgsl,
        "lib/libgsl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgslcblas,
        "lib/libgslcblas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MIToS.MSA
using MIToS.Information

const msa_long = read("../data/PF00089_aligned.fasta", FASTA)
const msa_wide = read("../data/PF16957_aligned.fasta", FASTA)

mapcolfreq!(entropy, msa_long, Probabilities{Float64,1,UngappedAlphabet}(ContingencyTable(Float64,Val{1},UngappedAlphabet())))
mapcolfreq!(entropy, msa_wide, Probabilities{Float64,1,UngappedAlphabet}(ContingencyTable(Float64,Val{1},UngappedAlphabet())))

println("[BENCH] Shannon entropy PF00089: ", @elapsed mapcolfreq!(entropy, msa_long, Probabilities{Float64,1,UngappedAlphabet}(ContingencyTable(Float64,Val{1},UngappedAlphabet()))))
println("[BENCH] Shannon entropy PF16957: ", @elapsed mapcolfreq!(entropy, msa_wide, Probabilities{Float64,1,UngappedAlphabet}(ContingencyTable(Float64,Val{1},UngappedAlphabet()))))

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
export AbstractDomain,
GenericDomain

export isWaveguide,
isBulkWaveguide,
isBackground,
isDefect,
isPC,
isPlanar,
isHalfSpace

abstract type AbstractForegroundBackgound end
struct Foreground<:AbstractForegroundBackgound end
struct Background<:AbstractForegroundBackgound end

abstract type AbstractWaveguideNot end
struct Waveguide<:AbstractWaveguideNot end
struct NotWaveguide<:AbstractWaveguideNot end

abstract type AbstractDomainType end
struct PlanarDomainType<:AbstractDomainType end
struct HalfspaceDomainType<:AbstractDomainType end
struct PCDomainType<:AbstractDomainType end
struct DefectDomainType<:AbstractDomainType end
struct GenericDomainType<:AbstractDomainType end

abstract type AbstractBulkAsymptote end
struct Bulk<:AbstractBulkAsymptote end
struct Asymptote<:AbstractBulkAsymptote end

abstract type AbstractDomain{TFB,TWN,TDT,TBA} end
struct PlanarWaveguide          <:AbstractDomain{Foreground,Waveguide   ,PlanarDomainType   ,Asymptote  } end
Base.show(io::IO,::PlanarWaveguide) = print(io,"Planar Waveguide")
struct PlanarWaveguideBackground<:AbstractDomain{Background,Waveguide   ,PlanarDomainType   ,Asymptote  } end
Base.show(io::IO,::PlanarWaveguideBackground) = print(io,"Planar Waveguide Background")
struct Halfspace                <:AbstractDomain{Foreground,Waveguide   ,HalfspaceDomainType,Asymptote  } end
Base.show(io::IO,::Halfspace) = print(io,"Halfspace")
struct PCWaveguide              <:AbstractDomain{Foreground,Waveguide   ,PCDomainType       ,Asymptote  } end
Base.show(io::IO,::PCWaveguide) = print(io,"PC Waveguides")
struct PCWaveguideBackground    <:AbstractDomain{Background,Waveguide   ,PCDomainType       ,Asymptote  } end
Base.show(io::IO,::PCWaveguideBackground) = print(io,"PC Waveguide Background")
struct BulkPlanarWaveguideX     <:AbstractDomain{Foreground,Waveguide   ,PlanarDomainType   ,Bulk       } end
Base.show(io::IO,::BulkPlanarWaveguideX) = print(io,"Bulk Planar Waveguide X")
struct BulkPlanarWaveguideY     <:AbstractDomain{Foreground,Waveguide   ,PlanarDomainType   ,Bulk       } end
Base.show(io::IO,::BulkPlanarWaveguideY) = print(io,"Bulk Planar Waveguide Y")
struct BulkPCWaveguideX         <:AbstractDomain{Foreground,Waveguide   ,PCDomainType       ,Bulk       } end
Base.show(io::IO,::BulkPCWaveguideX) = print(io,"Bulk PC Waveguide X")
struct BulkPCWaveguideY         <:AbstractDomain{Foreground,Waveguide   ,PCDomainType       ,Bulk       } end
Base.show(io::IO,::BulkPCWaveguideY) = print(io,"Bulk PC Waveguide Y")
struct BackgroundDomain         <:AbstractDomain{Background,NotWaveguide,GenericDomainType  ,Bulk       } end
Base.show(io::IO,::BackgroundDomain) = print(io,"Background")
struct DefectDomain             <:AbstractDomain{Foreground,NotWaveguide,DefectDomainType   ,Bulk       } end
Base.show(io::IO,::DefectDomain) = print(io,"Defect")
struct SiteDefect               <:AbstractDomain{Foreground,NotWaveguide,DefectDomainType   ,Bulk       } end
Base.show(io::IO,::SiteDefect) = print(io,"Site Defect")
struct LineDefect               <:AbstractDomain{Foreground,NotWaveguide,DefectDomainType   ,Bulk       } end
Base.show(io::IO,::LineDefect) = print(io,"Line Defect")
struct PCDomain                 <:AbstractDomain{Foreground,NotWaveguide,PCDomainType       ,Bulk       } end
Base.show(io::IO,::PCDomain) = print(io,"PC")
struct GenericDomain            <:AbstractDomain{Background,NotWaveguide,GenericDomainType  ,Bulk       } end
Base.show(io::IO,::GenericDomain) = print(io,"Generic")

isWaveguide(domain::T) where T<:AbstractDomain{_1,Waveguide} where _1= true
isWaveguide(domain) = false

isBulkWaveguide(domain::T) where T<:AbstractDomain{_1,Waveguide,_2,Bulk} where {_1,_2} = true
isBulkWaveguide(domain) = false

isBackground(domain::T) where T<:AbstractDomain{Background} = true
isBackground(domain) = false

isDefect(domain::T) where T<:AbstractDomain{_1,_2,DefectDomainType} where {_1,_2} = true
isDefect(domain) = false

isPC(domain::T) where T<:AbstractDomain{_1,_2,PCDomainType} where {_1,_2} = true
isPC(domain) = false

isPlanar(domain::T) where T<:AbstractDomain{_1,_2,PlanarDomainType} where {_1,_2} = true
isPlanar(domain) = false

isHalfSpace(domain::T) where T<:AbstractDomain{_1,_2,HalfspaceDomainType} where {_1,_2} = true
isHalfSpace(domain) = false


text:
######################################################################
# Voidables are wrappers around voidable values. The intention is to
# reference the Voidable rather than the wrapped value such that the
# wrapped value could be voided across all references.
# Essentially a more semantic implementation of `Ref{Union{Nothing, T}}`.
# -----
# Licensed under Apache License 2.0
module Voidables

export VoidError
struct VoidError <: Exception end

export Void, void
struct Void end
const void = Void()

export Voidable
mutable struct Voidable{T}
    value::Union{Void, T}
end
Voidable() = Voidable{Any}(void)
Voidable(::Void) = Voidable()
Voidable{T}() where T = Voidable{T}(void)

Base.convert(::Type{Voidable{T1}}, voidable::Voidable{T2}) where {T1, T2<:T1} = Voidable{T1}(voidable.value)
Base.convert(::Type{Voidable{T}}, voidable::Voidable{Any}) where T = Voidable{T}(void_convert(T, voidable.value))
void_convert(::Type, ::Void) = void
void_convert(::Type{T}, value::T) where T = value

function (Base.getindex(voidable::Voidable{T})::T) where T
    if voidable.value === void
        throw(VoidError())
    end
    voidable.value
end

Base.setindex!(voidable::Voidable{T}, value::T) where T = voidable.value = value

export isvoided
isvoided(voidable::Voidable) = voidable.value === void

export void!
void!(voidable::Voidable) = voidable.value = void

export voidable_param
(voidable_param(::Voidable{T})::Type{T}) where T = T

end # module Voidables


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

# PR #25567
Base.@deprecate_binding dSFMT DSFMT

# PR #21359

@deprecate srand(r::MersenneTwister, filename::AbstractString, n::Integer=4) srand(r, read!(filename, Vector{UInt32}(uninitialized, Int(n))))
@deprecate srand(filename::AbstractString, n::Integer=4) srand(read!(filename, Vector{UInt32}(uninitialized, Int(n))))
@deprecate MersenneTwister(filename::AbstractString)  srand(MersenneTwister(0), read!(filename, Vector{UInt32}(uninitialized, Int(4))))

function randjump(mt::MersenneTwister, jumps::Integer, jumppoly::AbstractString)
    depwarn("`randjump(rng, jumps, jumppoly::AbstractString)` is deprecated; use `randjump(rng, steps, jumps)` instead", :randjump)
    Base.Random._randjump(mt, DSFMT.GF2X(jumppoly), jumps)
end

@deprecate randjump(mt::MersenneTwister, jumps::Integer)  randjump(mt, big(10)^20, jumps)

@deprecate convert(::Type{UInt128},     u::UUID)     UInt128(u)
@deprecate convert(::Type{UUID}, s::AbstractString)  UUID(s)

# PR #25429
@deprecate rand(r::AbstractRNG, dims::Dims) rand(r, Float64, dims)
@deprecate rand(                dims::Dims) rand(Float64, dims)

# PR #25668
@deprecate RandomDevice(unlimited::Bool) RandomDevice(; unlimited=unlimited)

max_stars_repo_name: ashaywalke/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
export SingleGrid

using CausalityToolsBase
using SparseArrays

include("binvisits.jl")

"""
    SingleGrid(b::RectangularBinning, bc::String = "none", f::Real = 1.0) <: GridBasedTransferOperator

A simple estimator for approximating the transfer operator over a partition dictated by `b`.
`bc` is the boundary condition, and `f` is the allocation factor.

See also: [`RectangularBinning`](@ref).
"""
struct SingleGrid{B <: RectangularBinning} <: GridBasedTransferOperator
    b::B
    bc::String
    f::Real

    function SingleGrid(b::B, bc::String = "circular", f::Real = 1.0) where B <: RectangularBinning
        isboundarycondition(bc, "grid")  || error("Boundary condition '$bc' not valid.")
        new{B}(b, bc, f)
    end
end

Base.show(io::IO, g::SingleGrid) = print(io, "SingleGrid{$(g.b)}")

function transopergenerator(pts, method::SingleGrid)
    b = method.b

    # Calculate minima and edgelengths given the chosen grid
    mini, el = get_minima_and_edgelengths(pts, b)

    # Tie points to grid cells using integer tuples
    encoded_pts = encode(pts, mini, el)

    # Coordinates of the visited bins relative to grid.
    vb = [el .* pt .+ mini for pt in encoded_pts]

    # TODO: this is the expensive operation. Speed ups possible by
    # making versions of the groupslice functions that operate on
    # integers tuples instead of this matrix of integers?
    bv = get_binvisits(hcat(encoded_pts...,))

    init = (vb = vb, bv = bv)
    TransferOperatorGenerator(method, pts, init)
end

Base.show(io::IO, g::TransferOperatorGenerator) = print(io, "TransferOperatorGenerator{$(g.method)}")

function zerocols(M::AbstractArray{T, 2}) where T
    findall(sum(M[:, i] for i = 1:size(M, 2)) .== 0)
end

function (tog::TransferOperatorGenerator{<:SingleGrid})()
    boundary_condition = tog.method.bc

    alloc_frac = tog.method.f
    if !(0 < alloc_frac <= 1)
        error("allocation fraction needs to be on the interval (0, 1]")
    end

    bv, vb = getfield.(Ref(tog.init), (:bv, :vb))
    first_visited_by = bv.first_visited_by
    visitors = bv.visitors
    visits_whichbin = bv.visits_whichbin
    n_visited = length(first_visited_by)
    valid_boundary_conditions = [:none, :exclude, :circular, :invariantize]

    # Initialise transfer (Perron-Frobenius) operator as a sparse matrix
    # (keep track of indices and values in separate columns for now)
    n_possible_nonzero_entries = n_visited^2
    N = ceil(Int, n_possible_nonzero_entries / alloc_frac)

    I = zeros(Int, N)
    J = zeros(Int, N)
    P = zeros(N)

    # Preallocate target index for the case where there is only
    # one point of the orbit visiting a bin.
    target_bin_j::Int = 0
    n_visitsᵢ::Int = 0

    # Keep track of how many transitions we've considered.
    k = 0

    if boundary_condition == "circular"
        append!(visits_whichbin, [1])
    elseif boundary_condition == "random"
        append!(visits_whichbin, [rand(1:length(visits_whichbin))])
    end

    # Loop over the visited bins bᵢ
    for i in 1:n_visited
        # How many times is this bin visited?
        n_visitsᵢ = length(visitors[i])

        # If both conditions below are true, then there is just one
        # point visiting the i-th bin. If there is only one visiting point and
        # it happens to be the last, we skip it, because we don't know its
        # image.
        if n_visitsᵢ == 1 && !(i == visits_whichbin[end])
            k += 1
            # To which bin does the single point visiting bᵢ jump if we shift it one time step ahead along its orbit?
            target_bin_j = visits_whichbin[visitors[i][1] + 1][1]

            # We now know that exactly one point (the i-th) does the
            # transition from i to the target j.
            I[k] = i
            J[k] = target_bin_j
            P[k] = 1
        end
        # If more than one point of the orbit visits the i-th bin, we
        # identify the visiting points and track which bins bⱼ they end up
        # in after the forward linear map of the points.
        if n_visitsᵢ > 1
            timeindices_visiting_pts = visitors[i]

            # TODO: Introduce circular boundary condition. Simply excluding
            # might lead to a cascade of loosing points.

            # If bᵢ is the bin visited by the last point in the orbit, then
            # the last entry of `visiting_pts` will be the time index of the
            # last point of the orbit. In the next time step, that point will
            # have nowhere to go along its orbit (precisely because it is the
            # last data point). Thus, we exclude it.
            if i == visits_whichbin[end]
                #warn("Removing last point")
                n_visitsᵢ = length(timeindices_visiting_pts) - 1
                timeindices_visiting_pts = timeindices_visiting_pts[1:(end - 1)]
            end

            # To which boxes do each of the visitors to bᵢ jump in the next
            # time step?
            target_bins = visits_whichbin[timeindices_visiting_pts .+ 1]
            unique_target_bins = unique(target_bins)

            # Count how many points jump from the i-th bin to each of
            # the unique target bins, and use that to calculate the transition
            # probability from bᵢ to bⱼ.
            for j in 1:length(unique_target_bins)
                n_transitions_i_to_j = sum(target_bins .== unique_target_bins[j])

                k += 1
                I[k] = i
                J[k] = unique_target_bins[j]
                P[k] = n_transitions_i_to_j / n_visitsᵢ
            end
        end
    end

    # Combine indices and values into a sparse representation of the Perron-
    # Frobenius operator (transfer operator). Filter out the nonzero entries.
    # The number of rows in the matrix is given by the number of unique points
    # we have in the embedding.
    M = Array(sparse(I[1:k], J[1:k], P[1:k], n_visited, n_visited))

    # There may be boxes which are visited by points of the orbit, but not by
    # any image points.  # If the first bin is not visited, make every point starting in that bin
    # jump to the same bin with probability one.
    zc = zerocols(M)
    l = length(zc)
    if l > 0
        #@warn "There were $l all-zero columns.\nColumn(s) $(zc) is/are the culprit(s)."

        M[zc, 1] .= 1.0
    end

    params = NamedTuple()
    TransferOperatorApproximation(tog, M, params)
end

function transferoperator(pts, method::SingleGrid)
    tog = transopergenerator(pts, method)
    tog()
end


text:
using HorseML.Preprocessing
using HorseML.Classification
using HorseML.Tree
using CSV
import DataFrames

@testset "Classification" begin
    data = Matrix(dataloader("iris"))
    x, t = data[:, 1:4], data[:, 5]

    #The test for Decision Tree
    model = DecisionTree()
    @test_nowarn Tree.fit!(model, x, t)
    @test_throws DimensionMismatch Tree.fit!(model, x', t)
    @test_nowarn model(x)
    @test_nowarn MV("test.dot", model)

    #The test for Random forest
    model = RandomForest(10)
    @test_nowarn Tree.fit!(model, x, t)
    @test_nowarn model(x)
    paths = ["test$(string(i)).dot" for i in 1 : 10]
    @test_nowarn MV(paths, model)

    LE = LabelEncoder()
    OHE = OneHotEncoder()
    t = LE(t)
    t = OHE(t)

    #The test for Logistic Regression
    model = Logistic(alpha = 0.1)
    @test_nowarn Classification.fit!(model, x, t)
    @test_throws DimensionMismatch Classification.fit!(model, x', t)
    @test_nowarn model(x)

    #The test for Support Vector machine Classification(One-Vs-Rest)
    model = SVC()
    @test_nowarn Classification.fit!(model, x, t)
    @test_nowarn model(x)
end


text:
function rbase()
    return Distributions.sample([DNA_A, DNA_C, DNA_G, DNA_T])
end

function mutate_base(base::DNA)
    result = rbase()
    while result == base
        result = rbase()
    end
    return result
end

function mutate_seq(seq::DNASeq, n_diffs::Int)
    seq = copy(seq)
    positions = Base.rand(1:length(seq), n_diffs)
    for i in positions
        seq[i] = mutate_base(seq[i])
    end
    return seq
end

function random_codon()
    return (rbase(), rbase(), rbase())
end

function random_seq(n)
    return DNASeq(map(_ -> rbase(), 1:n))
end

# bounds for per-base error probabilities
const MIN_PROB = Prob(1e-10)
const MAX_PROB = Prob(0.5)

"""Add independent noise to each position in vector."""
function jitter_phred_domain(x::Vector{Prob},
                             phred_std::Prob)
    error = randn(length(x)) * phred_std / 10.0
    result = exp10.(log10.(x) + error)
    result[map(a -> a < MIN_PROB, result)] = MIN_PROB
    result[map(a -> a > MAX_PROB, result)] = MAX_PROB
    return result
end

function hmm_sample(sequence::DNASeq,
                    error_p::Vector{Prob},
                    errors::ErrorModel)
    errors = normalize(errors)
    codon = errors.codon_insertion > 0.0 || errors.codon_deletion > 0.0
    if codon && (errors.insertion > 0.0 || errors.deletion > 0.0)
        error("codon and non-codon indels are not both allowed")
    end
    sub_ratio = errors.mismatch
    ins_ratio = errors.insertion
    del_ratio = errors.deletion
    if codon
        ins_ratio = errors.codon_insertion
        del_ratio = errors.codon_deletion
    end
    final_seq = []
    final_error_p = Prob[]
    seqbools = Bool[]
    tbools = Bool[]
    skip = 0
    for i = 1:(length(sequence) + 1)
        p = (i > length(sequence) ? error_p[i - 1] : error_p[i])
        prev_p = (i == 1 ? error_p[1] : error_p[i - 1])
        # insertion between i-1 and i
        max_p = max(p, prev_p)
        ins_p = max_p * ins_ratio
        if codon
            ins_p /= 3.0
        end
        if Base.rand(Bernoulli(ins_p)) == 1
            if codon
                push!(final_seq, random_codon()...)
                push!(final_error_p, collect(Base.Iterators.repeated(max_p, 3))...)
                push!(seqbools, false, false, false)
            else
                push!(final_seq, rbase())
                push!(final_error_p, max_p)
                push!(seqbools, false)
            end
        end
        if i > length(sequence)
            break
        end

        # only skip after insertions, to ensure equal probability of
        # insertion and deletions
        if skip > 0
            skip -= 1
            continue
        end
        # deletion of i
        if codon
            if i > length(sequence) - 2
                del_p = 0.0
            else
                del_p = maximum(error_p[i:i + 2]) * del_ratio / 3.0
            end
        else
            del_p = p * del_ratio
        end
        if Base.rand(Bernoulli(del_p)) == 1
            # skip position i, and possibly the entire codon starting at i
            skip = codon ? 2 : 0
            append!(tbools, fill(false, skip + 1))
        else
            # mutation of position i
            if Base.rand(Bernoulli(p * sub_ratio)) == 1
                push!(final_seq, mutate_base(sequence[i]))
                push!(seqbools, false)
                push!(tbools, false)
            else
                push!(final_seq, sequence[i])
                push!(seqbools, true)
                push!(tbools, true)
            end
            push!(final_error_p, p)
        end
    end
    return DNASeq(final_seq), final_error_p, seqbools, tbools
end

function sample_reference(template::DNASeq,
                          error_rate::Prob,
                          errors::ErrorModel)
    errors = normalize(errors)
    if errors.insertion > 0.0 || errors.deletion > 0.0
        error("non-codon indels are not allowed in reference")
    end
    error_p = error_rate * ones(length(template))
    reference, _, _, _ = hmm_sample(template, error_p, errors)
    if length(reference) % 3 == 1
        # randomly remove a base
        idx = rand(1:length(reference))
        reference = DNASeq(reference[1:(idx-1)], reference[(idx+1):end])
    elseif length(reference) % 3 == 2
        # randomly add base
        idx = rand(1:(length(reference) + 1))
        reference = DNASeq(reference[1:(idx-1)], DNASeq([rbase()]), reference[idx:end])
    end
    return reference
end

function sample_from_template(template::DNASeq,
                              template_error_p::Vector{Prob},
                              errors::ErrorModel,
                              phred_scale::Float64,
                              actual_std::Float64,
                              reported_std::Float64)
    errors = normalize(errors)
    if errors.codon_insertion > 0.0 || errors.codon_deletion > 0.0
        error("codon indels are not allowed in sequences")
    end
    # add noise to simulate measurement error
    d = Exponential(phred_scale)
    base_vector = exp10.((-10.0 * log10.(template_error_p) + rand(d)) / (-10.0))
    jittered_error_p = jitter_phred_domain(base_vector,
                                           actual_std)

    seq, actual_error_p, sbools, tbools = hmm_sample(template,
                                                     jittered_error_p,
                                                     errors)

    # add noise to simulate quality score estimation error
    reported_error_p = jitter_phred_domain(actual_error_p,
                                           reported_std)
    phreds = p_to_phred(reported_error_p)
    return DNASeq(seq), actual_error_p, phreds, sbools, tbools
end

function sample_mixture(nseqs::Tuple{Int,Int}, len::Int, n_diffs::Int;
                        ref_error_rate::Prob=0.1,
                        ref_errors::ErrorModel=ErrorModel(10, 0, 0, 1, 0),
                        error_rate::Prob=0.01,
                        alpha::Float64=0.1,
                        phred_scale::Float64=1.5,
                        actual_std::Float64=3.0,
                        reported_std::Float64=1.0,
                        seq_errors::ErrorModel=ErrorModel(1, 5, 5))
    template1 = random_seq(len)
    template2 = mutate_seq(template1, n_diffs)
    templates = DNASeq[template1, template2]

    reference = sample_reference(template1,
                                 ref_error_rate,
                                 ref_errors)

    # generate template error rates from four-parameter Beta distribution
    beta = alpha * (error_rate - MAX_PROB) / (MIN_PROB - error_rate)
    error_dist = Beta(alpha, beta)
    template_error_p = rand(error_dist, len) * (MAX_PROB - MIN_PROB) + MIN_PROB

    seqs = DNASeq[]
    actual_error_ps = Vector{Prob}[]
    phreds = Vector{Phred}[]
    seqbools = Vector{Bool}[]
    tbools = Vector{Bool}[]
    for (t, n) in zip(templates, nseqs)
        for i = 1:n
            (seq, actual_error_p, phred, cb,
             db) = sample_from_template(t,
                                        template_error_p,
                                        seq_errors,
                                        phred_scale,
                                        actual_std,
                                        reported_std)
            push!(seqs, seq)
            push!(actual_error_ps, actual_error_p)
            push!(phreds, phred)
            push!(seqbools, cb)
            push!(tbools, db)
        end
    end
    return (DNASeq(reference),
            DNASeq[DNASeq(t) for t in templates],
            template_error_p, seqs, actual_error_ps, phreds,
            seqbools, tbools)
end

"""
    sample_sequences(nseqs, len; kwargs...)

Generate a template and sample simulated reads and Phred scores.

This function is meant for simple testing and benchmarking, and is not
meant to represent a realistic error model.

# Arguments:
- `nseqs::Int=3`: number of reads to generate

- `len::Int=90`: length of template

- `ref_error_rate::Prob=0.1`: reference error rate

- `ref_errors::ErrorModel=ErrorModel(10, 0, 0, 1, 1)`: reference error
  model

- `error_rate::Prob=0.01`: read error rate

- `alpha::Float64=0.1`: α parameter for beta distribution of
  per-base template error rates.

- `phred_scale::Float64=1.5`: λ parameter for exponential
  distribution of Phred error

- `actual_std::Float64=3.0`: σ^2 for true Gaussian errors in
  the Phred domain

- `reported_std::Float64=1.0`: σ^2 for Gaussian errors in the
  Phred domain

- `seq_errors::ErrorModel=ErrorModel(1, 5, 5)`: sequencing error model

# Returns:
- `reference::DNASeq`: reference sequence for `template`

- `template::DNASeq`: template sequence

- `t_p::Vector{Prob}`: template error probabilities

- `seqs::Vector{DNASeq}`: simulated reads

- `actual::Vector{Vector{Prob}}`: error probabilities

- `phreds::{Vector{Vector{Phred}}`: Phred values

- `seqbools::Vector{Vector{Bool}}`: `seqbools[i][j]` is `true` if
  `seqs[i][j]` was correctly sequenced from the template

- `tbools::Vector{Vector{Bool}}`: `tbools[i][j]` is `true` if
  `template[j]` was correctly sequenced in `seqs[i]`

"""
function sample_sequences(nseqs::Int=3,
                          len::Int=90;
                          ref_error_rate::Prob=0.1,
                          ref_errors::ErrorModel=ErrorModel(10, 0, 0, 1, 1),
                          error_rate::Prob=0.01,
                          alpha::Float64=0.1,
                          phred_scale::Float64=1.5,
                          actual_std::Float64=3.0,
                          reported_std::Float64=1.0,
                          seq_errors::ErrorModel=ErrorModel(1, 5, 5))
    (ref, templates, t_p, seqs, actual,
     phreds, cb, db) = sample_mixture((nseqs, 0), len, 0;
                                      ref_error_rate=ref_error_rate,
                                      ref_errors=ref_errors,
                                      error_rate=error_rate,
                                      alpha=alpha,
                                      phred_scale=phred_scale,
                                      actual_std=actual_std,
                                      reported_std=reported_std,
                                      seq_errors=seq_errors)
    return ref, templates[1], t_p, seqs, actual, phreds, cb, db
end

"""Write template into FASTA and sequences into FASTQ."""
function write_samples(filename, reference, template, template_error, seqs, phreds)
    template_phred = p_to_phred(template_error)
    write_fasta(string(filename, "-reference.fasta"), [reference])
    write_fastq(string(filename, "-template.fastq"), [template], Vector{Phred}[template_phred])
    write_fastq(string(filename, "-sequences.fastq"), seqs, phreds)
end

"""Read template from FASTA and sequences from FASTQ."""
function read_samples(filename)
    reference = read_fasta(string(filename, "-reference.fasta"))[1]
    template_seqs, template_phreds = read_fastq(string(filename, "-template.fastq"))
    template = template_seqs[1]
    template_error = phred_to_p(template_phreds[1])
    seqs, phreds = read_fastq(string(filename, "-sequences.fastq"))
    return reference, template, template_error, seqs, phreds
end


text:
# Unit tests for ab_integration.jl

@testset "ab_integration.jl" begin
    # TODO
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
function AbsolutePVCoordinates(arg0::AbsoluteDate, arg1::AbsolutePVCoordinates, arg2::AbsolutePVCoordinates)
    return AbsolutePVCoordinates((AbsoluteDate, AbsolutePVCoordinates, AbsolutePVCoordinates), arg0, arg1, arg2)
end

function AbsolutePVCoordinates(arg0::AbsoluteDate, arg1::jdouble, arg2::AbsolutePVCoordinates)
    return AbsolutePVCoordinates((AbsoluteDate, jdouble, AbsolutePVCoordinates), arg0, arg1, arg2)
end

function AbsolutePVCoordinates(arg0::AbsoluteDate, arg1::jdouble, arg2::AbsolutePVCoordinates, arg3::jdouble, arg4::AbsolutePVCoordinates)
    return AbsolutePVCoordinates((AbsoluteDate, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4)
end

function AbsolutePVCoordinates(arg0::AbsoluteDate, arg1::jdouble, arg2::AbsolutePVCoordinates, arg3::jdouble, arg4::AbsolutePVCoordinates, arg5::jdouble, arg6::AbsolutePVCoordinates)
    return AbsolutePVCoordinates((AbsoluteDate, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function AbsolutePVCoordinates(arg0::AbsoluteDate, arg1::jdouble, arg2::AbsolutePVCoordinates, arg3::jdouble, arg4::AbsolutePVCoordinates, arg5::jdouble, arg6::AbsolutePVCoordinates, arg7::jdouble, arg8::AbsolutePVCoordinates)
    return AbsolutePVCoordinates((AbsoluteDate, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates, jdouble, AbsolutePVCoordinates), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function AbsolutePVCoordinates(arg0::Frame, arg1::AbsoluteDate, arg2::FieldVector3D)
    return AbsolutePVCoordinates((Frame, AbsoluteDate, FieldVector3D), arg0, arg1, arg2)
end

function AbsolutePVCoordinates(arg0::Frame, arg1::AbsoluteDate, arg2::PVCoordinates)
    return AbsolutePVCoordinates((Frame, AbsoluteDate, PVCoordinates), arg0, arg1, arg2)
end

function AbsolutePVCoordinates(arg0::Frame, arg1::AbsoluteDate, arg2::Vector3D, arg3::Vector3D)
    return AbsolutePVCoordinates((Frame, AbsoluteDate, Vector3D, Vector3D), arg0, arg1, arg2, arg3)
end

function AbsolutePVCoordinates(arg0::Frame, arg1::AbsoluteDate, arg2::Vector3D, arg3::Vector3D, arg4::Vector3D)
    return AbsolutePVCoordinates((Frame, AbsoluteDate, Vector3D, Vector3D, Vector3D), arg0, arg1, arg2, arg3, arg4)
end

function AbsolutePVCoordinates(arg0::Frame, arg1::TimeStampedPVCoordinates)
    return AbsolutePVCoordinates((Frame, TimeStampedPVCoordinates), arg0, arg1)
end

function cross_product(::Type{PVCoordinates}, arg0::PVCoordinates, arg1::PVCoordinates)
    return jcall(PVCoordinates, "crossProduct", PVCoordinates, (PVCoordinates, PVCoordinates), arg0, arg1)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function estimate_velocity(::Type{PVCoordinates}, arg0::Vector3D, arg1::Vector3D, arg2::jdouble)
    return jcall(PVCoordinates, "estimateVelocity", Vector3D, (Vector3D, Vector3D, jdouble), arg0, arg1, arg2)
end

function get_acceleration(obj::PVCoordinates)
    return jcall(obj, "getAcceleration", Vector3D, ())
end

function get_angular_velocity(obj::PVCoordinates)
    return jcall(obj, "getAngularVelocity", Vector3D, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_date(obj::TimeStampedPVCoordinates)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_frame(obj::AbsolutePVCoordinates)
    return jcall(obj, "getFrame", Frame, ())
end

function get_momentum(obj::PVCoordinates)
    return jcall(obj, "getMomentum", Vector3D, ())
end

function get_position(obj::PVCoordinates)
    return jcall(obj, "getPosition", Vector3D, ())
end

function get_pv_coordinates(obj::AbsolutePVCoordinates)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, ())
end

function get_pv_coordinates(obj::AbsolutePVCoordinates, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::AbsolutePVCoordinates, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (Frame,), arg0)
end

function get_velocity(obj::PVCoordinates)
    return jcall(obj, "getVelocity", Vector3D, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function interpolate(::Type{AbsolutePVCoordinates}, arg0::Frame, arg1::AbsoluteDate, arg2::CartesianDerivativesFilter, arg3::Stream)
    return jcall(AbsolutePVCoordinates, "interpolate", AbsolutePVCoordinates, (Frame, AbsoluteDate, CartesianDerivativesFilter, Stream), arg0, arg1, arg2, arg3)
end

function interpolate(::Type{TimeStampedPVCoordinates}, arg0::AbsoluteDate, arg1::CartesianDerivativesFilter, arg2::Collection)
    return jcall(TimeStampedPVCoordinates, "interpolate", TimeStampedPVCoordinates, (AbsoluteDate, CartesianDerivativesFilter, Collection), arg0, arg1, arg2)
end

function interpolate(::Type{TimeStampedPVCoordinates}, arg0::AbsoluteDate, arg1::CartesianDerivativesFilter, arg2::Stream)
    return jcall(TimeStampedPVCoordinates, "interpolate", TimeStampedPVCoordinates, (AbsoluteDate, CartesianDerivativesFilter, Stream), arg0, arg1, arg2)
end

function interpolate(obj::AbsolutePVCoordinates, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", AbsolutePVCoordinates, (AbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function negate(obj::PVCoordinates)
    return jcall(obj, "negate", PVCoordinates, ())
end

function normalize(obj::PVCoordinates)
    return jcall(obj, "normalize", PVCoordinates, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function shifted_by(obj::AbsolutePVCoordinates, arg0::jdouble)
    return jcall(obj, "shiftedBy", AbsolutePVCoordinates, (jdouble,), arg0)
end

function to_derivative_structure_pv(obj::PVCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructurePV", FieldPVCoordinates, (jint,), arg0)
end

function to_derivative_structure_vector(obj::PVCoordinates, arg0::jint)
    return jcall(obj, "toDerivativeStructureVector", FieldVector3D, (jint,), arg0)
end

function to_string(obj::TimeStampedPVCoordinates)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::TimeStampedPVCoordinates, arg0::TimeScale)
    return jcall(obj, "toString", JString, (TimeScale,), arg0)
end

function to_taylor_provider(obj::AbsolutePVCoordinates)
    return jcall(obj, "toTaylorProvider", PVCoordinatesProvider, ())
end

function to_taylor_provider(obj::TimeStampedPVCoordinates, arg0::Frame)
    return jcall(obj, "toTaylorProvider", PVCoordinatesProvider, (Frame,), arg0)
end

function to_univariate_derivative1_vector(obj::PVCoordinates)
    return jcall(obj, "toUnivariateDerivative1Vector", FieldVector3D, ())
end

function to_univariate_derivative1pv(obj::PVCoordinates)
    return jcall(obj, "toUnivariateDerivative1PV", FieldPVCoordinates, ())
end

function to_univariate_derivative2_vector(obj::PVCoordinates)
    return jcall(obj, "toUnivariateDerivative2Vector", FieldVector3D, ())
end

function to_univariate_derivative2pv(obj::PVCoordinates)
    return jcall(obj, "toUnivariateDerivative2PV", FieldPVCoordinates, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
using InverseIsing
using Test

const tests =[
    "tests_inverseising",
    "tests_anneal",
    "tests_utils"
]

for t in tests
    @testset "Test $t" begin
        include("$t.jl")
    end
end

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
module YaoBlocksQASM

export convert_to_qasm

using YaoBlocks

include("qasm.jl")

end

max_stars_repo_name: Sov-trotter/YaoBlocksQASM

max_issues_repo_name: Sov-trotter/YaoBlocksQASM

max_forks_repo_name: Sov-trotter/YaoBlocksQASM


text:
function betweenness(g::AbstractStaticGraph)
    A = sparse_adjacency_matrix(g)
    expA = expm(A)
    ns = nodes(g)
    for v in ns
        i = node_index(v)
        row = A[i,:]
        column = A[:, i]
        A[i,:] = 0
        A[:,i] = 0
        B = (expA - expm(A)) / expA
    end
end

function betweenness(g::AbstractEvolvingGraph)

end


text:
# Extend MathOptInterface Test for Knitro.
#
# This file is largely adapted from :
# MathOptInterface/src/Tests/contconic.jl
# https://github.com/JuliaOpt/MathOptInterface.jl/blob/master/src/Test/contconic.jl
#
# LICENSE:
# MIT "expat" license
# Copyright (c) 2017: Miles Lubin and contributors Copyright (c) 2017: Google Inc.

using MathOptInterface

const MOI = MathOptInterface

const CONIC_OPTIMIZER = KNITRO.Optimizer(outlev=0, presolve=0, opttol=1e-8)
const CONIC_BRIDGED = MOIB.full_bridge_optimizer(CONIC_OPTIMIZER, Float64)

@testset "MOI SOCP 1" begin
    # Derived from MOI's problem SOC1
    model = CONIC_BRIDGED
    atol = 1e-4
    rtol = 1e-4

    # Problem SOC1
    # max 0x + 1y + 1z
    #  st  x            == 1
    #      x >= ||(y,z)||

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Zeros)
    @test MOI.supports_constraint(model, MOI.VectorOfVariables,MOI.SecondOrderCone)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x,y,z = MOI.add_variables(model, 3)

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0,1.0], [y,z]), 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    ceq = MOI.add_constraint(model,
                             MOI.VectorAffineFunction([MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))], [-1.0]), MOI.Zeros(1))
    vov = MOI.VectorOfVariables([x,y,z])
    csoc = MOI.add_constraint(model,
                              MOI.VectorAffineFunction{Float64}(vov), MOI.SecondOrderCone(3))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.Zeros}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.SecondOrderCone}()) == 1
    #= loc = MOI.get(model, MOI.ListOfConstraints()) =#
    #= @test length(loc) == 2 =#
    #= @test (MOI.VectorAffineFunction{Float64},MOI.Zeros) in loc =#
    #= @test (MOI.VectorAffineFunction{Float64},MOI.SecondOrderCone) in loc =#

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED

    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ √2 atol=atol rtol=rtol

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ 1 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1/√2 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), z) ≈ 1/√2 atol=atol rtol=rtol

    @test MOI.get(model, MOI.ConstraintPrimal(), ceq) ≈ [0.] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), csoc) ≈ [1., 1/√2, 1/√2] atol=atol rtol=rtol

    @test MOI.get(model, MOI.ConstraintDual(), ceq) ≈ [-√2] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintDual(), csoc) ≈ [√2, -1.0, -1.0] atol=atol rtol=rtol
end

MOI.empty!(CONIC_BRIDGED)

@testset "MOI SOCP 2" begin
    model = CONIC_BRIDGED
    atol = 1e-4
    rtol = 1e-4
    # Problem SOC2
    # min  x
    # s.t. y ≥ 1/√2
    #      x² + y² ≤ 1
    # in conic form:
    # min  x
    # s.t.  -1/√2 + y ∈ R₊
    #        1 - t ∈ {0}
    #      (t,x,y) ∈ SOC₃

    @test MOIU.supports_default_copy_to(model, #=copy_names=# false)
    @test MOI.supports(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
    @test MOI.supports(model, MOI.ObjectiveSense())
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Zeros)
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Nonpositives)
    @test MOI.supports_constraint(model, MOI.VectorAffineFunction{Float64},MOI.SecondOrderCone)

    MOI.empty!(model)
    @test MOI.is_empty(model)

    x,y,t = MOI.add_variables(model, 3)

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    cnon = MOI.add_constraint(model, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, y))], [1/√2]), MOI.Nonpositives(1))
    ceq = MOI.add_constraint(model, MOI.VectorAffineFunction([MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(-1.0, t))], [1.0]), MOI.Zeros(1))
    csoc = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([1,2,3], MOI.ScalarAffineTerm.(1.0, [t,x,y])), zeros(3)), MOI.SecondOrderCone(3))

    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64}, MOI.Nonpositives}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.Zeros}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64},MOI.SecondOrderCone}()) == 1

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED

    MOI.optimize!(model)

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED

    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ -1/√2 atol=atol rtol=rtol

    @test MOI.get(model, MOI.VariablePrimal(), x) ≈ -1/√2 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), y) ≈ 1/√2 atol=atol rtol=rtol
    @test MOI.get(model, MOI.VariablePrimal(), t) ≈ 1 atol=atol rtol=rtol

    # Getter of Nonnegatives constraints is currently broken!!
    @test_broken MOI.get(model, MOI.ConstraintPrimal(), cnon) ≈ [0.0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), ceq) ≈ [0.0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintPrimal(), csoc) ≈ [1., -1/√2, 1/√2] atol=atol rtol=rtol

    @test MOI.get(model, MOI.ConstraintDual(), cnon) ≈ [-1.0] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintDual(), ceq) ≈ [√2] atol=atol rtol=rtol
    @test MOI.get(model, MOI.ConstraintDual(), csoc) ≈ [√2, 1.0, -1.0] atol=atol rtol=rtol
end


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute()::Int
    result = 0
    for i ∈ 1:9
        j, n = 1, BigInteger(i)
        while length(n) == j
            result += 1
            j += 1
            n *= i
        end
    end
    return result
end

compute()

@benchmark compute()

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
# ---
# title: 575. Distribute Candies
# id: problem575
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table
# link: <https://leetcode.com/problems/distribute-candies/description/>
# hidden: true
# ---
#
# Alice has `n` candies, where the `ith` candy is of type `candyType[i]`. Alice
# noticed that she started to gain weight, so she visited a doctor.
#
# The doctor advised Alice to only eat `n / 2` of the candies she has (`n` is
# always even). Alice likes her candies very much, and she wants to eat the
# maximum number of different types of candies while still following the
# doctor's advice.
#
# Given the integer array `candyType` of length `n`, return _the **maximum**
# number of different types of candies she can eat if she only eats _`n / 2` _of
# them_.
#
#
#
# **Example 1:**
#
#
#
#     Input: candyType = [1,1,2,2,3,3]
#     Output: 3
#     Explanation: Alice can only eat 6 / 2 = 3 candies. Since there are only 3 types, she can eat one of each type.
#
#
# **Example 2:**
#
#
#
#     Input: candyType = [1,1,2,3]
#     Output: 2
#     Explanation: Alice can only eat 4 / 2 = 2 candies. Whether she eats types [1,2], [1,3], or [2,3], she still can only eat 2 different types.
#
#
# **Example 3:**
#
#
#
#     Input: candyType = [6,6,6,6]
#     Output: 1
#     Explanation: Alice can only eat 4 / 2 = 2 candies. Even though she can eat 2 candies, she only has 1 type.
#
#
#
#
# **Constraints:**
#
#   * `n == candyType.length`
#   * `2 <= n <= 104`
#   * `n` is even.
#   * `-105 <= candyType[i] <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#-------------------------------------------------------------------------------
# abstract classes for interpolators

export BasicInterpolator
export OneDimensionalInterpolator
export TwoDimensionalInterpolator

abstract type BasicInterpolator end
abstract type OneDimensionalInterpolator <: BasicInterpolator end
abstract type TwoDimensionalInterpolator <: BasicInterpolator end

function Base.show(io::IO, ϕ::P) where {P<:OneDimensionalInterpolator}
    print(io, "$P ∈ [$(ϕ.r.xa), $(ϕ.r.xb)]")
end

function Base.show(io::IO, ϕ::P) where {P<:TwoDimensionalInterpolator}
    print(io, "$P ∈ [$(ϕ.G.xa), $(ϕ.G.xb)], [$(ϕ.G.ya), $(ϕ.G.yb)]")
end

#-------------------------------------------------------------------------------
# cell finding stuff

export findcell

"""
    findcell(q, V, n)

Use bisection search to find the cell containing `q`, assuming `V` is a sorted vector of `n` coordinates. The returned integer is the index of the element in `V` immediately less than `q`. For example, if `findcell` returns 2, then `q ∈ [V[2],V[3])`. If `q` is less than every element in `V`, 1 is returned, indicating the first cell in `V`. If `q` is greater than every element in `V`, `length(V)-1` is returned, indicating the last cell in `V`.
"""
function findcell(q, V, n::Int64)::Int64
    #handle boundaries
    @inbounds (q <= V[1]) && return(1)
    @inbounds (q >= V[n]) && return(n - 1)
    #bisection search for the containing cell
    L = 0
    H = n
    while H - L > 1
        M = (H + L) ÷ 2
        @inbounds if V[M] > q
            H = M
        else
            L = M
        end
    end
    return L
end

incell(x, X, i::Int64)::Bool = @inbounds ((x >= X[i]) && (x <= X[i+1]))

function findcell(x, ϕ::OneDimensionalInterpolator)::Int64
    #check previous index used
    i::Int64 = ϕ.i[]
    xᵣ = ϕ.r.x
    incell(x, xᵣ, i) && return(i)
    #check adjacent cells
    n = ϕ.r.n
    i > 1 && incell(x, xᵣ, i-1) && return(i-1)
    i < n && incell(x, xᵣ, i+1) && return(i+1)
    #didn't work, find a fresh cell
    i = findcell(x, xᵣ, n)
    #store the index
    ϕ.i[] = i
    return i
end

function findcell(x, y, Φ::TwoDimensionalInterpolator)::NTuple{2,Int64}
    #check previous indices used
    i::Int64 = Φ.i[]
    j::Int64 = Φ.j[]
    xᵣ = Φ.G.x
    yᵣ = Φ.G.y
    incell(x, xᵣ, i) && incell(y, yᵣ, j) && return(i,j)
    #didn't work, find a fresh cell
    i = findcell(x, xᵣ, Φ.G.nx)
    j = findcell(y, yᵣ, Φ.G.ny)
    #store the indices
    Φ.i[] = i
    Φ.j[] = j
    return i, j
end

#-------------------------------------------------------------------------------
# types for handling boundaries

export AbstractBoundaries, NoBoundaries, WeakBoundaries, StrictBoundaries

upperbounderror(x, xb, axis::String) = error("Interpolation location $x outside upper interpolation limit $xb on the $axis axis")

lowerbounderror(x, xa, axis::String) = error("Interpolation location $x outside lower interpolation limit $xa on the $axis axis")

abstract type AbstractBoundaries end

#--------------------------------

"""
    NoBoundaries()

Performs no boundary checking. Allows interpolators to blindly extrapolate when possible.
"""
struct NoBoundaries <: AbstractBoundaries end

#totally empty function
function (B::NoBoundaries)(x...) end

#--------------------------------

"""
    WeakBoundaries()

Allows small overshoots at boundaries, but not large ones. Errors are only triggered when the interpolation coordinate is outside of a boundary and not close to it. At the lower boundar, for example, an error would be triggered when `(x < boundary) & !(x ≈ boundary)`.
"""
struct WeakBoundaries <: AbstractBoundaries end

function (B::WeakBoundaries)(x, xa, xb)
    ((x < xa) && !(x ≈ xa)) && lowerbounderror(x, xa, "first")
    ((x > xb) && !(x ≈ xb)) && upperbounderror(x, xb, "first")
end

function (B::WeakBoundaries)(x, xa, xb, y, ya, yb)
    ((x < xa) && !(x ≈ xa)) && lowerbounderror(x, xa, "first")
    ((x > xb) && !(x ≈ xb)) && upperbounderror(x, xb, "first")
    ((y < ya) && !(y ≈ ya)) && lowerbounderror(y, ya, "second")
    ((y > yb) && !(y ≈ yb)) && upperbounderror(y, yb, "second")
end

#--------------------------------

"""
    StrictBoundaries()

Triggers errors whenever interpolation coordinates are outside of the boundaries.
"""
struct StrictBoundaries <: AbstractBoundaries end

function (B::StrictBoundaries)(x, xa, xb)
    (x < xa) && lowerbounderror(x, xa, "first")
    (x > xb) && upperbounderror(x, xb, "first")
end

function (B::StrictBoundaries)(x, xa, xb, y, ya, yb)
    (x < xa) && lowerbounderror(x, xa, "first")
    (x > xb) && upperbounderror(x, xb, "first")
    (y < ya) && lowerbounderror(y, ya, "second")
    (y > yb) && upperbounderror(y, yb, "second")
end

#-------------------------------------------------------------------------------
# a base structure for 1d interpolations and some related functions

struct InterpolatorRange{T}
    #number of points
    n::Int64
    #grid/sample points
    x::Vector{T}
    #lowest sample value
    xa::T
    #highest sample value
    xb::T
    #values at sample points
    y::Vector{T}
end

function InterpolatorRange(x::AbstractVector, y::AbstractVector)
    #collect and copy
    T = promote_type(eltype(x), eltype(y))
    x = collect(T, x)
    y = collect(T, y)
    #check for basic problems
    rangecheck(x, y)
    #grid properties
    n = length(x)
    xa = minimum(x)
    xb = maximum(x)
    #construct
    InterpolatorRange{T}(n, x, xa, xb, y)
end

function rangecheck(x::AbstractVector, y::AbstractVector, minpts=2)
    nx, ny = length(x), length(y)
    @assert nx == ny "length of x ($nx) does not match length of y ($ny)"
    @assert nx >= minpts "must have at least $minpts points in the interpolation range"
    @assert all(diff(x) .> 0.0) "grid points must be in strictly ascending order"
end

function linstruct(I::Type, f::F, xa, xb, n::Int, boundaries::AbstractBoundaries) where {F}
    @assert n > 1 "cannot construct a length<=1 interpolation range"
    x = LinRange(xa, xb, n)
    y = f.(x)
    I(x, y, boundaries)
end

#-------------------------------------------------------------------------------
# a base structure for 2d grid interpolation and some related functions

struct InterpolatorGrid{T}
    #number of points along axis 1
    nx::Int64
    #number of points along axis 2
    ny::Int64
    #grid points along axis 1
    x::Vector{T}
    #grid points along axis 2
    y::Vector{T}
    #lowest value on axis 1
    xa::T
    #highest value on axis 1
    xb::T
    #lowest value on axis 2
    ya::T
    #highest value on axis 2
    yb::T
    #values at grid points
    Z::Array{T,2}
end

function InterpolatorGrid(x::AbstractVector, y::AbstractVector, Z::AbstractMatrix)
    #collect and copy
    T = promote_type(eltype(x), eltype(y), eltype(Z))
    x = collect(T, x)
    y = collect(T, y)
    Z = collect(T, Z)
    #check for basic grid problems
    gridcheck(x, y, Z)
    #boundaries
    xa, xb = minimum(x), maximum(x)
    ya, yb = minimum(y), maximum(y)
    #construct the object
    InterpolatorGrid{T}(length(x), length(y), x, y, xa, xb, ya, yb, Z)
end

function gridcheck(x::AbstractVector, y::AbstractVector, Z::AbstractMatrix, minpts=2)
    @assert (length(x) >= minpts) & (length(y) >= minpts) "must have at least $minpts points in each dimension of interpolation grid"
    @assert size(Z) == (length(x), length(y)) "dimensions mismatched, size(Z) = $(size(Z)) but length(x) = $(length(x)) and length(y) = $(length(y))"
    @assert (all(diff(x) .> 0.0) & all(diff(y) .>= 0.0)) "grid coordinates must be monotonically increasing without duplicates"
end

function linstruct(I::Type, f::F,
                   xa, xb, nx::Int,
                   ya, yb, ny::Int,
                   boundaries::AbstractBoundaries) where {F}
    @assert (nx > 1) & (ny > 1) "cannot use a length<=1 grid axis"
    xa, xb, ya, yb = promote(xa, xb, ya, yb)
    T = typeof(xa)
    x = LinRange(xa, xb, nx)
    y = LinRange(ya, yb, ny)
    X = x .* ones(T, ny)'
    Y = y' .* ones(T, nx)
    Z = f.(X, Y)
    I(x, y, Z, boundaries)
end


text:
cs_dict=Dict('A'=>:green,'C'=>:blue, 'G'=>:yellow, 'T'=>:red)

function print_emitters(state::Categorical) #print informative symbols from the state's emission distribution, if any
    order=Integer(log(4,length(state.p))-1)
    alphabet=CompoundAlphabet(ACGT,order)
    bits=log(2,length(state.p)) #higher order hmms express more information per symbol
    dummy_p=state.p.+10^-99 #prevent -Inf or NaN from zero probability symbols
    infoscore=(bits+sum(x*log(2,x) for x in dummy_p))
    infovec=[x*infoscore for x in dummy_p]
    infovec./=bits
    print("<<< ")
    for (symbol,symbol_prob) in enumerate(infovec)
        if symbol_prob >= .05
            str=string(alphabet.integers[symbol])
            if symbol_prob >= .7
                for char in str
                    printstyled(stdout, char; bold=true, color=cs_dict[char])
                end
            elseif .7 > symbol_prob >= .25
                for char in str
                    printstyled(stdout, char; color=cs_dict[char])
                end
            elseif .25 > symbol_prob
                for char in str
                    printstyled(stdout, lowercase(char); color=cs_dict[char])
                end
            end
            print(" ")
        end
    end
    println(">>>")
end


text:
# code for use with the MCTS.jl package
import MCTS

MCTS.node_tag(a::OPCSPAction) = "[$(a.next)]"

MCTS.node_tag(s::OPCSPBelief) = "[$(s.i)]"
MCTS.tooltip_tag(s::OPCSPBelief) = "[$(s.i), $(s.remaining)] $(s.dist.mean) | $(s.open)]"

type HeuristicActionGenerator <: MCTS.ActionGenerator
    solver::OPSolver
    rng::AbstractRNG
    _action_space::AbstractSpace
end
HeuristicActionGenerator(solver::OPSolver, rng::AbstractRNG) = HeuristicActionGenerator(solver, rng, OPCSPActionSpace())

function MCTS.next_action(gen::HeuristicActionGenerator, mdp::OPCSPBeliefMDP, s::OPCSPBelief, snode::MCTS.DPWStateNode)
    tail = tail_problem(mdp, s, prohibit_first=filter(i->i!=mdp.stop, map(a->a.next, keys(snode.A))))
    move = OPCSPAction(first_move(gen.solver, tail))
    if haskey(snode.A, move) # the best move was to the stop node, so just select something randomly
        move = rand(gen.rng, actions(mdp, s, gen._action_space))
    end
    return move
end

type PreSolvedActionGenerator <: MCTS.ActionGenerator
    presolved::Dict{Int,Int}
    secondary_solver::OPSolver
    rng::AbstractRNG
    _action_space
end
function PreSolvedActionGenerator(op::Union{OPCSPBeliefMDP,OPCSP};
                                  rng::AbstractRNG=MersenneTwister(),
                                  presolver::OPSolver=GurobiExactSolver(),
                                  secondary_solver::OPSolver=HeuristicSolver())
    d = Dict{Int,Int}()
    path = solve_op(presolver, op)
    for i in 1:length(path)-1
        d[path[i]] = path[i+1]
    end
    return PreSolvedActionGenerator(d, secondary_solver, rng, actions(op))
end

function MCTS.next_action(gen::PreSolvedActionGenerator, mdp::OPCSPBeliefMDP, s::OPCSPBelief, snode::MCTS.DPWStateNode)
    fromdict = false
    if haskey(gen.presolved, s.i)
        next = gen.presolved[s.i]
        if next in s.open && within_range(mdp, s.i, mdp.stop, next, s.remaining)
            move = OPCSPAction(gen.presolved[s.i])
            fromdict = true
        end
    end
    if !fromdict
        tail = tail_problem(mdp, s, prohibit_first=filter(i->i!=mdp.stop, map(a->a.next, keys(snode.A))))
        move = OPCSPAction(first_move(gen.secondary_solver, tail))
    end
    if haskey(snode.A, move) # the best move was to the stop node, so just select something randomly
        move = rand(gen.rng, actions(mdp, s, gen._action_space))
    end
    return move
end

# Aggregators
type OldVoronoiOPCSPAg <: MCTS.Aggregator
    radius::Float64
    anchors::Dict{Tuple{Int,Float64,IntSet},Vector{OPCSPBelief}}
end
OldVoronoiOPCSPAg(radius::Float64) = OldVoronoiOPCSPAg(radius, Dict{Tuple{Int,Float64,IntSet},Vector{OPCSPBelief}}())

function MCTS.initialize!(ag::OldVoronoiOPCSPAg, mdp::OPCSPBeliefMDP)
    ag.anchors = Dict{Tuple{Int, Float64, IntSet},Vector{OPCSPBelief}}()
end

function MCTS.assign(ag::OldVoronoiOPCSPAg, b::OPCSPBelief)
    found = false
    len = length(b.dist.mean)
    # @assert length(b.dist.mean) < 32
    # @assert b.open.limit == 256
    # @assert b.open.fill1s == false
    if haskey(ag.anchors, (b.i, b.remaining, b.open))
        similar = ag.anchors[(b.i, b.remaining, b.open)]
        local anchor
        bestdist = ag.radius
        for anch in similar
            dist = 0.0
            for i in 1:len
                @inbounds dist += abs(anch.dist.mean[i]-b.dist.mean[i])
                if dist > bestdist
                    break
                end
            end
            if dist <= bestdist
                found = true
                anchor = anch
                bestdist = dist
            end
        end
    else
        similar = ag.anchors[(b.i, b.remaining, b.open)] = Vector{OPCSPBelief}()
    end
    if !found
        push!(similar, deepcopy(b))
        return b
    end
    return anchor
end

type VoronoiOPCSPAg <: MCTS.Aggregator
    radius::Float64
    map::NearestStateMap{OPCSPState}
end
VoronoiOPCSPAg(radius::Float64) = VoronoiOPCSPAg(radius, NearestStateMap{OPCSPState}())

function MCTS.initialize!(ag::VoronoiOPCSPAg, mdp::OPCSPBeliefMDP)
    ag.map = NearestStateMap(OPCSPState, mdp)
end

function MCTS.assign(ag::VoronoiOPCSPAg, b::OPCSPBelief)
    if has_similar(ag.map, b)
        anchor, dist = get_nearest(ag.map, b)
        if dist <= ag.radius
            return anchor
        end
    end
    insert!(ag.map, b, mean(b))
    return b
end

type VoxelOPCSPAg <: MCTS.Aggregator
    diam::Float64
    # agstates::Dict{Tuple{Int,IntSet,Vector{Int}},Int}
    # next_id::Int
end
# VoxelOPCSPAg(diam::Float64) = VoxelOPCSPAg(diam, Dict{Tuple{Int,IntSet,Vector{Int}},Int}(), 1)

function MCTS.initialize!(ag::VoxelOPCSPAg, mdp::OPCSPBeliefMDP)
    # ag.agstates = Dict{Tuple{Int,IntSet,Vector{Int}},Int}()
    # ag.next_id = 1
end

function MCTS.assign(ag::VoxelOPCSPAg, b::OPCSPBelief)
    voxel = floor(Int, b.dist.mean/ag.diam)
    return (b.i, b.open, b.remaining, voxel)
end