text:
mutable struct Cluster
	points::Vector{Tuple{Int64,Int64}}
	boundary::Rect
	r::Int64
	x_bound::Int64
	y_bound::Int64
end

Cluster(r::Int64,x_bound::Int64,y_bound::Int64) = Cluster([], Rect([ [0,0] [0,0] ]),r,x_bound,y_bound)

function add_point_to_cluster(cluster::Cluster,p::Tuple{Int64,Int64})
	if p in cluster.points
		return false
	end

	push!(cluster.points,p)
	update_rect(cluster::Cluster,p::Tuple{Int64,Int64})
	return true
end

function update_rect(cluster::Cluster,p::Tuple{Int64,Int64})
	if length(cluster.points) == 1
		cluster.boundary = Rect( bounded_diff(p[1],cluster.r,0),
							  	 bounded_diff(p[2],cluster.r,0),
							 	 bounded_sum(p[1],cluster.r,cluster.x_bound),
							  	 bounded_sum(p[2],cluster.r,cluster.y_bound))
		return
	end

	x1,x2 = x(cluster.boundary)
	y1,y2 = y(cluster.boundary)

	if length(cluster.points) == 2
		cluster.boundary = Rect(x1,y1,
								 bounded_sum(p[1],cluster.r,cluster.x_bound),
								 bounded_sum(p[2],cluster.r,cluster.y_bound))
		return
	end

	x1 = min(x1, bounded_diff(p[1],cluster.r,0))
	x2 = max(x2, bounded_sum(p[1],cluster.r,cluster.x_bound))
	y1 = min(y1, bounded_diff(p[2],cluster.r,0))
	y2 = max(y2, bounded_sum(p[2],cluster.r,cluster.y_bound))

	cluster.boundary = Rect(x1,y1,x2,y2)
end

function is_inside_cluster(cluster::Cluster,p::Tuple{Int64,Int64})
	return inside(p,cluster.boundary)
end


text:
function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    as = [parse(Int, x) for x in split(readline())]
    ret = maximum(as) + 1

    pos = [[0] for _ in 1:ret]
    for (i, a) in enumerate(as)
        push!(pos[a+1], i)
    end
    for v in values(pos)
        push!(v, n+1)
    end

    for (k, v) in pairs(pos)
        for i in 1:length(v)-1
            if v[i+1] - v[i] > m
                ret = min(ret, k-1)
            end
        end
    end
    ret
end

println(solve())

ext: jl

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

# Based on the tutorial found at https://biojulia.net/Automa.jl/latest/#Overview-1

# Describe patterns in regular expression.
oct      = re"0o[0-7]+"
dec      = re"[-+]?[0-9]+"
# Hex numbers are actually present without the `0x` prefix, but the
# regexp without `0x` would conflict with `dec`.
hex      = re"[0-9A-Fa-f]+"
prefloat = re"[-+]?([0-9]+\.[0-9]*|[0-9]*\.[0-9]+)"
float    = prefloat | re.cat(prefloat | re"[-+]?[0-9]+", re"[eE][-+]?[0-9]+")
number   = oct | dec | hex | float
numbers  = re.cat(re.opt(number), re.rep(re" +" * number), re" *")

# # Register action names to regular expressions.
# number.actions[:enter] = [:mark]
# oct.actions[:exit]     = [:oct]
# dec.actions[:exit]     = [:dec]
# hex.actions[:exit]     = [:hex]
# float.actions[:exit]   = [:float]

open_paren = re"\("
close_paren = re"\)"
paren = open_paren | close_paren

# White space is correctly identified, but for some reason not
# captured.
whitespace = re"[ \n]+"

delim = paren | whitespace
delim.actions[:enter] = [:mark]
open_paren.actions[:exit] = [:open_paren]
close_paren.actions[:exit] = [:close_paren]
whitespace.actions[:exit] = [:whitespace]

str = re"\".*\""
str.actions[:enter] = [:mark]
str.actions[:exit] = [:str]

symbol = re"[a-zA-Z\*$%]"*re.rep(re"[^ ()]")
# symbol.actions[:enter] = [:mark]
# symbol.actions[:exit] = [:symbol]

num_or_sym = number | symbol
num_or_sym.actions[:enter] = [:mark]
oct.actions[:exit]     = [:oct]
dec.actions[:exit]     = [:dec]
hex.actions[:exit]     = [:hex]
float.actions[:exit]   = [:float]
symbol.actions[:exit] = [:symbol]

# kicad_lisp_tokens = number | symbol | str
kicad_lisp_tokens = num_or_sym | str

kicad_lisp = re.cat(re.opt(kicad_lisp_tokens),
                    re.rep(delim|delim*kicad_lisp_tokens),
                    re" *")

# # Compile a finite-state machine.
machine = Automa.compile(kicad_lisp)

# # This generates a SVG file to visualize the state machine.
# write("kicad_lisp.dot", Automa.machine2dot(machine))
# run(`dot -Tpng -o kicad_lisp.png kicad_lisp.dot`)

# Bind an action code for each action name.
actions = Dict(
    :mark  => :(mark = p),
    :oct   => :(emit(:oct)),
    :dec   => :(emit(:dec)),
    :hex   => :(emit(:hex)),
    :float => :(emit(:float)),
    :open_paren => :(emit(:open_paren)),
    :close_paren => :(emit(:close_paren)),
    :str => :(emit(:str)),
    :symbol => :(emit(:symbol)),
    :whitespace => :(emit(:whitespace))
)

# Generate a tokenizing function from the machine.
context = Automa.CodeGenContext()
@eval function tokenize(data::String)
    tokens = Tuple{Symbol,String}[]
    mark = 0
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    emit(kind) = push!(tokens, (kind, data[mark:p-1]))
    $(Automa.generate_exec_code(context, machine, actions))
    return tokens, cs == 0 ? :ok : cs < 0 ? :error : :incomplete
end

abstract type KObj end

mutable struct KNode <: KObj
    children::Vector{KObj}
end
KNode() = KNode(KObj[])
KNode(objs::KObj...) = KNode([objs...])

Base.copy(n::KNode) = KNode(copy(n.children))

function Base.show(io::IO, node::KNode)
    write(io, "(")
    col = 1
    max_col = 60
    for (i,c) in enumerate(node.children)
        i > 1 && write(io, " ")
        sc = split(string(c), "\n")
        if col + length(sc[1]) ≥ max_col
            write(io, "\n  ")
            col = 1
        end
        write(io, sc[1])
        for l in sc[2:end]
            write(io, "\n  ", l)
        end
        length(sc) != 1 && (col += length(sc[end]))
        if col >= max_col && i ∉ [1,length(node.children)]
            write(io, "\n  ")
            col = 1
        end
    end
    write(io, ")")
end

struct KSym <: KObj
    name::String
end
KSym(name::Symbol) = KSym(string(name))
Base.show(io::IO, s::KSym) = write(io, s.name)
macro ksym_str(name)
    KSym(name)
end

struct KStr <: KObj
    str::String
    KStr(str::String) = new(strip(str, '"'))
end
Base.show(io::IO, s::KStr) = write(io, "\"", s.str, "\"")

abstract type KNum <: KObj end

struct KInt{I<:Integer} <: KNum
    i::I
end
Base.show(io::IO, i::KInt) = write(io, string(i.i))
function Base.show(io::IO, i::KInt{<:Unsigned})
    if i.i == 0
        write(io, "0")
    else
        s = repr(i.i)
        write(io, lstrip(uppercase(s[3:end]), '0'))
    end
end

struct KFloat{F<:AbstractFloat} <: KNum
    f::F
end
Base.show(io::IO, f::KFloat) = write(io, string(f.f))

obj(o::Symbol) = KSym(o)
obj(o::String) = KStr(o)
obj(o::I) where {I<:Integer} = KInt(o)
obj(o::F) where {F<:AbstractFloat} = KFloat(o)
obj(v::Vector) = obj.(v)

function Parse(tokens::Vector{Tuple{Symbol,String}})
    nodes = KNode[]
    cur_node = KNode()
    for (token,value) in tokens
        token == :whitespace && continue
        if token == :open_paren
            push!(nodes, cur_node)
            new_node = KNode()
            push!(cur_node.children, new_node)
            cur_node = new_node
        elseif token == :close_paren
            cur_node = try
                pop!(nodes)
            catch
                error("Unbalanced parentheses")
            end
        else
            node = if token == :symbol
                KSym(value)
            elseif token == :str
                KStr(value)
            elseif token == :dec
                KInt(parse(Int, value))
            elseif token == :hex
                KInt(parse(UInt32, "0x$(value)"))
            elseif token == :float
                KFloat(parse(Float64, value))
            else
                throw(ArgumentError("Unknown token $(token)"))
            end
            push!(cur_node.children, node)
        end
    end
    cur_node.children
end

function Read(s::String)
    tokens,status = tokenize(s)
    status != :ok && error("Tokenization failed: $(status) [$(tokens)]")
    Parse(tokens)
end

Read(io::IO) = Read(read(io, String))

function traverse(visit::Function, tree::KNode)
    stack = KObj[tree]
    while !isempty(stack)
        node = pop!(stack)
        node isa KNode && append!(stack, node.children)
        visit(node)
    end
end

function traverse_nodes(visit::Function, tree::KNode)
    traverse(tree) do obj
        obj isa KNode && visit(obj)
    end
end

function traverse_nodes(visit::Function, tree::KNode, name::String)
    traverse_nodes(tree) do node
        fc = first(node.children)
        fc.name == name && visit(node)
    end
end

export @ksym_str


text:
module BenchFuzzy

import LinearAlgebra
using SyncBarriers
using BenchmarkTools
using UnPack: @unpack

function mul_naive!(
    y::AbstractVector,
    A::AbstractMatrix,
    x::AbstractVector,
    a = true,
    b = false,
)
    (axes(y, 1) == axes(A, 1) && axes(A, 2) == axes(x, 1)) || throw(DimensionMismatch())
    for i in axes(A, 1)
        s = zero(eltype(y))
        @simd for j in axes(A, 2)
            s += @inbounds A[i, j] * x[j]  # A is adjoint or transpose
        end
        y[i] = s * a + y[i] * b
    end
    return y
end

function sim_seq(model, niters, mul! = mul_naive!)
    @unpack f, g, A, B, C, x0, y0 = model
    x1 = copy(x0)
    x2 = similar(x0)
    y1 = copy(y0)
    y2 = similar(y0)
    for _ in 1:niters
        mul!(x2, A, x1)
        mul!(y2, C, x1)
        mul!(y2, B, y1, true, true)
        x2 .= f.(x2)
        y2 .= g.(y2)
        x1, x2 = x2, x1
        y1, y2 = y2, y1
    end
    return (x1, y1)
end

function sim_parallel_fuzzy(
    model,
    niters,
    bx,
    by,
    parallel_foreach,
    spin,
    mul! = mul_naive!,
)
    @assert length(bx) == length(by)
    @unpack f, g, A, B, C, x0, y0 = model
    x1 = copy(x0)
    x2 = similar(x0)
    y1 = copy(y0)
    y2 = similar(y0)
    xs = (x1, x2)
    ys = (y1, y2)

    ntasks = length(bx)
    xchunks = collect(Iterators.partition(eachindex(x0), cld(length(x0), ntasks)))
    ychunks = collect(Iterators.partition(eachindex(y0), cld(length(y0), ntasks)))
    @assert length(xchunks) == length(ychunks) == ntasks
    parallel_foreach(1:ntasks) do i
        for t in 1:niters
            local x1 , x2 , y1 , y2
            if isodd(t)
                x1, x2 = xs
                y1, y2 = ys
            else
                x2, x1 = xs
                y2, y1 = ys
            end
            @views begin
                xl = x2[xchunks[i]]
                yl = y2[ychunks[i]]
                Al = A[xchunks[i], :]
                Bl = B[ychunks[i], :]
                Cl = C[ychunks[i], :]
            end
            mul!(xl, Al, x1)
            t == 1 || depart!(by[i], spin)
            mul!(yl, Cl, x1)
            xl .= f.(xl)
            arrive!(bx[i])  # => x1 writable, x2 readable
            mul!(yl, Bl, y1, true, true)
            yl .= g.(yl)
            arrive!(by[i])  # => y1 writable, y2 readable
            depart!(bx[i], spin)
        end
    end

    return (x1, y1)
end

