text:
using DataStructures

function parseSnapData(fname::String)
    PP, T, PP_test, T_test, n_test = trainTestSplitSnapData(fname, 1.0)
    return PP, T
end

function trainTestSplitSnapData(fname::String, split_prop::Float64=0.8)
    N = countlines(open(fname))
    split_n = Int(round(N*split_prop))
    n_test = N - split_n

    f = open(fname)
    arrival_times = OrderedDict{String, Int}()
    degrees = OrderedDict{String, Int}()
    #scope
    T = nothing
    PP = nothing
    for (i, ln) in enumerate(eachline(f))
        a = split(ln)
        start = a[1]
        terminal = a[2]
        if ~haskey(arrival_times, start)
            arrival_times[start] = 2*i - 1
            degrees[start] = 1
        else
            degrees[start] += 1
        end
        if ~haskey(arrival_times, terminal)
            arrival_times[terminal] = 2*i
            degrees[terminal] = 1
        else
            degrees[terminal] += 1
        end
        if i == split_n
            T = collect(values(arrival_times))
            PP = collect(values(degrees))
        end
    end
    # T_test is a pure extension of T
    T_test = collect(values(arrival_times))
    # PP is a pure extension of PP
    PP_test = collect(values(degrees))
    return PP, T, PP_test, T_test, 2*split_n, 2*n_test
end

function generateInterarrivalTimes(TK::Char, N::Int, interarrival_dist::DiscreteDistribution)
  ia_dist = (x,y) -> interarrival_dist
  generateInterarrivalTimes(TK,N,ia_dist)
end

function generateInterarrivalTimes(TK::Char, N::Int, interarrival_dist::Function)
    """
    - `TK`: 'T' if `N` is the total number of observations;
        'K' if `N` is the total number of arrivals
    - `N`: number of arrival times to generate (modulated by `TK`)
    - `interarrival`: distribution object to generate i.i.d. interarrivals
    """

    # check function arguments
    TK != 'K' && TK != 'T' ? error("`TK` must be 'T' or 'K'") : nothing

    zero_shift = Int(minimum(interarrival_dist(1,1)) == 0)

    if TK == 'K'
      T = zeros(Int64,N)
      T[1] = 1 # first arrival time is always 1
      for j in 2:N
        T[j] = rand(interarrival_dist(T[j-1],j-1)) + T[j-1] + zero_shift
      end
      # return T
    else
      T = [1]
      j = 1
      while T[j] < N
        j += 1
        append!(T,rand(interarrival_dist(T[j-1],j-1)) + T[j-1] + zero_shift)
      end
      if T[end] > N
        pop!(T)
      end
      # return T
    end
    return T
end

function generateLabelSequence(N::Int, alpha::Float64,
        interarrival_dist::DiscreteDistribution)
    ia_dist = (x,y) -> interarrival_dist
    generateLabelSequence(N,alpha,ia_dist)
end

function generateLabelSequence(N::Int, alpha::Float64,
        interarrival_dist::Function)
    """
    - `N`: number of observations in the sequence
    - `alpha`: BNTL α
    - `interarrival_dist`: distribution object to generate interarrivals
    """
    Z = zeros(Int, N) # sequence of labels
    T = generateInterarrivalTimes('T', N, interarrival_dist)
    K = size(T,1) # number of clusters
    PP = zeros(Int, K) # arrival-ordered partition counts

    k = 0
    for n in 1:N
      if n <= T[end] && n == T[k+1]
        k += 1
        PP[k] = 1
        Z[n] = k
        k > K ? k = K : nothing
      else
        Z[n] = wsample(1:k, PP[1:k] .- alpha) # discounted size-biased sample
        PP[Z[n]] += 1
      end
    end
    return Z, PP, T
end

function generatePsis(T::Vector{Int},alpha::Float64)
    """
    - `T`: Arrival times
    - `alpha`: 'discount' parameter
    """
    K = size(T,1)
    Psi = zeros(Float64,K)
    Psi[1] = 1
    for j in 2:K
      Psi[j] = rand(Beta(1 - alpha, T[j] - 1 - (j-1)*alpha))
    end
    return Psi
end


text:
function move_to_filled(invest, result, client, invest_path)
    trades = client.get_my_trades(symbol=invest["order"]["symbol"], orderId=invest["order"]["orderId"])

    invest["status"] = "Done"
    invest["trades"] = trades
    invest["amount"] = result["executedQty"]
    invest["value"]  = result["cummulativeQuoteQty"]
    invest["order"]  = result
    datum = split(string(today()), "-")
    new_filename = joinpath(dirname(invest_path), "filled/" * datum[1] * "/" *
        datum[2] * "/" * invest["symbol"] * "_" * string(now()) * ".json")
    !isdir(dirname(new_filename)) ? mkpath(dirname(new_filename)) : ""
    rm(invest_path)
    open(new_filename, "w") do f
        JSON.print(f, invest, 4)
    end
end
export move_to_filled

ext: jl


text:
using Sockets: IPv4
is_ch_type(::Val{:IPv4})  = true
result_type(::Val{:IPv4})  = Vector{IPv4}

function read_col_data(sock::ClickHouseSock, num_rows::VarUInt, ::Val{:IPv4})
    tmp = Vector{IPv4}(undef, UInt64(num_rows))
    Base.read!(sock.io, tmp)
    return tmp
end

function write_col_data(sock::ClickHouseSock,
                        data::AbstractVector{IPv4}, ::Val{:IPv4})

    Base.write(sock.io, data)
end


text:
using Documenter
using EventTracker
using EventTrackerBase

makedocs(
    sitename = "EventTracker",
    format = Documenter.HTML(),
    modules = [EventTracker, EventTrackerBase],
    checkdocs = :export,  # `:all` contains many false-positives
    strict = true,
)

deploydocs(
    repo = "github.com/tkf/EventTracker.jl",
    push_preview = true,
)


text:
using LinearAlgebra

mutable struct Network

    n_exc::Int64
    n_neurons::Int64

    v::Array{Float64}
    u::Array{Float64}
    d::Array{Float64}
    a::Array{Float64}

    connection::Array{Int64,2} # WARNING read by columns !!!!
    s::Array{Float64,2} # WARNING read by columns !!!!
    sd::Array{Float64,2} # WARNING read by columns !!!!
    STDP::Array{Float64}

    exc::Array{Int64}
    inh::Array{Int64}

    I::Array{Float64}
    da::Float64
    arch::Array{Int64}
    #param::Dict

    function connect(n_neurons::Int64, n_exc::Int64, connectivity::Float64,
                     arch::Array{Int64})
        if arch == [0]
            connection = (rand(n_neurons,n_neurons) .< connectivity)
        else
            architr = [0; arch]
            connection = zeros(n_neurons,n_neurons)
            for i in 1:length(architr)-2
                si1 = sum(architr[1:i])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                ai2 = architr[i+1]
                ai3 = architr[i+2]
                connection[1+si1:si2,1+si2:si3] .= ones(ai2,ai3)
            end
            archinh = [0;arch[2:end-1]]
            for i in 1:length(archinh)-1
                sii1 = sum(archinh[1:i])
                sii2 = sum(archinh[1:i+1])
                si2 = sum(architr[1:i+1])
                si3 = sum(architr[1:i+2])
                a2 = arch[i+1]
                connection[n_exc+sii1+1:n_exc+sii2, si2+1:si3] .= ones(a2,a2) .- I(a2)
                connection[si2+1:si3, n_exc+sii1+1:n_exc+sii2] .= I(a2)
            end
        end
        exc = collect(1:n_exc)
        inh = collect(n_exc+1:n_neurons)
        return connection, exc, inh
    end

    function Network(arch::Array{Int64})
        n_exc = sum(arch)
        n_neurons = n_exc+sum(arch[2:end-1])
        v = -65.0*ones(n_neurons) #parameter
        u = 0.2 * v # parameter
        d = [(i<=n_exc) ? 8.0 : 2.0 for i in 1:n_neurons] # parameter
        a = [(i<=n_exc) ? 0.02 : 0.1 for i in 1:n_neurons] # parameter
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, 0.1, arch)
        connection = transpose(connection)
        s = zeros(n_neurons,n_neurons)
        # maybe set the Exc-> Inh conn weights to high value
        s[:,exc] .= 1.0 .* ones(n_neurons, n_exc) # parameter
        s[:,inh] .= -1.0 .* ones(n_neurons, n_inh) # paramter
        s .*= connection
        sd = zeros(n_neurons,n_neurons)
        STDP = zeros(n_neurons)
        I = Float64[]
        da = 0.0
        arch = arch
        new(n_exc,n_neurons,v,u,d,a,connection,s,sd,STDP,exc,inh,I,da,arch)
    end

    function Network(n_neurons::Int64, connectivity::Float64, n_exc::Int64)
        v = -65.0*ones(n_neurons)
        u = 0.2 * v
        d = [(i<=n_exc) ? 8.0 : 2.0 for i in 1:n_neurons]
        a = [(i<=n_exc) ? 0.02 : 0.1 for i in 1:n_neurons]
        n_inh = n_neurons-n_exc
        connection, exc, inh = connect(n_neurons, n_exc, connectivity, [0])
        connection = transpose(connection)
        s = zeros(n_neurons,n_neurons)
        # maybe set the Exc-> Inh conn weights to high value
        s[:,exc] .= 1.0 .* ones(n_neurons, n_exc) # parameter
        s[:,inh] .= -1.0 .* ones(n_neurons, n_inh) # paramter
        s .*= connection
        sd = zeros(n_neurons,n_neurons)
        STDP = zeros(n_neurons)
        I = Float64[]
        da = 0.0
        arch = [0]
        new(n_exc,n_neurons,v,u,d,a,connection,s,sd,STDP,exc,inh,I,da,arch)
    end
end

# %% functions

function input!(net::Network, input_spikes::Array{Int64})
    net.v[input_spikes] .= 40.0 # parameter
end

function spike!(net::Network)
    # spikes
    net.I = 13 .* (rand(length(net.STDP)) .- 0.5) # parameter
    spiked = findall(x->x>30,net.v) # parameter
    net.v[spiked] .= -65.0 # parameter
    net.u[spiked] .+= net.d[spiked]
    # increment current
    for neuron in spiked
        net.I .+= net.s[:,neuron]
    end
    #time step (Izhikevic model)
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I) # parametersss
    net.v .= @. net.v+0.5*((0.04*net.v+5)*net.v+140-net.u+net.I)
    net.u .= @. net.u+net.a*(0.2*net.v-net.u)   # parameters
    net.STDP[spiked] .= 0.1  # parameter
    # LTP & LTD
    for neuron in spiked
        net.sd[neuron,:] .+= 1.0 .* net.STDP #parameter
        net.sd[:,neuron] .-= 1.5 .* net.STDP # parameter
    end
    # time decay
    net.STDP .*= 0.95 # parameter
    net.da *= 0.995 # parameter
    spiked
end

