text:
using CompositeGrids: Grid
using StaticArrays

#
# generate a p grid for integration. dense on k and kf.
#

struct PGrid{T<:AbstractFloat, SIZE}
    k::T
    kF::T
    iscombined::Int
    grid::MVector{SIZE,T}
    size::Int
    head::T
    tail::T
    unilogs::SVector{4,Grid.UniLog{T}}

    function PGrid{T,SIZE}(k, kF, maxk, mink, MS) where {T<:AbstractFloat, SIZE}
        size = sum(MS)+1

        grid, segment,segindex = [],[],[]
        unilogs = []
        iscombined = 0

        if abs(k-kF)<2mink
            k1 = (k+kF)/2.0
            iscombined = 2
            M1 = MS[1]+MS[2]-1
            M2 = MS[3]+MS[4]-1

            g1 = Grid.UniLog{T}([0.0, k1], 1, mink, M1, 1, false, [false, true])
            push!(unilogs,g1)
            push!(unilogs,g1)
            for idx=g1.idx[1]:g1.idx[2]
                push!(grid, Grid._grid(g1, idx))
            end

            g2 = Grid.UniLog{T}([k1, maxk], M1+1, mink, M2, 1, true, [false, false])
            push!(unilogs,g2)
            push!(unilogs,g2)
            for idx=g2.idx[1]:g2.idx[2]
                push!(grid, Grid._grid(g2, idx))
            end

        elseif abs(k)<2mink
            iscombined = 1

            g0 = Grid.UniLog{T}([0.0, kF/2.0], 1, mink, MS[1]+MS[2]-1, 1, true, [false, true])
            push!(unilogs,g0)
            push!(unilogs,g0)
            for idx=g0.idx[1]:g0.idx[2]
                push!(grid, Grid._grid(g0, idx))
            end

            g1 = Grid.UniLog{T}([kF/2.0, kF], MS[1]+MS[2]+1, mink, MS[3]-1, 1, false, [false, true])
            push!(unilogs,g1)
            for idx=g1.idx[1]:g1.idx[2]
                push!(grid, Grid._grid(g1, idx))
            end

            g2 = Grid.UniLog{T}([kF, maxk], MS[1]+MS[2]+MS[3]+1, mink, MS[4]-1, 1, true, [false, false])
            push!(unilogs,g2)
            for idx=g2.idx[1]:g2.idx[2]
                push!(grid, Grid._grid(g2, idx))
            end

        else
            k1, k2 = minimum([k,kF]), maximum([k,kF])
            bounds = [[0.0, k1],[k1,(k1+k2)/2.0], [(k1+k2)/2.0, k2],[k2, maxk]]
            d2s = [false, true, false, true]
            isopen = [[false, true],[false, true],[false, true],[false, false]]
            init = 1

            for i in 1:4
                g = Grid.UniLog{T}(bounds[i], init, mink, MS[i]-1, 1, d2s[i], isopen[i])
                push!(unilogs, g)
                for idx=g.idx[1]:g.idx[2]
                    push!(grid, Grid._grid(g, idx))
                end
                init += MS[i]
            end

        end

        head, tail = grid[1], grid[end]
        return new{T,size}(k, kF, iscombined, grid, size, head, tail, unilogs)
    end

end

"""
   function pGrid(k, kF, maxk, mink, MS)

Create a logarithmic p Grid, which is densest near the Fermi momentum ``k_F`` and a given momentum ``k``

#Arguments:
- k: a given momentum
- kF: Fermi momentum
- maxK: the upper bound of the grid
- mink: the minimum interval of the grid
- MS: a list of # of points in each segments, the total size will be sum(MS)+1
"""
@inline function pGrid(k, kF, maxk, mink, MS)
    return PGrid{Float64, sum(MS)+1}(k, kF, maxk, mink, MS)
end

if abspath(PROGRAM_FILE) == @__FILE__
    pgrid = pGrid(0.001, 1.0, 10.0, 0.001, [4,4,4,4])
    println(pgrid.grid)
    pgrid = pGrid(0.999, 1.0, 10.0, 0.001, [4,4,4,4])
    println(pgrid.grid)
    pgrid = pGrid(2.0, 1.0, 10.0, 0.001, [4,4,4,4])
    println(pgrid.grid)
end


text:
immutable Cubic{BC<:Flag} <: Degree{3} end
Cubic{BC<:Flag}(::BC) = Cubic{BC}()

"""
Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
implemented here is defined as follows.

    y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))

where

    p(δx) = 1/6 * (1-δx)^3
    q(δx) = 2/3 - δx^2 + 1/2 δx^3

and the values `cX` for `X ∈ {m, _, p, pp}` are the pre-filtered coefficients.

For future reference, this expands out to the following polynomial:

    y_i(x) = 1/6 cm (1+i-x)^3 + c (2/3 - (x-i)^2 + 1/2 (x-i)^3) +
             cp (2/3 - (1+i-x)^2 + 1/2 (1+i-x)^3) + 1/6 cpp (x-i)^3

When we derive boundary conditions we will use derivatives `y_0'(x)` and
`y_0''(x)`
"""
Cubic

"""
`define_indices_d` for a cubic b-spline calculates `ix_d = floor(x_d)` and
`fx_d = x_d - ix_d` (corresponding to `i` `and `δx` in the docstring for
`Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`
"""
function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
    symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
    symixpp, symx, symfx = Symbol("ixpp_",d), Symbol("x_",d), Symbol("fx_",d)
    quote
        # ensure that all of ix_d, ixm_d, ixp_d, and ixpp_d are in-bounds no
        # matter the value of pad
        $symix = clamp(floor(Int, $symx), $(2-pad), size(itp,$d)+$(pad-2))
        $symfx = $symx - $symix
        $symix += $pad # padding for oob coefficient
        $symixm = $symix - 1
        $symixp = $symix + 1
        $symixpp = $symixp + 1
    end
end

"""
`define_indices_d`  for a cubic, periodic b-spline calculates `ix_d = floor(x_d)`
and `fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring entry
for `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`.

If any `ixX_d` for `x ∈ {m, p, pp}` (note: not `c_d`) should fall outside of
the data interval, they wrap around.
"""
function define_indices_d(::Type{BSpline{Cubic{Periodic}}}, d, pad)
    symix, symixm, symixp = Symbol("ix_",d), Symbol("ixm_",d), Symbol("ixp_",d)
    symixpp, symx, symfx = Symbol("ixpp_",d), Symbol("x_",d), Symbol("fx_",d)
    quote
        $symix = clamp(floor(Int, $symx), 1, size(itp,$d))
        $symfx = $symx - $symix
        $symixm = mod1($symix - 1, size(itp,$d))
        $symixp = mod1($symix + 1, size(itp,$d))
        $symixpp = mod1($symix + 2, size(itp,$d))
    end
end

padding{BC<:Flag}(::Type{BSpline{Cubic{BC}}}) = Val{1}()
padding(::Type{BSpline{Cubic{Periodic}}}) = Val{0}()

"""
In `coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that

    cm_d  = p(fx_d)
    c_d   = q(fx_d)
    cp_d  = q(1-fx_d)
    cpp_d = p(1-fx_d)

where `p` and `q` are defined in the docstring entry for `Cubic`, and
`fx_d` in the docstring entry for `define_indices_d`.
"""
function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
    symm, sym =  Symbol("cm_",d), Symbol("c_",d)
    symp, sympp = Symbol("cp_",d) ,Symbol("cpp_",d)
    symfx = Symbol("fx_",d)
    symfx_cub = Symbol("fx_cub_", d)
    sym_1m_fx_cub = Symbol("one_m_fx_cub_", d)
    quote
        $symfx_cub = cub($symfx)
        $sym_1m_fx_cub = cub(1-$symfx)
        $symm = SimpleRatio(1,6)*$sym_1m_fx_cub
        $sym  = SimpleRatio(2,3) - sqr($symfx) + SimpleRatio(1,2)*$symfx_cub
        $symp = SimpleRatio(2,3) - sqr(1-$symfx) + SimpleRatio(1,2)*$sym_1m_fx_cub
        $sympp = SimpleRatio(1,6)*$symfx_cub
    end
end

"""
In `gradient_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that

    cm_d  = p'(fx_d)
    c_d   = q'(fx_d)
    cp_d  = q'(1-fx_d)
    cpp_d = p'(1-fx_d)

where `p` and `q` are defined in the docstring entry for `Cubic`, and
`fx_d` in the docstring entry for `define_indices_d`.
"""
function gradient_coefficients{C<:Cubic}(::Type{BSpline{C}}, d)
    symm, sym, symp, sympp = Symbol("cm_",d), Symbol("c_",d), Symbol("cp_",d), Symbol("cpp_",d)
    symfx = Symbol("fx_",d)
    symfx_sqr = Symbol("fx_sqr_", d)
    sym_1m_fx_sqr = Symbol("one_m_fx_sqr_", d)
    quote
        $symfx_sqr = sqr($symfx)
        $sym_1m_fx_sqr = sqr(1 - $symfx)

        $symm  = -SimpleRatio(1,2) * $sym_1m_fx_sqr
        $sym   =  SimpleRatio(3,2) * $symfx_sqr     - 2 * $symfx
        $symp  = -SimpleRatio(3,2) * $sym_1m_fx_sqr + 2 * (1 - $symfx)
        $sympp =  SimpleRatio(1,2) * $symfx_sqr
    end
end

"""
In `hessian_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that

    cm_d  = p''(fx_d)
    c_d   = q''(fx_d)
    cp_d  = q''(1-fx_d)
    cpp_d = p''(1-fx_d)

where `p` and `q` are defined in the docstring entry for `Cubic`, and
`fx_d` in the docstring entry for `define_indices_d`.
"""
function hessian_coefficients{C<:Cubic}(::Type{BSpline{C}}, d)
    symm, sym, symp, sympp = Symbol("cm_",d), Symbol("c_",d), Symbol("cp_",d), Symbol("cpp_",d)
    symfx = Symbol("fx_",d)
    quote
        $symm  = 1 - $symfx
        $sym   = 3 * $symfx - 2
        $symp  = 1 - 3 * $symfx
        $sympp = $symfx
    end
end

function index_gen{C<:Cubic,IT<:DimSpec{BSpline}}(::Type{BSpline{C}}, ::Type{IT}, N::Integer, offsets...)
    if length(offsets) < N
        d = length(offsets)+1
        symm, sym, symp, sympp =  Symbol("cm_",d), Symbol("c_",d), Symbol("cp_",d), Symbol("cpp_",d)
        return :($symm * $(index_gen(IT, N, offsets...,-1)) + $sym * $(index_gen(IT, N, offsets..., 0)) +
                 $symp * $(index_gen(IT, N, offsets..., 1)) + $sympp * $(index_gen(IT, N, offsets..., 2)))
    else
        indices = [offsetsym(offsets[d], d) for d = 1:N]
        return :(itp.coefs[$(indices...)])
    end
end

# ------------ #
# Prefiltering #
# ------------ #

"""
`Cubic`: continuity in function value, first and second derivatives yields

    2/3 1/6
    1/6 2/3 1/6
        1/6 2/3 1/6
           ⋱  ⋱   ⋱
"""
function inner_system_diags{T,C<:Cubic}(::Type{T}, n::Int, ::Type{C})
    du = fill(convert(T, SimpleRatio(1, 6)), n-1)
    d = fill(convert(T, SimpleRatio(2, 3)), n)
    dl = copy(du)
    dl, d, du
end

"""
`Cubic{Flat}` `OnGrid` amounts to setting `y_1'(x) = 0` at `x = 1`.
Applying this condition yields

    -cm + cp = 0
"""
function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                   ::Type{Cubic{Flat}}, ::Type{OnGrid})
    dl, d, du = inner_system_diags(T, n, Cubic{Flat})
    d[1] = d[end] = -one(T)
    du[1] = dl[end] = zero(T)

    # Now Woodbury correction to set `[1, 3], [n, n-2] ==> 1`
    specs = WoodburyMatrices.sparse_factors(T, n, (1, 3, one(T)), (n, n-2, one(T)))

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

"""
`Cubic{Flat}`, `OnCell` amounts to setting `y_1'(x) = 0` at `x = 1/2`.
Applying this condition yields

    -9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0

or, equivalently,

    -9 cm + 11 c -3 cp + 1 cpp = 0

(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
"""
function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                   ::Type{Cubic{Flat}}, ::Type{OnCell})
    dl, d, du = inner_system_diags(T,n,Cubic{Flat})
    d[1] = d[end] = -9
    du[1] = dl[end] = 11

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> -3
    #    - [1, 4] and [n, n-3] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, T(-3)),
                                  (n, n-2, T(-3)),
                                  (1, 4, one(T)),
                                  (n, n-3, one(T))
                                  )

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

"""
`Cubic{Line}` `OnCell` amounts to setting `y_1''(x) = 0` at `x = 1/2`.
Applying this condition yields

    3 cm -7 c + 5 cp -1 cpp = 0

(Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
"""
function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                   ::Type{Cubic{Line}}, ::Type{OnCell})
    dl,d,du = inner_system_diags(T,n,Cubic{Line})
    d[1] = d[end] = 3
    du[1] = dl[end] = -7

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> -3
    #    - [1, 4] and [n, n-3] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, T(5)),
                                  (n, n-2, T(5)),
                                  (1, 4, -one(T)),
                                  (n, n-3, -one(T))
                                  )

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

