text:
dottosym(x) = x
dottosym(x::Expr) = Symbol(dottosym(x.args[1]), "###extractarray###", x.args[2].value)
function extract_array_symbol_from_ref!(ls::LoopSet, ex::Expr, offset1::Int)::Symbol
    ar = ex.args[1 + offset1]
    if isa(ar, Symbol)
        return ar
    elseif isa(ar, Expr) && ar.head === :(.)
        s = dottosym(ar)
        pushprepreamble!(ls, Expr(:(=), s, ar))
        return s
    else
        throw("Indexing into the following expression was not recognized: $ar")
    end
end

function ref_from_expr!(ls, ex, offset1::Int, offset2::Int)
    ar = extract_array_symbol_from_ref!(ls, ex, offset1)
    ar, @view(ex.args[2 + offset2:end])
end
ref_from_ref!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 0, 0)
ref_from_getindex!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 1, 1)
ref_from_setindex!(ls::LoopSet, ex::Expr) = ref_from_expr!(ls, ex, 1, 2)
function ref_from_expr!(ls::LoopSet, ex::Expr)
    if ex.head === :ref
        ref_from_ref!(ls, ex)
    else#if ex.head === :call
        f = first(ex.args)::Symbol
        f === :getindex ? ref_from_getindex!(ls, ex) : ref_from_setindex!(ls, ex)
    end
end

add_vptr!(ls::LoopSet, op::Operation) = add_vptr!(ls, op.ref)
add_vptr!(ls::LoopSet, mref::ArrayReferenceMeta) = add_vptr!(ls, mref.ref.array, vptr(mref))
using VectorizationBase: noaliasstridedpointer
function add_vptr!(ls::LoopSet, array::Symbol, vptrarray::Symbol, actualarray::Bool = true, broadcast::Bool = false)
    if !includesarray(ls, array)
        push!(ls.includedarrays, array)
        actualarray && push!(ls.includedactualarrays, array)
        if broadcast
            pushpreamble!(ls, Expr(:(=), vptrarray, Expr(:call, lv(:stridedpointer_for_broadcast), array)))
        else
            pushpreamble!(ls, Expr(:(=), vptrarray, Expr(:call, lv(:stridedpointer), array)))
            # pushpreamble!(ls, Expr(:(=), vptrarray, Expr(:call, lv(:noaliasstridedpointer), array)))
        end
    end
    nothing
end

@inline valsum() = Val{0}()
@inline valsum(::Val{M}) where {M} = Val{M}()
@generated valsum(::Val{M}, ::Val{N}) where {M,N} = Val{M+N}()
@inline valsum(::Val{M}, ::Val{N}, ::Val{K}, args...) where {M,N,K} = valsum(valsum(Val{M}(), Val{N}()), Val{K}(), args...)
@inline valdims(::Any) = Val{1}()
@inline valdims(::CartesianIndices{N}) where {N} = Val{N}()

function append_loop_valdims!(valcall::Expr, loop::Loop)
    if isstaticloop(loop)
        push!(valcall.args, :(Val{1}()))
    else
        push!(valcall.args, Expr(:call, lv(:valdims), loop_boundary(loop)))
    end
    nothing
end
function subset_vptr!(ls::LoopSet, vptr::Symbol, indnum::Int, ind, previndices, loopindex)
    subsetvptr = Symbol(vptr, "_subset_$(indnum)_with_$(ind)##")
    valcall = Expr(:call, Expr(:curly, :Val, 1))
    if indnum > 1
        offset = first(previndices) === DISCONTIGUOUS
        valcall = Expr(:call, lv(:valsum), valcall)
        for i ∈ 1:indnum-1
            if loopindex[i]
                append_loop_valdims!(valcall, getloop(ls, previndices[i+offset]))
            else
                for loopdep ∈ loopdependencies(ls.opdict[previndices[i+offset]])
                    append_loop_valdims!(valcall, getloop(ls, loopdep))
                end
            end
        end
    end
    # @show valcall
    indm1 = ind isa Integer ? ind - 1 : Expr(:call, :-, ind, 1)
    pushpreamble!(ls, Expr(:(=), subsetvptr, Expr(:call, lv(:subsetview), vptr, valcall, indm1)))
    subsetvptr
end

function addoffset!(ls, indices, offsets, loopedindex, loopdependencies, ind, offset)
    (typemin(Int8) ≤ offset ≤ typemax(Int8)) || return false
    push!(indices, ind);
    push!(offsets, offset % Int8)
    push!(loopedindex, true)
    push!(loopdependencies, ind)
    true
end
function addoffsetexpr!(ls, parents, indices, offsets, loopedindex, loopdependencies, reduceddeps, ind, offset, elementbytes)
    (typemin(Int8) ≤ offset ≤ typemax(Int8)) || return false
    parent = add_operation!(ls, gensym(:indexpr), ind, elementbytes, length(ls.loopsymbols))
    pushparent!(parents, loopdependencies, reduceddeps, parent)
    push!(indices, name(parent));
    push!(offsets, offset % Int8)
    push!(loopedindex, false)
    true
end

function checkforoffset!(
    ls::LoopSet, parents::Vector{Operation}, indices::Vector{Symbol}, offsets::Vector{Int8}, loopedindex::Vector{Bool},
    loopdependencies::Vector{Symbol}, reduceddeps::Vector{Symbol}, ind::Expr, elementbytes::Int
)
    ind.head === :call || return false
    f = first(ind.args)
    (((f === :+) || (f === :-)) && (length(ind.args) == 3)) || return false
    factor = f === :+ ? 1 : -1
    arg1 = ind.args[2]
    arg2 = ind.args[3]
    if arg1 isa Integer && isone(factor) # we want to return false when we're subtracting the index, e.g. A[3 - i]
        if arg2 isa Symbol && arg2 ∈ ls.loopsymbols
            addoffset!(ls, indices, offsets, loopedindex, loopdependencies, arg2, arg1 * factor)
        elseif arg2 isa Expr
            addoffsetexpr!(ls, parents, indices, offsets, loopedindex, loopdependencies, reduceddeps, arg2, arg1 * factor, elementbytes)
        else
            false
        end
    elseif arg2 isa Integer
        if arg1 isa Symbol && arg1 ∈ ls.loopsymbols
            addoffset!(ls, indices, offsets, loopedindex, loopdependencies, arg1, arg2 * factor)
        elseif arg1 isa Expr
            addoffsetexpr!(ls, parents, indices, offsets, loopedindex, loopdependencies, reduceddeps, arg1, arg2 * factor, elementbytes)
        else
            false
        end
    else
        false
    end
end

function move_to_last!(x, i)
    i == length(x) && return
    xᵢ = x[i]
    deleteat!(x, i)
    push!(x, xᵢ)
    nothing
end
# TODO: Make this work with Cartesian Indices
function repeated_index!(ls::LoopSet, indices::Vector{Symbol}, vptr::Symbol, indnum::Int, firstind::Int)
    # Move ind to last position
    vptrrepremoved = Symbol(vptr, "##ind##", firstind, "##repeated##", indnum, "##")
    f = Expr(:(.), Expr(:(.), :LoopVectorization, QuoteNode(:VectorizationBase)), QuoteNode(:double_index))
    fiv = Expr(:call, Expr(:curly, :Val, firstind - 1))
    siv = Expr(:call, Expr(:curly, :Val, indnum - 1))
    pushpreamble!(ls, Expr(:(=), vptrrepremoved, Expr(:call, f, vptr, fiv, siv)))
    vptrrepremoved
end