function learn!(net::Network)
    #increase the synaptic weigths
    net.s[:,net.exc] .+= (0.002+net.da) .* net.sd[:,net.exc]  # parameters
    net.s[:,net.exc] .= clamp.(net.s[:,net.exc], 0, 4)
    net.s .*= net.connection
    # time decay
    net.sd .*= 0.99 # parameter
end

# %% main functions

function step!(net::Network, train::Bool=false)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

function step!(net::Network, input_spikes::Array{Int64}, train::Bool=false)
    input!(net,input_spikes)
    spiked = spike!(net)
    if train
        learn!(net)
    end
    spiked
end

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
@testset "environment.jl" begin

    @testset "CPEnv{TS}" begin
        graph = Matrix(adjacency_matrix(random_regular_graph(6, 3)))
        ts = SeaPearl.DefaultTrajectoryState(GeometricFlux.FeaturedGraph(graph; nf=rand(3, 6)), 1)
        env = SeaPearl.CPEnv{SeaPearl.DefaultTrajectoryState}(
            .0,
            false,
            ts,
            collect(1:4),
            [1, 4],
            [true, false, false, true]
        )

        Random.seed!(0)
        actions = Set([agent(env) for I = 1:10])

        @test actions == Set([1,4])
    end
end


text:
"""
    annotations(strings::Vector{String}, positions::Vector{Point})

Plots an array of texts at each position in `positions`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Annotations, text, position) do scene
    default_theme(scene, Text)
end

function convert_arguments(::Type{<: Annotations},
                           strings::AbstractVector{<: AbstractString},
                           text_positions::AbstractVector{<: Point{N}}) where N
    return (map(strings, text_positions) do str, pos
        (String(str), Point{N, Float32}(pos))
    end,)
end

function plot!(plot::Annotations)
    # annotations are not necessary anymore with the different text behavior
    text!(plot, plot[1]; plot.attributes...)
    plot
end


text:
struct Generator
    𝐕ₜ::Set{Terminal}
    𝐕ₙ::Set{Nonterminal}
    𝐏::Production
    𝐀::Dict{Terminal,Function}
end

Generator(𝐆::Grammar) = Generator(𝐆.𝐕ₜ, 𝐆.𝐕ₙ, 𝐆.𝐏, Dict{Terminal,Function}())

function bind(g::Generator, sym::Terminal, f::Function)
    sym ∉ g.𝐕ₜ && throw(KeyError("$sym not in 𝐕ₜ"))
    g.𝐀[sym] = f
end

function bind(g::Generator, syms::Vector{Terminal}, fs::Vector{Function})
    foreach((s, a) -> bind(g, s, a), zip(syms, fs))
end

"广度优先搜索生成可能的句型"
function generate_bfs(g::Generator, token::Nonterminal;max_depth::Int64=2)
    queue = Statement[]
    used = Set{Statement}()
    result = Set{Statement}()
	isterminal = ∈(g.𝐕ₜ)
    for p in g.𝐏[token]
		if all(isterminal, p)
			push!(result, p)
		else
			push!(queue, p)
			push!(used, p)
		end
	end
    step = 0
    while !isempty(queue)
        num = length(queue)
        @info num
        for _ in 1:num
            p = popfirst!(queue)
            idx = findall(∈(g.𝐕ₙ), p)
            reverse!(idx)
            tmp = Statement[p,]
            for i in idx
                n = length(tmp)
                for j in 1:n
                    now = popfirst!(tmp)
                    for pp in g.𝐏[p[i]]
                        np=splice(now, i, pp)
                        if all(isterminal, np)
                            push!(result, np)
                        else
                            if np ∉ used && step < max_depth
                                push!(tmp, np)
                                push!(used, np)
                            end
                        end
                    end
                end
            end
            for np in tmp
                push!(queue, np)
            end
        end
        step = step + 1
	end
	true, result |> collect
end

"深度优先搜索产生可能的句型"
function generate_dfs(g::Generator, token::Nonterminal;max_depth::Int64=10)
    max_depth == 0 && return (false,nothing)
    ps = g.𝐏[token] |> collect
    for i in length(ps):-1:1
        remove_pi = false
        idx = findall(∈(g.𝐕ₙ), ps[i])
        isempty(idx) && continue
        nps = Statement[ps[i],]
        reverse!(idx)
        for j in idx
            (flag, res) = generate_dfs(g, ps[i][j], max_depth=max_depth - 1)
            if !flag
                remove_pi = true
                break
            end
            tmp = Statement[]
            while !isempty(nps)
                now = pop!(nps)
                for r in res
                   np=splice(now, j, r)
                   push!(tmp, np)
                end
            end
            nps = tmp
        end
        if remove_pi
            splice!(ps, i)
        else
            splice!(ps, i, unique!(nps))
        end
    end
    ifelse(isempty(ps), (false, nothing), (true, unique!(ps)))
end

function generate_ternamal(g::Generator, token::Token)
    f = get(g.𝐀, token, () -> string(token))
    f()
end

function generate(g::Generator, token::Token;max_depth::Int64=10,generator::Function=generate_dfs)
    token ∈ g.𝐕ₜ && return Statement[generate_ternamal(g, token),]
    (flag, states) = generator(g, token;max_depth)
    if !flag
        @warn "Generate fail."
        return [[string(:ϵ),]]
    end
    return map(states) do s
        join(map(x -> generate_ternamal(g, x), s), " ")
    end
end

max_stars_repo_name: yeruoforever/SimpleGrammar

max_issues_repo_name: yeruoforever/SimpleGrammar

max_forks_repo_name: yeruoforever/SimpleGrammar


text:
include("boundary_conditions_adapted.jl")
"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::ModelSetup)
apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::Union{ModelSetup},
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = state⁻.ρu

    return nothing
end

function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::Union{ModelSetup},
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0, -0])'

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::ModelSetup)
apply no slip boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for U
"""
@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::Union{ModelSetup},
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = -state⁻.ρu

    return nothing
end

@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::Union{ModelSetup},
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = -state⁻.ρu
    gradflux⁺.ν∇u = gradflux⁻.ν∇u

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Insulating, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for κ∇θ
"""
@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Union{Insulating, CoupledSecondaryBoundary, CoupledPrimaryBoundary},
    ::Union{ModelSetup},
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ

    return nothing
end

@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Union{Insulating, CoupledSecondaryBoundary, CoupledPrimaryBoundary},
    ::Union{ModelSetup},
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ
    gradflux⁺.κ∇θ = n⁻ * -0

    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Metatheory
using Metatheory.Classic
using Metatheory.EGraphs
using Metatheory.Library
using Metatheory.Util

using DataStructures
using Test

Metatheory.options.verbose = false
Metatheory.options.printiter = false

@metatheory_init ()

Test.FallbackTestSet(desc) = Test.FallbackTestSet()

ts = Test.FallbackTestSet

@testset ts "Metatheory Tests" begin
   @timev begin
      include("test_egraphs.jl")
      include("test_ematch.jl")
      include("test_ematch_assertions.jl")
      include("test_analysis.jl")
      include("test_extraction.jl")
      include("test_dynamic_ematch.jl")
      include("test_mu.jl")
      include("test_boson.jl")
      include("test_while_interpreter.jl")
      include("test_reductions.jl")
      include("test_taylor.jl")
      include("test_while_superinterpreter.jl")
      include("test_cat.jl")
      include("test_inequality.jl")
      include("test_custom_types.jl")
      include("test_multipat.jl")
      include("test_patallterm.jl")
      # use cases
      include("fib/test_fibonacci.jl")
      include("logic/test_calculational_logic.jl")
      include("logic/test_logic.jl")
      include("cas/test_infer.jl")
      include("cas/test_cas.jl")

      # exported consistency test
      for m ∈ [Metatheory, Metatheory.Util, Metatheory.Classic, Metatheory.EGraphs, Metatheory.EGraphs.Schedulers]
         for i ∈ propertynames(m)
            xxx = getproperty(m, i)
         end
      end
   end
end


text:
#!/usr/bin/env julia

using Couchzilla
using Base.Test

username = ENV["COUCH_USER"]
password = ENV["COUCH_PASS"]
host     = ENV["COUCH_HOST_URL"]

geo_database = ""
if haskey(ENV, "COUCH_GEO_DATABASE")
  geo_database = ENV["COUCH_GEO_DATABASE"]
end

database = "juliatest-$(Base.Random.uuid4())"
geo_database = "crimes"
cl = Client(username, password, host)
db, created = createdb(cl, database)

function testname(name::AbstractString)
  "\n\033[1m$name\033[0m" # Shell escape for bold text
end

try
  include("db_meta_tests.jl")
  include("crud_tests.jl")
  include("mango_tests.jl")
  include("replication_tests.jl")
  include("attachment_tests.jl")
  include("alldocs_tests.jl")
  include("view_tests.jl")
  include("retry_tests.jl")
  include("geo_tests.jl")
  include("apikey_tests.jl")

  # Include new tests above this line
  include("delete_db_tests.jl")
catch TestSetException
  println("\n[INFO] Deleting test database on failure")
  deletedb(cl, database)
end


text:
Base.@kwdef struct LinearAnalysis{I}
    npoints::Int=200
    dropcollinear::Bool=false
    interval::I=automatic
    level::Float64=0.95
end

function add_intercept_column(x::AbstractVector{T}) where {T}
    mat = similar(x, float(T), (length(x), 2))
    fill!(view(mat, :, 1), 1)
    copyto!(view(mat, :, 2), x)
    return mat
end

# TODO: add multidimensional version
function (l::LinearAnalysis)(input::ProcessedLayer)
    output = map(input) do p, n
        x, y = p
        weights = get(n, :weights, similar(x, 0))
        default_interval = length(weights) > 0 ? nothing : :confidence
        interval = l.interval === automatic ? default_interval : l.interval
        # FIXME: handle collinear case gracefully
        lin_model = GLM.lm(add_intercept_column(x), y; wts=weights, l.dropcollinear)
        x̂ = range(extrema(x)..., length=l.npoints)
        pred = GLM.predict(lin_model, add_intercept_column(x̂); interval, l.level)
        return if !isnothing(interval)
            # TODO: fix this in GLM https://github.com/JuliaStats/GLM.jl/pull/467
            ŷ, lower, upper = map(vec, pred) # GLM prediction returns matrices
            (x̂, ŷ), (; lower, upper)
        else
            ŷ = pred
            (x̂, ŷ), (;)
        end
    end
    default_plottype = isempty(output.named) ? Lines : LinesFill
    plottype = Makie.plottype(output.plottype, default_plottype)
    return ProcessedLayer(output; plottype)
end

"""
    linear(; interval=automatic, level=0.95, dropcollinear=false, npoints=200)

Compute a linear fit of `y ~ 1 + x`. An optional named mapping `weights` determines the weights.
Use `interval` to specify what type of interval the shaded band should represent,
for a given coverage `level` (the default `0.95` equates `alpha = 0.05`).
Valid values of `interval` are `:confidence`, to delimit the uncertainty of the predicted
relationship, and `:prediction`, to delimit estimated bounds for new data points.
Use `interval = nothing` to only compute the line fit, without any uncertainty estimate.
By default, this analysis errors on singular (collinear) data. To avoid that,
it is possible to set `dropcollinear=true`.
`npoints` is the number of points used by Makie to draw the shaded band.
"""
linear(; options...) = transformation(LinearAnalysis(; options...))


text:
using WordTokenizers

"""
    tokenize_document(path)