"""
`Cubic{Line}` `OnGrid` amounts to setting `y_1''(x) = 0` at `x = 1`. Applying this
condition gives:

    1 cm -2 c + 1 cp = 0
"""
function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
                                   ::Type{Cubic{Line}}, ::Type{OnGrid})
    dl,d,du = inner_system_diags(T,n,Cubic{Line})
    d[1] = d[end] = 1
    du[1] = dl[end] = -2

    # now need Woodbury correction to set :
    #    - [1, 3] and [n, n-2] ==> 1
    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, one(T)),
                                  (n, n-2, one(T)),
                                  )

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

"""
`Cubic{Periodic}` `OnGrid` closes the system by looking at the coefficients themselves
as periodic, yielding

    c0 = c(N+1)

where `N` is the number of data points.
"""
function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                ::Type{Cubic{Periodic}}, ::Type{GT})
    dl, d, du = inner_system_diags(T,n,Cubic{Periodic})

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, n, du[1]),
                                  (n, 1, dl[end])
                                  )

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

"""
`Cubic{Free}` `OnGrid` and `Cubic{Free}` `OnCell` amount to requiring an extra
continuous derivative at the second-to-last cell boundary; this means
`y_1'''(2) = y_2'''(2)`, yielding

    1 cm -3 c + 3 cp -1 cpp = 0
"""
function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
                                                ::Type{Cubic{Free}}, ::Type{GT})
    dl, d, du = inner_system_diags(T,n,Cubic{Periodic})

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, n, du[1]),
                                  (n, 1, dl[end])
                                  )

    Woodbury(lufact!(Tridiagonal(dl, d, du), Val{false}), specs...), zeros(TC, n)
end

ext: jl


text:
include("./util.jl")

using Test: @test

function line2nums(line::AbstractString)::Array{Int8,1}
    [
        if c == '.'
            0
        elseif c == '>'
            1
        elseif c == 'v'
            2
        else
            throw("Invalid character $c")
        end
        for c in line
    ]
end

parse(s::Vector)::Array{Int8,2} =
    transpose(hcat(map(line2nums, s)...))

nums = parse(util.input_lines(25))

@enum Direction east south

function next_idx(idx::CartesianIndex, direction::Direction, max::Int64)::CartesianIndex
    if direction == east
        next = idx[2] + 1
        if next > max
            next = 1
        end
        CartesianIndex(idx[1], next)
    else
        next = idx[1] + 1
        if next > max
            next = 1
        end
        CartesianIndex(next, idx[2])
    end
end

function step(curr::Array{Int8,2})::Bool
    east_shift = hcat(curr[:, 2:end], curr[:, 1:1])
    east_movers = findall((east_shift .== 0) .* (curr .== 1))
    east_max = size(curr)[2]
    for idx in east_movers
        target = next_idx(idx, east, east_max)
        curr[idx] = 0
        curr[target] = 1
    end

    south_shift = vcat(curr[2:end, :], curr[1:1, :])
    south_movers = findall((south_shift .== 0) .* (curr .== 2))
    south_max = size(curr)[1]
    for idx in south_movers
        target = next_idx(idx, south, south_max)
        curr[idx] = 0
        curr[target] = 2
    end
    length(south_movers) > 0 || length(east_movers) > 0
end

count = 1

test_input = split("""v...>>.vv>
.vv>>.vv..
>>.>v>...v
>>v>>.>.v.
v>v.vv.v..
>.>>..v...
.vv..>.>v.
v.v..>>v.v
....v..v.>""", '\n')

test_nums = parse(test_input)

test_count = 1

while step(test_nums)
    global test_count += 1
end

@test test_count == 58

n_p1 = copy(nums)

while step(n_p1)
    global count += 1
end

println("Part 1: $count")


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This code was ported from a MATLAB implementation available at
# http://exolete.com/code/lbm/
# with permission of the author.
# Original MATLAB implementation is copyright (c) 2006 Iain Haslam.

using ParallelAccelerator
using DocOpt

function computeF(density, c_squ, u_squ, t, u, v)
    t * density * (1 + u / c_squ + 0.5 * (v / c_squ) * (v / c_squ) - u_squ / (2 * c_squ))
end

@acc function boltzmann_kernel(numactivenodes, F, G, ON, UX, X1)
    avu = 1.
    prevavu = 1.
    ts = 0
    deltaU = 1e-7
    t1 = 4/9
    t2 = 1/9
    t3 = 1/36
    c_squ = 1/3
    omega = 1.0
    while (ts<4000 && 1e-10<abs((prevavu-avu)/avu)) || ts<100
        ux_sum = 0.
        runStencil(F[1], F[2], F[3], F[4], F[5], F[6], F[7], F[8], F[9],
                   G[1], G[2], G[3], G[4], G[5], G[6], G[7], G[8], G[9], ON, UX, X1,
                   1, :oob_wraparound) do F1, F2, F3, F4, F5, F6, F7, F8, F9, G1, G2, G3, G4, G5, G6, G7, G8, G9, ON, UX, X1
                       # Propagate
                       f4 = F4[1, -1]
                       f3 = F3[0,-1]
                       f2 = F2[-1,-1]
                       f5 = F5[1,0]
                       f1 = F1[-1,0]
                       f6 = F6[1,1]
                       f7 = F7[0,1]
                       f8 = F8[-1,1]
                       f9 = F9[0,0]
                       density = f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9
                       ux = (f1 + f2 + f8 - f4 - f5 - f6) / density
                       uy = (f2 + f3 + f4 - f6 - f7 - f8) / density
                       if X1[0,0] == 1
                           ux = ux + deltaU # Increase inlet pressure
                       end
                       if ON[0,0] == 1
                           ux = 0.
                           uy = 0.
                           density = 0.
                       end
                       UX[0,0] = ux
                       u_squ = ux * ux + uy * uy
                       u_c2  = ux + uy
                       u_c4  = uy - ux
                       u_c6  = - u_c2
                       u_c8  = - u_c4
                       # Calculate equilibrium distribution: stationary
                       feq9 = t1 * density * (1. - u_squ / (2. * c_squ))
                       # nearest-neighbours
                       feq1 = computeF(density, c_squ, u_squ, t2, ux, ux)
                       feq3 = computeF(density, c_squ, u_squ, t2, uy, uy)
                       feq5 = computeF(density, c_squ, u_squ, t2, -ux, ux)
                       feq7 = computeF(density, c_squ, u_squ, t2, -uy, uy)
                       # next-nearest neighbours
                       feq2 = computeF(density, c_squ, u_squ, t3, u_c2, u_c2)
                       feq4 = computeF(density, c_squ, u_squ, t3, u_c4, u_c4)
                       feq6 = computeF(density, c_squ, u_squ, t3, u_c6, u_c6)
                       feq8 = computeF(density, c_squ, u_squ, t3, u_c8, u_c8)
                       if ON[0,0] == 1
                           G5[0,0] = f1
                           G6[0,0] = f2
                           G7[0,0] = f3
                           G8[0,0] = f4
                           G1[0,0] = f5
                           G2[0,0] = f6
                           G3[0,0] = f7
                           G4[0,0] = f8
                       else
                           G1[0,0] = omega * feq1 + (1. - omega) * f1
                           G2[0,0] = omega * feq2 + (1. - omega) * f2
                           G3[0,0] = omega * feq3 + (1. - omega) * f3
                           G4[0,0] = omega * feq4 + (1. - omega) * f4
                           G5[0,0] = omega * feq5 + (1. - omega) * f5
                           G6[0,0] = omega * feq6 + (1. - omega) * f6
                           G7[0,0] = omega * feq7 + (1. - omega) * f7
                           G8[0,0] = omega * feq8 + (1. - omega) * f8
                       end
                       G9[0,0] = omega * feq9 + (1. - omega) * f9
                       return
                   end
        X = G
        G = F
        F = X
        prevavu = avu
        avu = sum(UX) / numactivenodes
        ts = ts + 1
    end
    return F
end

function boltzmann(nx, ny)
    density = 1.0
    F  = [ fill!(zeros(nx, ny), density/9) for i = 1:9 ]
    G  = deepcopy(F)
    X1 = [ x == 1 ? 1 : 0 for x = 1:nx, y = 1:ny ]
    UX = zeros(nx, ny)
    ON = zeros(Int, nx, ny)
    ON[rand(nx, ny) .> 0.7] = 1 # extremely porous random domain
    numactivenodes = sum(1 .- ON)
    boltzmann_kernel(numactivenodes, F, G, ON, UX, X1)
end

function main()
    doc = """boltzmann.jl

2D lattice Boltzmann fluid flow model.

Usage:
  boltzmann.jl -h | --help
  boltzmann.jl [--nx=<nx>] [--ny=<ny>]

Options:
  -h --help  Show this screen.
  --nx=<nx>  Specify size of simulation grid in the x dimension [default: 200].
  --ny=<ny>  Specify size of simulation grid in the y dimension [default: 200].
"""
    arguments = docopt(doc)
    nx = parse(Int, arguments["--nx"])
    ny = parse(Int, arguments["--ny"])

    srand(0)

    println("grid size: ", nx, " x ", ny)

    tic()
    F = boltzmann(1, 1)
    println("SELFPRIMED ", toq())

    tic()
    F = boltzmann(nx, ny)
    println("SELFTIMED ", toq())

    println("checksum: ", sum(F[1]), " ", sum(F[2]), " ", sum(F[3]), " ", sum(F[4]), " ", sum(F[5]), " ", sum(F[6]), " ", sum(F[7]), " ", sum(F[8]), " ", sum(F[9]))
end

main()


text:
using Test
using QCMaterial

@testset "core.fermionopertor" begin
    op1 = FermionOperator("1^")
    op2 = FermionOperator("1")

    @test op1 == op1
    @test op1 != op2
    @test op1 * op2 == FermionOperator("1^ 1")

    @test FermionOperator("", 0.0) == FermionOperator()
end

@testset "core.of_qubit_operator" begin
    op = OFQubitOperator("X1 X6", 1.0) + OFQubitOperator("Y2 Y3", 2.0)
    @test get_n_qubit(op) == 6
    @test get_term_count(op) == 2
    @test terms_dict(op)[((1,"X"), (6,"X"))] == 1.0
    @test terms_dict(op)[((2,"Y"), (3,"Y"))] == 2.0
end

@testset "core.qulacs_quantum_circuit" begin
    n_qubit = 2
    state = QulacsQuantumState(n_qubit, 0b00)
    # Initial state |00>
    #  (1) H gate to 1-th qubit
    #  (2) CNOT gate (control 1, target 2)
    #  (3) X gate to 1-th qubit
    # Final state (|01> + |10>)/sqrt(2)
    circuit = QulacsQuantumCircuit(n_qubit)
    add_H_gate!(circuit, 1)
    add_CNOT_gate!(circuit, 1, 2)
    add_X_gate!(circuit, 1)
    update_quantum_state!(circuit, state)

    @test get_vector(state) ≈ [0, 1, 1, 0]/sqrt(2)
end

@testset "core.qulacs_quantum_state" begin
    n_qubit = 2

    # Initialize to 0b0
    state = QulacsQuantumState(n_qubit)
    state2 = QulacsQuantumState(n_qubit, 0b0)
    @test get_vector(state) == get_vector(state2)

    # Initialize to 0b01
    state = QulacsQuantumState(n_qubit)
    set_computational_basis!(state, 0b01)
    state2 = QulacsQuantumState(n_qubit, 0b01)
    @test get_vector(state) == get_vector(state2)

    # Test inner_product
    state = QulacsQuantumState(n_qubit, 0b01)
    @test inner_product(state, state) ≈ 1.0
end

@testset "core.add_parametric_circuit_using_generator" begin
    n_qubit = 4
    a, i = 1, 1
    crr, ann = 1, 0
    @assert 1 <= 2*a <= n_qubit && 1 <= 2*i <= n_qubit
    generator = jordan_wigner(FermionOperator([(a, crr), (i, ann)], 1.0))

    circuit = UCCQuantumCircuit(n_qubit)
    add_parametric_circuit_using_generator!(circuit, generator, 0.0)
    add_parametric_circuit_using_generator!(circuit, generator, 1.0)
    @test circuit.theta_offsets == [
        (2, 0, [0.0]),
        (2, 2, [0.0])
    ]
end

@testset "core.qubit_operator" begin
    op1 = OFQubitOperator("Z1", 1.0)
    op2 = OFQubitOperator("Z1", 1.0im)
    op3 = OFQubitOperator("X1", 1.0)

    state0 = create_hf_state(1, 0) # Create |0>
    state1 = create_hf_state(1, 1) # Create |1>

    @test is_hermitian(op1)
    @test !is_hermitian(op2)
    @test get_expectation_value(op1, state1) == -1.0
    @test get_transition_amplitude(op3, state1, state0) == 1.0

    @test OFQubitOperator("X1") * OFQubitOperator("Y1") == OFQubitOperator("Z1", 1.0im)
    @test OFQubitOperator("X1 Y2") * OFQubitOperator("Y1 Z2") == OFQubitOperator("Z1 X2", -1.0)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from Euler angles to direction cosine matrices.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/angle_to_dcm.jl