function array_reference_meta!(ls::LoopSet, array::Symbol, rawindices, elementbytes::Int, var::Union{Nothing,Symbol} = nothing)
    vptrarray = vptr(array)
    add_vptr!(ls, array, vptrarray) # now, subset
    indices = Symbol[]
    offsets = Int8[]
    loopedindex = Bool[]
    parents = Operation[]
    loopdependencies = Symbol[]
    reduceddeps = Symbol[]
    loopset = ls.loopsymbols
    ninds = 1
    for ind ∈ rawindices
        if ind isa Integer # subset
            vptrarray = subset_vptr!(ls, vptrarray, ninds, ind, indices, loopedindex)
            length(indices) == 0 && push!(indices, DISCONTIGUOUS)
        elseif ind isa Expr
            #FIXME: position (in loopnest) wont be length(ls.loopsymbols) in general
            if !checkforoffset!(ls, parents, indices, offsets, loopedindex, loopdependencies, reduceddeps, ind, elementbytes)
                parent = add_operation!(ls, gensym(:indexpr), ind, elementbytes, length(ls.loopsymbols))
                pushparent!(parents, loopdependencies, reduceddeps, parent)
                push!(indices, name(parent));
                push!(offsets, zero(Int8))
                push!(loopedindex, false)
            end
            ninds += 1
        elseif ind isa Symbol
            if ind ∈ loopset
                ind_prev_index = findfirst(isequal(ind), indices)
                if isnothing(ind_prev_index)
                    push!(indices, ind); ninds += 1
                    push!(offsets, zero(Int8))
                    push!(loopedindex, true)
                    push!(loopdependencies, ind)
                else
                    move_to_last!(indices, ind_prev_index)
                    move_to_last!(offsets, ind_prev_index)
                    move_to_last!(loopedindex, ind_prev_index)
                    move_to_last!(loopdependencies, ind_prev_index)
                    vptrarray = repeated_index!(ls, indices, vptrarray, ninds, ind_prev_index + (first(indices) === DISCONTIGUOUS))
                    makediscontiguous!(indices)
                end
            else
                indop = get(ls.opdict, ind, nothing)
                if indop !== nothing  && !isconstant(indop)
                    pushparent!(parents, loopdependencies, reduceddeps, indop)
                    push!(indices, name(indop)); ninds += 1
                    push!(offsets, zero(Int8))
                    push!(loopedindex, false)
                else
                    vptrarray = subset_vptr!(ls, vptrarray, ninds, ind, indices, loopedindex)
                    length(indices) == 0 && push!(indices, DISCONTIGUOUS)
                end
            end
        else
            throw("Unrecognized loop index: $ind.")
        end
    end
    # (length(parents) != 0 && first(indices) !== Symbol("##DISCONTIGUOUSSUBARRAY##")) && pushfirst!(indices, Symbol("##DISCONTIGUOUSSUBARRAY##"))
    mref = ArrayReferenceMeta(ArrayReference( array, indices, offsets ), loopedindex, vptrarray)
    ArrayReferenceMetaPosition(mref, parents, loopdependencies, reduceddeps, isnothing(var) ? Symbol("") : var )
end
function tryrefconvert(ls::LoopSet, ex::Expr, elementbytes::Int, var::Union{Nothing,Symbol} = nothing)::Tuple{Bool,ArrayReferenceMetaPosition}
    ya, yinds = if ex.head === :ref
        ref_from_ref!(ls, ex)
    elseif ex.head === :call
        f = first(ex.args)
        if f === :getindex
            ref_from_getindex!(ls, ex)
        elseif f === :setindex!
            ref_from_setindex!(ls, ex)
        else
            return false, NOTAREFERENCEMP
        end
    else
        return false, NOTAREFERENCEMP
    end
    true, array_reference_meta!(ls, ya, yinds, elementbytes, var)
end


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function matchmom(μ,s,av,va, minvar,maxvar)
    newb = clamp(inv(1.0/va - 1.0/s),minvar,maxvar)
    newa = av + newb*(av-μ)/s
    isnan(newa) || isnan(newb) && @warn("a = $newa b = $newb")
    return newa, newb
end

ext: jl


text:
# game of life kinda thing: 5-tile rules.

function parseinput(f)
    initial_input = match(r"[#.]+", readline(f)).match
    readline(f)
    (collect(initial_input),
     Dict(collect(l[1]) => l[2][1]
          for l
          in split.(strip.(readlines(f)), Ref(" => "))))
end

test_start, test_rules = parseinput(IOBuffer("""
initial state: #..#.#..##......###...###

...## => #
..#.. => #
.#... => #
.#.#. => #
.#.## => #
.##.. => #
.#### => #
#.#.# => #
#.### => #
##.#. => #
##.## => #
###.. => #
###.# => #
####. => #
"""))

start, rules = open(parseinput, "12.input")

function pad!(state, starti; n=3)
    # need at least 3 '.' on the front and  on the end
    front_pad = max(0, n + 1 - findfirst(isequal('#'), state))
    starti -= front_pad
    for _ in 1:front_pad
        pushfirst!(state, '.')
    end
    back_pad = n-(length(state)-findlast(isequal('#'), state))
    for _ in 1:back_pad
        push!(state, '.')
    end
    state, starti
end

function generation!(next_state, state, rules)
    for i in 3:(lastindex(state)-2)
        next_state[i] = get(rules, view(state, i-2:i+2), '.')
    end
    next_state
end

function grow(start, rules; n=20)
    state, starti = pad!(copy(start), 0, n=n+3)
    next_state = fill!(similar(state), '.')
    println(String(state))
    for gen in 1:n
        state, next_state = generation!(next_state, state, rules), fill!(state, '.')
        println(String(state))
    end
    state, starti
end

star1(state, starti) = sum(findall(isequal('#'), state) .- 1 .+ starti)

star1(grow(test_start, test_rules; n=20)...)
star1(grow(start, rules; n=20)...)

# star 2: 50 billion generations!! there's gotta be some kind of shortcut...  if
# you look at the patterns visually it's clear that at some point (around
# generation 101) the game reaches a repeating state that's a bunch of "##..."
# repeated some number of time, and the whole thing moves to the right 1 square
# each time.  so need to just figure out the score for one generation and the
# number of live cells...

state100, starti100 = grow(start, rules; n=100)
score = star1(state100, starti100)
nlive = sum(isequal('#'), state100)

# score will increase by $nlive every generation...so score(t) = score(100) + (t-100)*nlive
score + (50_000_000_000-100)*nlive


text:
function removebracketspair(string)
    for brackets in ["()", "[]", "{}"]
        string = replace(string, brackets => "")
    end
    return string
end

function matching_brackets(string)
    string = replace(string, r"[^\[\](){}]" => "")
    tmp = removebracketspair(string)
    while length(tmp) < length(string)
        string = tmp
        tmp = removebracketspair(string)
    end
    return isempty(string)
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Metida
using  Test, CSV, DataFrames, StatsBase, PrettyTables, MixedModels, Distributions
path    = dirname(@__FILE__)


text:
using Test
using PersistenceDiagrams

@testset "Construction" begin
    for (start, stop) in ((0, 2), (ℯ, π)), len in (2, 13, 1001)
        curve = PersistenceCurve(identity, sum, start, stop; length=len)

        @testset "moving from `start` to `stop` by `step` yields `len` steps" begin
            t = curve.start
            count = 0
            while Float64(t) < curve.stop
                t += curve.step
                count += 1
            end
            @test count == len
        end

        @testset "show" begin
            @test sprint(show, curve) ==
                "PersistenceCurve(identity, sum, $(Float64(start)), $(Float64(stop)); " *
                  "length=$len, normalize=false, integrate=true)"
        end
    end

    @testset "ranges are selected correctly when learning parameters" begin
        diagram_1 = PersistenceDiagram([(0, 1), (0.5, 1), (0.5, 0.6), (1, 1.5), (0.5, Inf)])
        diagram_2 = PersistenceDiagram([(2, Inf)])

        @testset "for diagrams with infinite intervals" begin
            bc = PersistenceCurve((_...) -> 1.0, sum, [diagram_1, diagram_2])
            @test bc.start == 0
            @test bc.stop == Float64(2 + bc.step)
            @test bc.length == 10
            @test length(bc(diagram_1)) == 10
            @test length(bc(diagram_2)) == 10
        end

        @testset "for diagrams without infinite intervals" begin
            finite = [filter(isfinite, diagram_1), filter(isfinite, diagram_2)]
            bc = PersistenceCurve((_...) -> 1.0, sum, finite)
            @test bc.start == 0
            @test bc.stop == 1.5
            @test bc.length == 10
            @test length(bc(diagram_1)) == 10
            @test length(bc(diagram_2)) == 10
            @test all(iszero(bc(diagram_2)))
        end
    end
end

@testset "BettiCurve" begin
    @testset "With basic constructor" begin
        diagram = PersistenceDiagram([(0, 1), (0.5, 1), (0.5, 0.6), (1, 1.5), (0.5, Inf)])
        bc_1 = BettiCurve(0, 2; length=2)
        @test bc_1(diagram) == [2.1, 1.5]

        bc_2 = BettiCurve(0, 2; length=20)
        @test bc_2(diagram) == [1, 1, 1, 1, 1, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1]

        bc_3 = BettiCurve(0, 2; length=10)
        @test bc_3(diagram) == [1, 1, 2.5, 3, 3, 2, 2, 1.5, 1.0, 1.0]
    end

    @testset "With learned parameters" begin
        diagram_1 = PersistenceDiagram([(0, 1), (0.5, 1), (0.5, 0.6), (1, 1.5), (0.5, Inf)])
        diagram_2 = PersistenceDiagram([(2, Inf)])
        @testset "Last bucket only contains infinite intervals." begin
            for len in (10, 20, 40)
                bc = BettiCurve([diagram_1]; length=len)
                @test bc(diagram_1)[end] == 1
                @test bc(diagram_1)[end - 1] == 2
                @test bc.length == len
                @test length(bc(diagram_1)) == len
                @test length(bc(diagram_2)) == len
            end
        end
    end
end

