text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.TurbulenceClosures: VerticallyImplicit, Vertical, Horizontal
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VerticalVorticityField
using Printf
using GLMakie
using JLD2

grid = LatitudeLongitudeGrid(size = (64, 1024, 4),
                             longitude = (-180, 180),
                             latitude = (30, 50),
                             z = (-1kilometers, 0),
                             halo = (3, 3, 3))

const Lz = grid.Lz

# Uncomment to put a bump in the grid:
#=
const width = 10 # degrees
bump(λ, φ) = - Lz * (1 - 0.5 * exp(-λ^2 / 2width^2))
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bump))
=#

diffusive_closure = ScalarDiffusivity(ν = 1e-2, κ = 1e-2,
                                      isotropy = Vertical(),
                                      time_discretization = VerticallyImplicit())

horizontal_closure = ScalarBiharmonicDiffusivity(ν=1e6, κ=1e6, isotropy = Horizontal())

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    coriolis = HydrostaticSphericalCoriolis(),
                                    buoyancy = BuoyancyTracer(),
                                    closure = (diffusive_closure, horizontal_closure),
                                    tracers = (:b, :c),
                                    momentum_advection = VectorInvariant(),
                                    tracer_advection = WENO5(),
                                    free_surface = ImplicitFreeSurface())

# Initial condition: a baroclinically unstable situation!
ramp(φ, δφ, φ₀) = min(max(0, (φ - φ₀)/δφ + 1/2), 1)

# Parameters
N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

φ₀ = 40 # degrees
δφ = 1 # degrees
δz = 400 # m
δc = 2δφ
δb = δφ * M²
ϵb = 1e-2 * δb # noise amplitude

bᵢ(λ, φ, z) = N² * z + δb * ramp(φ, δφ, φ₀) + ϵb * randn()
cᵢ(λ, φ, z) = exp(-(φ - φ₀)^2 / 2δc^2) * exp(-(z + Lz/4)^2 / 2δz^2)

set!(model, b=bᵢ, c=cᵢ)

simulation = Simulation(model, Δt=20minutes, stop_time=30days)

print_progress(sim) = @printf("i: %d, t: %s, wall time: %s, max(u): %6.3e m s⁻¹ \n",
                              iteration(sim), prettytime(sim), prettytime(sim.run_wall_time),
                              maximum(abs, sim.model.velocities.u))

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

b, c = model.tracers
u, v, w = model.velocities
ζ = VerticalVorticityField(model)

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ζ, b, c),
                                                      schedule = TimeInterval(1hour),
                                                      field_slicer = FieldSlicer(k=grid.Nz),
                                                      prefix = "baroclinic_adjustment",
                                                      force = true)

run!(simulation)

file = jldopen("baroclinic_adjustment.jld2")
iterations = parse.(Int, keys(file["timeseries/t"]))
ζ = [file["timeseries/ζ/$i"][:, :, 1] for i in iterations]
b = [file["timeseries/b/$i"][:, :, 1] for i in iterations]
c = [file["timeseries/c/$i"][:, :, 1] for i in iterations]
t = [file["timeseries/t/$i"] for i in iterations]
close(file)

Nt = length(t)

fig = Figure()

axζ = Axis(fig[1, 1], aspect=:data)
axb = Axis(fig[2, 1], aspect=:data)
axc = Axis(fig[3, 1], aspect=:data)

slider = Slider(fig[2, :], range=1:Nt, startvalue=1)
n = slider.value

xζ, yζ, zζ = nodes((Face, Face, Center), grid)
xc, yc, zc = nodes((Center, Center, Center), grid)

ζⁿ = @lift ζ[$n]
bⁿ = @lift b[$n]
cⁿ = @lift c[$n]

hmζ = heatmap!(axζ, xζ, yζ, ζⁿ)
hmb = heatmap!(axb, xc, yc, bⁿ)
hmc = heatmap!(axc, xc, yc, cⁿ)

# Colorbar(fig[2, 1], xζ, yζ, hmζ, horizontal=true)
# Colorbar(fig[2, 2], xc, yc, hmb, horizontal=true)
# Colorbar(fig[2, 3], xc, yc, hmc, horizontal=true)

display(fig)

ext: jl


text:
export choice_probabilities, choice_probabilities!

"""
Compute the choice probabilities using an overflow safe algorithm.

This is just a wrapper around [`choice_probabilities!`](@ref) which computes them in-place.
"""
function choice_probabilities(utilities::AbstractMatrix{<:AbstractFloat})
    probabilities = similar(utilities)
    choice_probabilities!(probabilities, utilities)
    return probabilities
end

"""
Compute the choice probabilities in-place using an overflow safe algorithm.

Based on the SoftmaxThreePassReload function from
Dukhan and Ablavatski (2020), The Two-Pass Softmax Algorithm (https://arxiv.org/abs/2001.04438).
"""
function choice_probabilities!(output::AbstractMatrix{T}, utilities::AbstractMatrix{T}) where T<:AbstractFloat
    @assert size(output) == size(utilities) "Output matrix must have same size as utilities matrix."
    J, I = size(utilities)

    @inbounds for i = 1:I
        max_u = T(0)
        for j = 1:J
            if utilities[j, i] > max_u
                max_u = utilities[j, i]
            end
        end

        sigma = @fastmath exp(-max_u)

        @simd ivdep for j = 1:J
            output[j, i] = @fastmath exp(utilities[j, i] - max_u)
            # SIMD makes left-to-right summation both faster and more accurate
            # https://discourse.julialang.org/t/when-shouldnt-we-use-simd/18276/14
            sigma += output[j, i]
        end

        sigma_inv = 1 / sigma
        @simd ivdep for j = 1:J
            output[j, i] *= sigma_inv
        end
    end
end


text:
using Revise
using JDF, DataFrames
# a = gf()
# savejdf("c:/data/a.jdf", a)
@time a = loadjdf("c:/data/large.jdf")

# 2G file
# @time a = CSV.read(
#     "C:/data/Performance_All/Performance_2010Q3.txt",
#     delim = '|',
#     header = false
# );

@time b = Array(a[:, :Column2]);

@time res = type_compress(b);

all(b .== res)

# ne = findfirst(b .!= b1)
#
# b[ne], res[ne]

using RLEVectors, StatsBase
ar = rle(res)

ar[2]  .= cumsum(ar[2])
rlev= RLEVector(ar...)

all(b .== rlev)

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
function ingredients(path::String)
	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol(basename(path))
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
import Random
using AbstractAlgebra

function test_params_nf(rng, nvariables, exps, nterms,
                        npolys, grounds, coeffssize, orderings)

    for n in nvariables
        for e in exps
            for nt in nterms
                for np in npolys
                    for gr in grounds
                        for ord in orderings
                            for csz in coeffssize
                                set = Groebner.generate_set(
                                    n, e, nt, np, csz, rng, gr, ord
                                )
                                isempty(set) && continue
                                gb = Groebner.groebner(set)
                                f = rand(set)

                                mynf = Groebner.normalform(gb, f)
                                # TODO : add more tests
                                @test mynf == 0
                                if mynf != 0
                                    @error "Beda!" n e nt np gr ord
                                    println("Rng ", rng)
                                    println("Set ", set)
                                    println("Gb ", gb)
                                    error("exit tests")
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