# =======================================

# Functions: angle_to_dcm
# -----------------------

@testset "Euler angles => DCM (Float64)" begin
    T = Float64

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to DCM.
        D = angle_to_dcm(ea)
        @test eltype(D) === T

        # Create the DCM by composing the rotations.
        rot_seq_str = string(rot_seq)
        D₁ = create_rotation_matrix(ea.a1, Symbol(rot_seq_str[1]))
        D₂ = create_rotation_matrix(ea.a2, Symbol(rot_seq_str[2]))
        D₃ = create_rotation_matrix(ea.a3, Symbol(rot_seq_str[3]))
        De = D₃ * D₂ * D₁

        # Compare.
        @test D ≈ De
    end
end

@testset "Euler angles => DCM (Float32)" begin
    T = Float32

    for rot_seq in valid_rot_seqs
        # Sample Euler angles.
        ea = EulerAngles(_rand_ang(T), _rand_ang(T), _rand_ang(T), rot_seq)

        # Convert to DCM.
        D = angle_to_dcm(ea)
        @test eltype(D) === T

        # Create the DCM by composing the rotations.
        rot_seq_str = string(rot_seq)
        D₁ = create_rotation_matrix(ea.a1, Symbol(rot_seq_str[1]))
        D₂ = create_rotation_matrix(ea.a2, Symbol(rot_seq_str[2]))
        D₃ = create_rotation_matrix(ea.a3, Symbol(rot_seq_str[3]))
        De = D₃ * D₂ * D₁

        # Compare.
        @test D ≈ De
    end
end

@testset "Euler angles => DCM (Errors)" begin
    @test_throws ArgumentError angle_to_dcm(1, 2, 3, :ZZY)
end

# Functions: smallangle_to_dcm
# ----------------------------

@testset "Small Euler angles => DCM (Float64)" begin
    T = Float64

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Create the rotation matrix without normalization.
    D = smallangle_to_rot(θx, θy, θz, normalize = false)
    @test eltype(D) === T

    # Expected result.
    De = DCM(
          1, +θz, -θy,
        -θz,   1, +θx,
        +θy, -θx,   1
    )'

    @test De ≈ D

    # Create the rotation matrix with normalization.
    D = smallangle_to_rot(θx, θy, θz, normalize = true)
    @test eltype(D) === T

    @test orthonormalize(De) ≈ D
end

@testset "Small Euler angles => DCM (Float32)" begin
    T = Float32

    # Sample three small angles.
    θx = _rand_ang(T) * T(0.001)
    θy = _rand_ang(T) * T(0.001)
    θz = _rand_ang(T) * T(0.001)

    # Create the rotation matrix without normalization.
    D = smallangle_to_rot(θx, θy, θz, normalize = false)
    @test eltype(D) === T

    # Expected result.
    De = DCM(
          1, +θz, -θy,
        -θz,   1, +θx,
        +θy, -θx,   1
    )'

    @test De ≈ D

    # Create the rotation matrix with normalization.
    D = smallangle_to_rot(θx, θy, θz, normalize = true)
    @test eltype(D) === T

    @test orthonormalize(De) ≈ D
end


text:
module ExpectationStubs
using MacroTools
export @stub, @expect, Stub, all_expectations_used, @used, @usecount

include("donotcare.jl")
include("vectordict.jl")

struct SyntaxError <: Exception
end

"""
   ExpectationValueMismatchError

Similar to a julia MethodError
This is thrown if a call was made on a stub,
but the values/types did not match those that were expected.
"""
struct ExpectationValueMismatchError <: Exception
    name
    argvals
end

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

struct Stub{name} <: Function
    expectations::VectorDict{Any, Any}
    calls::Vector{Any}
end
Stub(name)=Stub{name}(VectorDict{Any, Any}(), Vector{Any}())

function(stub::Stub{name})(args...) where {name}
    if !haskey(stub.expectations, args)
        throw(ExpectationValueMismatchError(name, args))
    end
    push!(stub.calls, args)

    stub.expectations[args]
end

"""
    usecount(stub:Stub, arg)

Returns how many times the `stub` has been called with a matching argument
"""
function usecount(stub::Stub, arg)
    net = 0
    for call_record in stub.calls
        if call_record == arg
            net+=1
        end
    end
    net
end

"""
    @stub(name)

Declares that you will be making a stub function called `name`.
The key difference between a stub and a julia function,
is that stubs can have values for arguements declared,
which must match, as well as having types which must match.
Also they are simpler, in that they can not varying their return result based on argument.

Once you have declared a stub,
you should use `@expect` to declare what it should respond to.

Calling methods that do not exist (with values provided), will result in Errors.
This is intentional, as the stub exists to check that your function is only called with the arguements that you say are valid.
"""
macro stub(name)
    quote
        $(esc(name))=Stub($(QuoteNode(name)))
    end
end

"""
    onlyesc(v)

Like `esc`, except it ignores things that can not be escaped.
Eg literals.

Internal use
"""
onlyesc(v::Any)=v
onlyesc(v::Union{Symbol,Expr})=esc(v)

"""
    split_vals_and_sig(argsexpr)

Takes an expression from a function definition's args
eg `:([a::Int, ::Integer])`
and breaks it down into the values, and the types.
When value is not given it subsitutes `DoNotCare()`,
when type is not given it subsitutes `typeof(val)`
Returns a tuple of each, with all things in expression/symbol form.

Internal use
"""
function split_vals_and_sig(argsexpr)
    vals=Expr(:tuple)
    sig=Expr(:tuple)
    for term in argsexpr
        if @capture(term, v_::s_)
            push!(vals.args, onlyesc(v))
            push!(sig.args, s)
        elseif @capture(term, ::s_)
            push!(vals.args, :(DoNotCare{$(esc(s))}()))
            push!(sig.args, s)
        elseif @capture(term, v_)
            push!(vals.args, onlyesc(v))
            push!(sig.args, :(typeof($(onlyesc(v)))))
        else
            throw(SyntaxError())
        end
    end
    vals, sig
end

"""
    @expect(defn)

Prepares a stub, lets it know to expect a function call matching the `defn`,
and to return the result.

Has several forms, which align to julia function declations.
Examples:

- `@expect(foo(1, 1.5)=10)` this prepares the stub `foo`, to return 10, if it gets given the input `(1, 1.5)`
- `@expect(foo(3, ::Int)=20)` this prepares the stub `foo`, to return 20, if it gets given the input with first arguement 3, and second argument any `Int`
- `@expect(foo(5, ::Any)=30)` this prepares the stub `foo`, to return 30, if it gets given the input with first arguement 5, and second argument any type
- `@expect(foo(Any, ::Any)=40)` this prepares the stub `foo`, to return 40, if it is given 2 arguements

Notes that you can not for the same stub both the declare that it is bound by type (i.e. that it can take anything of a given type),
and declare that it is bound by value, for the same parameter.

Note that the result can not depend on the arguments of the function.
This is intentional, as it is there to keep your Stubs simple and to the point.
So you don't end up needing to test your tests.

Currently does not support KWArgs.
"""
macro expect(defn)
    @capture(defn, name_(args__)=ret_) || throw(SyntaxError())
    argvals, sig = split_vals_and_sig(args)
    quote
        ###########################################################
        # Check not already regeistered
        if haskey($(esc(name)).expectations, $(argvals))
            @warn "Expectation already set" name=$(Meta.quot(name)) argvals=$(argvals) sig=$(sig)
        end

        # actually register the value
        $(esc(name)).expectations[$(argvals)] = $(esc(ret))
    end |> unblock |>  MacroTools.striplines
end

"""
    all_expectations_used(stub::Stub)

Checks that every expectation setup for the stub was actually used.
It is good to have this as a sanity check at the end of your test script using the stub.
"""
function all_expectations_used(stub::Stub)
    all(usecount.(Ref(stub), keys(stub.expectations)) .> 0)
end

"""
    @used(defn)

Returns true if a particular expectionation was used.
Syntax is similar to `@expect`

Like `@expect` the key can be a value or a type.
It does not have to match to the one that was used in `@expect`

Normally one would used this inside a test:

- `@test @used(foo(1, 1.5))` test that `foo` was called with `(1, 1.5)`
- `@test @used(foo(3, ::Int))` test the `foo` was called with 3 and some Int
- `@test !@used(foo(5, ::Any))` test that `foo` was never called with the first arg 5 (and exactly 2 args)
- `@test !@used(foo(Any, ::Any))` test tjat `foo` was never alled with 2 args
"""
macro used(defn)
    name, argvals, sig = name_argvals_and_sig(defn)
    :(usecount($(esc(name)), $(argvals)) > 0)
end

"""
    @usecount(defn)

Like `@used` but returns the count of how many times the particular expectation
was called.

Syntax is the same.
"""
macro usecount(defn)
    name, argvals, sig = name_argvals_and_sig(defn)
   :(usecount($(esc(name)), $(argvals)))
end

"""
    name_argvals_and_sig(defn)
Does the extraction of these as is needed for `@used` and `@usedcounts`
"""
function name_argvals_and_sig(defn)
    @capture(defn, name_(args__)) || throw(SyntaxError())
    argvals, sig = split_vals_and_sig(args)
    name, argvals, sig
end

end #Module


text:
using XLA, IRTools, Test
using IRTools: IR, argument!, xcall
using XLA: XTuple, GetTupleElement, Sub, Mul, Gt, While, compile

ir = IR()
x = argument!(ir, Int)
n = argument!(ir, Int)
xnr = push!(ir, xcall(XTuple(), x, n, 1))
cond = let ir = IR()
  xnr = argument!(ir)
  n = push!(ir, xcall(GetTupleElement(1), xnr))
  push!(ir, xcall(Gt(), n, 0))
  ir
end
cond = push!(ir, cond)
body = let ir = IR()
  xnr = argument!(ir)
  x = push!(ir, xcall(GetTupleElement(0), xnr))
  n = push!(ir, xcall(GetTupleElement(1), xnr))
  r = push!(ir, xcall(GetTupleElement(2), xnr))
  n = push!(ir, xcall(Sub(), n, 1))
  r = push!(ir, xcall(Mul(), x, r))
  xnr = push!(ir, xcall(XTuple(), x, n, r))
  ir
end
body = push!(ir, body)
xnr = push!(ir, xcall(While(), cond, body, xnr))
push!(ir, xcall(GetTupleElement(2), xnr))
pow = compile(ir)

@test pow(2, 3)[1] == 2^3
@test pow(5, 10)[1] == 5^10


text:
mutable struct RegularizedSVD{T} <: MatrixFactorization{T}
    P::Array
    Q::Array
    preference::Persa.Preference{T}
    users::Int
    items::Int
end

const RSVD = RegularizedSVD

function RegularizedSVD(dataset::Persa.Dataset, features::Int)
    (users, items) = size(dataset)

    P = rand(users, features)
    Q = rand(items, features)

    return RegularizedSVD(P, Q, dataset.preference, Persa.users(dataset), Persa.items(dataset))
end

Persa.predict(model::RegularizedSVD, user::Int, item::Int) = model.P[user, :]' * model.Q[item, :]

function objective(model::RegularizedSVD, dataset::Persa.Dataset, λ::Float64)
    total = 0

    for (u, v, r) in dataset
        total += (r - model[u, v])^2
        total += λ * (norm(model.P[u,:])^2 + norm(model.Q[v,:])^2)
    end

    return total
end

function update!(model::RegularizedSVD, dataset::Persa.Dataset, γ::Float64, λ::Float64)

    idx = shuffle(1:length(dataset))

    for i = 1:length(dataset)
        (u, v, r) = dataset[idx[i]]

        e = r - Persa.predict(model, u, v)

        P = model.P[u,:]
        Q = model.Q[v,:]

        model.P[u,:] += γ * (e .* Q .- λ .* P)
        model.Q[v,:] += γ * (e .* P .- λ .* Q)
    end
end


text:
"""
The X Server reported a connection error.
"""
struct ConnectionError <: Exception
    msg
    code
end

"""
A request to the X Server returned with an error.
"""
struct RequestError <: Exception
    msg
end

"""
Error when flushing a connection to the X Server.
"""
struct FlushError <: Exception
    code
end


text:
module ComradeBase

using DocStringExtensions
using ChainRulesCore
using FLoops
using StaticArrays
using StructArrays

export  visibility, intensitymap, intensitymap!,
        StokesVector, CoherencyMatrix, evpa, m̆, SingleStokes,
        flux, fov, imagepixels, pixelsizes, IntensityMap

include("interface.jl")
include("images/images.jl")
#include("threaded.jl")

#include("polarized.jl")

end