function sim_parallel(model, niters, b, parallel_foreach, spin, mul! = mul_naive!)
    @unpack f, g, A, B, C, x0, y0 = model
    x1 = copy(x0)
    x2 = similar(x0)
    y1 = copy(y0)
    y2 = similar(y0)
    xs = (x1, x2)
    ys = (y1, y2)

    ntasks = length(b)
    xchunks = collect(Iterators.partition(eachindex(x0), cld(length(x0), ntasks)))
    ychunks = collect(Iterators.partition(eachindex(y0), cld(length(y0), ntasks)))
    @assert length(xchunks) == length(ychunks) == ntasks
    parallel_foreach(1:ntasks) do i
        for t in 1:niters
            local x1 , x2 , y1 , y2
            if isodd(t)
                x1, x2 = xs
                y1, y2 = ys
            else
                x2, x1 = xs
                y2, y1 = ys
            end
            @views begin
                xl = x2[xchunks[i]]
                yl = y2[ychunks[i]]
                Al = A[xchunks[i], :]
                Bl = B[ychunks[i], :]
                Cl = C[ychunks[i], :]
            end
            mul!(xl, Al, x1)
            xl .= f.(xl)
            mul!(yl, Cl, x1)
            mul!(yl, Bl, y1, true, true)
            yl .= g.(yl)
            cycle!(b[i], spin)
        end
    end

    return (x1, y1)
end

function sim_parallel_nobarrier(model, niters, ntasks, parallel_foreach, mul! = mul_naive!)
    @unpack f, g, A, B, C, x0, y0 = model
    x1 = copy(x0)
    x2 = similar(x0)
    y1 = copy(y0)
    y2 = similar(y0)
    xs = (x1, x2)
    ys = (y1, y2)

    xchunks = collect(Iterators.partition(eachindex(x0), cld(length(x0), ntasks)))
    ychunks = collect(Iterators.partition(eachindex(y0), cld(length(y0), ntasks)))
    @assert length(xchunks) == length(ychunks) == ntasks
    for t in 1:niters
        parallel_foreach(1:ntasks) do i
            local x1 , x2 , y1 , y2
            if isodd(t)
                x1, x2 = xs
                y1, y2 = ys
            else
                x2, x1 = xs
                y2, y1 = ys
            end
            @views begin
                xl = x2[xchunks[i]]
                yl = y2[ychunks[i]]
                Al = A[xchunks[i], :]
                Bl = B[ychunks[i], :]
                Cl = C[ychunks[i], :]
            end
            mul!(xl, Al, x1)
            xl .= f.(xl)
            mul!(yl, Cl, x1)
            mul!(yl, Bl, y1, true, true)
            yl .= g.(yl)
        end
    end

    return (x1, y1)
end

function parallel_foreach_static(f, xs)
    Threads.@threads :static for x in xs
        f(x)
    end
end

function parallel_foreach_dynamic(f, xs)
    tasks = empty!(Vector{Task}(undef, length(xs)))
    for x in xs
        t = Threads.@spawn f(x)
        push!(tasks, t)
    end
    foreach(wait, tasks)
end

function random_model(n, m, g = 1.5)
    At = (g / √n) .* rand(n, n)
    Bt = (g / √m) .* rand(m, m)
    Ct = (g / √n) .* rand(n, m)
    return (
        A = At',
        B = Bt',
        C = Ct',
        x0 = randn(n),
        y0 = randn(m),
        f = tanh,
        g = tanh,
        # ...
    )
end

const CACHE = Ref{Any}(nothing)

function setup(;
    n = 2^9,
    m = 2^10,
    niters = 1000,
    ntasks = Threads.nthreads(),
    spin = nothing,
    nbranches = 2,
)
    @debug "BenchFuzzy.setup: n=$n m=$m niters=$niters ntasks=$ntasks spin=$spin nbranches=$nbranches"

    CACHE[] = random_model(n, m)

    suite = BenchmarkGroup()

    suite["seq"] = @benchmarkable sim_seq(CACHE[], $niters)
    suite["blas"] = @benchmarkable sim_seq(CACHE[], $niters, LinearAlgebra.mul!)

    for (label, barrier) in [
        # Fuzzy barriers:
        ("tree", TreeBarrier{nbranches}),
        ("flat-tree", FlatTreeBarrier{nbranches}),
        ("centralized", CentralizedBarrier),
    ]
        s1 = suite[label] = BenchmarkGroup()
        s1["static"] = @benchmarkable sim_parallel_fuzzy(
            CACHE[],
            $niters,
            $barrier($ntasks),
            $barrier($ntasks),
            parallel_foreach_static,
            $spin,
        )
        s1["dynamic"] = @benchmarkable sim_parallel_fuzzy(
            CACHE[],
            $niters,
            $barrier($ntasks),
            $barrier($ntasks),
            parallel_foreach_dynamic,
            $spin,
        )
    end

    for (label, barrier) in [
        # Non-fuzzy barriers:
        ("flat-tree-no-fuzzy", FlatTreeBarrier{nbranches}),
        ("centralized-no-fuzzy", CentralizedBarrier),
        ("dissemination", DisseminationBarrier),
    ]
        s1 = suite[label] = BenchmarkGroup()
        s1["static"] = @benchmarkable sim_parallel(
            CACHE[],
            $niters,
            $barrier($ntasks),
            parallel_foreach_static,
            $spin,
        )
        s1["dynamic"] = @benchmarkable sim_parallel(
            CACHE[],
            $niters,
            $barrier($ntasks),
            parallel_foreach_dynamic,
            $spin,
        )
    end

    let s1 = suite["nobarrier"] = BenchmarkGroup()
        s1["static"] = @benchmarkable sim_parallel_nobarrier(
            CACHE[],
            $niters,
            $ntasks,
            parallel_foreach_static,
        )
        s1["dynamic"] = @benchmarkable sim_parallel_nobarrier(
            CACHE[],
            $niters,
            $ntasks,
            parallel_foreach_static,
        )
    end

    return suite
end

function clear()
    CACHE[] = nothing
end

end  # module


text:
using DataStructures
using StatsBase
using PyPlot
using ResumableFunctions
using SimJulia
using Random, Distributions

# Simulation run time and warm-up
sim_duration = 5000
warm_up = 1000

# Average time between patients arriving
inter_arrival_time = 10

# Number of doctors in ED
number_of_docs = 2

# Time between audits
audit_interval = 100

# Average time patients spends being *treated* in ED
appointment_time_mean = 18
appointment_time_sd = 18

# Priority queue for waiting patients
patient_queue = PriorityQueue()

# Vectors used to store audit results
audit_time = []
audit_patients_in_ED = []
audit_patients_waiting = []
audit_patients_waiting_p1 = []
audit_patients_waiting_p2 = []
audit_patients_waiting_p3 = []
audit_reources_used = []

# Proportion of patients in priorities 1,2,3
priority_dist = [0.3, 0.3, 0.4]

# Set up counter for number of a ptients entering simulation
patient_count = 0

# Set up running counts of patients waiting (total and by priority)
patients_waiting = 0
patients_waiting_by_priority = [0, 0, 0]

# Set up patient type

mutable struct Patient
    #=
    Patient type
    =#

    uid # Unique patient ID
    priority # Priority 1-3 (1=highest)
    consulation_time # Time needed with doc
    queuing_time # Record time spent queing
    time_in # Record time arrival
    time_see_doc # Record time see doc
    time_out # Record exit time
end

function patient_init(env)
    #=
    Create new patient.
    =#

    # Set references to global parameters
    global patient_count
    global priority_dist
    global appointment_time_mean
    global appointment_time_sd

    # Increment patient count
    patient_count += 1
    uid = patient_count

    # Sample priority

    cum_probs = cumsum(priority_dist)
    rand_num = rand()
    if rand_num <= cum_probs[1]
        p_priority = 1
    elseif rand_num <= cum_probs[2]
        p_priority = 2
    else
        p_priority = 3
    end

    #= Sample Sample time (assume log normal distribution). Note: log parameters
    are estimates from mean and sd. In a real sim fit σ & μ directly to data =#
    m = appointment_time_mean
    s = appointment_time_sd

    σ_log = s/m * 0.85
    μ_log = log(m) - σ_log^2/2
    d = LogNormal(μ_log, σ_log)
    consult_time = rand(d)

    # Set inital times (these will be updated in sim)
    queuing_time = 0
    time_in = now(env)
    time_see_doc = 0
    time_out = 0

    patient = Patient(uid, p_priority, consult_time, queuing_time, time_in,
        time_see_doc, time_out)

    return patient

end

# Generate new arrivals

@resumable function new_admission(env::Environment)

    #=
    Continuous loop of new arrivals:
    1) Sample interval until next arrival from exponential distribution
    2) Wait for required interval
    3) Call a new patient process
    =#

    x = 0
    # Set references to global parameters
    global inter_arrival_time

    while true
        # Get time to next patient arrival from exponential distribution
        next_p = rand(Exponential(inter_arrival_time))
        # Wait for next patient arrival
        @yield timeout(env, next_p)
        # Generate patient
        patient = patient_init(env)
        x += 1
        println(x)
        # Call patient pathway process
        # p = @process patient_pathway(env, patient)
    end
end

# Main

# Set up simulation environment
sim = Simulation()
# Set up starting processes (pass environment to process)
@process new_admission(sim)
# Run simulation (pass simulation time)
run(sim, warm_up + sim_duration)

println("end")


text:
"""
    dsblend(din, PARAM, iflag)

forward and adjoint blend operator

# Arguments
- `din`: input data
- `PARAM`: parameters
- `iflag`: 1=forward operator
          -1=adjoint operator

"""
function dsblend(din,PARAM,iflag)
    tt=PARAM[1];
    Nt=PARAM[2];
    Ns=PARAM[3];
    dt=PARAM[4];
    (nt,ns)=size(din);
    tn=floor.(Int,tt/dt) .+1;
    mn=maximum(tn);
    if iflag==1
        out=zeros((nt+mn),1);
        for i=1:ns
            out[tn[i]:(tn[i]+nt-1)]=out[tn[i]:(tn[i]+nt-1)].+din[:,i];
        end
    end
    if iflag==-1
        out=zeros(Nt,Ns);
        for j=1:Ns
            out[1:Nt,j]=out[1:Nt,j]+din[tn[j]:(tn[j]+Nt-1)];
        end
    end
    return out

end

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
"""
HotSideRadiationModel:
Implementation of types to model the incident radiation on the "hot" side of
the scrap bulk.
"""
abstract type HotSideRadiationModel end

"""
SimpleRadiationModel - For a simple radiation model:

dot{Q}_rad = sigma * epsilon * A * (T - T_inf)

"""
struct SimpleRadiationModel <: HotSideRadiationModel
    T_inf::Float64
end

### S2S Model types

"""
S2Ssurface:
Represents a general surface in the s2s model.
"""
struct S2Ssurface{T1 <: Int, T2 <: Real}
    name::String
    id::T1
    T::T2
    area::T2
    epsilon::T2
end

# """
# S2SRadiationModel
# To use a more sophisticated radiation model, where externally calculated view
# factors can be used, to model the influence an envoirement with different hot
# spots in front of the scrap bulk.
# """
mutable struct S2SRadiationModel{T1 <: Integer, T2 <: Real} <: HotSideRadiationModel
    front_surface_id::T1
    n_surfaces::T1
    surfaces::Vector{S2Ssurface{T1, T2}}
    vf_mat::Matrix{T2}
    has_initialized::Bool
    A::Matrix{T2}
    B::Vector{T2}
    # G::Vector{T2}
    # Qrad::Vector{T2}
end

function set_S2SRadiationModel(surf_names::Vector{String},  vf_mat::Matrix{T1},
T_surfaces::Vector{T1}, A_surfaces::Vector{T1}, epsilon_surfaces::Vector{T1},
front_surface_id::T2 = 1) where {T1 <: Real, T2 <: Integer}

    @assert size(vf_mat,1) == size(vf_mat,2)

    n = size(vf_mat,1)

    @assert n == length(surf_names)
    @assert n == length(T_surfaces)
    @assert n == length(A_surfaces)
    @assert n == length(epsilon_surfaces)

    surfaces = Vector{S2Ssurface{T2, T1}}(undef, n)
    for i = 1:n
        surfaces[i] = S2Ssurface(surf_names[i], i, T_surfaces[i], A_surfaces[i],
                epsilon_surfaces[i])
    end

    A = zeros(n, n)
    B = zeros(n)

    return S2SRadiationModel(front_surface_id, n, surfaces, vf_mat, false, A, B) # , G, Qrad)
end