@testset "normalform random stress tests" begin

    rng = Random.MersenneTwister(42)

    # :degrevlex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 4:5]
    nterms     = [1:1, 1:2, 4:5]
    npolys     = [1:1, 1:3, 4:5]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:degrevlex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for normal form"
    test_params_nf(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :lex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:lex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for normal form"
    test_params_nf(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

    # :deglex finite case tests
    nvariables = [2, 3, 4]
    exps       = [1:2, 2:4, 2:3]
    nterms     = [1:1, 1:2, 2:3]
    npolys     = [1:1, 1:3, 2:3]
    grounds    = [GF(1031), GF(2^31 - 1)]
    coeffssize = [3, 1000, 2^31 - 1]
    orderings  = [:deglex]
    p = prod(map(length, (nvariables, exps, nterms, npolys, grounds, orderings, coeffssize)))
    @info "producing $p $(orderings[1]) tests for normal form"
    test_params_nf(rng, nvariables, exps, nterms, npolys, grounds, coeffssize, orderings)

end


text:
function find_zero_indices(a)
  eachindex(a)'[a .== 0]
end

function find_nonzero_indices(a)
  eachindex(a)'[a .!= 0]
end

"""
`INCIDENCE_TO_INCIDENT_LIST_ARRAY` - takes in M edges on N nodes in NxM sparse format and returns an 3xN array where the N-th column includes the indices of the nodes connected (one-direction) to node N
 * Currently assumes that there's at most 3 points, and furthermore that the nonzeros alternate
 * 0 in the incidence_array means that there's less than 3 points
"""
function incidence_to_incident_list_array(E)
  incident_list_array = zeros(Int64, 3, size(E, 1))
  incident_list_counts = ones(Int64, size(E, 1))
  rows = rowvals(E)
  vals = nonzeros(E)

  ind = 1
  # for each edge
  for j in 1:size(E, 2)
    	node_first = rows[ind]
    	node_second = rows[ind+1]
	incident_list_count = incident_list_counts[node_first]
	incident_list_array[incident_list_count, node_first] = node_second
	incident_list_counts[node_first] = incident_list_count+1
	ind = ind+2
  end
  return incident_list_array
end

# checks if p -> p(pa,pb) is the same as p -> pc -> pca(b,c) and writes the triangle
function check_and_write_triangles!(triangles, cur_triangle, p, pa, pb, pc, pca, pcb, pcc)
  if p == 0 return cur_triangle end
  if pc == 0 return cur_triangle end
  if pa != 0
	if pa == pca || pa == pcb || pa == pcc
	triangles[1, cur_triangle] = p
	triangles[2, cur_triangle] = pc
	triangles[3, cur_triangle] = pa
	cur_triangle += 1
      end
  end
  if pb != 0
	if pb == pca || pb == pcb || pb == pcc
	triangles[1, cur_triangle] = p
	triangles[2, cur_triangle] = pc
	triangles[3, cur_triangle] = pb
	cur_triangle += 1
      end
  end
  return cur_triangle
end
"""
`INCIDENT_LIST_ARRAY_TO_TRIANGLES` - Convert incidence array to 3xN list triangle node index vectors
"""
function incident_list_array_to_triangles(incident_list_array)
  triangles = zeros(Int64, 3, size(incident_list_array, 2) * 10)
  max_edges = size(incident_list_array, 1)
  cur_triangle = 1
  # absolutely awful, but really fast
  # for each node at j
  for p in 1:size(incident_list_array, 2)
    # the up to 3 points that have edges from j
    	pa = incident_list_array[1, p]
    	pb = incident_list_array[2, p]
    	pc = incident_list_array[3, p]

      if pa != 0
	# the up to 9 points that have edges from pa, pb, pc
	paa = incident_list_array[1, pa]; pab = incident_list_array[2, pa]; pac = incident_list_array[3, pa]
      else
	paa = 0; pab = 0; pac = 0
      end

      if pb != 0
	pba = incident_list_array[1, pb]; pbb = incident_list_array[2, pb]; pbc = incident_list_array[3, pb]
      else
	pba = 0; pbb = 0; pbc = 0
      end

      if pc != 0
	pca = incident_list_array[1, pc]; pcb = incident_list_array[2, pc]; pcc = incident_list_array[3, pc]
      else
	pca = 0; pcb = 0; pcc = 0
      end

	cur_triangle = check_and_write_triangles!(triangles, cur_triangle, p, pa, pb, pc, pca, pcb, pcc)
	cur_triangle = check_and_write_triangles!(triangles, cur_triangle, p, pa, pc, pb, pba, pbb, pbc)
	cur_triangle = check_and_write_triangles!(triangles, cur_triangle, p, pb, pc, pa, paa, pab, pac)
  end
  return triangles[:, 1:(cur_triangle - 1)]
end

"""
`INCIDENCE_TO_TRIANGLES` - Convert node-edge incidence matrix into triangle
array (triples of node indices defining all triangles in the mesh). The triangle
array will be used when applying piecewise affine transforms through the
meshwarp function.

See Mesh type documentation for definition of an incidence matrix.

```
triangles = incidence_to_triangles(E)
```

* E: node-edge incidence matrix
* triangles: 3xN array of ints, each column represents the indices of three nodes
    in one triangle.
"""
function incidence_to_triangles(E)
	incident_list_array = incidence_to_incident_list_array(E)
	return incident_list_array_to_triangles(incident_list_array)
end

#= THESE FUNCTIONS ARE SLOW
"""
`INCIDENCE_TO_DICT` - Create dictionary of node sets connected to indexing node
"""
function incidence_to_dict(D)
  D = abs.(D)
  node_dict = Dict()

  for i = 1:size(D,1)
    j = 1
    while D[i,j] == 0
      j += 1
    end
    if !(j in keys(node_dict))
      node_dict[j] = Set{Int64}()
    end
    k = j+1
    while D[i,k] == 0
      k += 1
    end
    push!(node_dict[j], k)
  end
  return node_dict
end

"""
`DICT_TO_TRIANGLES` - Convert node dict to Nx3 list triangle node index vectors
"""
function dict_to_triangles(node_dict)
  triangles = Array{Int64}(0, 3)
  for a in sort(collect(keys(node_dict)))
    setA = node_dict[a]
    for b in sort(collect(setA))
      if !(b in keys(node_dict))
        continue
      end
      setB = node_dict[b]
      setC = intersect(setA, setB)
      for c in sort(collect(setC))
        triangles = vcat(triangles, [a b c])
      end
    end
  end
  return triangles
end

"""
`INCIDENCE_TO_TRIANGLES` - Convert edge-node incidence matrix into triangle
array (triples of node indices defining all triangles in the mesh). The triangle
array will be used when applying piecewise affine transforms through the
meshwarp function.

See Mesh type documentation for definition of an incidence matrix.

```
triangles = incidence_to_triangles(D)
```

* D: edge-node incidence matrix
* triangles: Nx3 array of ints, each row represents the indices of three nodes
    in one triangle.
"""
function incidence_to_triangles(D)
  node_dict = incidence_to_dict(D)
  return dict_to_triangles(node_dict)
end
=#


text:
@testset "Cross validating ooHMM on CTB..." begin
        @time ctb = ChTreebank()
        #@time ctb = ZhuXian()

        nb = 10
        batches = KongYiji.foldbatch(ctb, nb)
        tbs = KongYiji.HmmScoreTable[]
        @showprogress 1 "Cross Validating ooHMM..." for (tr, te) in batches
                base = KongYiji.HMM(tr)
                hmm = KongYiji.HMM(;base=base)
                lm = LM(tr; rnnType=:lstm, numLayers=1, bidirectional=false)
                oohmm = hmm + lm
                x = KongYiji.rawsents(te)
                z = KongYiji.wordsents(te)
                @time y = oohmm(x; recover=false, withpos=false)
                push!(tbs, KongYiji.HmmScoreTable(z, y))
                println(tbs[end])
        end
        hmmscoretable = sum(tbs)
        println(hmmscoretable)
end


text:
# https://github.com/spencerlyon2/PlotlyJS.jl

const _plotlyjs_attr        = _plotly_attr
const _plotlyjs_seriestype  = _plotly_seriestype
const _plotlyjs_style       = _plotly_style
const _plotlyjs_marker      = _plotly_marker
const _plotlyjs_scale       = _plotly_scale

# --------------------------------------------------------------------------------------

function add_backend_string(::PlotlyJSBackend)
    """
    if !Plots.is_installed("PlotlyJS")
        Pkg.add("PlotlyJS")
    end
    if !Plots.is_installed("Rsvg")
        Pkg.add("Rsvg")
    end
    import Blink
    Blink.AtomShell.install()
    """
end

function _initialize_backend(::PlotlyJSBackend; kw...)
    @eval begin
        import PlotlyJS
        export PlotlyJS
    end

    # # override IJulia inline display
    # if isijulia()
    #     IJulia.display_dict(plt::AbstractPlot{PlotlyJSBackend}) = IJulia.display_dict(plt.o)
    # end
end

# ---------------------------------------------------------------------------

function _create_backend_figure(plt::Plot{PlotlyJSBackend})
    PlotlyJS.plot()
end

function _series_added(plt::Plot{PlotlyJSBackend}, series::Series)
    syncplot = plt.o
    pdicts = plotly_series(plt, series)
    for pdict in pdicts
        typ = pop!(pdict, :type)
        gt = PlotlyJS.GenericTrace(typ; pdict...)
        PlotlyJS.addtraces!(syncplot, gt)
    end
end

function _series_updated(plt::Plot{PlotlyJSBackend}, series::Series)
    xsym, ysym = (ispolar(series) ? (:t,:r) : (:x,:y))
    kw = KW(xsym => (series.d[:x],), ysym => (series.d[:y],))
    z = series[:z]
    if z != nothing
        kw[:z] = (transpose_z(series, series[:z].surf, false),)
    end
    PlotlyJS.restyle!(
        plt.o,
        findfirst(plt.series_list, series),
        kw
    )
end

# ----------------------------------------------------------------

function _update_plot_object(plt::Plot{PlotlyJSBackend})
    pdict = plotly_layout(plt)
    syncplot = plt.o
    w,h = plt[:size]
    PlotlyJS.relayout!(syncplot, pdict, width = w, height = h)
end

# ----------------------------------------------------------------

function _show(io::IO, ::MIME"image/svg+xml", plt::Plot{PlotlyJSBackend})
    show(io, MIME("text/html"), plt.o)
end

function plotlyjs_save_hack(io::IO, plt::Plot{PlotlyJSBackend}, ext::String)
    tmpfn = tempname() * "." * ext
    PlotlyJS.savefig(plt.o, tmpfn)
    write(io, read(open(tmpfn)))
end
_show(io::IO, ::MIME"image/png", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "png")
_show(io::IO, ::MIME"application/pdf", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "pdf")
_show(io::IO, ::MIME"image/eps", plt::Plot{PlotlyJSBackend}) = plotlyjs_save_hack(io, plt, "eps")

function _display(plt::Plot{PlotlyJSBackend})
    display(plt.o)
end


text:
using MOCNeutronTransport
# Reactor-Physics-Benchmarks-Handbook/2015/LWR/CROCUS-LWR-RESR-001
# Model
# ----------------------------------------------------------------------------------------------
gmsh.initialize()
uo2_fuel_entities = Int32[]
um_fuel_entities = Int32[]
gap_entities = Int32[]
clad_entities = Int32[]
absorber_entities = Int32[]
# The external radius of the water is 65 cm
# How much water is needed to count as infinite reflector?
# What is the minimum distance from fuel to boundary?
# Don't want to simulate what we don't have to. Or could make the mesh towards the problem
# boundary really large to minimize computations on a larger problem.
# bb_apothem = 65.0
bb_apothem = 11.2*5
bb = (0.0, 2*bb_apothem, 0.0, 2*bb_apothem)

# UO2 pins
# ----------------------------------------------------------------------------------------------
# Fuel diameter = 10.52 mm (pg. 34)
# External cladding diameter = 12.6 mm (pg. 34)
# Cladding thickness = 0.85 mm (pg. 34)
# The document does not mention the gap thickness, so I believe we are left to compute that quantity
# Gap thickness = (External cladding diameter - 2(Clad thickness) - Fuel diameter)/2
# Gap thickness = (12.6 - 2(0.85) - 10.52)/2 = 0.19 mm
# Therefore
#   r_fuel = 10.52/2 = 5.26 mm = 0.526 cm
#   r_gap = 5.26 + 0.19 mm = 0.545 cm
#   r_clad = 12.6/2 = 6.3 mm = 0.63 cm
r_fuel = 0.526
r_gap = 0.545
r_clad = 0.63
pitch = 1.8370 # pg.34
# Sweet from left to right for each column in the uo2 pins.
# Place pins bottom to top for each row
for i = -11:11
    if i == 0
        continue
    end
    j_max = [11, 11, 11, 9, 9, 9, 6, 6, 6, 3, 3]
    for j = (-j_max[abs(i)] + 1):j_max[abs(i)]
        push!(uo2_fuel_entities,
              gmsh.model.occ.add_disk(i*pitch - sign(i)*pitch/2 + bb_apothem,
                                      j*pitch - pitch/2 + bb_apothem,
                                      0, r_fuel, r_fuel))
        push!(gap_entities,
              gmsh.model.occ.add_disk(i*pitch - sign(i)*pitch/2 + bb_apothem,
                                      j*pitch - pitch/2 + bb_apothem,
                                      0, r_gap,  r_gap))
        push!(clad_entities,
              gmsh.model.occ.add_disk(i*pitch - sign(i)*pitch/2 + bb_apothem,
                                      j*pitch - pitch/2 + bb_apothem,
                                      0, r_clad, r_clad))
    end
end
# U metal pins
# ----------------------------------------------------------------------------------------------
# Fuel diameter = 17.0 mm (pg. 34)
# External cladding diameter = 19.3 mm (pg. 34)
# Cladding thickness = 0.975 mm (pg. 34)
# Gap thickness = (External cladding diameter - 2(Clad thickness) - Fuel diameter)/2
# Gap thickness = (19.3 - 2(0.975) - 17.0)/2 = 0.175 mm
# Therefore
#   r_fuel = 17.0/2 = 8.5 mm = 0.85 cm
#   r_gap = 8.5 + 0.175 mm = 0.8675 cm
#   r_clad = 19.3/2 = 9.65 mm = 0.965 cm
# Note that the distribution of pins is not quarter symmetric.
r_fuel_m = 0.85
r_gap_m = 0.8675
r_clad_m = 0.965
pitch_m = 2.917 # pg.34
for i = -10:10
    if i == 0
        continue
    end
    j_max = [10, 10, 10, 9, 9, 8, 8, 7, 5, 3]
    j_min = [ 8,  8,  7, 7, 5, 5, 3, 1, 1, 1]
    for j = j_min[abs(i)]:j_max[abs(i)]
        push!(uo2_fuel_entities,
              gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      0, r_fuel_m, r_fuel_m))
        push!(gap_entities,
              gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      0, r_gap_m,  r_gap_m))
        push!(clad_entities,
              gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                      0, r_clad_m, r_clad_m))
    end
end
for i = -10:10
    if i == 0
        continue
    end
    j_max = [10, 10, 10, 9, 9, 9, 8, 7, 6, 4]
    j_min = [ 8,  8,  7, 7, 5, 5, 3, 1, 1, 1]
    for j = (j_min[abs(i)] - 1):(j_max[abs(i)] - 1)
        # Empty tube
        if (i == 6 && j == 5) || (i == -6 && j == 5)
            push!(gap_entities,
                  gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          -sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          0, r_gap_m,  r_gap_m))
            push!(clad_entities,
                  gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          -sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          0, r_clad_m, r_clad_m))
        else
            push!(uo2_fuel_entities,
                  gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          -sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          0, r_fuel_m, r_fuel_m))
            push!(gap_entities,
                  gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          -sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          0, r_gap_m,  r_gap_m))
            push!(clad_entities,
                  gmsh.model.occ.add_disk(i*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          -sign(i)*j*pitch_m - sign(i)*pitch_m/2 + bb_apothem,
                                          0, r_clad_m, r_clad_m))
        end
    end
end
# Absorber rod
# ----------------------------------------------------------------------------------------------
# Absorber diameter = 6.0 mm (pg. 32)
# Clad diameter = 8.0 mm (pg. 32)
push!(absorber_entities,
      gmsh.model.occ.add_disk(bb_apothem, bb_apothem, 0, 0.3, 0.3))
push!(clad_entities,
      gmsh.model.occ.add_disk(bb_apothem, bb_apothem, 0, 0.4, 0.4))