text:
function edit_line(
        line, direction_restriction::Vec2f0, clampto, window;
        knob_scale = 9f0,
        knob_color = RGBA{Float32}(0.7, 0.7, 0.7, 1.0),
        kw_args...
    )
    mouse_hover = mouse2id(window)
    inds = reinterpret(Cuint, collect(IterTools.partition(Cuint(0):Cuint(length(line)-1) , 2, 1)))
    line_robj = visualize(
        line, :linesegment;
        indices = inds,
        kw_args...
    ).children[]
    point_gpu = line_robj[:vertex]
    points = visualize(
        (Circle{Float32}(Point2f0(0), knob_scale), point_gpu);
        color = knob_color,
        kw_args...
    )
    point_robj = points.children[]
    gpu_position = point_robj[:position]
    gpu_color = point_robj[:color]
    gpu_scale = point_robj[:scale]

    m2id = mouse2id(window)
    ids = (point_robj.id, line_robj.id)
    isoverpoint = droprepeats(const_lift(is_same_id, m2id, ids))

    @materialize mouse_buttons_pressed, mouseposition = window.inputs

    key_pressed = const_lift(
        GLAbstraction.singlepressed,
        mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT
    )
    mousedragg = GLAbstraction.dragged(
        mouseposition, key_pressed, isoverpoint
    )
    T = Point2f0
    startvalue = (0, 0, Point2f0(0), Point2f0(0))
    sig = foldp(startvalue, mousedragg) do v0, dragg
        id, index, p0, np = v0
        if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
            id, index = value(m2id)
            if id==point_robj.id && length(gpu_position) >= index
                p0 = gpu_position[index]
            end
        else
            if id==point_robj.id && length(gpu_position) >= index
                np = p0 + T(dragg).*T(direction_restriction)
                np = T(np[1], clamp(np[2], clampto...))
                gpu_position[index] = np
            end
        end
        return id, index, p0, np
    end

    Context(line_robj, point_robj), map(x->(x[2], x[4]), sig)
end

function widget(
        line::Vector{Point{2,T}}, window;
        direction_restriction = Vec2f0(1),
        clampto = (-Inf, Inf),
        kw_args...
    ) where T
    vis, sig = edit_line(
        line, direction_restriction, clampto, window;
        kw_args...
    )
    pos_gpu = vis.children[2][:position]
    line_s = map(sig) do _
        gpu_data(pos_gpu)
    end
    vis, line_s
end

function mouse_drag_diff(past, drag_index)
    index_past, _, drag_past = past
    drag, index = drag_index
    index, drag-drag_past, drag
end

function dragg_gpu(v0, dragg)
    if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
        id, index = value(m2id)
        if id==point_robj.id && length(gpu_position) >= index
            p0 = gpu_position[index]
        else
            p0 = v0[3]
        end
    else
        id, index, p0 = v0
        if id==point_robj.id && length(gpu_position) >= index
            gpu_position[index] = Point2f0(p0) + Point2f0(dragg)
        end
    end
    return id, index, p0
end

struct ClampFunctor{T}
    a::T
    b::T
end

clampU8(x::RGBA{T}) where {T} = RGBA{T}(ntuple(i->clamp(getfield(x, i), 0.,1.), Val{4})...)
channel_color(channel, value) = RGBA{Float32}(ntuple(i->i==channel ? value : 0.0f0, Val{3})..., 1f0)
function c_setindex(color::RGBA{T}, val, channel) where T
    v = clamp.(val, 0, 1)
    RGBA{T}(
        1==channel ? v : comp1(color),
        2==channel ? v : comp2(color),
        3==channel ? v : comp3(color),
        4==channel ? v : alpha(color)
    )
end
function edit_color(tex, buff, index_value, channel, maxval)
    index, value = index_value
    if checkbounds(Bool, tex, index)
        color = c_setindex(buff[index], value[2]/maxval, channel) # we need buff, since getindex is very slow for textures
        buff[index] = color
        tex[index] = color
    end
    nothing
end

function widget(colormap::VecTypes{T}, window;
        area = (300, 30),
        slider_colors = (
            RGBA{Float32}(0.78125,0.1796875,0.41796875),
            RGBA{Float32}(0.41796875,0.78125,0.1796875),
            RGBA{Float32}(0.1796875,0.41796875,0.78125),
            RGBA{Float32}(0.9,0.9,0.9)
        ),
        knob_scale = 9f0,
        kw_args...
    ) where T<:Colorant
    colors = map(GLAbstraction.gl_promote(T), to_cpu_mem(value(colormap)))
    N = length(colors)
    color_tex = GLAbstraction.gl_convert(Texture, colormap)
    @assert colors == to_cpu_mem(color_tex)
    scale = Point2f0(area) .* Point2f0(1, 6)
    dir_restrict = Vec2f0(0,1)
    vis = ntuple(Val{4}) do i
        c_channel = Point2f0[
            Point2f0(x, getfield(c, i)) .* scale for (x, c) in zip(linspace(0,1,N), colors)
        ]
        c_i, diff = edit_line(
            c_channel, dir_restrict,
            (0, scale[2]), window, color=slider_colors[i],
            knob_scale = knob_scale
        )
        preserve(const_lift(edit_color, color_tex, colors, diff, i, scale[2]))
        c_i
    end
    tex = visualize(
        color_tex;
        is_fully_opaque=false,
        primitive=SimpleRectangle{Float32}(0, scale[2]+6, scale[1], 10),
        kw_args...
    )
    Context(tex, vis...), color_tex
end


text:
"""
$SIGNATURES

Find footnotes refs and defs and eliminate the ones that don't verify the
appropriate regex. For a footnote ref: `\\[\\^[\\p{L}0-0]+\\]` and
`\\[\\^[\\p{L}0-0]+\\]:` for the def.
"""
function validate_footnotes!(tokens::Vector{Token})
    fn_refs = Vector{Token}()
    rm      = Int[]
    for (i, τ) in enumerate(tokens)
        τ.name == :FOOTNOTE_REF || continue
        # footnote ref [^1]:
        m = match(FN_DEF_PAT, τ.ss)
        if !isnothing(m)
            if !isnothing(m.captures[1])
                # it's a def
                tokens[i] = Token(:FOOTNOTE_DEF, τ.ss)
            end
            # otherwise it's a ref, leave as is
        else
            # delete
            push!(rm, i)
        end
    end
    deleteat!(tokens, rm)
    return nothing
end

"""
$SIGNATURES

Verify that a given string corresponds to a well formed html entity.
"""
function validate_html_entity(ss::AS)
    !isnothing(match(HTML_ENT_PAT, ss))
end

"""
$(SIGNATURES)

Given a candidate header block, check that the opening `#` is at the start of a
line, otherwise ignore the block.
"""
function validate_headers!(tokens::Vector{Token})::Nothing
    isempty(tokens) && return
    s = str(tokens[1].ss) # does not allocate
    rm = Int[]
    for (i, τ) in enumerate(tokens)
        τ.name in MD_HEADER_OPEN || continue
        # check if it overlaps with the first character
        fromτ = from(τ)
        fromτ == 1 && continue
        # otherwise check if the previous character is a linereturn
        prevc = s[prevind(s, fromτ)]
        prevc == '\n' && continue
        push!(rm, i)
    end
    deleteat!(tokens, rm)
    return
end

"""
$(SIGNATURES)

Keep track of link defs. Deactivate any blocks that is within the span of
a link definition (see issue #314).
"""
function validate_and_store_link_defs!(blocks::Vector{OCBlock})::Nothing
    isempty(blocks) && return
    rm = Int[]
    parent = str(blocks[1])
    for (i, β) in enumerate(blocks)
        if β.name == :LINK_DEF
            # incremental backward look until we find a `[` or a `\n`
            # if `\n` (or start of string) first, discard
            ini  = prevind(parent, from(β))
            k    = ini
            char = '\n'
            while k ≥ 1
                char = parent[k]
                char ∈ ('[','\n') && break
                k = prevind(parent, k)
            end
            if char == '['
                # redefine the full block
                ftk = Token(:FOO,subs(""))
                # we have a [id]: lk add it to PAGE_LINK_DEFS
                id = subs(parent, nextind(parent, k), ini)
                # issue #266 in case there's formatting in the link
                id = fd2html(id, internal=true)
                id = replace(id, r"^<p>"=>"")
                id = replace(id, r"<\/p>\n$"=>"")
                lk = β |> content |> strip |> string
                PAGE_LINK_DEFS[id] = lk
                # replace the block by a full one so that it can be fully
                # discarded in the process of md blocks
                blocks[i] = OCBlock(:LINK_DEF, ftk=>ftk, subs(parent, k, to(β)))
            else
                # discard
                push!(rm, i)
            end
        end
    end
    isempty(rm) || deleteat!(blocks, rm)
    # cleanup: deactivate any block that it's in the span of a link def
    spans = [(from(ld), to(ld)) for ld in blocks if ld.name == :LINK_DEF]
    if !isempty(spans)
        for (i, block) in enumerate(blocks)
            # if it's in one of the span, discard
            curspan = (from(block), to(block))
            # early stopping if before all ld or after all of them
            curspan[2] < spans[1][1]   && continue
            curspan[1] > spans[end][2] && continue
            # go over each span break as soon as it's in
            for span in spans
                if span[1] < curspan[1] && curspan[2] < span[2]
                    push!(rm, i)
                    break
                end
            end
        end
    end
    isempty(rm) || deleteat!(blocks, rm)
    return nothing
end

"""
    find_double_brace_blocks(tokens)

Find `{{` and `}}` then form double brace blocks `{{ ... }}` (in markdown).
"""
function find_double_brace_blocks(tokens)
    # keep track of `{{` and `}}`
    db_tokens = Token[]
    # find `{` or `}` if next one exists and is the same
    # then make the first one a double with an increased substring
    # and remove the second one
    head = 1
    while head < length(tokens)
        thead = tokens[head]
        if thead.name in (:LXB_OPEN, :LXB_CLOSE)
            tcand = tokens[head+1]
            nhead = nextind(str(thead), to(thead))
            if from(tcand) == nhead && tcand.name == thead.name
                name = ifelse(thead.name == :LXB_OPEN, :DB_OPEN, :DB_CLOSE)
                ss   = subs(str(thead), from(thead), to(tcand))
                push!(db_tokens, Token(name, ss, 0))
                head += 1
            end
        end
        head += 1
    end
    # Now go back over the list of tokens and form blocks
    # these are non-nestable so it's easy.
    dbb = OCBlock[]
    head = 1
    while head < length(db_tokens)
        thead = db_tokens[head]
        if thead.name == :DB_OPEN
            hnext = findfirst(h -> db_tokens[h].name == :DB_CLOSE,
                              head+1:length(db_tokens))
            if isnothing(hnext)
                throw(OCBlockError("I found the opening token '{{' but not " *
                                   "the corresponding closing token '}}'.",
                                   context(thead)))
            end
            hnext += head
            ocb = OCBlock(:DOUBLE_BRACE, thead => db_tokens[hnext])
            push!(dbb, ocb)
            head = hnext
        end
        head += 1
    end
    return dbb
end

"""
$SIGNATURES

Check that a `---` or `***` or `___` starts at the beginning of a line and is preceded by an empty line, if that's not the case, discard it.
"""
function find_hrules!(tokens::Vector{Token})
    keep = Int[]
    rm   = Int[]
    for (i, τ) in enumerate(tokens)
        τ.name == :HORIZONTAL_RULE || continue
        # check if it has the right format
        if !(startswith(τ.ss, "---") ||
             startswith(τ.ss, "___") ||
             startswith(τ.ss, "***")) || length(unique(τ.ss)) > 1
            push!(rm, i)
            continue
        end
        # check if it's at the start of the string  or
        # if it's preceded by an empty line, in which case mark it as
        # horizontal rule, leave all other cases to be dealt with by Julia's
        # Markdown parser.
        s, k = str(τ), from(τ)
        if !(k == 1 || s[prevind(s, k, 2):prevind(s, k)] == "\n\n")
            push!(rm, i)
            continue
        end
        push!(keep, i)
    end
    hrules = tokens[keep]
    deleteat!(tokens, sort(vcat(rm, keep)))
    return hrules
end


text:
module EasyArgs

export get_arg, ARG_DICT, check_unused_args

const ARG_DICT = Dict{Any,String}()
const USED_ARGS = Set{Any}()

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@compiler_options"))
    @eval Base.Experimental.@compiler_options compile=min optimize=0 infer=false
end

function __init__()
    _parse_args!(ARG_DICT, ARGS)
end

function _parse_args!(output, input)
    position = 1
    for curr in input
        if startswith(curr, "--")
            # Long named argument
            sp = split(curr[3:end], '='; limit=2)
            name = sp[1]
            value = get(sp, 2, "")
            _put_arg!(output, name, value)
        elseif startswith(curr, "-")
            name = curr[2:end]
            _put_arg!(output, name[1:1], name[nextind(name, 1):end])
        else
            # Positional argument
            _put_arg!(output, position, curr)
            position += 1
        end
    end
    return output
end

function _put_arg!(dict, key, value)
    if haskey(dict, key)
        @warn "Duplicate argument \"$key\". Only the first occurenece will be used"
    else
        dict[key] = value
    end
end

"""
    get_arg(arg, default)

Get the command line argument `arg`. If `arg` was not given, the `default` value is returned.
If the argument exists, its value is parsed to be the same type as `default`.
"""
function get_arg(key, default::Bool)
    if haskey(ARG_DICT, key)
        if ARG_DICT[key] == ""
            push!(USED_ARGS, key)
            return true
        else
            return _get_arg(key, Bool)
        end
    else
        return default
    end