Helper function to prepare a text file for Word2Vec.  Takes in a filesystem path
as a String and returns an Array of tokenized sentences.
"""
function tokenize_document(path::T) where T<: AbstractString
    tokenized = open(path) do file
        text = replace(read(file, String), r"\n+|\r+|\s+" => " ")
        sentences = tokenize.(split_sentences(text))
        tokenized = Array{Array{AbstractString, 1}, 1}()
        for sentence in sentences
            no_punct = Array{AbstractString, 1}()
            for word in sentence
                occursin(r"[!\"#$%&\'()*+,-./:;<=>?@\[\\\]^_`\{\|\}~]+", word) ?
                                                continue : push!(no_punct, word)
            end
            push!(tokenized, lowercase.(no_punct))
        end
        tokenized
    end
    return tokenized
end


text:
#
#
#In the 20×20 grid below, four numbers along a diagonal line have been marked in red.

A = [
08 02 22 97 38 15 00 40 00 75 04 05 07 78 52 12 50 77 91 08
49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48 04 56 62 00
81 49 31 73 55 79 14 29 93 71 40 67 53 88 30 03 49 13 36 65
52 70 95 23 04 60 11 42 69 24 68 56 01 32 56 71 37 02 36 91
22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80
24 47 32 60 99 03 45 02 44 75 33 53 78 36 84 20 35 17 12 50
32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70
67 26 20 68 02 62 12 20 95 63 94 39 63 08 40 91 66 49 94 21
24 55 58 05 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72
21 36 23 09 75 00 76 44 20 45 35 14 00 61 33 97 34 31 33 95
78 17 53 28 22 75 31 67 15 94 03 80 04 62 16 14 09 53 56 92
16 39 05 42 96 35 31 47 55 58 88 24 00 17 54 24 36 29 85 57
86 56 00 48 35 71 89 07 05 44 44 37 44 60 21 58 51 54 17 58
19 80 81 68 05 94 47 69 28 73 92 13 86 52 17 77 04 89 55 40
04 52 08 83 97 35 99 16 07 97 57 32 16 26 26 79 33 27 98 66
88 36 68 87 57 62 20 72 03 46 33 67 46 55 12 32 63 93 53 69
04 42 16 73 38 25 39 11 24 94 72 18 08 46 29 32 40 62 76 36
20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74 04 36 16
20 73 35 29 78 31 90 01 74 31 49 71 48 86 81 16 23 57 05 54
01 70 54 71 83 51 54 69 16 92 33 48 61 43 52 01 89 19 67 48
]

#The product of these numbers is 26 × 63 × 78 × 14 = 1788696.
#
#What is the greatest product of four adjacent numbers in the same direction (up, down, left, right, or diagonally) in the 20×20 grid?

println(size(A))

maxval = typemin(Int64)
xi, xj, di = 0, 0, '.'
# horizontal
for i in 1:20
	for j in 1:20
		# horizontal
		if j <= 17
			x = prod(A[i, j:j+3])
			if x > maxval
				println(A[i, j:j+3])
				maxval, xi, xj, di = x, i, j, '>'
			end
		end
		# vertical
		if i <= 17
			x = prod(A[i:i+3, j])
			if x > maxval
				println(A[i:i+3, j])
				maxval, xi, xj, di = x, i, j, 'v'
			end
		end
		# top left -> bottom right
		if i <= 17 && j <= 17
			x = prod(diag(A[i:i+3, j:j+3]))
			if x > maxval
				println(diag(A[i:i+3, j:j+3]))
				maxval, xi, xj, di = x, i, j, '\\'
			end
		end
		# top right -> bottom left
		if i <= 17 && j <= 17
			x = prod(diag(rotr90(A[i:i+3, j:j+3])))
			if x > maxval
				println(diag(rotr90(A[i:i+3, j:j+3])))
				maxval, xi, xj, di = x, i, j, '/'
			end
		end
	end
end
print(maxval, ' ', xi, ' ', xj, ' ', di)
print(A[xi:xi+3, xj:xj+3])

# 89*94*97*87
# 70600674

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
push!(LOAD_PATH,"./")
using Fractal
using Images

fspecstr="0x08e961 0xe08e96 0xbc11d2 (0.00709466,-0.664701) zoom9.23647e-06 pow0.561485 affine1"
filename="output.png"
width=4096
height=4096
num_supersamples=2

camspec=parse(CameraSpec,fspecstr)
imagespec=ImageSpec(width, height, num_supersamples)

fracspec=FractalSpec(create_fractal("mandelbrot"))
fracgen=FractalGenerator(camspec,imagespec,fracspec)
@time img=generate_fractal(fracgen)
save(filename,img)

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal


text:
## Loading of inverse-model stuff

"""
Returns a default theta0 initialization dict.
"""
function theta0_dict_defaults(pm, n_1d_btilde, n_1d_fsl, n_1d_temp; with_sigma=false)
    out = OrderedDict(:btilde=> zeros(F,n_1d_btilde), # NOTE in m/y (error)
                      :fsl=> zeros(F,n_1d_fsl), # error
                      :temp=> zeros(F,n_1d_temp), # error
                      :dist_exp=> pm.dist_exp,
                       :iv_h_exp=> pm.iv_h_exp,
                      :iv_dist_exp1=> pm.iv_dist_exp1,
                      :iv_dist_exp2=> pm.iv_dist_exp2,
                      )
    if with_sigma
        out[:sigma_h_model] = pm.sigma_h_model
        out[:sigma_iv_model] = pm.sigma_iv_model
    end
    return out
end

"""
The default ball radius to be used in the emcee sampler initial conditions.
The size is chosen to be in line with the range of allowed values in priors.jl.
"""
const ball_radius_default = Dict(:btilde=> 1.4, # corresponds to sigma of 1 for dhdt and 1 for bdot
                                 :fsl=> 0.5,
                                 :temp=> 15.0,
                                 :dist_exp=> 0.8,
                                 :iv_h_exp=> 15.0, # @assert iv_h_exp>=0
                                 :iv_dist_exp1=> 0.4, # @assert 0<=iv_dist_exp1<=1
                                 :iv_dist_exp2=> 15.0, # @assert iv_dist_exp2>0
                                 :sigma_h_model => 100.0,
                                 :sigma_iv_model => 100.0,
                                 )
"""
    pmcmc_dict_defaults(theta0, runtyp=[:test,:mid,:prod][1])

Return, hopefully, sensible defaults for the construction of MCMCNum type with

    MCMCNum(pmcmc_dict_defaults(theta0, runtyp)...; other_kws)

The second argument allows to choose between quick-test, slower-test and production run
settings.
"""
function pmcmc_dict_defaults(theta0, runtyp)
    out = Dict(:sigma_ppdf => 0.001 + zeros(theta0.th0),
               :ball_radius => [ball_radius_default[name] for (name, inds) in theta0.names for i in inds]
               )
    dof = length(theta0.th0)
    return if runtyp==:test
        nchains = dof+2
        merge(out, Dict(:niter=>11*nchains,
                        :nthin => 1,
                        :nburnin=>10,
                        :nchains=>nchains))
    elseif runtyp==:testmid
        nchains = max(2*(dof+2), 30)
        merge(out, Dict(:niter=>10^4,
                        :nburnin=>9000,
                        :nchains=>nchains))
    elseif runtyp==:prodlow
        nchains = max(2*(dof+2), 50)
        merge(out, Dict(:niter=>10^5,
                        :nthin => 20,
                        :nburnin => 7*10^4,
                        :nchains=>nchains))
    elseif runtyp==:prodmid
        nchains = max(2*(dof+2), 50)
        merge(out, Dict(:niter=>5*10^5,
                        :nthin => 20,
                        :nburnin => 4*10^5,
                        :nchains=>nchains))
    elseif runtyp==:prodhigh
        nchains = min(10*dof, max(dof+2, 50))
        merge(out, Dict(:niter=>10^6,
                        :nchains=>nchains,
                        :nthin => 50,
                        :nburnin => 7*10^5))
    elseif runtyp==:prod
        nchains = min(10*dof, max(dof+2, 50))
        merge(out, Dict(:niter=>3*10^5,
                        :nchains=>nchains,
                        :nthin => 20,
                        :nburnin => 2*10^5))
    else
        error("Unknown option `runtyp`==$runtyp")
    end
end

"""
$(SIGNATURES)

Initialize all inverse-model parameters.

