text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Test of the function `_pager_key_processing`.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Key processing" begin
    pagerd = _create_pagerd("")

    # Down
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :down)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 11
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 0
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Shift down
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :down, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 15
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 2
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 12
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    # Up
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :up)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 9
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 1
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 1
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Shift up
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :up, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 5
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 2
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 0
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 1
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    # Right
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :right)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 11
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 0
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Right
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :right)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 11
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 0
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Alt right
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :right, alt = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 30
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 0
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Shift right
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :right, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 20
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 2
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 12
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    # Left
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :left)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 9
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 1
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 1
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Alt left
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :left, alt = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 1
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 1
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 1
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Shift left
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :left, shift = true)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 15
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 5
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 4
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 2
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 1
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    # End
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :end)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 30
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 0
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 10
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Home
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :home)

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 1
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 1
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 0
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 1
    @test pagerd.start_col == 10
    @test pagerd.redraw    == false
    @test pagerd.event     == nothing

    # Page up
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :pageup)

    pagerd.redraw          = false
    pagerd.start_row       = 40
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 20
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 10)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 21
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 0
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 1
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    # Page down
    # ==========================================================================

    k = TerminalPager.Keystroke(value = :pagedown)

    pagerd.redraw          = false
    pagerd.start_row       = 20
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 30
    pagerd.columns_cropped = 20
    pagerd.display_size    = (10, 20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 29
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing

    pagerd.redraw          = false
    pagerd.start_row       = 10
    pagerd.start_col       = 10
    pagerd.lines_cropped   = 10
    pagerd.columns_cropped = 20
    pagerd.display_size    = (20,20)

    TerminalPager._pager_key_process!(pagerd, k)

    @test pagerd.start_row == 20
    @test pagerd.start_col == 10
    @test pagerd.redraw    == true
    @test pagerd.event     == nothing
end


text:
"""
Computes the rain and snow advection (down) tendency
"""
function compute_precipitation_advection_tendencies(edmf, grid, state, gm)
    FT = eltype(grid)

    ρ0_c = center_ref_state(state).ρ0
    tendencies_pr = center_tendencies_precipitation(state)
    prog_pr = center_prog_precipitation(state)
    aux_tc = center_aux_turbconv(state)

    # helper to calculate the rain velocity
    # TODO: assuming gm.W = 0
    # TODO: verify translation
    term_vel_rain = aux_tc.term_vel_rain
    term_vel_snow = aux_tc.term_vel_snow

    q_rai = prog_pr.q_rai
    q_sno = prog_pr.q_sno

    If = CCO.DivergenceF2C()
    RB = CCO.RightBiasedC2F(; top = CCO.SetValue(FT(0)))
    ∇ = CCO.DivergenceF2C(; bottom = CCO.Extrapolate())
    wvec = CC.Geometry.WVector

    # TODO - some positivity limiters are needed
    @. aux_tc.qr_tendency_advection = ∇(wvec(RB(ρ0_c * q_rai * term_vel_rain))) / ρ0_c
    @. aux_tc.qs_tendency_advection = ∇(wvec(RB(ρ0_c * q_sno * term_vel_snow))) / ρ0_c

    @. tendencies_pr.q_rai += aux_tc.qr_tendency_advection
    @. tendencies_pr.q_sno += aux_tc.qs_tendency_advection
    return
end

"""
Computes the tendencies to θ_liq_ice, qt and qr due to rain evaporation
"""
function compute_precipitation_sink_tendencies(grid, state, gm, TS::TimeStepping)
    param_set = parameter_set(gm)
    Δt = TS.dt
    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    aux_gm = center_aux_grid_mean(state)
    aux_tc = center_aux_turbconv(state)
    prog_gm = center_prog_grid_mean(state)
    prog_pr = center_prog_precipitation(state)
    tendencies_pr = center_tendencies_precipitation(state)

    @inbounds for k in real_center_indices(grid)
        q_tot_gm = prog_gm.q_tot[k]
        T_gm = aux_gm.T[k]
        # When we fuse loops, this should hopefully disappear
        ts = TD.PhaseEquil_pTq(param_set, p0_c[k], T_gm, q_tot_gm)
        q = TD.PhasePartition(ts)

        # TODO - move limiters elsewhere
        rain_evap =
            -min(
                prog_pr.q_rai[k] / Δt,
                -CM1.evaporation_sublimation(param_set, rain_type, q, prog_pr.q_rai[k], ρ0_c[k], T_gm),
            )

        aux_tc.qr_tendency_evap[k] = rain_evap
        aux_tc.qs_tendency_melt[k] = 0.0
        aux_tc.qs_tendency_dep_sub[k] = 0.0

        aux_tc.qt_tendency_precip_sinks[k] = -rain_evap
        aux_tc.θ_liq_ice_tendency_precip_sinks[k] = θ_liq_ice_helper(ts, -rain_evap)

        tendencies_pr.q_rai[k] += rain_evap
    end
    return
end

ext: jl


text:
import Base: ==, @get!

macro bf(ex)
  @capture(ex, x_ = [w__]) && return :(words[$(Expr(:quote, x))] = @bf [$(esc.(w)...)])
  @capture(ex, [xs__])
  xs = [isexpr(x, :$) ? esc(x.args[1]) :
        isexpr(x, Symbol) ? Expr(:quote, x) :
        @capture(x, [w__]) ? :(Quote(@bf [$(esc.(w)...)])) :
        esc(x) for x in xs]
  :(Word([$(xs...)]))
end

struct Native
  code::String
end

struct Word
  code::Vector{Any}
end

Base.getindex(w::Word, i::AbstractArray) = Word(w.code[i])
Base.endof(w::Word) = endof(w.code)

struct Context
  icode::Vector{Any}
  cache::Dict{Symbol,Any}
end

Context() = Context([], Dict())

literal(ctx, x) = literal(x)
literal(x) = x

const cwords = Dict{Symbol,Any}()

lower(ctx, x) = x

lower_(ctx, w::Word) = Word([lower(ctx, w[1:end-1]).code..., lower(ctx, w.code[end])])

function lower(ctx, w::Word)
  w′ = Word([])
  for i = length(w.code):-1:1
    if haskey(cwords, w.code[i])
      unshift!(w′.code, cwords[w.code[i]](ctx, w[1:i]))
      return w′
    else
      unshift!(w′.code, literal(ctx, w.code[i]))
    end
  end
  return w′
end

const words = Dict{Symbol,Any}()

inline(ctx, x) = x
inline(ctx, w::Symbol) = @get!(ctx.cache, w, inline(ctx, words[w]))
inline(ctx, w::Word) = Word(inline.(ctx, lower(ctx, w).code))

function flatten(w::Word)
  w′ = Word([])
  for w in w.code
    w isa Word ? append!(w′.code, flatten(w).code) : push!(w′.code, w)
  end
  return w′
end

function compile_static(ctx, x)
  c = flatten(inline(ctx, x))
  for x in c.code
    x isa Native || error("Couldn't compile $x")
  end
  join(map(c -> c.code, c.code))
end

@bf inc!   = [Native("+")]
@bf dec!   = [Native("-")]
@bf left!  = [Native("<")]
@bf right! = [Native(">")]
@bf read!  = [Native(",")]
@bf write! = [Native(".")]
@bf debug! = [Native("#")]

struct Quote
  code::Vector{Any}
end

Quote(w::Word) = Quote(w.code)
Word(q::Quote) = Word(q.code)

struct Flip
  code::Any
end

flip(x) = Flip(x)
flip(x::Flip) = x.code
flip(w::Native) = w.code == ">" ? Native("<") : w.code == "<" ? Native(">") : w

lower(ctx, w::Flip) = Flip(lower(ctx, w.code))
inline(ctx, w::Flip) = Word(flip.(flatten(inline(ctx, w.code)).code))

cwords[:while!] = function (ctx, w::Word)
  if length(w.code) >= 2 && w.code[end-1] isa Quote
    lower(ctx, Word([w[1:end-2], Native("["), Word(w.code[end-1].code), Native("]")]))
  else
    lower_(ctx, w)
  end
end

repeated(w, i) = Word([w for _ = 1:i])

step!(n) = repeated(n > 0 ? :right! : :left!, abs(n))

@bf pass = []

@bf stack! = [[step!(2)], while!, step!(2), [step!(2)], while!, step!(-2)]
@bf rstack! = [Flip(:stack!)]

@bf stackswitch! = [dec!, left!,
                    [dec!, left!, rstack!, left!, inc!, right!, stack!, right!], while!,
                    left!, rstack!, left!, left!, inc!]

@bf rpush = [stackswitch!, stack!]
@bf rpop = [rstack!, Flip(:stackswitch!)]

function move!(locs...; mode = :inc!)
  @assert !any(x -> x == 0, locs)
  step = Word([@bf [step!(l), $mode, step!(-l)] for l in locs])
  @bf [[dec!, $step], while!]
end

@bf reset! = [[dec!], while!]

@bf drop = [dec!, left!, reset!, left!]

if!(t, f) = @bf [left!, [right!, $t, dec!], while!, right!,
                 [$f, dec!, right!], while!,
                 left!, inc!]

cwords[:if!] = function (ctx, w::Word)
  length(w.code) ≥ 3 && w.code[end-1] isa Quote && w.code[end-2] isa Quote ||
    return lower_(ctx, w)
  t, f = map(i -> @bf([drop, Word(w.code[end-i]), 0]), (2, 1))
  lower(ctx, @bf [w[1:end-3], if!(t, f), drop])
end

function bytecode(ctx::Context, code)
  n = findfirst(ctx.icode, code)
  n == 0 && (push!(ctx.icode, code); n = length(ctx.icode))
  return n
end

literal(ctx, w::Quote) = literal(bytecode(ctx, Word(w)))

struct Call end

words[:call] = Call()

function partition(w::Word)
  w′ = Word([])
  cur = Word([])
  for w in w.code
    if w == Call()
      push!(cur.code, :rpush)
      push!(w′.code, cur)
      cur = Word([])
    else
      push!(cur.code, w)
    end
  end
  !isempty(cur.code) && push!(w′.code, cur)
  return w′
end

function compile_dynamic(ctx::Context, w::Word)
  bytecode(ctx, w)
  is, i = [], 1
  while i ≤ length(ctx.icode)
    q = ctx.icode[i]
    q = partition(flatten(inline(ctx, q)))
    code = if q isa Word && length(q.code) > 1
      Word(reverse([bytecode(ctx, w) for w in q.code]))
    else
      Flip(@bf [stack!, $q, rstack!])
    end
    push!(is, @bf [1, -, if!(:pass, @bf [drop, $code, 0])])
    i += 1
  end
  interp = @bf [[is..., drop], while!, rstack!]
  compile_static(ctx, @bf [1, stackswitch!, Flip(interp)])
end

isdynamic(x) = false
isdynamic(x::Symbol) = x == :call || (haskey(words, x) && isdynamic(words[x]))
isdynamic(w::Word) = any(isdynamic, w.code)

compile(w::Word) = (isdynamic(w) ? compile_dynamic : compile_static)(Context(), w)

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
using Statistics
using StatsBase:sample
using Base.Threads:@threads
using BenchmarkTools
using Revise
using EvoTrees

n_obs = Int(1e6)
n_vars = 100
n_bins = 255
𝑖 = collect(1:n_obs);
𝑗 = collect(1:n_vars);
X_bin = rand(UInt8, n_obs, n_vars);

function split_set_chunk!(left, right, block, bid, X_bin, feat, cond_bin, offset, chunk_size, lefts, rights, bsizes)
    left_count = 0
    right_count = 0
    @inbounds for i in eachindex(block)
        @inbounds if X_bin[block[i], feat] <= cond_bin
            left_count += 1
            left[offset + chunk_size * (bid - 1) + left_count] = block[i]
        else
            right_count += 1
            right[offset + chunk_size * (bid - 1) + right_count] = block[i]
            # right[offset + chunk_size * (bid - 1) + length(block) - right_count] = block[i]
            # right_count += 1
        end
    end
    lefts[bid] = left_count
    rights[bid] = right_count
    bsizes[bid] = length(block)
    return nothing
end

function split_set_threads!(out, left, right, 𝑖, X_bin::Matrix{S}, feat, cond_bin, offset, chunk_size=2^14) where {S}

    left_count = 0
    right_count = 0
    iter = Iterators.partition(𝑖, chunk_size)
    nblocks = length(iter)
    lefts = zeros(Int, nblocks)
    rights = zeros(Int, nblocks)
    bsizes = zeros(Int, nblocks)

    @sync for (bid, block) in enumerate(iter)
        Threads.@spawn split_set_chunk!(left, right, block, bid, X_bin, feat, cond_bin, offset, chunk_size, lefts, rights, bsizes)
    end

    left_sum = sum(lefts)
    left_cum = 0
    right_cum = 0
    @inbounds for bid in 1:nblocks
        view(out, offset + left_cum + 1:offset + left_cum + lefts[bid]) .= view(left, offset + chunk_size * (bid - 1) + 1:offset + chunk_size * (bid - 1) + lefts[bid])
        view(out, offset + left_sum + right_cum + 1:offset + left_sum + right_cum + rights[bid]) .= view(right, offset + chunk_size * (bid - 1) + 1:offset + chunk_size * (bid - 1) + rights[bid])
        left_cum += lefts[bid]
        right_cum += rights[bid]
    end
    return (view(out, offset + 1:offset + sum(lefts)), view(out, offset + sum(lefts)+1:offset + length(𝑖)))
end

iter = Iterators.partition(rand(5), 3)
for i in enumerate(iter)
    println(i)
end

n = Int(1e6)
nvars = 100
nbins = 64
𝑖 = collect(1:n);
𝑗 = collect(1:nvars);
X_bin = reshape(sample(UInt8.(1:nbins), n * nvars), n, nvars);
𝑖 = sample(𝑖, Int(5e5), replace=false, ordered=true);
child_bool = zeros(Bool, length(𝑖));
left = similar(𝑖)
right = similar(𝑖)
out = similar(𝑖)

offset = 0
feat = 15
cond_bin = 32
@time l, r = split_set_threads!(out, left, right, 𝑖, X_bin, feat, cond_bin, offset, 2^14);
@btime split_set_threads!($out, $left, $right, $𝑖, $X_bin, $feat, $cond_bin, $offset, 2^14);
@code_warntype split_set_1!(left, right, 𝑖, X_bin, feat, cond_bin, offset)

offset = 0
feat = 15
cond_bin = 32
lid1, rid1 = split_set_threads!(out, left, right, 𝑖, X_bin, feat, cond_bin, offset)
offset = 0
feat = 14
cond_bin = 12
lid2, rid2 = split_set_threads!(out, left, right, lid1, X_bin, feat, cond_bin, offset)
offset = + length(lid1)
feat = 14
cond_bin = 12
lid3, rid3 = split_set_threads!(out, left, right, rid1, X_bin, feat, cond_bin, offset)

lid1_ = deepcopy(lid1)


text:
using ParametricPolyhedra
using Base.Test

# equilateral triangle
a = ParametricTriangle(1.0,1.0,1.0,pi/3,pi/3,pi/3)
@test typeof(a) == ParametricTriangle{Float64}
@test !isnull(a)
@test isvalid(a,atol=eps(1.0))

# right triangle
a = ParametricTriangle(sqrt(2),1.0,1.0,pi/2,pi/4,pi/4)
@test typeof(a) == ParametricTriangle{Float64}
@test !isnull(a)
@test isvalid(a,atol=eps(1.0))

# Not valid (angles == 5pi/4)
a = ParametricTriangle(sqrt(2),1.0,1.0,pi/2,pi/4,pi/2)
@test !isvalid(a,atol=eps(1.0))

e1 = edges(pi/3,pi/3,pi/3)
@test e1 == (sin(pi/3), sin(pi/3), sin(pi/3))
e1 = edges(pi/3,pi/3,pi/3,2)
@test e1 == (2*sin(pi/3), 2*sin(pi/3), 2*sin(pi/3))
e1 = edges(3*pi/5,pi/5,pi/5)
@test e1 == (0.9510565162951536,0.5877852522924731,0.5877852522924731)

# missing value, fill tests
a = ParametricTriangle(Nullable{Float64}(),1.0,1.0,pi/2,pi/4,pi/4)
@test isnull(a)
af = fill(a)
@test isvalid(af,atol=eps(10.0))
@test isapprox(get(af.a),sqrt(2))

a = ParametricTriangle(Nullable{Float64}(),Nullable{Float64}(),Nullable{Float64}(),
                       pi/2,pi/4,pi/4)
@test isnull(a)
af = fill(a)
@test isvalid(af,atol=eps(10.0))


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/functions.jl
# Defines variable info and function types.
#############################################################################

include(joinpath(@__DIR__, "nonlinear\\auxiliary_variables.jl"))

"""
$(TYPEDEF)

An abstract super-type used for representing constraints built by EAGO's backend.
"""
abstract type AbstractEAGOConstraint end

"""

Computes the lower interval bound for `AbstractEAGOConstraint` representing an
inequality constraint.
"""
function lower_interval_bound end

"""

Computes a tuple representing the lower and upper interval bounds for a
`AbstractEAGOConstraint` representing an equality constraint.
"""
function interval_bound end

"""

Eliminate fixed variables by rearrangment or restructuring of `AbstractEAGOConstraint`.
"""
function eliminate_fixed_variables! end

###
###
### Affine function
###
###

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine inequality.
"""
mutable struct AffineFunctionIneq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end
const AFI = AffineFunctionIneq

AffineFunctionIneq() = AffineFunctionIneq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionIneq(f::SAF, s::LT)
    terms = map(x -> (x.coefficient, x.variable.value), f.terms)
    AffineFunctionIneq(terms, f.constant - s.upper, length(f.terms))
end
function AffineFunctionIneq(f::SAF, s::GT)
    terms = map(x -> (-x.coefficient, x.variable.value), f.terms)
    AffineFunctionIneq(terms, s.lower - f.constant, length(f.terms))
end
function AffineFunctionIneq(f::VI; is_max = false)
    if is_max
        return AffineFunctionIneq(Tuple{Float64,Int}[(-1.0, f.value)], 0.0, 1)
    end
    AffineFunctionIneq(Tuple{Float64,Int}[(1.0, f.value)], 0.0, 1)
end

"""
$(TYPEDEF)

Current only used for bound tightening. Stores a representation
of an affine equality.
"""
mutable struct AffineFunctionEq <: AbstractEAGOConstraint
    terms::Vector{Tuple{Float64,Int}}
    constant::Float64
    len::Int
end
const AFE = AffineFunctionEq

AffineFunctionEq() = AffineFunctionEq(Tuple{Float64,Int}[], 0.0, 0)
function AffineFunctionEq(func::SAF, set::ET)
    terms = map(x -> (x.coefficient, x.variable_index.value), func.terms)
    return AffineFunctionEq(terms, func.constant - set.value, length(func.terms))
end

###
### Parsing definitions
###

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{AffineFunctionIneq,
                                                                                AffineFunctionEq}
    deleted_count = 0
    i = 1
    while i + deleted_count <= f.len
        coeff, indx = @inbounds f.terms[i]
        variable_info = @inbounds v[indx]
        if variable_info.is_fixed
            f.constant += coeff*variable_info.lower_bound
            deleteat!(f.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count
    return nothing
end

###
###
### Quadratic function
###
###

"""
$(TYPEDEF)

Stores a general quadratic inequality constraint with a buffer.
"""
mutable struct BufferedQuadraticIneq <: AbstractEAGOConstraint
    func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
const BQI = BufferedQuadraticIneq

"""
$(TYPEDEF)

Stores a general quadratic equality constraint with a buffer.
"""
mutable struct BufferedQuadraticEq <: AbstractEAGOConstraint
    func::SQF
    minus_func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
const BQE = BufferedQuadraticEq

#=
mutable struct BufferedConvexQuadratic <: AbstractEAGOConstraint
    func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
=#

function create_buffer_dict(func::SQF)

    buffer = Dict{Int, Float64}()

    for term in func.quadratic_terms
        buffer[term.variable_1.value] = 0.0
        buffer[term.variable_2.value] = 0.0
    end

    for term in func.affine_terms
        buffer[term.variable.value] = 0.0
    end

    return buffer
end

BufferedQuadraticIneq() = BufferedQuadraticIneq(SQF(SQT[], SAT[], 0.0), Dict{Int, Float64}(), SAF(SAT[], 0.0), 0)

function BufferedQuadraticIneq(func::SQF, set::LT)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc = copy(func)
    cfunc.constant -= set.upper

    return BufferedQuadraticIneq(cfunc, buffer, saf, len)
end

function BufferedQuadraticIneq(func::SQF, set::GT)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc = MOIU.operate(-, Float64, func)
    cfunc.constant += set.lower

    return BufferedQuadraticIneq(cfunc, buffer, saf, len)
end

BufferedQuadraticEq() = BufferedQuadraticEq(SQF(SQT[], SAT[], 0.0), SQF(SQT[], SAT[], 0.0), Dict{Int, Float64}(), SAF(SAT[], 0.0), 0)

function BufferedQuadraticEq(func::SQF, set::ET)

    buffer = create_buffer_dict(func)
    saf = SAF([SAT(0.0, VI(k)) for k in keys(buffer)], 0.0)
    len = length(buffer)
    cfunc1 = copy(func)
    cfunc1.constant -= set.value
    cfunc2 = MOIU.operate(-, Float64, func)
    cfunc2.constant += set.value

    return BufferedQuadraticEq(cfunc1, cfunc2, buffer, saf, len)
end

#=
function BufferedConvexQuadratic(f::BufferedQuadraticIneq)
    BufferedConvexQuadratic(copy(f.func), copy(f.buffer), copy(f.saf), f.len)
end
=#

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{BufferedQuadraticIneq,
                                                                                BufferedQuadraticIneq}
    deleted_count = 0
    i = 1
    while i + deleted_count <= f.len
        term = f.sqf.terms[i]
        variable_info_1 = v[term.variable_1.value]
        variable_info_2 = v[term.variable_2.value]
        if variable_info_1.is_fixed && variable_info_2.is_fixed
            f.sqf.constant += coeff*variable_info_1.lower_bound*variable_info_2.lower_bound
            deleteat!(f.sqf.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count

    return nothing
end

###
###
### Cone
###
###

"""
$(TYPEDEF)

Stores a second-order cone with a buffer.
"""
mutable struct BufferedSOC <: AbstractEAGOConstraint
    variables::VECOFVAR
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end

function BufferedSOC(func::VECOFVAR, set::SOC)
    len = length(func.variables)
    buffer = Dict{Int, Float64}([(variable.value, 0.0) for variable in func.variables])
    saf = SAF(fill(SAT(0.0, VI(1)), len), 0.0)
    return BufferedSOC(copy(func), buffer, saf, len)
end

include("nonlinear/nonlinear.jl")

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module SDFResultViewer

export add_λ_units, unit_l, unit_t, unit_E, unit_B, unit_L, pₑ,
    batch_plot,
    mean_Lx_plot, Lx_section_plot,
    phase_space_summary_plot,
    average_linear_density_x, npart_plot,
    section_widget,
    compare_E_slice_with_analytic, similar_E, similar_laser

using SDFResults
using PICDataStructures
using PICDataStructures: dir_to_idx
using PICAnalysisTools
using LaserTypes
using StaticArrays
using Rotations
using RecursiveArrayTools: recursive_bottom_eltype
using Unitful
using Unitful: superscript
using PhysicalConstants.CODATA2018: c_0, ε_0, m_e, e
using Statistics, StatsBase
using ProgressLogging
# plots
import Plots
import PyPlot
# widgets
using WGLMakie
using JSServe
using JSServe: DOM, Slider
using Markdown

include("units.jl")
include("batch.jl")
include("plots/angular_momentum.jl")
include("plots/fields.jl")
include("plots/phase_space.jl")
include("plots/qed.jl")
include("plots/statistics.jl")
include("widgets.jl")
include("utils.jl")

end


text: # This file should contain all spectroscopy specifc stuff like correlators, fits, effective masses, ...


text:
"""
Spin-weighted spherical harmonics and associated functions and operators.
"""
module Sphere

#setprecision(2048) # BigFloat precision in bits

import FastGaussQuadrature as FGQ
import Jacobi: jacobi

export Y_vals, cos_vals, sin_vals, swal, swal_laplacian_matrix, swal_raising_matrix, swal_lowering_matrix, swal_filter_matrix, angular_matrix_mult!

function convertBF(v)
   return return v#parse(Float64,"$v")
end

"""
number of l angular values given ny theta values

num_l = ny - 2*(max_s+max_m), where max_s and max_m
are hard coded in Sphere.jl

num_l(
   ny::Int64
   )::Int64
"""
const max_s = 2
const max_m = 6

function num_l(ny::Int64)::Int64
   @assert ny>(2*max_s+2*max_m)+4
   return ny-2*max_s-2*max_m
end

"""
Computes Gauss-Legendre points (y) over the interval [-1,1].

Y_vals(
   ny::Int64
   )::Vector{Float64}
"""
function Y_vals(ny::Int64)::Vector{Float64}

   roots, weights= FGQ.gausslegendre(ny)

   return roots
end

"""
Computes inner product over the interval [-1,1] using
Gauss quadrature:

return: ∫ dy v1 conj(v2)

inner_product(
   v1::Vector{T},
   v2::Vector{T},
   )::T where T<:Number
"""
function inner_product(
      v1::Array{T,1},
      v2::Array{T,1}
   )::T where T<:Number
   ny = size(v1)[1]

   roots, weights= FGQ.gausslegendre(ny)

   s = 0.0
   for j=1:ny
      s += weights[j]*v1[j]*conj(v2[j])
   end

   return s
end

"""
Computes cos(y) at Gauss-Legendre points over interval [-1,1].

cos_vals(
   ny::Int64
   )::Vector{Float64}
"""
function cos_vals(ny::Int64)::Vector{Float64}

   roots, weights= FGQ.gausslegendre(ny)

   return [-pt for pt in roots]
end

"""
Computes sin(y) at Gauss-Legendre points over interval [-1,1].

sin_vals(
   ny::Int64
   )::Vector{Float64}
"""
function sin_vals(ny::Int64)::Vector{Float64}

   roots, weights= FGQ.gausslegendre(ny)

   return [convert(Float64,
                   sqrt(1.0-convertBF(pt))*sqrt(1.0+convertBF(pt))
                  )
           for pt in roots
          ]
end

"""
Computes the spin-weighted associated Legendre function Y^s_{lm}(y).

swal(
   spin::Int64,
   m_ang::Int64,
   l_ang::Int64,
   y::Float64
   )::Float64
"""
function swal(
      spin::Int64,
      m_ang::Int64,
      l_ang::Int64,
      y::Float64
   )::Float64
   @assert l_ang>=abs(m_ang)

   al = abs(m_ang-spin)
   be = abs(m_ang+spin)
   @assert((al+be)%2==0)
   n = l_ang - (al+be)/2

   if n<0
      return convert(Float64,0)
   end

   norm = sqrt(
      (2*n+al+be+1)*(convertBF(2)^(-al-be-1.0))
   *  factorial(convertBF(n+al+be))/factorial(convertBF(n+al))
   *  factorial(convertBF(n      ))/factorial(convertBF(n+be))
   )
   norm *= convertBF(-1)^(max(m_ang,-spin))

   return convert(
      Float64,
      norm
      *(convertBF(1-y)^(al/2.))
      *(convertBF(1+y)^(be/2.))
      *jacobi(convertBF(y),convertBF(n),convertBF(al),convertBF(be))
   )
end

"""
Computes matrix swal^s_{lm}(y) at Gauss-Legendre points,
over a grid of l values (and fixed m).

swal_vals(
   ny::Int64,
   spin::Int64,
   m_ang::Int64,
   )::Matrix{Float64}
"""
function swal_vals(
      ny::Int64,
      spin::Int64,
      m_ang::Int64,
   )::Matrix{Float64}

   roots, weights= FGQ.gausslegendre(ny)

   nl = num_l(ny)

   vals = zeros(Float64,ny,nl)

   lmin = max(abs(spin),abs(m_ang))

   for k=1:nl
      l_ang = k-1+lmin
      for j=1:length(roots)
         vals[j,k] = swal(spin,m_ang,l_ang,roots[j])
      end
   end
   return vals
end

"""
Computes matrix to compute spin-weighted spherical harmonic laplacian.

swal_laplacian_matrix(
   ny::Int64,
   spin::Int64,
   m_ang::Int64
   )::Matrix{Float64}
"""
function swal_laplacian_matrix(
      ny::Int64,
      spin::Int64,
      m_ang::Int64
   )::Matrix{Float64}

   yv, wv = FGQ.gausslegendre(ny)

   nl   = num_l(ny)
   lmin = max(abs(spin),abs(m_ang))

   lap = zeros(Float64,ny,ny)

   for j=1:ny
      for i=1:ny
         for k=1:nl
            l = k-1+lmin
            lap[j,i] -= (l-spin)*(l+spin+1.0)*(swal(spin,m_ang,l,yv[i])*
                                               swal(spin,m_ang,l,yv[j])
                                              )
         end
         lap[j,i] *= wv[j]
      end
   end

   return lap
end

"""
Computes matrix to raise spin-weighted spherical harmonics.

swal_raising_matrix(
   ny::Int64,
   spin::Int64,
   m_ang::Int64
   )::Matrix{Float64}
"""
function swal_raising_matrix(
      ny::Int64,
      spin::Int64,
      m_ang::Int64
   )::Matrix{Float64}

   yv, wv = FGQ.gausslegendre(ny)

   nl   = num_l(ny)
   lmin = max(abs(spin), abs(spin+1), abs(m_ang))

   raise = zeros(Float64,ny,ny)

   for (i,yi) in enumerate(yv)
      for (j,yj) in enumerate(yv)

         for k=1:nl
            l = k-1+lmin
            raise[j,i] += (
                  sqrt((l-spin)*(l+spin+1.0))
                  *
                  wv[j]*swal(spin,m_ang,l,yj)
                  *
                  swal(spin+1,m_ang,l,yi)
                 )
         end
      end
   end

   return raise
end

"""
Computes matrix to lower spin-weighted spherical harmonics.

swal_lowering_matrix(
   ny::Int64,
   spin::Int64,
   m_ang::Int64
   )::Matrix{Float64}
"""
function swal_lowering_matrix(
      ny::Int64,
      spin::Int64,
      m_ang::Int64
   )::Matrix{Float64}

   yv, wv = FGQ.gausslegendre(ny)

   nl   = num_l(ny)
   lmin = max(abs(spin), abs(spin-1), abs(m_ang))

   lower = zeros(Float64,ny,ny)

   for (i,yi) in enumerate(yv)
      for (j,yj) in enumerate(yv)

         for k=1:nl
            l = k-1+lmin
            lower[j,i] += (
               -
               sqrt((l+spin)*(l-spin+1.0))
               *
               wv[j]*swal(spin,m_ang,l,yj)
               *
               swal(spin-1,m_ang,l,yi)
              )
         end
      end
   end

   return lower
end

"""
Computes matrix to compute low pass filter.

swal_filter_matrix(
   ny::Int64,
   spin::Int64,
   m_ang::Int64
   )::Matrix{Float64}

Multiply the matrix on the left: v[i]*M[i,j] -> v[j]
"""
function swal_filter_matrix(
      ny::Int64,
      spin::Int64,
      m_ang::Int64
   )::Matrix{Float64}

   yv, wv = FGQ.gausslegendre(ny)

   nl = num_l(ny)

   filter = zeros(Float64,ny,ny)

   lmin = max(abs(spin),abs(m_ang))

   for j=1:ny
      for i=1:ny
         for l=lmin:(nl-1+lmin)
            filter[j,i] += exp(-30.0*(l/(nl-1.0))^10)*(
                           swal(spin,m_ang,l,yv[i])*
                           swal(spin,m_ang,l,yv[j])
                        )
         end
         filter[j,i] *= wv[j]
      end
   end

   return filter
end

"""
Matrix multiplication in the angular direction

angular_matrix_mult!(
   f_m,
   f,
   mat
   )
"""
function angular_matrix_mult!(
      f_m,
      f,
      mat
   )
   nx, ny = size(f)
   for j=1:ny
      for i=1:nx
         f_m[i,j] = 0.0
         for k=1:ny
            f_m[i,j] += f[i,k]*mat[k,j]
         end
      end
   end
   return nothing
end

end


text:
module FeaturedFuncs

using ..PerFlagFuncs
using ..UtilityFuncs

export Flag
export FlagSet

export get_flat
export get_level
export get_level_flat
export remove_flat

"""
```
get_flat(s::String,
                  flags_start::Vector{String},
                  flags_stop::Vector{String},
                  inclusive::Bool = true)
```
Gets the substring based on the start and stop flag
vectors, and a Bool which determines whether the flags
themselves should be returned or not.

This function will grab the inner-most string, assuming
that you do not have multiple start flags before reaching
a corresponding stop flag.

```
julia> s = "Some text... {GRAB THIS}, some more text {GRAB THIS TOO}..."
"Some text... {GRAB THIS}, some more text {GRAB THIS TOO}..."

julia> L = BetweenFlags.get_flat(s, ["{"], ["}"])
2-element Array{String,1}:
 "{GRAB THIS}"
 "{GRAB THIS TOO}"
```
"""
function get_flat(s::String,
                  flags_start::Vector{String},
                  flags_stop::Vector{String},
                  inclusive::Bool = true)
  L = Vector{String}()
  L_start = [m for flag_start in flags_start for m in find_next_iter(s, flag_start)]
  L_stop  = [m for flag_stop  in flags_stop  for m in find_next_iter(s, flag_stop )]
  sort!(L_start)
  sort!(L_stop)
  if L_start==[] || L_stop==[]
    return Vector{String}([""])
  end
  (L_start, L_stop) = get_alternating_consecutive_vector(L_start, L_stop)
  for (i_start, i_stop) in zip(L_start, L_stop)
    b, m, a = substring_decomp_by_index(s, i_start, i_stop, flags_start, flags_stop, inclusive)
    push!(L, m)
  end
  return L
end

"""
```
get_level_flat(s::String,
               flags_start::Vector{String},
               flags_stop::Vector{String},
               inclusive::Bool = true)
```
get_level_flat gets the substring based on the `flags_start`
and `flags_stop` vectors, and a Bool which determines
whether the flags themselves should be returned or not.

This function will grab the outer-most string by ignoring
stop flags when multiple start flags occur before stop flags.

```
julia> using BetweenFlags

julia> s = "Some text... {GRAB {THIS}}, some more text {GRAB THIS TOO}..."
"Some text... {GRAB {THIS}}, some more text {GRAB THIS TOO}..."

julia> L = BetweenFlags.get_level_flat(s, ["{"], ["}"])
2-element Array{String,1}:
 "{GRAB {THIS}}"
 "{GRAB THIS TOO}"
```
"""
function get_level_flat(s::String,
                        flags_start::Vector{String},
                        flags_stop::Vector{String},
                        inclusive::Bool = true)
  L = Vector{String}()
  L_start = [m for flag_start in flags_start for m in find_next_iter(s, flag_start)]
  L_stop  = [m for flag_stop  in flags_stop  for m in find_next_iter(s, flag_stop )]
  L_start = unique(L_start)
  L_stop = unique(L_stop)
  if L_start==[] || L_stop==[]
    return Vector{String}([""])
  end
  sort!(L_start)
  sort!(L_stop)
  level = compute_level_total(s, flags_start, flags_stop)
  (L_start, L_stop) = get_alternating_consecutive_vector(L_start, L_stop, level)
  for (i_start, i_stop) in zip(L_start, L_stop)
    b, m, a = substring_decomp_by_index(s, i_start, i_stop, flags_start, flags_stop, inclusive)
    push!(L, m)
  end
  return L
end

"""
```
get_level(s::String,
          outer_flags::FlagSet,
          inner_flags::Vector{FlagSet},
          inclusive::Bool = true)
```
This is the featured function of BetweenFlags.

Gets the substring based on the outer and inner flag
sets, and a Bool which determines whether the flags
themselves should be returned or not.

To see an example of this function in action, go to
BetweenFlags/test/runtests.jl.
"""
function get_level(s::String,
                   outer_flags::FlagSet,
                   inner_flags::Vector{FlagSet},
                   inclusive::Bool = true)
  L = Vector{String}()
  flag_set_all = vcat([outer_flags], inner_flags)
  outer_flags_start = outer_flags.start.trigger
  outer_flags_stop = outer_flags.stop.trigger
  inner_flags_start = [y for x in inner_flags for y in x.start.trigger]
  inner_flags_stop  = [y for x in inner_flags for y in x.stop.trigger]

  flags_start = vcat(outer_flags_start, inner_flags_start)
  flags_stop  = vcat(outer_flags_stop , inner_flags_stop)
  L_start = [m for flag_start in flags_start for m in find_next_iter(s, flag_start)]
  L_stop  = [m for flag_stop  in flags_stop  for m in find_next_iter(s, flag_stop )]
  L_start = unique(L_start)
  L_stop = unique(L_stop)
  if L_start==[] || L_stop==[]
    return Vector{String}([""])
  end
  sort!(L_start)
  sort!(L_stop)
  level_total = compute_level_total(s, flags_start, flags_stop)
  level_per_flags = compute_level_per_flag(s, level_total, flag_set_all)
  level_outer = level_per_flags[outer_flags.ID]
  (L_start, L_stop) = get_alternating_consecutive_vector(L_start, L_stop, level_total, level_outer, s)
  for (i_start, i_stop) in zip(L_start, L_stop)
    b, m, a = substring_decomp_by_index(s, i_start, i_stop, flags_start, flags_stop, inclusive)
    push!(L, m)
  end
  return L
end

"""
```
remove_flat(s::String,
            flags_start::Vector{String},
            flags_stop::Vector{String},
            inclusive::Bool = true)
```
Removes text between flags.

The remove function is fundamentally different from get_level because
the string, `s`, in `get_level` does not change, whereas it does in
`remove_flat`. Therefore, the indexes found must, either be translated
by the number of removed characters in the correct location, or the
entire function must be called recursively. Alternatively, the
strings/flags can be removed in reverse order, preserving the output
string, which is what is done here. In addition, the return type of
`remove_flat` is a String.
"""
function remove_flat(s::String,
                     flags_start::Vector{String},
                     flags_stop::Vector{String},
                     inclusive::Bool = true)::String
  if !get_remaining_flags(s, flags_start, flags_stop)
    return s
  end
  same_flags = all([x==y for x in flags_start for y in flags_stop])
  L_start = [m for flag_start in flags_start for m in find_next_iter(s, flag_start)]
  L_stop  = [m for flag_stop  in flags_stop  for m in find_next_iter(s, flag_stop )]
  sort!(L_start)
  sort!(L_stop)
  (L_start, L_stop) = get_alternating_consecutive_vector(L_start, L_stop)
  if same_flags
    L_start = L_start[1:2:end]
    L_stop = L_stop[1:2:end]
  end
  s_new = s
  for (i_start, i_stop) in zip(reverse(L_start), reverse(L_stop))
    b, m, a = substring_decomp_by_index(s_new, i_start, i_stop, flags_start, flags_stop, inclusive)
    s_new = string(b, a)
  end
  return s_new
end

end


text:
#cd("/home/jchang/image/result/")
#
#print("Date: ")
#date=chomp(readline())
##date="20190315"
#@show date
#
#files=readdir()
#filelist=filter(x->occursin("_1--b.tif",x),files)
#filelist=filter(x->occursin(@sprintf("%s",date),x),filelist)
#@show filelist
#
#print("Filename: ")
#file=chomp(readline())
#@show file
#filelist=filter(x->occursin("--b.tif",x),files)
#filelist=filter(x->occursin(@sprintf("%s",file),x),filelist)
#@show filelist
#
#
#print("Frame start: ")
#z_start=chomp(readline())
#z_start=parse(Int,z_start)
#
#print("Frame end: ")
#z_end=chomp(readline())
#z_end=parse(Int,z_end)
#
#print("Is bidirectional(Y:1,N:0): ")
#bi=chomp(readline())
#bi=parse(Int,bi)
#
#filename=filelist[1][1:end-4]
#
#img = load(@sprintf("%s.tif", filename))
#img = convert(Array{N0f16}, img)
#
#(x,y,z) = size(img)
#
#
#@showprogress "Building time series file..." for s in z_start:z_end
#
#	if bi == 1
#		img_re = zeros(Int(floor(x/2)),y,size(filelist,1)*2)
#
#		for k in 1:size(filelist, 1)
#			filename=filelist[k][1:end-4]
#
#			img = load(@sprintf("%s.tif", filename))
#			img = convert(Array{N0f16}, img)
#
#			img_re[:,:,k*2-1] = img[1:Int(floor(x/2)),:,s]
#			img_re[end:-1:1,:,k*2] = img[end-Int(floor(x/2)):end-1,:,s]
#		end
#
#	else
#		img_re = zeros(x,y,size(filelist,1))
#
#		for k in 1:size(filelist, 1)
#			filename=filelist[k][1:end-4]
#
#			img = load(@sprintf("%s.tif", filename))
#			img = convert(Array{N0f16}, img)
#
#			img_re[:,:,k] = img[:,:,s]
#		end
#	end
#
#	img_re = Gray.(convert.(Normed{UInt16,16},img_re))
#	FileIO.save(File{format"TIFF"}(@sprintf("%s-%s-t.tif", file,s)), img_re)
#end
#
#cd("/home/jchang/image/")


text:
    # Functions and methods for cascade computation

    """
    `Cascade.computeSteps(scheme::Cascade.PhotonIonizationScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested transition amplitudes as well as PhotoIonization.Line's, etc. for all pre-specified
            (photo-) ionizing steps of the cascade. When compared with standard computations of these atomic processes, however,
            the amount of output is largely reduced and often just printed into the summary file. A set of
            data::Cascade.PhotoIonData  is returned.
    """
    function computeSteps(scheme::Cascade.PhotonIonizationScheme, comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        data = Cascade.PhotoIonData[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")

        # Perform computations in turn for all photon energies
        for  photonEnergy in  scheme.photonEnergies
            nt = 0;   st = 0;   linesP = PhotoIonization.Line[]
            settings = PhotoIonization.Settings(PhotoIonization.Settings(), photonEnergies=[Defaults.convertUnits("energy: from atomic", photonEnergy)])
            println("\n\n* Perform ionizing computations for the photon energy $photonEnergy  [a.u.]")
            if  printSummary   println(iostream, "\n\n* Perform ionizing computations for the photon energy $photonEnergy  [a.u.]")   end
            #
            for  step  in  stepList
                st = st + 1
                nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
                println("\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc ionizing lines (without selection rules): ")
                if  printSummary   println(iostream, "\n* $st) Perform $(string(step.process)) amplitude computations for " *
                                                     "up to $nc ionizing lines (without selection rules): ")   end

                if      step.process == Basics.Photo()
                    newLines = PhotoIonization.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                                   settings, output=true, printout=true)
                    append!(linesP, newLines);    nt = length(linesP)
                else   error("Unsupported atomic process for cascade computations.")
                end
                println("     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                        "to a total of $nt $(string(step.process)) ionizing lines." )
                if  printSummary   println(iostream, "\n*    Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, " *
                                                     "giving now rise to a total of $nt $(string(step.process)) ionizing lines." )   end
            end
            #
            push!(data, Cascade.PhotoIonData(photonEnergy, linesP))
        end

        return( data )
    end

    """
    `Cascade.determineSteps(scheme::Cascade.PhotonIonizationScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, ionizedList::Array{Cascade.Block,1})`
        ... determines all steps::Cascade.Step's that need to be computed for this ionizing cascade. It cycles through all processes
            of the given ionizing scheme and selects all pairs of initial and ionized blocks due to the selected photon energies and
            cascade approaches.             It also checks that the (averaged) energies of the configuration support a photoionization
            energetically.             A stepList::Array{Cascade.Step,1} is returned
    """
    function determineSteps(scheme::Cascade.PhotonIonizationScheme, comp::Cascade.Computation,
                            initialList::Array{Cascade.Block,1}, ionizedList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            for  initialBlock in initialList
                for  ionizedBlock in ionizedList
                    for  process  in  comp.scheme.processes
                        if      process == Basics.Photo()
                            if  initialBlock.NoElectrons == ionizedBlock.NoElectrons + 1
                                settings = PhotoIonization.Settings()   ##x settings, photonEnergies=[1.0], printBefore=false)
                                push!( stepList, Cascade.Step(process, settings, initialBlock.confs, ionizedBlock.confs,
                                                                                initialBlock.multiplet, ionizedBlock.multiplet) )
                            end
                        else    error("stop a")
                        end
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.generateConfigurationsForPhotoionization(multiplets::Array{Multiplet,1},  maxPhotoElectrons::Int64,
                                                      maxPhotonEnergy::Float64, nm::Nuclear.Model)`
        ... generates all possible (photoionized) configurations with upto maxPhotoElectrons less electrons and whose averaged energy
            is NOT higher than maxPhotoEnergy above of the given configuration. In the present version, no shake contributions are taken
            into account. A Tuple(initialConfList::Array{Configuration,1}, confList::Array{Configuration,1}) is returned.
    """
    function generateConfigurationsForPhotoionization(multiplets::Array{Multiplet,1}, maxPhotoElectrons::Int64, maxPhotonEnergy::Float64,
                                                      nm::Nuclear.Model)
        # Determine all (different) configurations from multiplets
        newconfList = Configuration[]
        for mp  in  multiplets
            confList = Basics.extractNonrelativisticConfigurations(mp.levels[1].basis)
            for  conf in confList   if  conf in newconfList   nothing   else   push!(newconfList, conf)      end      end
        end
        en     = Float64[];   for conf in confList    push!(en, -Semiempirical.estimate("binding energy", round(Int64, nm.Z), conf))   end
        maxen  = maximum(en)
        shList = Basics.generate("shells: ordered list for NR configurations", newconfList)

        # Collect all configuration with up to maxPhotoElectrons less electrons
        initialConfList = deepcopy(newconfList);
        allnewconfList  = Configuration[];    confList = deepcopy(newconfList);    newconfList = Configuration[]
        for ne = 1:maxPhotoElectrons
            for  sh in shList
                for conf in confList
                    shells = deepcopy(conf.shells);
                    if shells[sh] > 0   shells[sh] = shells[sh] - 1;     push!(newconfList, Configuration(shells, conf.NoElectrons - 1))  end
                end
            end
            append!(allnewconfList, newconfList)
            confList = deepcopy(newconfList);    newconfList = Configuration[]
        end

        # Exclude double configurations as well as those with too high average energy
        ##x confList = Basics.excludeDoubles(allnewconfList)
        confList = unique(allnewconfList)
        newconfList = Configuration[]
        for  conf in confList
            enconf = -Semiempirical.estimate("binding energy", round(Int64, nm.Z), conf)
            if enconf - maxen < maxPhotonEnergy   push!(newconfList, conf)    end
        end

        return( (initialConfList, newconfList)  )
    end

    """
    `Cascade.perform(scheme::PhotonIonizationScheme, comp::Cascade.Computation)`
        ... to set-up and perform a cascade computation that starts from a given set of initial configurations xor initial multiplets
            and proceeds via photoionizing processes to configurations with a given No of photoelectrons. The results of all individual
            photoionization steps are comprised into (output) data::PhotoIonData, although all data are only printed and nothing
            is returned.

    `Cascade.perform(scheme::PhotonIonizationScheme, comp::Cascade.Computation; output=true, outputToFile::Bool=true)`
        ... to perform the same but to return the complete output in a dictionary that is written to disk and can be used in subsequent
            cascade simulation. The particular output also depends on the specifications of the cascade.
    """
    function perform(scheme::PhotonIonizationScheme, comp::Cascade.Computation; output::Bool=false, outputToFile::Bool=true)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplets if initial configurations are given
        if  comp.initialConfigs != Configuration[]
            basis      = Basics.performSCF(comp.initialConfigs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplet  = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
            multiplets = [Multiplet("initial states", multiplet.levels)]
        else
            multiplets = comp.initialMultiplets
        end
        # Print out initial configurations and levels
        Cascade.displayLevels(stdout, multiplets, sa="initial ")
        if  printSummary   Cascade.displayLevels(iostream, multiplets, sa="initial ")                            end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        maxen = maximum(comp.scheme.photonEnergies) + 20
        wa  = Cascade.generateConfigurationsForPhotoionization(multiplets, comp.scheme.maxPhotoElectrons, maxen, comp.nuclearModel)
        wb1 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[1], sa="(initial part of the) ionizing ")
        wb2 = Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa[2], sa="(generated part of the) ionizing ")
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        wc1 = Cascade.generateBlocks(comp, wb1, basis.orbitals, sa="for the ionizing cascade ")
        wc2 = Cascade.generateBlocks(comp, wb2, basis.orbitals, sa="for the ionizing cascade ", printout=false)
        Cascade.displayBlocks(stdout, wc1, sa="for the (initial part of the) ionizing cascade ");
        Cascade.displayBlocks(stdout, wc2, sa="for the (generated part of the) ionizing cascade ")
        if  printSummary   Cascade.displayBlocks(iostream, wc1, sa="for the (initial part of the) ionizing cascade ")
                           Cascade.displayBlocks(iostream, wc2, sa="for the (generated part of the) ionizing cascade ")    end
        # Determine, modify and compute the transition data for all steps, ie. the PhotoIonization.Line's, etc.
        gMultiplets = Multiplet[];     for block in wc2  push!(gMultiplets, block.multiplet)    end
        we = Cascade.determineSteps(scheme, comp, wc1, wc2)
        Cascade.displaySteps(stdout, we, sa="ionizing ")
        if  printSummary   Cascade.displaySteps(iostream, we, sa="ionizing ")    end
        wf   = Cascade.modifySteps(we)
        data = Cascade.computeSteps(scheme, comp, wf)
        if output
            results = Base.merge( results, Dict("name"                       => comp.name) )
            results = Base.merge( results, Dict("cascade scheme"             => comp.scheme) )
            results = Base.merge( results, Dict("initial multiplets:"        => multiplets) )
            results = Base.merge( results, Dict("generated multiplets:"      => gMultiplets) )
            results = Base.merge( results, Dict("photo-ionizing line data:"  => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "zzz-cascade-ionizing-computations-" * string(Dates.now())[1:13] * ".jld"
            println("\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                    "\n   results = JLD.load(''$filename'')    ... to load the results back from file.")
            if  printSummary   println(iostream, "\n* Write all results to disk; use:\n   JLD.save(''$filename'', results) \n   using JLD " *
                                                 "\n   results = JLD.load(''$filename'')    ... to load the results back from file." )      end
            JLD2.@save filename results
        end
        ## return( results )
        return( results )
    end


text:
#folder = "CUTEst_infeasible"
folder = "CUTEst"
#folder = "CUTEst_low_tol"
#folder = "netlib_infeasible"

if folder == "CUTEst"
  mode = :Optimal
  data = :CUTEst
elseif folder == "CUTEst_low_tol"
  mode = :Optimal
  data = :CUTEst_low_tol
elseif folder == "CUTEst_infeasible"
  mode = :primal_infeasible
  data = :CUTEst
elseif folder == "netlib_infeasible"
  mode = :primal_infeasible
  data = :netlib
end

line_styles = Dict("IPOPT" => "-", "one-phase" => "--")
line_colors = Dict("IPOPT" => "black", "one-phase" => "black")

function get_CUTEst_results()
  results = Dict{String, Dict{String,problem_summary2}}()

  #GLOBAL mode;
  if mode == :Optimal
    if data == :CUTEst
    #results["IPOPT"] = load("../results/one_phase/sept_3_corrections/summary.jld", "summary")
      #results["one-phase"] = cps(load("../results/one_phase/Dec24/summary.jld", "summary"))
      results["one-phase"] = cps(load("../results/one_phase/Jan4/CUTEst/summary.jld", "summary"))
      #results["one-phase"] = cps(load("../results/one_phase/test_knitro_init7_Jan1/summary.jld", "summary"))
      #results["IPOPT"] = cps(load("../results/one_phase/test_knitro_init8_Jan2/summary.jld", "summary"))
      #results["one-phase"] = cps(load("../results/one_phase/test_meh_init_Dec30/summary.jld", "summary"))

      #results["IPOPT"] = convert_JuMP(cps(load("../results/ipopt/plain/summary.jld", "summary")))
    elseif data == :CUTEst_low_tol
      results["one-phase"] = cps(load("../results/one_phase/Jan4/LowTol/summary.jld", "summary"))
      results["IPOPT"] = convert_JuMP(cps(load("../results/ipopt/very_low_tol/summary.jld", "summary")))
    end
  elseif mode == :primal_infeasible
    if data == :CUTEst
      #results["one-phase"] = cps(load("../results/one_phase/infeas-Dec20/summary.jld", "summary"))
      results["one-phase"] = cps(load("../results/one_phase/Jan4/infeas/summary.jld", "summary"))
      results["IPOPT"] = convert_JuMP(cps(load("../results/ipopt/infeas-Dec15/summary.jld", "summary")))
    elseif data == :netlib
      results["one-phase"] = cps(load("../results/one_phase/Jan4/infeas-netlib/summary.jld", "summary"))
      results["IPOPT"] = convert_JuMP(cps(load("../results/ipopt/Jan4/infeas-netlib/summary.jld", "summary")))
    end
  else
    error("")
  end

  #=min_num = 0
  max_num = 10000
  lrg_problems = lrg_problem_func(min_num, max_num)
  problem_list = CUTEst.select(custom_filter=lrg_problems);
  problem_list = convert(Array{String,1},problem_list);

  results = restrict_problems(results, problem_list)=#

  overlapping_results = overlap(results)

  return overlapping_results
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
using MacroTools: @forward

#####
# Utils
#####

using IntervalSets

Random.rand(s::Union{Interval,Array{<:Interval}}) = rand(Random.GLOBAL_RNG, s)

function Random.rand(rng::AbstractRNG, s::Interval)
    rand(rng) * (s.right - s.left) + s.left
end

#####
# WorldSpace
#####

#####
# ZeroTo
#####

export ZeroTo

"""
Similar to `Base.OneTo`. Useful when wrapping third-party environments.
"""
struct ZeroTo{T<:Integer} <: AbstractUnitRange{T}
    stop::T
    ZeroTo{T}(n) where {T<:Integer} = new(max(zero(T) - one(T), n))
end

ZeroTo(n::T) where {T<:Integer} = ZeroTo{T}(n)

Base.show(io::IO, r::ZeroTo) = print(io, "ZeroTo(", r.stop, ")")
Base.length(r::ZeroTo{T}) where {T} = T(r.stop + one(r.stop))
Base.first(r::ZeroTo{T}) where {T} = zero(r.stop)

function getindex(v::ZeroTo{T}, i::Integer) where {T}
    Base.@_inline_meta
    @boundscheck ((i >= 0) & (i <= v.stop)) || throw_boundserror(v, i)
    convert(T, i)
end

#####
# Space
#####

#####
# Generate README
#####

gen_traits_table(envs) = gen_traits_table(stdout, envs)

function gen_traits_table(io, envs)
    trait_dict = Dict()
    for f in RLBase.env_traits()
        for env in envs
            if !haskey(trait_dict, f)
                trait_dict[f] = Set()
            end
            t = f(env)
            if f == StateStyle
                if t isa Tuple
                    for x in t
                        push!(trait_dict[f], nameof(typeof(x)))
                    end
                else
                    push!(trait_dict[f], nameof(typeof(t)))
                end
            else
                push!(trait_dict[f], nameof(typeof(t)))
            end
        end
    end

    println(io, "
")

    print(io, "Traits	")
    for i in 1:length(envs)
        print(io, " $(i) 	")
    end

    for k in sort(collect(keys(trait_dict)), by = nameof)
        vs = trait_dict[k]
        print(io, "
  $(nameof(k)) 	")
        for (i, v) in enumerate(vs)
            if i != 1
                print(io, "
 ")
            end
            print(io, " $(v) 	")
            for env in envs
                if k == StateStyle && k(env) isa Tuple
                    ss = k(env)
                    if v in map(x -> nameof(typeof(x)), ss)
                        print(io, " ✔ 	")
                    else
                        print(io, " 	 ")
                    end
                else
                    if nameof(typeof(k(env))) == v
                        print(io, " ✔ 	")
                    else
                        print(io, " 	 ")
                    end
                end
            end
            println(io, "
")
        end
    end

    println(io, "
")

    print(io, "
")
    for env in envs
        println(
            io,
            "
  $(nameof(env))
",
        )
    end
    print(io, "
")
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.7.2
#     language: julia
#     name: julia-1.7
# ---

using ArgoData, DataFrames, CSV, NCDatasets, Plots

# ## Input File
#
# 1. Choose one Argo float (e.g., `wmo=6900900`)
# 1. Download file to temporary folder (if needed).
# 1. Set up file names and parameters (`meta`)
# 1. Open file for lazy access (`data`)

# +
wmo=6900900

url0="https://data-argo.ifremer.fr/dac/coriolis/"
input_url=url0*"/$(wmo)/$(wmo)_prof.nc"
input_file=joinpath(tempdir(),"$(wmo)_prof.nc")
output_file=joinpath(tempdir(),"$(wmo)_MITprof.nc")

!isfile(input_file) ? fil=Downloads.download(input_url,input_file) : nothing

meta=ArgoTools.meta(input_file,output_file)

fil=joinpath(tempdir(),"ar_greylist.txt")
isfile(fil) ? greylist=GDAC.grey_list(fil) : greylist=""

meta["greylist"]=greylist

data=Dataset(input_file)
# -

# ## Read Profile
#
# Let's read the first T/S profile from the input file (`prof`) and set up the data structure (`prof_std`) that will be used when interpolating to standard levels.

m=1
prof=ArgoTools.GetOneProfile(data,m);

nz=length(meta["z_std"])
prof_std=ArgoData.ProfileStandard(nz);

# +
#for verification, record intermediate step:
prof_step0=(T=prof.T, S=prof.S, P=prof.pressure)

#display date
data["JULD"][m]
# -

# ## Process Profile
#
# 1. Conversions, as specified in `meta`.
# 1. Interpolation to standard depth levels.

# +
#various conversions (as specified in `meta`):
ArgoTools.prof_convert!(prof,meta)

#for verification, record intermediate step:
prof_step1=(T=prof.T, S=prof.S, D=prof.depth);
# -

#interpolate to standard depth levels
ArgoTools.prof_interp!(prof,prof_std,meta)

#first series of tests
ArgoTools.prof_test_set1!(prof,prof_std,meta)

# ## Verification / CI

#verification / CI
k=findall((!ismissing).(prof_step0.P))[50]
println(prof_step0.P[k]==937.0)

prof_step0.T[k]

#verification / CI
println(prof_step1.T[k]==2.8032383695780343)

# ## Verification / Visual

scatter(prof_step0.S,prof_step0.T,title="Temperature-Salinity relationship",
    label="original levels",xlabel="psu",ylabel="degree C")
scatter!(prof_std.S,prof_std.T,color=:red,
    label="standard levels",legend = :topright)

# + {"cell_style": "center"}
scatter(prof_step1.T,-prof_step1.D,title="Potential Temperature",
    label="original levels",xlabel="degree C",ylabel="m")
scatter!(prof_std.T,-meta["z_std"],color=:red,ylims=(-2000.0,0.0),
    label="standard levels",legend = :bottomright)

# + {"cell_style": "center"}
scatter(prof_step1.S,-prof_step1.D,title="Salinity",
    label="original levels",xlabel="psu",ylabel="m")
scatter!(prof_std.S,-meta["z_std"],color=:red,ylims=(-2000.0,0.0),
    label="standard levels",legend = :topright)
# -


text:
#=
Indexer struct for array type entities

An indexer has a symbol, range and current value
The range can be numbers or other indexers.
It can be a min,max pair, or an array of values.
=#

mutable struct Indexer
    symbol::Symbol   # symbol used in expressions
    range::Array     # Range can be a min/max pair, or an array of values.
    value::Int       # Current value (may not be needed?)
end


text:
function paraview_collection(filename::AbstractString)
    # Initialise Paraview collection file (extension .pvd).
    # filename: filename with or without the extension (.pvd).
    xpvd = XMLDocument()
    xroot = create_root(xpvd, "VTKFile")
    set_attribute(xroot, "type", "Collection")
    set_attribute(xroot, "version", "1.0")
    if IS_LITTLE_ENDIAN
        set_attribute(xroot, "byte_order", "LittleEndian")
    else
        set_attribute(xroot, "byte_order", "BigEndian")
    end
    set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    new_child(xroot, "Collection")
    return CollectionFile(xpvd, add_extension(filename, ".pvd"))
end

function paraview_collection_load(filename::AbstractString)
    pvd_filename = add_extension(filename, ".pvd")
    xpvd = parse_file(pvd_filename)
    pvd = paraview_collection(filename)
    xroot = root(pvd.xdoc)
    xMBDS = find_element(xroot, "Collection")
    # Iterate the child elements and only add
    # attributes considered by pvd.
    # This also preserves the formatting of the resulting
    # file as the empty textnodes created during loading
    # are removed.
    for c in child_elements(find_element(root(xpvd), "Collection"))
        xDataSet = new_child(xMBDS, "DataSet")
        set_attribute(xDataSet, "timestep", attribute(c, "timestep"))
        set_attribute(xDataSet, "part", attribute(c, "part"))
        set_attribute(xDataSet, "file", attribute(c, "file"))
    end
    return pvd
end

function collection_add_timestep(pvd::CollectionFile, datfile::VTKFile,
                                 t::Real)
    xroot = root(pvd.xdoc)
    xMBDS = find_element(xroot, "Collection")
    xDataSet = new_child(xMBDS, "DataSet")
    fname = relpath(abspath(datfile.path), first(splitdir(abspath(pvd.path))))
    set_attribute(xDataSet, "timestep", string(t))
    set_attribute(xDataSet, "part", "0")
    set_attribute(xDataSet, "file", fname)
    append!(pvd.timeSteps, vtk_save(datfile))
    return
end

function vtk_save(pvd::CollectionFile)
    # Saves paraview collection file (.pvd).
    outfiles = [pvd.path; pvd.timeSteps]::Vector{String}
    if isopen(pvd)
        save_file(pvd.xdoc, pvd.path)
        close(pvd)
    end
    return outfiles
end


text:
# test for primitives.jl
using Test
#include("primitives.jl")

@testset "vars and const" begin
 @test isvar(:x, [:y,:x]) == true
 @test isvar(:y, [:y,:x]) == true
 @test isvar(:z, [:y,:x]) == false
 @test isvar(12, [:y,:x]) == false
 @test isvar(23.1, [:y,:x]) == false
 @test isvar(0, [:y,:x]) == false
 @test isvar([:a :b], [:y,:x]) == false
 @test isvar([], [:y,:x]) == false

 @test isevar(:Abc, [:Abc, :x]) == true
 @test isevar(:Abc, [:abc, :x]) == false
 @test isevar(:dbc, [:abc, :x]) == false

 @test isconst(:x, [:y,:x]) == false
 @test isconst(:a, [:y,:x]) == true
 @test isconst(12, [:y,:x]) == true
 @test isconst(23.1, [:y,:x]) == true
 @test isconst(0, [:y,:x]) == true
 @test isvar([:a :b], [:y,:x]) == false
 @test isvar([], [:y,:x]) == false

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
ix = 1.0

z = "My String"

# http://docs.julialang.org/en/release-0.2/manual/integers-and-floating-point-numbers/
for T = {Int8,Int16,Int32,Int64,Int128,Uint8,Uint16,Uint32,Uint64,Uint128}
  println("$(lpad(T,6)): [$(typemin(T)),$(typemax(T))]")
end

println(bits(10))
println(bits(10.0))

onepointtwo = 1.1 + 0.1
onepointtworounded = with_rounding(RoundDown) do
  1.1 + 0.1
end

println(onepointtworounded)
println(onepointtwo)

println(BigInt("123456789012345678901234567890") + 1)

rational = 1//2
println(rational, " ", typeof(rational))

# "end" keyword?

str = "Hello"
println(str[end])

# Regular expressions

r = r"^\s*(?:#|$)"
println(typeof(r))
m = ismatch(r, "not a comment")
println("match:$m")

m = match(r"(a|b)(c)?(d)", "ad")
println(m)

first, second, third = m.captures
println(first, ",", second, ",", third)

# Byte array literal
println(typeof(b"DATA\xff\u2200"))

# Functions

function f1(x, y)
  x + y
end

f2(x,y) = x + y

println(f1(3, 4))
println(f2(5, 6))
println(apply(f1, 5, 6))

# Operators are functions
println(+(1, 2))

# anonymous function
println(map(x -> x^x, [1, 2, 3, 4]))

# Multiple return values
function tuplize(x)
  x, x
end

tx = tuplize(3)
println(tx)
println(typeof(tx))

# varargs

function vararghello(a, y, x...)
  println(typeof(x))
  println(x)
end

vararghello(1, 2)
vararghello(1, 2, 3)
vararghello(1, 2, 3, 4)

# Keyword args
function kwfunc(x; foo="bar")
  println(foo)
end

kwfunc(1)
kwfunc(2, foo="BAR")
println(kwfunc)

# Block syntax

sq = map([1, 2, 3]) do x
  x*x
end

println(sq)

println(begin
  1
  2
  3
end)

# Loop and range
for i = 1:3
  println("i:$i")
end

for i = [4,5,6]
  println("i:$i")
end

# Exception

try
  sqrt(-1)
catch e
  println("catch:$e")
finally
  println("finally.")
end

function producer()
  produce("start")
  for n=1:4
    produce(2n)
  end
  produce("stop")
end

p = Task(producer)

while true
  c = consume(p)
  println(c)
  if c == "stop"
    break
  end
end

# This works How iteration is modeled?
for i = Task(producer)
  println(i)
end

funcs = map([1,2,3]) do i
  let value = i*2
    () -> println("Value is: $i, $value")
  end
end

for f = funcs; f(); end

# Types

type Foo
  bar
  baz::Int
  qux::Float64
end

foo = Foo(1, 2, 3.0)
println(foo, " ", typeof(foo))

immutable type ImmFoo
  bar
  baz::Int
  qux::Float64
end

println(ImmFoo(1, 2, 3.0))

IntOrString = Union(Int, String)

# ios::IntOrString # Doesn't work. Why?
isstring = "Hello"::IntOrString
isint = 1::IntOrString

abstract Pointy{T}

type Point{T <: Real} <: Pointy{T}
  x::T
  y::T
end

println(typeof(Point(1,2)))
println(Point{Int64} <: Point)

println(isa(1, Integer))

inttype = Int32
while inttype != Any
  println(inttype)
  inttype = super(inttype)
end

#println(methods(+))
methods(+)

type OrderedPair
  x::Real
  y::Real

  # Inner constructor
  OrderedPair(x,y) = x > y ? error("out of order") : new(x,y)
end

# Outer constructor
OrderedPair() = OrderedPair(0, 0)

op0 = OrderedPair()
println(op0)

op1 = OrderedPair(1, 2)
println(op1)

try
  op2 =  OrderedPair(2, 1)
  println(op2)
catch e
  println(e)
end

type ParametricP{T<:Real}
  x::T
  y::T
end

println(ParametricP(1, 2))
println(ParametricP{Float64}(1.0, 2.0))

# Type promotion - Doesn't work. Why?
# Seems like buggy or error prone. https://github.com/JuliaLang/julia/issues/3571
# Even though the promotion concept seems one of the core of this language.
#  http://docs.julialang.org/en/release-0.2/manual/conversion-and-promotion/
#ParametricT(x::Real, y::Real) = ParametricT(promote(x, y)...)
#println(ParametricT(3, 4.0))

module MyModule
modVar = 123
end

println(MyModule.modVar)

# Metaprogramming

exp = :(1 + (2 * Rational(3, 4)))
println(exp, typeof(exp))
println(eval(exp))

sym = :x
println(sym, typeof(sym))

for op = (:+, :*, :&, :|, :$)
  #eval(:(($op)(a,b,c) = ($op)(($op)(a,b),c)))
  @eval ($op)(a,b,c) = ($op)(($op)(a,b),c)
end

macro printvalnexp(ex)
  :(println($(string(ex)), "->", $ex))
end

@printvalnexp 1 + 2 + 3

macro zerox()
  esc(:(x = 0))
end

function usezerox()
  x = 1
  @zerox
  println(x)  # is zero
end

usezerox()

# List comprehension
const x = rand(8)
noise = [ 0.25*x[i-1] + 0.5*x[i] + 0.25*x[i+1] for i=2:length(x)-1 ]
println(noise)

# broadcast
a = rand(2,1)
b = rand(1,2)
c = broadcast(+, a, b)
println(c)
println(typeof(c))

# IO

write(STDOUT, "Hello")

cmd = `echo hello`
println(run(cmd), " ", typeof(cmd))

# Command pipelining
run(`echo world` & `echo hello` |> `sort -r`)

max_stars_repo_name: omo/hello

max_issues_repo_name: omo/hello

max_forks_repo_name: omo/hello


text:
### Global surface temperature map

for i = 1:length(casenames)

    fig, ax = plt[:subplots](1, 1, figsize=(12,8))

#=
    m = basemap.Basemap(projection="merc",llcrnrlat=-80,urcrnrlat=80,
                llcrnrlon=0,urcrnrlon=360,lat_ts=20,resolution="c",lon_0=180.0, ax=ax)

    lab_lons = collect(0:30:360)
    lab_lats = collect(-90:30:90)

    m[:drawcoastlines](linewidth=1.5)
    m[:fillcontinents](color="#888888", lake_color="aqua")
    m[:drawmeridians](lab_lons, labels=repeat([true, ], outer=(length(lab_lons),)))
    m[:drawparallels](lab_lats, labels=repeat([true, ], outer=(length(lab_lats),)))
=#
    casename = casenames[i]
    cmap = plt[:get_cmap]("bwr")
    cmap[:set_over]("r")
    cmap[:set_under]("b")
    cntrs = range(-30, step=5, stop=30)
    Dataset("$nc_file_dir/$casename.h1.nc", "r") do ds

        TREFHT = mean(mreplace(ds["TREFHT"][:, :, t_rng]), dims=(3,))[:, :, 1] .- 273.15
        #plt[:contourf](lon, lat, mask'[:,:], [.45,.55], cmap=cmap, antialiased=false, extend="both")
        #plt[:contourf](lon, lat, mask'[:,:], [.45,.55], cmap=cmap, antialiased=false, extend="both")
        cmapping = plt[:contourf](lon, lat, TREFHT'[:,:], cntrs, cmap=cmap, antialised=true, extend="both")
        plt[:contour](lon, lat, mask'[:,:], [.5], colors="#000000")
        #cs = plt[:contour](lon, lat, TREFHT'[:,:], cntrs, colors="k")
        #plt[:clabel](cs, fmt="%d")
        cba = plt[:colorbar](cmapping, ax=ax)

        ax[:set_title](format("Mean surface temperature ({:d} years): $casename", years))
        ax[:set_xlabel]("Lat [deg]")
        ax[:set_ylabel]("Lon [deg]")
        cba[:set_label]("Temperature [degC]")
    end

    fig[:savefig]("$(casename)_map_surface_temperature.png", dpi=200)

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
function modes(values)
    dict = Dict() # Values => Number of repetitions
    modesArray = typeof(values[1])[] # Array of the modes so far
    max = 0 # Max of repetitions so far

    for v in values
        # Add one to the dict[v] entry (create one if none)
        if v in keys(dict)
            dict[v] += 1
        else
            dict[v] = 1
        end

        # Update modesArray if the number of repetitions
        # of v reaches or surpasses the max value
        if dict[v] >= max
            if dict[v] > max
                empty!(modesArray)
                max += 1
            end
            append!(modesArray, [v])
        end
    end

    return modesArray
end

println(modes([1,3,6,6,6,6,7,7,12,12,17]))
println(modes((1,1,2,4,4)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Advent of Code 2017 - day 4
println("Advent of Code 2017, day 4")

f = open("input4")
input = readlines(f)
close(f)

# Part 1

function no_duplicates(l)
  length(l) == length(unique(l))
end

valid1 = length([l for l in input if no_duplicates(split(l))])

println("part 1: ", valid1)

# Part 2

function no_anagram(l)
  no_duplicates(map(w -> sort(collect(w)), l))
end

valid2 = length([l for l in input if no_anagram(split(l))])

println("part 2: ", valid2)

max_stars_repo_name: venanzio/AoC

max_issues_repo_name: venanzio/AoC

max_forks_repo_name: venanzio/AoC


text:
function translational_symmetry_table(dims)
    nv = prod(dims)
    lat = reshape(collect(1:nv), Tuple(dims))

    iso_table = Vector{Vector{Int}}()
    if length(dims) == 1
        for i=0:dims[1]-1
            push!(iso_table, vec(circshift(lat, i)))
        end
    elseif length(dims) == 2
        for i=0:dims[1]-1
            for j=0:dims[2]-1
                push!(iso_table, vec(circshift(lat, [i,j])))
            end
        end
    elseif length(dims) == 3
        for i=0:dims[1]-1
            for j=0:dims[2]-1
                for k=0:dims[3]-1
                    push!(iso_table, vec(circshift(lat, [i,j,k])))
                end
            end
        end
    else
        @error "not supported"
    end
    iso_table
end

function all_isomorph(graph)
    isomorphs = LightGraphs.Experimental.all_isomorph(graph, graph)

    iso_table = Vector{Vector{Int}}()
    for iso=isomorphs
        map_table = zeros(Int, nv(graph))
        for (src,dst)=iso
            map_table[src] = dst
        end
        push!(iso_table, map_table)
    end
    iso_table
end


text:
Base.@kwdef mutable struct SyntheticCalibrationDataType
    BMRB_setup_label = "BMRB-700" # for use with the NMRData package. See that package's readme for details.

    save_dir = "/home/roy/MEGAsync/outputs/NMR/combined/"
    project_name = "cal01"

    target_names = Vector{String}(undef, 0)

    # mixing weights for each experiment.
    weights = ones(Float64, 0)

    solvent_ppm_guess = 4.7 # guess where the solvent ppm is for each BMRB experiment.
    solvent_window_ppm = 0.1 # allowed solvent shift, in ppm.

end

# taken from an example script in the NMRDataSetup repo.
function setupsyntheticcalibrationdata(config::SyntheticCalibrationDataType)

    ## get BMRB experiment paths from the NMRData package. You can specify your own array of experiments to mix.
    paths = Vector{String}(undef, 0)
    for i = 1:length(config.target_names)
        push!( paths, NMRData.getBMRBexperimentpath(config.BMRB_setup_label, config.target_names[i])[1] )
    end
    @assert all(isdir.(paths)) # check for bad paths.

    ## mix.
    s_t, ν_0ppm, fs, SW,
    s_t_set, ν_0ppm_set, fs_set, SW_set = NMRDataSetup.assemblemixture(paths; weights = config.weights)

    ## store.
    save_folder_path = joinpath(config.save_dir, config.project_name)
    isdir(save_folder_path) || mkdir(save_folder_path); # make save folder if it doesn't exist.

    save_path = joinpath(save_folder_path, "$(config.project_name).bson")
    BSON.bson(save_path,
    s_t = s_t,
    molecule_names = config.target_names,
    fs = fs,
    SW = SW,
    ν_0ppm = ν_0ppm,
    s_t_set = s_t_set,
    ν_0ppm_set = ν_0ppm_set,
    fs_set = fs_set,
    SW_set = SW_set)

    return s_t, ν_0ppm, fs, SW,
        s_t_set, ν_0ppm_set, fs_set, SW_set
end


text:
"""
Useful for a basis for the search and test subspace.
"""
mutable struct SubSpace{matT <: StridedMatrix, matViewT <: StridedMatrix, vecViewT <: StridedVector}
    basis::matT
    all::matViewT
    curr::vecViewT
    prev::matViewT
end

SubSpace(V::AbstractMatrix) = SubSpace(V, view(V, :, 1 : 1), view(V, :, 1), view(V, :, 1 : 0))

function resize!(V::SubSpace, size::Int)
    V.prev = view(V.basis, :, 1 : size - 1)
    V.all = view(V.basis, :, 1 : size)
    V.curr = view(V.basis, :, size)

    V
end

"""
Holds a small projected matrix = W'AV and a view `curr` to the currently active part.
"""
mutable struct ProjectedMatrix{matT <: StridedMatrix, matViewT <: StridedMatrix}
    matrix::matT
    curr::matViewT
end

ProjectedMatrix(M::StridedMatrix) = ProjectedMatrix(M, view(M, 1 : 0, 1 : 0))

function resize!(M::ProjectedMatrix, size::Int)
    M.curr = view(M.matrix, 1 : size, 1 : size)
end

"""
`Q` is the pre-allocated work space of Schur vector
`R` will be the upper triangular factor

`Q` has the form [q₁, ..., qₖ, qₖ₊₁, qₖ₊₂, ..., q_pairs] where q₁ ... qₖ are already converged
and thus locked, while qₖ₊₁ is the active Schur vector that is converging, and the remaining
columns are just garbage data.

It is very useful to have views for these columns:
`locked` is the matrix view [q₁, ..., qₖ]
`active` is the column vector qₖ₊₁
`all` is the matrix view of the non-garbage part [q₁, ..., qₖ, qₖ₊₁]
`locked` is the number of locked Schur vectors
"""
mutable struct PartialSchur{matT <: StridedMatrix, vecT <: AbstractVector, matViewT <: StridedMatrix, vecViewT <: StridedVector}
    Q::matT
    values::vecT

    locked::matViewT
    active::vecViewT
    all::matViewT

    num_locked::Int
end

PartialSchur(Q, numT) = PartialSchur(
    Q,
    numT[],
    view(Q, :, 1 : 0), # Empty view initially
    view(Q, :, 1),
    view(Q, :, 1 : 1),
    0
)

function lock!(pschur::PartialSchur)
    pschur.num_locked += 1
    pschur.locked = view(pschur.Q, :, 1 : pschur.num_locked)

    # Don't extend beyond the max number of Schur vectors
    if pschur.num_locked < size(pschur.Q, 2)
        pschur.all = view(pschur.Q, :, 1 : pschur.num_locked + 1)
        pschur.active = view(pschur.Q, :, pschur.num_locked + 1)
    end

    pschur
end

mutable struct PartialGeneralizedSchur{subT <: SubSpace, vecT <: AbstractVector}
    Q::subT
    Z::subT
    alphas::vecT
    betas::vecT
end

PartialGeneralizedSchur(Q, Z, numT) = PartialGeneralizedSchur(
    SubSpace(Q),
    SubSpace(Z),
    numT[],
    numT[]
)

@inline function resize!(pschur::PartialGeneralizedSchur, size::Int)
    resize!(pschur.Q, size)
    resize!(pschur.Z, size)
end

function shrink!(temporary, subspace::SubSpace, combination::StridedMatrix, dimension)
    tmp = view(temporary, :, 1 : dimension)
    mul!(tmp, subspace.all, combination)
    copyto!(subspace.all, tmp)
    resize!(subspace, dimension)
end

function shrink!(M::ProjectedMatrix, replacement, dimension)
    copyto!(view(M.matrix, 1 : dimension, 1 : dimension), replacement)
    resize!(M, dimension)
end

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
import Logging
using LoggingExtras, Dates

"""
    initialise_logging(true, "/logs/", 1000, dev)

Initialise logging object.

# Examples
```julia-repl
julia> initialise_logging()
julia> @warn("It is bad")
┌ Warning: 2021-06-21 11:47:44 It is bad
└ @ Main REPL[166]:1

julia> @info("normal stuff")
[ Info: 2021-06-21 11:47:44 normal stuff

julia> @error("THE WORSE THING")
┌ Error: 2021-06-21 11:47:44 THE WORSE THING
└ @ Main REPL[168]:1

julia> @debug("it is chill")
┌ Debug: 2021-06-21 11:47:46 it is chill
└ @ Main REPL[169]:1
```
"""
function initialise_logging(log_to_file::Bool, log_path::String, log_level::Int32, app_env::String)
    date_format = "yyyy-mm-dd HH:MM:SS"
    log_level = log_level |> Logging.LogLevel

    logger =  if log_to_file
                isdir(log_path) || mkpath(log_path)
                LoggingExtras.TeeLogger(
                    LoggingExtras.FileLogger(joinpath(log_path, "$(app_env)-$(Dates.today()).log"), always_flush = true, append = true),
                    Logging.ConsoleLogger(stdout, log_level)
                )
                else
                Logging.ConsoleLogger(stdout, log_level)
                end

    timestamp_logger(logger) = LoggingExtras.TransformerLogger(logger) do log
        merge(log, (; message = "$(Dates.format(Dates.now(), date_format)) $(log.message)"))
    end

    LoggingExtras.TeeLogger(LoggingExtras.MinLevelLogger(logger, log_level)) |> timestamp_logger |> global_logger

    nothing
end


text:
# Old idea: use numerical solver for stochastic PDE
# with periodic BCs as a proposal.
# New idea: use periodic GP kernel.
include("Utils.jl")

@doc """
A class defining a trajectory in 1-D space
with a GP prior distribution and periodic boundary
conditions.
""" ->
type Track
	# Number of points in the trajectory
	num_points::Int64

	# Hyperparameter
	length_scale::Float64

	# Covariance matrix and its cholesky decomposition
	C::Array{Float64, 2}
	L::Array{Float64, 2}

	# The trajectory itself
	y::Array{Float64, 1}
end

@doc """
Constructor. Input: num_dimensions, num_points
""" ->
function Track(num_points::Int64)
	return Track(num_points, 0.3,
						Array(Float64, (num_points, num_points)),
						Array(Float64, (num_points, num_points)),
						Array(Float64, (num_points, )))
end

@doc """
Calculate covariance matrix
""" ->
function calculate_C!(track::Track)
	pinv = 1.0/track.num_points
	l2inv = 1.0/track.length_scale^2

	# Fill covariance matrix
	for(j in 1:track.num_points)
		for(i in 1:track.num_points)
			track.C[i, j] = exp(-2*sin(pi*abs(i-j)*pinv)^2*l2inv)
			if(i == j)
				# Small term for positive definiteness
				track.C[i, j] += 1E-6
			else
				track.C[j, i] = track.C[i, j]
			end
		end
	end

	# Do Cholesky decomposition
	track.L = chol(track.C)'

	return nothing
end

@doc """
Generate realisation of trajectory
""" ->
function generate_y!(track::Track)
	# Generate realisation
	track.y = track.L*randn(track.num_points)
	return nothing
end

@doc """
Metropolis proposal
""" ->
function perturb!(track::Track; theta=nothing)
	if(rand() <= 0.5)
		# Make backup of old track
		y_old = copy(track.y)

		# Generate totally new track
		generate_y!(track)

		# Step size
		if(theta == nothing)
			theta = 2*pi*10.0^(-6.0*rand())
		end
		cos_theta = cos(theta)
		sin_theta = sin(theta)

		# AR(1) proposal
		for(i in 1:track.num_points)
			track.y[i] = cos_theta*y_old[i] + sin_theta*track.y[i]
		end
	else
		track.length_scale += 2*randh()
		track.length_scale = mod(track.length_scale, 2.0)
		n = \(track.L, track.y)
		calculate_C!(track)
		track.y = track.L*n
	end

	return 0.0
end

max_stars_repo_name: eggplantbren/MogTrack

max_issues_repo_name: eggplantbren/MogTrack

max_forks_repo_name: eggplantbren/MogTrack


text:
struct GBWeightedEdge{W} <: Graphs.AbstractEdge{UInt64}
    src::UInt64
    dst::UInt64
    weight::W
end

# Not sure about the default to Float64 here, but it's probably the least bad.
GBWeightedEdge(t::NTuple{2}) = GBWeightedEdge(t..., one(Float64))
GBWeightedEdge(t::NTuple{3}) = GBWeightedEdge(t...)
GBWeightedEdge(p::Pair) = GBWeightedEdge(p.first, p.second, one(Float64))
GBWeightedEdge{W}(p::Pair) where W = GBWeightedEdge(UInt64(p.first), UInt64(p.second), one(W))
GBWeightedEdge{W}(t::NTuple{3}) where W = GBWeightedEdge(UInt64(t[1]), UInt64(t[2]), W(t[3]))
GBWeightedEdge{W}(t::NTuple{2}) where W = GBWeightedEdge(UInt64(t[1]), UInt64(t[2]), one(W))
GBWeightedEdge(x, y) = GBWeightedEdge(x, y, one(Float64))

Graphs.src(e::GBWeightedEdge) = e.src
Graphs.dst(e::GBWeightedEdge) = e.dst
weight(e::GBWeightedEdge) = e.weight

#Printing
Base.show(io::IO, e::GBWeightedEdge) = print(io, "Edge $(e.src) => $(e.dst) with weight $(e.weight)")

Base.Tuple(e::GBWeightedEdge) = (src(e), dst(e), weight(e))

Graphs.reverse(e::GBWeightedEdge) = GBWeightedEdge(dst(e), src(e), weight(e))

==(e1::GBWeightedEdge, e2::GBWeightedEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2))
==(e1::GBWeightedEdge, e2::AbstractEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2))
==(e1::AbstractEdge, e2::GBWeightedEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2))


text:
export MinimalDominatingSet

struct MinimalDominatingSet end

"""
    dominating_set(g, MinimalDominatingSet(); seed=-1)

Find a set of vertices that consitute a dominating set (all vertices in `g` are either adjacent to a vertex
in the set or is a vertex in the set) and it is not possible to delete a vertex from the set
without sacrificing the dominating property.

### Implementation Notes
Initially, every vertex is in the dominating set.
In some random order, we check if the removal of a vertex from the set will destroy the
dominating property. If no, the vertex is removed from the dominating set.

### Performance
Runtime: ``\\mathcal{O}(|V|+|E|)``
Memory: ``\\mathcal{O}(|V|)``

### Optional Arguments
- If `seed >= 0`, a random generator is seeded with this value.
"""
function dominating_set(
    g::AbstractGraph{T},
    alg::MinimalDominatingSet;
    seed::Int=-1
    ) where T <: Integer

    nvg = nv(g)
    in_dom_set = trues(nvg)
    length_ds = Int(nvg)
    dom_degree = degree(g)
    @inbounds @simd for v in vertices(g)
        dom_degree[v] -= (has_edge(g, v, v) ? 1 : 0)
    end

    for v in randperm(getRNG(seed), nvg)
    	(dom_degree[v] == 0) && continue #It is not adjacent to any dominating vertex
    	#Check if any vertex is depending on v to be dominated
        dependent = findfirst(u -> !in_dom_set[u] && dom_degree[u] <= 1, neighbors(g, v))

        (dependent != nothing) && continue
        in_dom_set[v] = false
        length_ds -= 1
        dom_degree[neighbors(g, v)] .-= 1
    end

    return LightGraphs.findall!(in_dom_set, Vector{T}(undef, length_ds))
end


text:
#
# oipost.jl --
#
# Post-processing of OI-FITS data-block definitions.  Must be the last
# file included by the main OIFITS.jl source file.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015-2019: Éric Thiébaut.
#
#------------------------------------------------------------------------------

# Automatically define getters from all fields of a data-block.
for (dbname, dbtype) in _DATABLOCKS
    for symb in _FIELDS[dbname]
        eval(Meta.parse("get_$symb(db::$dbtype) = db.contents[:$symb]"))
    end
end

# Define getters which rely on indirections.
get_eff_wave(db::Union{OIVis,OIVis2,OIT3}) = db.ins[:eff_wave]
get_eff_band(db::Union{OIVis,OIVis2,OIT3}) = db.ins[:eff_band]

@doc """
Assuming `db` is an instance of `OITarget`, `OIVis`, `OIVis2` or `OIT3`, then:
```
    get_target_id(db)
```
yields the "TARGET_ID" column of `db` which is an array of integers.
""" get_target_id