end
function get_arg(key, default::T) where {T}
    if haskey(ARG_DICT, key)
        return _get_arg(key, T)
    else
        return default
    end
end
function get_arg(key, T::Type)
    if !haskey(ARG_DICT, key)
        if key isa Int
            throw(ArgumentError("required positional argument $key not given"))
        else
            throw(ArgumentError("required argument `$key` not given"))
        end
    else
        return _get_arg(key, T)
    end
end

function _get_arg(key, T::Type)
    push!(USED_ARGS, key)
    if T <: AbstractString
        return T(ARG_DICT[key])
    else
        return parse(T, ARG_DICT[key])
    end
end

function check_unused_args(; error=false)
    unused = Pair[]
    for (k, v) in ARG_DICT
        (k ∉ USED_ARGS) && push!(unused, k => v)
    end
    if !isempty(unused)
        msg = "the following arguments were given, but not used: "
        if error
            Base.error(msg, "`", join(unused, "`, `"), "`")
        else
            @warn msg unused
        end
    end
end

end


text:
# Vector shuffles

function shufflevector_instrs(N, T, I, two_operands)
    typ = llvmtype(T)
    vtyp2 = vtyp1 = "<$N x $typ>"
    M = length(I)
    vtyp3 = "<$M x i32>"
    vtypr = "<$M x $typ>"
    mask = "<" * join(map(x->string("i32 ", x), I), ", ") * ">"
    instrs = String[]
    v2 = two_operands ? "%1" : "undef"
    push!(instrs, "%res = shufflevector $vtyp1 %0, $vtyp2 $v2, $vtyp3 $mask")
    push!(instrs, "ret $vtypr %res")
    return M, String[], instrs
end

@generated function shufflevector(v1::Vec{N,T}, v2::Vec{N,T},
                                  ::Val{I}) where {N,T,I}
    M, decls, instrs = shufflevector_instrs(N, T, I, true)
    quote
        $(Expr(:meta, :inline))
        Vec{$M,T}(Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$M,T},
            Tuple{Vec{N,T}, Vec{N,T}},
            v1, v2))
    end
end
@inline function shufflevector(v1::AbstractStructVec{N,T}, v2::AbstractStructVec{N,T}, ::Val{I}) where {N,T,I}
    SVec(shufflevector(extract_data(v1), extract_data(v2), Val(I)))
end

@generated function shufflevector(v1::Vec{N,T}, ::Val{I}) where {N,T,I}
    M, decls, instrs = shufflevector_instrs(N, T, I, false)
    quote
        $(Expr(:meta, :inline))
        Base.llvmcall($((join(decls, "\n"), join(instrs, "\n"))),
            Vec{$M,T},
            Tuple{Vec{$N,T}},
            v1)
    end
end
@inline function shufflevector(v1::AbstractStructVec{N,T}, ::Val{I}) where {N,T,I}
    SVec(shufflevector(extract_data(v1), Val(I)))
end

@inline rotate_vector_left(v::AbstractSIMDVector{N}) where {N} = shufflevector(v, Val(ntuple(i -> i, Val(N))))


text:
function smolyak_find_bracket(x::Array{T,1},nodes::Array{T,2}) where T<:AbstractFloat

  n = size(nodes,1)
  num_points = 2^size(nodes,2)

  dist_vector    = zeros(n)
  index_vector_k = Array{Int64}(num_points)
  dist_vector_k  = Array{Float64}(num_points)

  for i = 1:n
    dist_vector[i] = norm(x-nodes[i,:])
  end

  (dist,index)       = findmin(dist_vector)
  dist_vector_k[1]   = dist
  index_vector_k[1]  = index
  dist_vector[index] = Inf

  bool_vector = collect(x .< nodes[index,:])'
  counter = 1

  while counter < num_points

    (dist,index) = findmin(dist_vector)
    trial_vector = collect(x .< nodes[index,:])'

    inclusion_check = 0
    for i = 1:size(bool_vector,1)

      if sum(trial_vector .== bool_vector[i,:]) == size(bool_vector,2)
        inclusion_check += 1
      end

    end

    if inclusion_check == 0
      counter += 1
      dist_vector_k[counter]  = dist
      index_vector_k[counter] = index
      dist_vector[index]    = Inf
      bool_vector = [bool_vector; trial_vector]
    else
      dist_vector[index]    = Inf
    end

  end

  return index_vector_k,dist_vector_k

end

function smolyak_bracket_evaluate(y::AbstractArray{T,1},nodes::Array{T,2},x::Array{T,1}) where T<:AbstractFloat

  (index_vector_k,dist_vector_k) = smolyak_find_bracket(x,nodes)
  k = size(index_vector_k,1)

  weights = Array{T}(k)
  for i = 1:k
    weights[i] = 1.0/(dist_vector_k[i]+eps(T))
  end

  weights = weights./sum(weights)

  y_hat = zero(T)

  for i = 1:k
    y_hat += y[index_vector_k[i]]*weights[i]
  end

  return y_hat

end


text:
using ConstraintSolver
CS = ConstraintSolver

function main(filename; benchmark=false)
    com = CS.init()

    lines = readlines(filename)
    num_colors = 0
    x = nothing
    for line in lines
        parts = split(line, " ")
        if parts[1] == "p"
            num_colors = parse(Int, parts[3])
            x = Vector{CS.Variable}(undef, num_colors)
            for i=1:num_colors
                x[i] = add_var!(com, 1, num_colors)
            end
        elseif parts[1] == "e"
            f = parse(Int, parts[2])
            t = parse(Int, parts[3])
            add_constraint!(com, x[f] != x[t])
        end
    end
    println("num_colors: ", num_colors)

    set_objective!(com, :Min, CS.vars_max(x))

    status = solve!(com)

    if !benchmark
        println("status: ", status)
        println("objective: ", com.best_sol)
        com.info
    end
end


text:
# # Useful patterns
#
# This page includes some useful patterns using Transducers.jl.

using Transducers

using Test                                                             #src

# ## Flattening nested objects using `MapCat`
#
# ### Simple `MapCat` usage
#
# Consider a vector of "objects" (here just `NamedTuple`s) which in
# turn contain a vector of objects:

#! format: off
nested_objects = [
    (a = 1,  b = [(c = 2,  d = 3),  (c = 4,  d = 5)]),
    (a = 10, b = [(c = 20, d = 30), (c = 40, d = 50)]),
];
#! format: on

# We can flatten this into a table by using `Map` inside `MapCat`:

using TypedTables
astable(xs) = copy(Table, xs)  # using `TypedTables` for a nice display

table1 = nested_objects |> MapCat() do x
    x.b |> Map() do b  # not `MapCat`
        (a = x.a, b...)
    end
end |> astable
@test table1 == Table(a = [1, 1, 10, 10], c = [2, 4, 20, 40], d = [3, 5, 30, 50])

# (Note that the transducer used inside `MapCat` is `Map`, not `MapCat`)

# ### Nested `MapCat`
#
# This pattern can handle more nested objects:

#! format: off
more_nested_objects = [
    (a = 1,  b = [(c = 2,  d = [(e = 3,  f = 4),  (e = 4,  f = 5)]),
                  (c = 6,  d = [])]),
    (a = 10, b = [(c = 20, d = [(e = 30, f = 40), (e = 40, f = 50)])]),
];
#! format: on

# By using nested `MapCat` (except for the "inner most" processing
# which uses `Map` since there is nothing to con**cat**enate):

table3 =
    more_nested_objects |> MapCat() do x
        x.b |> MapCat() do b
            b.d |> Map() do d
                (a = x.a, c = b.c, d...)
            end
        end
    end |> astable

# ### Comparison with iterator comprehension
#
# As a comparison, here is how to do it with iterator comprehension

#! format: off
rows = (
    (a = x.a, c = b.c, d...)
    for x in more_nested_objects
    for b in x.b
    for d in b.d
)
@assert Table(collect(rows)) == table3
#! format: on

# For a simple flattening and mapping, iterator comprehension as above
# perhaps is the simplest solution.
#
# Note that Transducers.jl works well with iterator comprehensions.
# Transducers.jl-specific entry points like [`foldxl`](@ref) converts
# iterator comprehensions to transducers internally.
# [`eduction`](@ref) can be used to explicitly do this conversion:

@assert astable(eduction(rows)) == table3

# ### Complex `MapCat` example
#
# For more complex processing that requires intermediate variables,
# the iterator comprehension does not work well.  Fortunately, it is
# easy to use intermediate variables with transducers:

ans = begin
    more_nested_objects |>
    MapCat() do x
        a2 = x.a * 2
        x.b |> MapCat() do b
            a2_plus_c = a2 + b.c
            b.d |> Map() do d
                c_plus_e = b.c + d.e
                c_plus_f = b.c + d.f
                (a2_plus_c = a2_plus_c, c_plus_e = c_plus_e, c_plus_f = c_plus_f)
            end
        end
    end |>
    astable
end

@testset begin
    @test ans == Table(
        a2_plus_c = [4, 4, 40, 40],
        c_plus_e = [5, 6, 50, 60],
        c_plus_f = [6, 7, 60, 70],
    )
end

# ### `MapCat` with `zip`
#
# Note also that `MapCat` can be combined with arbitrary iterator
# combinators such as `zip`

ans = begin
    [(a = 1:3, b = 'x':'z'), (a = 1:4, b = 'i':'l')] |>
    MapCat() do x
        zip(x.a, x.b)
    end |>
    MapSplat((a, b) -> (a = a, b = b)) |>
    astable
end

@testset begin
    @test ans == Table(a = [1, 2, 3, 1, 2, 3, 4], b = ['x', 'y', 'z', 'i', 'j', 'k', 'l'])
end

# ### `MapCat` with `Iterators.product`
#
# ... and `product`

ans = begin
    [(a = 1:3, b = 'x':'z'), (a = 1:4, b = 'i':'l')] |>
    MapCat() do x
        Iterators.product(x.a, x.b)
    end |>
    Enumerate() |>
    Filter(x -> x[1] % 5 == 0) |>  # include only every five item
    MapSplat((n, (a, b)) -> (n = n, a = a, b = b)) |>
    astable
end

@testset begin
    @test ans ==
          Table(n = [5, 10, 15, 20, 25], a = [2, 1, 2, 3, 4], b = ['y', 'i', 'j', 'k', 'l'])
end

# ## "Missing value" handling with `KeepSomething`

# Transducers.jl has a generic filtering such as [`Filter`](@ref) as
# well as type-based filtering such as [`NotA`](@ref) and
# [`OfType`](@ref).  These transducers can be used to filter out
# "missing values" represented as `missing` or `nothing`.
#
# [`KeepSomething`](@ref) is a transducer that is useful for working
# on `Union{Nothing,Some{T}}`.  It filters out `nothing` and yield
# itmes after applying `something`.

@test begin
    [nothing, 1, Some(nothing), 2, 3] |> KeepSomething(identity) |> collect
end == [1, nothing, 2, 3]

# Thus, `KeepSomething` works well with any tools that operate on
# `Union{Nothing,Some{T}}`.  Here is an example of using it with
# [Maybe.jl](https://github.com/tkf/Maybe.jl).  Consider a vector of
# heterogeneous dictionaries with varying set of keys:

heterogeneous_objects = [
    Dict(:a => 1, :b => Dict(:c => 2)),
    Dict(:a => 1),                          # missing key
    Dict(:a => 1, :b => Dict()),            # missing key
    Dict(:b => Dict(:c => 2)),              # missing key
    Dict(:a => 10, :b => Dict(:ccc => 20)), # alternative key name
];

# Using `@something` and `@?` macros from Maybe.jl, we can convert
# this to a regular table quite easily:

using Maybe
using Maybe: @something

ans = begin
    heterogeneous_objects |>
    KeepSomething() do x
        c = @something {       # (1)
            @? x[:b][:c];      # (2)
            @? x[:b][:ccc];    # (3)
            return;            # (4)
        }
        @? (a = x[:a], c = c)  # (5)
    end |>
    astable
end

@testset begin
    @test ans == Table(a = [1, 10], c = [2, 20])
end

# In this example, for each dictionary `x`, the body of the `do` block
# works as follows:
#
# * (1) Try to extract the item `c`.
#   * (2) First, try to get it from `x[:b][:c]`.
#   * (3) If `x[:b][:c]` doesn't exist, try `x[:b][:ccc]` next.
#   * (4) If both `x[:b][:c]` and `x[:b][:ccc]` do not exist, return
#     `nothing`.  `KeepSomething` will filter out this entry.
# * (5) Try to extract the item `a` from `x[:a]`.
#   * If this does not exist, the whole expression wrapped by `@?`
#     evaluates to `nothing`.  This, in turn, will be filtered out by
#     `KeepSomething`.
#   * If `x[:a]` exists, `@? (a = x[:a], c = c)` evaluates to `Some((a
#     = value_of_a, c = value_of_c))`.  The `Some` wrapper is
#     unwrapped by `something` called by `KeepSomething`.
#
# For more information, see the tutorial in [Maybe.jl
# documentation](https://tkf.github.io/Maybe.jl/dev/).

# ## Multiple outputs