# Materials
# ----------------------------------------------------------------------------------------------
gmsh.model.occ.synchronize()
p = gmsh.model.add_physical_group(2, uo2_fuel_entities)
gmsh.model.set_physical_name(2, p, "MATERIAL_UO2")
p = gmsh.model.add_physical_group(2, gap_entities)
gmsh.model.set_physical_name(2, p, "MATERIAL_GAP")
p = gmsh.model.add_physical_group(2, clad_entities)
gmsh.model.set_physical_name(2, p, "MATERIAL_CLAD")
p = gmsh.model.add_physical_group(2, absorber_entities)
gmsh.model.set_physical_name(2, p, "MATERIAL_ABSORBER")
gmsh_group_preserving_fragment(gmsh.model.get_entities(2),
                               gmsh.model.get_entities(2);
                               material_hierarchy = ["MATERIAL_UO2",
                                                     "MATERIAL_ABSORBER",
                                                     "MATERIAL_GAP",
                                                     "MATERIAL_CLAD"])

# Overlay Grid
# ---------------------------------------------------------------------------------------------
# The preferable locations to place a grid division are the horizontal or vertical lines that do
# not intersect any pins. This way, we can potentially capture one pin per coarse cell, and easily
# show pin powers.
#
# If we define the center of the model at (0, 0)
# The valid x or y locations for the UO2 pins in the positive x, y quadrant are:
# [0, pitch/2 - r_clad)
# ( [(2i-1)/2]pitch + r_clad, [(2(i+1)-1)/2]pitch - r_clad)
# (23pitch/2 + r_clad, ∞)
# i = 1, 2, ..., 11
# The valid x or y locations for the U metal pins are
# [0, pitch_m/2 - r_clad_m)
# ( [(2i-1)/2]pitch_m + r_clad_m, [(2(i+1)-1)/2]pitch_m - r_clad_m)
# (21pitch_m/2 + r_clad_m, ∞)
# i = 1, 2, ..., 10
# Not that due to symmetry, -1*interval is also valid.
# We may compute the intersections of these sets.
uo2_sets = [(0.0, pitch/2 - r_clad)]
for i = 1:11
    push!(uo2_sets, ( ((2i-1)/2)pitch + r_clad, ((2(i+1)-1)/2)pitch - r_clad))
end
push!(uo2_sets, (23pitch/2 + r_clad, Inf))

um_sets = [(0.0, pitch_m/2 - r_clad_m)]
for i = 1:10
    push!(um_sets, ( ((2i-1)/2)pitch_m + r_clad_m, ((2(i+1)-1)/2)pitch_m - r_clad_m))
end
push!(um_sets, (21pitch_m/2 + r_clad_m, Inf))

set_intersections = Tuple{Float64, Float64}[]
for (uo2_min, uo2_max) in uo2_sets
    for (um_min, um_max) in um_sets
        # um set intersects on the left
        if um_min ≤ uo2_min ≤ min(um_max, uo2_max)
            push!(set_intersections, (uo2_min, min(um_max, uo2_max)))
        # um set intersects on the right
        elseif um_min ≤ uo2_max ≤ um_max
            push!(set_intersections, (max(um_min, uo2_min), uo2_max))
        end
    end
end
# Visualize the intervals in 1d
#for (ymin, ymax) in set_intersections
#    if ymax == Inf
#        ymax = bb_apothem
#    end
#    gmsh.model.occ.add_rectangle(bb_apothem, ymin + bb_apothem, 0.0, 30.0, ymax - ymin)
#end
#gmsh.model.occ.synchronize()
#
# Note that we may also use intervals of length 2*value, and divide the problem using an odd
# number of divisions
#
# Possible RT modules ± ϵ. Check the intervals for full range of valid values.
# 1 x 1 - (-65, 65)
# 2 x 2 - (-65, 0, 65)
# 3 x 3 - (-60, -20,  20, 60)
# 4 x 4 - (-65, -32.5, 0, 32.5, 65)
# 5 x 5 - (-56, -33.6, -11.2, 11.2, 33.6, 56) & (-73, -43.8, -14.6, 14.6, 43.8, 73)
# 6 x 6 - None
# 7 x 7 - None
# 8 x 8 - None
# ... None
nx = [5]
ny = [5]
grid_tags = gmsh_overlay_rectangular_grid(bb, "MATERIAL_MODERATOR", nx, ny)
gmsh.fltk.run()

gmsh.finalize()

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
# https://docs.julialang.org/en/v1/base/collections/#Base.Pair
karty = [Pair(7,"7"),Pair(8,"8"),Pair(9,"9"),Pair(10,"10"),Pair(11,"dolník"),Pair(12,"horník"),Pair(13,"kráľ"),Pair(14,"eso")]
farby_kariet = ["zeleň", "žaluď", "červeň", "guľa"]

# https://discourse.julialang.org/t/cleanest-way-to-generate-all-combinations-of-n-arrays/20127
# help?> Iterators.product
vsetky_karty = collect(Iterators.product(karty, farby_kariet))

najhodnotnejsia_karta = vsetky_karty[1,1:4]

max_stars_repo_name: LRjunior/kartova-hra-sedma

max_issues_repo_name: LRjunior/kartova-hra-sedma

max_forks_repo_name: LRjunior/kartova-hra-sedma


text:
using AuditoryStimuli
using Test
using DSP
using StatsBase
using Statistics
using Logging
using Plots
using Unitful
using SampledSignals
using Pipe
using Images

Fs = 48000

ENV["JULIA_DEBUG"] = "all"

@testset "Offline Stimuli" begin

    @testset "Generator Functions" begin
    # ==================================

        @testset "One hit signal generation" begin

            @testset "Bandpass Noise" begin

                # Test different constructors
                bn = bandpass_noise(Fs * 30, 2, 300, 700, Fs)

                # Test data is actuall filtered
                for lower_bound = 500:500:1500
                    for upper_bound = 2000:500:3000

                        bn = bandpass_noise(Fs * 30, 2, lower_bound, upper_bound, Fs)
                        @test size(bn, 1) == Fs * 30
                        spec = welch_pgram(bn[:, 1], fs=Fs)

                        val, idx_lb = findmin(abs.(freq(spec) .- lower_bound))
                        val, idx_bl = findmin(abs.(freq(spec) .- (lower_bound - 250)))
                        @test (amp2db(power(spec)[idx_lb]) - amp2db(power(spec)[idx_bl])) > 10

                        val, idx_ub = findmin(abs.(freq(spec) .- upper_bound))
                        val, idx_bu = findmin(abs.(freq(spec) .- (upper_bound + 250)))
                        @test (amp2db(power(spec)[idx_ub]) - amp2db(power(spec)[idx_bu])) > 10
                    end
                end
            end
        end

    end

    @testset "Modifier Functions" begin
    # ==================================

        @testset "Filter Signals" begin

            @testset "Bandpass Butterworth" begin

                @testset "Abstract Arrays" begin

                    # Test data is actuall filtered
                    for lower_bound = 500:500:1500
                        for upper_bound = 2000:500:3000

                            x = randn(Fs*30, 2)
                            bn = bandpass_filter(x, lower_bound, upper_bound, Fs)

                            @test size(bn, 1) == Fs * 30
                            @test size(bn, 2) == 2

                            for channel = 1:2
                                spec = welch_pgram(bn[:, channel], fs=Fs)

                                val, idx_lb = findmin(abs.(freq(spec) .- lower_bound))
                                val, idx_bl = findmin(abs.(freq(spec) .- (lower_bound - 250)))
                                @test (amp2db(power(spec)[idx_lb]) - amp2db(power(spec)[idx_bl])) > 10

                                val, idx_ub = findmin(abs.(freq(spec) .- upper_bound))
                                val, idx_bu = findmin(abs.(freq(spec) .- (upper_bound + 250)))
                                @test (amp2db(power(spec)[idx_ub]) - amp2db(power(spec)[idx_bu])) > 10
                            end
                        end
                    end
                end

                @testset "Sampled Signals" begin

                    source = CorrelatedNoiseSource(Float64, 48000, 2, 1, 0.1)
                    a = read(source, 48000)
                    b = bandpass_filter(a, 300u"Hz", 700u"Hz")

                    @test typeof(b) == typeof(a)
                    @test typeof(b) == SampleBuf{Float64,2}

                end

            end

        end

        @testset "Modulate Signals" begin

            @testset "Amplitude Modulation" begin

                @testset "Abstract Arrays" begin

                    for modulation_frequency = 1:1:10
                        x = randn(Fs, 1)
                        @test_nowarn amplitude_modulate(x, modulation_frequency, Fs)
                    end

                    for modulation_frequency = 1.3:1:10
                        x = randn(Fs, 1)
                        amplitude_modulate(x, modulation_frequency, Fs)
                        @test_logs (:warn, "Not a complete modulation") amplitude_modulate(x, modulation_frequency, Fs)
                    end
                end

                @testset "Sampled Signals" begin

                    source = CorrelatedNoiseSource(Float64, 48000, 2, 1, 0.1)
                    a = read(source, 48000)
                    b = amplitude_modulate(a, 20u"Hz")

                    @test typeof(b) == typeof(a)
                    @test typeof(b) == SampleBuf{Float64,2}

                end

            end

            @testset "ITD Modulation" begin

                @testset "Abstract Arrays" begin

                    source = CorrelatedNoiseSource(Float64, Fs, 2, 0.3, 0.99)
                    cn = read(source, Fs * 1)
                    bn = bandpass_filter(cn, 300, 700, Fs)
                    mn = amplitude_modulate(bn, 40, Fs)
                    im = ITD_modulate(mn, 8, 24, -24, Fs)

                    source = CorrelatedNoiseSource(Float64, Fs, 2, 0.3, 0.99)
                    cn = read(source, Fs * 1)
                    bn = bandpass_filter(cn, 300, 700, Fs)
                    mn = amplitude_modulate(bn, 40, Fs)
                    im = ITD_modulate(mn, 8, 48, -48, Fs)
                end
            end

        end

        @testset "RMS" begin

            @testset "Abstract Arrays" begin

                for desired_rms = 01:0.1:1
                    bn = bandpass_noise(Fs * 30, 2, 300, 700, Fs)
                    bn = set_RMS(bn, desired_rms)
                    @test rms(bn) ≈ desired_rms
                end
            end

            @testset "Sampled Signals" begin

                for desired_rms = 01:0.1:1
                    source = CorrelatedNoiseSource(Float64, 48000, 2, 1, 0.5)
                    bn = read(source, Fs * 30)
                    bn = set_RMS(bn, desired_rms)
                    @test rms(bn) ≈ desired_rms
                end
            end

        end

        @testset "Ramps" begin

            @testset "Ramp on" begin

                for ramp_length = [1, 2]
                    bn = bandpass_noise(Fs * 5, 2, 300, 700, Fs)
                    bn = ramp_on(bn, Fs * ramp_length)
                    @test rms(bn[1:Fs, :]) < rms(bn[2*Fs:3*Fs, :])
                end
            end

            @testset "Ramp off" begin

                for ramp_length = [1, 2]
                    bn = bandpass_noise(Fs * 5, 2, 300, 700, Fs)
                    bn = ramp_off(bn, Fs * ramp_length)
                    @test rms(bn[end-Fs:end, :]) < rms(bn[2*Fs:3*Fs, :])
                end
            end
        end

        @testset "ITD" begin

            for desired_itd = -100:10:100

                source = CorrelatedNoiseSource(Float64, Fs, 2, 0.3, 0.9)
                cn = read(source, Fs * 5)
                bn = bandpass_filter(cn, 300, 700, Fs)
                bn = set_ITD(bn, desired_itd)
                lags = round.(Int, -150:1:150)
                c = crosscor(bn[:, 2], bn[:, 1], lags)
                x, idx = findmax(c)
                @test lags[idx] == desired_itd
            end
        end

    end

    @testset "Plotting" begin
    # =======================

        @testset "SpectroTempral" begin

        source = CorrelatedNoiseSource(Float64, Fs, 2, 0.3, 0.99)
        cn = read(source, Fs * 1)
        bn = bandpass_filter(cn, 300, 700, Fs)
        mn = amplitude_modulate(bn, 40, Fs)
        im = ITD_modulate(mn, 8, 24, -24, Fs)
        p = PlotSpectroTemporal(im, 48000)
        @test isa(p, Plots.Plot) == true

        source = NoiseSource(Float64, 48000, 2, 0.1)
        sink = DummySampleSink(Float64, 48000, 2)
        some_data = read(source, 3.0u"s"/frames)
        @pipe read(source, 3.0u"s"/frames) |>  write(sink, _)
        p = PlotSpectroTemporal(sink)
        @test isa(p, Plots.Plot) == true
        p = plot(sink)
        @test isa(p, Plots.Plot) == true

        p = plot_cross_correlation(sink)
        @test isa(p, Plots.Plot) == true
        p = plot_cross_correlation(sink, lags=0.1u"s")
        @test isa(p, Plots.Plot) == true
        p = plot_cross_correlation(some_data)
        @test isa(p, Plots.Plot) == true

        end
    end