@testset "Landscape(s) and Slihuette" begin
    diagram = PersistenceDiagram([(3, 9), (4, 6), (5, 11)])

    @testset "(k + 1)-th landscape is below the k-th landscape" begin
        for k in 1:3
            scape_1 = Landscape(k, 0, 12; length=12)
            scape_2 = Landscape(k + 1, 0, 12; length=12)
            @test all(scape_2(diagram) .≤ scape_1(diagram))
        end
    end

    @testset "for a large enough k, a landscape is zero everywhere" begin
        scape = Landscape(4, 0, 12; length=12)
        @test all(iszero, scape(diagram))
    end

    @testset "columns in Landscapes are equal to Landscape" begin
        scapes = Landscapes(4, 0, 12; length=12)
        result = scapes(diagram)
        @test size(result) == (12, 4)
        for i in 1:4
            scape = Landscape(i, 0, 12; length=12)
            @test result[:, i] == scape(diagram)
        end

        @test sprint(show, scapes) == "Landscapes(4, 0.0, 12.0; length=12)"
    end

    @testset "a silhuette is equivalent to the sum of landscapes" begin
        scape_res = zeros(24)
        for k in 1:3
            scape_res .+= Landscape(k, 0, 12; length=24)(diagram)
        end
        @test Silhuette(0, 12; length=24)(diagram) == scape_res
    end

    @test_throws ArgumentError Landscape(-1, 0, 12)
    @test_throws ArgumentError Landscapes(-1, 0, 12)
end

@testset "sum-based time independent curves" begin
    diagram = PersistenceDiagram([(1, 3), (2, 4.5)])

    for constructor in (BettiCurve, Life, Midlife, LifeEntropy, MidlifeEntropy)
        @testset "values are correct for $(nameof(constructor))" begin
            curve = constructor(0, 5; length=5)
            result = curve(diagram)
            x1 = curve.fun(diagram[1], diagram, nothing)
            x2 = curve.fun(diagram[2], diagram, nothing)
            @test result[1] == 0.0
            @test result[2] == x1
            @test result[3] == x1 + x2
            @test result[4] == x2
            @test result[5] == x2 / 2
        end
    end
end

@testset "PDThresholding" begin
    @testset "is always positive" begin
        diagram = PersistenceDiagram([(3, 9), (4, 6), (5, 11)])
        @test all(PDThresholding([diagram])(diagram) .≥ 0)
    end

    @testset "is highest in the middle for a long interval" begin
        diagram = PersistenceDiagram([(ℯ, π)])
        result = PDThresholding([diagram]; length=11)(diagram)
        @test all(result[end ÷ 2 + 1] .≥ result)
    end

    @testset "is unchanged with multiple equal intervals" begin
        diagram_1 = PersistenceDiagram([(ℯ, π)])
        diagram_2 = PersistenceDiagram([(ℯ, π), (ℯ, π)])
        pdt = PDThresholding([diagram_1, diagram_2])
        @test pdt(diagram_1) == pdt(diagram_2)
    end
end

@testset "Normalization" begin
    diagram = PersistenceDiagram([
        (0.0, 1.0),
        (0.0, 1.0),
        (ℯ, π),
        (1.0, 1.5),
        (1.1, 1.6),
        (1.2, 1.7),
        (1.3, 1.8),
        (1.0, ℯ),
        (1.0, π),
        (√5, ℯ),
    ])

    for constructor in (BettiCurve, Life, Midlife, LifeEntropy, MidlifeEntropy)
        @testset "$(nameof(constructor))" begin
            curve_reg = constructor([diagram])
            curve_norm = constructor([diagram]; normalize=true)
            @testset "all values in normalized curve are below 1" begin
                @test all(curve_norm(diagram) .≤ 1)
            end
        end
    end
    for constructor in (Landscape, Landscapes, Silhuette, PDThresholding)
        @testset "$(nameof(constructor))" begin
            @test_throws MethodError constructor([diagram]; normalize=true)
        end
    end
end