Returns
- theta0
- logposterior
- logprior, logprior_debug
- loglikelihood
- fwdm_fn
- pmcmc-dict-defaults
- fit_target_eff: what run-type was effectively chosen (depending on data availability)
"""
init_inverse(gb::Bands, pp, pm, pn; kwargs...) =
    init_inverse(gb.gl.gid, gb, pp, pm, pn; kwargs...)

function init_inverse(gid::GlacierID, gb::Bands, pp, pm, pn;
                      runtyp=[:test,:mid,:prod][1],
                      n_1d_btilde = 3, n_1d_fsl = 3, n_1d_temp = 1,
                      fit_target=FitTarget.h_iv,
                      theta0_dict=theta0_dict_defaults(pm, n_1d_btilde, n_1d_fsl, n_1d_temp),
                      check_pdf=true, # if true check and try to update if -Inf
                      do_grid_serach=true
                      )
    theta0_dict = deepcopy(theta0_dict)
    if pm.iv_nye
        pop!(theta0_dict, :iv_h_exp, nothing)
        pop!(theta0_dict, :iv_dist_exp1, nothing)
    else
        pop!(theta0_dict, :iv_dist_exp2, nothing)
    end

    ## Modify fit_target:
    if fit_target==FitTarget.h_iv
        # - if there is no IV data, then drop iv-part
        if gb.gl.iv isa BITEModel.SomeData{Void}
            fit_target = FitTarget.h
        end
        # - if there is no h data, then drop h-part
        if gb.gl.h isa BITEModel.SomeData{Void}
            fit_target = FitTarget.iv
        end
    end
    #  - if there is no data to fit to, the only fit the glacier length
    if  (fit_target==FitTarget.h && gb.gl.h isa BITEModel.SomeData{Void}) ||
        (fit_target==FitTarget.iv && gb.gl.iv isa BITEModel.SomeData{Void})
        fit_target = FitTarget.length
    end

    if fit_target==FitTarget.length || fit_target==FitTarget.prior
        pop!(theta0_dict, :sigma_h_model, nothing)
        pop!(theta0_dict, :sigma_iv_model, nothing)
    elseif fit_target==FitTarget.iv
        pop!(theta0_dict, :sigma_h_model, nothing)
    elseif fit_target==FitTarget.h
        pop!(theta0_dict, :sigma_iv_model, nothing)
    elseif  fit_target==FitTarget.h_iv
        nothing
    else
        error("Unknow $(fit_target)")
    end

    theta0 = Theta0(theta0_dict, gb, pp, pm, pn)

    logprior, logprior_debug = make_logprior(theta0)
    pn.test && Base.Test.@inferred logprior(theta0.th0)

    # forward model
    loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn = make_loglikelihood(theta0, pp, pm, pn, fit_target)
    pn.test && Base.Test.@inferred fwdm_fn(theta0.th0)
    pn.test && Base.Test.@inferred loglikelihood(theta0.th0)
    pn.test && Base.Test.@inferred loglikelihood1d(theta0.th0)

    logposterior, logposterior1d = make_logposterior(loglikelihood, loglikelihood1d, logprior)
    pn.test && Base.Test.@inferred logposterior(theta0.th0)

    # Try hard to get a btilde which works decently
    if check_pdf && logposterior(theta0.th0)[1]==-Inf || do_grid_serach
        orig = theta0_dict[:btilde]

        # Note: only fit length (through prior for calving glaciers and through 1d likelihood otherwise)

        if length(orig)!=3
            # Do an abbreviated grid search
            best = -Inf
            delta_btilde_best = zeros(orig)
            for delta_btilde = logspace(0, 1.0, 10) - 0.95
                for plusminus = [-1,1] # add or subtract delta_btilde
                    for ind=0:2^(length(orig))-1
                        # index vector
                        vec = [parse.(Bool, c) for c in bin(ind, length(orig))]
                        theta0_dict[:btilde] = orig + plusminus * delta_btilde .* vec
                        theta0 = Theta0(theta0_dict, gb, pp, pm, pn)
                        lp = logposterior1d(theta0.th0)[1]
                        if lp>best
                            best = lp
                            delta_btilde_best .= plusminus * delta_btilde .* vec
                        end
                    end
                end
            end
        else
            # Do a grid search
            best = -Inf
            delta_btilde_best = [0.0, 0.0, 0.0]
            nn = 5
            delta_btildes = [-(logspace(0, 1.0, nn) - 0.95); 0 ; logspace(0, 1.0, nn) - 0.95]
            for di in delta_btildes
                for dj in delta_btildes
                    for dk in delta_btildes
                        theta0_dict[:btilde] = orig .+ [di, dj, dk]
                        theta0 = Theta0(theta0_dict, gb, pp, pm, pn)
                        lp = logposterior1d(theta0.th0)[1]
                        if lp>best
                            best = lp
                            delta_btilde_best = [di, dj, dk]
                        end
                    end
                end
            end
        end
        # use the delta_btilde_best
        theta0_dict[:btilde] = orig .+ delta_btilde_best
        theta0 = Theta0(theta0_dict, gb, pp, pm, pn)

        @assert logposterior(theta0.th0)[1]>-Inf "Could not find a theta0 such that logposterior>-Inf"
    end

    return (theta0, logposterior, logposterior1d, logprior, logprior_debug,
            loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
            pmcmc_dict_defaults(theta0, runtyp), fit_target)
end


text:
# Examples used in docs, tests, and benchmarking
module Examples
using Parameters, Symbolics
using DifferentiableStateSpaceModels

include("rbc.jl")
include("sgu.jl")
include("FVGQ20.jl")
end

# Macro to load example can't be in that module.  Needs to be macro due to world-age issues

"""
Return model from a stored module, creating as necessary .  If the module is
- already included, it just returns the `Module`
- saved to disk but not already included, it does `include`
- not saved to disk, it saves it, includes it

The specification of a `model_generator` has a sepcification like `H, mod_vals = model_generator()`
where `H` is an symbolics expression the `mod_vals` can be splatted into the `PerturbationModel` constructor.

The name of the stored `module` is the name of `model_generator` function and the order.

$(SIGNATURES)

