text:
module MiscellaneousBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat
using Compat.Dates

const SUITE = BenchmarkGroup()

###########################################################################
# Splatting penalties (issue #13359)

g = addgroup!(SUITE, "splatting", ["array", "getindex"])

@noinline function perf_splatting_(A, xs...)
    A[xs...]
end
function perf_splatting(A, n, xs...)
    s = zero(eltype(A))
    for i = 1:n
        s += perf_splatting_(A, xs...)
    end
    return s
end

g[(3,3,3)] = @benchmarkable perf_splatting($(samerand(3,3,3)), 100, 1, 2, 3)

###########################################################################
# crossover from x + y + ... to afoldl (issue #13724)

function perf_afoldl(n, k)
    s = zero(eltype(k))
    @assert length(k) >= 12
    for i = 1:n
        s += k[1] + k[2] + k[3] + k[4] + k[5] + 2 * k[6] + k[7] + k[8] + k[9] + k[10] + k[11] + k[12]
    end
    return s
end

g = addgroup!(SUITE, "afoldl", ["+", "getindex"])
g["Int"] = @benchmarkable perf_afoldl(100, $(zeros(Int, 20)))
g["Float64"] = @benchmarkable perf_afoldl(100, $(zeros(Float64, 20)))
g["Complex{Float64}"] = @benchmarkable perf_afoldl(100, $(zeros(Complex{Float64}, 20)))

###########################################################################
# repeat function (issue #15553)

g = addgroup!(SUITE, "repeat", ["array"])
g[200, 24, 1] = @benchmarkable repeat($(collect(1:200)), inner=$[24], outer=$[1])
g[200, 1, 24] = @benchmarkable repeat($(collect(1:200)), inner=$[1], outer=$[24])

###########################################################################
# bitshift operators (from #18135)

function perf_bitshift(r, n)
    s = zero(eltype(r))
    for i in r
        s += i<<n
    end
    return s
end

g = addgroup!(SUITE, "bitshift", ["range"])
g["Int", "Int"] = @benchmarkable perf_bitshift($(1:1000), Int(3))
g["Int", "UInt"] = @benchmarkable perf_bitshift($(1:1000), UInt(3))
g["UInt", "UInt"] = @benchmarkable perf_bitshift($(UInt(1):UInt(1000)), UInt(3))
g["UInt32", "UInt32"] = @benchmarkable perf_bitshift($(UInt32(1):UInt32(1000)), UInt32(3))

###########################################################################
# Integer, Float64, and Date (#18000) parsing

if !hasmethod(parse, Tuple{Type{DateTime}, AbstractString})
    Base.parse(::Type{DateTime}, s::AbstractString) = DateTime(s)
end

function perf_parse(result::AbstractVector{T}, strings::AbstractVector) where T
    @assert length(result) == length(strings)
    for i = 1:length(strings)
        @inbounds result[i] = parse(T, strings[i])
    end
    return result
end

g = addgroup!(SUITE, "parse", ["DateTime"])
if VERSION <= v"0.7.0-DEV.3986"
    datestr = map(string, range(DateTime("2016-02-19T12:34:56"),Dates.Millisecond(123),200))
else
    datestr = map(string, range(DateTime("2016-02-19T12:34:56"), step = Dates.Millisecond(123), length = 200))
end
g["DateTime"] = @benchmarkable perf_parse($(similar(datestr, DateTime)), $datestr)
g["Int"] = @benchmarkable perf_parse($(Vector{Int}(undef, 1000)), $(map(string, 1:1000)))
g["Float64"] = @benchmarkable perf_parse($(Vector{Float64}(undef, 1000)), $(map(string, 1:1000)))

###########################################################################
# Julia language components (parser, etc.)

# horner-like nested expression with n levels: 1*(x + 2*(x + 2*(x + 3* ...
# ... written as a string so that we can also benchmark parsing of this function.
nestedexpr_str = """
function nestedexpr(n)
    ex = :x
    for i = n:-1:1
        ex = :(\$i * (x + \$ex))
    end
    return ex
end"""
include_string(@__MODULE__, nestedexpr_str)

if VERSION >= v"0.7.0-DEV.2437"
    const _parse = Meta.parse
else
    const _parse = Base.parse
end

g = addgroup!(SUITE, "julia")
g["parse", "array"] = @benchmarkable _parse($("[" * "a + b, "^100 * "]"))
g["parse", "nested"] = @benchmarkable _parse($(string(nestedexpr(100))))
g["parse", "function"] = @benchmarkable _parse($nestedexpr_str)
g["macroexpand", "evalpoly"] = @benchmarkable macroexpand(@__MODULE__, $(Expr(:macrocall, Symbol("@evalpoly"), 1:10...)))

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

# Issue #12165

struct FloatingPointDatatype
    class::UInt8
    bitfield1::UInt8
    bitfield2::UInt8
    bitfield3::UInt8
    size::UInt32
    bitoffset::UInt16
    bitprecision::UInt16
    exponentlocation::UInt8
    exponentsize::UInt8
    mantissalocation::UInt8
    mantissasize::UInt8
    exponentbias::UInt32
end

h5type(::Type{Float16}) =
    FloatingPointDatatype(0x00, 0x20, 0x0f, 0x00, UInt32(2), 0x0000, UInt16(16), UInt8(10), 0x05, 0x00, UInt32(10), 0x0000000f)
h5type(::Type{Float32}) =
    FloatingPointDatatype(0x00, 0x20, 0x1f, 0x00, UInt32(4), 0x0000, UInt16(32), UInt8(23), 0x08, 0x00, UInt32(23), 0x0000007f)
h5type(::Type{Float64}) =
    FloatingPointDatatype(0x00, 0x20, 0x3f, 0x00, UInt32(8), 0x0000, UInt16(64), UInt8(52), 0x0b, 0x00, UInt32(52), 0x000003ff)

struct UnsupportedFeatureException <: Exception end

function jltype(dt::FloatingPointDatatype)
    if dt == h5type(Float64)
        return 64
    elseif dt == h5type(Float32)
        return 32
    elseif dt == h5type(Float16)
        return 16
    else
        throw(UnsupportedFeatureException())
    end
end

x_16 = fill(h5type(Float16), 1000000)
x_32 = fill(h5type(Float32), 1000000)
x_64 = fill(h5type(Float64), 1000000)

function perf_jltype(x)
    y = 0
    for i = 1:length(x)
        y += jltype(x[i])
    end
    y
end

g = addgroup!(SUITE, "issue 12165")
g["Float16"] = @benchmarkable perf_jltype($x_16)
g["Float32"] = @benchmarkable perf_jltype($x_32)
g["Float64"] = @benchmarkable perf_jltype($x_64)

#########################################################################
# issue #18129

function perf_cheapest_insertion_18129(distmat::Matrix{T}, initpath::Vector{Int}) where {T<:Real}
    check_square(distmat, "Distance matrix passed to cheapest_insertion must be square.")

    n = size(distmat, 1)
    path = copy(initpath)

    # collect cities to visited
    visitus = setdiff(collect(1:n), initpath)

    # helper for insertion cost
    # tour cost change for inserting node k after the node at index after in the path
    function inscost(k, after)
        return distmat[path[after], k] +
              distmat[k, path[after + 1]] -
              distmat[path[after], path[after + 1]]
    end

    counter = 0
    while !isempty(visitus)
        bestCost = Inf
        bestInsertion = (-1, -1)
        for k in visitus
            for after in 1:(length(path) - 1) # can't insert after end of path
                counter += 1
                c = inscost(k, after)
                if c < bestCost
                    bestCost = c
                    bestInsertion = (k, after)
                end
            end
        end
        # bestInsertion now holds (k, after)
        # insert into path, remove from to-do list
        k, after = bestInsertion
        insert!(path, after + 1, k)
        visitus = setdiff(visitus, k)
    end

    return (path, pathcost(distmat, path))
end

###
# helpers
###

# make sure a passed distance matrix is a square
function check_square(m, msg)
    if size(m, 1) != size(m, 2)
        error(msg)
    end
end

# helper for readable one-line path costs
# optionally specify the bounds for the subpath we want the cost of
# defaults to the whole path
# but when calculating reversed path costs can help to have subpath costs
function pathcost(distmat::Matrix{T}, path::Vector{Int}, lb::Int = 1, ub::Int = length(path)) where {T<:Real}
    cost = zero(T)
    for i in lb:(ub - 1)
        @inbounds cost += distmat[path[i], path[i+1]]
    end
    return cost
end

dm = samerand(Float64, 300, 300)
SUITE["18129"] = @benchmarkable perf_cheapest_insertion_18129($dm, $([1, 1]))

###############################################################################
# issue #20517

function perf_dsum_20517(A::Matrix)
    z = zero(A[1,1])
    n = size(A,1)
    B = Vector{typeof(z)}(undef, n)

    @inbounds for j in 1:n
        @static if VERSION < v"0.7.0-beta.81"
            B[j] = mapreduce(k -> A[j,k]*A[k,j], +, z, 1:j)
        else
            B[j] = mapreduce(k -> A[j,k]*A[k,j], +, 1:j; init=z)
        end
    end
    B
end

A = samerand(127,127)
SUITE["20517"] = @benchmarkable perf_dsum_20517($A)

###############################################
# issue # 23042

struct Foo_23042{T<:Number, A<:AbstractMatrix{T}}
    data::A
end

Foo_23042(data::AbstractMatrix) = Foo_23042{eltype(data), typeof(data)}(data)

function perf_copy_23042(a, b)
    for i in 1:length(a.data)
        @inbounds a.data[i] = b.data[i]
    end
    a
end

g = addgroup!(SUITE, "23042")

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    b = samerand(T, 128, 128)
    a = similar(b)
    g[string(T)] = @benchmarkable perf_copy_23042($(Foo_23042(a)), $(Foo_23042(b)))
end

###############################################
# zip iterator

g = addgroup!(SUITE, "iterators", ["zip"])
for N in (1,1000), M in 1:4
    X = zip(Iterators.repeated(1:N, M)...)
    g["zip($(join(fill("1:$N", M), ", ")))"] = @benchmarkable collect($X)
end

####################################################
# Allocation elision stumped by conditional #28226 #
# Note, not fixed when this benchmark was written  #
####################################################

function perf_colwise_alloc!(r, a, b)
    @inbounds for j = 1:size(a,2)
        r[j] = evaluate_cond(view(a, :, j), view(b, :, j))
    end
    r
end

@inline function evaluate_cond(a, b)
    length(a) == 0 && return 0.0 # comment out and 0.7 is super fast
    @inbounds begin
        s = 0.0
        @simd for I in eachindex(a, b)
            ai = a[I]
            bi = b[I]
            s += abs2(ai - bi)
        end
        return s
    end
end

function perf_colwise_noalloc!(r, a, b)
    @inbounds for j = 1:size(a,2)
        r[j] = evaluate_nocond(view(a, :, j), view(b, :, j))
    end
    r
end

@inline function evaluate_nocond(a, b)
    @inbounds begin
        s = 0.0
        @simd for I in eachindex(a, b)
            ai = a[I]
            bi = b[I]
            s += abs2(ai - bi)
        end
        return s
    end
end

z = zeros(41); A = rand(2, 41); B = rand(2, 41);
g = addgroup!(SUITE, "allocation elision view")
g["conditional"] = @benchmarkable perf_colwise_alloc!($z, $A, $B)
g["no conditional"] = @benchmarkable perf_colwise_noalloc!($z, $A, $B)

####################################################
# Fastmath infererence large number of args #22275 #
####################################################

function f2(a,b,c,d,e,f,g,h,j,k,l,m,n,o,p)
    aidx = eachindex(a)
    @fastmath for i in aidx
        @inbounds a[i] = b[i]+c*(d*e[i]+f*g[i]+h*j[i]+k*l[i]+m*n[i]+o*p[i])
    end
end
a = rand(10)
b = rand(10)
c = 0.1
d = 0.1
e = rand(10)
f = 0.1
g = rand(10)
h = 0.1
j = rand(10)
k = 0.1
l = rand(10)
m = 0.1
n = rand(10)
o = 0.1
p = rand(10)

SUITE["fastmath many args"] = @benchmarkable f2($a,$b,$c,$d,$e,$f,$g,$h,$j,$k,$l,$m,$n,$o,$p)

##############################################################
# Performance and typing of 6+ dimensional generators #21058 #
##############################################################
perf_g6() = sum([+(a,b,c,d,e,f) for a in 1:4, b in 1:4, c in 1:4, d in 1:4, e in 1:4, f in 1:4])

SUITE["perf highdim generator"] = @benchmarkable perf_g6()

end


text:
#-------------------------------------------------------------------------------
# This function creates an instance of MimiRICE2010, given user specifications.
#-------------------------------------------------------------------------------

# Load packages.
using Mimi
using DataFrames
using CSVFiles
using MimiRICE2010

include("new_components/updated_welfare_rice.jl")

# Load necessary model and data files.
un_population = DataFrame(load(joinpath(@__DIR__, "..", "data", "UN_population_rice_regions.csv"), skiplines_begin=3))