end

# =======================
# Stream Processing
# =======================

@testset "Online Stimuli" begin

    @testset "Generators" begin

        @testset "NoiseSource generator" begin

            # Test instansiation
            source = NoiseSource(Float64, 48.0u"kHz", 2)
            @test source.samplerate == 48000
            @test source.samplerate == samplerate(source)

            source = NoiseSource(Float64, 48u"kHz", 2)
            @test source.samplerate == 48000

            source = NoiseSource(Float64, 48000u"Hz", 2)
            @test source.samplerate == 48000

            source = NoiseSource(Float64, 48000.1u"Hz", 2)
            @test source.samplerate == 48000.1

            source = NoiseSource(Float64, 44100, 1)
            @test source.samplerate == 44100

            source = NoiseSource(Float64, 48000, 2)
            @test source.samplerate == 48000

            # Test read
            a = read(source, 48000)
            @test size(a) == (48000, 2)
            @test std(a) ≈ 1 atol = 0.01

            a = read(source, 1u"s")
            @test size(a) == (48000, 2)
            @test std(a) ≈ 1 atol = 0.01

            # Test result
            for deviation = 0.1:0.1:1.3
                source = NoiseSource(Float64, 48000, 1, deviation)
                a = read(source, 48000*3)
                @test std(a) ≈ deviation atol = 0.01
            end
        end

        @testset "CorrelatedNoiseSource generator" begin

            # Test instansiation
            source = CorrelatedNoiseSource(Float64, 48.0u"kHz", 2, 1, 0.1)
            @test source.samplerate == 48000
            @test source.samplerate == samplerate(source)

            source = CorrelatedNoiseSource(Float64, 48u"kHz", 2, 1, 1)
            @test source.samplerate == 48000

            source = CorrelatedNoiseSource(Float64, 46u"Hz", 2, 1, 1)
            @test source.samplerate == 46

            source = CorrelatedNoiseSource(Float64, 44100.0u"Hz", 2, 1, 1)
            @test source.samplerate == 44100

            # Test read
            source = CorrelatedNoiseSource(Float64, 48000, 2, 1, 0.1)
            a = read(source, 48000)
            @test size(a) == (48000, 2)
            @test std(a) ≈ 1 atol = 0.01

            # Test behaviour
            for deviation = 0.1:0.2:1.3
                for correlation = 0.0:0.2:1
                    source = CorrelatedNoiseSource(Float64, 48000, 2, deviation, correlation)
                    a = read(source, 3u"s")
                    @test std(a) ≈ deviation atol = 0.025
                    @test cor(a.data)[2, 1] ≈ correlation atol = 0.025
                    @test interaural_coherence(a.data) ≈ correlation atol = 0.025
                end
            end
        end

        @testset "Harmonic Complex" begin

            # Test instansiation
            source = SinusoidSource(Float64, 48000u"Hz", 2)
            @test source.samplerate == 48000
            source = SinusoidSource(Float64, 48000.7u"Hz", [100, 324, 55, 999])
            @test source.samplerate == 48000.7
            source = SinusoidSource(Float64, 48000, 300:300:3000)

            source = SinusoidSource(Float64, 48000, 2000)
            a = read(source, 48000)
            @test size(a) == (48000, 1)

            freqs = collect(200:200:2400.0)
            source = SinusoidSource(Float64, 48000, freqs)
            a = read(source, 48000)
            b = welch_pgram(vec(a.data), fs=a.samplerate)
            maxs_cart = Images.findlocalmaxima(power(b))
            maxs = [idx[1] for idx in maxs_cart]
            maxs = maxs[power(b)[maxs] .> 0.02]
            @test freq(b)[maxs] == freqs

        end

    end

    @testset "Modifiers" begin

        @testset "Input / Output" begin

            desired_rms = 0.3
            for num_channels = 1:2
                for frames = [100, 50, 20]
                    source = NoiseSource(Float64, Fs, num_channels, desired_rms)
                    sink = DummySampleSink(Float64, 48000, num_channels)
                    for idx = 1:frames
                        @pipe read(source, 1.0u"s"/frames) |>  write(sink, _)
                    end
                    @test size(sink.buf, 1) == 48000
                    @test size(sink.buf, 2) == num_channels
                    @test rms(sink.buf) ≈ desired_rms  atol = 0.01
                end
            end

            source = NoiseSource(Float64, 96000u"Hz", 2, 0.1)
            @test source.samplerate == 96000
            source = NoiseSource(Float64, 96u"kHz", 2)
            @test source.samplerate == 96000
            sink = DummySampleSink(Float64, 48000u"Hz", 3)
            @test sink.samplerate == 48000
            sink = DummySampleSink(Float64, 48000.0u"Hz", 3)
            @test sink.samplerate == 48000
            sink = DummySampleSink(Float64, 48u"kHz", 4)
            @test sink.samplerate == 48000
            sink = DummySampleSink(Float64, 48.0u"kHz", 4)
            @test sink.samplerate == 48000
        end

        @testset "Amplification" begin

            # Test different ways of instanciating the modifier
            @test Amplification().target == 0
            @test Amplification().current == 0
            @test Amplification().change_limit == Inf
            @test Amplification().enable == true
            @test Amplification(1).target == 1.0
            @test Amplification(1, 3).target == 1.0
            @test Amplification(1, 3).current == 3.0
            @test Amplification(1, 3, 4, false).change_limit == 4.0
            @test Amplification(1, 3, 4, false).enable == false
            @test Amplification(target=3).target == 3
            @test Amplification(current=3).current == 3
            @test Amplification(change_limit=3).change_limit == 3
            @test Amplification(enable=false).enable == false
            @test Amplification(enable=false, current=2).enable == false
            @test Amplification(enable=false, current=2).current == 2

            # Ensure disabled does nothing
            source = NoiseSource(Float64, 48000, 2, 1)
            sink = DummySampleSink(Float64, 48000, 2)
            amp = Amplification(2, 3, 0.5, false)
            data = read(source, 1u"s")
            @pipe data |> modify(amp, _) |>  write(sink, _)
            @test sink.buf == data.data

            desired_rms = 0.3
            amp_mod = 0.1
            num_channels = 1

            @testset "Static" begin

                source = NoiseSource(Float64, Fs, num_channels, desired_rms)
                sink = DummySampleSink(Float64, 48000, num_channels)
                amp = Amplification(amp_mod, amp_mod, 0.005)

                for idx = 1:200
                    @pipe read(source, 0.01u"s") |> modify(amp, _) |>  write(sink, _)
                end
                @test size(sink.buf, 1) == 96000
                @test size(sink.buf, 2) == num_channels
                @test rms(sink.buf) ≈ desired_rms * amp_mod  atol = 0.01
            end

            @testset "Dynamic" begin

                source = NoiseSource(Float64, Fs, num_channels, desired_rms)
                sink = DummySampleSink(Float64, 48000, num_channels)
                amp = Amplification(target=amp_mod, current=amp_mod, change_limit=0.5)

                for idx = 1:200
                    if idx == 100
                        setproperty!(amp, :target, 1.0)
                    end
                    @pipe read(source, 0.01u"s") |> modify(amp, _) |>  write(sink, _)
                end
                @test size(sink.buf, 1) == 96000
                @test size(sink.buf, 2) == num_channels
                @test rms(sink.buf[1:48000]) ≈ desired_rms * amp_mod  atol = 0.01
                @test rms(sink.buf[48000:96000]) ≈ desired_rms atol = 0.01
            end
        end

        @testset "Filtering" begin

            desired_rms = 0.3
            amp_mod = 0.1
            for num_channels = 1:2

                @testset "Channels: $num_channels" begin

                    source = NoiseSource(Float64, Fs, num_channels, desired_rms)
                    sink = DummySampleSink(Float64, Fs, num_channels)

                    lower_bound = 1000
                    upper_bound = 4000

                    responsetype = Bandpass(lower_bound, upper_bound; fs=Fs)
                    designmethod = Butterworth(14)
                    zpg = digitalfilter(responsetype, designmethod)
                    f1 = DSP.Filters.DF2TFilter(zpg)
                    f2 = DSP.Filters.DF2TFilter(zpg)
                    filters = [f1]
                    if num_channels == 2; filters = [filters[1], f2]; end
                    bandpass = AuditoryStimuli.Filter(filters)

                    for idx = 1:500
                        @pipe read(source, 0.01u"s") |> modify(bandpass, _) |>  write(sink, _)
                    end
                    @test size(sink.buf, 1) == 48000 * 5
                    @test size(sink.buf, 2) == num_channels

                    for chan = 1:num_channels

                        spec = welch_pgram(sink.buf[:, chan], 12000, fs=Fs)

                        val, idx_lb = findmin(abs.(freq(spec) .- lower_bound))
                        val, idx_bl = findmin(abs.(freq(spec) .- (lower_bound - 500)))
                        @test (amp2db(power(spec)[idx_lb]) - amp2db(power(spec)[idx_bl])) > 10

                        val, idx_ub = findmin(abs.(freq(spec) .- upper_bound))
                        val, idx_bu = findmin(abs.(freq(spec) .- (upper_bound + 500)))
                        @test (amp2db(power(spec)[idx_ub]) - amp2db(power(spec)[idx_bu])) > 10
                    end

                    # Test turning filter off
                    setproperty!(bandpass, :enable, false)

                    for idx = 1:1500
                        @pipe read(source, 0.01u"s") |> modify(bandpass, _) |>  write(sink, _)
                    end
                    @test size(sink.buf, 1) == 48000 * 20
                    @test size(sink.buf, 2) == num_channels

                    for chan = 1:num_channels

                        spec = welch_pgram(sink.buf[48000* 5:48000*20, chan], 12000, fs=Fs)

                        val, idx_lb = findmin(abs.(freq(spec) .- lower_bound))
                        val, idx_bl = findmin(abs.(freq(spec) .- (lower_bound - 500)))
                        @test (amp2db(power(spec)[idx_lb]) - amp2db(power(spec)[idx_bl])) < 6

                        val, idx_ub = findmin(abs.(freq(spec) .- upper_bound))
                        val, idx_bu = findmin(abs.(freq(spec) .- (upper_bound + 500)))
                        @test (amp2db(power(spec)[idx_ub]) - amp2db(power(spec)[idx_bu])) < 6
                    end
                end
            end
        end

        @testset "Amplitude Modulation" begin

            for num_channels = 1:5

                @testset "Channels: $num_channels" begin

                    source = NoiseSource(Float64, Fs, num_channels)
                    sink = DummySampleSink(Float64, 48000, num_channels)
                    am = AmplitudeModulation(10u"Hz")

                    for idx = 1:10
                        @pipe read(source, 0.01u"s") |> modify(am, _) |>  write(sink, _)
                    end
                    @test size(sink.buf, 1) == 4800
                    @test size(sink.buf, 2) == num_channels

                    start_mod = Statistics.sum(abs.(sink.buf[1:1000, :]))
                    mid_mod = Statistics.sum(abs.(sink.buf[2000:3000, :]))
                    end_mod = Statistics.sum(abs.(sink.buf[3800:4800, :]))

                    @test mid_mod > start_mod
                    @test mid_mod > end_mod

                    # Test different ways of instanciating the modifier
                    @test AmplitudeModulation(1u"MHz").rate == 1000000u"Hz"
                    @test AmplitudeModulation(1u"kHz").rate == 1000u"Hz"
                    @test AmplitudeModulation(1u"mHz").rate == 0.001u"Hz"
                    @test AmplitudeModulation(10u"Hz").rate == 10u"Hz"
                    @test AmplitudeModulation(1.0u"Hz").rate == 1u"Hz"
                    @test AmplitudeModulation(10u"Hz", 0.0).rate == 10u"Hz"
                    @test AmplitudeModulation(10.3u"Hz", 0.0).rate == 10.3u"Hz"
                    @test AmplitudeModulation(10u"Hz", π).rate == 10u"Hz"
                    @test AmplitudeModulation(10u"Hz", π, 0.0).rate == 10u"Hz"
                    @test AmplitudeModulation(10u"Hz", π, 0.5).rate == 10u"Hz"
                    @test AmplitudeModulation(10.8u"Hz", π, 0.5).rate == 10.8u"Hz"
                    @test AmplitudeModulation(1.0u"kHz", π, 1.5).rate == 1000u"Hz"
                    @test AmplitudeModulation(1.0u"kHz", π, 1.5, false).enable == false
                    @test AmplitudeModulation(10u"Hz", π, 1.5).depth == 1.5
                    @test AmplitudeModulation(rate=3u"Hz").rate == 3.0u"Hz"
                    @test AmplitudeModulation(phase=3).phase == 3
                    @test AmplitudeModulation(depth=0.5).depth == 0.5

                    # Test defaults
                    @test AmplitudeModulation().rate == 0u"Hz"
                    @test AmplitudeModulation().phase == π
                    @test AmplitudeModulation().depth == 1
                    @test AmplitudeModulation().enable == true
                    @test AmplitudeModulation().time == 0

                    # Test bad instansiation
                    @test_logs (:error, "You must use units for modulation rate.")  AmplitudeModulation(33)

                    # Ensure disabled does nothing
                    sink = DummySampleSink(Float64, 48000, num_channels)
                    amp = AmplitudeModulation(enable=false)
                    data = read(source, 1u"s")
                    @pipe data |> modify(amp, _) |>  write(sink, _)
                    @test sink.buf == data.data

                end
            end
        end

        @testset "ITD" begin

            # Test instansiation
            @test TimeDelay().channel == 1
            @test TimeDelay(2).channel == 2
            @test TimeDelay(1, 22).delay == 22
            @test TimeDelay(1, 22, false).delay == 22
            @test TimeDelay(1, 22, false, ones(22, 1)).buffer == ones(22, 1)
            @test TimeDelay(delay=33, buffer=zeros(33, 1)).delay == 33
            @test TimeDelay(channel=33).channel == 33
            @test TimeDelay(enable=false, channel=3).channel == 3
            @test TimeDelay(enable=false, channel=3).enable == false
            @test TimeDelay(2, 0.5u"ms", samplerate=48u"kHz").delay == 24
            @test TimeDelay(2, 1u"ms", samplerate=48u"kHz").delay == 48
            @test_logs (:info, "Not rounded number of samples 23.04") TimeDelay(2, 0.48u"ms", samplerate=48u"kHz")

            # Test correct behaiour
            for desired_itd = -100:10:100

                source = CorrelatedNoiseSource(Float64, 48000, 2, 0.3, 1)
                if desired_itd >= 0
                    itd = TimeDelay(2, desired_itd)
                else
                    itd = TimeDelay(1, -desired_itd)
                end
                sink = DummySampleSink(Float64, 48000, 2)

                for idx = 1:10
                    @pipe read(source, 0.1u"s") |> modify(itd, _) |>  write(sink, _)
                end

                lags = round.(Int, -150:1:150)
                c = crosscor(sink.buf[:, 1], sink.buf[:, 2], lags)
                x, idx = findmax(c)
                @test lags[idx] == desired_itd
            end
        end
    end