"""
heat_flow_rad_front(rm::S2SRadiationModel,eps::Real, A::Real,T::Real)

heat flow calculation for scrap front cells with s2s modelß0
"""
function heat_flow_front_rm(rm::S2SRadiationModel, eps::Real, b::Real, A::Real, T::Real)::Real

    if !rm.has_initialized
        @inbounds for i1 = 1:rm.n_surfaces, i2 = 1:rm.n_surfaces
            if i1 == i2
                rm.A[i1, i2] = ((rm.surfaces[i1].epsilon - 1) * rm.vf_mat[i1, i2]) + 1
            else
                rm.A[i1, i2] = (rm.surfaces[i1].epsilon - 1) * rm.vf_mat[i1, i2]
            end
        end

        T_faces = zeros(rm.n_surfaces)

        for i = 1:rm.n_surfaces
            T_faces[i]  = rm.surfaces[i].T
        end

        T_faces[rm.front_surface_id] = T

        for ii = 1:rm.n_surfaces
            rm.B[ii] = sigma * rm.surfaces[ii].epsilon * (T_faces[ii]^4);
        end

        rm.has_initialized = true
    else
        rm.B[rm.front_surface_id] = sigma * rm.surfaces[rm.front_surface_id].epsilon * (T^4);
    end

    J = rm.A\rm.B;

    # what do we realy need for Q_other_to_front_surface_id
    # @inbounds for i1 = 1:rm.n_surfaces
    #     G_zw = 0
    #     for i2 = 1:rm.n_surfaces
    #         G_zw = G_zw + J[i2] * rm.vf_mat[i1, i2]
    #     end
    #     rm.G[i1] = G_zw
    # end

    # @inbounds for ii=1:rm.n_surfaces
    #     rm.Qrad[ii] = (J[ii] - rm.G[ii]) * rm.surfaces[ii].area
    # end

    # Qrad = -1 .* rm.Qrad[rm.front_surface_id]

    # G: incident radiation
    # J: outgoining including reflection

    G_front_face = 0.0
    for i2 = 1:rm.n_surfaces
        G_front_face = G_front_face + J[i2] * rm.vf_mat[rm.front_surface_id, i2]
    end

    Qrad = (J[rm.front_surface_id] - G_front_face) * rm.surfaces[rm.front_surface_id].area

    return - b * Qrad
end

"""
heat_flow_rad_front(rm::SimpleRadiationModel,eps::Real, A::Real,T::Real)

heat flow calculation for scrap front cells ...
"""
function heat_flow_front_rm(rm::SimpleRadiationModel, eps::Real, b::Real, A::Real, T::Real)::Real

    return b * hf_radiation_inf(eps, A, rm.T_inf, T)
end


text:
type Network
    g::LightGraphs.DiGraph
    osm_id::Vector{Int}                     # Graph ID (1:n) -> OSM ID
    node_id::Dict{Int,Int}                  # OSM ID -> Graph ID
    way_ids::Dict{Tuple{Int,Int},Vector{Int}}    # (Graph ID) -> osm.ways index
    distmx::SparseMatrixCSC{Float64, Int}   # Edge weights, indexed by Graph ID
end

"fetch all osm ids from the list of osm elements"
function osm_indices{T <: OSMElement}(osm_elements::Vector{T})
    osm_id = Array(Int, length(osm_elements))
    for i in 1:length(osm_elements)
        osm_id[i] = osm_elements[i].id
    end
    osm_id
end

"collect all the osm indices in the list of ways"
function collect_nodes(ways::Vector{Way})
    way_nodes = Set{Int}()
    for w in ways
        for n in w.nodes
            push!(way_nodes, n)
        end
    end
    collect(way_nodes)
end

"construct reverse mapping from osmid -> nodeid"
function node_indices(osm_id::Vector{Int})
    node_id = Dict{Int,Int}()
    sizehint!(node_id, length(osm_id))
    for i in 1:length(osm_id)
        node_id[osm_id[i]] = i
    end
    node_id
end

toradians(degree::Float64) = degree * π / 180.0
todegrees(radian::Float64) = radian * 180.0 / π

"distance between the two points in kilometres"
function distance(pt1::Tuple{Float64, Float64}, pt2::Tuple{Float64, Float64})
    dLat = toradians(pt2[2] - pt1[2])
    dLon = toradians(pt2[1] - pt1[1])
    lat1 = toradians(pt1[2])
    lat2 = toradians(pt2[2])
    a = sin(dLat/2)^2 + sin(dLon/2)^2 * cos(lat1) * cos(lat2)
    2.0 * atan2(sqrt(a), sqrt(1-a)) * 6373.0
end

"""
construct a road network from the osm data,

* According to the corresponding speeds (for each way) in km/hour.
* If there are multiple ways correpsonding to an edge, it assumes
  that the road classes corresponding to those ways are the same
* Works only with visible road_classes[highways] for now.
"""
function osmnetwork(osm::OSMData, road_speed::Dict{Int,Int}=SPEED_ROADS_URBAN)
    # filter down to roads along highways
    ways = filter(highway, osm.ways)
    ways = ways[map(visible, ways) & ~(map(services, ways))]
    road_class = map(roadway, ways)
    ways = ways[~(road_class .== 0)]
    speed = [road_speed[c] for c in road_class[~(road_class .== 0)]]
    @assert length(speed) == length(ways)

    # collect all nodes found along the roads
    nodes = osm.nodes
    numnodes = length(nodes)
    osm_ids = osm_indices(nodes)
    node_id = node_indices(osm_ids)

    # construct the graph
    edges, way_ids = construct_edges(node_id, ways)
    V = Float64[distance(nodes[e[1]].lonlat, nodes[e[2]].lonlat) / speed[way_ids[e][1]]
                for e in edges]
    edges = reinterpret(Int, edges)
    I = edges[1:2:end] # collect all start nodes
    J = edges[2:2:end] # collect all end nodes
    distmx = sparse(I, J, V, numnodes, numnodes)

    Network(LightGraphs.DiGraph(distmx), osm_ids, node_id, way_ids, distmx)
end

function construct_edges(node_id::Dict{Int,Int}, ways::Vector{Way})
    rev_ways = map(reverse, ways)
    directed_ways = map(oneway, ways)
    edges = Set{Tuple{Int,Int}}()
    way_ids = Dict{Tuple{Int,Int},Vector{Int}}()

    for i in 1:length(ways)
        way = ways[i]
        rev, nrev = rev_ways[i], !rev_ways[i]
        for n in 2:length(way.nodes)
            n0 = node_id[way.nodes[n-1]] # map osm_id -> node_id
            n1 = node_id[way.nodes[n]]
            startnode = n0*nrev + n1*rev # reverse the direction if need be
            endnode = n0*rev + n1*nrev

            edge = (startnode, endnode)
            push!(edges, edge)
            !haskey(way_ids, edge) && (way_ids[edge] = Vector{Int}())
            push!(way_ids[edge], i)

            if !directed_ways[i]
                edge = (endnode, startnode)
                push!(edges, edge)
                !haskey(way_ids, edge) && (way_ids[edge] = Vector{Int}())
                push!(way_ids[edge], i)
            end
        end
    end
    collect(edges), way_ids
end


text:
"""
    get_sha(obj, directory::AbstractString = "data", file::AbstractString = "states_current.csv")::String

Navigate the entries in the tree to find the SHA1 of the file at path (directory/file)
"""
function get_sha(obj, directory::AbstractString = "data", file::AbstractString = "states_current.csv")::String
    node = obj.node.tree.entries
    node_ = node[findfirst(elem -> elem.name == directory, node)].object.entries
    node_[findfirst(elem -> elem.name == file, node_)].oid
end
"""
    find_shas(obj::GitHubPersonalAccessToken,
              id::AbstractString = "MDEwOlJlcG9zaXRvcnkyNDY0MTE2MDc=",
              directory::AbstractString = "data",
              file::AbstractString = "states_current.csv",
              just_last::Bool = false)

Return slug for the repository and SHA1 for the file at directory/file (all versions).
"""
function find_shas(obj::GitHubPersonalAccessToken,
                   id::AbstractString = "MDEwOlJlcG9zaXRvcnkyNDY0MTE2MDc=",
                   directory::AbstractString = "data",
                   file::AbstractString = "states_current.csv",
                   just_last::Bool = false)
    # For testing
    # obj = opt.pat
    # id = "MDEwOlJlcG9zaXRvcnkyNDY0MTE2MDc="
    # directory = "data"
    # file = "states_current.csv"
    # just_last = false
    # The node ID is the repository with slug COVID19Tracking/covid-tracking-data
    # The file that has the information we want is at: data/states_current.csv

    # uva-bi-sdad/COVID_Tracking_Clean
    # id = "MDEwOlJlcG9zaXRvcnkyNTIwNDUwOTQ="
    # directory = "data"
    # file = "daily.tsv"
    # just_last = true
    vars = Dict("id" => id,
                "path" => "$directory/$file",
                # When first started publishing the quality checks
                "since" => "2020-03-20T23:00:07Z",
                "cursor" => nothing,
                "first" => just_last ? 1 : 100)
    # The initial request finds out how many commits are for the file
    # and the SHA1 for the last 100 commits affecting that file
    response = graphql(obj, "Magic", vars, GITHUB_API_QUERY = GITHUB_API_QUERY)
    json = JSON3.read(response.Data)
    # Get the slug just in case
    slug = json.data.node.nameWithOwner
    # We store the total count of commits for a final check
    total = json.data.node.defaultBranchRef.target.history.totalCount
    # We get the SHA1 for each of the file/versions
    shas = get_sha.(json.data.node.defaultBranchRef.target.history.edges, directory, file)
    just_last && return shas[1]
    # This strategy is valid while the total number of commits is below 1,000
    @assert total ≤ 1_000 "Code needs to be updated for more than 1,000 commits!"
    # If there are more than 100 commits we paginate
    while json.data.node.defaultBranchRef.target.history.pageInfo.hasNextPage
        response = graphql(obj,
                           "Magic",
                           merge(vars, Dict("cursor" => json.data.node.defaultBranchRef.target.history.pageInfo.endCursor)),
                           GITHUB_API_QUERY = GITHUB_API_QUERY)
        json = JSON3.read(response.Data)
        append!(shas, get_sha.(json.data.node.defaultBranchRef.target.history.edges, directory, file))
    end
    # Verify the total number of SHA1 matches the number of commits
    @assert slug == "COVID19Tracking/covid-tracking-data" "Repository has been renamed or moved!"
    @assert total == length(shas) "Some commit SHA1 were not collected!"
    shas
end
"""
    get_tbl(sha::AbstractString)

Uses the SHA1 to download from GitHub the version of [COVID19Tracking/covid-tracking-data/data/states_current.csv](https://github.com/COVID19Tracking/covid-tracking-data/blob/master/data/states_current.csv)
"""
function get_tbl(obj::GitHubPersonalAccessToken, sha::AbstractString)
    response = restful(obj, "repos/COVID19Tracking/covid-tracking-data/git/blobs/$sha")
    @assert response.status == 200
    json = JSON3.read(response.body)
    data = File(base64decode(json.content),
                # select = [:state, :positiveScore, :negativeScore, :negativeRegularScore, :commercialScore, :checkTimeEt, :dataQualityGrade]) |>
    ) |>
        DataFrame
        names(data)
    colnames = names(data)
    if :dataQualityGrade in colnames
        data = dropmissing!(data[!,[:state, :checkTimeEt, :dataQualityGrade]])
        data[!,:checkTimeEt] = ZonedDateTime.(string.("2020/", data.checkTimeEt, " America/New_York"),
                                              COVID_TRACKING_DT)
    elseif all(elem -> elem ∈ colnames, (:positiveScore, :negativeScore, :negativeRegularScore, :commercialScore))
        data = dropmissing!(data[!,[:state, :checkTimeEt, :positiveScore, :negativeScore, :negativeRegularScore, :commercialScore]])
        data[!,:positiveScore] = isone.(data.positiveScore)
        data[!,:negativeScore] = isone.(data.negativeScore)
        data[!,:negativeRegularScore] = isone.(data.negativeRegularScore)
        data[!,:commercialScore] = isone.(data.commercialScore)
        data[!,:checkTimeEt] = ZonedDateTime.(string.("2020/", data.checkTimeEt, " America/New_York"),
                                              COVID_TRACKING_DT)
        data[!,:dataQualityGrade] = get.(Ref(Dict(1 => "D", 2 => "C", 3 => "B", 4 => "A")),
                              data.positiveScore + data.negativeScore + data.negativeRegularScore + data.commercialScore,
                              missing)
        data = data[!,[:state, :checkTimeEt, :dataQualityGrade]]
    else
        data = DataFrame()
    end
    data
end
"""
    states_daily()

Uses the Covid Tracking API to return the states/daily table.
"""
function states_daily()
    response = request("GET", "https://covidtracking.com/api/states/daily.csv")
    @assert response.status == 200
    data = File(response.body) |>
        (tbl -> select(tbl, :state, :dateChecked,
                            :positive, :negative, :pending,
                            :hospitalizedCurrently, :hospitalizedCumulative,
                            :inIcuCurrently, :inIcuCumulative,
                            :onVentilatorCurrently, :onVentilatorCumulative,
                            :recovered, :death)) |>
        DataFrame
    data[!,:dateChecked] = ZonedDateTime.(data.dateChecked, DateFormat("yyyy-mm-ddTHH:MM:SSz"))
    rename!(data, :dateChecked => :checkts,
                  :hospitalizedCurrently => :hospitalized_currently,
                  :hospitalizedCumulative => :hospitalized_cumulative,
                  :inIcuCurrently => :icu_currently,
                  :inIcuCumulative => :icu_cumulative,
                  :onVentilatorCurrently => :ventilator_currently,
                  :onVentilatorCumulative => :ventilator_cumulative,
                  )