# Usually, reducers like `sum` and `collect` have one output.  However
# we can use [`TeeRF`](@ref) etc. to "fan-out" input items to multiple
# outputs.
#
# ### Multiple output vectors
#
# Here is an example of creating two output vectors of integers and
# symbols in one go:

ints, symbols =
    [1, :two, missing, 3, 4, :five, 6] |>
    Filter(!isequal(6)) |>
    foldxl(TeeRF(
        OfType(Int)'(push!!),    # push integers to a vector
        OfType(Symbol)'(push!!), # push symbols to a vector
    ))

@testset begin
    @test ints == [1, 3, 4]
    @test symbols == [:two, :five]
end

# Here, we use `TeeRF(rf₁, rf₂, ..., rfₙ)` to fan-out input items to
# multiple reducing functions.  To compose each reducing function, we
# use [`OfType`](@ref) transducer as reducing function transformation
# [`xf'(rf)`](@ref adjoint).
#
# ### Handling empty results
#
# Note that fold with `push!!` throws when the input is empty.  To
# obtain an empty vector when the input is empty or all filtered out,
# we need to specify `init`.
# [MicroCollections.jl](https://github.com/JuliaFolds/MicroCollections.jl)
# includes a library of collections useful as `init`.  Here, we can
# use `EmptyVector`:

using MicroCollections

ints, strings =
    [1, :two, missing, 3, 4, :five, 6] |>
    Filter(!isequal(6)) |>
    foldxl(TeeRF(
        OfType(Int)'(push!!),    # push integers to a vector
        OfType(String)'(push!!), # push strings to a vector (but there is no string)
    ); init = EmptyVector())

@testset begin
    @test ints == [1, 3, 4]
    @test strings == []
end

# ### Composed transducers with `TeeRF`
#
# Each reducing function passed to `TeeRF` can use arbitrary complex
# transducers.  Here is an example of filtering-in symbols and then
# map them to strings:

ints, strings =
    [1, :two, missing, 3, 4, :five, 6] |>
    Filter(!isequal(6)) |>
    foldxl(
        TeeRF(
            OfType(Int)'(push!!),
            opcompose(OfType(Symbol), Map(String))'(push!!),  # filter _then_ map
        );
        init = EmptyVector(),
    )

@testset begin
    @test ints == [1, 3, 4]
    @test strings == ["two", "five"]
end

# ### Nested `TeeRF`
#
# Each reducing function itself passed to `TeeRF` can even be composed
# using `TeeRF` (or other reducing function combinators; e.g.,
# [`ProductRF`](@ref)).  Here is an example of computing extrema on
# integers:

(imin, imax), strings =
    [1, :two, missing, 3, 4, :five, 6] |>
    Filter(!isequal(6)) |>
    foldxl(
        TeeRF(
            OfType(Int)'(TeeRF(max, min)),  # extrema on integers
            opcompose(OfType(Symbol), Map(String))'(push!!),  # filter _then_ map
        );
        init = ((typemin(Int), typemax(Int)), EmptyVector()),
    )

@testset begin
    @test imin == 4
    @test imax == 1
    @test strings == ["two", "five"]
end

# ### When input is a tuple: `ProductRF`
#
# [`ProductRF`](@ref) is like `TeeRF` but it expects that the input is
# already a tuple:

ints, io =
    [(1:3, 'x':'z'), nothing, (1:4, 'i':'l')] |>
    NotA(Nothing) |>
    foldxl(
        ProductRF(
            opcompose(Cat(), Filter(isodd))'(push!!),    # process 1:3 etc.
            Cat()'((io, char) -> (write(io, char); io)), # process 'x':'z' etc.
        );
        init = (EmptyVector(), IOBuffer()),
    );
#-

@test begin
    String(take!(io))
end == "xyzijkl"
#-

@test begin
    ints
end == [1, 3, 1, 3]
#-

# ### When input is a row: `DataTools.oncol`
#
# [`oncol`](https://juliafolds.github.io/DataTools.jl/dev/#DataTools.oncol)
# from [DataTools.jl](https://github.com/JuliaFolds/DataTools.jl) is
# like `ProductRF` but acts on `NamedTuple` (as well as any
# [Setfield.jl](https://github.com/jw3126/Setfield.jl)-compatible
# possibly nested objects).

using DataTools
@test begin
    foldxl(oncol(a = +, b = *), [(a = 1, b = 2), (a = 3, b = 4)])
end == (a = 4, b = 8)


text:
using WAV
using SpeechAugment: speedWav

# speed parameter
maxspeed = 1.2
minspeed = 0.8
# read a wav file
wav = wavread("xxx.wav")[1]
# speedup speech by speed parameter
y = speedWav(wav, rand()*(maxspeed-minspeed) + minspeed)

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
import Distributions
import DataFrames
import GLM
import Random
import StatsBase

n = 500
p = 10
X = hcat(ones(n), randn(n, p-1))

noise_coefs = zeros(p - 5)
β = vcat([1.5, 1.1, -10.5, 4.5, 1.9], noise_coefs)
η = X*β                                                   # linear predictor
pr = 1 ./(1 + exp.(-η))                                     # inv-logit

# simulate outcome variable
y = map(π -> rand(Distributions.Binomial(1, π)), pr)
mean(y)

df = DataFrames.DataFrame(X)
df[:y] = y

# keep only a sub-sample of positive cases
keep_prop = 0.05
pos_cases = findall(y .== 1)
keep_indcs = vcat(StatsBase.sample(pos_cases, Int(round(2 * keep_prop * length(pos_cases)))), findall(y .== 0))
df2 = df[keep_indcs, :]

sum(df2[:y])/n

GLM.glm(y ~ 1 + x2 + x3 + x4 + x5 + x6 + x7, df2, Binomial())

# Simulation of class-imbalanced data requires data drawn from
# multivariate distribution (here we use MV Normal).

n = 100
p = 10
# X = hcat(ones(n), randn(n, p-1))

Σ = [1.0 0.8 0.4 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
     0.8 1.0 0.6 0.2 0.0 0.0 0.0 0.0 0.0 0.0;
     0.4 0.6 1.0 0.6 0.4 0.0 0.0 0.0 0.0 0.0;
     0.0 0.2 0.6 1.0 0.8 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.4 0.8 1.0 0.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
     0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]

mvn = MvNormal(ones(p), Σ)

Random.seed!(round(Int, time()))

X = rand(mvn, n)'
cor(X)

ϵ = rand(Normal(2, 0.5), n)                                  # gives 10% minority class,
# ϵ = rand(Normal(-2.5, 0.5), n)                               # gives 5% minority class,
# ϵ = rand(Normal(4, 0.5), n)                                  # gives 15% minority class,

β = [-3, -4, 1, -5, -4, 0.0, 0.0, 0.0, 0.0, 0.0]
η = X*β + ϵ                                                    # linear predictor w/ error
pr = 1.0 ./ (1.0 + exp.(-η))                                    # inv-logit

# simulate outcome variable
y = map(π -> rand(Binomial(1, π)), pr)


text:
using SparseTransforms
using Test
using Pipe.Pipe
using LinearAlgebra

function test_reconstruct()
    query_method = :simple

    n = 4
    locs = [4, 6, 10, 15]
    strengths = Float64.([2, 4, 1, 1])
    signal = TestSignal(n, locs, strengths, 0.0)
    b = get_b(signal; method=query_method)
    Ms = get_Ms(signal.n, b; method=query_method)
    K = binary_ints(signal.n)

    select_froms = Array{Int64}[]
    for M in Ms
        selects = @pipe M' * K |> transpose |> Bool.(_) |> eachrow |> bin_to_dec.(_)
        push!(select_froms, selects)
    end

    cutoff = 2 * signal.noise_sd ^ 2 * (2 ^ (signal.n - b)) * (signal.n + 1) # noise threshold

    for method_name in all_methods["reconstruct"]
        println("Testing reconstruction with $method_name")
        if method_name == :nso
            delays_method = :nso
        else
            delays_method = :identity_like
        end
        D = get_D(signal.n; method=delays_method, num_delays=2*n)
        Us = Array{Float64}[]
        for M in Ms
            U, used_i = compute_delayed_subtransform(signal, M, D, fwht)
            push!(Us, U)
        end
        S = (-1) .^ (D * K)
        singleton_guesses = Int64[]
        for (i, (U, select_from)) in enumerate(zip(Us, select_froms))
            col_gen = U |> eachrow |> enumerate
            for (j, col) in col_gen
                if col⋅col > cutoff
                    selection = findall(==(j-1), select_from) # pick all the k such that M.T @ k = j
                    k, sgn = singleton_detection(
                        col;
                        method=method_name,
                        selection=selection,
                        S_slice=S[:, selection],
                        n=signal.n,
                        D=D
                    ) # find the best fit singleton
                    push!(singleton_guesses, k |> bin_to_dec)
                end
            end
        end
        @test singleton_guesses == [4, 6, 15, 6, 10, 15]
    end
end


text:
# Provides a structure for turning a mdp into a rpomdp where observations are the states of the original mdp

struct FullyObservablePOMDP{S, A} <: RPOMDP{S,A,S}
    mdp::MDP{S, A}
end

# observations are the state of the MDP
# The observation distribution is modeled by a SparseCat distribution with only one element
RPOMDPs.observations(rpomdp::FullyObservablePOMDP) = states(rpomdp.mdp)
RPOMDPs.n_observations(rpomdp::FullyObservablePOMDP) = n_states(rpomdp.mdp)
RPOMDPs.obs_index{S, A}(rpomdp::FullyObservablePOMDP{S, A}, o::S) = state_index(rpomdp.mdp, o)

RPOMDPs.convert_o(T::Type{V}, o, rpomdp::FullyObservablePOMDP) where {V<:AbstractArray} = convert_s(T, s, rpomdp.mdp)
RPOMDPs.convert_o(T::Type{S}, vec::V, rpomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, rpomdp.mdp)

function RPOMDPs.generate_o(rpomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG)
    return s
end

function RPOMDPs.observation(rpomdp::FullyObservablePOMDP, s, a)
    return SparseCat((s,), (1.,))
end

function RPOMDPs.observation(rpomdp::FullyObservablePOMDP, s, a, sp)
    return SparseCat((sp,), (1.,))
end

RPOMDPs.isterminal_obs{S,A}(problem::FullyObservablePOMDP{S,A}, o::S) = isterminal(rpomdp.mdp, o)

# inherit other function from the MDP type