end

@testset "Signal Metrics" begin

    @testset "Interaural Coherence" begin

        for correlation = 0.2:0.2:0.8
            source = CorrelatedNoiseSource(Float64, 48000, 2, 0.4, correlation)
            a = read(source, 3u"s")
            @test interaural_coherence(a.data) ≈ correlation atol = 0.025
            @test interaural_coherence(a.data, lags=100) ≈ correlation atol = 0.025
            @test interaural_coherence(a) ≈ correlation atol = 0.025
            @test interaural_coherence(a, lags=1u"s") ≈ correlation atol = 0.025
            @test interaural_coherence(a, lags=1.0u"s") ≈ correlation atol = 0.025
            @test interaural_coherence(a, lags=900u"ms") ≈ correlation atol = 0.025
        end

        # Test still works when ITD applied
        correlation = 0.5
        delay_samples = 100
        delay_time = (delay_samples / 48000)u"s"
        source = CorrelatedNoiseSource(Float64, 48000, 2, 0.4, correlation)
        itd = TimeDelay(1, 22)
        sink = DummySampleSink(Float64, 48000, 2)
        @pipe read(source, 3u"s") |> modify(itd, _) |> write(sink, _)
        @test interaural_coherence(sink) ≈ correlation atol = 0.025
        @test interaural_coherence(sink, lags=2*delay_time) ≈ correlation atol = 0.025

        # Whereas if the range of lags doesnt encompas the maximum the iac will be
        # less than the global maximum
        @test interaural_coherence(sink, lags=0.1*delay_time) < (correlation - 0.1)

    end
end


text:
# Latafat, Themelis, Patrinos, "Block-coordinate and incremental aggregated
# proximal gradient methods for nonsmooth nonconvex problems."
# arXiv:1906.10053 (2019).
#
# Latafat. "Distributed proximal algorithms for large-scale structured optimization"
# PhD thesis, KU Leuven, 7 2020.
#
# Mairal, "Incremental majorization-minimization optimization with application to
# large-scale machine learning."
# SIAM Journal on Optimization 25, 2 (2015), 829–855.
#
# Defazio, Domke, "Finito: A faster, permutable incremental gradient method
# for big data problems."
# In International Conference on Machine Learning (2014), pp. 1125-1133.
#

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms.IterationTools
# using ProximalAlgorithms: LBFGS, update!, mul!
using Printf
using Base.Iterators # to use take and enumerate functions
using Random
using StatsBase: sample
# using BregmanBC
using Flux

import ProximalOperators: prox!, prox
export solution, epoch_count

abstract type CIAO_iterable end #? what is this? where is used?

include("Finito_basic.jl")
include("Finito_LFinito.jl")
include("Finito_adaptive.jl")
include("Finito_LFinito_lbfgs.jl")
include("Finito_DLFinito.jl")
include("Finito_DLFinito_lbfgs.jl")
include("Finito_LFinito_lbfgs_adaptive.jl")
include("Finito_LFinito_lbfgs_adaptive_DNN.jl")

struct Finito{R<:Real}
    γ::Maybe{Union{Array{R},R}} #? this is defined in CIAOAlgorithms, where CIAOAlgorithms is called/used? #other questions in that file ...
    sweeping::Int8
    LFinito::Bool
    lbfgs::Bool
    memory::Int
    η::R
    β::R
    adaptive::Bool
    DeepLFinito::Tuple{Bool,Int, Int}
    minibatch::Tuple{Bool,Int}
    maxit::Int
    verbose::Bool
    freq::Int
    α::R
    tol::R
    tol_b::R
    DNN_training::Bool
    function Finito{R}(;
        γ::Maybe{Union{Array{R},R}} = nothing,
        sweeping = 1,
        LFinito::Bool = false,
        lbfgs::Bool = false,
        memory::Int = 6, # lbfgs memory
        η::R = 0.7,
        β::R = 1/50,
        adaptive::Bool = false,
        # DeepLFinito::Tuple{Bool,Int} = (false, 3),
        DeepLFinito::Tuple{Bool,Int, Int} = (false, 3, 3),
        minibatch::Tuple{Bool,Int} = (false, 1),
        maxit::Int = 10000,
        verbose::Bool = false,
        freq::Int = 10000,
        α::R = R(0.999), # R is type
        tol::R = R(1e-8),
        tol_b::R = R(1e-9),
        DNN_training::Bool = false,
    ) where {R}
        @assert γ === nothing || minimum(γ) > 0
        @assert maxit > 0
        @assert memory >= 0
        @assert tol > 0
        @assert tol_b > 0
        @assert freq > 0
        new(γ, sweeping, LFinito, lbfgs, memory, η, β, adaptive, DeepLFinito, minibatch, maxit, verbose, freq, α, tol, tol_b, DNN_training)
    end
end

function (solver::Finito{R})( # this is a function definition. if solver = Finito(; kwargs...), and then solver(x0;kwargs...) is called, this function will be executed.
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}

    stop(state) = false # the stopping function for halt function

    disp(it, state) = @printf "%5d | %.3e  \n" it state.hat_γ

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    if solver.LFinito
        if solver.DeepLFinito[1]
            iter = FINITO_DLFinito_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                solver.DeepLFinito[2],
                solver.DeepLFinito[3],
            )
        else
            iter = FINITO_LFinito_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
            )
        end
    elseif solver.lbfgs
        if solver.DeepLFinito[1]
            iter = FINITO_DFlbfgs_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
            LBFGS(x0, solver.memory),
            solver.DeepLFinito[2],
            solver.DeepLFinito[3]
        )
        elseif solver.adaptive
            iter = FINITO_lbfgs_adaptive_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.η,
                solver.β,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                LBFGS(x0, solver.memory),
                solver.adaptive
            )
        else
            iter = FINITO_lbfgs_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.β,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                LBFGS(x0, solver.memory),
            )
        end
    elseif solver.adaptive
        iter = FINITO_adaptive_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.tol,
            solver.tol_b,
            solver.sweeping,
            solver.α,
        )
    else
        iter = FINITO_basic_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    end

    iter = halt(iter, stop) #? where is halt defined?
    iter = take(iter, solver.maxit)
    iter = enumerate(iter)

    num_iters, state_final = nothing, nothing
    for (it_, state_) in iter  # unrolling the iterator (acts as tee and loop functions in the tutorial)
        # see https://docs.julialang.org/en/v1/manual/interfaces/index.html
        if solver.verbose && mod(it_, solver.freq) == 0
            disp(it_, state_)
        end
        num_iters, state_final = it_, state_
    end
    if solver.verbose && mod(num_iters, solver.freq) !== 0
        disp(num_iters, state_final)
    end # for the final iteration
    return solution(state_final), num_iters
end

"""
    Finito([γ, sweeping, LFinito, adaptive, minibatch, maxit, verbose, freq, tol, tol_b])

Instantiate the Finito algorithm for solving fully nonconvex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x) + g(x)

where `f_i` are smooth and `g` is possibly nonsmooth, all of which may be nonconvex.

If `solver = Finito(args...)`, then the above problem is solved with

	solver(x0, [F, g, N, L])

where F is an array containing f_i's, x0 is the initial point, and L is an array of
smoothness moduli of f_i's; it is optional in the adaptive mode or if γ is provided.

Optional keyword arguments are:
* `γ`: an array of N stepsizes for each coordinate
* `sweeping::Int` 1 for uniform randomized (default), 2 for cyclic, 3 for shuffled
* `LFinito::Bool` low memory variant of the Finito/MISO algorithm
* `adaptive::Bool` to activate adaptive smoothness parameter computation
* `minibatch::(Bool,Int)` to use batchs of a given size
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10000`), frequency of verbosity.
* `α::R` parameter where γ_i = αN/L_i
* `tol::Real` (default: `1e-8`), absolute tolerance for the adaptive case
* `tol_b::R` tolerance for the backtrack (default: `1e-9`)
"""

Finito(::Type{R}; kwargs...) where {R} = Finito{R}(; kwargs...) #? outer constructor? why is needed? where it is used? Type{R}?
Finito(; kwargs...) = Finito(Float64; kwargs...)