end


text:
const obsnames = ["Energy", "Energy^2", "Specific Heat",
                  "|Magnetization|", "Magnetization^2", "Magnetization^4",
                  "Susceptibility", "Connected Susceptibility", "Binder Ratio",
                 ]
const Ns = [2,8]
const J = 1.0
const Ts = [0.3, 1.0, 3.0, 10.0]

"""
Exact calculation for Ising model on `N` site fully connected network
"""
function exact(J, N, Ts)
    j = J/N
    Es = zeros(N+1)
    Ms = zeros(N+1)
    ns = zeros(Int,N+1)
    for Nup in 0:N
        Ndown = N - Nup
        n = binomial(N, Nup)
        push!(ns, n)
        E = j*(Nup*Ndown - 0.5(Nup*(Nup-1)) - 0.5(Ndown*(Ndown-1)))
        push!(Es, E)
        M = (Nup-Ndown)/N
        push!(Ms, M)
    end
    sorted = sortperm(Es, rev=true)
    E0 = Es[sorted[end]]

    mbetas = Ts .\ (-1.0)
    nT = length(Ts)
    res = Dict(n=>zeros(nT) for n in obsnames)
    res["Magnetization"] = zeros(nT)
    Zs = zeros(nT)
    for i in sorted
        E = Es[i] - E0
        zs = exp.(mbetas.*E) * ns[i]
        Zs .+= zs
        E += E0
        res["Energy"] .+= zs.*(E/N)
        res["Energy^2"] .+= zs.*(E/N)^2

        M = Ms[i]
        res["Magnetization"] .+= zs.*M
        res["|Magnetization|"] .+= zs.*abs(M)
        res["Magnetization^2"] .+= zs.*M^2
        res["Magnetization^4"] .+= zs.*M^4
    end
    res["Magnetization"] ./= Zs
    for n in obsnames
        res[n] ./= Zs
    end
    return res
end

betas = 1.0 ./ Ts

for N in Ns
    res = exact(J,N,Ts)
    res["Specific Heat"] = N.*(res["Energy^2"] .- res["Energy"].^2).*betas.^2
    res["Susceptibility"] = N.*(res["Magnetization^2"] .- res["Magnetization"].^2).*betas
    res["Connected Susceptibility"] = N.*(res["Magnetization^2"] .- res["|Magnetization|"].^2).*betas
    res["Binder Ratio"] = res["Magnetization^4"] ./ (res["Magnetization^2"].^2)
    open(@sprintf("J_%.1f__N_%d.dat", J, N), "w") do io
        for i in 1:length(Ts)
            @printf(io, "%.15f", Ts[i])
            for n in obsnames
                @printf(io, " %.15f", res[n][i])
            end
            println(io)
        end
    end
end


text:
module TestDataStream
    using Base.Test
    using DataArrays
    using DataFrames

    path = Pkg.dir("DataFrames",
                   "test",
                   "data",
                   "separators",
                   "sample_data.csv")

    ds = readstream(path, nrows = 1)
    n = start(ds)
    (df, n) = next(ds, n)
    @assert done(ds, n) == false
    (df, n) = next(ds, n)
    @assert done(ds, n) == false
    (df, n) = next(ds, n)
    @assert done(ds, n) == true

    path = Pkg.dir("DataFrames",
                   "test",
                   "data",
                   "scaling",
                   "10000rows.csv")

    ds = readstream(path, nrows = 100)
    n = start(ds)
    (df, n) = next(ds, n)
    @assert done(ds, n) == false

    ds = readstream(path, nrows = 5)

    for minibatch in ds
        @assert isa(minibatch, DataFrame)
    end

    # means = colmeans(ds)
    # @assert abs(means[1, 4] - (-0.005686449)) < 10e-4
    # @assert abs(means[1, 5] - 19.01197) < 10e-4

    # vars = colvars(ds)
    # @assert abs(vars[1, 4] - 0.98048) < 10e-4
    # @assert abs(vars[1, 5] - 0.989416) < 10e-4

    # (mins, maxs) = colranges(ds)
    # @assert abs(mins[1, 4] - (-4.33635)) < 10e-4
    # @assert abs(mins[1, 5] - 15.6219) < 10e-4
    # @assert abs(maxs[1, 4] - 3.86857) < 10e-4
    # @assert abs(maxs[1, 5] - 22.574) < 10e-4

    # covariances = cov(ds)
    # for i in 1:ncol(covariances)
    #     if isna(covariances[i, i])
    #         @assert isna(vars[1, i])
    #     else
    #         @assert abs(covariances[i, i] - vars[1, i]) < 10e-4
    #     end
    # end
    # @assert abs(covariances[4, 4] - 0.980479916) < 10e-4
    # @assert abs(covariances[4, 5] - 0.009823644) < 10e-4
    # @assert abs(covariances[5, 4] - 0.009823644) < 10e-4
    # @assert abs(covariances[5, 5] - 0.989415811) < 10e-4

    # correlations = cor(ds)
    # for i in ncol(correlations)
    #     for j in ncol(correlations)
    #         true_value = covariances[i, j] / sqrt(covariances[i, i] * covariances[j, j])
    #         if isna(true_value)
    #             @assert isna(correlations[i, j])
    #         else
    #             @assert abs(correlations[i, j] - true_value) < 10e-4
    #         end
    #     end
    # end

    # Deal with different delimiters
    path = Pkg.dir("DataFrames",
                   "test",
                   "data",
                   "separators",
                   "sample_data.csv")
    ds = readstream(path)
    for row in ds
        @assert size(row, 1) <= 1
    end

    path = Pkg.dir("DataFrames",
                   "test",
                   "data",
                   "separators",
                   "sample_data.tsv")
    ds = readstream(path)
    for row in ds
        @assert size(row, 1) <= 1
    end

    path = Pkg.dir("DataFrames",
                   "test",
                   "data",
                   "separators",
                   "sample_data.wsv")
    ds = readstream(path)
    for row in ds
        @assert size(row, 1) <= 1
    end

    # # DataFrame to DataStream conversion
    # df = DataFrame(A = 1:25)

    # ds = DataStream(df, 5)

    # for mini in ds
    #     @assert nrow(mini) <= 5
    # end

    # colmeans(ds)
end