RPOMDPs.states(rpomdp::FullyObservablePOMDP) = states(rpomdp.mdp)
RPOMDPs.actions(rpomdp::FullyObservablePOMDP) = actions(rpomdp.mdp)
RPOMDPs.transition{S,A}(rpomdp::FullyObservablePOMDP{S,A}, s::S, a::A) = transition(rpomdp.mdp, s, a)
RPOMDPs.initial_state_distribution(rpomdp::FullyObservablePOMDP) = initial_state_distribution(rpomdp.mdp)
RPOMDPs.initial_state(rpomdp::FullyObservablePOMDP, rng::AbstractRNG) = initial_state(rpomdp.mdp, rng)
RPOMDPs.generate_s(rpomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_s(rpomdp.mdp, s, a, rng)
RPOMDPs.generate_sr(rpomdp::FullyObservablePOMDP, s, a, rng::AbstractRNG) = generate_sr(rpomdp.mdp, s, a, rng)
RPOMDPs.reward{S,A}(rpomdp::FullyObservablePOMDP{S, A}, s::S, a::A) = reward(rpomdp.mdp, s, a)
RPOMDPs.isterminal(rpomdp::FullyObservablePOMDP, s) = isterminal(rpomdp.mdp, s)
RPOMDPs.discount(rpomdp::FullyObservablePOMDP) = discount(rpomdp.mdp)
RPOMDPs.n_states(rpomdp::FullyObservablePOMDP) = n_states(rpomdp.mdp)
RPOMDPs.n_actions(rpomdp::FullyObservablePOMDP) = n_actions(rpomdp.mdp)
RPOMDPs.state_index{S,A}(rpomdp::FullyObservablePOMDP{S,A}, s::S) = state_index(rpomdp.mdp, s)
RPOMDPs.action_index{S, A}(rpomdp::FullyObservablePOMDP{S, A}, a::A) = action_index(rpomdp.mdp, a)
RPOMDPs.convert_s(T::Type{V}, s, rpomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_s(T, s, rpomdp.mdp)
RPOMDPs.convert_s(T::Type{S}, vec::V, rpomdp::FullyObservablePOMDP) where {S,V<:AbstractArray} = convert_s(T, vec, rpomdp.mdp)
RPOMDPs.convert_a(T::Type{V}, a, rpomdp::FullyObservablePOMDP) where V<:AbstractArray = convert_a(T, a, rpomdp.mdp)
RPOMDPs.convert_a(T::Type{A}, vec::V, rpomdp::FullyObservablePOMDP) where {A,V<:AbstractArray} = convert_a(T, vec, rpomdp.mdp)


text:
using Test, AutoARIMA

@testset "Mathematical Properties" begin

    @test reduce(&, yule_walker(seriesF, 15) .≈ levinson_durbin(seriesF, 15))

    @testset "MA(q) forecast after q steps should be equal to the constant" begin
        μ = 0.1
        ma = MAModel{3}(μ,[0.2,0.3,-0.5],1.0)
        z = [0.2,0.5,0.1]
        f = Float64[]
        for i in 1:10
            zhat = forecast(ma,z)
            push!(z,zhat)
            push!(f,zhat)
        end
        @test all(f[4:end] .≈ μ)
    end

    @testset "BoxCox" begin
        z = rand(1000)
        @test z ≈ inv_boxcox(boxcox(z,0.1),0.1)
        @test z ≈ inv_boxcox(boxcox(z,0.0),0.0)
    end

end


text:
module ViewHelper

using Genie, Genie.Helpers, SearchLight, Genie.Router

export output_flash, book_cover, book_form_uri

function output_flash(params::Dict{Symbol,Any}) :: String
  ! isempty( flash(params) ) ? """
$(flash(params))
""" : ""
end

function book_cover(book)
  isempty(book.cover) ? "img/docs.png" : book.cover
end

function book_form_uri(book)
  ispersisted(book) ? link_to(:update_book, id = book.id |> get) : link_to(:create_book)
end

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
# This file was generated, do not modify it. # hide
using Plots
@views scatter(db.matrix[1, :], db.matrix[2, :], fmt=:png, size=(600, 600), color=:cyan, ma=0.3, a=0.3, ms=1, msw=0, label="")
for c in eachcol(I)
    R = db.matrix[:, c]
    @views scatter!(R[1, :], R[2, :], m=:diamond, ma=0.3, a=0.3, color=:auto, ms=2, msw=0, label="")
end

@views scatter!(Q.matrix[1, :], Q.matrix[2, :], color=:black, m=:star, ma=0.5, a=0.5, ms=4, msw=0, label="")

savefig(joinpath(@OUTPUT, "fig-2d-t1.png")) # hide

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
#
# errors.jl --
#
# Management of errors for the Julia interface to the C libraries of TAO, a
# Toolkit for Adaptive Optics.
#
#------------------------------------------------------------------------------
#
# This file is part of TAO software (https://git-cral.univ-lyon1.fr/tao)
# licensed under the MIT license.
#
# Copyright (C) 2018-2021, Éric Thiébaut.
#

# Structure to store a single error of the TAO Library.
struct ErrorInfo
    func::String
    code::Cint
    proc::Ptr{Cvoid}
end

"""
    TaoBindings.TaoError(errs)

yields an exception to be thrown in case of errors in TAO Library.
Argument `errs` is a vector of `TaoBindings.ErrorInfo` elements.

""" TaoError

struct TaoError <: Exception
    errs::Array{ErrorInfo}
end

"""
    TaoBindings.get_error_stack()

yields the address of the calling thread error stack.

""" get_error_stack

struct ErrorStack
    top::Ptr{Cvoid} # pointer to last error
end

# We use a per-thread cache to speed-up querying the address of the calling
# thread error stack. On my Intel Core i7-5500U: _call_get_error_stack() takes
# 26ns while get_error_stack() only takes 5ns.
const ERROR_STACK = Ptr{ErrorStack}[]
function get_error_stack()
    id = Threads.threadid()
    while length(ERROR_STACK) < id
        push!(ERROR_STACK, C_NULL)
    end
    @inbounds begin
        if ERROR_STACK[id] == C_NULL
            ERROR_STACK[id] = _call_get_error_stack()
        end
        return ERROR_STACK[id]
    end
end

_call_get_error_stack() =
    ccall((:tao_get_error_stack, taolib), Ptr{ErrorStack}, ())

"""
    TaoBindings.any_errors()

yields whether there are any errors in the calling thread error stack.

"""
any_errors() = any_errors(get_error_stack())
any_errors(ptr::Ptr{ErrorStack}) = (unsafe_load(ptr)).top != C_NULL

"""
    TaoBindings._check(ans::Status, warn=false)

or

    TaoBindings._check(success::Bool, warn=false)

deals with error reporting after calling one or several functions of the TAO
Library.  Argument `ans` is the status value returned by a call to a TAO
function returning a `tao_status_t`; the call is considered as successful if
`ans != TaoBindings.ERROR`.  Another possibility is to specify a boolean value
`success` to indicate whether the call(s) to the function(s) in the TAO Library
was (were) successful.  If there are errors, optional argument `warning`
(`false` by default) indicates whether to just print the error messages as
warnings or to throw an exception.

Note that, if the first argument indicates a success but there are some
unreported erros in the calling thread error stacks, these errors are printed
as warnings.

!!! warn
    All methods in `TaoBindings` which `ccall` one of the error prone fucntions
    of the TAO Library shall call `TaoBindings._check` before they return to
    avoid memory leaks.

"""
_check(ans::Status, warn::Bool=false) = _check(ans != ERROR, warn)

@inline function _check(success::Bool, warn::Bool=false)
    stack = get_error_stack()
    if !success || any_errors(stack)
        _something_went_wrong(success, stack, warn)
    end
    nothing
end

# Slow method to deal with errors/warnings.  Should be called by _check only if
# something unexpected occured.
@noinline function _something_went_wrong(success::Bool,
                                         stack::Ptr{ErrorStack},
                                         warn::Bool)
    if any_errors(stack)
        if success || warn
            _call_report_errors(stack, "WARNING ", "        ", "\n")
        else
            throw(TaoError(pop_errors(stack)))
        end
    elseif !success
        mesg = "failure detected in a call to a function of TAO Library without error information!"
        if warn
            println("WARNING ", mesg)
        else
            error(mesg)
        end
    end
end

"""
    TaoBindings.get_error_reason(err)

yields the error message associated to `err` which is either an instance of
`TaoBindings.ErrorInfo` or an integer error code.

"""
get_error_reason(err::ErrorInfo) = get_error_reason(err.code)
get_error_reason(code::Integer) =
    unsafe_string(ccall((:tao_get_error_reason, taolib), Ptr{UInt8},
                        (Cint,), code))

"""
    TaoBindings.get_error_name(err)

yields the literal name associated to `err` which is either an instance of
`TaoBindings.ErrorInfo` or an integer error code.

"""
get_error_name(err::ErrorInfo) = get_error_name(err.code)
get_error_name(code::Integer) =
    unsafe_string(ccall((:tao_get_error_name, taolib), Ptr{UInt8},
                        (Cint,), code))

"""
    TaoBindings.report_errors(first, other, last, sfx)

prints the errors of the calling thread to the standard error output and clear
the contents of the error stack.  Arguments `first`, `other` and `last` are the
prefixes to use for printing the first, subsequent and last errors.  Argument
`sfx` is the suffix to use for printing errors.

"""
report_errors() = report_errors("(TAO-ERROR) ",
                                "         ├─ ",
                                "         └─ ", "\n")
function report_errors(first::AbstractString, other::AbstractString,
                       last::AbstractString, sfx::AbstractString)
    _call_report_errors(get_error_stack(), first, other, last, sfx)
end

function _call_report_errors(stack::Ptr{ErrorStack}, first::AbstractString,
                             other::AbstractString, last::AbstractString,
                             sfx::AbstractString)
    ccall((:tao_report_errors_to_stream, taolib), Cvoid,
          (Ptr{Cvoid}, Ptr{Cvoid}, Cstring, Cstring, Cstring, Cstring),
          C_NULL, Ptr{ErrorStack}, first, other, last, sfx)
end

"""
    TaoBindings.discard_errors(stack=get_error_stack())

discards the contents of the calling thread error stack.

"""
discard_errors(stack::Ptr{ErrorStack} = get_error_stack()) =
    ccall((:tao_discard_errors, taolib), Cvoid, (Ptr{ErrorStack},), stack)

"""
    TaoBindings.pop_errors(stack=get_error_stack())

pops the errors out of the calling thread error stack and returns an array of
instances of `TaoBindings.ErrorInfo` (possibly empty if there are no errors).

"""
function pop_errors(stack::Ptr{ErrorStack} = get_error_stack())
    func = Ref{Ptr{UInt8}}(0)
    code = Ref{Cint}(0)
    proc = Ref{Ptr{Cvoid}}(0)
    errs = ErrorInfo[]
    while 0 != ccall((:tao_pop_error, taolib), Cint,
                     (Ptr{ErrorStack}, Ptr{Ptr{UInt8}}, Ptr{Cint},
                      Ptr{Ptr{Cvoid}}), stack, func, code, proc)
        funcname = func[] == C_NULL ? "" : unsafe_string(func[])
        push!(errs, ErrorInfo(funcname, code[], proc[]))
    end
    return errs
end

struct ErrorInfoWorkspace
    reason::Ref{Ptr{UInt8}}
    info::Ref{Ptr{UInt8}}
    buffer::Vector{UInt8}
    ErrorInfoWorkspace() = new(Ref{Ptr{UInt8}}(0),
                               Ref{Ptr{UInt8}}(0),
                               Array{UInt8}(undef, 20))
end

unsafe_retrieve_error_details(err::ErrorInfo) =
    unsafe_retrieve_error_details!(err, ErrorInfoWorkspace())
unsafe_retrieve_error_details!(err::ErrorInfo, wrk::ErrorInfoWorkspace) = begin
    GC.@preserve wrk begin # FIXME: preserve not necessary?
        ccall((:tao_retrieve_error_details, taolib), Cvoid,
              (Cint, Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}},
               Ptr{Ptr{Cvoid}}, Ptr{UInt8}),
              err.code, wrk.reason, wrk.info, err.proc, wrk.buffer)
        return (err.func, unsafe_string(wrk.reason[]),
                unsafe_string(wrk.info[]))
    end
end

function showerror(io::IO, err::ErrorInfo)
    func, reason, info = unsafe_retrieve_error_details(err)
    print(io, "TaoBindings.ErrorInfo: ", reason, " in function `",
          func, "` [", info, ']')
end

function showerror(io::IO, e::TaoError)
    wrk = ErrorInfoWorkspace()
    errs = e.errs
    prefix1= "TaoError: "
    prefix2= "          "
    for i in 1:length(errs)
        func, reason, info = unsafe_retrieve_error_details!(errs[i], wrk)
        print(io, (i == 1 ? prefix1 : prefix2), reason, " in function `",
              func, "` [", info, ']')
    end
end


text:
using Compat, DSPRemezWrap, Compat.Test

testfiles = [ "remez_wrap.jl" ]

for testfile in testfiles
    eval(@__MODULE__, :(@testset $testfile begin include($testfile) end))
end


text:
import Dierckx
import DSP: conv

const actionc_steer = 7
const actionc_accel = 5
const MINSPEED = -0.9
const MAXACCEL = 4.0 # meter / second^2
const MAXSTEER = pi/30  # in radians. This is 90 deg of steerwheel asumming steer factor of 15
const MAXSPEED = 14.0 # in meter / second
const YOLOCAR_MAX_VEL = 5.0 # meter / second
const YOLOCAR_MIN_VEL = 0.9 # meter / second
const YOLOCAR_MAX_ACCEL = 6.0 # meter / second^2
const DRAWING_SCALE = 0.3
const SAFETY = 1.3

include("draw.jl")
include("car.jl")
include("collisions.jl")

function gen_simple_road(n, width, len, curvy, smth, seed=1)
    #rand(seed)
    smth = smth*(1000/n)
    da = (rand(n + 102) .- 0.5) .* 2.0 .* pi .* curvy
    a = Float64[da[1]]
    for i in 2:length(da)
        v = a[end] + da[i]
        if abs(v)>pi/2
            v -= 2.0*da[i]
        end
        push!(a, v)
    end
    l = range(-100, stop=100, length=n)
    g = exp.(-l .* l / smth / smth)
    a = conv(g, a)[div(n,2):end-div(n,2)+1]
    a /= maximum(abs.(a))
    a *= pi/2

    a .-= a[50]
    dl = len / n
    dx = dl * cos.(a)
    dy = dl * sin.(a)
    xl = cumsum(dx)
    yl = cumsum(dy)
    a .+= pi / 2
    dx = width .* cos.(a)
    dy = width .* sin.(a)
    xr = xl + dx
    yr = yl + dy
    xr = xr[50:end-50]
    yr = yr[50:end-50]
    xl = xl[50:end-50]
    yl = yl[50:end-50]
    a = a[50:end-50]
    x0 = xr[1]
    y0 = 0.5*(yr[1].+yl[1])
    return xr.-x0, yr.-y0, xl.-x0, yl.-y0, 0.5*(xr.+xl).-x0, 0.5*(yr.+yl).-y0, a
end

struct MyState
    rcx::Float32
    rcy::Float32
    angle::Float32
    steer::Float32
    accel::Float32
    speed::Float32
    yolocar_dist::Float32
    yolocar_timestep::Int32
    times_passed_yolocar::Int32
end

MyState(car::Car, yolocar_dist, yolocar_timestep, times_passed) =
    MyState(car.rcx,car.rcy,car.angle,car.steer,car.accel,car.speed,yolocar_dist, yolocar_timestep, times_passed)

mutable struct MyEnv
    coll::Collider
    max_dist::Float64
    splx::Dierckx.Spline1D
    sply::Dierckx.Spline1D
    yolocar_width::Float64
    yolocar_height::Float64
    yolocar_vels::Vector{Float32}
    triangle_points::Vector{Tuple{Float32,Float32}}