"""
If `solver = Finito(args...)`, then

    itr = iterator(solver, x0, [F, g, N, L])

is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

function iterator( #? how it is called?
    solver::Finito{R},
    x0::Union{AbstractArray{C},Tp};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
    S = nothing,
    F_full = nothing,
    data = nothing,
    DNN_config::Tdnn
) where {R,C<:RealOrComplex{R},Tp,Tdnn}
    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    w0 = DNN_config()
    if solver.DNN_training
        iter = FINITO_lbfgs_adaptive_DNN_iterable(
            F,
            F_full,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.η,
            solver.β,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
            LBFGS(w0, solver.memory),
            solver.adaptive,
            solver.tol_b,
            data,
            DNN_config
        )
    elseif solver.LFinito
        if solver.DeepLFinito[1]
            iter = FINITO_DLFinito_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                solver.DeepLFinito[2],
                solver.DeepLFinito[3],
            )
        else
            iter = FINITO_LFinito_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
            )
        end
    elseif solver.lbfgs
        if solver.DeepLFinito[1]
            iter = FINITO_DFlbfgs_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
            LBFGS(x0, solver.memory),
            solver.DeepLFinito[2],
            solver.DeepLFinito[3]
        )
        elseif solver.adaptive
            iter = FINITO_lbfgs_adaptive_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.η,
                solver.β,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                LBFGS(x0, solver.memory),
                solver.adaptive,
                solver.tol_b,
                S,
            )
        else
            iter = FINITO_lbfgs_iterable(
                F,
                g,
                x0,
                N,
                L,
                solver.γ,
                solver.β,
                solver.sweeping,
                solver.minibatch[2],
                solver.α,
                LBFGS(x0, solver.memory),
            )
        end
    elseif solver.adaptive
        iter = FINITO_adaptive_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.tol,
            solver.tol_b,
            solver.sweeping,
            solver.α,
        )
    else
        iter = FINITO_basic_iterable(
            F,
            g,
            x0,
            N,
            L,
            solver.γ,
            solver.sweeping,
            solver.minibatch[2],
            solver.α,
        )
    end
    return iter
end

#### TODO
    # remove finito_adaptive


text:
## types

struct FoodWeb
    links::Matrix
    nicheValues::Vector{Float64}
    distribution::Vector{Int64}
    index::Vector{UnitRange}
end

## functions

function generateFoodWeb(S::Int, C::Float64)
    @assert (C <= 0.5 && C >= 0.0)
    β = 1.0 / (2.0 * C) - 1.0
    n = vcat(0.0,sort(rand(S - 1)))
    # n = sort(rand(S))
    r = n .* rand(Beta(1,β),S)
    c = (n .- (r ./ 2)) .* rand(S) .+ (r ./ 2)
    A = getAdjacency(n, r, c)
    while (isolatedSpecies(A))
        replaceIsolatedSpecies!(A, n, r, c, β)
        A = getAdjacency(n, r, c)
    end
    T = findall(sum(A, dims=2)[:] .== 0)
    B = findall(sum(A, dims=1)[:] .== 0)
    I = findall((sum(A, dims=2)[:] .> 0) .& (sum(A, dims=1)[:] .> 0))
    newInd = vcat(B,I,T)
    distri = [size(B,1), size(I,1), size(T,1)]
    ind = [1:distri[1],distri[1]+1:size(B,1) + size(I,1),size(B,1) + size(I,1) + 1:S]
    A = getAdjacency(n[newInd], r[newInd], c[newInd])
    return FoodWeb(A,n[newInd],distri,ind)
end

function isolatedSpecies(A::Matrix)
    return size(findall((sum(A, dims=2)[:] .== 0) .& (sum(A, dims=1)[:] .== 0)),1) > 0
end

function replaceIsolatedSpecies!(A, n, r, c, β)
    ind = findall((sum(A, dims=2)[:] .== 0) .& (sum(A, dims=1)[:] .== 0))
    # ind = findall(x -> x != 1, ind)
    for i in ind
        n[i] = rand()
        r[i] = n[i] * rand(Beta(1,β))
        c[i] = (n[i] - (r[i] / 2)) * rand().+ (r[i] / 2)
    end
end

function getAdjacency(n::Vector,r::Vector,c::Vector)
    A = zeros(size(n,1),size(n,1))
    for i = 1:size(n,1)
        for j = 2:size(n,1)
            A[i,j] = (n[i] <= c[j] + r[j] / 2) && (n[i] >= c[j] - r[j] / 2)
        end
    end
    return A - Diagonal(A)
end


text:
# ---------------------------------------------------
module CoreInterfaces
# ---------------------------------------------------

abstract type ICore end
abstract type IConfig end
abstract type ILog end

end


text:
using EngEconomics, Roots

# Given
aInit = -1000000
aAnnualSavings = 100000
aAnnualG = 50000
MARR = 0.10
aLife = 6

# Determine Payback Period


text:
export PassManager,
       add!, dispose

abstract type PassManager end
reftype(::Type{T}) where {T<:PassManager} = API.LLVMPassManagerRef

add!(pm::PassManager, pass::Pass) =
    API.LLVMAddPass(ref(pm), ref(pass))

dispose(pm::PassManager) = API.LLVMDisposePassManager(ref(pm))

#
# Module pass manager
#

export ModulePassManager, run!

@checked struct ModulePassManager <: PassManager
    ref::reftype(PassManager)
end

ModulePassManager() = ModulePassManager(API.LLVMCreatePassManager())

function ModulePassManager(f::Core.Function, args...)
    mpm = ModulePassManager(args...)
    try
        f(mpm)
    finally
        dispose(mpm)
    end
end

run!(mpm::ModulePassManager, mod::Module) =
    convert(Core.Bool, API.LLVMRunPassManager(ref(mpm), ref(mod)))

#
# Function pass manager
#

export FunctionPassManager,
       initialize!, finalize!, run!

@checked struct FunctionPassManager <: PassManager
    ref::reftype(PassManager)
end

FunctionPassManager(mod::Module) =
    FunctionPassManager(API.LLVMCreateFunctionPassManagerForModule(ref(mod)))

function FunctionPassManager(f::Core.Function, args...)
    fpm = FunctionPassManager(args...)
    try
        f(fpm)
    finally
        dispose(fpm)
    end
end

initialize!(fpm::FunctionPassManager) =
    convert(Core.Bool, API.LLVMInitializeFunctionPassManager(ref(fpm)))
finalize!(fpm::FunctionPassManager) =
    convert(Core.Bool, API.LLVMFinalizeFunctionPassManager(ref(fpm)))

run!(fpm::FunctionPassManager, f::Function) =
    convert(Core.Bool, API.LLVMRunFunctionPassManager(ref(fpm), ref(f)))


text:
function assign(isfix::BitArray{1})
    fix = 1
    n = length(isfix)
    fixcluster = vcat(fix, zeros(n))

    for i in 1:n

        if !isfix[i]
            fix += 1
        end

        fixcluster[i+1] = fix
    end

    fixcluster
end


text:
immutable ErrorBarGeometry <: Gadfly.GeometryElement
end

immutable XErrorBarGeometry <: Gadfly.GeometryElement
end

immutable YErrorBarGeometry <: Gadfly.GeometryElement
end

const errorbar = ErrorBarGeometry
const xerrorbar = XErrorBarGeometry
const yerrorbar = YErrorBarGeometry

function element_aesthetics(::ErrorBarGeometry)
    [:x, :y, :xmin, :xmax, :ymin, :ymax]
end

function element_aesthetics(::YErrorBarGeometry)
    [:x, :ymin, :ymax]
end

function element_aesthetics(::XErrorBarGeometry)
    [:y, :xmin, :xmax]
end

# Generate a form for the errorbar geometry.
#
# Args:
#   geom: errorbar geometry
#   theme: the plot's theme
#   aes: aesthetics
#
# Returns:
#   A compose Form.
#
function render(geom::ErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    # check for X and Y error bar aesthetics
    if isempty(Gadfly.undefined_aesthetics(aes, element_aesthetics(xerrorbar())...))
        xform = render(xerrorbar(), theme, aes)
    else
        xform = empty_form
    end
    if isempty(Gadfly.undefined_aesthetics(aes, element_aesthetics(yerrorbar())...))
        yform = render(yerrorbar(), theme, aes)
    else
        yform = empty_form
    end
    compose(canvas(units_inherited=true, order=3), xform, yform)
end

function render(geom::YErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.errorbar", aes,
                                     element_aesthetics(geom)...)
    Gadfly.assert_aesthetics_equal_length("Geom.errorbar", aes,
                                          element_aesthetics(geom)...)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    # What are the actual extents. We can'n use 1/6
    caplen = theme.errorbar_cap_length/2

    compose(
        canvas(units_inherited=true, order=3),
        combine(
            # top cap
            lines([[(x*cx - caplen, ymax), (x*cx + caplen, ymax)]
                   for (x, ymax) in zip(aes.x, aes.ymax)]...),

            # error bar
            lines([[(x*cx, ymax), (x*cx, ymin)]
                   for (x, ymin, ymax) in zip(aes.x, aes.ymin, aes.ymax)]...),

            # bottom cap
            lines([[(x*cx - caplen, ymin), (x*cx + caplen, ymin)]
                   for (x, ymin) in zip(aes.x, aes.ymin)]...)),

        stroke([theme.highlight_color(c) for c in aes.color]),
        linewidth(theme.line_width),
        aes.color_key_continuous == true ?
            svgclass("geometry") :
            svgclass([@sprintf("geometry color_%s",
                escape_id(aes.color_label([c])[1]))
                          for c in aes.color]))
end

function render(geom::XErrorBarGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.errorbar", aes,
                                     element_aesthetics(geom)...)
    Gadfly.assert_aesthetics_equal_length("Geom.errorbar", aes,
                                          element_aesthetics(geom)...)

    colored = aes.color != nothing
    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(ColorValue[theme.default_color])
    aes = inherit(aes, default_aes)

    # What are the actual extents. We can'n use 1/6
    caplen = theme.errorbar_cap_length/2

    compose(
        canvas(units_inherited=true, order=3),
        combine(
            # left cap
            lines([[(xmin, y*cy - caplen), (xmin, y*cy + caplen)]
                   for (xmin, y) in zip(aes.xmin, aes.y)]...),

            # error bar
            lines([[(xmin, y*cy), (xmax, y*cy)]
                   for (xmin, xmax, y) in zip(aes.xmin, aes.xmax, aes.y)]...),

            # right cap
            lines([[(xmax, y*cy - caplen), (xmax, y*cy + caplen)]
                   for (xmax, y) in zip(aes.xmax, aes.y)]...)),

        stroke([theme.highlight_color(c) for c in aes.color]),
        linewidth(theme.line_width),
        (aes.color_key_continuous == true || !colored) ?
            svgclass("geometry") :
            svgclass([@sprintf("geometry color_%s", escape_id(aes.color_label(c)[1]))
                      for c in aes.color]))
end


text:
"""

# ScriptEntry

Define processing requirements for chapter scripts

### Constructor
```julia
scriptentry(scriptfile;; nb=true, exe=true, doc=true)
```

### Required arguments
```julia
* `scriptfile::AbstractString`        : Script file
```

### Optional arguments
```julia
* `nb::Bool`      : Generate a notebook version in notebooks directory
* `exe::Bool`     : Execute the notebook (for testing or documentation purposes)
* `doc::Bool`     : Insert documention into Github pages
```

If exe = false and doc = true it is assumed the appropriate .md files
have been manually created and stored in docs/src/_nn_/... (Travis
will terminate if runs take too long).

"""
struct ScriptEntry
  scriptfile::AbstractString
  nb::Bool
  exe::Bool
  doc::Bool
end

"""

# scriptentry

Constructor for ScriptEntry objects.

### Constructor
```julia
scriptentry(scriptfile;; nb=true, exe=true, doc=true)
```

### Required arguments
```julia
* `scriptfile::AbstractString`        : Script file
```

### Optional arguments
```julia
* `nb::Bool`      : Generate a notebook version in notebooks directory
* `exe::Bool`     : Execute the notebook (for testing or documentation purposes)
* `doc::Bool`     : Insert documention into Github pages
```

If exe = false and doc = true it is assumed the appropriate .md files
have been manually created and stored in docs/src/_nn_/... (Travis
will terminate if runs take too long).