text:
function free_energy(xs::Vector)
  I1 = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*u(θs),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  I2 = pcubature(θs -> ρ(θs, xs)*sin(θs[2])*log(ρ(θs, xs)),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  return I1[1] + kT*I2[1] - kT*N*log(N);
end

function free_energy_smallomega(x::Vector)
  c, λ, α = x[:];
  I1 = pcubature(θs -> c*(1 - ω(θs) + α*cos(θs[1])*sin(θs[2]))*exp(λ*cos(θs[2]))*sin(θs[2])*u(θs),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  I2 = pcubature(θs -> c*(1 - ω(θs) + α*cos(θs[1])*sin(θs[2]))*exp(λ*cos(θs[2]))*sin(θs[2])*log(ρ(θs, xs)),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  return I1[1] + kT*I2[1] - kT*N*log(N);
end

function free_energy_smalllambda(x::Vector)
  c, λ, α = x[:];
  I1 = pcubature(θs -> c*(1 + λ*cos(θs[2]) + α*cos(θs[1])*sin(θs[2]))*exp(-ω(θs))*sin(θs[2])*u(θs),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
                 I2 = pcubature(θs -> c*(1 + λ*cos(θs[2]) + α*cos(θs[1])*sin(θs[2]))*exp(-ω(θs))*sin(θs[2])*log(ρ(θs, xs)),
                 _LOWER_BOUNDS_, _UPPER_BOUNDS_; reltol=RELTOL_INT, maxevals=MAXEVALS_INT);
  return I1[1] + kT*I2[1] - kT*N*log(N);
end

max_stars_repo_name: grasingerm/polymer-stats

max_issues_repo_name: grasingerm/polymer-stats

max_forks_repo_name: grasingerm/polymer-stats


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

@testset ExtendedTestSet "DocumentFolder" begin
    dir = Path(mktempdir())
    try
        write(joinpath(dir, "hi.md"), "Hello")
        rewriter = DocumentFolder(dir)

        docs = createsources!(rewriter)
        @test haskey(docs, p"hi.md")
        @test tag(docs[p"hi.md"]) == :body

        ch = Channel()
        fw = geteventsource(rewriter, ch)
        Pollen.start(fw)
        open(joinpath(dir, p"hi.md"), "w") do f
            write(f, "Hello")
        end
        event = take!(ch)
        @test event.name == p"hi.md"
        @test Pollen.gettext(event.doc) == "Hello"
        Pollen.stop(fw)
    catch e
        rethrow(e)
    finally
        rm(dir; recursive=true)
    end
end


text:
# Contiguous on-device arrays

export CuDeviceArray, CuDeviceVector, CuDeviceMatrix, ldg

## construction

"""
    CuDeviceArray(dims, ptr)
    CuDeviceArray{T}(dims, ptr)
    CuDeviceArray{T,N}(dims, ptr)
    CuDeviceArray{T,N,A}(dims, ptr)

Construct an `N`-dimensional dense CUDA device array with element type `T` wrapping a
pointer, where `N` is determined from the length of `dims` and `T` is determined from the
type of `ptr`. `dims` may be a single scalar, or a tuple of integers corresponding to the
lengths in each dimension). If the rank `N` is supplied explicitly as in `Array{T,N}(dims)`,
then it must match the length of `dims`. The same applies to the element type `T`, which
should match the type of the pointer `ptr`.
"""
CuDeviceArray

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're currently requiring a trailing pointer argument.

struct CuDeviceArray{T,N,A} <: DenseArray{T,N}
    ptr::LLVMPtr{T,A}
    maxsize::Int

    dims::Dims{N}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    # TODO: deprecate; put `ptr` first like CuArray
    CuDeviceArray{T,N,A}(dims::Dims{N}, ptr::LLVMPtr{T,A},
                         maxsize::Int=prod(dims)*sizeof(T)) where {T,A,N} =
        new(ptr, maxsize, dims)
end

const CuDeviceVector = CuDeviceArray{T,1,A} where {T,A}
const CuDeviceMatrix = CuDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
CuDeviceArray(dims::NTuple{N,Integer}, p::LLVMPtr{T,A})        where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceVector(dims::NTuple{1,Integer}, p::LLVMPtr{T,A})       where {T,A}   = CuDeviceVector{T,A}(dims, p)
CuDeviceMatrix(dims::NTuple{2,Integer}, p::LLVMPtr{T,A})       where {T,A}   = CuDeviceMatrix{T,A}(dims, p)
CuDeviceArray(len::Integer,            p::LLVMPtr{T,A})        where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceArray(m::Integer, n::Integer,  p::LLVMPtr{T,A})        where {T,A}   = CuDeviceMatrix{T,A}((m,n), p)
CuDeviceVector(len::Integer,           p::LLVMPtr{T,A})        where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceMatrix(m::Integer, n::Integer, p::LLVMPtr{T,A})        where {T,A}   = CuDeviceMatrix{T,A}((m,n), p)

# outer constructors, partially parameterized
CuDeviceArray{T}(dims::NTuple{N,Integer},     p::LLVMPtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceVector{T}(dims::NTuple{1,Integer},    p::LLVMPtr{T,A}) where {T,A}   = CuDeviceVector{T,A}(dims, p)
CuDeviceMatrix{T}(dims::NTuple{2,Integer},    p::LLVMPtr{T,A}) where {T,A}   = CuDeviceMatrix{T,A}(dims, p)
CuDeviceArray{T}(len::Integer,                p::LLVMPtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceArray{T}(m::Integer, n::Integer,      p::LLVMPtr{T,A}) where {T,A}   = CuDeviceMatrix{T,A}((m,n), p)
CuDeviceArray{T,N}(dims::NTuple{N,Integer},   p::LLVMPtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(dims, p)
CuDeviceVector{T}(len::Integer,               p::LLVMPtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceMatrix{T}(m::Integer, n::Integer,     p::LLVMPtr{T,A}) where {T,A}   = CuDeviceMatrix{T,A}((m,n), p)

# outer constructors, fully parameterized
CuDeviceArray{T,N,A}(dims::NTuple{N,Integer}, p::LLVMPtr{T,A}) where {T,A,N} = CuDeviceArray{T,N,A}(convert(Tuple{Vararg{Int}}, dims), p)
CuDeviceVector{T,A}(len::Integer,             p::LLVMPtr{T,A}) where {T,A}   = CuDeviceVector{T,A}((len,), p)
CuDeviceMatrix{T,A}(m::Integer, n::Integer,   p::LLVMPtr{T,A}) where {T,A}   = CuDeviceMatrix{T,A}((m,n), p)

## array interface

Base.elsize(::Type{<:CuDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::CuDeviceArray) = g.dims
Base.length(g::CuDeviceArray) = prod(g.dims)

Base.sizeof(x::CuDeviceArray) = Base.elsize(x) * length(x)

Base.pointer(x::CuDeviceArray{T,<:Any,A}) where {T,A} = Base.unsafe_convert(LLVMPtr{T,A}, x)
@inline function Base.pointer(x::CuDeviceArray{T,<:Any,A}, i::Integer) where {T,A}
    Base.unsafe_convert(LLVMPtr{T,A}, x) + Base._memory_offset(x, i)
end

typetagdata(a::CuDeviceArray{<:Any,<:Any,A}, i=1) where {A} =
  reinterpret(LLVMPtr{UInt8,A}, a.ptr + a.maxsize) + i - one(i)

## conversions

Base.unsafe_convert(::Type{LLVMPtr{T,A}}, x::CuDeviceArray{T,<:Any,A}) where {T,A} =
  x.ptr

## indexing intrinsics

# TODO: arrays as allocated by the CUDA APIs are 256-byte aligned. we should keep track of
#       this information, because it enables optimizations like Load Store Vectorization
#       (cfr. shared memory and its wider-than-datatype alignment)

@generated function alignment(::CuDeviceArray{T}) where {T}
    if Base.isbitsunion(T)
        _, sz, al = Base.uniontype_layout(T)
        al
    else
        Base.datatype_alignment(T)
    end
end

@inline function arrayref(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    if isbitstype(T)
        arrayref_bits(A, index)
    else #if isbitsunion(T)
        arrayref_union(A, index)
    end
end

@inline function arrayref_bits(A::CuDeviceArray{T}, index::Integer) where {T}
    align = alignment(A)
    unsafe_load(pointer(A), index, Val(align))
end

@inline @generated function arrayref_union(A::CuDeviceArray{T,<:Any,AS}, index::Integer) where {T,AS}
    typs = Base.uniontypes(T)

    # generate code that conditionally loads a value based on the selector value.
    # lacking noreturn, we return T to avoid inference thinking this can return Nothing.
    ex = :(Base.llvmcall("unreachable", $T, Tuple{}))
    for (sel, typ) in Iterators.reverse(enumerate(typs))
        ex = quote
            if selector == $(sel-1)
                ptr = reinterpret(LLVMPtr{$typ,AS}, data_ptr)
                unsafe_load(ptr, 1, Val(align))
            else
                $ex
            end
        end
    end

    quote
        selector_ptr = typetagdata(A, index)
        selector = unsafe_load(selector_ptr)

        align = alignment(A)
        data_ptr = pointer(A, index)

        return $ex
    end
end

@inline function arrayset(A::CuDeviceArray{T}, x::T, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    if isbitstype(T)
        arrayset_bits(A, x, index)
    else #if isbitsunion(T)
        arrayset_union(A, x, index)
    end
    return A
end

@inline function arrayset_bits(A::CuDeviceArray{T}, x::T, index::Integer) where {T}
    align = alignment(A)
    unsafe_store!(pointer(A), x, index, Val(align))
end

@inline @generated function arrayset_union(A::CuDeviceArray{T,<:Any,AS}, x::T, index::Integer) where {T,AS}
    typs = Base.uniontypes(T)
    sel = findfirst(isequal(x), typs)

    quote
        selector_ptr = typetagdata(A, index)
        unsafe_store!(selector_ptr, $(UInt8(sel-1)))

        align = alignment(A)
        data_ptr = pointer(A, index)

        unsafe_store!(reinterpret(LLVMPtr{$x,AS}, data_ptr), x, 1, Val(align))
        return
    end
end

@inline function const_arrayref(A::CuDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = alignment(A)
    unsafe_cached_load(pointer(A), index, Val(align))
end

## indexing

Base.IndexStyle(::Type{<:CuDeviceArray}) = Base.IndexLinear()

Base.@propagate_inbounds Base.getindex(A::CuDeviceArray{T}, i1::Integer) where {T} =
    arrayref(A, i1)
Base.@propagate_inbounds Base.setindex!(A::CuDeviceArray{T}, x, i1::Integer) where {T} =
    arrayset(A, convert(T,x)::T, i1)

# preserve the specific integer type when indexing device arrays,
# to avoid extending 32-bit hardware indices to 64-bit.
Base.to_index(::CuDeviceArray, i::Integer) = i

# Base doesn't like Integer indices, so we need our own ND get and setindex! routines.
# See also: https://github.com/JuliaLang/julia/pull/42289
Base.@propagate_inbounds Base.getindex(A::CuDeviceArray,
                                       I::Union{Integer, CartesianIndex}...) =
    A[Base._to_linear_index(A, to_indices(A, I)...)]
Base.@propagate_inbounds Base.setindex!(A::CuDeviceArray, x,
                                        I::Union{Integer, CartesianIndex}...) =
    A[Base._to_linear_index(A, to_indices(A, I)...)] = x

## const indexing

"""
    Const(A::CuDeviceArray)

Mark a CuDeviceArray as constant/read-only. The invariant guaranteed is that you will not
modify an CuDeviceArray for the duration of the current kernel.

This API can only be used on devices with compute capability 3.5 or higher.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N,AS} <: DenseArray{T,N}
    a::CuDeviceArray{T,N,AS}
end
Base.Experimental.Const(A::CuDeviceArray) = Const(A)

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
Base.@propagate_inbounds Base.getindex(A::Const, i1::Integer) = const_arrayref(A.a, i1)

# deprecated
Base.@propagate_inbounds ldg(A::CuDeviceArray, i1::Integer) = const_arrayref(A, i1)

## other

Base.show(io::IO, a::CuDeviceVector) =
    @printf(io, "%g-element device array at %p", length(a), Int(pointer(a)))
Base.show(io::IO, a::CuDeviceArray) =
    @printf(io, "%s device array at %p", join(a.dims, '×'), Int(pointer(a)))

Base.show(io::IO, mime::MIME"text/plain", a::CuDeviceArray) = show(io, a)

@inline function Base.iterate(A::CuDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end

function Base.reinterpret(::Type{T}, a::CuDeviceArray{S,N,A}) where {T,S,N,A}
  err = _reinterpret_exception(T, a)
  err === nothing || throw(err)

  if sizeof(T) == sizeof(S) # fast case
    return CuDeviceArray{T,N,A}(size(a), reinterpret(LLVMPtr{T,A}, a.ptr), a.maxsize)
  end

  isize = size(a)
  size1 = div(isize[1]*sizeof(S), sizeof(T))
  osize = tuple(size1, Base.tail(isize)...)
  return CuDeviceArray{T,N,A}(osize, reinterpret(LLVMPtr{T,A}, a.ptr), a.maxsize)
end


text:
using Test
using DataFrames
using PrettyTables

using ScottishTaxBenefitModel
using .FRSHouseholdGetter
using .ModelHousehold
using .BenefitGenerosity: initialise, to_set, change_status, adjust_disability_eligibility!
using .NonMeansTestedBenefits: calc_pip, calc_dla, calc_attendance_allowance
using .Definitions
using .RunSettings
using .STBIncomes

using .ExampleHelpers

settings = Settings()
ruksys = get_system( scotland = false )
scosys = get_system( scotland = true )

FRSHouseholdGetter.initialise( settings )

@testset "Loading Tests" begin

    # println( FRSHouseholdGetter.MODEL_HOUSEHOLDS.pers_map )
    ks = sort( collect( keys( FRSHouseholdGetter.MODEL_HOUSEHOLDS.pers_map )), lt=ModelHousehold.isless )
    for k in ks
        if k.data_year == 2015
            print( k )
            if k.id == 120150847501
                print("!!!!")
            end
            println()
        end
    end
    BenefitGenerosity.initialise( "$(MODEL_DATA_DIR)/disability/")

    for ben in [ATTENDANCE_ALLOWANCE, PERSONAL_INDEPENDENCE_PAYMENT_DAILY_LIVING,
        PERSONAL_INDEPENDENCE_PAYMENT_MOBILITY, DLA_SELF_CARE]
        for peeps in [-100_000, -10_000, -1000, 0, 1000, 10_000, 100_000 ]
            s = to_set( ben, peeps )
            println( "set for $peeps $ben = $(s)")
        end
    end
end

"""
Quicky thing to make a little set of targets for the testset below
"""
function makeset( r :: UnitRange ) :: Set
    s = Set{OneIndex}()
    for i in r
        push!( s, OneIndex( makePID(i),2018))
    end
    s
end

@enum Benvals noben lowben medben highben

@testset "Status Change Thingy" begin

    s = makeset(1:20)
    println( "s=$s")

    # in set, +ive => set to random low:high
    @test change_status(
        candidates = s,
        pid = makePID(1),
        change = 1,
        choices = [lowben, medben, highben],
        current_value = noben,
        disqual_value = noben ) in [lowben,medben,highben]
    # not in set - unchanged
    @test change_status(
        candidates = s,
        pid = makePID(100),
        change = 1,
        choices = [lowben, medben, highben],
        current_value = lowben,
        disqual_value = noben ) == lowben
    # in set - change -ive; disqual
    @test change_status(
        candidates = s,
        pid = makePID(10),
        change = -1,
        choices = [lowben, medben, highben],
        current_value = lowben,
        disqual_value = noben ) == noben
    # out of set - no change
    @test change_status(
        candidates = s,
        pid = makePID(100),
        change = -1,
        choices = [lowben, medben, highben],
        current_value = lowben,
        disqual_value = noben ) == lowben
    # no extras - unchanged even if in set
    @test change_status(
        candidates = s,
        pid = makePID(10),
        change = 0,
        choices = [lowben, medben, highben],
        current_value = lowben,
        disqual_value = noben ) == lowben
end

@testset "Add/Subtract from PIP aggregate tests" begin

    rc = @timed begin
        num_households,total_num_people,nhh2 = FRSHouseholdGetter.initialise( DEFAULT_SETTINGS )
    end
    println( "num_households=$num_households, num_people=$(total_num_people)")
    BenefitGenerosity.initialise( MODEL_DATA_DIR*"/disability/" )

    out = DataFrame(
        which = ["RUK","-50,000","0","50,000"],
        pip_d = zeros(4),
        pip_m = zeros(4),
        dla_d = zeros(4),
        dla_m = zeros(4),
        aa    = zeros(4)
    )
    changes = [0,-50_000,0,50_000]
    for sysno in 1:4
        outr = out[sysno,:]
        if sysno == 1
            sys = ruksys
        else
            sys = scosys
            sys.nmt_bens.pip.extra_people = changes[sysno]
            sys.nmt_bens.dla.extra_people = changes[sysno]
            sys.nmt_bens.attendance_allowance.extra_people = changes[sysno]

            adjust_disability_eligibility!( sys.nmt_bens )
            adjust_disability_eligibility!( sys.nmt_bens )
        end
        r = 0
        @time for hhno in 1:num_households
            hh = FRSHouseholdGetter.get_household( hhno )
            r += 1
            for (pid,pers) in hh.people
                (pip_d,pip_m) = calc_pip( pers, sys.nmt_bens.pip )
                if pip_d > 0
                    outr.pip_d += hh.weight
                end
                if pip_m > 0
                    outr.pip_m += hh.weight
                end
                (dla_d,dla_m) = calc_dla( pers, sys.nmt_bens.dla )
                if dla_d > 0
                    outr.dla_d += hh.weight
                end
                if dla_m > 0
                    outr.dla_m += hh.weight
                end
                aa = calc_attendance_allowance( pers, sys.nmt_bens.attendance_allowance )
                if aa > 0
                    outr.aa += hh.weight
                end
            end
        end
        # should be ~50k bigger/smaller for scottish system
    end
    pretty_table( out; formatters = ft_printf("%10.0f", [2,3,4,5,6]))
    dout =
    pretty_table( out; formatters = ft_printf("%10.0f", [2,3,4,5,6]))
end


text:
mutable struct _DrivenProblem{O, S <: WavePacketState, WP <: _WavePacket{S}} <: WaveguideProblem
	H::O
	Ls::Array{O}
	σ::O
	system_state

	ψᵢ::WP

	ts::StepRangeLen
end

function WaveguideProblem(sys, ts, disp_out::Bool = false)
	WaveguideProblem(sys, NoPacket(), ts, false)
end

function WaveguideProblem(sys, ψᵢ::_WavePacket, ts, disp_out::Bool = false)
	disp_out && @warn("Displacing the output ineffective without output modes")
	ts isa StepRangeLen || (ts = _toRange(ts))

	length(sys) == 3 && return _DrivenProblem(sys[1], typeof(sys[1])[], sys[2], sys[3], ψᵢ, ts)
	return _DrivenProblem(sys[1], sys[2], sys[3], sys[4], ψᵢ, ts)
end


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using EngEconomics, Plots, ColorSchemes, Colors
plotly()

uwEngColorScheme = ColorScheme([parse(Colorant, "#D0B4E7"), parse(Colorant, "#FFFFFF"), parse(Colorant, "#57058B"), parse(Colorant, "#8100B4")])

F = 50
A_base = 7
i_base = 0.10
sensitivityVector = [.90, 0.95, 1.0, 1.05, 1.10]

baseN = log(F / A_base * i_base + 1) / log(1 + i_base)

N_annualSavings = log.(F ./ (A_base .* sensitivityVector) .* i_base .+ 1) ./ log(1 + i_base)
N_interestRate  = log.(F / A_base .* (sensitivityVector * i_base) .+ 1) ./ log.(1 .+ (i_base .* sensitivityVector))

plot(sensitivityVector, N_annualSavings, label="Annual Savings Sensitivity",
	xlabel="Sensitivity (+/- 5%, 10%)", ylabel="Number of Years",
	palette=cgrad(uwEngColorScheme))
plot!(sensitivityVector, N_interestRate, label="Interest Rate Sensitivity")


text:
# General bayesian filtering and updates

using HistogramFilters
using LinearAlgebra
using Base.Iterators
using StatsBase
using Random

import HistogramFilters.generate_prior

export process_update, process_update!, measurement_update, measurement_update!,
  Filter, SparseFilter, get_data, sparsity, drop_below_threshold!

# Alias the histogram filters
const Filter = HistogramFilters.HistogramFilter
const SparseFilter = HistogramFilters.SparseHistogramFilter
const AbstractFilter = HistogramFilters.AbstractHistogramFilter

Filter(g::Grid) = Filter((1:g.width, 1:g.height))

# Provide a more general method to compute the states in a grid
get_states(g::Grid, _::AbstractFilter) = get_states(g)

get_states(::Grid, f::SparseFilter) = get_states(f)
function get_states(f::SparseFilter)
  vec = HistogramFilters.get_values(f)
  states = State[]

  # Use the existing machinery to produce the appropriate indices
  shape = CartesianIndices(size(f))
  map(x->shape[x].I, vec.nzind)
end

# Compute cartesian indices in sparse arrays
function index_to_state(f::SparseFilter, ind)
  shape = CartesianIndices(size(f))
  shape[ind].I
end

# Constructor with known initial state
function Filter(g::Grid, initial_state::State)
  data = zeros(dims(g))
  data[initial_state...] = 1

  Filter((1:g.width, 1:g.height), data)
end
function SparseFilter(g::Grid, initial_state::State; kwargs...)
  data = spzeros(prod(dims(g)))

  idx = LinearIndices(dims(g))
  data[idx[initial_state...]] = 1

  SparseFilter((1:g.width, 1:g.height), data; kwargs...)
end

# sample from the prior
function sample_state(grid::Grid, prior::Filter; rng=Random.GLOBAL_RNG)
  @views ind = sample(rng, 1:length(get_data(prior)),
                      Weights(get_data(prior)[:]))
  index_to_state(grid, ind)
end
function sample_state(grid::Grid, prior::SparseFilter; rng=Random.GLOBAL_RNG)
  vec = HistogramFilters.get_values(prior)

  # returns the linear index of the sampled state
  ind = sample(rng, vec.nzind, Weights(vec.nzval))

  # Turn this linear index into a proper state
  index_to_state(prior, ind)
end

# General process update
function process_update!(prior::Filter, transition_matrix)

  @inbounds @views mul!(get_buffer(prior)[:],
                        transition_matrix,
                        get_data(prior)[:])

  swap_buffer!(prior)

  prior
end
function process_update!(prior::SparseFilter, transition_matrix)
  # Note that the values and buffer are vectors that represent the flattened
  # state grid
  mul!(get_buffer(prior), transition_matrix, HistogramFilters.get_values(prior))

  swap_buffer!(prior)

  prior
end
function process_update(prior::AbstractFilter, transition_matrix)
  # Copy
  posterior = duplicate(prior)

  # Perform the update
  process_update!(posterior, transition_matrix)
end

# General measurement update
#
# Note: likelihoods should match the positions of the non-zeros for sparse
# matrices
function measurement_update!(prior::Filter, likelihoods::Vector{Float64})
  vals = HistogramFilters.get_values(prior)

  # update belief in place
  @inbounds @simd for ii in 1:length(vals)
    vals[ii] *= likelihoods[ii]
  end

  # normalize
  s = sum(vals)
  @inbounds @simd for ii in 1:length(vals)
    vals[ii] /= s
  end

  prior
end
function measurement_update!(prior::SparseFilter, likelihoods::Vector{Float64})
  vals = HistogramFilters.get_values(prior).nzval

  # update belief in place
  @inbounds @simd for ii in 1:length(vals)
    vals[ii] *= likelihoods[ii]
  end

  # normalize
  s = sum(vals)
  @inbounds @simd for ii in 1:length(vals)
    vals[ii] /= s
  end

  prior
end
function measurement_update(prior::AbstractFilter, xs...; kwargs...)
  # Copy
  posterior = duplicate(prior)

  # Perform the update
  measurement_update!(posterior, xs...; kwargs...)
end

# compute measurement update by computing likelihood
function measurement_update!(prior::AbstractFilter, xs...; kwargs...)
  likelihoods = compute_likelihoods(xs...; kwargs...)
  measurement_update!(prior, likelihoods)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# stitch images and estimate together

using ImageMagick
# using ImageShow
using Images, ImageView, ImageDraw

datadir = joinpath(ENV["HOME"],"data","racecar")
resultsparentdir = joinpath(datadir, "results")
imgfolder = "images"

currdirtime = "2018-11-07T01:36:52.274"
resultsdir = joinpath(resultsparentdir, "$(currdirtime)")

i = 1
for i in 1:10

tag_img = load( joinpath(resultsdir, imgfolder, "x$(i).jpg") )
top_img = load( joinpath(resultsdir, imgfolder, "img_$(i).png") )

rows = size(im1)[1]+size(im1t)[1]

end

#


text:
# Функционал, его градиент и минимизация

@doc raw"""
    J_q(uˢ::Matrix, ψˢ::Matrix,
        Xₙ::Array, N::Int,
        Tₘ::Vector, M::Int) -> Vector

Вычисляет градиент ``J_q(x) = \int_0^T u^s(x, t) \psi^s(x, t) dt ``.
"""
function J_q(uˢ::Matrix, ψˢ::Matrix,
             Xₙ::Matrix, N::Int,
             Tₘ::Vector, M::Int,
             X_q::Vector = [NaN, NaN])

    size(Xₙ) == (N+1, M+1) ||
    throw(ArgumentError("size(Xₙ) == $(size(Xₙ)), N == $(N), M == $(M) " *
                        "Массив Xₙ должен иметь размерность (N+1, M+1)."))

    isDynamicMesh = ! isapprox(Xₙ[:, 1], Xₙ[:, end])
    # Когда мы понимаем, что у нас динамическая сетка,
    # мы должны убедиться, что нам передали необходимую вспомогательную сетку.
    if isDynamicMesh
        if all(isnan.(X_q))
            throw(ArgumentError("Xₙ передан в виде матрицы, значит подразумевается," *
                                "\nчто используется динамическая сетка. Это требует" *
                                "\nпередачи вспомогательной сетки `X_q`, для ``q``"))
        end
    end
    @assert length(Tₘ) == M+1

    # заполним рабочие массивы специальным образом
    # Если сетка динамическая,
    # Нам нужна интерполяция решений на вспомогательную сетку `X_q`
    if isDynamicMesh
        J_q = zero(X_q);                # Градиент будет определен на вспомогательной сетке
        u = zeros(length(X_q), M+1);
        ψ = zeros(length(X_q), M+1);
        for m in 1:M+1
            X = Xₙ[:, m];
            uspl = Spline1D(X, uˢ[:, m]);
            ψspl = Spline1D(X, ψˢ[:, m]);
            u[:, m] = uspl(X_q);
            ψ[:, m] = ψspl(X_q);
        end
    else                        # Если все вычислялось на статических сетках
        u = uˢ;                 # То будем работать с исходными массивами
        ψ = ψˢ;
        J_q = zeros(N+1);       # Градиент определяется на сетке Xₙ
    end

    #  Здесь уже не важно, на какой сетке определены `u`, `ψ`
    #  Вычислим интеграл по формуле трапеций.
    #  ``\int_0^T u^s(x, t) \psi^s(x, t) dt ``
    for m in 1:M
            τ = Tₘ[m+1] - Tₘ[m]
            J_q  += (( u[:, m] .* ψ[:, m] ) + ( u[:, m+1] .* ψ[:, m+1] )).* τ / 2
    end

    return -J_q
end

@doc raw"""
    J(uˢ::Matrix, Xₙ::Array, N::Int,
      Tₘ::Vector, M::Int,
      f1::Vector, f2::Vector,
      qₙˢ::Vector, α::Real = 0.0) -> Real

Вычисляет функционал `` J(\mathbf{x}) = \int_0^T \left( u(f_1(t), t; q^s) - f_2(t) \right)^2 +
\alpha \int_0^1 q^2(x) dx ``.

Использует [`f2`](@ref) для вычисления ``u(f_1(t), t)`` и после по формуле трапеций.
"""
function J(uˢ::Matrix, Xₙ::Array, N::Int,
           Tₘ::Vector, M::Int,
           f1::Vector, f2::Vector,
           qₙˢ::Vector, α::Real = 0.0)
    J = 0.0
    ξ = zeros(M+1);

    ξ = NonLinearReactionAdvectionDiffusionWithFrontData.f2(f1, uˢ, Xₙ, N, M);

    for m in 1:M
            τ = Tₘ[m+1] - Tₘ[m];
            J += ( (ξ[m] - f2[m])^2 + (ξ[m+1] - f2[m+1])^2) * τ / 2
    end

    return J
end

@doc raw"""
    minimize(q₀::Vector, u₀::Vector,
             ulₘ::Vector, urₘ::Vector,
             Xₙ, N,
             Tₘ::Vector, M,
             ε,
             f1_data::Vector, f2_data::Vector;
             S::Int = 10,
             β::Real = 0.01,
             w::Real = 0.0001,
             create_mesh::Function = x -> [NaN, NaN])  -> Vector, Vector, Matrix

Вернет `qˢ` на оригинальной сетке `Xₙ[:,1]`.
"""
function minimize(q₀::Vector, u₀::Vector,
                  ulₘ::Vector, urₘ::Vector,
                  Xₙ, N,
                  Tₘ::Vector, M,
                  ε,
                  f1_data::Vector, f2_data::Vector;
                  S::Int = 10,
                  β::Real = 0.01,
                  w::Real = 0.0001,
                  create_mesh::Function = x -> [NaN, NaN],
                  showProgress = false)

    ψ₀ = zero(Xₙ);
    ψl = zero(Tₘ);
    ψr = zero(Tₘ);

    isDynamicMesh = all(isnan.(create_mesh(Xₙ[end - div(end,2)])))
    @info isDynamicMesh ? "Используем оригинальную сетку Xₙ для q₀" : "Используем отдельную сетку для q₀"

    # Создадим интерполяционный объект для q
    # Он остается постоянным, вдоль всего цикла минимизации
    # Сетка Xₙ — статическая, определенная только для начального
    # приближения, то, что она двигается — дело сугубо внутрннее
    # для прямой и сопряженной задачи.
    qspl = Spline1D(Xₙ, q₀)
    k = 100;                                 # Кол-во интервалов в вспомогательной сетке
    if isDynamicMesh
        X_q = Xₙ;
    else
        X_q = [ first(Xₙ) + n * (last(Xₙ) - first(Xₙ))/Float64(k) for n in 0:k]
    end
    # вектор сеточных значений q на вспомогательной сетке.
    q_aux = qspl(X_q);

    # Проверим, что нам передали q₀ именно на сетке Xₙ
    if ! isapprox(qspl.(Xₙ), q₀)
        throw(ArgumentError("Сеточные значения q₀ определены на какой-то другой сетке, а не Xₙ"))
    end

    #' ## Подготовка к итерационному процессу
    qˢ = q₀;                    # Приближение вектора `q` на текущей итерации.
    uˢ = zeros(N+1, M+1);       # Решение прямой задачи на текущем шаге итерации.
    ψˢ = zeros(N+1, M+1);       # Решение сопряженной задачи на текущем шаге итерации.
    J_values = zeros(S);        # История значений градиента.
    Q_values = zeros(N+1, S);   # История вектора `q`.

    p = showProgress ? Progress(S, 5, "Iterating minimization loop S=$(S)... ") : nothing;

    for s in 1:S

        showProgress && next!(p);

        Q_values[:, s] .= qˢ; # Сохраним вектор q на текущей итерации

        # Нахождения решения прямой задачи внутри итерационного числа
        # в аргументах отличается лишь вектором `qˢ`.
        # Решение `u` зависит от него как от параметра.
        uˢ, XXˢ, TPˢ = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qˢ, create_mesh = create_mesh);

        # Решаем сопряженную задачу, на вход которой подадим
        # `uˢ`      — Параметр, изменяющийся внутри цикла
        # `f1, f2`  — Параметры, сформированные единожды до цикла

        ψˢ = solve_adjoint(ψ₀, XXˢ, N, Tₘ, M, ε, ψl, ψr, qˢ, uˢ, f1_data, f2_data, w=w);

        J_values[s] = J(uˢ, XXˢ, N, Tₘ, M, f1_data, f2_data, qˢ);

        ∇J = J_q(uˢ, ψˢ[:,end:-1:1], XXˢ, N, Tₘ, M, X_q);

        # Делаем шаг минимизации
        q_aux = q_aux - β * ∇J

        qs_spl = Spline1D(X_q, q_aux)
        qˢ = qs_spl(Xₙ);

    end

    return (qˢ, J_values, Q_values)
end


text:
"""
    @test_cases [ntests = 100] [maxtests = 1000] [logto = ""] [working_mode = :none] [argument_data, proposition]

    Takes proposition depending on some variables, generates samples of their values
    and tests the proposition. Prints number of samples for which the proposition evaluated to `true`.

## Arguments

Macro accepts a tuple of expressions that describe either statement of a test
or specify optional keyword arguments.

Statement of a test must be a sequence of comma-separated expressions, where
all but last expressions are read as `argument_data` and last statement is
`proposition` depending on these arguments.

* `argument_data` consists of arbitrarily many expressions containing
information about variables, their types and conditions on them. Expressions
must have one of the following forms:

  * `a<x::T<b` or `a>x::T>b`, where middle variable `x` must be a symbol;
  `T`  must be a type; bounds `a` & `b` can be arbitrary expressions;
  * `a<<x::T<b`, `a<x::T<<b`, `a>>x::T>b` or `a>x::T>>b`,
  similar to the one above;
  * `x::T`, where `T` is a julia type for which function `custom_generator(x::T...)`
  is specified (see ~/examples/new-type-generator.jl for more details).
  * For inequalities shorthand `x` instead of `x::T` may be used, e.g. `a<x<b`,
  then type of variable defaults to `Float64`.

* `proposition`: is the expression to be tested. Two types of propositions are allowed:

  * `proposition` of form `condition-->statement`, where `condition` and
  `statement` are boolean-valued expressions
  * `proposition` is boolean-valued expression itself

Optional keyword arguments include the following:

* `ntests = 100`: number of tests that must be performed.

* `maxtests = 1000`: maximal number of attempts to generate arguments
satisfying `condition`

* `logto = "path_to_file"`: if provided then parsed as a path to file,
where log of tests will be recorded.

* `working_mode = :none`: accepts one of three symbols: `:none`, :test_formany`,
  `:test_exists`. For modes different from `:none` see documentation for `@test_formany`
  and `@test_exists`.

## Description

Macro `@test_cases` provides a convenient syntax for testing a boolean-valued
proposition multiple times against randomly generated input. Output prints number
of cases for which the proposition was satisfied.

The macro attempts to generate a sample of values for all variables identified in
`argument_data`. Value of every variable is generated according to its type `T`
and rules specified for this type in function `custom_generator(x::T...)`.
Generated input is used to evaluate the statement of the test.

If `proposition` is just a boolean-valued expression, code simply checks whether
`prop` holds for generated values.

If `proposition` is of form `condition-->statement` (see arguments list),
then test checks whether implication holds. Namely, the test passes if whenever
`condition` evaluates to `true`, then `statement` also evaluates to `true`.

In the latter setting the test breaks and returns error, if the number of attempts to
generate argument values satisfying `condition` exceeds `maxtests` parameter before
succesfully performing `ntests` checks of `proposition`. Thus, some caution is
required to avoid vacuous conditions. Failure to generate `ntests` values
satisfying the `condition` in `maxtests` attempts results in error.

Default for optional keyword arguments is the following: `ntests = 100, maxtests = 1000`.
If only `ntests` parameter is specified, `maxtests = 10*ntests`. If `logto = "path_to_file"`
is given, log of test will be written to corresponding file,for example:

* `@test_cases ntests = 10 maxtests = 1000  1<x<10, proposition logto =  "./test_log.csv"`

Working mode argument has default `working_mode = :none` and is designed to
refer to @test_cases` code from other macros `@test_formany` and `@test_exists`.
Macro calls `@test_cases working_mode = :test_formany argument_data, proposition`
and `@test_formany argument_data, proposition` are fully equivalent.

## Examples

```julia
# example with default genenerator that draws values from uniform distribution
@test_cases ntests = 10000 0<x<100, x>=50
Test passed in 5035/10000 simulations.

# simulated values can be recorded, using `logto` optional field
## theoretical expectation .5*log(2) +.5 ~ 0.847...
@test_cases ntests = 10000 1<x<1, 0<y<1, x*y <.5 logto =  "./test_log.csv"
Test passed in 8399/10000 simulations.

#elementary example of testing cdf properties using Checkers
type gaussian_rv
	val::Float64
end

function Checkers.custom_generator{T<:gaussian_rv}(var_type::Type{T}, mu::Number, sigma::Number)
	return (size) -> gaussian_rv(mu + sigma*randn())
end

@test_cases x::gaussian_rv, 0, 1, x.val<1.96 ntests = 10000
Test passed in 9777/10000 simulations.
```
"""
macro test_cases(exprs...)
    # everything is embedded in try/catch block and if it can't be compiled,
    # cathch-part returns Error type object from Base.test
    output_expr = try

    #  identify pseduo-keyword arguments
    maxtests = 0; ntests = 0; logto = ""
    # working mode
    working_mode = :none # can be :test_formany or :test_exists

    prop = :(); cond = :()

    # parsing expression
	var_data = Array(Any, 3, 1)
    for ex in exprs
        if isa(ex, Expr) && ex.head == :tuple  # parsing statements about
            var_data = parse_argument_data(ex) # variables, condition and
            if ex.args[end].head == :-->       # proposition
                prop = ex.args[end].args[2]
                cond = ex.args[end].args[1]
            else
                prop = ex.args[end]
                cond = true
            end
        elseif isa(ex, Expr) && ex.head == :(=) # parsing pseduo-keyword arguments
            if ex.args[1] == :ntests
                ntests = ex.args[2]
            elseif ex.args[1] == :maxtests
                maxtests = ex.args[2]
            elseif ex.args[1] == :logto
                logto = ex.args[2]
            elseif ex.args[1] == :mode
                if ex.args[2] in (:test_formany, :test_exists, :none)
                    working_mode = ex.args[2]
                else
                    error("Unrecognized mode $(ex.args[2]).")
                end
            else
                error("Invalid macro input $ex.")
            end
        else
            error("Invalid macro input $ex.")
        end
    end

    # setting default parameters for ntests and maxtests
    if ntests == 0
        ntests = 100
    end
    if maxtests == 0
        maxtests = 10 * ntests
    end

    # creating a block expression `generate_values` and expression values
    # that stores references to generated variable
    num_of_vars = size(var_data,2)
    generate_values = Expr(:block)
    values = Expr(:vect)
    for i in 1:num_of_vars
        params = Expr(:call, :custom_generator, esc(var_data[2,i]))
        for p in var_data[3,i]
            push!(params.args, esc(p))
        end
        next_expr = :($(esc(var_data[1,i]))= $params(div(n,2)+3))
        push!(generate_values.args, next_expr)
        push!(values.args, :($(esc(var_data[1,i]))))
    end
    # defining innermost expression in the loop running tests
    inner_ex = quote
        num_good_args += 1
        tmp = $(esc(prop))::Bool
        num_passes += tmp

        # possibly log here

        # possibly record break_values and break if mode = test_exists/test_formany

        if num_good_args >= $ntests
            break
        end
    end

    if logto != ""
        insert!(inner_ex.args, length(inner_ex.args) - 1,
         quote
             writedlm(log_file, transpose(push!(convert(Vector{Any},$values),tmp)), ",")
         end
        )
    end

    if working_mode == :test_formany
        insertion_ex = quote
            if !tmp
                break_values = $values
                break_test = true
                break
            end
        end

        insert!(inner_ex.args, length(inner_ex.args) - 1, insertion_ex)
    end

    if working_mode == :test_exists
        insertion_ex = quote
            if tmp
                break_values = $values
                break_test = true
                break
            end
        end

        insert!(inner_ex.args, length(inner_ex.args) - 1, insertion_ex)
    end

    # defining output expression
    output_expr = quote
        num_good_args = 0
        num_passes = 0
        break_values = []
        break_test = false

        # possibly open file to log values

        for n in 1:$maxtests
            $generate_values

            if $(esc(cond))
                $inner_ex
            end
        end

        if !break_test && num_good_args < $ntests
            nt = $ntests
            error("Found only $num_good_args/$nt values satisfying given condition.")
        end

        # if logging close file here

        # if mode == none print the result
        # else specify result in Pass/Fail format
    end

    # depending on pseudo-args modify the output_expr
    if working_mode == :none
        show_result_expr = quote
            nt = $ntests
            print_with_color(:green, "Test passed in $num_passes/$nt simulations.\n")
        end
        push!(output_expr.args, show_result_expr)
    end

    if logto!= ""
        logging_expr = quote
            log_file = open($logto,"a")
            writedlm(log_file,
                reshape([string(s) for s in $(var_data[1,:])],(1,$num_of_vars)),",")
        end

        unshift!(output_expr.args, logging_expr)
        push!(output_expr.args, quote close(log_file) end)
    end

    if working_mode == :test_formany
        insertion_ex = quote
            result = !break_test ?
                 Pass(:test,$(Expr(:quote, exprs)), nothing, nothing) :
                 Fail(:test,
                    $(Expr(:quote, exprs)),
                    [string(s[1])*" = "*string(s[2]) for s in zip($(var_data[1,:]), break_values)],
                    nothing)

            if isa(result, Fail)
                println(result.data)
            end

            record(get_testset(), result)
        end

        push!(output_expr.args, insertion_ex)
    end

    if working_mode == :test_exists
        insertion_ex = quote
            result = break_test ?
                Pass(:test,
                    $(Expr(:quote, exprs)),
                    [string(s[1])*" = "*string(s[2]) for s in zip($(var_data[1,:]), break_values)],
                    nothing) :
                Fail(:test,$(Expr(:quote, exprs)), nothing, nothing)

            if isa(result, Pass)
                println(result.data)
            end

            record(get_testset(), result)
        end

        push!(output_expr.args, insertion_ex)
    end

    output_expr = Expr(:try, output_expr, :err,
        quote
            result = Error(:test_error, $(Expr(:quote, exprs)), err, catch_backtrace())
            record(get_testset(), result)
        end
        )

    output_expr

    # if something of the above didn't work, will return error
    catch err
        output_expr = quote
        	result = Error(:test_error, $(Expr(:quote, exprs)), $err, catch_backtrace())
            record(get_testset(), result)
        end
    end

    return output_expr
end


text:
let

    @info("Testing FileMTimeEvent")
    e = FileMTimeEvent()
    file = joinpath(tempdir(), "FileMTimeEvent.test")
    touch(file)
    update!(e ,file)

    events_count = 0
    trigger_at = [3, 4, 5]
    for it in 1:10
        if has_event!(e ,file)
            println(basename(file), " changed!!!")
            events_count += 1
        end

        if it in trigger_at
            touch(file)
        end
        sleep(0.1)
    end
    @test events_count == length(trigger_at)

    return nothing
end


text:
function linear_kernel(x1::Vector, x2::Vector)
  return dot(x1, x2)::Float64
end

function polynomial_kernel(x1::Vector, x2::Vector, d::Integer=2)
  return ((1 + dot(x1, x2)) ^ d)::Float64
end

# function gaussian_kernel(x1::Vector, x2::Vector, sigma::Real=1.0)
#   # k = exp(-1/(2*sigma^2)  ||x1 - x2||^2 )
#   return exp( -norm(x1 - x2)^2 / (2*sigma^2) )::Float64
# end

function gaussian_kernel(x1::Vector, x2::Vector, gamma::Real=1.0)
  # k = exp(-gamma  ||x1 - x2||^2 )
  return exp( -gamma * norm(x1 - x2)^2 )::Float64
end

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
# Utils for every script to load kayak experimental data

function getPaths(expID::String , expType::String ; currloc::String="local" , trialID::Int=-1)
    if currloc == "local"
        topdir = joinpath(ENV["HOME"],"data", "kayaks")
    elseif currloc == "lj"
        topdir = joinpath("/media","data1","data","kayaks")
    end

    if expType == "range"
        expstr = "rangeOnly";
    elseif expType  == "dynsas"
        expstr = expType
    elseif expType == "pointsas"
        expstr = expType
    elseif expType == "dynsource"
        expstr = expType
    end

    if expID == "dock"
        trialstr = "20_gps_pos";
        gtjldpath = ""
    elseif expID == "drift"
        trialstr = "08_10_parsed";
        if trialID == -1
            error("Specify Trial ID: 1-3")
        else
            gtjldpath = joinpath(topdir,trialstr,"exp$(trialID).jld")
        end
    end

    rangedir = joinpath(topdir, expstr*"_"*trialstr)
    datadir = joinpath(topdir,trialstr)
    savedir = joinpath(topdir,expstr)

    return datadir, rangedir, savedir, gtjldpath
end

function loadNav(wIn::Array, allpaths)
    thisIter = 1:wIn[end]-wIn[1]
    navalltemp = zeros(length(thisIter),2);
    for i in thisIter
        navfile = joinpath(allpaths[1],"nav$(wIn[i]).csv")
        navalltemp[i,:] = readdlm(navfile,',',Float64,'\n')
    end
    return navalltemp
end

function loadRanges(wIn::Array, allpaths)
    thisIter = 1:wIn[end]-wIn[1]
    rangealltemp = zeros(7501,2,length(thisIter));
    for i in thisIter
        rfile = joinpath(allpaths[2],"range$(wIn[i]).txt")
        rangealltemp[:,:,i] = readdlm(rfile,',',Float64,'\n')
    end
    return rangealltemp
end

function loadGT(wIn::Array, expID, allpaths)
    if expID == "dock"
        igttemp = [17.0499;1.7832];
    elseif expID == "drift"
        tmpdict = load(allpaths[4])
        igtFull = tmpdict["icarus_gt"]
        windowstart = tmpdict["ibegin"];
        windowend = tmpdict["iend"];

        if length(wIn) > (windowend - windowstart)
            error("Specified trial window is too long")
        else
            igttemp = igtFull[wIn[1]-windowstart:length(wIn),:]
        end
    end
    return igttemp
end


text:
#push!(LOAD_PATH,utf8("T:\\Research\\sdelap") )
using SDELab2 # SDELab2.jl should be findable via LOAD_PATH
println("run_FHN.jl")

#
# FitzHugh-Nagumo example from
#
# H. Alzubaidi and T. Shardlow. Numerical simulations of SDEs and SPDEs.
# In: Stochastic Methods in Neuroscience. Ed. by C. Laing and G. Lord. OUP,
# 2009. Chap. 12, pp. 344–366.
# doi: 10.1093/acprof:oso/9780199235070.003.0012 .
#
type fhn_params
  D
  Mu
  a
  b
  gamma
end

function drift_f(y,p::fhn_params)
  d=round(Int,length(y)/2)
  A=SymTridiagonal(-2*ones(d),ones(d-1))
  u=y[1:d] # fasts
  v=y[(d+1):end]
  F=(u.*(1-u).*(u-p.a))-v
  B=p.D*((d+1)^2)*A*u
  z1=B+F+p.Mu
  z2=p.b*(u-p.gamma*v)
  return[z1;z2]
end
#
#
function diff_g_mat(u,dw,sigma)
  d=length(dw)
  Q=sqrt(d+1) # orthognoal matrix does not change distribution
  return   [sigma*Q*dw;zeros(d)]
end
#
function diff_g_vecs(sig,m)
  A=Array(Function,m)
  for i=1:m
    A[i]=(u->1)
  end
  # do not use Milstein methods (no point for additive!)
  return A
end
# dimensions
d=9
d_all=2*d
noise_type=2 # diagonal
# params
p=fhn_params(0.01,0.5,0.05,0.008,0.5)
sig=0.05
# initial data
a0_n = -.1/(exp(1.)-1.); b0_n=0.125
v0=a0_n/(a0_n+b0_n)
y0=[zeros(d);v0*ones(d)]
#
t0=0
tf=200
#
fcn=set_fcn(d_all,d, x->drift_f(x,p),(x,dw)->diff_g_mat(x,dw,sig),diff_g_vecs(sig,d), noise_type)
opt=set_opt(0.01,"SIBDF")
opt["do_plotting"]=true
opt["output_plot_type"]="path_plot"
opt["output_sel"]=1:d
using PyPlot
println("Start sde_strong_soln")
tic();t,y,W=sde_strong_solution(fcn, linspace(t0,tf,100), y0, opt);run_time=toc()

#


text:
using Documenter, BEAST

makedocs(
    clean=false,
    sitename="BEAST Documentation")

deploydocs(
    # deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/krcools/BEAST.jl.git",
    # julia  = "1.0",
)


text:
module PreProcessing

end