end

function action_ix_to_action(::MyEnv, ix::Int32)
    R = div(ix-1, actionc_accel) + 1
    C = ix-(R-1)*actionc_accel
    steer = range(-MAXSTEER, stop=MAXSTEER, length=actionc_steer)[R]
    accel = range(-MAXACCEL, stop=MAXACCEL, length=actionc_accel)[C]
    (steer, accel)
end

function initialize_simple_road(xl,yl,xr,yr,xc,yc,v=-1.0)
    segs = Segment[]
    triangle_points = Tuple{Float64,Float64}[]
    ix = 1
    for i in 1:(length(xr)-1)
        push!(segs,Segment(xl[i],yl[i],xl[i+1],yl[i+1],ix))
        push!(triangle_points, (xr[i+1],yr[i+1]))
        ix +=1
        push!(segs,Segment(xr[i],yr[i],xr[i+1],yr[i+1],ix))
        push!(triangle_points, (xl[i],yl[i]))
        ix += 1
    end
    push!(segs,Segment(xr[1],yr[1],xl[1],yl[1],ix))
    speed = if v<0.0 0.01 + rand()*MAXSPEED else v end
    av = if v<0.0
        Car(2.0 + 2.0*(rand()-0.5),2.0*(rand()-0.5), randn()*0.08,0.0,0.0,speed)
    else
        Car(2.0,0.0, 0.0,0.0,0.0,speed)
    end
    col = Collider(71.0)
    addp!(col, segs)

    yolocar_dfac = 0.1 + rand()*0.8 # right vs left weighting
    path_x = xl.*yolocar_dfac + xr.*(1.0-yolocar_dfac)
    path_y = yl.*yolocar_dfac + yr.*(1.0-yolocar_dfac)
    dx = diff(path_x)
    dy = diff(path_y)
    d = [0; cumsum(sqrt.(dx.*dx + dy.*dy))]
    splx = Dierckx.Spline1D(d, path_x)
    sply = Dierckx.Spline1D(d, path_y)

    yolocar_dist = 40.0 + randn()*7.0
    yolocar_width = 1.6 + randn()*0.1
    yolocar_height = 1.5 + randn()*0.1

    yolocar_vels = Float64[rand()*YOLOCAR_MAX_VEL]
    for _ in 1:50000
        a = 2*(rand()-0.5)*YOLOCAR_MAX_ACCEL
        dv = 0.1*a
        if yolocar_vels[end]+dv>YOLOCAR_MAX_VEL || yolocar_vels[end]+dv<YOLOCAR_MIN_VEL
            dv = -dv
        end
        push!(yolocar_vels, yolocar_vels[end]+dv)
    end
    l = range(-120,stop=120,length=210)
    yolocar_vels = conv(yolocar_vels,exp.(-l.*l/1000))/200.0*18.0/3.7
    yolocar_vels = yolocar_vels[1000:end-1000]

    # scaling yolo car vel between max and min:
    yolocar_vels .-= minimum(yolocar_vels)
    yolocar_vels ./= maximum(yolocar_vels)
    yolocar_vels .*= (YOLOCAR_MAX_VEL - YOLOCAR_MIN_VEL)
    yolocar_vels .+= YOLOCAR_MIN_VEL

    yolocar_timestep = 1
    times_passed::Int32 = 0
    return MyState(av,yolocar_dist,yolocar_timestep, times_passed), MyEnv(col,d[end],splx,sply,yolocar_width,yolocar_height,yolocar_vels, triangle_points)
end

function initialize_simple_road(v=-1.0)
    width = 4.9 + rand()*2.1
    # a very long route (100km)
    xr,yr,xl,yl, xc,yc, a = gen_simple_road(round(Int64,100000/5), width, 500000.0/5, 0.1, 5.7, -1.0);
    return initialize_simple_road(xl,yl,xr,yr,xc,yc,v)
end

function iscolliding(av::Car, borders::Collider, safety_fac=1.3)
    for seg in getsegs(borders,av.cx,av.cy)
        test_coll(av, seg, safety_fac) && return true
    end
    return false
end

function my_calc_reward(av::Car, borders::Collider, yolo_x, yolo_y, safety=1.3)
    av.rcx < 0.0 && return 0.0, 1 # don't let it escape the road by reversing
    av.speed < MINSPEED && return 0.0, 1
    av.speed > MAXSPEED && return 0.0, 1
    iscolliding(av, borders, safety) && return 0.0, 2
    dx = av.cx - yolo_x
    dy = av.cy - yolo_y
    d2 = dx*dx + dy*dy
    d2 < 1.5*1.5*safety*safety && return 0.0, 3 # colliding with yolo car
    _,_,d2_1,_,_ = get_beam_coll(av.fcx,av.fcy, getsegs(borders,av.cx,av.cy), -pi/2+av.angle, 50.0)
    _,_,d2_2,_,_ = get_beam_coll(av.fcx,av.fcy, getsegs(borders,av.cx,av.cy),  pi/2+av.angle, 50.0)
    l = sqrt(d2_1)
    r = sqrt(d2_2)
    return (2.0*(av.speed>0.0) + max(0.0, av.speed)^1.2)*(min(l,r)/max(l,r))^1.5, 0
end

function calc_reward(av::Car, borders::Collider, yolo_x, yolo_y, safety=1.3)
    reward, result = my_calc_reward(av, borders, yolo_x, yolo_y, safety)
    return reward, result>0
end

getcar(state::MyState) =
    Car(state.rcx,state.rcy,state.angle,state.steer,state.accel,state.speed)

function get_sensors(env::MyEnv, state::MyState)
    car = getcar(state)
    spd = max(0.0, car.speed)
    bkg = floor(UInt8, 200.0*spd/MAXSPEED)
    img = ones(UInt8, 84,84) * bkg
    c = cos(-car.angle+pi/2)
    s = sin(-car.angle+pi/2)
    # plotting the road
    coll = env.coll
    segs = getsegs(coll, car.cx, car.cy)
    for seg in segs
        seg.ix > length(env.triangle_points) && continue
        # move and rotate the segment...
        x1 = seg.x1 - car.cx
        y1 = seg.y1 - car.cy
        x1,y1 = x1*c-y1*s, x1*s+y1*c
        x2 = seg.x2 - car.cx
        y2 = seg.y2 - car.cy
        x2,y2 = x2*c-y2*s, x2*s+y2*c
        xt = env.triangle_points[seg.ix][1] - car.cx
        yt = env.triangle_points[seg.ix][2] - car.cy
        xt,yt = xt*c-yt*s, xt*s+yt*c

        # to image coordinates:
        x1 = floor(Int64, x1/DRAWING_SCALE) + 42
        x2 = floor(Int64, x2/DRAWING_SCALE) + 42
        xt = floor(Int64, xt/DRAWING_SCALE) + 42
        y1 = floor(Int64, y1/DRAWING_SCALE)
        y2 = floor(Int64, y2/DRAWING_SCALE)
        yt = floor(Int64, yt/DRAWING_SCALE)

        # plot!
        draw_triangle!(img, xt,yt, x2,y2, x1,y1, 255)
        line!(img, x1,y1, x2,y2, 255)
    end

    # drawing the yolo-car
    xc = env.splx(state.yolocar_dist)
    yc = env.sply(state.yolocar_dist)

    # move and rotate
    xc -= car.cx
    yc -= car.cy
    xc,yc = xc*c-yc*s, xc*s+yc*c

    # to image coordinates:
    xc = floor(Int64, xc/DRAWING_SCALE) + 42
    yc = floor(Int64, yc/DRAWING_SCALE)

    # plot!
    maxdv = YOLOCAR_MAX_VEL - MINSPEED
    mindv = YOLOCAR_MIN_VEL - MAXSPEED

    spd = env.yolocar_vels[state.yolocar_timestep]-car.speed
    spd -= mindv
    spd /= maxdv - mindv
    spd += 0.1
    spd *= 200.0
    if spd < 0.0
        @show spd
        @show mindv
        @show maxdv
        @show env.yolocar_vels[state.yolocar_timestep]
        @show car.speed
    end
    yolo_color = floor(UInt8, spd)
    for x in (xc-2):(xc+2), y in (yc-2):(yc+2)
        x < 1 && continue
        y < 1 && continue
        x > 84 && continue
        y > 84 && continue
        img[y,x] = yolo_color
    end

    img[1,39:45] .= bkg
    img[2,40:44] .= bkg
    img[3,41:43] .= bkg
    img[4,42:42] .= bkg

    # plot current steering
    steer_color = floor(UInt8, (0.05+pi/30 + car.steer)/(pi/15) * 200.0)
    for y in 1:84
        img[y,1:4] .= steer_color
        img[y,5] = 228
        img[y,81:84] .= steer_color
        img[y,80] = 228
    end

    return img
end

function sim!(env::MyEnv, state::MyState, action::Tuple{Float64,Float64}, stepc=5, safety=SAFETY)
    reward = 0.0
    done = false
    steer, accel = action
    av = getcar(state)
    yolocar_dist::Float32 = state.yolocar_dist
    yolocar_timestep::Int32 = state.yolocar_timestep
    times_passed::Int32 = state.times_passed_yolocar
    dt = 0.1
    for j in 1:stepc
        # moving the yolo car!
        yolocar_timestep += 1
        yolocar_dist += dt*env.yolocar_vels[yolocar_timestep]

        av = move(av,dt,steer,accel)

        # check if we passed the yolo car
        # (if yes, we should spawn a new car)
        c = cos(-av.angle+pi/2)
        s = sin(-av.angle+pi/2)
        xc = env.splx(state.yolocar_dist)
        yc = env.sply(state.yolocar_dist)
        xc -= av.cx
        yc -= av.cy
        if xc*xc + yc*yc < 100.0*100.0
            # we are cloose to the car, passing it was possible within a short timescale
            xc,yc = xc*c-yc*s, xc*s+yc*c
            xc = floor(Int64, xc/DRAWING_SCALE) + 42
            yc = floor(Int64, yc/DRAWING_SCALE)
            if yc < 2
                # passed, will no be drawn on screen
                yolocar_dist += 70 + randn()*7
                times_passed += 1
            end
        end

        # calculating reward
        yolo_x::Float32 = env.splx(yolocar_dist)
        yolo_y::Float32 = env.sply(yolocar_dist)
        imre, done = calc_reward(av, env.coll, yolo_x, yolo_y, safety)
        reward += imre
        done && break
    end
    if done
        reward = 0.0
    end
    state = MyState(av,yolocar_dist,yolocar_timestep, times_passed)
    return state, reward, done
end


text:
using StackOverflow

"""
    getopenquestions(page::Int = 1, pagesize::Int = 1, fromdate::String = "", todate::String = "",
                            order::String = "desc", sort::String = "activity", tags::String = "Julia",
                            site::String = "stackoverflow")

Will return questions that have no marked answer given a tag. `pagesize` is the number of questions it will return,
`page` is the page on which it takes questions from, and returns `pagesize` number of questions.

Possible order values: "desc" or "asc"
Possible sort values: "activity", "votes", "creation", "hot", "week", or "month"
"""
function getopenquestions(page::Int = 1, pagesize::Int = 1, fromdate::String = "", todate::String = "",
                                order::String = "desc", sort::String = "activity", tags::String = "Julia",
                                site::String = "stackoverflow")
    tags = replace(tags, " " => "%20")

    if fromdate != "" && todate != ""
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions/unanswered?page=$page&pagesize=$pagesize&fromdate=$fromdate&todate=$todate&order=$order&sort=$sort&tagged=$tags&site=$site")
    else
        r = HTTP.request("GET", "https://api.stackexchange.com/2.2/questions/unanswered?page=$page&pagesize=$pagesize&order=$order&sort=$sort&tagged=$tags&site=$site")
    end

    json = StackOverflow.convert_HTTP_Response_To_JSON(r)

    questionholder = Vector()
    for (k, v) in json
        if occursin("items", k)
            for item in v
                link = ""
                view_count = 0
                creation_date = 0
                is_answered = false
                owner = Dict()
                last_activity_date = 0
                score = 0
                accepted_answer_id = 0
                question_id = 0
                title = 0
                answer_count = 0
                tags = []
                title = ""

                for (key, value) in item
                    if occursin("link", key)
                        link = value
                    elseif occursin("view_count", key)
                        view_count = value
                    elseif occursin("creation_date", key)
                        creation_date = value
                    elseif occursin("is_answered", key)
                        is_answered = value
                    elseif occursin("owner", key)
                        owner = value
                    elseif occursin("last_activity_date", key)
                        last_activity_date = value
                    elseif occursin("score", key)
                        score = value
                    elseif occursin("accepted_answer_id", key)
                        accepted_answer_id = value
                    elseif occursin("question_id", key)
                        question_id = value
                    elseif occursin("tags", key)
                        tags = value
                    elseif occursin("title", key)
                        title = value
                    elseif occursin("answer_count", key)
                        answer_count = value
                    end
                    # println(key)
                end
                question = StackOverflow.Question(link, view_count, creation_date, is_answered, owner,
                    last_activity_date, score, accepted_answer_id, question_id, tags, title,
                    answer_count)

                push!(questionholder, question)
            end

        end
    end
    return questionholder
end