"""
scriptentry(scriptfile; nb=true, exe=true, doc=true) = ScriptEntry(scriptfile, nb, exe, doc)

script_dict = DataStructures.OrderedDict{AbstractString, Vector{ScriptEntry}}(
  "00" => [
    scriptentry("clip-01-03.jl"),
    scriptentry("clip-04-05.jl")
  ],
  "02" => [
    scriptentry("clip-01-02.jl"),
    scriptentry("clip-03-05.jl"),
    scriptentry("clip-06-07.jl"),
    scriptentry("clip-08.jl"),
    scriptentry("m2.1s.jl")
  ],
  "03" => [
    scriptentry("clip-01.jl"),
    scriptentry("clip-02-05.jl"),
    scriptentry("clip-05s.jl"),
    scriptentry("clip-06-10.jl"),
    scriptentry("clip-11-16.jl"),
    scriptentry("clip-17s.jl")
  ],
  "04" => [
    scriptentry("m4.1s.jl"),
    scriptentry("clip-01-06.jl"),
    scriptentry("clip-07-13s.jl"),
    scriptentry("clip-14-20.jl"),
    scriptentry("clip-21-23.jl"),
    scriptentry("clip-24-29s.jl"),
    scriptentry("clip-30s.jl"),
    scriptentry("clip-38s.jl"),
    scriptentry("clip-43s.jl"),
    scriptentry("clip-45-47s.jl"),
    scriptentry("clip-48-54s.jl")
  ],
  "05" => [
    scriptentry("clip-01-05s.jl"),
    scriptentry("clip-06-10s.jl")
   ],
  "09" => [
    scriptentry("clip-01.jl"),
    scriptentry("clip-02.jl"),
    scriptentry("clip-03d.jl")
  ]
);


text:
function e005(N)
    return prod_of_factors(N)
end

function prod_of_factors(n)
    factors = zeros(Int8, n+1)
    for i = 2:n
        for (f, mul) in Channel(c -> get_prime_factors(c, i))
            if factors[f] < mul
                factors[f] = mul
            end
        end
    end
    prod = 1
    for i = 2:n
        prod *= i^factors[i]
    end
    return prod
end

function solution005()
    N = 20
    answer = e005(N)
    println("The answer to problem 005 is $answer")
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
@testset "Simulator: Hopper (2D)" begin

    s = get_simulation("hopper_2D", "flat_2D_lc", "flat")
    model = s.model
    env = s.env

    ref_traj = get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/hopper_2D/gaits/gait_in_place.jld2"),
        load_type = :joint_traj)
    T = ref_traj.H
    h = ref_traj.h

    for t = 1:T
    	r = ContactImplicitMPC.residual(model, env, ref_traj.z[t], ref_traj.θ[t], 0.0)
    	@test norm(r, Inf) < 1.0e-5
    end

    # initial conditions
    q0 = SVector{model.dim.q}(ref_traj.q[1])
    q1 = SVector{model.dim.q}(ref_traj.q[2])

    # simulator
    sim = ContactImplicitMPC.simulator(s, q0, q1, h, T,
        p = ContactImplicitMPC.open_loop_policy([SVector{model.dim.u}(ut) for ut in ref_traj.u]),
        ip_opts = ContactImplicitMPC.InteriorPointOptions(
    		r_tol = 1.0e-8, κ_tol = 1.0e-8, solver = :lu_solver),
        sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

    # simulate
    status = ContactImplicitMPC.simulate!(sim, verbose = false)
    @test status
    @test norm(ref_traj.q[end] - sim.traj.q[end], Inf) < 1.0e-3
end


text:
using LinearAlgebra, Plots, JLD2
using Flux, Flux.Zygote, Optim, DiffEqFlux
using KitBase, ProgressMeter
import KitML
include("math.jl")

cd(@__DIR__)
include("math.jl")
model = KitML.load_model("best_model.h5"; mode = :tf)

begin
    # space
    x0 = -3.5
    x1 = 3.5
    y0 = -3.5
    y1 = 3.5
    nx = 100
    ny = 100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 5.0
    cfl = 0.7
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # moments
    L = 1
    ne = GetBasisSize(L,3)
    phi = zeros(Float32, ne, nx, ny)
    α = zeros(Float32, ne, nx, ny)
    m = ComputeSphericalBasisKarth(points,1,3)
end

function is_absorb(x::T, y::T) where {T<:Real}
    cds = Array{Bool}(undef, 11) # conditions

    cds[1] = -2.5<x<-1.5 && 1.5<y<2.5
    cds[2] = -2.5<x<-1.5 && -0.5<y<0.5
    cds[3] = -2.5<x<-1.5 && -2.5<y<-1.5
    cds[4] = -1.5<x<-0.5 && 0.5<y<1.5
    cds[5] = -1.5<x<-0.5 && -1.5<y<-0.5
    cds[6] = -0.5<x<0.5 && -2.5<y<-1.5
    cds[7] = 0.5<x<1.5 && 0.5<y<1.5
    cds[8] = 0.5<x<1.5 && -1.5<y<-0.5
    cds[9] = 1.5<x<2.5 && 1.5<y<2.5
    cds[10] = 1.5<x<2.5 && -0.5<y<0.5
    cds[11] = 1.5<x<2.5 && -2.5<y<-1.5

    if any(cds) == true
        return true
    else
        return false
    end
end

# particle
begin
    σs = zeros(Float64, nx, ny)
    σa = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if is_absorb(pspace.x[i, j], pspace.y[i, j])
            σs[i, j] = 0.0
            σa[i, j] = 10.0
        else
            σs[i, j] = 1.0
            σa[i, j] = 0.0
        end
    end
    σt = σs + σa
    σq = zeros(Float64, nx, ny)
    for i = 1:nx, j = 1:ny
        if -0.5<pspace.x[i, j]<0.5 && -0.5<pspace.y[i, j]<0.5
            σq[i, j] = 1.0 / (4.0 * π)
        else
            σq[i, j] = 0.0
        end
    end
end

#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], σq')

for j = 1:nx
    for i = 1:ny
        phi[1, i, j] = 1e-6
    end
end

global t = 0.0
flux1 = zeros(ne, nx + 1, ny)
flux2 = zeros(ne, nx, ny + 1)
αT = zeros(Float32, nx*ny, ne)
phiT = zeros(Float32, nx*ny, ne)

anim = @animate for iter = 1:100
    println("Iteration $(iter)")

    # neural closure of the system
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :, :][:]
    end
    αT .= KitML.neural_closure(model, phiT)
    @inbounds for i = 1:ne
        α[i, :, :] .= reshape(αT[:, i], nx, ny)
    end
    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            phi[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)
        end
    end

    # flux
    fη1 = zeros(nq)
    @inbounds for j = 1:ny
        for i = 2:nx
            KitBase.flux_kfvs!(fη1, KitBase.maxwell_boltzmann_dual.(α[:, i-1, j]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 1], dt)

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    @inbounds for i = 1:nx
        for j = 2:ny
            KitBase.flux_kfvs!(fη2, KitBase.maxwell_boltzmann_dual.(α[:, i, j-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 2], dt)

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (weights .* fη2))
            end
        end
    end

    # update
    @inbounds for j = 2:ny-1
        for i = 2:nx-1
            for q = 1:1
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (σs[i, j] * phi[q, i, j] - σt[i, j] * phi[q, i, j]) * dt +
                    σq[i, j] * dt * 100
            end

            for q = 2:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy +
                    (-σt[i, j] * phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
    contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])
end

cd(@__DIR__)
gif(anim, "lattice_mn_pure_neural.gif")
#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
# Benchmark adopted from the house_heating.ipynb example
#

using Simulate, Random, Distributions, Statistics, BenchmarkTools
import Dates.now

const Th = 40     # temperature of heating fluid
const R = 1e-6    # thermal resistance of room insulation
const α = 2e6     # represents thermal conductivity and capacity of the air
const β = 3e-7    # represents mass of the air and heat capacity
const Δt = 1//60  # update every minute

mutable struct House
    Te::Float64   # environment temperature
    Tr::Float64   # room temperature
    heating::Bool # heating on or off
    η::Float64    # efficiency factor reducing R if doors or windows are open
    rng::MersenneTwister       # random number generator
    datTr::Vector{Float64}     # room temperature data
end

Δte(t, t1, t2) = cos((t-10)*π/12) * (t2-t1)/2  # change of a sinusoidal Te

function Δtr(h::House)                         # calculate room temperature change
    Δqc = (h.Tr - h.Te)/(R * h.η)                # cooling rate
    Δqh = h.heating ? α * (Th - h.Tr) : 0    # heating rate
    return β * (Δqh - Δqc)
end

function setTemperatures!(h::House, t1=8, t2=20)               # change the temperatures
    h.Te += Δte(tau(), t1, t2) * 2π/1440 + rand(h.rng, Normal(0, 0.1))
    h.Tr += Δtr(h) * Δt
    push!(h.datTr, h.Tr)       # append room temperature data
end

function switch!(h::House, t1=20, t2=23)                       # simulate the thermostat
    if h.Tr ≥ t2
        h.heating = false
        event!(Fun(switch!, h, t1, t2), Fun((h, x)-> h.Tr ≤ x, h, t1))  # setup a conditional event
        # event!(Fun(switch, t1, t2), @val :Tr :≤ t1)  # setup a conditional event
    elseif h.Tr ≤ t1
        h.heating = true
        event!(Fun(switch!, h, t1, t2), Fun((h, x)-> h.Tr ≥ x, h, t2))  # setup a conditional event
        # event!(Fun(switch, t1, t2), @val :Tr :≥ t2)  # setup a conditional event
    end
end

function people!(clk::Clock, h::House)
    delay!(clk, 6 + rand(Normal(0, 0.5)))         # sleep until around 6am
    sleeptime = 22 + rand(Normal(0, 0.5))    # calculate bed time
    while tau(clk)%24 < sleeptime
        h.η = rand()                         # open door or window
        delay!(clk, 0.1 * rand(Normal(1, 0.3)))   # for some time
        h.η = 1.0                            # close it again
        delay!(clk, rand())
    end
end

function setup()
    reset!(𝐶)                                      # reset the clock
    global house = House(11, 20, false, 1.0, MersenneTwister(122), Float64[])
    Random.seed!(1234)

    for i in 1:2                              # put 2 people in the house
        process!(Prc(i, people!, house))       # run processes
    end
    periodic!(Fun(setTemperatures!, house), Δt)  # set sampling function
    switch!(house)                                     # start the thermostat
end

# first run to compile
setup()
run!(𝐶, 48)

# take measurements
println(now())
setup()
@time info = run!(𝐶, 48)
println(info)
println("measurements:", length(house.datTr), " mean_tr=", mean(house.datTr))

t = run(@benchmarkable run!(𝐶, 48) setup=setup() evals=1 seconds=15.0 samples=100)


text:
function write_pdb(data::Data, name_file::AbstractString)
    io = open(name_file, "w")
    write_pdb_CRYST1(io, data)
    write_pdb_ATOM(io, data)
    write_pdb_END(io)
    close(io)
end

function write_pdb_CRYST1(io, data::Data)
    box_size = data.data_basic.box_size[:, 2] - data.data_basic.box_size[:, 1]
    box_tilt = data.data_basic.box_tilt
    x = sqrt(box_size[1]^2 + box_tilt[1]^2 + box_tilt[2]^2)
    y = sqrt(box_size[2]^2 + box_tilt[3]^2)
    z = box_size[3]

    @printf(io, "CRYST1")
    @printf(io, "%9.3f", x)
    @printf(io, "%9.3f", y)
    @printf(io, "%9.3f", z)
    @printf(io, "%7.2f", 90)
    @printf(io, "%7.2f", 90)
    @printf(io, "%7.2f", 90)
    @printf(io, "\n")
end

function write_pdb_ATOM(io, data::Data_Unit)
    atom_vec = data.vec_atom
    atom_name = data.data_str.atom_name
    res_name = data.data_str.res_name
    if isa(data.data_str, Structure_Wat)
        type_name = "WT1"
        chain_id = "W"
    elseif isa(data.data_str, Structure_Ion)
        type_name = "ION"
        chain_id = "I"
    else
        type_name = "U0"
        chain_id = "U"
    end
    for (id_atom, atom) in enumerate(atom_vec)
        @printf(io, "%-6s", "ATOM")                            #  1 -  6        Record name   "ATOM  "
        @printf(io, "%5d", atom.atom)                   #  7 - 11        Integer       serial       Atom  serial number.
        @printf(io, "%2s%-3s", " ", uppercase(atom_name[atom.typ]))
        #@printf(io, "%5s", atom_name[atom.typ])        # 13 - 16        Atom          name         Atom name.
        #@printf(io, " ")                                # 17             Character     altLoc       Alternate location indicator.
        @printf(io, "%1s%-4s", " ", res_name)                    # 18 - 20        Residue name  resName      Residue name.
        @printf(io, "%1s", chain_id)                   # 22             Character     chainID      Chain identifier.
        @printf(io, "%4d", atom.mol)                    # 23 - 26        Integer       resSeq       Residue sequence number.
        @printf(io, "%4s", " ")                        # 27             AChar         iCode        Code for insertion of residues.
        @printf(io, "%8.3f", atom.coord[1])             # 31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
        @printf(io, "%8.3f", atom.coord[2])             # 39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
        @printf(io, "%8.3f", atom.coord[3])             # 47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
        @printf(io, "%6.2f", 1.0)                       # 55 - 60        Real(6.2)     occupancy    Occupancy.
        @printf(io, "%6.2f", 0.0)                       # 61 - 66        Real(6.2)     tempFactor   Temperature  factor.
        @printf(io, "%6s", " ")
        @printf(io, "%-3s%3s", type_name, uppercase(atom_name[atom.typ][1:end-1]))# 73 - 78        LString(2)    element      Element symbol, right-justified.
        # 79 - 80        LString(2)    charge       Charge  on the atom.
        @printf(io, "\n")
    end
end

function write_pdb_ATOM(io, data::Data_Sum)
    atom_vec = data.vec_atom
    atom_name = []
    res_name = []
    type_name = []
    chain_name = []
    chain_id = []
    for (id, str) in enumerate(data.vec_str)
        atom_name = vcat(atom_name, str.atom_name)
        res_name = vcat(res_name, str.res_name) # Residue name
        append!(chain_id, ones(Int64, str.num_atom_types)*id)
        if isa(str, Structure_Wat)
            type_name = vcat(type_name, "WT1")
            chain_name = vcat(chain_name, "W")
        elseif isa(str, Structure_Ion)
            type_name = vcat(type_name, "ION")
            chain_name = vcat(chain_name, "I")
        else
            type_name = vcat(type_name, "U0")
            chain_name = vcat(chain_name, "U")
        end
    end
    for (id_atom, atom) in enumerate(atom_vec)
        @printf(io, "%-6s", "ATOM")                                         #  1 -  6        Record name   "ATOM  "
        @printf(io, "%5d", atom.atom)                                       #  7 - 11        Integer       serial       Atom  serial number.
        @printf(io, "%2s%-3s", " ", uppercase(atom_name[atom.typ]))         # 13 - 16        Atom          name         Atom name.
        @printf(io, "%1s%-4s", " ",res_name[chain_id[atom.typ]])            # 17 - 20        Residue name  resName      Residue name.
        @printf(io, "%1s", chain_name[chain_id[atom.typ]])                  # 22             Character     chainID      Chain identifier.
        @printf(io, "%4d", atom.mol)                                        # 23 - 26        Integer       resSeq       Residue sequence number.
        @printf(io, "%4s", " ")                                             # 27             AChar         iCode        Code for insertion of residues.
        @printf(io, "%8.3f", atom.coord[1])                                 # 31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
        @printf(io, "%8.3f", atom.coord[2])                                 # 39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
        @printf(io, "%8.3f", atom.coord[3])                                 # 47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
        @printf(io, "%6.2f", 0.5)                                           # 55 - 60        Real(6.2)     occupancy    Occupancy.
        @printf(io, "%6.2f", 300)                                           # 61 - 66        Real(6.2)     tempFactor   Temperature  factor.
        @printf(io, "%6s", " ")
        @printf(io, "%-3s%3s", type_name[chain_id[atom.typ]], uppercase(atom_name[atom.typ][1:end-1]))# 73 - 78        LString(2)    element      Element symbol, right-justified.
        #@printf(io, "%9s%6s", type_name[chain_id[atom.typ]], atom_name[atom.typ])# 73 - 78        LString(2)    element      Element symbol, right-justified.
        # 79 - 80        LString(2)    charge       Charge  on the atom.
        @printf(io, "\n")
    end
end

function write_pdb_END(io)
    @printf(io, "END")
end


text:
function edit(style::Style, dict::Dict, customization::Dict{Symbol,Any})
  screen         = customization[:screen]
  yposition      = float32(screen.area.value.h)

  glypharray     = Array(Romeo.GLGlyph{Uint16}, 0)
  visualizations = RenderObject[]
  xgap = 60f0
  ygap = 20f0
  lineheight = 24f0
  currentline = 0
  aabb = AABB(Vec3(0), Vec3(0))
  i = 0
  for (name,value) in dict
    if method_exists(edit, (typeof(value),))
        currentline       += int(abs(aabb.max[2]-aabb.min[2])/lineheight) + i
        yposition         -= lineheight*2
        i = 3
        append!(glypharray, Romeo.GLGlyph{Uint16}[Romeo.GLGlyph(c, currentline, k, 0) for (k,c) in enumerate(string(name)*":")])
        visual, signal     = Romeo.edit(value, style, screen=screen)
        aabb               = visual.boundingbox(visual)

        translatm          = translationmatrix(Vec3(xgap - aabb.min[1], yposition - aabb.max[2], 0))
        visual[:model]     = translatm * visual[:model]
        dict[name]         = signal
        yposition          -= abs(aabb.max[2]-aabb.min[2]) + lineheight
        push!(visualizations, visual)
    end
  end
  labels = visualize(glypharray, screen=screen, color=rgba(0.19, 0.70,0.88,1.0), model=translationmatrix(Vec3(30, float32(screen.area.value.h), 0)))
  push!(visualizations, labels)
  visualizations
end
function edit(style::Style, obj::RenderObject, customization::Dict{Symbol,Any})
  edit(style, obj.uniforms, customization)
end


text:
function penalize!{T<:Integer,U<:String}(sd::Array{Array{Resource{T},1},1},
                                         tp::TProblem{T,U})
    avail = BitArray{1}[]
    for dim in 2:-1:1
        push!(avail, bitpack(map(isavailable, sd[dim])))
    end
    for dim in 1:2, r in sd[dim]
        if r.quant == 0
            r.l = r.m = r.p = r.q = 0
            continue
        end
        r.l == 0 || !avail[dim][r.l] || !avail[dim][r.m] || continue
        rsort = filter(x->avail[dim][x], vec(slicedim(tp.tsort[dim],dim,r.i)))
        rcost = vec(slicedim(tp.toc, dim, r.i))[rsort]
        if length(rsort) == 1
            r.l = r.m = rsort[1]
            r.p = r.q = rcost[1]
        else
            r.l, r.m = rsort[1:2]
            r.p = rcost[2] - rcost[1]
            r.q = rcost[1]
        end
    end
    nothing
end

function vogel{T<:Integer,U<:String}(tp::TProblem{T,U})
    sdcnt = collect(size(tp.toc))
    sol = spzeros(T, sdcnt[1], sdcnt[2])
    sd = Array{Resource{T},1}[]
    for dim in 1:2
        push!(sd, [Resource(dim, i, tp.sd[dim][i]) for i in 1:sdcnt[dim]])
    end
    while any(map(isavailable, sd[1])) && any(map(isavailable, sd[2]))
        penalize!(sd, tp)
        a = maximum([sd[1], sd[2]])
        b = sd[rem1(a.dim+1,2)][a.l]
        if a.dim == 2 # swap to make a supply and b demand
            a, b = b, a
        end
        expend = min(a.quant, b.quant)
        sol[a.i, b.i] = expend
        a.quant -= expend
        b.quant -= expend
    end
    return sol
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#TODO convert to real MVC, unify interface, plugin style, fix grid interface
#TODO additional abstraction layer encapsulating an element plus data and handler

"""
    ctrl(D, x)