# Create a function to construct an updated version of RICE2010.
function create_rice(ρ::Float64, η::Float64, remove_negishi::Bool)

    # ---------------------------------------------
    # Create MimiRICE2010 model and set parameters.
    # ---------------------------------------------

    # Load base version of MimiRICE2010.
    m = MimiRICE2010.get_model()
    #include(joinpath(@__DIR__, "..", "src", "new_components\\Emissions_SlowDown_rice.jl"))
    #include(joinpath(@__DIR__, "..", "src", "new_components\\CatastrophicDamages.jl"))
    #replace_comp!(m, emissions, :emissions, reconnect=true) #unblock this for political feasibility run; need to lower initial guess in cost-min run
    #replace_comp!(m, damages, :damages, reconnect=true)  #unblock this for catastrophic damages output

    # Set savings rate to 25.8% for all regions and time periods.
    update_param!(m, :S, ones(60,12) .* 0.258)

    # Set population to updated UN projctions.
    update_param!(m, :l, convert(Matrix, un_population))

    # Replace welfare component if switching off Negishi weights and then set necessary parameters.
    if remove_negishi == true
        # We cannot use `reconnect` here because the new welfare component does not have all the same
        # parameters as the old one
        replace!(m, :welfare=>welfare, reconnect=false)

        # Set the welfare function parameters
        set_param!(m, :welfare, :ρ, ρ)
        set_param!(m, :welfare, :η, η)

        # Connect the population parameter in the new component to the existing
        # population model level parameter
        connect_param!(m, :welfare, :pop, :l)

        # Connect the welfare component with the net economy component
        connect_param!(m, :welfare, :cpc, :neteconomy, :CPC)
    end

    # Return user-specified model.
    return m
end

ext: jl

max_stars_repo_name: Environment-Research/Utilitarianism

max_issues_repo_name: Environment-Research/Utilitarianism

max_forks_repo_name: Environment-Research/Utilitarianism


text:
## runtests.jl : run all test files

using GaussianRandomFields
using Suppressor
using Plots
using Test

# test indexsets
include("test_covariance_functions.jl")
include("test_cholesky.jl")
include("test_spectral.jl")
include("test_karhunen_loeve.jl")
include("test_circulant_embedding.jl")
include("test_gaussian_random_fields.jl")
include("test_separable.jl")
include("test_fem.jl")
include("test_show.jl")
include("test_plot.jl")


text:
using DecisionTree
using Statistics