# Examples
```julia-repl
julia> m = @include_example_module(Examples.rbc_observables_benchmark)  # stores or loads `.function_cache/rbc_observables.jl`
"""
macro include_example_module(model_generator)
    return quote
        function_name = "$(Base.nameof($(model_generator)))"  # strips out module name
        local model_name = "$($model_generator)"
        local model_filename = "$(model_name).jl"
        local model_cache_path = joinpath(DifferentiableStateSpaceModels.default_model_cache_location(),
                                          model_filename)
        if !isfile(model_cache_path)
            local H, mod_vals = $model_generator()
            make_perturbation_model(H; model_name, mod_vals...)
            include(model_cache_path)
        elseif !isdefined(Main, Symbol(model_name))
            include(model_cache_path)
        end
        local mod = getfield(Main, Symbol(model_name))
        DifferentiableStateSpaceModels.PerturbationModel(mod)
    end |> esc
end


text:
using Plots, StatsPlots
using DataFrames, FileIO, CSVFiles
# gr()
pyplot()

palette = [:red, :green, :blue, :cyan, :black]
noise_levels = [20, 40, 60, 80, 100]

root = isdir("results") ? "." : "benchmark"
result_dir = joinpath(root, "results")
experiments = filter!(readdir(result_dir, join=true)) do expdir
    isdir(expdir) || return false
    any(endswith(".csv"), readdir(expdir))
end

default(
    titlefont = (15, "times"),
    legendfontsize = 10,
    size=(1000, 700),
    guidefont = (15, :black),
    dpi=300,
    # tickfont = (10, :orange),
    # guide = "x",
    # framestyle = :zerolines,
    yminorgrid = true
)

for expdir in experiments
    expname = basename(expdir)
    infofile = joinpath(expdir, "versioninfo.txt")
    @assert isfile(infofile) "versioninfo.txt does not exist"
    m = match(r"\d+", filter!(contains("JULIA_NUM_THREADS"), readlines(infofile))[1])
    @assert !isnothing(m) "versioninfo.txt file must contain JULIA_NUM_THREADS info"
    num_threads = m.match

    julia_results = map(noise_levels) do σ
        julia_file = joinpath(expdir, "julia_$σ.csv")
        df = DataFrame(load(julia_file))
        accumulate!(+, df[!, :runtime], df[!, :runtime])
        df
    end
    matlab_results = map(noise_levels) do σ
        matlab_file = joinpath(expdir, "matlab_$σ.csv")
        df = DataFrame(load(matlab_file))
        accumulate!(+, df[!, :runtime], df[!, :runtime])
        df
    end

    p1 = plot(
        title="$expname (n=$num_threads)",
        ylabel="Time (seconds) in log scale",
        xlabel="PSNR(dB)",
        legend=:bottomright,
        xlim=(10, 50),
    )
    for i in 1:length(noise_levels)
        julia_df = julia_results[i]
        matlab_df = matlab_results[i]
        c = palette[i]
        σ = noise_levels[i]

        plot!(p1, matlab_df[!, :psnr], matlab_df[!, :runtime], marker=:circle, label="MATLAB(σ=$σ)", yscale=:log10, c=c, line=:dash)
        plot!(p1, julia_df[!, :psnr], julia_df[!, :runtime], marker=:square, label="Julia(σ=$σ)", yscale=:log10, c=c, line=:solid)
    end
    p1

    julia_runtime = map(julia_results) do df
        last(df[!, :runtime])
    end
    matlab_runtime = map(matlab_results) do df
        last(df[!, :runtime])
    end

    p2 = plot(
        title="Runtime",
        xlabel="Noise levels",
        ylabel="Runtime (seconds)",
        legend=:topleft
    )
    ctg = repeat(["Julia", "MATLAB"], inner=5)
    nam = repeat(noise_levels, outer = 2)
    groupedbar!(p2, nam, [julia_runtime matlab_runtime]; group=ctg)

    p3 = plot(
        title="Runtime (Scaled)",
        ylabel="runtime ratio",
        xlabel="Noise levels",
        legend=:topleft
    )
    groupedbar!(p3, nam, [ones(size(matlab_runtime)) matlab_runtime./julia_runtime]; group=ctg)

    p = plot(p1, plot(p2, p3, layout=(2, 1)), layout=(1, 2))
    outfile = joinpath(result_dir, basename(expdir)*".png")
    savefig(p, outfile)
end


text:
"""
    $SIGNATURES

Internal function to apply a vector of `Drawing` objects contained in an `Axes` container in a GLE
context. The `axidx` and `figid` help keep track of where individual drawings belong to which is
useful when writing auxiliary files containing the drawing data.
"""
function apply_drawings!(g::GLE, drawings::Vector{<:Drawing}, figid::String, axidx::Int)
    # element counter to have an index over objects drawn
    el_cntr = 1
    for drawing ∈ drawings
        el_cntr = apply_drawing!(g, drawing, el_cntr, figid, axidx)
    end
    # this is recuperated in `apply_axes!` and further processed by `apply_legend!`
    return nothing
end

"""
    $SIGNATURES

Internal function to generate a path to an auxiliary file storing drawing data for the current
axes of the current figure. The argument `h` is the hash of what has to be written (so that if
the exact same command is issued, the file is not re-written). The `axidx` and `figid` are there
to make the auxiliary file unique and easy to delete after use.
"""
auxpath(h::UInt, fid::String, axid::Int) = joinpath(GP_ENV["TMP_PATH"], "$(fid)_$(axid)_$h.csv")

####
#### Apply a Scatter2D object
####

"""
    $SIGNATURES

Internal function to apply a `Drawing` object `obj` in a GLE context `g` with current legend
entries `leg_entries` (possibly empty), current element counter `el_cntr` and `axidx` and `figid`
are used to make the auxiliary file name unique.
"""
function apply_drawing!(g::GLE, scatter::Scatter2D, el_cntr::Int, figid::String, axidx::Int)
    # write data to a temporary CSV file
    faux = auxpath(hash(scatter.data), figid, axidx)
    dir  = splitdir(faux)[1]
    # don't rewrite if it's the exact same zipper
    isfile(faux) || csv_writer(faux, scatter.data, scatter.hasmissing)

    #
    # GLE syntax is:
    #
    #   data datafile.dat d1=c1,c2
    #   d1 line blue
    #

    for k ∈ 1:scatter.nobj
        # (1) indicate what data to read
        "\n\tdata \"$faux\" d$(el_cntr)=c1,c$(k+1)" |> g
        # if no color has been specified, assign one according to the PALETTE
        if !isdef(scatter.linestyles[k].color)
            cc = mod(el_cntr, GP_ENV["SIZE_PALETTE"])
            (cc == 0) && (cc = GP_ENV["SIZE_PALETTE"])
            scatter.linestyles[k].color = GP_ENV["PALETTE"][cc]
        end
        # (2) line and marker description
        # build a tuple with the current buffer and the legend entry buffer
        if scatter.linestyles[k].lstyle != -1
            # Line plot
            "\n\td$el_cntr" |> g
            "line" |> g; apply_linestyle!(g, scatter.linestyles[k])
            # if a marker color is specified and different than the line
            # color we need to have a special subroutine for GLE
            mcol = false
            mscale = 1
            if isdef(scatter.markerstyles[k].color) &&
                    (scatter.markerstyles[k].color != scatter.linestyles[k].color)
                mcol = true
                f = Figure(figid; _noreset=true)
                add_sub_marker!(f, scatter.markerstyles[k])
                mscale = f.textstyle.hei
            end
            apply_markerstyle!(g, scatter.markerstyles[k], mcol=mcol, mscale=mscale)
        else
            # Scatter plot; if there's no specified marker color,
            # take the default line color
            if !isdef(scatter.markerstyles[k].color)
                scatter.markerstyles[k].color = scatter.linestyles[k].color
            end
            "\n\td$el_cntr" |> g
            apply_markerstyle!(g, scatter.markerstyles[k])
        end
        el_cntr += 1
    end
    return el_cntr
end

####
#### Apply a Fill2D object
####

function apply_drawing!(g::GLE, fill::Fill2D, el_cntr::Int, figid::String, axidx::Int)
    # write data to a temporary CSV file
    faux = auxpath(hash(fill.data), figid, axidx)
    isfile(faux) || csv_writer(faux, fill.data, false)

    #
    # GLE syntax is:
    #
    #   data datafile.dat d1=c1,c2 d2=c1,c3
    #   fill d1,d2 color rgb(1,1,1) xmin 0 xmax 1
    #

    "\n\tdata \"$faux\" d$(el_cntr)=c1,c2 d$(el_cntr+1)=c1,c3" |> g
    "\n\tfill d$(el_cntr),d$(el_cntr+1)" |> g
    "color $(col2str(fill.fillstyle.fill))"    |> g
    isdef(fill.xmin) && "xmin $(fill.xmin)"    |> g
    isdef(fill.xmax) && "xmax $(fill.xmax)"    |> g

    el_cntr += 2
    return el_cntr
end

####
#### Apply a Hist2D object
####

function apply_drawing!(g::GLE, hist::Hist2D, el_cntr::Int, figid::String, axidx::Int)
    # write data to a temporary CSV file
    faux = auxpath(hash(hist.data), figid, axidx)
    isfile(faux) || csv_writer(faux, hist.data, hist.hasmissing)

    #
    # GLE syntax is:
    #
    #   data datafile.dat d1
    #   let d2 = hist d1 from xmin to xmax bins nbins
    #   bar d2 width width$ fill col$ color col2$ pattern pat$ horiz
    #

    # if no color has been specified, assign one according to the PALETTE
    if !isdef(hist.barstyle.color)
        if hist.barstyle.fill == c"white"
            cc = mod(el_cntr, GP_ENV["SIZE_PALETTE"])
            (cc == 0) && (cc = GP_ENV["SIZE_PALETTE"])
            hist.barstyle.color = GP_ENV["PALETTE"][cc]
        else
            hist.barstyle.color = c"white" # looks nicer than black
        end
    end

    # (1) indicate what data to read
    "\n\tdata \"$faux\" d$(el_cntr)" |> g

    # (2) hist description
    minx, maxx = hist.range
    "\n\tlet d$(el_cntr+1) = hist d$(el_cntr)" |> g
    "from $minx to $maxx" |> g
    el_cntr += 1

    # number of bins
    nobs = hist.nobs
    nbauto = (nobs == 0 ? 1 : ceil(Integer, log2(nobs))+1) # sturges, see StatsBase
    bins = isdef(hist.bins) ? hist.bins : nbauto
    "bins $bins" |> g

    # (3) compute appropriate scaling
    width   = (maxx - minx) / bins
    scaling = 1.0
    hist.scaling == "probability" && (scaling /= nobs)
    hist.scaling == "pdf"         && (scaling /= (nobs * width))
    "\n\tlet d$(el_cntr) = d$(el_cntr)*$scaling" |> g

    # (4) apply histogram
    "\n\tbar d$(el_cntr) width $width" |> g

    # apply styling
    apply_barstyle!(g, hist.barstyle)
    hist.horiz && "horiz" |> g

    return el_cntr+1
end

####
#### Apply Bar2D
####

function apply_drawing!(g::GLE, bar::Bar2D, el_cntr::Int, figid::String, axidx::Int)
    # write data to a temporary CSV file
    faux = auxpath(hash(bar.data), figid, axidx)
    isfile(faux) || csv_writer(faux, bar.data, bar.hasmissing)

    # >>>>>>>>>>>>>>>>
    # general GLE syntax is:
    # (1) data datafile.dat d1 d2 ...
    # (2) bar d1,d2,d3 fill color_ color color_
    # (STACKED, NO GROUP)
    # bar d1 ...
    # bar d2 from d1 ...
    # XXX (STACKED, GROUP) ---> see first whether this has a use case...
    # could solve it by specifying what gets stacked stack=[(1, 2), ...] but
    # clunky (and tbf not great visualisation anyway...)
    # bar d1,d2 ...
    # bar d3,d4 from d1,d2 ...
    # <<<<<<<<<<<<<<<<

    # if no color has been specified, assign one according to the PALETTE
    for c ∈ eachindex(bar.barstyles)
        if !isdef(bar.barstyles[c].color)
            if bar.barstyles[c].fill == c"white"
                cc = mod(el_cntr+c-1, GP_ENV["SIZE_PALETTE"])
                (cc == 0) && (cc = GP_ENV["SIZE_PALETTE"])
                bar.barstyles[c].color = GP_ENV["PALETTE"][cc]
            else
                bar.barstyles[c].color = c"white"
            end
        end
    end

    nbars = bar.nobj

    # (1) indicate what data to read "data file d1 d2 d3..."
    "\n\tdata \"$faux\"" |> g
    prod("d$(el_cntr+i-1) " for i ∈ (1:nbars)) |> g

    # (2) non stacked (or single barset)
    if nbars==1 || !bar.stacked
        # bar d1,d2,d3
        "\n\tbar $(svec2str(("d$(el_cntr+i-1)" for i ∈ 1:nbars)))" |> g
        isdef(bar.bwidth) && "width $(bar.bwidth)" |> g
        # apply bar styles
        apply_barstyles_nostack!(g, bar.barstyles)
        bar.horiz && "horiz" |> g

    # (2) stacked
    else
        # first base bar
        "\n\tbar d$(el_cntr)" |> g
        isdef(bar.bwidth) && "width $(bar.bwidth)" |> g
        apply_barstyle!(g, bar.barstyles[1])
        bar.horiz && "horiz" |> g
        # bars stacked on top
        for i ∈ 2:nbars
            "\n\tbar d$(el_cntr+i-1) from d$(el_cntr+i-2)" |> g
            apply_barstyle!(g, bar.barstyles[i])
        end
    end

    return el_cntr + nbars
end


text:
using BasisFunctions, Test, DoubleFloats

for T in (Float64,LargeFloat)

    @testset "Discrete Gauss OPS Rule" begin
        for d in (Legendre{T}(10), Laguerre(10,rand(T)), Hermite{T}(10), ChebyshevT{T}(10),ChebyshevU{T}(10),Jacobi(10,rand(T),rand(T)),Fourier{T}(10))
            g = interpolation_grid(d)
            gr = gauss_rule(d)
            @test g ≈ points(gr)
            @test isorthogonal(d, gr)
            @test iscompatible(d, g)
            @test gram(d, gr) isa DiagonalOperator

            test_orthogonality_orthonormality(d, true, d isa Fourier, gr)

            @test BasisFunctions.quadweights(gauss_rule(d),measure(d)) ≈ ones(T, 10)
        end

        B = Fourier(10)
        d = B^2
        g = interpolation_grid(d)
        gr = gauss_rule(d)
        @test g ≈ points(gr)
        @test isorthogonal(d, gr)
        @test iscompatible(d, g)
        @test gram(d, gr) isa TensorProductOperator && all(map(x->isa(x, DiagonalOperator), components(gram(d,gr))))
        @test BasisFunctions.quadweights(gauss_rule(d),measure(d)) ≈ BasisFunctions.OuterProductArray(map(ones, axes(B^2))...)

    end
end


text:
#=
Kinematics
=#

using SpaceTimeFields
import SpaceTimeFields: Abstract1DProfile, >>, ConstantProfile, d_dt

struct Constant{C <: Complex, A <: Real} <: Kinematics
    ċ::C
    α̇::A
end
Constant(ċ, α̇) = Constant(complex(ċ), α̇)
(c::Constant{C})(t) where C = zero(C), c.ċ, zero(C), 0.0, c.α̇, 0.0
show(io::IO, c::Constant) = print(io, "Constant (ċ = $(c.ċ), α̇ = $(c.α̇))")

"""
    Pitchup(U₀,a,K,α₀,t₀,Δα,ramp=EldredgeRamp) <: Kinematics

Kinematics describing a pitch-ramp motion (horizontal translation with rotation)
starting at time ``t_0`` about an axis at `a` (expressed relative to the centroid, in the ``\\tilde{x}``
  direction in the body-fixed coordinate system), with translational velocity `U₀`
in the inertial ``x`` direction, initial angle ``\\alpha_0``, and angular
change ``\\Delta\\alpha``. The optional ramp argument is assumed to be
given by the smooth ramp `EldredgeRamp`, though this can be changed to
`ColoniusRamp`.
"""
struct Pitchup <: Kinematics
    "Freestream velocity"
    U₀::Float64
    "Axis of rotation, relative to the plate centroid"
    a::Float64

    "Non-dimensional pitch rate ``K = \\dot{\\alpha}_0\\frac{c}{2U_0}``"
    K::Float64

    "Initial angle of attack"
    α₀::Float64
    "Nominal start of pitch up"
    t₀::Float64

    "Total pitching angle"
    Δα::Float64

    α::Abstract1DProfile
    α̇::Abstract1DProfile
    α̈::Abstract1DProfile
end

function Pitchup(U₀, a, K, α₀, t₀, Δα, ramp=EldredgeRamp)
    Δt = 0.5Δα/K
    p = ConstantProfile(α₀) + 2K*((ramp >> t₀) - (ramp >> (t₀ + Δt)))
    ṗ = d_dt(p)
    p̈ = d_dt(ṗ)

    Pitchup(U₀, a, K, α₀, t₀, Δα, p, ṗ, p̈)
end

function (p::Pitchup)(t)
    α = p.α(t)
    α̇ = p.α̇(t)
    α̈ = p.α̈(t)

    c = p.U₀*t - p.a*exp(im*α)
    ċ = p.U₀ - p.a*im*α̇*exp(im*α)
    if (t - p.t₀) > p.Δα/p.K
        c̈ = 0.0im
    else
        c̈ = p.a*exp(im*α)*(α̇^2 - im*α̈)
    end

    return c, ċ, c̈, α, α̇, α̈
end

function show(io::IO, p::Pitchup)
    print(io, "Pitch-up kinematics with rate K = $(p.K)")
end

"""
    Oscillation(Ux,Uy,α̇₀,ax,ay,Ω,Ax,Ay,ϕx,ϕy,α₀,Δα,ϕα) <: Kinematics

Set 2-d oscillatory kinematics. This general constructor sets up motion of a
rotational axis, located at `ax`, `ay` (expressed relative to the body centroid, in
  a body-fixed coordinate system). The rotational axis motion is described by

``
x(t) = U_x t + A_x\\sin(\\Omega t - \\phi_x), \\quad y(t) = U_y t + A_y\\sin(\\Omega t - \\phi_y),
 \\quad \\alpha(t) = \\alpha_0 + \\dot{\\alpha}_0 t +
 \\Delta\\alpha \\sin(\\Omega t - \\phi_{\\alpha})
 ``

"""
struct Oscillation <: Kinematics
    "Translational velocity in x"
    Ux::Float64

    "Translational velocity in y"
    Uy::Float64

    "Constant rotation rate"
    α̇₀::Float64

    "Axis of pitch rotation, in body coordinates"
    ax::Float64
    ay::Float64

    "Angular frequency"
    Ω::Float64

    "Amplitude of oscillatory motion in x"
    Ax::Float64

    "Amplitude of oscillatory motion in y"
    Ay::Float64

    "Phase of x motion (in radians)"
    ϕx::Float64

    "Phase of y motion (in radians)"
    ϕy::Float64

    "Mean angle of attack"
    α₀::Float64

    "Amplitude of oscillatory rotation (in radians)"
    Δα::Float64

    "Phase of rotational motion (in radians)"
    ϕα::Float64

    px::Abstract1DProfile
    ṗx::Abstract1DProfile
    p̈x::Abstract1DProfile

    py::Abstract1DProfile
    ṗy::Abstract1DProfile
    p̈y::Abstract1DProfile

    α::Abstract1DProfile
    α̇::Abstract1DProfile
    α̈::Abstract1DProfile
end

function Oscillation(Ux, Uy, α̇₀, ax, ay, Ω, Ax, Ay, ϕx, ϕy, α₀, Δα, ϕa)
    Ω > 0.0 || error("Frequency must be positive and non-zero")
    px = Ax*(Sinusoid(Ω) >> (ϕx/Ω))
    ṗx = d_dt(px)
    p̈x = d_dt(ṗx)
    py = Ay*(Sinusoid(Ω) >> (ϕy/Ω))
    ṗy = d_dt(py)
    p̈y = d_dt(ṗy)
    α = ConstantProfile(α₀) + Δα*(Sinusoid(Ω) >> (ϕa/Ω))
    α̇ = d_dt(α)
    α̈ = d_dt(α̇)
    Oscillation(Ux, Uy, α̇₀, ax, ay, Ω, Ax, Ay, ϕx, ϕy, α₀, Δα, ϕa, px, ṗx, p̈x, py, ṗy, p̈y, α, α̇, α̈)
end

function (p::Oscillation)(t)
    α = p.α̇₀*t + p.α(t)
    α̇ = p.α̇₀   + p.α̇(t)
    α̈ = p.α̈(t)

    a = complex(p.ax,p.ay)
    U = complex(p.Ux,p.Uy)
    aeiα = a*exp(im*α)
    c = U*t + p.px(t) + im*p.py(t) - aeiα
    ċ = U + p.ṗx(t) + im*p.ṗy(t) - im*α̇*aeiα
    c̈ = p.p̈x(t) + im*p.p̈y(t) + aeiα*(α̇^2 - im*α̈)

    return c, ċ, c̈, α, α̇, α̈
end

function show(io::IO, p::Oscillation)
    println(io, "Oscillatory kinematics with")
    println(io, "     Steady velocity U = ($(p.Ux),$(p.Uy))")
    println(io, "     Ref angle α₀ = $(p.α₀)")
    println(io, "     Mean rotation rate α̇₀ = $(p.α̇₀)")
    println(io, "     Pitch axis (rel. to centroid) a = ($(p.ax),$(p.ay))")
    println(io, "     Frequency Ω = $(p.Ω)")
    println(io, "     x amplitude Ax, phase lag ϕx = ($(p.Ax), $(p.ϕx))")
    println(io, "     y amplitude Ay, phase lag ϕy = ($(p.Ay), $(p.ϕy))")
    println(io, "     α amplitude Δα, phase lag ϕα = ($(p.Δα), $(p.ϕα))")
end

"""
    PitchHeave(U₀,a,Ω,α₀,Δα,ϕp,A,ϕh)

Create oscillatory pitching and heaving kinematics of a pitch axis at location ``a`` (expressed
relative to the centroid in the ``\\tilde{x}`` direction of the body-fixed coordinate system),
of the form (in inertial coordinates)

``
x(t) = U_0, \\quad y(t) = A\\sin(\\Omega t - \\phi_h),
 \\quad \\alpha(t) = \\alpha_0 + \\Delta\\alpha \\sin(\\Omega t - \\phi_p)
 ``
"""
PitchHeave(U₀, a, Ω, α₀, Δα, ϕp, A, ϕh) = Oscillation(U₀, 0, 0, a, 0, Ω, 0, A, 0, ϕh, α₀, Δα, ϕp)

"""
    OscillationXY(Ux,Uy,Ω,Ax,ϕx,Ay,ϕy)

Set oscillatory kinematics in the ``x`` and ``y`` directions, of the form

``
x(t) = U_x t + A_x \\sin(\\Omega t - \\phi_x), \\quad y(t) = U_y t + A_y \\sin(\\Omega t - \\phi_y)
``
"""
OscillationXY(Ux,Uy,Ω,Ax,ϕx,Ay,ϕy) = Oscillation(Ux, Uy, 0, 0, 0, Ω, Ax, Ay, ϕx, ϕy, 0, 0, 0)

"""
    OscillationX(Ux,Ω,Ax,ϕx)

Set oscillatory kinematics in the ``x`` direction, of the form

``
x(t) = U_x t + A_x \\sin(\\Omega t - \\phi_x),
``
"""
OscillationX(Ux,Ω,Ax,ϕx) = Oscillation(Ux, 0, 0, 0, 0, Ω, Ax, 0, ϕx, 0, 0, 0, 0)

"""
    OscillationY(Uy,Ω,Ay,ϕy)

Set oscillatory kinematics in the ``y`` direction, of the form

``
y(t) = U_y t + A_y \\sin(\\Omega t - \\phi_y)
``
"""
OscillationY(Uy,Ω,Ay,ϕy) = Oscillation(0, Uy, 0, 0, 0, Ω, 0, Ay, 0, ϕy, 0, 0, 0)

"""
    RotationalOscillation(ax,ay,Ω,α₀,α̇₀,Δα,ϕα)

Set oscillatory rotational kinematics about an axis located at `ax`, `ay`
(expressed relative to the body centroid, in a body-fixed coordinate system), of the form

``
\\alpha(t) = \\alpha_0 + \\dot{\\alpha}_0 t +
\\Delta\\alpha \\sin(\\Omega t - \\phi_{\\alpha})
``
"""
RotationalOscillation(ax,ay,Ω,α₀,α̇₀,Δα,ϕα) = Oscillation(0, 0, α̇₀, ax, ay, Ω, 0, 0, 0, 0, α₀, Δα, ϕα)

"""
    RotationalOscillation(Ω,Δα,ϕα)

Set oscillatory rotational kinematics about the centroid of the form

``
\\alpha(t) = \\Delta\\alpha \\sin(\\Omega t - \\phi_{\\alpha})
``
"""
RotationalOscillation(Ω,Δα,ϕα) = RotationalOscillation(0,0,Ω,0,0,Δα,ϕα)

abstract type Switch end
abstract type SwitchOn <: Switch end
abstract type SwitchOff <: Switch end

"""
    SwitchedKinematics <: Kinematics

Modulates a given set of kinematics between simple on/off states. The velocity
specified by the given kinematics is toggled on/off.

# Fields
$(FIELDS)
"""
struct SwitchedKinematics{S <: Switch} <: Kinematics

    "time at which the kinematics should be turned on"
    t_on :: Float64

    "time at which the kinematics should be turned off"
    t_off :: Float64

    "kinematics to be followed in the on state"
    kin :: Kinematics

    off :: Kinematics

    SwitchedKinematics(t_on,t_off,kin) = t_on > t_off ?
            new{SwitchOn}(t_on,t_off,kin,RigidBodyMotions.Constant(0,0)) :
            new{SwitchOff}(t_on,t_off,kin,RigidBodyMotions.Constant(0,0))
end

# note that these do not introduce impulsive changes into the derivatives
(p::SwitchedKinematics{SwitchOn})(t) = t <= p.t_on ? p.off(t) : p.kin(t-p.t_on)

(p::SwitchedKinematics{SwitchOff})(t) = t <= p.t_off ? p.kin(t-p.t_on) : p.off(t)


text:
module CentroidLinearClassifierImpl

export train_model

using LinearAlgebra
using ..Clusters, ..LinearClassifiers

"""
    train_model(positive_data, negative_data)

Trains a binary classification model.

The returned model has a field `classify`, a function that classifies a point
and returns a scalar value, where positive and negative results indicate that
the point is positively or negatively classified, respectively.
"""
function train_model(positive_data::AbstractMatrix{<:Number},
                     negative_data::AbstractMatrix{<:Number})::CentroidLinearClassifier
    @assert size(positive_data, 2) == size(negative_data, 2)

    positive_data_centroid = compute_centroid(positive_data)
    negative_data_centroid = compute_centroid(negative_data)

    normal_vector = positive_data_centroid .- negative_data_centroid
    midpoint = (positive_data_centroid .+ negative_data_centroid) ./ 2

    offset = -1 * dot(normal_vector, midpoint)

    decision_function =
        function (x)
            return dot(normal_vector, x) + offset
        end

    return CentroidLinearClassifier(decision_function, normal_vector, offset)
end

end # module

max_stars_repo_name: brianxie/automaton-edification

max_issues_repo_name: brianxie/automaton-edification

max_forks_repo_name: brianxie/automaton-edification


text:
for i in 1:10
    print(i)
    if i%5 == 0
        println()
        continue
    end
    print(", ")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApimachineryPkgRuntimeRawExtension <: SwaggerModel
    Raw::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: Raw

    function IoK8sApimachineryPkgRuntimeRawExtension(;Raw=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgRuntimeRawExtension, Symbol("Raw"), Raw)
        setfield!(o, Symbol("Raw"), Raw)
        o
    end
end # type IoK8sApimachineryPkgRuntimeRawExtension

const _property_map_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,Symbol}(Symbol("Raw")=>Symbol("Raw"))
const _property_types_IoK8sApimachineryPkgRuntimeRawExtension = Dict{Symbol,String}(Symbol("Raw")=>"Vector{UInt8}")
Base.propertynames(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }) = collect(keys(_property_map_IoK8sApimachineryPkgRuntimeRawExtension))
Swagger.property_type(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApimachineryPkgRuntimeRawExtension[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgRuntimeRawExtension[property_name]

function check_required(o::IoK8sApimachineryPkgRuntimeRawExtension)
    (getproperty(o, Symbol("Raw")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgRuntimeRawExtension }, name::Symbol, val)
    if name === Symbol("Raw")
    end
end


text:
const use_local_dependencies = Ref(false)
const use_local_plotlyjs = Ref(false)

function _init_ijulia_plotting()
    # IJulia is more stable with local file
    use_local_plotlyjs[] = plotly_local_file_path[] === nothing ? false : isfile(plotly_local_file_path[])

    ENV["MPLBACKEND"] = "Agg"
end

"""
Add extra jupyter mimetypes to display_dict based on the plot backed.

The default is nothing, except for plotly based backends, where it
adds data for `application/vnd.plotly.v1+json` that is used in
frontends like jupyterlab and nteract.
"""
_ijulia__extra_mime_info!(plt::Plot, out::Dict) = out

function _ijulia__extra_mime_info!(plt::Plot{PlotlyJSBackend}, out::Dict)
    out["application/vnd.plotly.v1+json"] = Dict(
        :data => plotly_series(plt),
        :layout => plotly_layout(plt)
    )
    out
end

function _ijulia__extra_mime_info!(plt::Plot{PlotlyBackend}, out::Dict)
    out["application/vnd.plotly.v1+json"] = Dict(
        :data => plotly_series(plt),
        :layout => plotly_layout(plt)
    )
    out
end

function _ijulia_display_dict(plt::Plot)
    output_type = Symbol(plt.attr[:html_output_format])
    if output_type == :auto
        output_type = get(_best_html_output_type, backend_name(plt.backend), :svg)
    end
    out = Dict()
    if output_type == :txt
        mime = "text/plain"
        out[mime] = sprint(show, MIME(mime), plt)
    elseif output_type == :png
        mime = "image/png"
        out[mime] = base64encode(show, MIME(mime), plt)
    elseif output_type == :svg
        mime = "image/svg+xml"
        out[mime] = sprint(show, MIME(mime), plt)
    elseif output_type == :html
        mime = "text/html"
        out[mime] = sprint(show, MIME(mime), plt)
        _ijulia__extra_mime_info!(plt, out)
    else
        error("Unsupported output type $output_type")
    end
    out
end


text:
type AggressivenessBelief <: BehaviorBelief
    gen::CorrelatedIDMMOBIL
    physical::MLPhysicalState
    particles::Vector{Vector{Float64}} # First index is the position in physical.cars
    weights::Vector{Vector{Float64}}   # Second index is the particle number
end

function rand(rng::AbstractRNG,
                 b::AggressivenessBelief,
                # s::MLState=MLState(b.physical, Array(CarState, length(b.physical.cars))))
                 s::MLState=MLState(b.physical, Array{CarState}(length(b.physical.cars))))
    return rand(rng, b, zeros(length(b.physical.cars)), s)
end

function rand(rng::AbstractRNG,
                b::AggressivenessBelief,
                sample_noises::Vector{Float64},
                s::MLState=MLState(b.physical, Array{CarState}(length(b.physical.cars))))
                #s::MLState=MLState(b.physical, Array(CarState, length(b.physical.cars))))

    s.x = b.physical.x
    s.t = b.physical.t
    resize!(s.cars, length(b.physical.cars))
    for i in 1:length(s.cars)
        particle = sample(rng, b.particles[i], WeightVec(b.weights[i]))
        nudged = max(min(particle+sample_noises[i]*randn(rng),1.0),0.0)
        s.cars[i] = CarState(b.physical.cars[i], create_model(b.gen, nudged))
    end
    s.terminal = b.physical.terminal
    return s
end

# action(p::Policy, b::AggressivenessBelief) = action(p, b.physical)

function most_likely_state(b::AggressivenessBelief)
    #s = MLState(b.physical, Array(CarState, length(b.physical.cars)))
    s = MLState(b.physical, Array{CarState}(length(b.physical.cars)))
    for i in 1:length(s.cars)
        ml_ind = indmax(b.weights[i])
        behavior = create_model(b.gen, b.particles[i][ml_ind])
        s.cars[i] = CarState(b.physical.cars[i], behavior)
    end
    return s
end
agg_means(b::AggressivenessBelief) = [mean(b.particles[i], WeightVec(b.weights[i])) for i in 1:length(b.particles)]
function agg_stds(b::AggressivenessBelief)
    means = agg_means(b)
    #stds = Array(Float64, length(b.physical.cars))
    stds = Array{Float64}(length(b.physical.cars))
    for i in 1:length(b.physical.cars)
        stds[i] = sqrt(sum(b.weights[i].*(b.particles[i]-means[i]).^2)/sum(b.weights[i]))
    end
    return stds
end

function weights_from_particles!(b::AggressivenessBelief,
                                 problem::NoCrashProblem,
                                 o::MLPhysicalState,
                                 particles,
                                 p::WeightUpdateParams)

    b.physical = o
    resize!(b.weights, length(o.cars))
    resize!(b.particles, length(o.cars))
    for i in 1:length(o.cars)
        # make sure we're not going to be allocating a bunch of memory in the loop
        #if isdefined(b.particles, i)
        if isassigned(b.particles, i)
            sizehint!(b.particles[i], length(particles))
            resize!(b.particles[i], 0)
        else
            #b.particles[i] = Array(Float64, length(particles))
            b.particles[i] = Array{Float64}(length(particles))
            resize!(b.particles[i], 0)
        end
        if isassigned(b.weights, i)
            sizehint!(b.weights[i], length(particles))
            resize!(b.weights[i], 0)
        else
            #b.weights[i] = Array(Float64, length(particles))
            b.weights[i] = Array{Float64}(length(particles))
            resize!(b.weights[i], 0)
        end
    end
    for (j, sp) in enumerate(particles)
        isp = 1
        io = 1
        while io <= length(o.cars) && isp <= length(sp.cars)
            co = o.cars[io]
            csp = sp.cars[isp]
            if co.id == csp.id
                if abs(co.x-csp.x) < 0.2*problem.dmodel.phys_param.lane_length
                    @assert isa(csp.behavior, IDMMOBILBehavior)
                    a = csp.behavior.p_idm.a
                    dt = problem.dmodel.phys_param.dt
                    veld = TriangularDist(csp.vel-a*dt/2.0, csp.vel+a*dt/2.0, csp.vel)
                    proportional_likelihood = Distributions.pdf(veld, co.vel)
                    if co.y == csp.y
                        push!(b.particles[io], aggressiveness(b.gen, csp.behavior))
                        push!(b.weights[io], proportional_likelihood)
                    elseif abs(co.y - csp.y) <= 1.1
                        push!(b.particles[io], aggressiveness(b.gen, csp.behavior))
                        push!(b.weights[io], p.wrong_lane_factor*proportional_likelihood)
                    end # if greater than one lane apart, do nothing
                end
                io += 1
                isp += 1
            elseif co.id < csp.id
                io += 1
            else
                @assert co.id > csp.id
                isp += 1
            end
        end
    end

    return b
end

mutable struct AggressivenessUpdater <: Updater #{AggressivenessBelief}
    problem::Nullable{NoCrashProblem}
    nb_sims::Int
    p_resample_noise::Float64
    resample_noise_factor::Float64
    params::WeightUpdateParams
    rng::AbstractRNG
end
function set_problem!(u::AggressivenessUpdater, p::Union{POMDP,MDP})
    u.problem = Nullable{NoCrashProblem}(p)
end
function set_rng!(u::AggressivenessUpdater, rng::AbstractRNG)
    u.rng = rng
end

function update(up::AggressivenessUpdater,
                b_old::AggressivenessBelief,
                a::MLAction,
                o::MLPhysicalState,
                b_new::AggressivenessBelief=AggressivenessBelief(CorrelatedIDMMOBIL(
                                                    get(up.problem).dmodel.behaviors), o,
                                                    Array{Vector{Float64}}(length(o.cars)),
                                                    Array{Vector{Float64}}(length(o.cars))))
                                                    #Array(Vector{Float64}, length(o.cars)),
                                                    #Array(Vector{Float64}, length(o.cars))))

    #particles = Array(MLState, up.nb_sims)
    particles = Array{MLState}(up.nb_sims)
    stds = max(agg_stds(b_old), 0.01)
    for i in 1:up.nb_sims
        #=
        if rand(up.rng) < up.p_resample_noise
            s = rand(up.rng, b_old, up.resample_noise_factor*stds)
        else
            s = rand(up.rng, b_old)
        end
        =#
        s = rand(up.rng, b_old)
        particles[i] = generate_s(get(up.problem), s, a, up.rng)
    end

    weights_from_particles!(b_new, get(up.problem), o, particles, up.params)

    for i in 1:length(o.cars)
        if isempty(b_new.weights[i])
            # println("car $i has empty weights")
            b_new.particles[i] = rand(up.rng, up.nb_sims)
            b_new.weights[i] = ones(up.nb_sims)
        end
    end

    return b_new
end

function initialize_belief(up::AggressivenessUpdater, distribution)
    gen = CorrelatedIDMMOBIL(get(up.problem).dmodel.behaviors)
    states = [rand(up.rng, distribution) for i in 1:up.nb_sims]
    #particles = Array(Array{Float64}, length(first(states).cars))
    particles = Array{Array{Float64}}(length(first(states).cars))
    #weights = Array(Array{Float64}, length(first(states).cars))
    weights = Array{Array{Float64}}(length(first(states).cars))#
    for i in 1:length(first(states).cars)
        #particles[i] = Array(Float64, length(states))
        particles[i] = Array{Float64}(length(states))
        weights[i] = ones(length(states))
        for (j,s) in enumerate(states)
            particles[i][j] = aggressiveness(gen, s.cars[i].behavior)
        end
    end
    return AggressivenessBelief(gen, MLPhysicalState(first(states)), particles, weights)
end

initialize_belief(up::AggressivenessUpdater, d::AggressivenessBelief) = d

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
function DataFrames.rename!(df, ::DataFile{Val{Symbol("C&D")}})
    namemap = Dict("WE(1).Potential (V)"=>"Potential (V)",
                   "Time (s)"=>"Other Time (s)",
                   "Corrected time (s)"=>"Time (s)")
    rename!(df, namemap)
end

function find_pair(datafile, list, ending)
    findfirst(f -> filevalue(f) == filevalue(datafile) &&
                   foldervalue(f) == foldervalue(datafile) &&
                   endswith(f.filename, ending),
              list)
end

function add_CD(filename)
    parts = rsplit(filename, '.', limit=2)
    parts[1][1:end-1] * "CD." * parts[2]
end

function process_data(::Val{Symbol("C&D")}, data; insert_D, continue_col)
    done = Vector{Int}()
    for (i,f) in enumerate(data["C&D"])
        if i in done
            continue
        else
            push!(done, i)

            df = read_file(f)

            valid_idx = setdiff(axes(data["C&D"], 1), done)
            pair_idx = find_pair(f, data["C&D"], endswith(f.filename, "_C") ? "_D" : "_C")

            if isnothing(pair_idx)
                if endswith(f.filename, "_D")
                    pushfirst(df, insert_D)
                    write_file(f, df, ';')
                end
                continue
            end

            f_pair = data["C&D"][pair_idx]
            push!(done, pair_idx)

            df_pair = read_file(f_pair)
            df_C, df_D = endswith(f.filename, "_C") ? (df, df_pair) : (df_pair, df)
            df_CD, df_D = postprocess(f, df_C, df_D, insert_D, continue_col)

            new_name = add_CD(f.savename)
            mergedf = DataFile{Val{Symbol("C&D")}}(f.filename, new_name, f.units,
                f.legend_units, f.idx)

            write_file(mergedf, df_CD, ';')
            write_file(endswith(f.filename, "_D") ? f : f_pair, df_D, ';')
        end
    end

    return nothing
end

function pushfirst(df, value)
    types = eltype.(eachcol(df))
    line = (Vector{types[i]}([v]) for (i,v) in enumerate(value))
    to_add = DataFrame(;zip(propertynames(df), line)...)
    append!(to_add, df)
end

function postprocess(datafile::DataFile{Val{Symbol("C&D")}}, df_C, df_D, value, cont_col)
    # Insert values in _D file
    df_D_mod = pushfirst(df_D, value)
    # Append column with types
    df_C[!, :Type] .= "C"
    df_D[!, :Type] .= "D"
    # Add last value of cont_col form _C to _D values
    last_time = df_C[end, cont_col]
    df_D[!, cont_col] .+= last_time
    # Append DataFrames
    append!(df_C, df_D)
    return df_C, df_D_mod
end


text:
# ---
# title: 744. Find Smallest Letter Greater Than Target
# id: problem744
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Binary Search
# link: <https://leetcode.com/problems/find-smallest-letter-greater-than-target/description/>
# hidden: true
# ---
#
# Given a list of sorted characters `letters` containing only lowercase letters,
# and given a target letter `target`, find the smallest element in the list that
# is larger than the given target.
#
# Letters also wrap around. For example, if the target is `target = 'z'` and
# `letters = ['a', 'b']`, the answer is `'a'`.
#
# **Examples:**
#
#
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "a"
#     Output: "c"
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "c"
#     Output: "f"
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "d"
#     Output: "f"
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "g"
#     Output: "j"
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "j"
#     Output: "c"
#
#     Input:
#     letters = ["c", "f", "j"]
#     target = "k"
#     Output: "c"
#
#
# **Note:**
#
#   1. `letters` has a length in range `[2, 10000]`.
#   2. `letters` consists of lowercase letters, and contains at least 2 unique letters.
#   3. `target` is a lowercase letter.
#
#
## @lc code=start
using LeetCode

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


text:
"""
    cholUnblocked!(A, Val{:L})

Overwrite the lower triangle of `A` with its lower Cholesky factor.

The name is borrowed from [https://github.com/andreasnoack/LinearAlgebra.jl]
because these are part of the inner calculations in a blocked Cholesky factorization.
"""
function cholUnblocked! end

function cholUnblocked!(D::Diagonal{T}, ::Type{Val{:L}}) where {T<:AbstractFloat}
    map!(sqrt, D.diag, D.diag)
    D
end

function cholUnblocked!(A::StridedMatrix{T}, ::Type{Val{:L}}) where {T<:BlasFloat}
    n = LinearAlgebra.checksquare(A)
    if n == 1
        A[1] < zero(T) && throw(PosDefException(1))
        A[1] = sqrt(A[1])
    elseif n == 2
        A[1] < zero(T) && throw(PosDefException(1))
        A[1] = sqrt(A[1])
        A[2] /= A[1]
        (A[4] -= abs2(A[2])) < zero(T) && throw(PosDefException(2))
        A[4] = sqrt(A[4])
    else
        _, info = LAPACK.potrf!('L', A)
        iszero(info) || throw(PosDefException(info))
    end
    A
end

function cholUnblocked!(D::UniformBlockDiagonal, ::Type{Val{:L}})
    for f in D.facevec
        cholUnblocked!(f, Val{:L})
    end
    D
end


text:
# Signal to Revise.jl that this file should be tracked as a package
#__revise_mode__ = :eval

# Types
abstract type AbstractGrid end

struct Grid{N, T} <: AbstractGrid
    r₀::SVector{N, T}
    l::SVector{N, T}
    Δl::SVector{N, T}
    size::SVector{N, Int}
    periodicity::SVector{N, Bool}

    function Grid(lower, upper, size, periodicity)
        return Grid(lower, upper, SVector(size), periodicity)
    end

    function Grid(lower, upper, size::SVector{N, <: Integer}, periodicity) where {N}
        l = SVector(upper .- lower)
        Δ = l ./ size
        T = eltype(Δ)
        return new{N, T}(SVector(lower), l, Δ, size, SVector(periodicity))
    end
end

function Base.show(io::IO, g::Grid{N, T}) where {N, T}
    context = IOContext(io, :compact => true)
    join(io, g.size, '×')
    print(io, "-element Grid{", N, ", ", T, "}:\n ")
    itr = Iterators.Stateful(zip(g.r₀, g.l))
    for (lo, l) in itr
        print(context, '[', lo, " .. ", lo + l, ']')
        if Base.peek(itr) !== nothing
            print(io, " × ")
        end
    end
    return nothing
end

function circ_index(r, r₀, l, Δl, m)
    n = floor(Int, mod(r - r₀, l) / Δl) + 1
    return ifelse(n > m, 1, n)
end

function Base.getindex(g::Grid{N}, r::Union{Real, NTuple{N, Real}, SVector{N}}) where {N}
    return CartesianIndex(Tuple(circ_index.(r, g.r₀, g.l, g.Δl, g.size)))
end

Base.size(g::Grid) = Tuple(g.size)
Base.ndims(::Grid{N}) where {N} = N


text:
using Luxor

function drawmatrix(A::Matrix;
    cellsize = (10, 10))
    table = Table(size(A)..., cellsize...)
    used = Set()
    for i in CartesianIndices(A)
        r, c = Tuple(i)
        if A[r, c] ∈ used
            sethue("orange")
        else
            sethue("purple")
            push!(used, A[r, c])
        end
        text(string(A[r, c]), table[r, c],
            halign = :center,
            valign = :middle)
        sethue("white")
        box(table, r, c, :stroke)
    end
end

function draw_matrix()
    A = rand(1:99, 5, 8)
    @draw begin
        background("black")
        fontsize(30)
        setline(0.5)
        sethue("white")
        drawmatrix(A, cellsize = 10 .* size(A))
    end
end

itworks() = print("It works !")

struct MatrixChunks
    X::Matrix{<:AbstractFloat}
    n::Int  # size of matrix
end

function MatrixChunks(X::Matrix{<:AbstractFloat})
    (n1, n2) = size(X)
    if n1 != n2
        error("Matrix should be square")
    end
    return MatrixChunks(X, n1)
end

function Base.iterate(iter::MatrixChunks)
    element = view(iter.X, 1:2, 1:2)
    return (element, (1, 1))
end

function Base.iterate(iter::MatrixChunks, state)
    if state[1] + 1 < iter.n
        x = state[1] + 1
        y = state[2]
    else
        if state[2] + 1 < iter.n
            x = 1
            y = state[2] + 1
        else
            return nothing
        end
    end

    element = view(iter.X, y:y+1, x:x+1)

    return (element, (x, y))
end

function iterate_iso_field(a::Matrix{<:AbstractFloat})
    for y = 1:size(a, 1)-1
        for x = 1:size(a, 2)-1
            s = view(a, y:y+1, x:x+1)
            display((x, y))
            display(s)
        end
    end
end

function generate_iso_field()
    a = randn(4, 4)

    display(a)

    iterate_iso_field(a)

    println()

    for n in MatrixChunks(a)
        display(n)
    end
end

max_stars_repo_name: vansweej/MetaBalls

max_issues_repo_name: vansweej/MetaBalls

max_forks_repo_name: vansweej/MetaBalls


text:
using Images

url = "https://atmire.com/dspace-labs3/bitstream/handle/123456789/7618/earth-map-huge.jpg?sequence=1"

const savedir = joinpath(tempdir(), "Images")
if !isdir(savedir)
    mkdir(savedir)
end

const savename = joinpath(savedir, "earth-map-huge.jpg")
if !isfile(savename)
    download(url, savename)
end

img = imread(savename)


text:
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cartesian})::T where {T <: SSDFloat}
    return wp(pt.x, pt.y, pt.z)::T
end
@inline function get_interpolation(wp::Interpolations.Extrapolation{T, 3}, pt::CartesianPoint{T}, S::Val{:cylindrical})::T where {T <: SSDFloat}
    p::CylindricalPoint{T} = CylindricalPoint(pt)
    return wp(p.r, p.φ, p.z)::T
end

function add_signal!(signal::Vector{T}, path::Vector{CartesianPoint{T}}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    @inbounds for i in eachindex(signal)
        signal[i] += get_interpolation(wp, path[i], S)::T * charge
    end
    nothing
end
function add_signal!(signal::Vector{T}, path::DriftPath{T}, charge::T, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    add_signal!(signal, path.e_path, -charge, wp, S) # electrons induce negative charge
    add_signal!(signal, path.h_path,  charge, wp, S)
    nothing
end
function add_signal!(signal::Vector{T}, paths::Vector{DriftPath{T}}, charges::Vector{T}, wp::Interpolations.Extrapolation{T, 3}, S::Union{Val{:cylindrical}, Val{:cartesian}})::Nothing where {T <: SSDFloat}
    for ipath in eachindex(paths)
        add_signal!(signal, paths[ipath], charges[ipath], wp, S)
    end
    nothing
end

function signal_contributions_from_drift_paths(drift_paths, energy_depositions::AbstractVector{T}, Wpot_interp::Interpolations.Extrapolation{T,3}) where T<:Real
    charge_signal_e = zeros(T,size(drift_paths[1].e_path,1))
    charge_signal_h = zeros(T,size(drift_paths[1].e_path,1))
    @inbounds for i in eachindex(drift_paths[1].e_path)
        for (idp, drift_path) in enumerate(drift_paths)
            # charge_signal_e[i] = muladd(energy_depositions[idp] ,charge_signal_e[i], _get_wpot_at(Wpot_interp, drift_path.e_path[i]))
            # charge_signal_h[i] = muladd(energy_depositions[idp] ,charge_signal_h[i], _get_wpot_at(Wpot_interp, drift_path.h_path[i]))
            charge_signal_e[i] +=  _get_wpot_at(Wpot_interp, drift_path.e_path[i]) * energy_depositions[idp]
            charge_signal_h[i] +=  _get_wpot_at(Wpot_interp, drift_path.h_path[i]) * energy_depositions[idp]
        end
    end
    return charge_signal_e .* -1, charge_signal_h
end
@deprecate signal_contributions_from_drift_paths(drift_paths, energy_depositions::AbstractVector{T}, Wpot_interp::Interpolations.Extrapolation{T,3}) where {T <: SSDFloat} get_signal(dp::DriftPath{T}, energy::T, wp::Interpolations.Extrapolation{T,3,ITPT,IT,ET} where ET where IT where ITPT) where {T <: SSDFloat}