Dictionary `D`, Tuple `x=(name::String, ::DataType, init, v...)`; used by mkControlGrid
"""
function ctrl(D::Dict{AbstractString, Any}, x)
	local tmp = Void
	if x[2] == Int
		tmp = @SpinButton(x[4][1], x[4][end], step(x[4]))
		setproperty!(tmp, :value, x[3])
		signal_connect(tmp, "value_changed") do w
			D[x[1]] = getproperty(w, :value, x[2])
		end
		signal_emit(tmp, "value_changed", Void)
	elseif x[2] == Float64
		tmp = @Scale(:h, x[4][1], x[4][end], step(x[4]))
		local adj = @Adjustment(tmp)
		setproperty!(adj, :value, x[3])
		signal_connect(adj, "value_changed") do w
			D[x[1]] = getproperty(w, :value, x[2])
		end
		signal_emit(tmp, "value_changed", Void)
	elseif x[2] == Bool
		tmp = @CheckButton()
		setproperty!(tmp, :active, x[3])
		signal_connect(tmp, "toggled") do w
			D[x[1]] = getproperty(w, :active, x[2])
		end
		signal_emit(tmp, "toggled", Void)
	end

	if tmp == Void
		warn("No handling for type $(x[2]) defined.")
		return Void
	end

	D[x[1]] = x[3]

	setproperty!(tmp, :hexpand, true)
	setproperty!(tmp, :vexpand, false)
	return tmp
end

function button(handler::Function, label::AbstractString)
	B = @Button(label)
	signal_connect(B, "clicked") do w
		@schedule handler(w)
	end
	return B
end

"""
    mkControlGrid(D, C)
Creates a grid of controls with labels, handlers and encapsulated storage (Dictionary `D`)
`c` in `C` is Tuple `(name::String, ::DataType, init, v...)`
"""
function mkControlGrid(D::Dict{AbstractString, Any}, C)
	Control = map(C) do X
		reduce([], X) do A,a
			isa(a, Gtk.GtkWidget) && return [A; a]
			a == Void && return [A; Void]
			return [A; @Label(a[1]); ctrl(D, a)]
		end
	end

	G = @Grid()
	setproperty!(G, :column_spacing, 10)
	setproperty!(G, :row_spacing, 6)
	setproperty!(G, :margin, 6)
	setproperty!(G, :column_homogeneous, true)
	# setproperty!(G, :row_homogeneous, true)

	S = reduce((A,a)->lcm(A,length(a)), 1, Control)
	for i in 1:length(Control)
		l = S÷length(Control[i])
		for j in 1:length(Control[i])
			Control[i][j] == Void && continue
			G[(j-1)*l + (1:l), i] = Control[i][j]
		end
	end

	return G
end


text:
using Circuitscape

const FILE_CELLMAP_ASC = joinpath("data_files", "cellmap.asc")
const FILE_CELLMAP_ASC_GZ = joinpath("data_files", "cellmap.asc.gz")
const FILE_INCLUDE_MATRIX = joinpath("data_files", "include_matrix.txt")

@assert false == try Circuitscape._check_file_exists("cellmap.asc") end
@assert Circuitscape._check_file_exists(FILE_CELLMAP_ASC)

@assert Circuitscape.FILE_TYPE_AAGRID == Circuitscape._guess_file_type(FILE_CELLMAP_ASC)
@assert Circuitscape.FILE_TYPE_AAGRID == Circuitscape._guess_file_type(FILE_CELLMAP_ASC_GZ)
@assert Circuitscape.FILE_TYPE_INCL_PAIRS_AAGRID == Circuitscape._guess_file_type(FILE_INCLUDE_MATRIX)
@assert Circuitscape.FILE_TYPE_INCL_PAIRS == Circuitscape._guess_file_type(FILE_INCLUDE_MATRIX)

f = FileAAGrid(FILE_CELLMAP_ASC)
@assert f.hdr.ncols == 10
@assert f.hdr.nrows == 10
@assert f.hdr.nodata == -9999
@assert f.data_type == "float"
@assert f.data[1,1] == 1.0

fname = tempname()
write(fname, f)
f = FileAAGrid(fname)
rm(fname)
@assert f.hdr.ncols == 10
@assert f.hdr.nrows == 10
@assert f.hdr.nodata == -9999
@assert f.data_type == "float"
@assert f.data[1,1] == 1.0

fname = "$(tempname()).gz"
write(fname, f, compress=true)
f = FileAAGrid(fname)
rm(fname)
@assert f.hdr.ncols == 10
@assert f.hdr.nrows == 10
@assert f.hdr.nodata == -9999
@assert f.data_type == "float"
@assert f.data[1,1] == 1.0