function InvFeatureGen(X, y; n_trees::Int=200, min_interval::Int=3)
    n_samps, series_length = size(X)
    transform_xt = Array{Array{Float64,2},1}()
    n_intervals = floor(Int, sqrt(series_length))
    intervals = zeros(Int, n_trees, 3*n_intervals, 2)
    for i in range(1, stop = n_trees)
       transformed_x = Array{Float64,2}(undef, 3*n_intervals, n_samps)
       for j in range(1, stop = n_intervals)
           intervals[i,j,1] = rand(1:(series_length - min_interval))
           len = rand(1:(series_length - intervals[i,j,1]))
           if len < min_interval
               len = min_interval
           end
           intervals[i,j,2] = intervals[i,j,1] + len
           x = Array(1:len+1)
           Y = X[:, intervals[i,j,1]:intervals[i,j,2]]
           means = mean(Y, dims=2)
           stds =  std(Y, dims=2)
           slope = (mean(transpose(x).*Y, dims=2) -
                    mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
           transformed_x[3*j-2,:] =  means
           transformed_x[3*j-1,:] =  stds
           transformed_x[3*j,:]   =  slope
       end
           push!(transform_xt, transpose(transformed_x))
    end
    return transform_xt
end

function RandomForestClassifierTS(transform_xt)
    n_trees = length(transform_xt)
    model = DecisionTreeClassifier()
    forest = Array{DecisionTreeClassifier,1}()
    for i in range(1, stop=n_trees)
        mdl = deepcopy(model)
        fit!(mdl, transform_xt[i], y)
        push!(forest, mdl)
    end
    forest
end

function predict_f(forest, X) #takes only single instance input of Type "transform_xt"
    a, b = proba_predict(forest, X)
    a > b ? 1 : 2
end

function proba_predict(forest, X)
    a = 0
    for tree in forest
        if predict(tree, X) == 1
            a += 1
        end
    end
    (a, length(forest)-a)
end


text:
using ModernGL, GeometryTypes, GLAbstraction, GLWindow, Images

# Load our textures. See "downloads.jl" to get the images.
kitten = load(GLAbstraction.dir("tutorials", "images", "kitten.png"))
puppy  = load(GLAbstraction.dir("tutorials", "images", "puppy.png"))

# Create the window. This sets all the hints and makes the context current.
window = create_glcontext("Textures 2", resolution=(800,600))

vao = glGenVertexArrays()
glBindVertexArray(vao)

# The positions of the vertices in our rectangle
vertex_positions = Point{2,Float32}[(-0.5,  0.5),     # top-left
                                    ( 0.5,  0.5),     # top-right
                                    ( 0.5, -0.5),     # bottom-right
                                    (-0.5, -0.5)]     # bottom-left

# The texture coordinates of each vertex
vertex_texcoords = Vec2f0[(0, 0),
                          (1, 0),
                          (1, 1),
                          (0, 1)]

# Specify how vertices are arranged into faces
# Face{N,T,O} type specifies a face with N vertices, with index type
# T (you should choose UInt32), and index-offset O. If you're
# specifying faces in terms of julia's 1-based indexing, you should set
# O=0. (If you instead number the vertices starting with 0, set
# O=-1.)
elements = Face{3,UInt32,-1}[(0,1,2),          # the first triangle
                             (2,3,0)]          # the second triangle

# The vertex shader---note the `vert` in front of """
vertex_shader = vert"""
#version 150

in vec2 position;
in vec2 texcoord;

out vec2 Texcoord;

void main()
{
    Texcoord = texcoord;
    gl_Position = vec4(position, 0.0, 1.0);
}
"""

# The fragment shader
fragment_shader = frag"""
# version 150

in vec3 Color;
in vec2 Texcoord;

out vec4 outColor;

uniform sampler2D texKitten;
uniform sampler2D texPuppy;

void main()
{
    vec4 colKitten = texture(texKitten, Texcoord);
    vec4 colPuppy  = texture(texPuppy,  Texcoord);
    outColor = mix(colKitten, colPuppy, 0.5);
}
"""

# Link everything together, using the corresponding shader variable as
# the Dict key
bufferdict = Dict(:position=>Buffer(vertex_positions),
                  :texcoord=>Buffer(vertex_texcoords),
                  :texKitten=>Texture(kitten'),
                  :texPuppy=>Texture(puppy'),
                  :indexes=>indexbuffer(elements)) # special for element buffers

ro = std_renderobject(bufferdict,
                      LazyShader(vertex_shader, fragment_shader))

# Draw until we receive a close event
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT)
    GLAbstraction.render(ro)
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
using SafeTestsets
# TODO add description to the tests

## ===-===-===-===-===-===-===-===-
@safetestset "MatrixProcessing tests" begin
    include("matrixProcessing_tests.jl")
end

## ===-===-===-===-===-===-===-===-
@safetestset "MatrixOrganization tests" begin
    include("matrixOrganisation_tests.jl")
end

## ===-===-===-===-===-===-===-===-
@safetestset "BettiCurves tests" begin
    include("bettiCurves_tests.jl")
end

## ===-===-===-===-===-===-===-===-


text:
is_foggy() = rand(Bool)
is_dog_distracted() = rand() < 0.25
loot() = rand(3:13)
loot(crate) = rand(crate)

max_stars_repo_name: tripitakit/julia

max_issues_repo_name: tripitakit/julia

max_forks_repo_name: tripitakit/julia


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MPICH_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MPICH")
JLLWrappers.@generate_main_file("MPICH", UUID("7cb0a576-ebde-5e09-9194-50597f1243b4"))
end  # module MPICH_jll


text:
#=

In the following, ``uhat`` and ``tauhat`` are what the planner would
choose if he could reset at time t, ``uhatdif`` and ``tauhatdif`` are
the difference between those and what the planner is constrained to
choose.  The variable ``mu`` is the Lagrange multiplier associated with
the constraint at time t.

For more complete description of inputs and outputs see the website.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-21

References
----------

Simple port of the file examples/evans_sargent.py

http://quant-econ.net/hist_dep_policies.html

=#
using QuantEcon
using Optim
using PyPlot

type HistDepRamsey
    # These are the parameters of the economy
    A0::Real
    A1::Real
    d::Real
    Q0::Real
    tau0::Real
    mu0::Real
    bet::Real

    # These are the LQ fields and stationary values
    R::Matrix
    A::Matrix
    B::Matrix
    Q::Real
    P::Matrix
    F::Matrix
    lq::LQ
end

type RamseyPath
    y::Matrix
    uhat::Vector
    uhatdif::Vector
    tauhat::Vector
    tauhatdif::Vector
    mu::Vector
    G::Vector
    GPay::Vector
end

function HistDepRamsey(A0, A1, d, Q0, tau0, mu, bet)
    # Create Matrices for solving Ramsey problem
    R = [0.0  -A0/2  0.0    0.0
        -A0/2 A1/2   -mu/2  0.0
        0.0   -mu/2  0.0    0.0
        0.0   0.0    0.0    d/2]

    A = [1.0   0.0  0.0  0.0
         0.0   1.0  0.0  1.0
         0.0   0.0  0.0  0.0
         -A0/d A1/d 0.0  A1/d+1.0/bet]

    B = [0.0 0.0 1.0 1.0/d]'

    Q = 0.0

    # Use LQ to solve the Ramsey Problem.
    lq = LQ(Q, -R, A, B, bet=bet)

    P, F, _d = stationary_values(lq)

    HistDepRamsey(A0, A1, d, Q0, tau0, mu0, bet, R, A, B, Q, P, F, lq)
end

function compute_G(hdr::HistDepRamsey, mu)
    # simplify notation
    Q0, tau0, P, F, d, A, B = hdr.Q0, hdr.tau0, hdr.P, hdr.F, hdr.d, hdr.A, hdr.B
    bet = hdr.bet

    # Need y_0 to compute government tax revenue.
    u0 = compute_u0(hdr, P)
    y0 = vcat([1.0 Q0 tau0]', u0)

    # Define A_F and S matricies
    AF = A - B * F
    S = [0.0 1.0 0.0 0]' * [0.0 0.0 1.0 0]

    # Solves equation (25)
    Omega = solve_discrete_lyapunov(sqrt(bet) .* AF', bet .* AF' * S * AF)
    T0 = y0' * Omega * y0

    return T0[1], A, B, F, P
end

function compute_u0(hdr::HistDepRamsey, P::Matrix)
    # simplify notation
    Q0, tau0 = hdr.Q0, hdr.tau0

    P21 = P[4, 1:3]
    P22 = P[4, 4]
    z0 = [1.0 Q0 tau0]'
    u0 = -P22^(-1) .* P21*(z0)

    return u0[1]
end

function init_path(hdr::HistDepRamsey, mu0, T::Int=20)
    # Construct starting values for the path of the Ramsey economy
    G0, A, B, F, P = compute_G(hdr, mu0)

    # Compute the optimal u0
    u0 = compute_u0(hdr, P)

    # Initialize vectors
    y = Array(Float64, 4, T)
    uhat       = Array(Float64, T)
    uhatdif    = Array(Float64, T)
    tauhat     = Array(Float64, T)
    tauhatdif  = Array(Float64, T-1)
    mu         = Array(Float64, T)
    G          = Array(Float64, T)
    GPay       = Array(Float64, T)

    # Initial conditions
    G[1] = G0
    mu[1] = mu0
    uhatdif[1] = 0.0
    uhat[1] = u0
    y[:, 1] = vcat([1.0 hdr.Q0 hdr.tau0]', u0)

    return RamseyPath(y, uhat, uhatdif, tauhat, tauhatdif, mu, G, GPay)
end

function compute_ramsey_path!(hdr::HistDepRamsey, rp::RamseyPath)
    # simplify notation
    y, uhat, uhatdif, tauhat, = rp.y, rp.uhat, rp.uhatdif, rp.tauhat
    tauhatdif, mu, G, GPay = rp.tauhatdif, rp.mu, rp.G, rp.GPay
    bet = hdr.bet

    G0, A, B, F, P = compute_G(hdr, mu[1])

    for t=2:T
        # iterate government policy
        y[:, t] = (A - B * F) * y[:, t-1]

        # update G
        G[t] = (G[t-1] - bet*y[2, t]*y[3, t])/bet
        GPay[t] = bet.*y[2, t]*y[3, t]

        #=
        Compute the mu if the government were able to reset its plan
        ff is the tax revenues the government would receive if they reset the
        plan with Lagrange multiplier mu minus current G
        =#
        ff(mu) = abs(compute_G(hdr, mu)[1]-G[t])

        # find ff = 0
        mu[t] = optimize(ff, mu[t-1]-1e4, mu[t-1]+1e4).minimum
        temp, Atemp, Btemp, Ftemp, Ptemp = compute_G(hdr, mu[t])

        # Compute alternative decisions
        P21temp = Ptemp[4, 1:3]
        P22temp = P[4, 4]
        uhat[t] = (-P22temp^(-1) .* P21temp * y[1:3, t])[1]

        yhat = (Atemp-Btemp * Ftemp) * [y[1:3, t-1], uhat[t-1]]
        tauhat[t] = yhat[3]
        tauhatdif[t-1] = tauhat[t] - y[3, t]
        uhatdif[t] = uhat[t] - y[3, t]
    end

    return rp
end

function plot1(rp::RamseyPath)
    tt = 1:length(rp.mu)  # tt is used to make the plot time index correct.
    y = rp.y

    n_rows = 3
    fig, axes = subplots(n_rows, 1, figsize=(10, 12))

    subplots_adjust(hspace=0.5)
    for ax in axes
        ax[:grid]()
        ax[:set_xlim](0, 15)
    end

    bbox = (0., 1.02, 1., .102)
    legend_args = {:bbox_to_anchor => bbox, :loc => 3, :mode => "expand"}
    p_args = {:lw => 2, :alpha => 0.7}

    ax = axes[1]
    ax[:plot](tt, squeeze(y[2, :], 1), "b-", label="output"; p_args...)
    ax[:set_ylabel](L"$Q$", fontsize=16)
    ax[:legend](ncol=1; legend_args...)

    ax = axes[2]
    ax[:plot](tt, squeeze(y[3, :], 1), "b-", label="tax rate"; p_args...)
    ax[:set_ylabel](L"$\tau$", fontsize=16)
    ax[:set_yticks]((0.0, 0.2, 0.4, 0.6, 0.8))
    ax[:legend](ncol=1; legend_args...)

    ax = axes[3]
    ax[:plot](tt, squeeze(y[4, :], 1), "b-", label="first difference in output";
              p_args...)
    ax[:set_ylabel](L"$u$", fontsize=16)
    ax[:set_yticks]((0, 100, 200, 300, 400))
    ax[:legend](ncol=1; legend_args...)
    ax[:set_xlabel](L"time", fontsize=16)

end

function plot2(rp::RamseyPath)
    y, uhatdif, tauhatdif, mu = rp.y, rp.uhatdif, rp.tauhatdif, rp.mu
    G, GPay = rp.G, rp.GPay
    T = length(rp.mu)
    tt = 1:T  # tt is used to make the plot time index correct.
    tt2 = 1:T-1

    n_rows = 4
    fig, axes = subplots(n_rows, 1, figsize=(10, 16))

    subplots_adjust(hspace=0.5)
    for ax in axes
        ax[:grid](alpha=.5)
        ax[:set_xlim](-0.5, 15)
    end

    bbox = (0., 1.02, 1., .102)
    legend_args = {:bbox_to_anchor => bbox, :loc => 3, :mode => "expand"}
    p_args = {:lw => 2, :alpha => 0.7}

    ax = axes[1]
    ax[:plot](tt2, tauhatdif,
              label="time inconsistency differential for tax rate"; p_args...)
    ax[:set_ylabel](L"$\Delta\tau$", fontsize=16)
    ax[:set_yticks]((0.0, 0.4, 0.8, 1.2))
    ax[:legend](ncol=1; legend_args...)

    ax = axes[2]
    ax[:plot](tt, uhatdif,
              label=L"time inconsistency differential for $u$"; p_args...)
    ax[:set_ylabel](L"$\Delta u$", fontsize=16)
    ax[:set_yticks]((-3.0, -2.0, -1.0, 0.0))
    ax[:legend](ncol=1; legend_args...)

    ax = axes[3]
    ax[:plot](tt, mu, label="Lagrange multiplier"; p_args...)
    ax[:set_ylabel](L"$\mu$", fontsize=16)
    ax[:set_yticks]((2.34e-3, 2.43e-3, 2.52e-3))
    ax[:legend](ncol=1; legend_args...)

    ax = axes[4]
    ax[:plot](tt, G, label="government revenue"; p_args...)
    ax[:set_ylabel](L"$G$", fontsize=16)
    ax[:set_yticks]((9200, 9400, 9600, 9800))
    ax[:legend](ncol=1; legend_args...)

    ax[:set_xlabel](L"time", fontsize=16)

end

# Primitives
T    = 20
A0   = 100.0
A1   = 0.05
d    = 0.20
bet = 0.95

# Initial conditions
mu0  = 0.0025
Q0   = 1000.0
tau0 = 0.0

# Solve Ramsey problem and compute path
hdr = HistDepRamsey(A0, A1, d, Q0, tau0, mu0, bet)
rp = init_path(hdr, mu0, T)
compute_ramsey_path!(hdr, rp)  # updates rp in place
plot1(rp)
plot2(rp)

max_stars_repo_name: gsulis/structural

max_issues_repo_name: gsulis/structural

max_forks_repo_name: gsulis/structural


text:
"""
    BipartiteView{T<:Real} <: LightGraphs.SimpleGraphs.AbstractSimpleGraph{Int}

Represents a bipartite view of a hypergraph `h`.
Note this is a view - changes to the original hypergraph will be automatically reflected in the view.

**Constructors**

BipartiteView(::Hypergraph)

The bipartite view of a hypergraph is suitable for processing with the LightGraphs.jl package.
Several LightGraphs methods are provided for the compability.

"""
struct BipartiteView{T<:Real} <: LightGraphs.SimpleGraphs.AbstractSimpleGraph{Int}
    h::Hypergraph{T}
end

"""
  Return the number of vertices in a bipartite view `b` of a hypergraph.
"""
LightGraphs.nv(b::BipartiteView) = length(b.h.v2he)+length(b.h.he2v)

LightGraphs.vertices(b::BipartiteView) = Base.OneTo(LightGraphs.nv(b))

"""
  Return the number of edges in a bipartite view `b` of a hypergraph.
"""
LightGraphs.ne(b::BipartiteView) = sum(length.(b.h.v2he))

function LightGraphs.all_neighbors(b::BipartiteView, v::Integer)
    n1 = length(b.h.v2he)
    if v <= n1
      n1 .+ keys(b.h.v2he[v])
    else
      collect(keys(b.h.he2v[v-n1]))
    end
end

function LightGraphs.has_edge(b::BipartiteView, s, d)
    n1 = length(b.h.v2he)
    if s <= n1
        d > n1 && has_key(b.v2he[s], d - n1)
    else
        d <= n1 && has_key(b.h.he2v[s-n1], d)
    end
end

LightGraphs.has_vertex(b::BipartiteView, v::Integer) = 1 <= v <= LightGraphs.nv(b)

LightGraphs.outneighbors(b::BipartiteView, v::Integer) = LightGraphs.all_neighbors(b::BipartiteView, v)

LightGraphs.inneighbors(b::BipartiteView, v::Integer) = LightGraphs.all_neighbors(b::BipartiteView, v)

"""
    LightGraphs.SimpleGraph(b::BipartiteView)

Creates a `LightGraphs.SimpleGraph` representation of a `BipartiteView` b.

This creates a copy of the date. Note that the weights information is not stored
in the created `SimpleGraph`.
"""
function LightGraphs.SimpleGraph(b::BipartiteView)
    g = SimpleGraph(nv(b))
    for v in keys(b.h.v2he)
        for he in keys(b.h.v2he[v])
            add_edge!(g, v, length(b.h.v2he) + he)
        end
    end
    g
end

LightGraphs.is_directed(b::BipartiteView{T}) where T = false

LightGraphs.is_directed(::Type{BipartiteView{T}}) where T = false

Base.eltype(::BipartiteView{T}) where T = Int

"""
    shortest_path(b::BipartiteView,source::Int, target::Int)

Finds a single shortest path in a graph `b` between vertices
`source` and `target`.
Note that if several paths of the same length exist, only one
will be returned.

"""
function shortest_path(b::BipartiteView,source::Int, target::Int)
    checkbounds(b.h.v2he, source)
    checkbounds(b.h.v2he, target)
    dj = dijkstra_shortest_paths(b, source)
    enumerate_paths(dj)[target][1:2:end]
end

"""
    LightGraphs.SimpleGraphs.fadj(b::BipartiteView)

Generates an adjency list for this view of a hypergraph.
"""
function LightGraphs.SimpleGraphs.fadj(b::BipartiteView)
    res = Vector{Vector{Int}}(undef, LightGraphs.nv(b))

    h_nv = length(b.h.v2he)
    for i in 1:h_nv
       res[i] = h_nv .+ sort!(collect(keys(b.h.v2he[i])))
    end
    for i in 1:length(b.h.he2v)
        res[i+h_nv] = sort!(collect(keys(b.h.he2v[i])))
    end
    res
end

LightGraphs.SimpleGraphs.fadj(b::BipartiteView, v::Integer) = LightGraphs.all_neighbors(b,v)
LightGraphs.edges(b::BipartiteView) = LightGraphs.SimpleGraphs.SimpleEdgeIter(b)

LightGraphs.edgetype(b::BipartiteView{T}) where T = LightGraphs.SimpleGraphs.SimpleEdge{Int}

LightGraphs.zero(t::BipartiteView{T}) where T = BipartiteView(Hypergraph{T}(0,0))
LightGraphs.zero(::Type{BipartiteView{T}}) where T = BipartiteView(Hypergraph{T}(0,0))


text:
using ComoniconOptions
using Test

@testset "load options" begin
    options = read_options(joinpath(pkgdir(ComoniconOptions), "test"))
    display(options)
    @test options == ComoniconOptions.Comonicon(;
        name = "foo",
        install = ComoniconOptions.Install(;
            compile = "min",
        ),
        sysimg = ComoniconOptions.SysImg(;
            cpu_target = "native",
            precompile = ComoniconOptions.Precompile(;
                execution_file = ["deps/precopmile.jl"],
            ),
        ),
        download = ComoniconOptions.Download(;
            user = "Roger-luo",
            repo = "Foo.jl",
        ),
        application = ComoniconOptions.Application(;
            assets = ComoniconOptions.Asset[asset"PkgTemplate: templates", asset"assets/images"],
            incremental = true,
            filter_stdlibs = false,
        ),
    )
end

@testset "has_comonicon_toml" begin
    @test ComoniconOptions.has_comonicon_toml(Main) == false
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm
using NearestNeighbors
using PyCall

p = PyCall.pyimport("pyplasm")

filename = "/Users/paoluzzi/Documents/dev/Plasm.jl/test/svg/tile.svg"
V,EV = Plasm.svg2lar(filename)
Plasm.view(Plasm.numbering(.2)((V,[[[k] for k=1:size(V,2)], EV])))

function pointsinout(V,EV, n=10000)
	result = [Plasm.lar2hpc(V,EV)]
	classify = Lar.pointInPolygonClassification(V,EV)
	Threads.@threads for k=1:n
		queryPoint = [rand(),rand()]
		inOut = classify(queryPoint)
		# println("k = $k, queryPoint = $queryPoint, inOut = $inOut")
		if inOut=="p_in"
			push!(result, p["MK"](queryPoint))
		elseif inOut=="p_out"
			push!(result, Plasm.color("red")(p["MK"](queryPoint)))
		elseif inOut=="p_on"
			push!(result, Plasm.color("green")(p["MK"](queryPoint)))
		end
	end
	return result
end

result = @time pointsinout(V,EV);
Plasm.view(p["STRUCT"](result))


text:
#!/usr/bin/env julia
# Sets and Dictionaries
# https://teaching.healthtech.dtu.dk/22110/index.php/Sets_and_Dictionaries

const codons = Dict("TTT"=>"F", "TTC"=>"F", "TTA"=>"L", "TTG"=>"L",
    "TCT"=>"S", "TCC"=>"S", "TCA"=>"S", "TCG"=>"S",
    "TAT"=>"Y", "TAC"=>"Y", "TAA"=>"", "TAG"=>"",
    "TGT"=>"C", "TGC"=>"C", "TGA"=>"", "TGG"=>"W",
    "CTT"=>"L", "CTC"=>"L", "CTA"=>"L", "CTG"=>"L",
    "CCT"=>"P", "CCC"=>"P", "CCA"=>"P", "CCG"=>"P",
    "CAT"=>"H", "CAC"=>"H", "CAA"=>"Q", "CAG"=>"Q",
    "CGT"=>"R", "CGC"=>"R", "CGA"=>"R", "CGG"=>"R",
    "ATT"=>"I", "ATC"=>"I", "ATA"=>"I", "ATG"=>"M",
    "ACT"=>"T", "ACC"=>"T", "ACA"=>"T", "ACG"=>"T",
    "AAT"=>"N", "AAC"=>"N", "AAA"=>"K", "AAG"=>"K",
    "AGT"=>"S", "AGC"=>"S", "AGA"=>"R", "AGG"=>"R",
    "GTT"=>"V", "GTC"=>"V", "GTA"=>"V", "GTG"=>"V",
    "GCT"=>"A", "GCC"=>"A", "GCA"=>"A", "GCG"=>"A",
    "GAT"=>"D", "GAC"=>"D", "GAA"=>"E", "GAG"=>"E",
    "GGT"=>"G", "GGC"=>"G", "GGA"=>"G", "GGG"=>"G")

function translate_dna(dna)
    if length(dna) % 3 != 0
        println("Whoops")
    end
    peptide = ""
    for i = 1:3:length(dna)
        peptide *= codons[dna[i:i+2]]
    end
    peptide
end

function read_fasta(f)
    headlines = String[]
    seqs = String[]
    seq = ""
    for line = eachline(f)
        ln = chomp(line)
        if ln[1] == '>'
            if seq != ""
                push!(seqs, seq)
                seq = ""
            end
            push!(headlines, ln)
        else
            seq *= translate_dna(ln)
        end
    end
    push!(seqs, seq)
    headlines, seqs
end

open("data/aa7.fsa", "w") do outfile
    headlines, seqs = open(read_fasta, "data/dna7.fsa")
    for (headline, seq) in zip(headlines, seqs)
        println(outfile, "$headline Amino Acid Sequence")
        for i in 1:60:length(seq)
            if i + 59 > length(seq)
                println(outfile, seq[i:end])
            else
                println(outfile, seq[i:(i + 59)])
            end
        end
    end
end

# Ex 3
function read_as_set(f, colnum)
    extract_line = if colnum === 0
        ef1(ln) = string(chomp(ln))
    else
        ef2(ln) = string(split(ln, "\t")[colnum])
    end
    outset = Set{String}()
    for line = eachline(f)
        push!(outset, extract_line(line))
    end
    outset
end

function read_set_diff(f1, f2)
    set1 = read_as_set(f1, 0)
    set2 = read_as_set(f2, 2)
    setdiff(set1,set2)
end

println(read_set_diff("data/start10.dat", "data/res10.dat"))

# Ex 4
function read_as_count_dict(f)
    outdict = Dict{String, Int64}()
    for line = eachline(f)
        x = chomp(line)
        outdict[x] = get(outdict, x, 0) + 1
    end
    outdict
end

open("data/ex5.acc.count", "w") do f
    for (id, count) in sort(collect(open(read_as_count_dict, "data/ex5.acc")),
            by=x->x[2], rev=true)
        println(f, id, " ", count)
    end
end

# Ex 5
read_data() = try
    if length(ARGS) == 1
        filename = ARGS[1]
    else
        filename = "data/data1.gb"
    end
    open(_read_data, filename)
catch e
    msg = sprint(showerror, e)
    println(msg)
    exit(1)
end

function _read_data(f)
    MED = String[]
    SEQ_flag = AA_flag = CDS_flag = false
    SEQ = AA = CDS = ""
    CDS_range = ""

    for line in eachline(f)
        if startswith(line, "//")
            SEQ_flag = false
        end
        if SEQ_flag
            SEQ_re = findall(r"([atcg]+)", line)
            SEQ_parts = [line[x] for x in SEQ_re]
            SEQ *= join(SEQ_parts)
        else
            ACC_re = match(r"^ACCESSION\s+(\S+)$", line)
            DEF_re = match(r"^DEFINITION\s+(.+)$", line)
            ORG_re = match(r"^\s+ORGANISM\s+(.+)$", line)
            MED_re = match(r"^\s+MEDLINE\s+(\S+)$", line)

            if ACC_re !== nothing
                println("Accession number: ", ACC_re[1])
            elseif DEF_re !== nothing
                println("Definition: ", DEF_re[1])
            elseif ORG_re !== nothing
                println("Organism: ", ORG_re[1])
            elseif MED_re !== nothing
                push!(MED, MED_re[1])
            end
        end

        if startswith(line, "ORIGIN")
            SEQ_flag = true
        end

        CDS_re = match(r"^\s+CDS\s+(\S+)$", line)
        AA_re = match(r"^\s+/translation=(\S+)$", line)

        if CDS_flag
            CDS_data = strip(line)
            if startswith(CDS_data, "/")
                CDS_flag = false
            else
                CDS_range *= CDS_data
            end
        end
        if CDS_re !== nothing
            CDS_range *= CDS_re[1]
            CDS_flag = true
        end

        if AA_flag
            AA_data = match(r"^\s+(\S+)$", line)
            if AA_data !== nothing
                AA *= AA_data[1]
            else
                AA_flag = false
            end
        end
        if AA_re !== nothing
            AA *= AA_re[1]
            AA_flag = true
        end
    end

    CDS_list = [split(CDS_range[x], "..") for x in findall(r"(\d+\.\.\d+)", CDS_range)]

    for x in CDS_list
        CDS *= SEQ[parse(Int64, x[1]):parse(Int64, x[2])]
    end

    CDS_dict = Dict{String, Int64}()
    for i = 1:3:length(CDS)
        codon = CDS[i:(i+2)]
        CDS_dict[codon] = get(CDS_dict, codon, 0) + 1
    end

    println("MEDLINE article number(s): ", join(MED, ", "))
    println()
    println("Sequence: ", SEQ)
    println()
    println("Coding sequence (CDS): ", CDS)
    println()
    println("Amino acid sequence: ", AA)
    println()
    println("Codons: ")
    for (codon, count) in sort(collect(CDS_dict), by=x->x[2], rev=true)
        println(codon, ": ", count)
    end
end

read_data()

# Ex 6

read_authors() = try
    if length(ARGS) == 1
        filename = ARGS[1]
    else
        filename = "data/data1.gb"
    end
    open(_read_authors, "data/data1.gb")
catch e
    msg = sprint(showerror, e)
    println(msg)
    #exit(1)
end

function _read_authors(f)
    MED = String[]
    MED_flag = false

    for line in eachline(f)
        MED_reg = match(r"^\s+AUTHORS\s+(.*)$", line)

        if MED_flag
            MED_extra_reg = match(r"^\s+[A-Z]+\s", line)
            if MED_extra_reg === nothing
                MED[end] *= strip(line)
            else
                MED_flag = false
            end
        end

        if MED_reg !== nothing
            push!(MED, MED_reg[1])
            MED_flag = true
        end
    end

    MED = [replace(x, " and" => ", ") for x in MED]

    authors_arr = String[]
    for names in MED
        append!(authors_arr, [strip(name) for name in split(names, ", ")])
    end

    authors = Set(authors_arr)

    println(MED)
    println(authors)
    authors
end

read_authors()


text:
"""
    qubits(N::Int; mixed::Bool=false)

    qubits(sites::Vector{<:Index}; mixed::Bool=false)

Initialize qubits to:
- An MPS wavefunction `|ψ⟩` if `mixed = false`
- An MPO density matrix `ρ` if `mixed = true`
"""
qubits(N::Int; mixed::Bool=false) = qubits(siteinds("Qubit", N); mixed=mixed)

function qubits(sites::Vector{<:Index}; mixed::Bool=false)
  @warn "Method `qubits` is deprecated, use `productstate` or `productoperator` instead."
  ψ = productMPS(sites, "0")
  mixed && return MPO(ψ)
  return ψ
end

"""
    qubits(M::Union{MPS,MPO,LPDO}; mixed::Bool=false)

Initialize qubits on the Hilbert space of a reference state,
given as `MPS`, `MPO` or `LPDO`.
"""
qubits(M::Union{MPS,MPO,LPDO}; mixed::Bool=false) = qubits(hilbertspace(M); mixed=mixed)

"""
    qubits(N::Int, states::Vector{String}; mixed::Bool=false)

    qubits(sites::Vector{<:Index}, states::Vector{String};mixed::Bool = false)

Initialize the qubits to a given single-qubit product state.
"""
function qubits(N::Int, states::Vector{String}; mixed::Bool=false)
  return qubits(siteinds("Qubit", N), states; mixed=mixed)
end

function qubits(sites::Vector{<:Index}, states::Vector{String}; mixed::Bool=false)
  @warn "Method `qubits` is deprecated, use `productstate` or `productoperator` instead."
  N = length(sites)
  @assert N == length(states)

  ψ = productMPS(sites, "0")

  if N == 1
    s1 = sites[1]
    state1 = state(states[1])
    if eltype(state1) <: Complex
      ψ[1] = complex(ψ[1])
    end
    for j in 1:dim(s1)
      ψ[1][s1 => j] = state1[j]
    end
    mixed && return MPO(ψ)
    return ψ
  end

  # Set first site
  s1 = sites[1]
  l1 = linkind(ψ, 1)
  state1 = state(states[1])
  if eltype(state1) <: Complex
    ψ[1] = complex(ψ[1])
  end
  for j in 1:dim(s1)
    ψ[1][s1 => j, l1 => 1] = state1[j]
  end

  # Set sites 2:N-1
  for n in 2:(N - 1)
    sn = sites[n]
    ln_1 = linkind(ψ, n - 1)
    ln = linkind(ψ, n)
    state_n = state(states[n])
    if eltype(state_n) <: Complex
      ψ[n] = complex(ψ[n])
    end
    for j in 1:dim(sn)
      ψ[n][sn => j, ln_1 => 1, ln => 1] = state_n[j]
    end
  end

  # Set last site N
  sN = sites[N]
  lN_1 = linkind(ψ, N - 1)
  state_N = state(states[N])
  if eltype(state_N) <: Complex
    ψ[N] = complex(ψ[N])
  end
  for j in 1:dim(sN)
    ψ[N][sN => j, lN_1 => 1] = state_N[j]
  end

  mixed && return MPO(ψ)
  return ψ
end


text:
#!/usr/bin/julia
"""
Collection of input/output functions.
"""

function write_xyz(fname, bl, xyz)
    """Input:
    * bl: (N, 1) vector of names
    * xyz: (N, 3) matrix of positions
    """
    N = size(bl, 1)
    open(fname, "w") do f
        write(f, "$(N)\nbla\n")
        for i in 1:N
            write(f, "$(bl[i])\t$(xyz[i, 1])\t$(xyz[i, 2])\t$(xyz[i, 3])\n")
        end
    end
end

function read_xyz(fname)
    s = readlines(open(fname))[3:end]
    N = length(s)
    bl = zeros(Int, N)
    xyz = zeros(N, 3)
    for i in 1:N
        tmp = split(s[i])
        bl[i] = parse(Int, tmp[1])
        xyz[i, :] = map(x->parse(Float64, x), tmp[2:4])
    end
    bl, xyz
end

function read_config(fname)
    """Returns only positions and velocities
    not yet atom names."""
    s = readlines(open(fname))
    Ns = length(s)
    println("Number of lines: $Ns")
    vel_flag = parse(Int, split(s[2])[2])
    println("Velocity flag: $(vel_flag)")
    if vel_flag == 0
        N = Int((Ns - 2) / 2)
        println("Number of atoms: $N")
        mask_X = collect(4:2:Ns)
        X = zeros(N, 3)
        for i in 1:N
            tmp = split(s[mask_X[i]])
            X[i, :] = map(x->parse(Float64, x), tmp[1:3])
        end
        X
    elseif vel_flag == 1
        N = Int((Ns - 2) / 3)
        println("Number of atoms: $N")
        mask_X = collect(4:3:Ns)
        mask_V = collect(5:3:Ns)
        X = zeros(N, 3)
        V = zeros(N, 3)
        for i in 1:N
            tmp = split(s[mask_X[i]])
            X[i, :] = map(x->parse(Float64, x), tmp[1:3])
            tmp = split(s[mask_V[i]])
            V[i, :] = map(x->parse(Float64, x), tmp[1:3])
        end
        X, V
    end
end

max_stars_repo_name: petervanya/DPDsim

max_issues_repo_name: petervanya/DPDsim

max_forks_repo_name: petervanya/DPDsim


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

using Random

@testset "gcd/lcm" begin
    # All Integer data types take different code paths -- test all
    # TODO: Test gcd and lcm for BigInt.
    for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128)
        @test gcd(T(3)) === T(3)
        @test gcd(T(3), T(5)) === T(1)
        @test gcd(T(3), T(15)) === T(3)
        @test gcd(T(0), T(15)) === T(15)
        @test gcd(T(15), T(0)) === T(15)
        if T <: Signed
            @test gcd(T(-12)) === T(12)
            @test gcd(T(0), T(-15)) === T(15)
            @test gcd(T(-15), T(0)) === T(15)
            @test gcd(T(3), T(-15)) === T(3)
            @test gcd(T(-3), T(-15)) === T(3)
        end
        @test gcd(T(0), T(0)) === T(0)

        @test gcd(T(2), T(4), T(6)) === T(2)
        if T <: Signed
            @test gcd(T(2), T(4), T(-6)) === T(2)
            @test gcd(T(2), T(-4), T(-6)) === T(2)
            @test gcd(T(-2), T(4), T(-6)) === T(2)
            @test gcd(T(-2), T(-4), T(-6)) === T(2)
        end

        @test gcd(typemax(T), T(1)) === T(1)
        @test gcd(T(1), typemax(T)) === T(1)
        @test gcd(typemax(T), T(0)) === typemax(T)
        @test gcd(T(0), typemax(T)) === typemax(T)
        @test gcd(typemax(T), typemax(T)) === typemax(T)
        @test gcd(typemax(T), typemax(T)-T(1)) === T(1)     # gcd(n, n-1) = 1. n and n-1 are always coprime.

        if T <: Signed
            @test gcd(-typemax(T), T(1)) === T(1)
            @test gcd(T(1), -typemax(T)) === T(1)
            @test gcd(-typemax(T), T(0)) === typemax(T)
            @test gcd(T(0), -typemax(T)) === typemax(T)
            @test gcd(-typemax(T), -typemax(T)) === typemax(T)
            @test gcd(typemax(T), -typemax(T)) === typemax(T)
            @test gcd(-typemax(T), typemax(T)) === typemax(T)

            @test gcd(typemin(T), T(1)) === T(1)
            @test gcd(T(1), typemin(T)) === T(1)
            @test gcd(typemin(T), typemin(T)+T(1)) === T(1) # gcd(n, n+1) = 1. n and n+1 are always coprime.
            @test_throws OverflowError gcd(typemin(T), typemin(T))
            @test_throws OverflowError gcd(typemin(T), T(0))
            @test_throws OverflowError gcd(T(0), typemin(T))
        else
            # For Unsigned Integer types, -typemax(T) == 1.
            @test gcd(-typemax(T), T(1)) === T(1)
            @test gcd(T(1), -typemax(T)) === T(1)
            @test gcd(-typemax(T), T(0)) === T(1)
            @test gcd(T(0), -typemax(T)) === T(1)
            @test gcd(-typemax(T), -typemax(T)) === T(1)
            @test gcd(-typemax(T), typemax(T)) === T(1)
            @test gcd(typemax(T), -typemax(T)) === T(1)

            # For Unsigned Integer types, typemin(T) == 0.
            @test gcd(typemin(T), T(1)) === T(1)
            @test gcd(T(1), typemin(T)) === T(1)
            @test gcd(typemin(T), typemin(T)+T(1)) === T(1) # gcd(n, n+1) = 1. n and n+1 are always coprime.
            @test gcd(typemin(T), typemin(T)) === T(0)
            @test gcd(typemin(T), T(0)) === T(0)
            @test gcd(T(0), typemin(T)) === T(0)
        end

        @test lcm(T(0)) === T(0)
        @test lcm(T(2)) === T(2)
        @test lcm(T(2), T(3)) === T(6)
        @test lcm(T(3), T(2)) === T(6)
        @test lcm(T(4), T(6)) === T(12)
        @test lcm(T(6), T(4)) === T(12)
        @test lcm(T(3), T(0)) === T(0)
        @test lcm(T(0), T(3)) === T(0)
        @test lcm(T(0), T(0)) === T(0)
        if T <: Signed
            @test lcm(T(-12)) === T(12)
            @test lcm(T(0), T(-4)) === T(0)
            @test lcm(T(-4), T(0)) === T(0)
            @test lcm(T(4), T(-6)) === T(12)
            @test lcm(T(-4), T(-6)) === T(12)
        end

        @test lcm(T(2), T(4), T(6)) === T(12)
        @test lcm(T(2), T(4), T(0)) === T(0)
        if T <: Signed
            @test lcm(T(2), T(4), T(-6)) === T(12)
            @test lcm(T(2), T(-4), T(-6)) === T(12)
            @test lcm(T(-2), T(-4), T(-6)) === T(12)
            @test lcm(T(-2), T(0), T(-6)) === T(0)
        end

        @test lcm(typemax(T), T(1)) === typemax(T)
        @test lcm(T(1), typemax(T)) === typemax(T)
        @test lcm(typemax(T), T(0)) === T(0)
        @test lcm(T(0), typemax(T)) === T(0)
        @test lcm(typemax(T), typemax(T)) === typemax(T)
        @test_throws OverflowError lcm(typemax(T), typemax(T)-T(1)) # lcm(n, n-1) = n*(n-1). Since n and n-1 are always coprime.
        @test_throws OverflowError lcm(typemax(T), T(2))

        let x = isqrt(typemax(T))+T(1) # smallest number x such that x^2 > typemax(T)
            @test lcm(x, x) === x
            @test_throws OverflowError lcm(x, x+T(1))   # lcm(n, n+1) = n*(n+1). Since n and n+1 are always coprime.
        end

        if T <: Signed
            @test lcm(-typemax(T), T(1)) === typemax(T)
            @test lcm(T(1), -typemax(T)) === typemax(T)
            @test lcm(-typemax(T), T(0)) === T(0)
            @test lcm(T(0), -typemax(T)) === T(0)
            @test lcm(-typemax(T), -typemax(T)) === typemax(T)
            @test lcm(typemax(T), -typemax(T)) === typemax(T)
            @test lcm(-typemax(T), typemax(T)) === typemax(T)

            @test_throws OverflowError lcm(typemin(T), T(1))
            @test_throws OverflowError lcm(T(1), typemin(T))
            @test lcm(typemin(T), T(0)) === T(0)
            @test lcm(T(0), typemin(T)) === T(0)
            @test_throws OverflowError lcm(typemin(T), typemin(T)+T(1)) # lcm(n, n+1) = n*(n+1).
            @test_throws OverflowError lcm(typemin(T), typemin(T))
        else
            # For Unsigned Integer types, -typemax(T) == 1.
            @test lcm(-typemax(T), T(1)) === T(1)
            @test lcm(T(1), -typemax(T)) === T(1)
            @test lcm(-typemax(T), T(0)) === T(0)
            @test lcm(T(0), -typemax(T)) === T(0)
            @test lcm(-typemax(T), -typemax(T)) === T(1)
            @test lcm(-typemax(T), typemax(T)) === typemax(T)
            @test lcm(typemax(T), -typemax(T)) === typemax(T)

            # For Unsigned Integer types, typemin(T) == 0.
            @test lcm(typemin(T), T(1)) === lcm(T(0), T(1)) === T(0)
            @test lcm(T(1), typemin(T)) === T(0)
            @test lcm(typemin(T), T(0)) === T(0)
            @test lcm(T(0), typemin(T)) === T(0)
            @test lcm(typemin(T), typemin(T)) === T(0)
            @test lcm(typemin(T), typemin(T)+T(1)) === T(0)
        end
    end
    @test lcm(0x5, 3) == 15
    @test gcd(0xf, 20) == 5
    @test gcd(UInt32(6), Int8(-50)) == 2
    @test gcd(typemax(UInt), -16) == 1
end

@testset "gcd/lcm for arrays" begin
    # TODO: Test gcd and lcm for BigInt arrays.
    for T in (Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128)
        @test gcd(T[]) === T(0)
        @test gcd(T[3, 5]) === T(1)
        @test gcd(T[3, 15]) === T(3)
        @test gcd(T[0, 15]) === T(15)
        if T <: Signed
            @test gcd(T[-12]) === T(12)
            @test gcd(T[3,-15]) === T(3)
            @test gcd(T[-3,-15]) === T(3)
        end
        @test gcd(T[0, 0]) === T(0)

        @test gcd(T[2, 4, 6]) === T(2)
        @test gcd(T[2, 4, 3, 5]) === T(1)

        @test lcm(T[]) === T(1)
        @test lcm(T[2, 3]) === T(6)
        @test lcm(T[4, 6]) === T(12)
        @test lcm(T[3, 0]) === T(0)
        @test lcm(T[0, 0]) === T(0)
        if T <: Signed
            @test lcm(T[-2]) === T(2)
            @test lcm(T[4, -6]) === T(12)
            @test lcm(T[-4, -6]) === T(12)
        end

        @test lcm(T[2, 4, 6]) === T(12)
    end
end

@testset "gcdx" begin
    # TODO: Test gcdx for BigInt.
    for T in (Int8, Int16, Int32, Int64, Int128)
        @test gcdx(T(5), T(12)) === (T(1), T(5), T(-2))
        @test gcdx(T(5), T(-12)) === (T(1), T(5), T(2))
        @test gcdx(T(-5), T(12)) === (T(1), T(-5), T(-2))
        @test gcdx(T(-5), T(-12)) === (T(1), T(-5), T(2))
        @test gcdx(T(-25), T(-4)) === (T(1), T(-1), T(6))
    end
    x, y = Int8(-12), UInt(100)
    d, u, v = gcdx(x, y)
    @test x*u + y*v == d
end

@testset "gcd/lcm/gcdx for custom types" begin
    struct MyRational <: Real
        val::Rational{Int}
    end
    Base.promote_rule(::Type{MyRational}, T::Type{<:Real}) = promote_type(Rational{Int}, T)
    (T::Type{<:Real})(x::MyRational) = T(x.val)

    @test gcd(MyRational(2//3), 3) == gcd(2//3, 3) == gcd(Real[MyRational(2//3), 3])
    @test lcm(MyRational(2//3), 3) == lcm(2//3, 3) == lcm(Real[MyRational(2//3), 3])
    @test gcdx(MyRational(2//3), 3) == gcdx(2//3, 3)
end

@testset "invmod" begin
    @test invmod(6, 31) === 26
    @test invmod(-1, 3) === 2
    @test invmod(1, -3) === -2
    @test invmod(-1, -3) === -1
    @test invmod(0x2, 0x3) === 0x2
    @test invmod(2, 0x3) === UInt(2)
    @test invmod(0x8, -3) === -1
    @test_throws DomainError invmod(0, 3)

    # For issue 29971
    @test invmod(UInt8(1), typemax(UInt8))  === 0x01
    @test invmod(UInt16(1), typemax(UInt16)) === 0x0001
    @test invmod(UInt32(1), typemax(UInt32)) === 0x0000_0001
    @test invmod(UInt64(1), typemax(UInt64)) === 0x0000_0000_0000_0001

    for T in (UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128, BigInt)
        @test invmod(T(3), T(124))::T == 83
    end

    for T in (Int8, UInt8)
        for x in typemin(T):typemax(T)
            for m in typemin(T):typemax(T)
                if m != 0 && try gcdx(x, m)[1] == 1 catch _ true end
                    y = invmod(x, m)
                    @test mod(widemul(y, x), m) == mod(1, m)
                    @test div(y, m) == 0
                else
                    @test_throws DomainError invmod(x, m)
                end
            end
        end
    end
end

@testset "powermod" begin
    @test powermod(2, 3, 5) == 3
    @test powermod(2, 3, -5) == -2

    @test powermod(2, 0, 5) == 1
    @test powermod(2, 0, -5) == -4

    @test powermod(2, -1, 5) == 3
    @test powermod(2, -2, 5) == 4
    @test powermod(2, -1, -5) == -2
    @test powermod(2, -2, -5) == -1
end

@testset "nextpow/prevpow" begin
    @test nextpow(2, 3) == 4
    @test nextpow(2, 4) == 4
    @test nextpow(2, 7) == 8
    @test_throws DomainError nextpow(0, 3)
    @test_throws DomainError nextpow(3, 0)

    @test prevpow(2, 3) == 2
    @test prevpow(2, 4) == 4
    @test prevpow(2, 5) == 4
    @test_throws DomainError prevpow(0, 3)
    @test_throws DomainError prevpow(0, 3)
end

@testset "ndigits/ndigits0z" begin
    @testset "issue #8266" begin
        @test ndigits(-15, base=10) == 2
        @test ndigits(-15, base=-10) == 2
        @test ndigits(-1, base=10) == 1
        @test ndigits(-1, base=-10) == 2
        @test ndigits(2, base=10) == 1
        @test ndigits(2, base=-10) == 1
        @test ndigits(10, base=10) == 2
        @test ndigits(10, base=-10) == 3
        @test ndigits(17, base=10) == 2
        @test ndigits(17, base=-10) == 3
        @test ndigits(unsigned(17), base=-10) == 3

        @test ndigits(146, base=-3) == 5
    end
    @testset "ndigits with base power of 2" begin
        @test ndigits(17, base = 2) == 5
        @test ndigits(123, base = 4) == 4
        @test ndigits(64, base = 8) == 3
        @test ndigits(8436, base = 16) == 4
        @test ndigits(159753, base = 32) == 4
        @test ndigits(3578951, base = 64) == 4
    end
    let (n, b) = rand(Int, 2)
        -1 <= b <= 1 && (b = 2) # invalid bases
        @test ndigits(n) == ndigits(big(n)) == ndigits(n, base=10)
        @test ndigits(n, base=b) == ndigits(big(n), base=b)
    end

    for b in -1:1
        @test_throws DomainError ndigits(rand(Int), base=b)
    end
    @test ndigits(Int8(5)) == ndigits(5)

    # issue #19367
    @test ndigits(Int128(2)^64, base=256) == 9

    # test unsigned bases
    @test ndigits(9, base=0x2) == 4
    @test ndigits(0x9, base=0x2) == 4

    # ndigits is defined for Bool
    @test iszero([Base.ndigits0z(false, b) for b in [-20:-2;2:20]])
    @test all(n -> n == 1, Base.ndigits0z(true, b) for b in [-20:-2;2:20])
    @test all(n -> n == 1, ndigits(x, base=b) for b in [-20:-2;2:20] for x in [true, false])

    # issue #29148
    @test ndigits(typemax(UInt64), base=-2) == ndigits(big(typemax(UInt64)), base=-2)
    for T in Base.BitInteger_types
        n = rand(T)
        b = -rand(2:100)
        @test ndigits(n, base=b) == ndigits(big(n), base=b)
    end

end

primitive type BitString128 128 end

@testset "bin/oct/dec/hex/bits" begin
    @test string(UInt32('3'), base = 2) == "110011"
    @test string(UInt32('3'), pad = 7, base = 2) == "0110011"
    @test string(3, base = 2) == "11"
    @test string(3, pad = 2, base = 2) == "11"
    @test string(3, pad = Int32(2), base = Int32(2)) == "11"
    @test string(3, pad = typemin(Int128) + 3, base = 0x2) == "11"
    @test string(3, pad = 3, base = 2) == "011"
    @test string(-3, base = 2) == "-11"
    @test string(-3, pad = 3, base = 2) == "-011"

    @test string(9, base = 8) == "11"
    @test string(-9, base = 8) == "-11"
    @test string(-9, base = 8, pad = 5) == "-00011"
    @test string(-9, base = 8, pad = Int32(5)) == "-00011"

    @test string(121, base = 10) == "121"
    @test string(121, base = 10, pad = 5) == "00121"
    @test string(121, base = 10, pad = 5) == "00121"

    @test string(12, base = 16) == "c"
    @test string(-12, pad = 3, base = 16) == "-00c"
    @test string(-12, pad = Int32(3), base = Int32(16)) == "-00c"

    @test string(5, pad = 7, base = 2) == "0000101"

    @test bitstring(Int16(3)) == "0000000000000011"
    @test bitstring('3') == "00110011000000000000000000000000"
    @test bitstring(1035) == (Int == Int32 ? "00000000000000000000010000001011" :
        "0000000000000000000000000000000000000000000000000000010000001011")
    @test bitstring(Int128(3)) == "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011"
    @test bitstring(reinterpret(BitString128, Int128(3))) == "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000011"
end

@testset "digits/base" begin
    @test digits(5, base = 3) == [2, 1]
    @test digits(5, pad = 3) == [5, 0, 0]
    @test digits(5, pad = Int32(3)) == [5, 0, 0]
    # The following have bases powers of 2, but don't enter the fast path
    @test digits(-3, base = 2) == -[1, 1]
    @test digits(-42, base = 4) == -[2, 2, 2]

    @test_throws DomainError string(5, base = typemin(Int128) + 10)

    @testset "digits/base with bases powers of 2" begin
        @test digits(4, base = 2) == [0, 0, 1]
        @test digits(5, base = Int32(2), pad=Int32(3)) == [1, 0, 1]
        @test digits(42, base = 4) == [2, 2, 2]
        @test digits(321, base = 8) == [1, 0, 5]
        @test digits(0x123456789abcdef, base = 16) == 15:-1:1
        @test digits(0x2b1a210a750, base = 64) == [16, 29, 10, 4, 34, 6, 43]
        @test digits(0x02a01407, base = Int128(1024)) == [7, 5, 42]
    end

    @testset "digits/base with negative bases" begin
        @testset "digits(n::$T, base = b)" for T in (Int, UInt, BigInt, Int32, UInt32)
            @test digits(T(8163), base = -10) == [3, 4, 2, 2, 1]
            if !(T<:Unsigned)
                @test digits(T(-8163), base = -10) == [7, 7, 9, 9]
            end
            if T !== BigInt
                b = rand(-32:-2)
                for n = T[rand(T), typemax(T), typemin(T)]
                    # issue #29183
                    @test digits(n, base=b) == digits(signed(widen(n)), base=b)
                end
            end
        end
        @test [string(n, base = b)
               for n = [-10^9, -10^5, -2^20, -2^10, -100, -83, -50, -34, -27, -16, -7, -3, -2, -1,
                        0, 1, 2, 3, 4, 7, 16, 27, 34, 50, 83, 100, 2^10, 2^20, 10^5, 10^9]
               for b = [-2, -3, -7, -10, -60]] ==
                   ["11000101101001010100101000000000", "11211100201202120012",
                    "144246601121", "1000000000", "2hANlK", "111000111010100000",
                    "122011122112", "615462", "100000", "1XlK", "1100000000000000000000",
                    "11000202101022", "25055043", "19169584", "59Hi", "110000000000",
                    "12102002", "3005", "1036", "Iu", "11101100", "121112", "1515",
                    "1900", "2K", "11111101", "120011", "1651", "97", "2b", "11010010",
                    "2121", "1616", "50", "1A", "100010", "2202", "51", "46", "1Q",
                    "100101", "1000", "41", "33", "1X", "110000", "1102", "35", "24",
                    "1i", "1001", "1202", "10", "13", "1r", "1101", "10", "14", "17",
                    "1v", "10", "11", "15", "18", "1w", "11", "12", "16", "19", "1x", "0",
                    "0", "0", "0", "0", "1", "1", "1", "1", "1", "110", "2", "2", "2",
                    "2", "111", "120", "3", "3", "3", "100", "121", "4", "4", "4",
                    "11011", "111", "160", "7", "7", "10000", "211", "152", "196", "G",
                    "1101111", "12000", "146", "187", "R", "1100110", "12111", "136",
                    "174", "Y", "1110110", "11022", "101", "150", "o", "1010111", "10002",
                    "236", "123", "1xN", "110100100", "10201", "202", "100", "1xe",
                    "10000000000", "2211011", "14012", "19184", "1h4",
                    "100000000000000000000", "2001112212121", "162132144", "1052636",
                    "1uqiG", "1101001101111100000", "21002022201", "1103425", "1900000",
                    "SEe", "1001100111011111101111000000000", "120220201100111010001",
                    "44642116066", "19000000000", "1xIpcEe"]
    end
end

@testset "leading_ones and count_zeros" begin
    @test leading_ones(UInt32(Int64(2) ^ 32 - 2)) == 31
    @test leading_ones(1) == 0
    @test leading_zeros(Int32(1)) == 31
    @test leading_zeros(UInt32(Int64(2) ^ 32 - 2)) == 0

    @test count_zeros(Int64(1)) == 63
end

@testset "factorial" begin
    @test factorial(3) == 6
    @test factorial(Int8(3)) === 6
    @test_throws DomainError factorial(-3)
    @test_throws DomainError factorial(Int8(-3))
end

@testset "isqrt" begin
    @test isqrt(4) == 2
    @test isqrt(5) == 2
    @test isqrt(Int8(4)) === Int8(2)
    @test isqrt(Int8(5)) === Int8(2)
end

@testset "issue #4884" begin
    @test isqrt(9223372030926249000) == 3037000498
    @test isqrt(typemax(Int128)) == parse(Int128,"13043817825332782212")
    @test isqrt(Int128(typemax(Int64))^2-1) == 9223372036854775806
    @test isqrt(0) == 0
    for i = 1:1000
        n = rand(UInt128)
        s = isqrt(n)
        @test s*s <= n
        @test (s+1)*(s+1) > n
        n = rand(UInt64)
        s = isqrt(n)
        @test s*s <= n
        @test (s+1)*(s+1) > n
    end
end

# issue #9786
let ptr = Ptr{Cvoid}(typemax(UInt))
    for T in (Int, Cssize_t)
        @test T(ptr) == -1
        @test ptr == Ptr{Cvoid}(T(ptr))
        @test typeof(Ptr{Float64}(T(ptr))) == Ptr{Float64}
    end
end

# issue #15911
@inferred string(1)

# issue #22837
for b in [-100:-2; 2:100;]
    @test Base.ndigits0z(0, b) == 0
end

max_stars_repo_name: barucden/julia

max_issues_repo_name: KlausC/julia

max_forks_repo_name: KlausC/julia


text:
function derivative(fn)
    return function(x)
        # pick a small value for h
        h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
        # floating point arithmetic gymnastics
        xph = x + h
        dx = xph - x
        # evaluate f at x + h
        f1 = fn(xph)
        # evaluate f at x
        f0 = fn(x)
        # divide the difference by h
        # return (fn(xph) - fn(x)) / dx
        return (f1 - f0) / dx
    end

# 1st approximation:
#    return function(x)
#        h = 0.00001
#        f1 = fn(x + h)
#        f0 = fn(x)
#        return (f1 - f0)/h
#    end
end


text:
using LinearAlgebra

struct NullLattice <: Lattice end

mutable struct Branching
  k::Vector{Int}
  probs::Vector{Vector{Float64}}
  kMin::Int
  jMin::Int
  kMax::Int
  jMax::Int
end

function Branching()
  probs = Vector{Vector{Float64}}(undef, 3)
  probs[1] = zeros(0)
  probs[2] = zeros(0)
  probs[3] = zeros(0)
  return Branching(zeros(Int, 0), probs, typemax(Int), typemax(Int), typemin(Int), typemin(Int))
end

mutable struct TrinomialTree{S <: StochasticProcess} <: AbstractTree
  process::S
  timeGrid::TimeGrid
  dx::Vector{Float64}
  branchings::Vector{Branching}
  isPositive::Bool
end

function TrinomialTree(process::S, timeGrid::TimeGrid, isPositive::Bool = false) where {S <: StochasticProcess}
  x0 = process.x0
  dx = zeros(length(timeGrid.times))
  nTimeSteps = length(timeGrid.times) - 1
  jMin = 0
  jMax = 0
  branchings = Vector{Branching}(undef, nTimeSteps)

  for i = 1:nTimeSteps
    t = timeGrid.times[i]
    dt = timeGrid.dt[i]

    # Variance must be independent of x
    v2 = variance(process, t, 0.0, dt)
    v = sqrt(v2)
    dx[i+1] = v * sqrt(3.0)

    branching = Branching()

    @simd for j =jMin:jMax
      @inbounds x = x0 + j * dx[i]
      m = expectation(process, t, x, dt)
      @inbounds temp = round(Int, floor((m - x0) / dx[i+1] + 0.5))

      if isPositive
        @inbounds while (x0 + (temp - 1) * dx[i + 1] <= 0)
          temp += 1
        end
      end

      @inbounds e = m - (x0 + temp * dx[i + 1])
      e2 = e * e
      e3 = e * sqrt(3.0)

      p1 = (1.0 + e2 / v2 - e3 / v) / 6.0
      p2 = (2.0 - e2 / v2) / 3.0
      p3 = (1.0 + e2 / v2 + e3 / v) / 6.0

      add!(branching, temp, p1, p2, p3)
    end

    @inbounds branchings[i] = branching # check if we need copy

    jMin = branching.jMin
    jMax = branching.jMax
  end

  return TrinomialTree{S}(process, timeGrid, dx, branchings, isPositive)
end

function rebuild_tree!(tt::TrinomialTree, timeGrid::TimeGrid)
  x0 = tt.process.x0
  dx = zeros(length(timeGrid.times))
  nTimeSteps = length(timeGrid.times) - 1
  jMin = 0
  jMax = 0
  branchings = Vector{Branching}(undef, nTimeSteps)

  for i = 1:nTimeSteps
    t = timeGrid.times[i]
    dt = timeGrid.dt[i]

    # Variance must be independent of x
    v2 = variance(tt.process, t, 0.0, dt)
    v = sqrt(v2)
    dx[i+1] = v * sqrt(3.0)

    branching = Branching()

    @simd for j =jMin:jMax
      @inbounds x = x0 + j * dx[i]
      m = expectation(tt.process, t, x, dt)
      @inbounds temp = round(Int, floor((m - x0) / dx[i+1] + 0.5))

      if tt.isPositive
        @inbounds while (x0 + (temp - 1) * dx[i + 1] <= 0)
          temp += 1
        end
      end

      @inbounds e = m - (x0 + temp * dx[i + 1])
      e2 = e * e
      e3 = e * sqrt(3.0)

      p1 = (1.0 + e2 / v2 - e3 / v) / 6.0
      p2 = (2.0 - e2 / v2) / 3.0
      p3 = (1.0 + e2 / v2 + e3 / v) / 6.0

      add!(branching, temp, p1, p2, p3)
    end

    @inbounds branchings[i] = branching # check if we need copy

    jMin = branching.jMin
    jMax = branching.jMax
  end
  tt.dx = dx
  tt.timeGrid = timeGrid
  tt.branchings = branchings

  return tt
end

mutable struct TreeLattice1D{T <: TreeLattice} <: TreeLattice
  tg::TimeGrid
  impl::T
  statePrices::Vector{Vector{Float64}}
  n::Int
  statePricesLimit::Int
end

function TreeLattice1D(tg::TimeGrid, n::Int, impl::T) where {T <: TreeLattice}
  statePrices = Vector{Vector{Float64}}(undef, 1)
  statePrices[1] = ones(1)

  statePricesLimit = 1

  return TreeLattice1D{T}(tg, impl, statePrices, n, statePricesLimit)
end

function get_grid(lat::TreeLattice1D, t::Float64)
  i = return_index(lat.tg, t)
  grid = zeros(get_size(lat.impl, i))

  @simd for j in eachindex(grid)
    @inbounds grid[j] = get_underlying(lat.impl, i, j)
  end

  return grid
end

mutable struct TreeLattice2D{T <: TreeLattice} <: TreeLattice
  tg::TimeGrid
  impl::T
  statePrices::Vector{Vector{Float64}}
  n::Int
  statePricesLimit::Int
  tree1::TrinomialTree
  tree2::TrinomialTree
  m::Matrix{Float64}
  rho::Float64
end

function TreeLattice2D(tree1::TrinomialTree, tree2::TrinomialTree, correlation::Float64, impl::T) where {T <: TreeLattice}
  tg = tree1.timeGrid
  statePrices = Vector{Vector{Float64}}(undef, 1)
  statePrices[1] = ones(1)

  statePricesLimit = 1
  branch_num = branches(TrinomialTree)
  n = branch_num ^ 2
  m = zeros(branch_num, branch_num)
  rho = abs(correlation)

  if correlation < 0.0 && branch_num == 3
    m[1,1] = -1.0
    m[2,1] = -4.0
    m[3,1] =  5.0
    m[1,2] = -4.0
    m[2,2] =  8.0
    m[3,2] = -4.0
    m[1,3] =  5.0
    m[2,3] = -4.0
    m[3,3] = -1.0
  else
    m[1,1] =  5.0
    m[2,1] = -4.0
    m[3,1] = -1.0
    m[1,2] = -4.0
    m[2,2] =  8.0
    m[3,2] = -4.0
    m[1,3] = -1.0
    m[2,3] = -4.0
    m[3,3] =  5.0
  end

  return TreeLattice2D{T}(tg, impl, statePrices, n, statePricesLimit, tree1, tree2, m, rho)
end

get_size(tr::TreeLattice2D, i::Int) = get_size(tr.tree1, i) * get_size(tr.tree2, i)

function descendant(tr::TreeLattice2D, i::Int, idx::Int, branch::Int)
  modulo = get_size(tr.tree1, i)
  new_idx = idx - 1
  new_branch = branch - 1

  index1 = (new_idx % modulo) + 1
  index2 = (round(Int, floor(new_idx / modulo))) + 1

  branch1 = (new_branch % branches(TrinomialTree)) + 1
  branch2 = (round(Int, floor(new_branch / branches(TrinomialTree)))) + 1

  modulo = get_size(tr.tree1, i + 1)
  return ((descendant(tr.tree1, i, index1, branch1) - 1) + (descendant(tr.tree2, i, index2, branch2) - 1) * modulo) + 1
end

function probability(tr::TreeLattice2D, i::Int, idx::Int, branch::Int)
  modulo = get_size(tr.tree1, i)
  new_idx = idx - 1
  new_branch = branch - 1

  index1 = (new_idx % modulo) + 1
  index2 = (round(Int, floor(new_idx / modulo))) + 1

  branch1 = (new_branch % branches(TrinomialTree)) + 1
  branch2 = (round(Int, floor(new_branch / branches(TrinomialTree)))) + 1

  # println("$(index1) $(index2) $(branch1) $(branch2) $(idx) $(branch) $(modulo)")

  prob1 = probability(tr.tree1, i, index1, branch1)
  prob2 = probability(tr.tree2, i, index2, branch2)

  # println("$(prob1) $(prob2)")

  return prob1 * prob2 + tr.rho * tr.m[branch2,branch1] / 36.0 # this 36 could depend on the branches(TrinomialTree)
end

function get_state_prices!(t::TreeLattice, i::Int)
  if i > t.statePricesLimit
    compute_state_prices!(t, i)
  end

  return t.statePrices[i]
end

function compute_state_prices!(t::TreeLattice, until::Int)
  @inbounds @simd for i = t.statePricesLimit:until - 1
    push!(t.statePrices, zeros(get_size(t.impl, i + 1)))
    for j = 1:get_size(t.impl, i)
      disc = discount(t.impl, i, j)
      statePrice = t.statePrices[i][j]
      for l = 1:t.n
        t.statePrices[i + 1][descendant(t.impl, i, j, l)] += statePrice * disc * probability(t.impl, i, j, l)
      end
    end
  end

  t.statePricesLimit = until

  return t
end

function initialize!(lattice::TreeLattice, asset::DiscretizedAsset, t::Float64)
  #i = findfirst(lattice.tg.times .>= t)
  i = return_index(lattice.tg, t)
  set_time!(asset, t)
  reset!(asset, get_size(lattice.impl, i))
end

function rollback!(lattice::TreeLattice, asset::DiscretizedAsset, t::Float64)
  partial_rollback!(lattice, asset, t)
  adjust_values!(asset)

  return asset
end

function partial_rollback!(lattice::TreeLattice, asset::DiscretizedAsset, t::Float64)
  from = asset.common.time

  if QuantLib.Math.is_close(from, t)
    return
  end

  # iFrom = findfirst(lattice.tg.times .>= from)
  # iTo = findfirst(lattice.tg.times .>= t)
  iFrom = return_index(lattice.tg, from)
  iTo = return_index(lattice.tg, t)

  @simd for i = iFrom-1:-1:iTo
    newVals = zeros(get_size(lattice.impl, i))
    step_back!(lattice, i, asset.common.values, newVals)
    @inbounds asset.common.time = lattice.tg.times[i]
    asset.common.values = newVals
    # println("newVals before: ", newVals)
    if i != iTo
      adjust_values!(asset)
    end
    # println("newVals after: ", newVals)
  end

  return asset
end

function present_value(lattice::TreeLattice, asset::DiscretizedAsset)
  i = findfirst(lattice.tg.times .>= asset.common.time)
  return dot(asset.common.values, get_state_prices!(lattice, i))
end

function step_back!(lattice::TreeLattice, i::Int, vals::Vector{Float64}, newVals::Vector{Float64})
  # pragma omp parallel for
  for j = 1:get_size(lattice.impl, i)
    val = 0.0
    @simd for l = 1:lattice.n
      @inbounds val += probability(lattice.impl, i, j, l) * vals[descendant(lattice.impl, i, j, l)]
    end
    @inbounds val *= discount(lattice.impl, i, j)
    @inbounds newVals[j] = val
  end

  return newVals
end

get_size(b::Branching) = b.jMax - b.jMin + 1

function add!(branch::Branching, k::Int, p1::Float64, p2::Float64, p3::Float64)
  push!(branch.k, k)
  push!(branch.probs[1], p1)
  push!(branch.probs[2], p2)
  push!(branch.probs[3], p3)

  # maintain invariants
  branch.kMin = min(branch.kMin, k)
  branch.jMin = branch.kMin - 1
  branch.kMax = max(branch.kMax, k)
  branch.jMax = branch.kMax + 1

  return branch
end

descendant(b::Branching, idx::Int, branch::Int) = b.k[idx] - b.jMin - 1 + branch
probability(b::Branching, idx::Int, branch::Int) = b.probs[branch][idx]
branches(::Type{TrinomialTree}) = 3
get_size(t::TrinomialTree, i::Int) = i == 1 ? 1 : get_size(t.branchings[i-1])

function get_underlying(t::TrinomialTree, i::Int, idx::Int)
  if i == 1
    return t.process.x0
  else
    @inbounds return t.process.x0 + (t.branchings[i - 1].jMin + (idx - 1)) * t.dx[i]
  end
end

descendant(t::TrinomialTree, i::Int, idx::Int, branch::Int) = descendant(t.branchings[i], idx, branch)

probability(t::TrinomialTree, i::Int, j::Int, b::Int) = probability(t.branchings[i], j, b)


text:
using DifferentialEquations
using ProgressMeter
using Printf
using MATLAB
include("escape_exit_function.jl")

max_hit=5
t_end=1e3
n_iter_P=1000
n_iter_Q=2001
width=.7
height=1
ArrP=range(-width,stop=width,length=n_iter_P)
ArrQ=range(-height,stop=height,length=n_iter_Q)
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/Escape/"
H=range(.21,stop=0.24,length=21)

@showprogress 1 "Computing..." for Energy in H

println(Energy)

Q_0=zeros(0)
P_0=zeros(0)

Q_1=zeros(0)
P_1=zeros(0)

Q_2=zeros(0)
P_2=zeros(0)

Q_3=zeros(0)
P_3=zeros(0)

Q_4=zeros(0)
P_4=zeros(0)

Q_5=zeros(0)
P_5=zeros(0)
#
# @showprogress 1 "Computing..." for Energy in H

h=replace(@sprintf("%.13f",Energy),"."=>"_")
file_name=location*"Escape_"*h*".fig"
for j=1:n_iter_P
    for k=1:n_iter_Q
        K=escape_exit_function(ArrQ[k],ArrP[j], Energy, t_end, max_hit)
        if K==0
            push!(Q_0,ArrQ[k])
            push!(P_0,ArrP[j])
        elseif K==1
            push!(Q_1,ArrQ[k])
            push!(P_1,ArrP[j])
        elseif K==2
            push!(Q_2,ArrQ[k])
            push!(P_2,ArrP[j])
        elseif K==3
            push!(Q_3,ArrQ[k])
            push!(P_3,ArrP[j])
        elseif K==4
            push!(Q_4,ArrQ[k])
            push!(P_4,ArrP[j])
        elseif K==5
            push!(Q_5,ArrQ[k])
            push!(P_5,ArrP[j])
        # elseif K==-1
        #     push!(Q_bound,ArrQ[k])
        #     push!(P_bound,ArrP[j])
        end
    end
end

size=3
#
# mat"axis([ -$width,$width,-$height,$height ])"
mat"figure(); hold on;"
mat"axis([ -$width,$width,-$height,$height ])"
# mat"plot($P_bound,$Q_bound ,'k.','MarkerSize',10)"
mat"plot($P_0,$Q_0 ,'b.','MarkerSize',$size)"
mat"plot($P_1,$Q_1 ,'r.','MarkerSize',$size)"
mat"plot($P_2,$Q_2 ,'g.','MarkerSize',$size)"
mat"plot($P_4,$Q_4 ,'c.','MarkerSize',$size)"
mat"plot($P_5,$Q_5 ,'y.','MarkerSize',$size)"
mat"savefig($file_name)"
mat"close"
 end


text:
using FiniteDifferences: Forward, Backward, Central, Nonstandard, add_tiny

@testset "Methods" begin

    @testset "add_tiny" begin
        @test add_tiny(convert(Float64, 5)) isa Float64
        @test add_tiny(convert(Float32, 5)) isa Float32
        @test add_tiny(convert(Float16, 5)) isa Float16

        @test add_tiny(convert(Int, 5)) isa Float64
        @test add_tiny(convert(UInt, 5)) isa Float64
        @test add_tiny(convert(Bool, 1)) isa Float64
    end

    # The different approaches to approximating the gradient to try.
    methods = [forward_fdm, backward_fdm, central_fdm]

    # The different floating point types to try and the associated required relative
    # tolerance.
    types = [Float32, Float64]

    # The different functions to evaluate (.f), their first derivative at 1 (.d1),
    # and second derivative at 1 (.d2).
    foos = [
        (f=sin, d1=cos(1), d2=-sin(1)),
        (f=exp, d1=exp(1), d2=exp(1)),
        (f=abs2, d1=2, d2=2),
        (f=x -> sqrt(x + 1), d1=0.5 / sqrt(2), d2=-0.25 / 2^(3/2)),
    ]

    # Test all combinations of the above settings, i.e. differentiate all functions using
    # all methods and data types.
    @testset "foo=$(foo.f), method=$m, type=$(T)" for foo in foos, m in methods, T in types

        @testset "method-order=$order" for order in [1, 2, 3]
            @test m(order, 0; adapt=2)(foo.f, T(1)) isa T
            @test m(order, 0; adapt=2)(foo.f, T(1)) == T(foo.f(1))
        end

        @test m(10, 1)(foo.f, T(1)) isa T
        @test m(10, 1)(foo.f, T(1)) ≈ T(foo.d1)

        @test m(10, 2; bound=1)(foo.f, T(1)) isa T
        if T == Float64
            @test m(10, 2; bound=1)(foo.f, T(1)) ≈ T(foo.d2)
        end
    end

    # Integration test to ensure that Integer-output functions can be tested.
    @testset "Integer Output" begin
        @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12)
    end

    @testset "Adaptation improves estimate" begin
        @test forward_fdm(5, 1)(log, 0.001; adapt=0) ≈ 969.2571703
        @test forward_fdm(5, 1)(log, 0.001; adapt=1) ≈ 1000
    end

    @testset "Limiting step size" begin
        @test !isfinite(central_fdm(5, 1)(abs, 0.001; max_step=0))
        @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0
    end

    @testset "Accuracy at high orders, with high adapt" begin
        # Regression test against issues with precision during computation of _coeffs
        # see https://github.com/JuliaDiff/FiniteDifferences.jl/issues/64

        @test fdm(central_fdm(9, 5), exp, 1.0, adapt=4) ≈ exp(1) atol=1e-7

        poly(x) = 4x^3 + 3x^2 + 2x + 1
        @test fdm(central_fdm(9, 3), poly, 1.0, adapt=4) ≈ 24 atol=1e-11
    end

    @testset "Printing FiniteDifferenceMethods" begin
        @test sprint(show, central_fdm(2, 1)) == """
            FiniteDifferenceMethod:
              order of method:       2
              order of derivative:   1
              grid:                  [-1, 1]
              coefficients:          [-0.5, 0.5]
            """
        m, _ = fdm(central_fdm(2, 1), sin, 1, Val(true))
        report = sprint(show, m)
        regex_float = r"[\d\.\+-e]+"
        regex_array = r"\[([\d.+-e]+(, )?)+\]"
        @test occursin(Regex(join(map(x -> x.pattern,
            [
                r"FiniteDifferenceMethod:",
                r"order of method:", r"\d+",
                r"order of derivative:", r"\d+",
                r"grid:", regex_array,
                r"coefficients:", regex_array,
                r"roundoff error:", regex_float,
                r"bounds on derivatives:", regex_float,
                r"step size:", regex_float,
                r"accuracy:", regex_float,
                r""
            ]
        ), r"\s*".pattern)), report)
    end

    @testset "Breaking deprecations" begin
        @test_throws ErrorException fdm([1,2,3], 4)  # Custom grids need Nonstandard
        for f in (forward_fdm, backward_fdm, central_fdm)
            @test_throws ErrorException f(2, 1; M=1)  # Old kwarg, now misplaced
            @test_throws ErrorException f(2, 1, Val(true))  # Ask fdm for reports instead
        end
    end

    @testset "Types" begin
        @testset "$T" for T in (Forward, Backward, Central)
            @test T(5, 1)(sin, 1; adapt=4) ≈ cos(1)
            @test_throws ArgumentError T(3, 3)
            @test_throws ArgumentError T(3, 4)
            @test_throws ArgumentError T(40, 5)
            @test_throws ArgumentError T(5, 1)(sin, 1; adapt=200)
            @test_throws ArgumentError T(5, 1)(sin, 1; eps=0.0)
            @test_throws ArgumentError T(5, 1)(sin, 1; bound=0.0)
        end
        @testset "Nonstandard" begin
            @test Nonstandard([-2, -1, 1], 1)(sin, 1) ≈ cos(1)
            @test_throws ArgumentError Nonstandard([-2, -1, 1], 1)(sin, 1; adapt=2)
        end
    end
end


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export serve

#------------------------------------------------------------------------------
# External Server API
#------------------------------------------------------------------------------

function serve (portno::Integer)
    info("[server] Listening on port $portno")
    srv = listen(portno)
    @async accept_loop(srv)
end

#------------------------------------------------------------------------------
# Internal Functions
#------------------------------------------------------------------------------

function accept_loop (srv)
    while true
        s = accept(srv)
        info("[server] Accepted connection, handling request")
        @async handle_request(s)
    end
end

max_issues_repo_name: esl/FsServer

max_forks_repo_name: esl/FsServer


text:
# Base Quality
# ============
#
# A representation of positions-specific integer quality scores, as in FASTQ.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"A `QualityEncoding` value holds a set of encoding compatible with a string"
bitstype 16 QualityEncoding

Base.convert(::Type{QualityEncoding}, nt::UInt16) = reinterpret(QualityEncoding, nt)
Base.convert(::Type{UInt16}, nt::QualityEncoding) = reinterpret(UInt16, nt)

@compat function Base.:|(a::QualityEncoding, b::QualityEncoding)
    return convert(QualityEncoding, convert(UInt16, a) | convert(UInt16, b))
end

@compat function Base.:&(a::QualityEncoding, b::QualityEncoding)
    return convert(QualityEncoding, convert(UInt16, a) & convert(UInt16, b))
end

"`QualityEncoding` indicating no compatible encodings."
const EMPTY_QUAL_ENCODING      = convert(QualityEncoding, UInt16(0))

"Sanger (Phred+33) quality score encoding"
const SANGER_QUAL_ENCODING     = convert(QualityEncoding, UInt16(0b00001))

"Solexa (Solexa+64) quality score encoding"
const SOLEXA_QUAL_ENCODING     = convert(QualityEncoding, UInt16(0b00010))

"Illumina 1.3 (Phred+64) quality score encoding"
const ILLUMINA13_QUAL_ENCODING = convert(QualityEncoding, UInt16(0b00100))

"Illumina 1.5 (Phred+64) quality score encoding"
const ILLUMINA15_QUAL_ENCODING = convert(QualityEncoding, UInt16(0b01000))

"Illumina 1.8 (Phred+33) quality score encoding"
const ILLUMINA18_QUAL_ENCODING = convert(QualityEncoding, UInt16(0b10000))

"`QualityEncoding` indicating all known encodings are compatible."
const ALL_QUAL_ENCODINGS =
    SANGER_QUAL_ENCODING | SOLEXA_QUAL_ENCODING | ILLUMINA13_QUAL_ENCODING |
    ILLUMINA15_QUAL_ENCODING | ILLUMINA18_QUAL_ENCODING

# Index into this with `trailing_zeros(encoding) + 1`
const qual_encoding_ranges = [
    (UInt8('!'), UInt8('~'), Int8(0) ),  # SANGER
    (UInt8(';'), UInt8('~'), Int8(-5)),  # SOLEXA
    (UInt8('@'), UInt8('~'), Int8(0) ),  # ILLUMINA13
    (UInt8('B'), UInt8('~'), Int8(2) ),  # ILLUMINA15
    (UInt8('!'), UInt8('~'), Int8(0) )   # ILLUMINA18
]

# Build an encoding lookup table
const compatible_qual_encoding = fill(EMPTY_QUAL_ENCODING, length('!':'~'))

for (encoding_num, (first, last, start_score)) in enumerate(qual_encoding_ranges)
    encoding = convert(QualityEncoding, UInt16(1) << (encoding_num - 1))
    for c in first:last
        compatible_qual_encoding[c - UInt8('!') + 1] |= encoding
    end
end

# When a quality string has multiple compatible encodings, we choose the first
# compatible alphabet in this list.
const preferred_quality_encodings = [
    ILLUMINA15_QUAL_ENCODING,
    ILLUMINA13_QUAL_ENCODING,
    SOLEXA_QUAL_ENCODING,
    SANGER_QUAL_ENCODING,
    ILLUMINA18_QUAL_ENCODING,
]

"""
Infer the encoding of a FASTQ quality string.

# Arguments
  * `data`: sequence data in a string
  * `start`: first position to consider in data
  * `stop`: last position to consider in data
  * `encodings`: valid encodings
  * `default`: if there are multiple compatible encodings, default
    to this one if it's compatible.

# Returns
A pair of `QualityEncoding`s. The first given the chosen encoding, the second
giving the set of compatible encodings.
"""
function infer_quality_encoding(data::Vector{UInt8}, start, stop,
                                encodings::QualityEncoding=ALL_QUAL_ENCODINGS,
                                default::QualityEncoding=EMPTY_QUAL_ENCODING)
    encodings = ALL_QUAL_ENCODINGS
    @inbounds for i in start:stop
        c = data[i]
        if '!' <= c <= '~'
            encodings &= compatible_qual_encoding[c - UInt8('!') + 1]
        else
            error("Character $(convert(Char, c)) is not compatible with any known quality encoding.")
        end

        if count_ones(convert(UInt16, encodings)) <= 0
            break
        end
    end

    if count_ones(convert(UInt16, encodings)) == 0
        error("String is not compatible with any known sequence type.")
    elseif default != EMPTY_QUAL_ENCODING && (encodings & default) != EMPTY_ALPHABET
        return (default, encodings)
    else
        for encoding in preferred_quality_encodings
            if encoding & encodings != EMPTY_QUAL_ENCODING
                return (encoding, encodings)
            end
        end
    end
    return (default, encodings)
end

"""
Decode a quality string in place into integer Phred scores.

# Arguments:
  * `encoding::QualityEncoding`: how the quality scores are encoded.
  * `input::Vector{UInt8}`: character data to decode
  * `output:::Vector{UInt8}`: Phred score output vector, assumed to be
    of length `stop - start + 1`.
  * `start`: First position in `input` to decode.
  * `stop`: Last position with in `input` to decode.

# Returns
`output`
"""
function decode_quality_string!(encoding::QualityEncoding, input::Vector{UInt8},
                                output::Vector{Int8}, start=1,
                                stop=min(length(output), length(input)))
    if length(output) != stop - start + 1
        resize!(output, stop - start + 1)
    end

    encoding_num = trailing_zeros(convert(UInt16, encoding)) + 1
    first, last, startqual = qual_encoding_ranges[encoding_num]

    fill!(output, startqual - first)
    for i in 1:(stop - start + 1)
        @inbounds output[i] += input[start + i - 1]
    end

    return output
end

"""
Decode a quality string in place into integer Phred scores.

# Arguments:
  * `encoding::QualityEncoding`: how the quality scores are encoded.
  * `input::Vector{UInt8}`: character data to decode
  * `start`: First position in `input` to decode.
  * `stop`: Last position with in `input` to decode.

# Returns
A `Vector{UInt8}` of length `stop - start + 1` containing integer Phred scores.
"""
function decode_quality_string(encoding::QualityEncoding, input::Vector{UInt8},
                               start=1, stop=length(input))
    output = Array(Int8, stop - start + 1)
    return decode_quality_string!(encoding, input, output, start, stop)
end

"""
Encode a Phred quality score vector in place into a quality string.

# Arguments:
  * `encoding::QualityEncoding`: how the quality scores are encoded.
  * `input::Vector{UInt8}`: character data to decode
  * `output:::Vector{UInt8}`: Phred score output vector, assumed to be
    of length `stop - start + 1`.
  * `start`: First position in `input` to decode.
  * `stop`: Last position with in `input` to decode.

# Returns
`output`
"""
function encode_quality_string!(encoding::QualityEncoding, input::Vector{Int8},
                                output::Vector{UInt8}, start=1,
                                stop=min(length(output), length(input)))
    @inbounds begin
        if length(output) != stop - start + 1
            resize!(output, stop - start + 1)
        end

        encoding_num = trailing_zeros(convert(UInt16, encoding)) + 1
        first, last, startqual = qual_encoding_ranges[encoding_num]
        for (i, j) in enumerate(start:stop)
            c = input[j]
            output[i] =  (c + first) - startqual
        end
        return output
    end
end

"""
Encode a vector of integer Phred scores into a quality string

# Arguments:
  * `encoding::QualityEncoding`: how the quality scores are encoded.
  * `input::Vector{UInt8}`: character data to decode
  * `start`: First position in `input` to decode.
  * `stop`: Last position with in `input` to decode.

# Returns
A `Vector{UInt8}` of length `stop - start + 1` containing ASCII-encoded phred
scores
"""
function encode_quality_string(encoding::QualityEncoding, input::Vector{Int8},
                               start=1, stop=length(input))
    output = Array(UInt8, stop - start + 1)
    return encode_quality_string!(encoding, input, output, start, stop)
end


text:
function mouseconnectivitycache()
    mcc = mouse_connectivity_cache.MouseConnectivityCache(manifest_file=mouseconnectivitymanifest)
end
export MouseConnectivityCache

function getannotationvolume()
    vol, info = mouseconnectivitycache().get_annotation_volume()
end

function getstructuremask(structureid::Number)
    strucutreid = Int(structureid)
    mask, info =  mouseconnectivitycache().get_structure_mask(structureid)
end

function gettemplatevolume()
    template, template_info =  mouseconnectivitycache().get_template_volume()
end

function getreferencespace()
    reference, reference_info = mouseconnectivitycache().get_reference_space()
end


text:
struct Mlrda
    fm
    lev::AbstractVector
    ni::AbstractVector
end

"""
    mlrda(X, y, weights = ones(size(X, 1)))
Discrimination based on multple linear regression (MLR-DA).
* `X` : X-data.
* `y` : Univariate class membership.
* `weights` : Weights of the observations.

The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a multiple linear regression (MLR) is run between the X-data and and each column
of Ydummy, returning predictions of the dummy variables (object `posterior`
retunerd by fuction `predict`).
These predictions can be  considered as unbounded (i.e. eventuall outside of [0, 1])
estimates of the class memberships.
For a given observation, the final class prediction is the class corresponding
to the dummy variable for which the prediction is the highest.
"""
function mlrda(X, y, weights = ones(size(X, 1)))
    z = dummy(y)
    fm = mlr(X, z.Y, weights)
    Mlrda(fm, z.lev, z.ni)
end

function predict(object::Mlrda, X)
    X = ensure_mat(X)
    m = size(X, 1)
    zp = predict(object.fm, X).pred
    z =  mapslices(argmax, zp; dims = 2)
    pred = reshape(replacebylev(z, object.lev), m, 1)
    (pred = pred, posterior = zp)
end


text:
# Julia wrapper for header: /usr/include/scip/cons_setppc.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrSetppc(scip)
    ccall((:SCIPincludeConshdlrSetppc, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateConsSetpart(scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSetpart, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSetpart(scip, cons, name, nvars, vars)
    ccall((:SCIPcreateConsBasicSetpart, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}), scip, cons, name, nvars, vars)
end

function SCIPcreateConsSetpack(scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSetpack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSetpack(scip, cons, name, nvars, vars)
    ccall((:SCIPcreateConsBasicSetpack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}), scip, cons, name, nvars, vars)
end

function SCIPcreateConsSetcover(scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSetcover, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, nvars, vars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSetcover(scip, cons, name, nvars, vars)
    ccall((:SCIPcreateConsBasicSetcover, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Cint, Ptr{Ptr{SCIP_VAR}}), scip, cons, name, nvars, vars)
end

function SCIPaddCoefSetppc(scip, cons, var)
    ccall((:SCIPaddCoefSetppc, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{SCIP_VAR}), scip, cons, var)
end

function SCIPgetNVarsSetppc(scip, cons)
    ccall((:SCIPgetNVarsSetppc, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetVarsSetppc(scip, cons)
    ccall((:SCIPgetVarsSetppc, libscip), Ptr{Ptr{SCIP_VAR}}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetTypeSetppc(scip, cons)
    ccall((:SCIPgetTypeSetppc, libscip), SCIP_SETPPCTYPE, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualsolSetppc(scip, cons)
    ccall((:SCIPgetDualsolSetppc, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetDualfarkasSetppc(scip, cons)
    ccall((:SCIPgetDualfarkasSetppc, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetRowSetppc(scip, cons)
    ccall((:SCIPgetRowSetppc, libscip), Ptr{SCIP_ROW}, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetNFixedonesSetppc(scip, cons)
    ccall((:SCIPgetNFixedonesSetppc, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end

function SCIPgetNFixedzerosSetppc(scip, cons)
    ccall((:SCIPgetNFixedzerosSetppc, libscip), Cint, (Ptr{SCIP_}, Ptr{SCIP_CONS}), scip, cons)
end


text:
using SimLib
using XXZNumerics
using Test
using Distributed
using Random

to_file(func, path; mode="w") = open(path, mode) do out
    redirect_stdout(out) do
        func()
    end
end

log_contains(path, str) = any(l -> contains(l, str), eachline(path))

const PREFIX = tempname()
mkdir(PREFIX)
@show PREFIX

Random.seed!(5)

@show pwd()

@testset "SimLib.jl" begin
    # Write your tests here.
    include("positions.jl")
    include("ed.jl") # depends on generating positions
    include("ensembles.jl") # should be after ed.jl to reuse data
end


text:
# LSAC JUNE 2007 Form 8LSN75
# https://www.lsac.org/sites/default/files/legacy/docs/default-source/jd-docs/sampleptjune.pdf

# Questions 1-5

#A company employee generates a series of five-digit product
#codes in accordance with the following rules:
#The codes use the digits 0, 1, 2, 3, and 4, and no others.

julia> using AbstractLogic
Start the repl in command prompt by typing `=`.

# unique forces each value to be different
a, b, c, d, e  ∈  0:4 || unique [clear]

#The second digit has a value exactly twice that of the first digit.
b == a*2

#The value of the third digit is less than the value of the fifth digit
c < e

#1. If the last digit of an acceptable product code is 1...
# (which of the following has must be true N out of N cases).

#(A) first digit is 2
prove: e == 1 ==> a==2
# Returns true: 6 of 6 possible.

#(B) second digit is 0
prove: e == 1 ==> b==0
# False: 5 out of 6

#(C) third digit is 3
prove: e == 1 ==> c==3
# False: 5 out of 6

#(D) fourth digit is 4
prove: d == 1 ==> d==4
# False: 5 out of 6

#(E) fourth digit is 0
prove: d == 1 ==> d==0
# False: 5 out of 6

#2. Which one of the following must be true about any
#acceptable product code?
range

#(A) The digit 1 appears in some position before the digit 2.
prove: {{i}} == 1 &&& {{>i}} == 2 {{1}}
# False: only 3 of 6

#(B) The digit 1 appears in some position before the digit 3.
prove: {{i}} == 1 &&& {{>i}} == 3 {{1}}
# False: only 5 of 6

#(C) The digit 2 appears in some position before the digit 3.
prove: {{i}} == 2 &&& {{>i}} == 3 {{1}}
# True! 6 of 6

#(D) The digit 3 appears in some position before the digit 0.
prove: {{i}} == 0 &&& {{>i}} == 3 {{1}}
# False: only 5 of 6

#(E) The digit 4 appears in some position before the digit 3.
prove: {{i}} == 4 &&& {{>i}} == 3 {{1}}
# False: only 4 of 6

#3. If the third digit of an acceptable product code is not 0, which one of the
# following must be true?

#(A) The second digit of the product code is 2.
prove: c != 0 ==> b == 2
# False: only 5 of 6

#(B) The third digit of the product code is 3.
prove: c != 0 ==> b == 2
# False: only 5 of 6

#(C) The fourth digit of the product code is 0.
prove: c != 0 ==> d == 0
# True! 6 of 6

#(D) The fifth digit of the product code is 3. <true>
prove: c != 0 ==> e == 3
# False: 5 of 6

#(E) The fifth digit of the product code is 1.
prove: c != 0 ==> e == 1
# False: 4 of 6

#4. Any of the following pairs could be the third and fourth digits,
# respectively, of an acceptable product
#code, EXCEPT:

# (A) 0, 1
any: c == 0 &&& d == 1
# Possible 1 of 6

# (B) 0, 3
any: c == 0 &&& d == 3
# Possible 2 of 6

#(C) 1, 0
any: c == 1 &&& d == 0
# Possible 1 of 6

#(D) 3, 0
any: c == 3 &&& d == 0
# Possible 1 of 6

#(E) 3, 4
any: c == 3 &&& d == 4
# False! 0 of 6

#5. Which one of the following must be true about any acceptable product code?
show

#(A) There is exactly one digit between the digit 0 and the digit 1.
prove: {{i}} == 1|0 ==> {{i+2}} == 1|0
# False 1 of 6 outcomes

#(B) There is exactly one digit between the digit 1 and the digit 2.
prove: {{i}} == 1|2 ==> {{i+2}} == 1|2
# False 0 of 6 outcomes

#(C) There are at most two digits between the digit 1 and the digit 3.
prove: {{i}} == 1|3 ==> {{i+4}} != 1,3
# False 5 out of 6 outcomes true

#(D) There are at most two digits between the digit 2 and the digit 3.
prove: {{i}} == 2|3 ==> {{i+4}} != 2,3
# False 4 out of 6 outcomes feasible

#(E) There are at most two digits between the digit 2 and the digit 4.
prove: {{i}} == 2|4 ==> {{i+4}} != 2,4
# True 6 out of 6 outcomes feasible