text:
int_rules_1_1_1_6 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.1.6*P(x)*(a+b*x)^m*(c+d*x)^n*(e+f*x)^p =#
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => Antiderivative(~Px * (~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m * (~e + ~f * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ(~m, ~n) && (IntegerQ(~m) || GtQ(~a, 0) && GtQ(~c, 0)))))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~n * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => (((~a + ~b * ~x) ^ FracPart(~m) * (~c + ~d * ~x) ^ FracPart(~m)) / (~a * ~c + ~b * ~d * (~x) ^ 2) ^ FracPart(~m)) * Antiderivative(~Px * (~a * ~c + ~b * ~d * (~x) ^ 2) ^ ~m * (~e + ~f * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && (EqQ(~b * ~c + ~a * ~d, 0) && (EqQ(~m, ~n) && Not(IntegerQ(~m)))))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => Antiderivative(PolynomialQuotient(~Px, ~a + ~b * ~x, ~x) * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && EqQ(PolynomialRemainder(~Px, ~a + ~b * ~x, ~x), 0))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand(~Px * (~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && IntegersQ(~m, ~n))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => With([Qx = PolynomialQuotient(~Px, ~a + ~b * ~x, ~x), R = PolynomialRemainder(~Px, ~a + ~b * ~x, ~x)], (~b * R * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * ExpandToSum((((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * Qx + ~a * ~d * ~f * R * (~m + 1)) - ~b * R * (~d * ~e * (~m + ~n + 2) + ~c * ~f * (~m + ~p + 2))) - ~b * ~d * ~f * R * (~m + ~n + ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && ILtQ(~m, -1))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~m * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => With([Qx = PolynomialQuotient(~Px, ~a + ~b * ~x, ~x), R = PolynomialRemainder(~Px, ~a + ~b * ~x, ~x)], (~b * R * (~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f)) + (1 / ((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f))) * Antiderivative((~a + ~b * ~x) ^ (~m + 1) * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * ExpandToSum((((~m + 1) * (~b * ~c - ~a * ~d) * (~b * ~e - ~a * ~f) * Qx + ~a * ~d * ~f * R * (~m + 1)) - ~b * R * (~d * ~e * (~m + ~n + 2) + ~c * ~f * (~m + ~p + 2))) - ~b * ~d * ~f * R * (~m + ~n + ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~n, ~p], ~x) && (PolyQ(~Px, ~x) && (LtQ(~m, -1) && IntegersQ(2 * ~m, 2 * ~n, 2 * ~p)))
    @apply_utils Antiderivative(~Px * (~(a') + ~(b') * ~x) ^ ~(m') * (~(c') + ~(d') * ~x) ^ ~(n') * (~(e') + ~(f') * ~x) ^ ~(p'), ~x) => With([q = Expon(~Px, ~x), k = Coeff(~Px, ~x, Expon(~Px, ~x))], (k * (~a + ~b * ~x) ^ ((~m + q) - 1) * (~c + ~d * ~x) ^ (~n + 1) * (~e + ~f * ~x) ^ (~p + 1)) / (~d * ~f * (~b) ^ (q - 1) * (~m + ~n + ~p + q + 1)) + (1 / (~d * ~f * (~b) ^ q * (~m + ~n + ~p + q + 1))) * Antiderivative((~a + ~b * ~x) ^ ~m * (~c + ~d * ~x) ^ ~n * (~e + ~f * ~x) ^ ~p * ExpandToSum((~d * ~f * (~b) ^ q * (~m + ~n + ~p + q + 1) * ~Px - ~d * ~f * k * (~m + ~n + ~p + q + 1) * (~a + ~b * ~x) ^ q) + k * (~a + ~b * ~x) ^ (q - 2) * (((~a) ^ 2 * ~d * ~f * (~m + ~n + ~p + q + 1) - ~b * (~b * ~c * ~e * ((~m + q) - 1) + ~a * (~d * ~e * (~n + 1) + ~c * ~f * (~p + 1)))) + ~b * (~a * ~d * ~f * (2 * (~m + q) + ~n + ~p) - ~b * (~d * ~e * (~m + q + ~n) + ~c * ~f * (~m + q + ~p))) * ~x), ~x), ~x) <-- NeQ(~m + ~n + ~p + q + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~m, ~n, ~p], ~x) && PolyQ(~Px, ~x)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.1 Linear/1.1.1.6 P(x) (a+b x)^m (c+d x)^n (e+f x)^p.jl


text:
"""
    reduce(op, ss::AbstractArray{<:ShiftedVector}, args...; default = missing, filter = t->true, dropmissing = true)

Align all vectors in `ss`. For each of the indices in `args`, extract the collection of
elements corresponding to that index and reduce using `op`. `op` can be any binary reduction
function. Indices for which the iterable is empty will return `default`.

# Examples

```jldoctest reduce
julia> v = [1, 3, 5, 9, 6, 7];

julia> ss = ShiftedArray.((v,), [-2, -4])
2-element Array{ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}},1}:
 Union{Int64, Missing}[5, 9, 6, 7, missing, missing]
 Union{Int64, Missing}[6, 7, missing, missing, missing, missing]

julia> reduce(+, ss, 1:2)
2-element Array{Int64,1}:
 11
 16

julia> reduce(+, ss, -5:2)
8-element Array{Any,1}:
   missing
   missing
  1
  3
  6
 12
 11
 16
```
"""
function reduce end

"""
    mapreduce(g, op, ss::AbstractArray{<:ShiftedVector}, args...; default = missing, filter = t->true, dropmissing = true)

Align all vectors in `ss`. For each of the indices in `args`, extract the collection of
elements corresponding to that index, apply `g` and reduce using `op`. `op` can be any binary reduction
function. Indices for which the iterable is empty will return `default`.

# Examples

```jldoctest mapreduce
julia> v = [1, 3, 5, 9, 6, 7];

julia> ss = ShiftedArray.((v,), [-2, -4])
2-element Array{ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}},1}:
 Union{Int64, Missing}[5, 9, 6, 7, missing, missing]
 Union{Int64, Missing}[6, 7, missing, missing, missing, missing]

julia> mapreduce(t -> t^2, +, ss, 1:2)
2-element Array{Int64,1}:
 61
 130

julia> mapreduce(t -> t^2, +, ss, -5:2)
8-element Array{Any,1}:
    missing
    missing
   1
   9
  26
  90
  61
 130
```
"""
function mapreduce end

"""
    reduce_vec(f, ss::AbstractArray{<:ShiftedVector}, args...; default = missing, filter = t->true, dropmissing = true)

Align all vectors in `ss`. For each of the indices in `args`, extract the collection of
elements corresponding to that index and reduce using `f`. `f` is any function that
takes an iterable as input and outputs a scalar, such as `mean`. Indices for which the
iterable is empty will return `default`.

# Examples

```jldoctest reduce_vec
julia> v = [1, 3, 5, 9, 6, 7];

julia> ss = ShiftedArray.((v,), [-2, -4])
2-element Array{ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}},1}:
 Union{Int64, Missing}[5, 9, 6, 7, missing, missing]
 Union{Int64, Missing}[6, 7, missing, missing, missing, missing]

julia> reduce_vec(mean, ss, 1:2)
2-element Array{Float64,1}:
 5.5
 8.0

julia> reduce_vec(mean, ss, -5:2)
8-element Array{Any,1}:
  missing
  missing
 1.0
 3.0
 3.0
 6.0
 5.5
 8.0
```
"""
function reduce_vec end

"""
    mapreduce_vec(g, f, ss::AbstractArray{<:ShiftedVector}, args...; default = missing, filter = t->true, dropmissing = true)

Align all vectors in `ss`. For each of the indices in `args`, extract the collection of
elements corresponding to that index, apply `g` and reduce using `f`. `f` is any function that
takes an iterable as input and outputs a scalar, such as `mean`. Indices for which the
iterable is empty will return `default`.

# Examples

```jldoctest mapreduce_vec
julia> v = [1, 3, 5, 9, 6, 7];

julia> ss = ShiftedArray.((v,), [-2, -4])
2-element Array{ShiftedArrays.ShiftedArray{Int64,Missing,1,Array{Int64,1}},1}:
 Union{Int64, Missing}[5, 9, 6, 7, missing, missing]
 Union{Int64, Missing}[6, 7, missing, missing, missing, missing]

julia> mapreduce_vec(log, mean, ss, 1:2)
2-element Array{Float64,1}:
 1.7006
 2.07157

julia> mapreduce_vec(log, mean, ss, -5:2)
8-element Array{Any,1}:
  missing
  missing
 0.0
 1.09861
 0.804719
 1.64792
 1.7006
 2.07157
```
"""
function mapreduce_vec end

mapreduce_vec(g, f, itr) = f(g(x) for x in itr)

for (_mapreduce, mapreduce, reduce) in [(:_mapreduce, :mapreduce, :reduce), (:_mapreduce_vec, :mapreduce_vec, :reduce_vec)]
    @eval begin
        # auxiliary method to filter away missings and return a default if no values are found
        function ($_mapreduce)(g, op, itr; default = missing, filter = t->true, dropmissing = true)
            missingfree_itr = dropmissing ? skipmissing(itr) : itr
            filtered_itr = Iterators.filter(filter, missingfree_itr)
            isempty(filtered_itr) ? default : ($mapreduce)(g, op, filtered_itr)
        end

        # align the shifted arrays and apply _mapreduce
        function ($mapreduce)(g, f, ss::AbstractArray{<:ShiftedArray{<:Any, <:Any, N}}, args::Vararg{<:AbstractArray, N}; kwargs...) where{N}
            inds = Base.product(args...)
            [($_mapreduce)(g, f, (s[CartesianIndex(i)] for s in ss); kwargs...) for i in inds]
        end

        # define corresponding reduce methods
        ($reduce)(op, ss::AbstractArray{<:ShiftedArray{<:Any, <:Any, N}}, args::Vararg{<:AbstractArray, N}; kwargs...) where{N} =
            ($mapreduce)(identity, op, ss, args...; kwargs...)
    end
end


text:
using CSV

function stan_generate_quantities(m::SampleModel, id::Int64; kwargs...)

  local fname

  cmd = ``
  if isa(m, SampleModel)
    # Handle the model name field for unix and windows
    cmd = `$(m.exec_path)`

    # Sample() specific portion of the model
    cmd = `$cmd generate_quantities`

    # Fitted_params is required
    fname = "$(m.output_base)_chain_$id.csv"
    cmd = `$cmd fitted_params=$fname`

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file[id])
      fname = m.data_file[id]
      cmd = `$cmd data file=$fname`
    end

    fname = "$(m.output_base)_generated_quantities_$id.csv"
    cmd = `$cmd output file=$fname`
  end

  cd(m.tmpdir) do
    run(pipeline(cmd, stdout="$(m.output_base)_generated_quantities_$id.log"))
  end

  CSV.read(fname, delim=",", comment="#")

end


text:
export POP!, PUSH!, COPYPUSH!, COPYPOP!, GLOBAL_STACK

const GLOBAL_STACK = []

############# global stack operations ##########
@inline function PUSH!(x)
    push!(GLOBAL_STACK, x)
    cleared(x)
end

@inline function POP!(x::T) where T
    @invcheck x cleared(x)
    loaddata(T, pop!(GLOBAL_STACK))
end

UNSAFE_PUSH!(args...) = PUSH!(args...)
@inline function UNSAFE_POP!(x::T) where T
    loaddata(T, pop!(GLOBAL_STACK))
end

############# local stack operations ##########
@inline function PUSH!(stack, x)
    push!(stack, x)
    stack, cleared(x)
end

@inline function POP!(stack, x::T) where T
    @invcheck x cleared(x)
    stack, loaddata(T, pop!(stack))
end

@inline function UNSAFE_POP!(stack, x::T) where T
    stack, loaddata(T, pop!(stack))
end

cleared(x) = zero(x)
cleared(x::AbstractArray{T,N}) where {T,N} = similar(x, zeros(Int, N)...)
cleared(x::AbstractVector{T}) where {T} = empty(x)

loaddata(::Type{T}, x::T) where T = x

@dual POP! PUSH!
@dual UNSAFE_POP! UNSAFE_PUSH!

############# copied push/pop stack operations ##########
@inline function COPYPUSH!(stack, x)
    push!(stack, copy(x))
    stack, x
end

@inline function COPYPOP!(stack, x)
    y = pop!(stack)
    @invcheck x y
    stack, x
end

UNSAFE_COPYPUSH!(args...) = COPYPUSH!(args...)
@inline function UNSAFE_COPYPOP!(stack, x)
    y = pop!(stack)
    stack, y
end

@inline function COPYPUSH!(x)
    push!(GLOBAL_STACK, x)
    x
end

@inline function COPYPOP!(x)
    y = pop!(GLOBAL_STACK)
    @invcheck x y
    x
end

@inline function UNSAFE_COPYPOP!(x)
    y = pop!(GLOBAL_STACK)
    y
end

@dual COPYPOP! COPYPUSH!
@dual UNSAFE_COPYPOP! UNSAFE_COPYPUSH!

loaddata(::Type{T}, x::TX) where {T<:NullType, TX} = T(x)


text:
MLT_PATH = joinpath(splitpath(@__DIR__)[1:end-1]...)
push!(LOAD_PATH, MLT_PATH)

using Documenter, MLToolkit

makedocs(;
    sitename = "MLToolkit Documentation",
    pages = [
        "Home" => "index.md",
        "Modules" => [
            "plots.md",
        ]
    ]
)

deploydocs(
    repo = "https://github.com/xukai92/MLToolkit.jl.git",
)


text:
##############################################################################################
module TestIterativeSmootherExps
##############################################################################################
##############################################################################################
# imports and exports
using DataAssimilationBenchmarks.SmootherExps
using JLD2, Statistics

##############################################################################################
##############################################################################################
# run and analyze the IEnKS for state estimation with the Lorenz-96 model
# arguments are
# time_series, method, seed, nanl, lag, shift, mda, obs_un, obs_dim, γ, N_ens, infl = args

function run_sda_smoother_state_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "ienks-transform", 0, 3500, 10, 1, false,
                1.0, 40, 1.00, 21, 1.02)
        iterative_state(args)
        true
    catch
        false
    end
end

function analyze_sda_smoother_state_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/ienks-transform/")
        data = load(path * "ienks-transform_L96_state_seed_0000_diff_0.000_sysD_40_obsD_40_" *
                    "obsU_1.00_gamma_001.0_nanl_03500_tanl_0.05_h_0.05_lag_010_shift_001_" *
                    "mda_false_nens_021_stateInfl_1.02.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if mean(filt_rmse[501:end]) < 0.2 && mean(post_rmse[501:end]) < 0.10
            true
        else
            false
        end
    catch
        false
    end
end

function run_mda_smoother_state_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "ienks-transform", 0, 3500, 10, 1, true,
                1.0, 40, 1.00, 21, 1.02)
        iterative_state(args)
        true
    catch
        false
    end
end

function analyze_mda_smoother_state_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/ienks-transform/")
        data = load(path * "ienks-transform_L96_state_seed_0000_diff_0.000_sysD_40_obsD_40_" *
                    "obsU_1.00_gamma_001.0_nanl_03500_tanl_0.05_h_0.05_lag_010_shift_001_" *
                    "mda_true_nens_021_stateInfl_1.02.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if mean(filt_rmse[501:end]) < 0.2 && mean(post_rmse[501:end]) < 0.10
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# run and analyze the IEnKS for joint state-parameter estimation with the Lorenz-96 model
# time_series, method, seed, nanl, lag, shift, mda, obs_un, obs_dim, γ,
# param_err, param_wlk, N_ens, state_infl, param_infl = args

function run_sda_smoother_param_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "ienks-transform", 0, 3500, 10, 1, false, 1.0, 40, 1.0,
                0.10, 0.0010, 21, 1.02, 1.0)
        iterative_param(args)
        true
    catch
        false
    end
end

function analyze_sda_smoother_param_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/ienks-transform/")
        data = load(path * "ienks-transform_L96_param_seed_0000_diff_0.000_sysD_41_obsD_40_" *
                    "obsU_1.00_gamma_001.0_paramE_0.10_paramW_0.0010_nanl_03500_tanl_0.05_" *
                    "h_0.05_lag_010_shift_001_mda_false_nens_021_stateInfl_1.02_" *
                    "paramInfl_1.00.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]
        para_rmse = data["param_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if (mean(filt_rmse[501:end]) < 0.2) && (mean(post_rmse[501:end]) < 0.10) &&
            (mean(para_rmse[501:end]) < 0.01)
            true
        else
            false
        end
    catch
        false
    end
end

function run_mda_smoother_param_L96()
    try
        path = joinpath(@__DIR__, "../src/data/time_series/")
        L96_test_ts = "L96_time_series_seed_0000_dim_40_diff_0.000_F_08.0" *
                      "_tanl_0.05_nanl_05000_spin_1500_h_0.050.jld2"

        args = (path * L96_test_ts, "ienks-transform", 0, 3500, 10, 1, true, 1.0, 40, 1.0,
                0.10, 0.0010, 21, 1.02, 1.0)
        iterative_param(args)
        true
    catch
        false
    end
end

function analyze_mda_smoother_param_L96()
    try
        # test if the filter RMSE for standard simulation falls below adequate threshold
        path = joinpath(@__DIR__, "../src/data/ienks-transform/")
        data = load(path * "ienks-transform_L96_param_seed_0000_diff_0.000_sysD_41_obsD_40_" *
                    "obsU_1.00_gamma_001.0_paramE_0.10_paramW_0.0010_nanl_03500_tanl_0.05_" *
                    "h_0.05_lag_010_shift_001_mda_true_nens_021_stateInfl_1.02_" *
                    "paramInfl_1.00.jld2")
        filt_rmse = data["filt_rmse"]
        post_rmse = data["post_rmse"]
        para_rmse = data["param_rmse"]

        # note, we use a small burn-in to reach more regular cycles
        if (mean(filt_rmse[501:end]) < 0.2) && (mean(post_rmse[501:end]) < 0.10) &&
            (mean(para_rmse[501:end]) < 0.01)
            true
        else
            false
        end
    catch
        false
    end
end

##############################################################################################
# end module

end

max_stars_repo_name: cgrudz/DataAssimilationBenchmarks

max_issues_repo_name: cgrudz/DataAssimilationBenchmarks


text:
using DataFrames, Dates, Parameters

# Abstract type for dispatch
abstract type TreatmentPanel end

# Types for number of treatment units and periods
abstract type UnitTreatmentType end
abstract type TreatmentDurationType end
struct SingleUnitTreatment <: UnitTreatmentType end
struct MultiUnitSimultaneousTreatment <: UnitTreatmentType end
struct MultiUnitStaggeredTreatment <: UnitTreatmentType end
struct ContinuousTreatment <: TreatmentDurationType end
struct StartEndTreatment <: TreatmentDurationType end

# BalancedPanel will have an N×T matrix of treatment assigment and outcomes
"""
    BalancedPanel{UnitTreatmentType, TreatmentDurationType}

A TreatmentPanel in which all N treatment units are observed for the same T periods.

The object is constructed from a `DataFrame` which contains all outcome and covariate data. The
constructor requires passing information on which column in the `DataFrame` holds the subject and
time period identifiers, as well as information on the timing of treatment. Treatments are specified
as a `Pair` of either a `String` or `Symbol` identifying the unit treated, and a `Date` or
`Int` value, identifying the treatment period. Where treatments have an end point, the
treatment period os specified as a `Tuple` of either `Date` or `Int`, indicating the first
and last period of treatment. Where individual units have multiple treatment periods,
these are specified as a `Pair{Union{String, Symbol}, Vector{Union{Int, Date}}}` of treatment
unit identifier and vector of timings. Finally, single continuous, single periodic, and
multiple period treatments can be generalized to multiple treated units.
The following table provides an overview of the types of treatment pattern supported:

|                 |  Only starting point        |   Start and end point                     |   Multiple start & end points                       |
|-----------      |------------------------     |-------------------------                  |------------------------------                       |
| **one unit**         |  Pair{String, Date}         |   Pair{String, Tuple{Date, Date}}         |  Pair{String}, Vector{Tuple{Date, Date}}}           |
| **multiple units**  |  Vector{Pair{String, Date}} |   Vector{Pair{String, Tuple{Date, Date}}} |  Vector{Pair{String}, Vector{Tuple{Date, Date}}}}   |

Currently, only single treatment unit and continuous treatment is supported.
"""
@with_kw struct BalancedPanel{UTType, TDType} <: TreatmentPanel where UTType <: UnitTreatmentType where TDType <: TreatmentDurationType
    N::Int64
    T::Int64
    W::Union{Matrix{Bool}, Matrix{Union{Missing, Bool}}}
    ts::Vector{T1} where T1 <: Union{Date, Int64}
    is::Vector{T2} where T2 <: Union{Symbol, String, Int64}
    Y::Matrix{Float64}
end

# Constructor based on DatFrame and treatment assignment in pairs
function BalancedPanel(df::DataFrame, treatment_assignment::Vector{Pair{NType, TType}};
    id_var = nothing,
    t_var = nothing,
    outcome_var = nothing,
    sort_inplace = false) where NType where TType <: Union{Date, Int64}

    ### SANITY CHECKS ###

    # Check relevant info has been provided
    !isnothing(outcome_var) || error(ArgumentError(
        "Please specify outcome_var, the name of the column in your dataset holding the "*
        "outcome variable of interest in your dataset."
        ))
    !isnothing(id_var) || error(ArgumentError(
            "Please specify id_var, the name of the column in your data set holding the "*
            "identifier of your units of observation (panel dimension)."
            ))
    !isnothing(t_var) || error(ArgumentError(
            "Please specify t_var, the name of the column in your dataset holding the "*
            "time dimension."
        ))
    !isnothing(treatment_assignment) || error(ArgumentError(
        "Please specify treatment assignment, the identifier of the treated unit(s) in your "*
        "dataset and associated start date(s) of treatment."
        ))

    # Ensure columns exist in data
    f = in(names(df))
    f(string(id_var)) || throw("Error: ID variable $id_var is not present in the data.")
    f(string(t_var)) || throw("Error: Time variable $t_var is not present in the data.")
    f(string(outcome_var)) || throw("Error: ID variable $outcome_var is not present in the data.")

    # Get all units and time periods
    is = sort(unique(df[!, id_var])); i_set = Set(is)
    ts = sort(unique(df[!, t_var])); t_set = Set(ts)

    for tp ∈ treatment_assignment
        in(tp[1], i_set) || throw("Error: Treatment assignment $tp provided, but $(tp[1]) is not in the list of unit identifiers $id_var")
        in(tp[2], t_set) || throw("Error: Treatment assignment $tp provided, but $(tp[2]) is not in the list of time identifiers $t_var")
    end

    # Sort data if necessary, in place if required
    df = ifelse(issorted(df, [id_var, t_var]), df,
                                               ifelse(sort_inplace, sort!(df, [id_var, t_var]),
                                                                    sort(df, [id_var, t_var])))

    # Dimensions
    N = length(is)
    T = length(ts)

    # Treatment matrix
    W = [false for i ∈ 1:N, t ∈ 1:T]

    for tp ∈ treatment_assignment
        i_id = findfirst(==(tp[1]), is)
        t_id = findfirst(==(tp[2]), ts)
        W[i_id, t_id:end] .= true
    end

    # Outcome matrix
    Y = zeros(size(W))

    for (row, i) ∈ enumerate(is), (col, t) ∈ enumerate(ts)
        try
            Y[row, col] = only(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), outcome_var])
        catch ArgumentError
            throw("$(nrow(df[(df[!, id_var] .== i) .& (df[!, t_var] .== t), :])) outcomes present in the data for unit $i in period $t")
        end
    end

    # Determine UnitTreatmentType and TreatmentDurationType
    uttype = if length(treatment_assignment) == 1
        SingleUnitTreatment
    else
        if all(==(treatment_assignment[1][2]), last.(treatment_assignment))
            MultiUnitSimultaneousTreatment
        else
            MultiUnitStaggeredTreatment
        end
    end

    tdtype = ContinuousTreatment

    BalancedPanel{uttype, tdtype}(N, T, W, ts, is, Y)
end

# Constructor for single treatment
function BalancedPanel(df::DataFrame, treatment_assignment::Pair{NType, TType};
    id_var = nothing, t_var= nothing, outcome_var = nothing, sort_inplace = false) where NType where TType

    BalancedPanel(df, [treatment_assignment]; id_var = id_var,
        t_var = t_var, outcome_var = outcome_var, sort_inplace = sort_inplace)

end

# UnblancedPanel - N observations but not all of them for T periods

#!# Not yet implemented

# Utility functions
function treated_ids(x::BalancedPanel{SingleUnitTreatment, T}) where T
    for i ∈ 1:x.N
        for t ∈ 1:x.T
            if x.W[i, t]
                return i
            end
        end
    end
end

function treated_labels(x::BalancedPanel{SingleUnitTreatment, T}) where T
    x.is[treated_ids(x)]
end

function first_treated_period_ids(x::BalancedPanel{SingleUnitTreatment, T}) where T
    findfirst(x.W[treated_ids(x), :])
end

function first_treated_period_labels(x::BalancedPanel{SingleUnitTreatment, T}) where T
    x.ts[first_treated_period_ids(x)]
end

function length_T₀(x::BalancedPanel{SingleUnitTreatment, T}) where T
    first_treated_period_ids(x) - 1
end

function length_T₁(x::BalancedPanel{SingleUnitTreatment, T}) where T
    x.T - first_treated_period_ids(x) + 1
end


text:
using ProgrammingProjects
using Documenter

DocMeta.setdocmeta!(ProgrammingProjects, :DocTestSetup, :(using ProgrammingProjects); recursive=true)

makedocs(;
    modules=[ProgrammingProjects],
    authors="Gil Junqueira",
    repo="https://github.com/gjunqueira-sys/ProgrammingProjects.jl/blob/{commit}{path}#{line}",
    sitename="ProgrammingProjects.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://gjunqueira-sys.github.io/ProgrammingProjects.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/gjunqueira-sys/ProgrammingProjects.jl",
    devbranch="main",
)


text:
# Test 2D modeling
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

parsed_args = parse_commandline()

# Set parallel if specified
nw = parsed_args["parallel"]
if nw > 1 && nworkers() < nw
   addprocs(nw-nworkers() + 1; exeflags=["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
end

@everywhere using JOLI
@everywhere using JUDI, LinearAlgebra, Test, Distributed

### Model
model, model0, dm = setup_model(parsed_args["tti"], parsed_args["nlayer"]; n=(101, 101), d=(10., 10.))
q, srcGeometry, recGeometry, info = setup_geom(model; nsrc=2, tn=500f0)
dt = srcGeometry.dt[1]

# Modeling operators
println("Generic modeling and misc test with ", parsed_args["nlayer"], " layers and tti: ", parsed_args["tti"])

ftol = 1f-5
######################## WITH DENSITY ############################################

cases = [(true, true), (true, false), (false, true), (false, false)]

@testset "Generic tests with limit_m = $(limit_m)  and save to disk = $(disk)" for (limit_m, disk)=cases
	# Options structures
	opt = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0,
					file_path=pwd(),	# path to files
					file_name="shot_record")	# saves files as file_name_xsrc_ysrc.segy

	opt0 = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0,
					file_path=pwd(),	# path to files
					file_name="smooth_shot_record")	# saves files as file_name_xsrc_ysrc.segy

	optJ = Options(save_data_to_disk=disk, limit_m=limit_m, buffer_size=100f0,
					file_path=pwd(),	# path to files
					file_name="linearized_shot_record")	# saves files as file_name_xsrc_ysrc.segy

	# Setup operators
	Pr = judiProjection(info, recGeometry)
	F = judiModeling(info, model; options=opt)
	F0 = judiModeling(info, model0; options=opt0)
	Ps = judiProjection(info, srcGeometry)

	# Combined operator Pr*F*adjoint(Ps)
	Ffull = judiModeling(model, srcGeometry, recGeometry)
	J = judiJacobian(Pr*F0*adjoint(Ps),q; options=optJ) # equivalent to J = judiJacobian(Ffull,q)

	# Nonlinear modeling
	d1 = Pr*F*adjoint(Ps)*q	# equivalent to d = Ffull*q
	dfull = Ffull*q
	@test isapprox(get_data(d1), dfull, rtol=ftol)

	qad = (Ps*adjoint(F))*adjoint(Pr)*d1
	qfull = adjoint(Ffull)*d1
	@test isapprox(qad, qfull, rtol=ftol)

	# fwi objective function
	f, g = fwi_objective(model0, q, d1; options=opt)
	f, g = fwi_objective(model0, subsample(q,1), subsample(d1,1); options=opt)

	# Subsampling
	for inds=[2, [1, 2]]
		dsub = subsample(dfull, inds)
		qsub = subsample(q, inds)
			Fsub = subsample(F, inds)
		Jsub = subsample(J, inds)
		Ffullsub = subsample(Ffull, inds)
		Pssub = subsample(Ps, inds)
		Prsub = subsample(Pr, inds)
		ds1 = Ffullsub*qsub
		ds2 = Prsub * Fsub * adjoint(Pssub) *qsub
		@test isapprox(ds1, get_data(ds2), rtol=ftol)
		@test isapprox(ds1, dsub, rtol=ftol)
		@test isapprox(get_data(ds2), dsub, rtol=ftol)
	end

	# vcat, norms, dot
	dcat = [d1, d1]
	@test isapprox(norm(d1)^2, .5f0*norm(dcat)^2)
	@test isapprox(dot(d1, d1), norm(d1)^2)
end

############################# Full wavefield ############################################

@testset "Basic judiWavefield modeling tests" begin
	opt = Options(dt_comp=dt)
	F = judiModeling(info, model; options=opt)
	Fa = judiModelingAdjoint(info, model; options=opt)
	Ps = judiProjection(info, srcGeometry)
	Pr = judiProjection(info, recGeometry)

	# Return wavefields
	u = F * adjoint(Ps) * q

	# Adjoint from data
	dobs = Pr*F*u
	v = Fa*(adjoint(Pr)*dobs)

	a = dot(u, v)
	b = dot(dobs, dobs)
	@printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", b, a, a/b - 1)
	@test isapprox(a, b, rtol=1f-4)

	# Wavefields as source + return wavefields
	u2 = F*u
	v2 = Fa*v

	a = dot(u2, v)
	b = dot(v2, u)
	@printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", a, b, a/b - 1)
	@test isapprox(a, b, rtol=1f-4)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
function init_local_gaussian(center, dx, dy, dz, H)
    X = Data.Array([
        2 * exp(
            -1.0 * (
                (x_g(ix, dx, H) + dx / 2 - center[1])^2 +
                (y_g(iy, dy, H) + dy / 2 - center[2])^2 +
                (z_g(iz, dz, H) + dz / 2 - center[3])^2
            ),
        ) for ix = 1:size(H, 1), iy = 1:size(H, 2), iz = 1:size(H, 3)
    ])
    return X
end

function apply_boundary_conditions!(H, coords, dims)
    if coords[1] == 1
        H[1, :, :] .= 0.0
    end
    if coords[2] == 1
        H[:, 1, :] .= 0.0
    end
    if coords[3] == 1
        H[:, :, 1] .= 0.0
    end

    if coords[1] == dims[1]
        H[end, :, :] .= 0.0
    end
    if coords[2] == dims[2]
        H[:, end, :] .= 0.0
    end
    if coords[3] == dims[3]
        H[:, :, end] .= 0.0
    end
end

function dist_norm_L2(Rh, comm_cart)
    sq_residual = [sum(Rh .^ 2)]  # this is local
    MPI.Allreduce!(sq_residual, +, comm_cart)
    return sqrt(sq_residual[1])
end

function linear_interpolate_3D(H, dx)
    dy = dz = dx

    ix, iy, iz = Int.((LOCATION_OF_INTEREST.÷ dx) .+ 1)

    # manual trilinear interpolation
    # see https://en.wikipedia.org/wiki/Trilinear_interpolation#Alternative_algorithm
    x0, x1 = ix*dx+dx/2, (ix+1)*dx+dx/2
    y0, y1 = iy*dy+dy/2, (iy+1)*dy+dy/2
    z0, z1 = iz*dz+dz/2, (iz+1)*dz+dz/2

    cvec = [H[ix, iy, iz], H[ix+1, iy, iz], H[ix, iy+1, iz], H[ix+1, iy+1, iz], H[ix, iy, iz+1], H[ix+1, iy, iz+1], H[ix, iy+1, iz+1], H[ix+1, iy+1, iz+1]]

    xvec = repeat([x0;x1], 4)
    yvec = repeat([y0; y0; y1; y1], 2)
    zvec = repeat(repeat([z0], 4), 2)

    M = hcat(ones(8), xvec, yvec, zvec, xvec.*yvec, xvec.*zvec, yvec.*zvec, xvec.*yvec.*zvec)

    interp_value = try
        avec = M \ cvec
        x, y, z = LOCATION_OF_INTEREST
        xvec_ = [1 x y z x*y x*z y*z x*z*z]
        interp_value = xvec_ * avec
        interp_value
    catch e
        interp_value = H[ix, iy, iz]
        interp_value
    end
end


text:
##
using Plots
using FFTW
using DSP

##
# Number of points
N = 2^14 - 1
# Sample period
Ts = 1 / (1.1 * N)
# Start time
t0 = 0
tmax = t0 + N * Ts
# time coordinate
t = t0:Ts:tmax

# signal
signal = sin.(2π * 60 .* t) # sin (2π f t)

# Fourier Transform of it
F = fft(signal) |> fftshift
w = fftfreq(length(t), 1.0/Ts) |> fftshift

# plots
time_domain = plot(t, signal, title = "Signal")
freq_domain = plot(w, abs.(F), title = "Spectrum", xlim=(-1000, +1000))
plot(time_domain, freq_domain, layout = 2)

##
f(x,y) = x + y

function g(x,y)
    x * y
end


text:
#
# supportedAxes(::AbstractBackend) = [:left]
# supportedTypes(::AbstractBackend) = []
# supportedStyles(::AbstractBackend) = [:solid]
# supportedMarkers(::AbstractBackend) = [:none]
# supportedScales(::AbstractBackend) = [:identity]
# subplotSupported(::AbstractBackend) = false
# stringsSupported(::AbstractBackend) = false
#
# supportedAxes() = supportedAxes(backend())
# supportedTypes() = supportedTypes(backend())
# supportedStyles() = supportedStyles(backend())
# supportedMarkers() = supportedMarkers(backend())
# supportedScales() = supportedScales(backend())
# subplotSupported() = subplotSupported(backend())
# stringsSupported() = stringsSupported(backend())

# --------------------------------------------------------------------------------------
#
#
# supportedArgs(::GadflyBackend) = [
#     :annotation,
#     :background_color, :foreground_color, :color_palette,
#     :group,
#     :label,
#     :linetype,
#     :seriescolor, :seriesalpha,
#     :linecolor, :linestyle, :linewidth, :linealpha,
#     :markershape, :markercolor, :markersize, :markeralpha,
#     :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
#     :fillrange, :fillcolor, :fillalpha,
#     :bins,
#     :n, :nc, :nr, :layout,
#     :smooth,
#     :title, :windowtitle, :show, :size,
#     :x, :xlabel, :xlims, :xticks, :xscale, :xflip,
#     :y, :ylabel, :ylims, :yticks, :yscale, :yflip,
#     # :z, :zlabel, :zlims, :zticks, :zscale, :zflip,
#     :z,
#     :tickfont, :guidefont, :legendfont,
#     :grid, :legend, :colorbar,
#     :marker_z, :levels,
#     :xerror, :yerror,
#     :ribbon, :quiver,
#     :orientation,
#   ]
# supportedAxes(::GadflyBackend) = [:auto, :left]
# supportedTypes(::GadflyBackend) = [
#         :none, :line, :path, :steppre, :steppost, :sticks,
#         :scatter, :hist2d, :hexbin, :hist,
#         :bar, :box, :violin, :quiver,
#         :hline, :vline, :contour, :shape
#     ]
# supportedStyles(::GadflyBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
# supportedMarkers(::GadflyBackend) = vcat(_allMarkers, Shape)
# supportedScales(::GadflyBackend) = [:identity, :ln, :log2, :log10, :asinh, :sqrt]
# subplotSupported(::GadflyBackend) = true
#
#
# # --------------------------------------------------------------------------------------

#
# supportedArgs(::ImmerseBackend) = supportedArgs(GadflyBackend())
# supportedAxes(::ImmerseBackend) = supportedAxes(GadflyBackend())
# supportedTypes(::ImmerseBackend) = supportedTypes(GadflyBackend())
# supportedStyles(::ImmerseBackend) = supportedStyles(GadflyBackend())
# supportedMarkers(::ImmerseBackend) = supportedMarkers(GadflyBackend())
# supportedScales(::ImmerseBackend) = supportedScales(GadflyBackend())
# subplotSupported(::ImmerseBackend) = true
#
# # --------------------------------------------------------------------------------------

#
#
# supportedArgs(::PyPlotBackend) = [
#     :annotation,
#     :background_color, :foreground_color, :color_palette,
#     :background_color_legend, :background_color_inside, :background_color_outside,
#     :foreground_color_legend, :foreground_color_grid, :foreground_color_axis,
#         :foreground_color_text, :foreground_color_border,
#     :group,
#     :label,
#     :linetype,
#     :seriescolor, :seriesalpha,
#     :linecolor, :linestyle, :linewidth, :linealpha,
#     :markershape, :markercolor, :markersize, :markeralpha,
#     :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
#     :fillrange, :fillcolor, :fillalpha,
#     :bins,
#     :n, :nc, :nr, :layout,
#     :smooth,
#     :title, :windowtitle, :show, :size,
#     :x, :xlabel, :xlims, :xticks, :xscale, :xflip,
#     :y, :ylabel, :ylims, :yticks, :yscale, :yflip,
#     :axis, :yrightlabel,
#     :z, :zlabel, :zlims, :zticks, :zscale, :zflip,
#     :z,
#     :tickfont, :guidefont, :legendfont,
#     :grid, :legend, :colorbar,
#     :marker_z, :levels,
#     :xerror, :yerror,
#     :ribbon, :quiver,
#     :orientation,
#     :overwrite_figure,
#     :polar,
#     :normalize, :weights, :contours, :aspect_ratio
#   ]
# supportedAxes(::PyPlotBackend) = _allAxes
# supportedTypes(::PyPlotBackend) = [
#         :none, :line, :path, :steppre, :steppost, :shape,
#         :scatter, :hist2d, :hexbin, :hist, :density,
#         :bar, :sticks, :box, :violin, :quiver,
#         :hline, :vline, :heatmap, :pie,
#         :contour, :contour3d, :path3d, :scatter3d, :surface, :wireframe
#     ]
# supportedStyles(::PyPlotBackend) = [:auto, :solid, :dash, :dot, :dashdot]
# supportedMarkers(::PyPlotBackend) = vcat(_allMarkers, Shape)
# supportedScales(::PyPlotBackend) = [:identity, :ln, :log2, :log10]
# subplotSupported(::PyPlotBackend) = true
#
#
# # --------------------------------------------------------------------------------------

#
#
# supportedArgs(::GRBackend) = [
#     :annotation,
#     :background_color, :foreground_color, :color_palette,
#     :group,
#     :label,
#     :linetype,
#     :seriescolor, :seriesalpha,
#     :linecolor, :linestyle, :linewidth, :linealpha,
#     :markershape, :markercolor, :markersize, :markeralpha,
#     :markerstrokewidth, :markerstrokecolor, :markerstrokealpha,
#     :fillrange, :fillcolor, :fillalpha,
#     :bins,
#     :n, :nc, :nr, :layout,
#     :smooth,
#     :title, :windowtitle, :show, :size,
#     :x, :xlabel, :xlims, :xticks, :xscale, :xflip,
#     :y, :ylabel, :ylims, :yticks, :yscale, :yflip,
#     :axis, :yrightlabel,
#     :z, :zlabel, :zlims, :zticks, :zscale, :zflip,
#     :z,
#     :tickfont, :guidefont, :legendfont,
#     :grid, :legend, :colorbar,
#     :marker_z, :levels,
#     :xerror, :yerror,
#     :ribbon, :quiver,
#     :orientation,
#     :overwrite_figure,
#     :polar,
#   ]
# supportedAxes(::GRBackend) = _allAxes
# supportedTypes(::GRBackend) = [:none, :line, :path, :steppre, :steppost, :sticks,
#                                :scatter, :hist2d, :hexbin, :hist, :density, :bar,
#                                :hline, :vline, :contour, :heatmap, :path3d, :scatter3d, :surface,
#                                :wireframe, :ohlc, :pie]
# supportedStyles(::GRBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
# supportedMarkers(::GRBackend) = vcat(_allMarkers, Shape)
# supportedScales(::GRBackend) = [:identity, :log10]
# subplotSupported(::GRBackend) = true
#
#
# # --------------------------------------------------------------------------------------

#
#
# supportedArgs(::QwtBackend) = [
#     :annotation,
#     :axis,
#     :background_color,
#     :linecolor,
#     :color_palette,
#     :fillrange,
#     :fillcolor,
#     :foreground_color,
#     :group,
#     :label,
#     :layout,
#     :legend,
#     :seriescolor, :seriesalpha,
#     :linestyle,
#     :linetype,
#     :linewidth,
#     :markershape,
#     :markercolor,
#     :markersize,
#     :n,
#     :bins,
#     :nc,
#     :nr,
#     :pos,
#     :smooth,
#     :show,
#     :size,
#     :title,
#     :windowtitle,
#     :x,
#     :xlabel,
#     :xlims,
#     :xticks,
#     :y,
#     :ylabel,
#     :ylims,
#     :yrightlabel,
#     :yticks,
#     :xscale,
#     :yscale,
#   ]
# supportedTypes(::QwtBackend) = [:none, :line, :path, :steppre, :steppost, :sticks, :scatter, :hist2d, :hexbin, :hist, :bar, :hline, :vline]
# supportedMarkers(::QwtBackend) = [:none, :auto, :rect, :ellipse, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5, :star8, :hexagon]
# supportedScales(::QwtBackend) = [:identity, :log10]
# subplotSupported(::QwtBackend) = true
#
#
# # --------------------------------------------------------------------------------------

#
# supportedArgs(::UnicodePlotsBackend) = [
#     # :annotation,
#     # :args,
#     # :axis,
#     # :background_color,
#     # :linecolor,
#     # :fill,
#     # :foreground_color,
#     :group,
#     # :heatmap_c,
#     # :kwargs,
#     :label,
#     # :layout,
#     :legend,
#     :seriescolor, :seriesalpha,
#     :linestyle,
#     :linetype,
#     # :linewidth,
#     :markershape,
#     # :markercolor,
#     # :markersize,
#     # :markerstrokewidth,
#     # :markerstrokecolor,
#     # :markerstrokestyle,
#     # :n,
#     :bins,
#     # :nc,
#     # :nr,
#     # :pos,
#     # :reg,
#     # :ribbon,
#     :show,
#     :size,
#     :title,
#     :windowtitle,
#     :x,
#     :xlabel,
#     :xlims,
#     # :xticks,
#     :y,
#     :ylabel,
#     :ylims,
#     # :yrightlabel,
#     # :yticks,
#     # :xscale,
#     # :yscale,
#     # :xflip,
#     # :yflip,
#     # :z,
#   ]
# supportedAxes(::UnicodePlotsBackend) = [:auto, :left]
# supportedTypes(::UnicodePlotsBackend) = [:none, :line, :path, :steppre, :steppost, :sticks, :scatter, :hist2d, :hexbin, :hist, :bar, :hline, :vline]
# supportedStyles(::UnicodePlotsBackend) = [:auto, :solid]
# supportedMarkers(::UnicodePlotsBackend) = [:none, :auto, :ellipse]
# supportedScales(::UnicodePlotsBackend) = [:identity]
# subplotSupported(::UnicodePlotsBackend) = true
#
#
#
#
# # --------------------------------------------------------------------------------------

#
# supportedArgs(::WinstonBackend) = [
#     :annotation,
#     # :args,
#     # :axis,
#     # :background_color,
#     :linecolor,
#     :color_palette,
#     :fillrange,
#     :fillcolor,
#     # :foreground_color,
#     :group,
#     # :heatmap_c,
#     # :kwargs,
#     :label,
#     # :layout,
#     :legend,
#     :seriescolor, :seriesalpha,
#     :linestyle,
#     :linetype,
#     :linewidth,
#     :markershape,
#     :markercolor,
#     :markersize,
#     # :markerstrokewidth,
#     # :markerstrokecolor,
#     # :markerstrokestyle,
#     # :n,
#     :bins,
#     # :nc,
#     # :nr,
#     # :pos,
#     :smooth,
#     # :ribbon,
#     :show,
#     :size,
#     :title,
#     :windowtitle,
#     :x,
#     :xlabel,
#     :xlims,
#     # :xticks,
#     :y,
#     :ylabel,
#     :ylims,
#     # :yrightlabel,
#     # :yticks,
#     :xscale,
#     :yscale,
#     # :xflip,
#     # :yflip,
#     # :z,
#   ]
# supportedAxes(::WinstonBackend) = [:auto, :left]
# supportedTypes(::WinstonBackend) = [:none, :line, :path, :sticks, :scatter, :hist, :bar]
# supportedStyles(::WinstonBackend) = [:auto, :solid, :dash, :dot, :dashdot]
# supportedMarkers(::WinstonBackend) = [:none, :auto, :rect, :ellipse, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5]
# supportedScales(::WinstonBackend) = [:identity, :log10]
# subplotSupported(::WinstonBackend) = false
#
#
# # --------------------------------------------------------------------------------------

#
#
# supportedArgs(::BokehBackend) = [
#     # :annotation,
#     # :axis,
#     # :background_color,
#     :linecolor,
#     # :color_palette,
#     # :fillrange,
#     # :fillcolor,
#     # :fillalpha,
#     # :foreground_color,
#     :group,
#     # :label,
#     # :layout,
#     # :legend,
#     :seriescolor, :seriesalpha,
#     :linestyle,
#     :linetype,
#     :linewidth,
#     # :linealpha,
#     :markershape,
#     :markercolor,
#     :markersize,
#     # :markeralpha,
#     # :markerstrokewidth,
#     # :markerstrokecolor,
#     # :markerstrokestyle,
#     # :n,
#     # :bins,
#     # :nc,
#     # :nr,
#     # :pos,
#     # :smooth,
#     # :show,
#     :size,
#     :title,
#     # :windowtitle,
#     :x,
#     # :xlabel,
#     # :xlims,
#     # :xticks,
#     :y,
#     # :ylabel,
#     # :ylims,
#     # :yrightlabel,
#     # :yticks,
#     # :xscale,
#     # :yscale,
#     # :xflip,
#     # :yflip,
#     # :z,
#     # :tickfont,
#     # :guidefont,
#     # :legendfont,
#     # :grid,
#     # :surface,
#     # :levels,
#   ]
# supportedAxes(::BokehBackend) = [:auto, :left]
# supportedTypes(::BokehBackend) = [:none, :path, :scatter] #,:steppre, :steppost, :sticks, :hist2d, :hexbin, :hist, :bar, :hline, :vline, :contour]
# supportedStyles(::BokehBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
# supportedMarkers(::BokehBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5] #vcat(_allMarkers, Shape)
# supportedScales(::BokehBackend) = [:identity, :ln] #, :ln, :log2, :log10, :asinh, :sqrt]
# subplotSupported(::BokehBackend) = false
#
#
# # --------------------------------------------------------------------------------------
#
# supportedArgs(::PlotlyBackend) = [
#     :annotation,
#     # :axis,
#     :background_color,
#     :color_palette,
#     :fillrange,
#     :fillcolor,
#     :fillalpha,
#     :foreground_color,
#     :group,
#     :label,
#     :layout,
#     :legend,
#     :seriescolor, :seriesalpha,
#     :linecolor,
#     :linestyle,
#     :linetype,
#     :linewidth,
#     :linealpha,
#     :markershape,
#     :markercolor,
#     :markersize,
#     :markeralpha,
#     :markerstrokewidth,
#     :markerstrokecolor,
#     :markerstrokestyle,
#     :n,
#     :bins,
#     :nc,
#     :nr,
#     # :pos,
#     # :smooth,
#     :show,
#     :size,
#     :title,
#     :windowtitle,
#     :x,
#     :xlabel,
#     :xlims,
#     :xticks,
#     :y,
#     :ylabel,
#     :ylims,
#     # :yrightlabel,
#     :yticks,
#     :xscale,
#     :yscale,
#     :xflip,
#     :yflip,
#     :z,
#     :marker_z,
#     :tickfont,
#     :guidefont,
#     :legendfont,
#     :grid,
#     :levels,
#     :xerror,
#     :yerror,
#     :ribbon,
#     :quiver,
#     :orientation,
#     :polar,
#   ]
# supportedAxes(::PlotlyBackend) = [:auto, :left]
# supportedTypes(::PlotlyBackend) = [:none, :line, :path, :scatter, :steppre, :steppost,
#                                    :hist2d, :hist, :density, :bar, :contour, :surface, :path3d, :scatter3d,
#                                    :pie, :heatmap] #,, :sticks, :hexbin, :hline, :vline]
# supportedStyles(::PlotlyBackend) = [:auto, :solid, :dash, :dot, :dashdot]
# supportedMarkers(::PlotlyBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross,
#                                      :pentagon, :hexagon, :octagon, :vline, :hline] #vcat(_allMarkers, Shape)
# supportedScales(::PlotlyBackend) = [:identity, :log10] #, :ln, :log2, :log10, :asinh, :sqrt]
# subplotSupported(::PlotlyBackend) = true
# stringsSupported(::PlotlyBackend) = true
#
#
# # --------------------------------------------------------------------------------------
#
# supportedArgs(::PlotlyJSBackend) = [
#     :annotation,
#     # :axis,
#     :background_color,
#     :color_palette,
#     :fillrange,
#     :fillcolor,
#     :fillalpha,
#     :foreground_color,
#     :group,
#     :label,
#     :layout,
#     :legend,
#     :seriescolor, :seriesalpha,
#     :linecolor,
#     :linestyle,
#     :linetype,
#     :linewidth,
#     :linealpha,
#     :markershape,
#     :markercolor,
#     :markersize,
#     :markeralpha,
#     :markerstrokewidth,
#     :markerstrokecolor,
#     :markerstrokestyle,
#     :n,
#     :bins,
#     :nc,
#     :nr,
#     # :pos,
#     # :smooth,
#     :show,
#     :size,
#     :title,
#     :windowtitle,
#     :x,
#     :xlabel,
#     :xlims,
#     :xticks,
#     :y,
#     :ylabel,
#     :ylims,
#     # :yrightlabel,
#     :yticks,
#     :xscale,
#     :yscale,
#     :xflip,
#     :yflip,
#     :z,
#     :marker_z,
#     :tickfont,
#     :guidefont,
#     :legendfont,
#     :grid,
#     :levels,
#     :xerror,
#     :yerror,
#     :ribbon,
#     :quiver,
#     :orientation,
#     :polar,
#   ]
# supportedAxes(::PlotlyJSBackend) = [:auto, :left]
# supportedTypes(::PlotlyJSBackend) = [:none, :line, :path, :scatter, :steppre, :steppost,
#                                    :hist2d, :hist, :density, :bar, :contour, :surface, :path3d, :scatter3d,
#                                    :pie, :heatmap] #,, :sticks, :hexbin, :hline, :vline]
# supportedStyles(::PlotlyJSBackend) = [:auto, :solid, :dash, :dot, :dashdot]
# supportedMarkers(::PlotlyJSBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross,
#                                      :pentagon, :hexagon, :octagon, :vline, :hline] #vcat(_allMarkers, Shape)
# supportedScales(::PlotlyJSBackend) = [:identity, :log10] #, :ln, :log2, :log10, :asinh, :sqrt]
# subplotSupported(::PlotlyJSBackend) = true
# stringsSupported(::PlotlyJSBackend) = true
#
# # --------------------------------------------------------------------------------------
#
# supportedArgs(::GLVisualizeBackend) = [
#     # :annotation,
#     # :axis,
#     # :background_color,
#     # :color_palette,
#     # :fillrange,
#     # :fillcolor,
#     # :fillalpha,
#     # :foreground_color,
#     # :group,
#     # :label,
#     # :layout,
#     # :legend,
#     # :linecolor,
#     # :linestyle,
#      :linetype
#     #  :seriescolor, :seriesalpha,
#     # :linewidth,
#     # :linealpha,
#     # :markershape,
#     # :markercolor,
#     # :markersize,
#     # :markeralpha,
#     # :markerstrokewidth,
#     # :markerstrokecolor,
#     # :markerstrokestyle,
#     # :n,
#     # :bins,
#     # :nc,
#     # :nr,
#     # :pos,
#     # :smooth,
#     # :show,
#     # :size,
#     # :title,
#     # :windowtitle,
#     # :x,
#     # :xlabel,
#     # :xlims,
#     # :xticks,
#     # :y,
#     # :ylabel,
#     # :ylims,
#     # :yrightlabel,
#     # :yticks,
#     # :xscale,
#     # :yscale,
#     # :xflip,
#     # :yflip,
#     # :z,
#     # :tickfont,
#     # :guidefont,
#     # :legendfont,
#     # :grid,
#     # :surface
#     # :levels,
#   ]
# supportedAxes(::GLVisualizeBackend) = [:auto, :left]
# supportedTypes(::GLVisualizeBackend) = [:surface] #, :path, :scatter ,:steppre, :steppost, :sticks, :heatmap, :hexbin, :hist, :bar, :hline, :vline, :contour]
# supportedStyles(::GLVisualizeBackend) = [:auto, :solid] #, :dash, :dot, :dashdot, :dashdotdot]
# supportedMarkers(::GLVisualizeBackend) = [:none, :auto, :ellipse] #, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5] #vcat(_allMarkers, Shape)
# supportedScales(::GLVisualizeBackend) = [:identity] #, :log, :log2, :log10, :asinh, :sqrt]
# subplotSupported(::GLVisualizeBackend) = false
#
# # --------------------------------------------------------------------------------------
#
# supportedArgs(::PGFPlotsBackend) = [
#     # :annotation,
#     # :axis,
#      :background_color,
#     # :color_palette,
#     # :fillrange,
#      :fillcolor,
#      :fillalpha,
#     # :foreground_color,
#     # :group,
#     # :label,
#     # :layout,
#     # :legend,
#      :seriescolor, :seriesalpha,
#      :linecolor,
#      :linestyle,
#      :linetype,
#      :linewidth,
#      :linealpha,
#      :markershape,
#      :markercolor,
#      :markersize,
#      :markeralpha,
#     # :markerstrokewidth,
#      :markerstrokecolor,
#      :markerstrokestyle,
#     # :n,
#     # :bins,
#     # :nc,
#     # :nr,
#     # :pos,
#     # :smooth,
#     # :show,
#     # :size,
#      :title,
#     # :windowtitle,
#      :x,
#      :xlabel,
#      :xlims,
#     # :xticks,
#      :y,
#      :ylabel,
#      :ylims,
#     # :yrightlabel,
#     # :yticks,
#      :xscale,
#      :yscale,
#      :xflip,
#      :yflip,
#      :z,
#      :zscale,
#     # :tickfont,
#     # :guidefont,
#     # :legendfont,
#      :grid,
#     # :surface
#     # :levels,
#   ]
# supportedAxes(::PGFPlotsBackend) = [:auto, :left]
# supportedTypes(::PGFPlotsBackend) = [:path, :path3d, :scatter, :line, :steppre, :stepmid, :steppost, :hist, :bar, :hist2d, :sticks, :ysticks, :xsticks, :contour] #  :hexbin, :hline, :vline,]
# supportedStyles(::PGFPlotsBackend) = [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
# supportedMarkers(::PGFPlotsBackend) = [:none, :auto, :ellipse, :rect, :diamond, :utriangle, :dtriangle, :cross, :xcross, :star5, :pentagon] #vcat(_allMarkers, Shape)
# supportedScales(::PGFPlotsBackend) = [:identity, :log, :ln, :log2, :log10] # :asinh, :sqrt]
# subplotSupported(::PGFPlotsBackend) = false


text:
using EchogramProcessing
using Test

a = rand(100, 100)

medianfilter(a, 9, 9)

meanfilter(a, 9, 9)

max_stars_repo_name: EchoJulia/EchogramProcessing

max_issues_repo_name: EchoJulia/EchogramProcessing

max_forks_repo_name: EchoJulia/EchogramProcessing