text:
#=
Popis
  Funkcia vypocita diameter grafu.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String
  components : parameter pre prepinanie graf/max_komponent :: Bool

Navratova hodnota
  : diameter maximalneho komponentu :: Float64
  g_diameter : diameter grafu :: Float64

Poznamky
  Funkcia podla hodnoty parametra components vypocita hodnotu
  pre cely graf alebo pre najvacsi komponent.
=#

function graph_diameter(g::GenericGraph; dist_key::AbstractString="weight", components::Bool=false)

  if components
    comp_diameter = Dict{Array{ExVertex,1}, Float64}()
    comp_length = Int64[]

    for comp in connected_components(g)
      comp_diameter[comp] = maximum([vertex_eccentricity(v.index, g) for v in comp])
      push!(comp_length, length(comp))
    end

    max_comp = maximum(comp_length)

    return [d[2] for d=filter!(x -> length(x[1]) == max_comp, collect(zip(keys(comp_diameter), values(comp_diameter))))]
  else
    g_diameter = [vertex_eccentricity(v.index, g) for v in vertices(g)]

    return maximum(g_diameter)
  end

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
@testset "PeriodicSelector" begin
    st = Date(2019, 1, 1)
    ed = Date(2019, 2, 1)
    date_range = st:Day(1):ed

    @testset "1 week period, 1 day stride" begin

        selector = PeriodicSelector(Week(1))
        @test selector isa DateSelector

        result = partition(date_range, selector)
        @test sort(vcat(result...)) == date_range
        @test isempty(intersect(result...))

        @test all(isequal(Week(1)), diff(result.holdout))
    end

    @testset "2 week period, 5 day stride" begin
        selector = PeriodicSelector(Week(2), Day(5))
        result = partition(date_range, selector)
        @test sort(vcat(result...)) == date_range
        @test isempty(intersect(result...))

        expected_holdout = [
            st:Day(1):st + Day(4)...,
            st + Week(2):Day(1):st + Week(2) + Day(4)...,
            st + Week(4):Day(1):st + Week(4) + Day(3)...,
        ]
        @test diff(result.holdout) == Day.([1, 1, 1, 1, 10, 1, 1, 1, 1])
    end

    @testset "1 week period, 1 day stride, 2 day offset" begin
        selector = PeriodicSelector(Week(1), Day(1), Day(2))

        result = partition(date_range, selector)
        @test sort(vcat(result...)) == date_range
        @test isempty(intersect(result...))

        expected_holdout = [st+Day(1):Week(1):ed...]
        @test result.holdout == expected_holdout
    end

    @testset "BiWeekly" begin
        st1 = Date(1900, 1, 1)
        ed1 = Date(1900, 2, 1)
        date_range1 = st1:Day(1):ed1

        selector = PeriodicSelector(Week(2), Week(1))

        result = partition(date_range1, selector)
        @test isempty(intersect(result...))

        expected_holdout = [
            st1:Day(1):st1+Day(6)...,
            st1+Week(2):Day(1):st1+Week(3)-Day(1)...,
            st1+Week(4):Day(1):ed1...,
        ]

        expected_validation = [
            st1+Week(1):Day(1):st1+Week(2)-Day(1)...,
            st1+Week(3):Day(1):st1+Week(4)-Day(1)...,
        ]

        @test result.holdout == expected_holdout
        @test result.validation == expected_validation
    end

    @testset "Day of week" begin
        @testset "Weekends as holdout" begin
            selector = PeriodicSelector(Week(1), Day(2), Day(5))
            validation, holdout = partition(date_range, selector)

            @test isequal(
                sort(unique(dayname.(validation))),
                ["Friday", "Monday", "Thursday", "Tuesday", "Wednesday"]
            )
            @test sort(unique(dayname.(holdout))) == ["Saturday", "Sunday"]
        end

        @testset "Week days as holdout" begin
            selector = PeriodicSelector(Week(1), Day(5))
            validation, holdout = partition(date_range, selector)

            @test sort(unique(dayname.(validation))) == ["Saturday", "Sunday"]
            @test isequal(
                sort(unique(dayname.(holdout))),
                ["Friday", "Monday", "Thursday", "Tuesday", "Wednesday"]
            )
        end
    end

    @testset "stride, offset, period wrong domain" begin
        @test_throws DomainError PeriodicSelector(Day(0))
        @test_throws DomainError PeriodicSelector(Day(7), Day(0))

        @test_throws DomainError PeriodicSelector(Hour(24))
        @test_throws MethodError PeriodicSelector(Day(7), Hour(48))
        @test_throws MethodError PeriodicSelector(Day(7), Day(3), Hour(72))
    end

    @testset "errors if stride > period" begin
        @test_throws ArgumentError PeriodicSelector(Day(2), Day(3))
        @test_throws ArgumentError PeriodicSelector(Week(1), Day(8))
    end

    @testset "Invarient to date_range (period $period, stride $stride)" for
        (period, stride) in ((Week(2), Week(1)), (Day(30), Day(10)), (Day(5), Day(2)))

        @testset "offset $offset" for offset in Day.((0, 1, 2, 3))
            selector = PeriodicSelector(period, stride, offset)
            initial_range = Date(2029, 1, 1):Day(1):Date(2030, 1, 31)
            initial_sets = map(Set, partition(initial_range, selector))

            later_range = initial_range + Day(20)  # 20 days later.
            later_sets = map(Set, partition(later_range, selector))

            @test initial_sets.holdout ∩ later_sets.validation == Set()
            @test initial_sets.validation ∩ later_sets.holdout == Set()
        end
    end
end

ext: jl


text:
using Test
using GtkExtensions, Gtk

@testset "TextIter" begin

    import Gtk: GtkTextIter, mutable

    w = GtkWindow()
    b = GtkTextBuffer()
    b.text[String] = "test"
    v = GtkTextView(b)

    push!(w,v)
    showall(w)

    it = GtkTextIter(b,1)
end

max_stars_repo_name: jonathanBieler/GtkExtensions

max_issues_repo_name: jonathanBieler/GtkExtensions

max_forks_repo_name: jonathanBieler/GtkExtensions


text:
using UKCarbonIntensityData
using Documenter

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

makedocs(;
    modules=[UKCarbonIntensityData],
    authors="Ian Goddard",
    repo="https://github.com/ianlmgoddard/UKCarbonIntensityData.jl/blob/{commit}{path}#{line}",
    sitename="UKCarbonIntensityData.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ianlmgoddard.github.io/UKCarbonIntensityData.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md"
    ],
    checkdocs=:exports,
    strict=false,
)

deploydocs(;
    repo="github.com/Ianlmgoddard/UKCarbonIntensityData.jl",
    devbranch="main",
    #push_preview=true,
)


text:
@testset "transformedkernel" begin
    rng = MersenneTwister(123456)
    x = rand(rng)*2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    s = rand(rng)
    s2 = rand(rng)
    v = rand(rng, 3)
    k = SqExponentialKernel()
    kt = TransformedKernel(k,ScaleTransform(s))
    ktard = TransformedKernel(k,ARDTransform(v))
    @test kt(v1, v2) == transform(k, ScaleTransform(s))(v1, v2)
    @test kt(v1, v2) == transform(k, s)(v1,v2)
    @test kt(v1, v2) ≈ k(s * v1, s * v2) atol=1e-5
    @test ktard(v1, v2) ≈ transform(k, ARDTransform(v))(v1, v2) atol=1e-5
    @test ktard(v1, v2) == transform(k,v)(v1, v2)
    @test ktard(v1, v2) == k(v .* v1, v .* v2)
    @test transform(kt, s2)(v1, v2) ≈ kt(s2 * v1, s2 * v2)
    @test KernelFunctions.kernel(kt) == k
    @test repr(kt) == repr(k) * "\n\t- " * repr(ScaleTransform(s))

    @testset "kernelmatrix" begin
        rng = MersenneTwister(123456)

        Nx = 5
        Ny = 4
        D = 3

        k = SqExponentialKernel()
        t = ScaleTransform(randn(rng))
        kt = TransformedKernel(k, t)

        @testset "$(typeof(x))" for (x, y) in [
            (ColVecs(randn(rng, D, Nx)), ColVecs(randn(rng, D, Ny))),
            (RowVecs(randn(rng, Nx, D)), RowVecs(randn(rng, Ny, D))),
        ]
            @test kernelmatrix(kt, x, y) ≈ kernelmatrix(k, map(t, x), map(t, y))

            @test kernelmatrix(kt, x) ≈ kernelmatrix(k, map(t, x))

            @test kerneldiagmatrix(kt, x) ≈ kerneldiagmatrix(k, map(t, x))

            tmp = Matrix{Float64}(undef, length(x), length(y))
            @test kernelmatrix!(tmp, kt, x, y) ≈ kernelmatrix(kt, x, y)

            tmp_square = Matrix{Float64}(undef, length(x), length(x))
            @test kernelmatrix!(tmp_square, kt, x) ≈ kernelmatrix(kt, x)

            tmp_diag = Vector{Float64}(undef, length(x))
            @test kerneldiagmatrix!(tmp_diag, kt, x) ≈ kerneldiagmatrix(kt, x)
        end
    end
    test_ADs(x->transform(SqExponentialKernel(), x[1]), rand(1))# ADs = [:ForwardDiff, :ReverseDiff])
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# License); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# AS IS BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# runtime detection of compile time features in the native library

module MXRuntime

using ..mx

export LibFeature
export feature_list, isenabled

# defined in include/mxnet/c_api.h
struct LibFeature
  _name::Ptr{Cchar}
  enabled::Bool
end

function Base.getproperty(x::LibFeature, p::Symbol)
  (p == :name) && return unsafe_string(getfield(x, :_name))
  getfield(x, p)
end

Base.show(io::IO, x::LibFeature) =
  print(io, ifelse(x.enabled, "✔", "✖"), " ", x.name)

"""
    feature_list()

Check the library for compile-time features.
The list of features are maintained in libinfo.h and libinfo.cc
"""
function feature_list()
  ref = Ref{Ptr{LibFeature}}(C_NULL)
  s = Ref{Csize_t}(C_NULL)
  @mx.mxcall(:MXLibInfoFeatures, (Ref{Ptr{LibFeature}}, Ref{Csize_t}), ref, s)
  unsafe_wrap(Array, ref[], s[])
end

"""
    isenabled(x::Symbol)::Bool

Returns the given runtime feature is enabled or not.

```julia-repl
julia> mx.isenabled(:CUDA)
false

julia> mx.isenabled(:CPU_SSE)
true
```

See also `mx.feature_list()`.
"""
isenabled(x::Symbol) =
  any(feature_list()) do i
    Symbol(i.name) == x && i.enabled
  end

end  # module MXRuntime

using .MXRuntime


text:
using TensorBoardLogger, Logging
using Flux: gradient, update!, WeightDecay, gpu
using Juno
using Profile
using Statistics: norm
using LinearAlgebra: rmul!
import Flux.Optimise: apply!
using Statistics: norm, mean, std
using Flux: ADAMW, ADAM, RADAM, Momentum, Optimiser
lg=TBLogger("./TesorBoardLog/run", min_level=Logging.Info)

## ---

# python impl uses clip_grad_norm_ to clip the norms.  It calculates the total norm of all gradients then does:
# clip_coef = max_norm / (total_norm + 1e-6) to clip all the gradients
#  See https://pytorch.org/docs/stable/_modules/torch/nn/utils/clip_grad.html#clip_grad_norm_
function clipGlobalNorm!(gs, ps, max=1.0)
	sumnorm = sum(norm(gs[p])  for p in ps)
	if sumnorm > max
		coef = max/sumnorm
		for p in ps
				rmul!(gs[p], coef)
		end
	end
end

## ---
mutable struct SelectiveWeightDecay
	wd::WeightDecay
	noDecayList::IdDict
end

function apply!(o::SelectiveWeightDecay, x, Δ)
	length(o.noDecayList)
	if haskey(o.noDecayList, x)
		return Δ
	end
	apply!(o.wd, x,  Δ)
end

## ---

function cosRate( lr=6e-4, warmuptoken=512*20, len_train_dataset=1115394, blocksize=128, repeat=1)
    finalTok=repeat*len_train_dataset*blocksize

	function cosRateFun(token)
		if token < warmuptoken
			lr_mult = token / max(1, warmuptoken)
		else
	    	progress = (token - warmuptoken) / max(1, finalTok - warmuptoken)
	    	lr_mult = max(0.1, 0.5 * (1.0 + cos(pi * progress)))
		end
    	lr_mult*lr
	end
	cosRateFun
end

updateLR(opt::ADAM, lr::Float64) = opt.eta=lr
updateLR(opt::Optimiser, lr::Float64) = for o in opt.os  updateLR(o, lr) end
updateLR(o, ::Float64)=nothing

defaultValidate()=0
function train!(loss, ps, data, opt, validate=defaultValidate, logger=lg;rateDecay=cosRate(), repeat=1)
	tokenCnt = 0
	for batchCnt in 1:repeat
	  # training_loss is declared local so it will be available for logging outside the gradient calculation.
	  local training_loss
	  ctr=0
	  normList=[]
	  tokensInBatch(b)=size(b[1])[1] * size(b[1])[2]   #number of char in sentences, and number of sentences
	  numberOfBatches=size(data)
	  @progress for next in data
		tokenCnt = tokenCnt + tokensInBatch(next)
		lr=rateDecay(tokenCnt)
		updateLR(opt, lr)
		d = next |> gpu
		@show ctr, numberOfBatches
	    gs = gradient(ps) do
	      training_loss = loss(d...)
	      # Code inserted here will be differentiated, unless you need that gradient information
	      # it is better to do the work outside this block.
	      return training_loss
	    end

		normBeforeAdj=sum(norm(gs[p])  for p in ps)
		push!(normList, normBeforeAdj)
		clipGlobalNorm!(gs, ps)
		#normAfterAdj=sum(norm(gs[p])  for p in ps)

	    # Insert whatever code you want here that needs training_loss, e.g. logging.
	    # logging_callback(training_loss)
	    # Insert what ever code you want here that needs gradient.
	    # E.g. logging with TensorBoardLogger.jl as histogram so you can see if it is becoming huge.

		if ctr % 10 == 0

			##TODO:::: SHOW AVE gradient norm over the last 10 runs
			@show mean(normList), std(normList)
			normList=[]
			@show ctr, training_loss
			@show training_accuracy = validate()
			with_logger(logger) do
			   @info "train" loss=training_loss
			   @info "train" accuracy=training_accuracy
			end
		end
		ctr+=1
	    update!(opt, ps, gs)
	    # Here you might like to check validation set accuracy, and break out to do early stopping.
	  end
  end
end

max_stars_repo_name: daryoush/minGPT

max_issues_repo_name: daryoush/minGPT

max_forks_repo_name: daryoush/minGPT


text:
"Check whether a set of coefficients is compatible with the function set."
compatible_coefficients(dict::Dictionary, coefficients) =
    length(dict) == length(coefficients)

export Expansion
"""
An `Expansion` describes a function using its expansion coefficients in a certain
dictionary.

An expansion acts both like an array, the array of coefficients, and like a
function.

Parameters:
- D is the dictionary.
- C is the type of the expansion coefficients
"""
struct Expansion{S,T,D<:Dictionary{S,T},C}
    dictionary      ::  D
    coefficients    ::  C

    function Expansion{S,T,D,C}(dict::D, coefficients::C) where {S,T,D<:Dictionary{S,T},C}
        @assert compatible_coefficients(dict, coefficients)
        new(dict, coefficients)
    end
end

Expansion(dict::Dictionary) = Expansion(dict, zeros(dict))
Expansion(dict::Dictionary, coef) =
    Expansion{domaintype(dict),codomaintype(dict),typeof(dict),typeof(coef)}(dict, coef)

export Expansion1d, Expansion2d, Expansion3d, Expansion4d
# Warning: not all 2d function sets have SVector{2,T} type, they could have (S,T) type
const Expansion1d{S <: Number,T,D,C} = Expansion{S,T,D,C}
const Expansion2d{S <: Number,T,D,C} = Expansion{SVector{2,S},T,D,C}
const Expansion3d{S <: Number,T,D,C} = Expansion{SVector{3,S},T,D,C}
const Expansion4d{S <: Number,T,D,C} = Expansion{SVector{4,S},T,D,C}

@forward Expansion.dictionary domaintype, prectype, dimension,
                                size, Span, support, interpolation_grid, numtype,
                                ncomponents

@forward Expansion.coefficients length, eachindex, getindex, setindex!
codomaintype(e::Expansion) = promote_type(codomaintype(dictionary(e)), eltype(coefficients(e)))
isreal(e::Expansion) = isreal(e.dictionary) && isreal(eltype(coefficients(e)))

export expansion
"""
    expansion(dict::Dictionary, coefficients)

An expansion of a dictionary with given coefficients.
"""
expansion(dict::Dictionary, coefficients) =
    Expansion(dict, native_coefficients(dict, coefficients))

similar(e::Expansion, coefficients) = Expansion(dictionary(e), coefficients)

dictionary(e::Expansion) = e.dictionary
coefficients(e::Expansion) = e.coefficients

export random_expansion
random_expansion(dict::Dictionary) = Expansion(dict, rand(dict))

# For expansions of composite types, return a Expansion of a subdict
component(e::Expansion, i) = Expansion(component(e.dictionary, i), component(e.coefficients, i))
components(e::Expansion) = map(Expansion, components(e.dictionary), components(e.coefficients))

# This indirect call enables dispatch on the type of the dict of the expansion
(e::Expansion)(x; options...) =
    call_expansion(e, dictionary(e), coefficients(e), x; options...)
(e::Expansion)(x, y; options...) =
    call_expansion(e, dictionary(e), coefficients(e), SVector(x, y); options...)
(e::Expansion)(x, y, z; options...) =
    call_expansion(e, dictionary(e), coefficients(e), SVector(x, y, z); options...)
(e::Expansion)(x, y, z, t; options...) =
    call_expansion(e, dictionary(e), coefficients(e), SVector(x, y, z, t); options...)
(e::Expansion)(x, y, z, t, u...; options...) =
    call_expansion(e, dictionary(e), coefficients(e), SVector(x, y, z, t, u...); options...)

call_expansion(e::Expansion, dict::Dictionary, coefficients, x; options...) =
    eval_expansion(dict, coefficients, x; options...)

function differentiate(e::Expansion, order = difforder(dictionary(e)); options...)
    op = differentiation(codomaintype(e), dictionary(e), order; options...)
    Expansion(dest(op), apply(op,e.coefficients))
end

function antidifferentiate(e::Expansion, order = 1; options...)
    op = antidifferentiation(codomaintype(e), dictionary(e); options...)
    Expansion(dest(op), apply(op,e.coefficients))
end

Base.broadcast(e::Expansion, grid::AbstractGrid) = eval_expansion(dictionary(e), coefficients(e), grid)

# Shorthands for partial derivatives
∂x(f::Expansion) = differentiate(f, 1, 1)
∂y(f::Expansion) = differentiate(f, 2, 1)
∂z(f::Expansion) = differentiate(f, 3, 1)
# Shorthands for partial integrals
∫∂x(f::Expansion) = antidifferentiate(f, 1, 1)
∫∂y(f::Expansion) = antidifferentiate(f, 2, 1)
∫∂z(f::Expansion) = antidifferentiate(f, 3, 1)
# little helper function
ei(dim,i, coefficients) = tuple((coefficients*Matrix{Int}(I, dim, dim)[:,i])...)
# we allow the differentiation of one specific variable through the var argument
differentiate(f::Expansion, var, order) = differentiate(f, ei(dimension(f), var, order))
antidifferentiate(f::Expansion, var, order) = antidifferentiate(f, ei(dimension(f), var, order))

# To be implemented: Laplacian (needs multiplying functions)
## Δ(f::Expansion)

# This is just too cute not to do: f' is the derivative of f. Then f'' is the second derivative, and so on.
adjoint(f::Expansion) = differentiate(f)

∫(f::Expansion) = antidifferentiate(f)

roots(f::Expansion) = roots(dictionary(f), coefficients(f))

# Delegate generic operators
for op in (:extension, :restriction, :evaluation, :approximation, :transform)
    @eval $op(src::Expansion, dest::Expansion) = $op(promote_type(codomaintype(src),codomaintype(dest)), dictionary(src), dictionary(dest))
end

for op in (:interpolation, :approximation)
    @eval $op(s::Expansion) = $op(dictionary(s))
end

differentiation(e::Expansion; options...) = differentiation(codomaintype(e), dictionary(e); options...)
differentiation(e::Expansion, order; options...) = differentiation(codomaintype(e), dictionary(e), order; options...)

show(io::IO, ::MIME"text/plain", fun::Expansion) = show_setexpansion(io, fun, dictionary(fun))
function show_setexpansion(io::IO, fun::Expansion, fs::Dictionary)
    println(io, "A ", dimension(fun), "-dimensional Expansion with ", length(coefficients(fun)), " degrees of freedom.")
    println(io, "Basis: ", repr(fs))
end

# Invoke split_interval on the set and compute the coefficients such that
# the piecewise function agrees with the original one.
split_interval(s::Expansion, x) = Expansion(split_interval_expansion(dictionary(s), coefficients(s), x)...)

##############################
# Arithmetics with expansions
##############################

# Arithmetics are only possible when the basis type is equal.
iscompatible(d1::D, d2::D) where {D<:Dictionary} = true
iscompatible(d1::Dictionary, d2::Dictionary) = false

function promote_length(e1::Expansion, e2::Expansion)
    if length(e1) == length(e2)
        e1, e2
    elseif length(e1) < length(e2)
        extension(e1, e2) * e1, e2
    else
        e1, extension(e2, e1) * e2
    end
end

function (+)(e1::Expansion, e2::Expansion)
    @assert iscompatible(dictionary(e1),dictionary(e2))
    f1, f2 = promote_length(e1, e2)
    Expansion(dictionary(f1), coefficients(f1)+coefficients(f2))
end

function (-)(e1::Expansion, e2::Expansion)
    @assert iscompatible(dictionary(e1),dictionary(e2))
    f1, f2 = promote_length(e1, e2)
    Expansion(dictionary(f1), coefficients(f1)-coefficients(f2))
end

(-)(e::Expansion) = similar(e, -coefficients(e))

function (*)(s1::Expansion, s2::Expansion)
    @assert iscompatible(dictionary(s1),dictionary(s2))
    (mset,mcoefficients) = (*)(dictionary(s1),dictionary(s2),coefficients(s1),coefficients(s2))
    Expansion(mset,mcoefficients)
end

(*)(op::DictionaryOperator, e::Expansion) = apply(op, e)

for op in (:+, :-, :*)
    @eval Base.$op(a::Number, e::Expansion) = Expansion(dictionary(e), $op(a, coefficients(e)))
    @eval Base.$op(e::Expansion, a::Number) = Expansion(dictionary(e), $op(coefficients(e), a))
end

apply(op::DictionaryOperator, e::Expansion) = Expansion(dest(op), op * coefficients(e))

iterate(e::Expansion) = iterate(coefficients(e))

iterate(e::Expansion, state) = iterate(coefficients(e), state)

Base.collect(e::Expansion) = coefficients(e)

Base.BroadcastStyle(e::Expansion) = Base.Broadcast.DefaultArrayStyle{dimension(e)}()

+(f1::Function, f2::Expansion) = (x->f1(x)+f2(x))
+(f1::Expansion, f2::Function) = (x->f1(x)+f2(x))


text:
module SparseMatrixCSRTests

    using SparseArrays
    using LinearAlgebra
    using Gridap.Algebra
    using Test

    maxnz=5
    maxrows=5
    maxcols=5
    maxrowsorcols=7
    int_types=(Int32,Int)
    float_types=(Float32,Float64)
    Bi_types=(0,1)

    for Ti in int_types
        for Tv in float_types
            for Bi in Bi_types
                I = Vector{Ti}()
                J = Vector{Ti}()
                V = Vector{Tv}()
                for (ik, jk, vk) in zip(rand(1:maxrows, maxnz), rand(1:maxcols, maxnz), rand(1:Tv(maxnz), maxnz-1))
                    push_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,ik,jk,vk)
                end
                push_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows,maxcols,maxnz)
                finalize_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows,maxcols)
                CSC = sparse(I, J, V, maxrows,maxcols)
                CSR = sparsecsr(SparseMatrixCSR{Bi,Tv,Ti},I, J, V,maxrows,maxcols)

                @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},1,1)
                @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},1,2)
                @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},2,1)

                _l = 10
                _I, _J, _V = allocate_coo_vectors(SymSparseMatrixCSR{Bi,Tv,Ti},_l)
                @test length(_I) == _l
                @test length(_J) == _l
                @test length(_V) == _l
                @test eltype(_I) == Ti
                @test eltype(_J) == Ti
                @test eltype(_V) == Tv

                V2 = copy(V)
                V2 .= 0
                CSR2 = sparsecsr(SparseMatrixCSR{Bi,Tv,Ti},I, J, V2,maxrows,maxcols)
                copy_entries!(CSR2,CSR)
                @test CSR2 == CSR
                copy_entries!(CSR2,CSR2)
                @test CSR2 == CSR

                show(devnull, CSR);

                @test CSC == CSR
                @test nnz(CSC) == count(!iszero, CSC) == nnz(CSR) == count(!iszero, CSR)

                @test hasrowmajororder(CSR) == true
                @test hascolmajororder(CSR) == false
                @test getptr(CSR)           == CSR.rowptr
                @test getindices(CSR)       == colvals(CSR)

                TCSC = sparse(J, I, V, maxrows, maxcols)
                TCSR = sparsecsr(SparseMatrixCSR{Bi}, J, I, V, maxrows, maxcols)

                @test size(CSC)==size(CSR)==reverse(size(TCSC))==reverse(size(TCSC))

                @test [nzrange(CSC,col) for col in 1:size(CSC,2)] == [nzrange(TCSR,row) for row in 1:size(TCSR,1)]
                @test [nzrange(CSR,row) for row in 1:size(CSR,1)] == [nzrange(TCSC,col) for col in 1:size(TCSC,2)]

                @test nonzeros(CSC) == nonzeros(TCSR) && nonzeros(CSR) == nonzeros(TCSC)

                ICSC,JCSC,VCSC= findnz(CSC)
                ICSR,JCSR,VCSR= findnz(CSR)

                @test sort(ICSC)==sort(ICSR) && sort(JCSC)==sort(JCSR) && sort(VCSC)==sort(VCSR)

                v = rand(size(CSC)[2])
                @test CSC*v == CSR*v

                for cBi in Bi_types
                    if Bi == cBi
                        @test convert(SparseMatrixCSR{Bi}, CSR) === CSR
                    else
                        CSRC = convert(SparseMatrixCSR{cBi}, CSR)
                        @test CSRC == CSR
                        @test CSRC !== CSR
                    end
                end

                for cTi in int_types
                    for cTv in float_types
                        if (Ti,Tv) == (cTi,cTv)
                            @test convert(SparseMatrixCSR{Bi,Tv,Ti}, CSR) === CSR
                        else
                            CSRC = convert(SparseMatrixCSR{Bi,cTv,cTi}, CSR)
                            @test CSRC == CSR
                            @test CSRC !== CSR
                        end
                    end
                end

                vold = getindex(CSR,maxrows,maxcols)
                add_entry!(CSR,1,maxrows,maxcols,+)
                @test getindex(CSR,maxrows,maxcols) == vold+1

                fill_entries!(CSR,0)
                @test all(x->x==0, nonzeros(CSR))

            end

            function test_rectangular(Ti::Type,Tv::Type,nrows,ncols)
               I = Vector{Ti}()
               J = Vector{Ti}()
               V = Vector{Tv}()
               for (ik, jk, vk) in zip(rand(1:nrows, maxnz), rand(1:ncols, maxnz), rand(1:Tv(maxnz), maxnz-1))
                 push_coo!(SparseMatrixCSR{1,Tv,Ti},I,J,V,ik,jk,vk)
               end
               push_coo!(SparseMatrixCSR{1,Tv,Ti},I,J,V,nrows,ncols,maxnz)
               finalize_coo!(SparseMatrixCSR{1,Tv,Ti},I,J,V,nrows,ncols)
               CSC = sparse(I, J, V, nrows,ncols)
               CSR = sparsecsr(I, J, V, nrows,ncols)
               @test CSR == CSC
           end
           test_rectangular(Ti,Tv,maxrowsorcols,maxcols)
           test_rectangular(Ti,Tv,maxrowsorcols,maxcols)
        end
    end
end


text:
using MicroscopyTools, Documenter

 # set seed fixed for documentation
DocMeta.setdocmeta!(MicroscopyTools, :DocTestSetup, :(using MicroscopyTools); recursive=true)
makedocs(modules = [MicroscopyTools],
         sitename = "MicroscopyTools.jl",
         pages = Any[
            "MicroscopyTools.jl" => "index.md",
            "Binning" => "binning.md",
         ]
        )

deploydocs(repo = "github.com/JuliaMicroscopy/MicroscopyTools.jl.git", devbranch="main")

max_stars_repo_name: JuliaMicroscopy/MicroscopyTools

max_issues_repo_name: JuliaMicroscopy/MicroscopyTools


text:
import StatsBase: sample

"""
    sample(N::T, n::Int64) where {T<:AbstractUnipartiteNetwork}

Samples a sub-network from a unipartite network. `n` is the number of species to
have in the sampled network. This functions makes *no* attempt to ensure that
the network is not degenerate, or even has a single interaction. This is the
recommended way to sample a unipartite network.
"""
function sample(N::T, n::Int64) where {T<:AbstractUnipartiteNetwork}
    @assert n <= richness(N)
    ns = StatsBase.sample(species(N), n, replace=false)
    return N[ns]
end

"""
    sample(N::T, n::Tuple{Int64}) where {T<:AbstractUnipartiteNetwork}

Same as `sample`, but work when called with `(n,)` instead of a species number.
This is an accepted way to sample a unipartite network.
"""
function sample(N::T, n::Tuple{Int64}) where {T<:AbstractUnipartiteNetwork}
    return sample(N, n[1])
end

"""
    sample(N::T, n::Tuple{Int64,Int64}) where {T<:AbstractUnipartiteNetwork}

Same as `sample` but called with `(n,n)` instead of a species number. Note that
this will fail if the size requested is not square. This is not a really good
way to sample a unipartite network.
"""
function sample(N::T, n::Tuple{Int64,Int64}) where {T<:AbstractUnipartiteNetwork}
    @assert first(n) == last(n)
    return sample(N, n[1])
end

"""
    sample(N::T, n::Tuple{Int64}) where {T<:AbstractBipartiteNetwork}

Same as `sample` but with a single species number given as `(n,)`, to return a
bipartite network of equal richness on both sides. This is not a very good way
to sample a bipartite network.
"""
function sample(N::T, n::Tuple{Int64}) where {T<:AbstractBipartiteNetwork}
    @assert n[1] <= richness(N; dims=1)
    @assert n[1] <= richness(N; dims=2)
    return sample(N, (n[1],n[1]))
end

"""
    sample(N::T, n::Int64) where {T<:AbstractBipartiteNetwork}

Same thing as `sample` but with a single species number, to return a bipartite
network of equal richness on both sides. This is not a very good way to sample a
bipartite network.
"""
function sample(N::T, n::Int64) where {T<:AbstractBipartiteNetwork}
    @assert n <= richness(N; dims=1)
    @assert n <= richness(N; dims=2)
    return sample(N, (n,n))
end

"""
    sample(N::T, n::Tuple{Int64,Int64}) where {T<:AbstractBipartiteNetwork}

Samples a sub-network from a bipartite network. `n` is the size of the network
to return, *i.e.* number of top and bottom species. This functions makes *no*
attempt to ensure that the network is not degenerate, or even has a single
interaction.

This is the recommended way to sample a bipartite network.
"""
function sample(N::T, n::Tuple{Int64,Int64}) where {T<:AbstractBipartiteNetwork}
    @assert n[1] <= richness(N; dims=1)
    @assert n[2] <= richness(N; dims=2)
    ns1 = StatsBase.sample(species(N; dims=1), n[1], replace=false)
    ns2 = StatsBase.sample(species(N; dims=2), n[2], replace=false)
    return N[ns1,ns2]
end


text:
module CellMLToolkit

using MathML

using SymbolicUtils: FnType, Sym, operation
using ModelingToolkit
using EzXML

include("cellml.jl")

"""
    reads a CellML path or io and returns an ODESystem
"""
function read_cellml(path, tspan)
    xml = readxml(path)
    ml = CellModel(xml, process_cellml_xml(xml))
    ODEProblem(ml, tspan)
end

"""
    parses a CellML XML string and returns an ODESystem
"""
function parse_cellml(xmlstr::AbstractString, tspan)
    xml = parsexml(xmlstr)
    ml = CellModel(xml, process_cellml_xml(xml))
    ODEProblem(ml, tspan)
end

##############################################################################

export CellModel, ODEProblem
export read_cellml, parse_cellml
export list_params, list_initial_conditions, list_states, update_list!
export eval_string, readxml

struct CellModel
    xml::EzXML.Document
    sys::ODESystem
end

getxml(ml::CellModel) = ml.xml
getsys(ml::CellModel) = ml.sys

"""
    constructs a CellModel struct for the CellML model defined in path
"""
function CellModel(path::AbstractString; dependency=true)
    xml = readxml(path)
    CellModel(xml, process_cellml_xml(xml))

    # fall-back option:
    # s = eval_string(xml)
    # eval(Meta.parse(s))
    # sys = ODESystem(eqs)
    # CellModel(xml, sys)
end

list_states(ml::CellModel) = getsys(ml).states
list_params(ml::CellModel) = list_params(getxml(ml))

"""
    extracts the initial conditions from ml.xml and sorts them
    according to ml.sys.states
"""
function list_initial_conditions(ml::CellModel)
    states = list_states(ml)
    u0 = list_initial_conditions(getxml(ml))
    map(x -> value_of(x, u0), states)
end

function value_of(x, u0)
    name = string(x)
    i = findfirst(x -> string(first(x)) == name, u0)
    return last(u0[i])
end

import ModelingToolkit.ODEProblem

"""
    ODEProblem constructs an ODEProblem from a CellModel
"""
function ODEProblem(ml::CellModel, tspan;
        jac=false, level=1, p=list_params(ml), u0=list_initial_conditions(ml))
    ODEProblem(ml.sys, u0, tspan, p; jac=jac)
end

"""
    update_list! updates the value of an item in an initial value list (either p or u0)
    l is the list
    sym is a Symbol pointing to the variable (we also pass a string)
    val is the new value
"""
function update_list!(l, sym::Symbol, val)
    for (i,x) in enumerate(l)
        if first(x).name == sym
            l[i] = (first(x) => val)
            return
        end
    end
    parse_error(e, "param not found: $name")
end

update_list!(l, name::AbstractString, val) = update_list!(l, Symbol(name), val)

end # module


text:
@testset "check" begin

	using TensorCast: _check!, index_store, size_store

	@check!  throw=true empty

	A = rand(3)
	@check! A[i]

	@test_throws ErrorException _check!(:(A[i,j])) # parse-time error
	@test_throws ErrorException _check!(:(A[z]))

	@check! size=true
	@check! A[j]

	B = rand(4)
	@test_throws ErrorException @check! B[j] # run-time error

	sq(x) = x .^ 2
	@check! sq(B)[j] # these should be ignored
	@check! B[(i,j)]
	@check! B[-z]
	@check! B[2]
	@test_broken length(keys(index_store)) == 2
	@check! info

    @test_throws ErrorException _check!(:(B[i,2]))
    @test_throws ErrorException _check!(:(B[i,j\k]))

    E = rand(5)
    @check! E[abc]

end
@testset "shape! reduce! cast!" begin

	# using TensorCast: _shape, _reduce, _cast

	# @check!  throw=true size=false # empty

	# A = rand(3)
	# C = rand(3)
	# @shape! A[j] = C[j]

	# @test_throws ErrorException _shape( :( A[z] = C[z] ); icheck=true)
	# @test_throws ErrorException _reduce( :( A[z] := sum(i) ), :( D[z,i] ); icheck=true)

end
@testset "einsum!" begin

    @check!  throw=true  empty

    AB = rand(4,4)
    XY = rand(4,4)

    @reduce! A[a] := sum(b) AB[a,b]
    @cast! ABC[a,b,c] := A[a] + AB[b,c]

    # @einsum! AC[a,c] := AB[a,b] * AB[b,c]
    # @test AC == AB * AB

end


text:
using LinearAlgebra, LowRankApprox, Formatting, ImageFiltering, Images, Plots, ImageMagick, Colors, TestImages, ImageView, ImageTransformations

I=load("rgb.jpg");
X=imresize(I, ratio=1/3);
(p,q)=size(X);
A=channelview(X);
R=Array{Float64}(A[1,:,:]);
G=Array{Float64}(A[2,:,:]);
B=Array{Float64}(A[3,:,:]);

using DSP, Wavelets
R1=dwt(R, wavelet(WT.db2), 1);
R2=dwt(R, wavelet(WT.db2), 2);
R3=dwt(R, wavelet(WT.db2), 3);
R4=dwt(R, wavelet(WT.db2), 4);
G1=dwt(G, wavelet(WT.db2), 1);
G2=dwt(G, wavelet(WT.db2), 2);
G3=dwt(G, wavelet(WT.db2), 3);
G4=dwt(G, wavelet(WT.db2), 4);
B1=dwt(B, wavelet(WT.db2), 1);
B2=dwt(B, wavelet(WT.db2), 2);
B3=dwt(B, wavelet(WT.db2), 3);
B4=dwt(B, wavelet(WT.db2), 4);

R1a=zeros(Float64, p, q);
R2a=zeros(Float64, p, q);
R3a=zeros(Float64, p, q);
R4a=zeros(Float64, p, q);
R1a[1:Int(p/2),1:Int(q/2)]=R1[1:Int(p/2),1:Int(q/2)];
R2a[1:Int(p/4),1:Int(q/4)]=R2[1:Int(p/4),1:Int(q/4)];
R3a[1:Int(p/8),1:Int(q/8)]=R3[1:Int(p/8),1:Int(q/8)];
R4a[1:Int(p/16),1:Int(q/16)]=R4[1:Int(p/16),1:Int(q/16)];
R1d=R1-R1a;
R2d=R2-R2a;
R3d=R3-R3a;
R4d=R4-R4a;
G1a=zeros(Float64, p, q);
G2a=zeros(Float64, p, q);
G3a=zeros(Float64, p, q);
G4a=zeros(Float64, p, q);
G1a[1:Int(p/2),1:Int(q/2)]=G1[1:Int(p/2),1:Int(q/2)];
G2a[1:Int(p/4),1:Int(q/4)]=G2[1:Int(p/4),1:Int(q/4)];
G3a[1:Int(p/8),1:Int(q/8)]=G3[1:Int(p/8),1:Int(q/8)];
G4a[1:Int(p/16),1:Int(q/16)]=G4[1:Int(p/16),1:Int(q/16)];
G1d=G1-G1a;
G2d=G2-G2a;
G3d=G3-G3a;
G4d=G4-G4a;
B1a=zeros(Float64, p, q);
B2a=zeros(Float64, p, q);
B3a=zeros(Float64, p, q);
B4a=zeros(Float64, p, q);
B1a[1:Int(p/2),1:Int(q/2)]=B1[1:Int(p/2),1:Int(q/2)];
B2a[1:Int(p/4),1:Int(q/4)]=B2[1:Int(p/4),1:Int(q/4)];
B3a[1:Int(p/8),1:Int(q/8)]=B3[1:Int(p/8),1:Int(q/8)];
B4a[1:Int(p/16),1:Int(q/16)]=B4[1:Int(p/16),1:Int(q/16)];
B1d=B1-B1a;
B2d=B2-B2a;
B3d=B3-B3a;
B4d=B4-B4a;

RR1a=idwt(R1a, wavelet(WT.db2), 1);
RR2a=idwt(R2a, wavelet(WT.db2), 2);
RR3a=idwt(R3a, wavelet(WT.db2), 3);
RR4a=idwt(R4a, wavelet(WT.db2), 4);
RR1d=idwt(R1d, wavelet(WT.db2), 1);
RR2d=idwt(R2d, wavelet(WT.db2), 2);
RR3d=idwt(R3d, wavelet(WT.db2), 3);
RR4d=idwt(R4d, wavelet(WT.db2), 4);
GG1a=idwt(G1a, wavelet(WT.db2), 1);
GG2a=idwt(G2a, wavelet(WT.db2), 2);
GG3a=idwt(G3a, wavelet(WT.db2), 3);
GG4a=idwt(G4a, wavelet(WT.db2), 4);
GG1d=idwt(G1d, wavelet(WT.db2), 1);
GG2d=idwt(G2d, wavelet(WT.db2), 2);
GG3d=idwt(G3d, wavelet(WT.db2), 3);
GG4d=idwt(G4d, wavelet(WT.db2), 4);
BB1a=idwt(B1a, wavelet(WT.db2), 1);
BB2a=idwt(B2a, wavelet(WT.db2), 2);
BB3a=idwt(B3a, wavelet(WT.db2), 3);
BB4a=idwt(B4a, wavelet(WT.db2), 4);
BB1d=idwt(B1d, wavelet(WT.db2), 1);
BB2d=idwt(B2d, wavelet(WT.db2), 2);
BB3d=idwt(B3d, wavelet(WT.db2), 3);
BB4d=idwt(B4d, wavelet(WT.db2), 4);

A1a=zeros(3,p,q);
A1d=zeros(3,p,q);
A2a=zeros(3,p,q);
A2d=zeros(3,p,q);
A3a=zeros(3,p,q);
A3d=zeros(3,p,q);
A4a=zeros(3,p,q);
A4d=zeros(3,p,q);

A1a[1,:,:]=RR1a;
A1a[2,:,:]=GG1a;
A1a[3,:,:]=BB1a;
A1d[1,:,:]=RR1d;
A1d[2,:,:]=GG1d;
A1d[3,:,:]=BB1d;
A2a[1,:,:]=RR2a;
A2a[2,:,:]=GG2a;
A2a[3,:,:]=BB2a;
A2d[1,:,:]=RR2d;
A2d[2,:,:]=GG2d;
A2d[3,:,:]=BB2d;
A3a[1,:,:]=RR3a;
A3a[2,:,:]=GG3a;
A3a[3,:,:]=BB3a;
A3d[1,:,:]=RR3d;
A3d[2,:,:]=GG3d;
A3d[3,:,:]=BB3d;
A4a[1,:,:]=RR4a;
A4a[2,:,:]=GG4a;
A4a[3,:,:]=BB4a;
A4d[1,:,:]=RR4d;
A4d[2,:,:]=GG4d;
A4d[3,:,:]=BB4d;

X1a=colorview(RGB, A1a);
X1d=colorview(RGB, A1d);
X2a=colorview(RGB, A2a);
X2d=colorview(RGB, A2d);
X3a=colorview(RGB, A3a);
X3d=colorview(RGB, A3d);
X4a=colorview(RGB, A4a);
X4d=colorview(RGB, A4d);

plot(X,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 0: Original RGB Image")
savefig("rgb_db2_gbb2_level0.png")

plot(X1a,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 1 Approximation (db2)")
savefig("rgb_db2_gbb2_level1_a.png")

plot(X1d,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 1 Detail (db2)")
savefig("rgb_db2_gbb2_level1_d.png")

plot(X2a,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 2 Approximation (db2)")
savefig("rgb_db2_gbb2_level2_a.png")

plot(X2d,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 2 Detail (db2)")
savefig("rgb_db2_gbb2_level2_d.png")

plot(X3a,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 3 Approximation (db2)")
savefig("rgb_db2_gbb2_level3_a.png")

plot(X3d,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 3 Detail (db2)")
savefig("rgb_db2_gbb2_level3_d.png")

plot(X4a,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 4 Approximation (db2)")
savefig("rgb_db2_gbb2_level4_a.png")

plot(X4d,
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 4 Detail (db2)")
savefig("rgb_db2_gbb2_level4_d.png")


text:
using DynamicalSystems

# Observational noise
x = ones(10, 2) .+ rand(Normal(0, 0.1), 10, 2)
@test add_observational_noise!(x, 10) isa typeof(x)
@test add_observational_noise!(Dataset(x), 10) isa Dataset

#####################
# Continuous systems
#####################

# Rössler-Lorenz
# --------------
@test RosslerLorenzUnidir() isa RosslerLorenzUnidir
@test ContinuousDynamicalSystem(RosslerLorenzUnidir()) isa ContinuousDynamicalSystem
@test rand(RosslerLorenzUnidir) isa RosslerLorenzUnidir
@test trajectory(RosslerLorenzUnidir(), 100) isa Dataset


text:
using DashBootstrapComponents, DashHtmlComponents

spinners = html_div([
    dbc_spinner(color="primary", type="grow"),
    dbc_spinner(color="secondary", type="grow"),
    dbc_spinner(color="success", type="grow"),
    dbc_spinner(color="warning", type="grow"),
    dbc_spinner(color="danger", type="grow"),
    dbc_spinner(color="info", type="grow"),
    dbc_spinner(color="light", type="grow"),
    dbc_spinner(color="dark", type="grow"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using LatPhysReciprocal
using Test

@testset "LatPhysReciprocal.jl" begin
    # Write your own tests here.
end


text:
module OceanGrids

using Unitful
using Unitful: °
using Interpolations
using LinearAlgebra, SparseArrays
using NearestNeighbors
using Inpaintings

abstract type OceanGrid end

"""
    OceanRectilinearGrid

Ocean rectilinear grid.
Here I use the slightly abuse the term rectilinear,
because the grid supposedly represents layers around the (curved) globe.
Here "rectilinear" means that the grid is defined entirely by just its
`lat`, `lon`, `depth`, `δlat`, `δlon`, and `δdepth` vectors.
For example, the Ocean Circulation Inverse Model (OCIM) version 1
grid is rectilinear in that sense.
"""
struct OceanRectilinearGrid <: OceanGrid
    lat          # °
    lon          # °
    depth        # m
    δlat         # °
    δlon         # °
    δdepth       # m
    lat_3D       # °
    lon_3D       # °
    depth_3D     # m
    δy           # m
    δx_3D        # m
    δy_3D        # m
    δz_3D        # m
    volume_3D    # m³
    depth_top    # m
    depth_top_3D # m
    A_2D         # m²
    wet3D::BitArray
    nlon
    nlat
    ndepth
    nboxes
end

"""
    OceanCurvilinearGrid

Ocean curvilinear grid.
The grid is curvilinear if it can be accessed using Cartesian Indices.
I.e., it maps to a rectilinear grid.
This type of grid typically cannot be inferred from only
`lat`, `lon`, `depth`, `δlat`, `δlon`, and `δdepth` vectors.
Instead, the full 3-dimensional arrays
`lat_3D`, `lon_3D`, `D`epth_3D`, `δx_3D`, `δy_3D`, and `δz_3D`
are required.
For exaample, the displaced-pole grid of the POP model is curcvilinear.
"""
struct OceanCurvilinearGrid <: OceanGrid
    lat_3D       # °
    lon_3D       # °
    depth_3D     # m
    δx_3D        # m
    δy_3D        # m
    δz_3D        # m
    volume_3D    # m³
    depth_top    # m
    depth_top_3D # m
    A_2D         # m²
    wet3D::BitArray
    nlon
    nlat
    ndepth
    nboxes
end

const TU = AbstractArray{<:Quantity}
const T1 = AbstractArray{<:Number}

"""
    OceanGrid(elon::T, elat::U, edepth::V; R=6371.0u"km") where {T<:TU, U<:TU, V<:TU}

Returns an `OceanRectilinearGrid` with boxes whose edges are defined by the
`elon`, `elat`, and `edepth` vectors.
The globe radius can be changed with the keyword `R` (default value 6371 km)
The 3D array of wet boxes, `wet3D` is true everywhere by default.
"""
function OceanGrid(elon::TU, elat::TU, edepth::TU; R=6371.0u"km")
    unitful_data, unitless_data = generate_grid_data_no_wet3D(elon, elat, edepth, R)
    #
    nlat, nlon, ndepth = unitless_data[[2,1,3]]
    wet3D = trues(nlat, nlon, ndepth)
    return OceanRectilinearGrid(unitful_data..., wet3D, unitless_data...)
end

"""
    OceanGrid(elon::T, elat::U, edepth::V, wet3D::BitArray; R=6371.0u"km")

Same as `OceanGrid(elon, elat, edepth; R)` but with specified `wet3D`.
"""
function OceanGrid(elon::TU, elat::TU, edepth::TU, wet3D::BitArray; R=6371.0u"km")
    unitful_data, unitless_data = generate_grid_data_no_wet3D(elon, elat, edepth, R)
    return OceanRectilinearGrid(unitful_data..., wet3D, unitless_data...)
end

function generate_grid_data_no_wet3D(elon::TU, elat::TU, edepth::TU, R)
    R = R |> u"m"
    edepth = edepth .|> u"m"
    elat = elat .|> °
    edon = elon .|> °
    # δ objects
    nlon, nlat, ndepth = length(elon) - 1, length(elat) - 1, length(edepth) - 1
    nboxes = nlon * nlat * ndepth
    lon = edges_to_centers(elon)
    lat = edges_to_centers(elat)
    depth = edges_to_centers(edepth)
    δlon, δlat, δdepth = diff(elon), diff(elat), diff(edepth)
    # depth objects
    δz_3D = repeat(reshape(δdepth, (1,1,ndepth)), outer=(nlat,nlon,1))
    depth_3D = repeat(reshape(depth, (1,1,ndepth)), outer=(nlat,nlon,1))
    depth_top = cumsum(δdepth) - δdepth
    depth_top_3D = repeat(reshape(depth_top, (1,1,ndepth)), outer=(nlat,nlon,1))
    # lat objects
    lat_3D = repeat(reshape(lat, (nlat,1,1)), outer=(1,nlon,ndepth))
    δy = R * δlat ./ 360°
    δy_3D = repeat(reshape(δy, (nlat,1,1)), outer=(1,nlon,ndepth))
    # lon objects
    lon_3D = repeat(reshape(lon, (1,nlon,1)), outer=(nlat,1,ndepth))
    # For δx_3D, first calculate the area and then infer the equivalent distance in meters
    A_2D = R^2 * abs.(sin.(elat[1:end-1]) - sin.(elat[2:end])) * ustrip.(δlon .|> u"rad")'
    A_3D = repeat(A_2D, outer=(1, 1, ndepth))
    δx_3D = A_3D ./ δy_3D
    # volume
    volume_3D = δx_3D .* δy_3D .* δz_3D
    unitful_data = (lat, lon, depth, δlat, δlon, δdepth, lat_3D, lon_3D, depth_3D,
                    δy, δx_3D, δy_3D, δz_3D, volume_3D, depth_top, depth_top_3D, A_2D)
    unitless_data = (nlon, nlat, ndepth, nboxes)
    return unitful_data, unitless_data
end

# Convert units in case not provided
function OceanGrid(elon::TU, elat::TU, edepth::T1; R=6371.0u"km")
    @warn "`edepth` provided without units — assuming meters"
    return OceanGrid(elon, elat, edepth * u"m"; R=R)
end
function OceanGrid(elon::T1, elat::TU, edepth::TU; R=6371.0u"km")
    @warn "`elon` provided without units — assuming degrees"
    return OceanGrid(elon * °, elat, edepth; R=R)
end
function OceanGrid(elon::TU, elat::T1, edepth::TU; R=6371.0u"km")
    @warn "`elat` provided without units — assuming degrees"
    return OceanGrid(elon, elat * °, edepth; R=R)
end
function OceanGrid(elon::T1, elat::T1, edepth::TU; R=6371.0u"km")
    @warn "`elat` and `elon` provided without units — assuming degrees"
    return OceanGrid(elon * °, elat * °, edepth; R=R)
end
function OceanGrid(elon::T1, elat::T1, edepth::T1; R=6371.0u"km")
    @warn "`elat`, `elon`, and `edepth` provided without units — assuming degrees and meters"
    return OceanGrid(elon * °, elat * °, edepth * u"m"; R=R)
end
function OceanGrid(elon::TU, elat::T1, edepth::T1; R=6371.0u"km")
    @warn "`elat` and `edepth` provided without units — assuming and meters"
    return OceanGrid(elon, elat * °, edepth * u"m"; R=R)
end
function OceanGrid(elon::T1, elat::TU, edepth::T1; R=6371.0u"km")
    @warn "`elon` and `edepth` provided without units — assuming degrees and meters"
    return OceanGrid(elon * °, elat, edepth * u"m"; R=R)
end

edges_to_centers(x::Vector) = 0.5 * (x[1:end-1] + x[2:end])
edges_to_centers(x::AbstractRange) = x[1:end-1] .+ 0.5step(x)

"""
    OceanGrid(nlat::Int, nlon::Int, ndepth::Int)

Returns a regularly spaced `OceanRectilinearGrid` with size `nlat`, `nlon`, and `ndepth`.
"""
function OceanGrid(nlat::Int, nlon::Int, ndepth::Int)
    elat = range(-90,90,length=nlat+1) * °
    elon = range(0,360,length=nlon+1) * °
    edepth = range(0,3682,length=ndepth+1) * u"m"
    return OceanGrid(elon, elat, edepth)
end

function Base.show(io::IO, g::OceanGrid)
    print("OceanGrid of size $(g.nlat)×$(g.nlon)×$(g.ndepth) (lat×lon×depth)")
end

"""
    OceanGridBox

Ocean grid box.
Each grid can be looped over, where the `OceanGridBox` are the elements of the grid.
This useful to investigate specific boxes in the grid.
"""
struct OceanGridBox
    I
    lat          # °
    lon          # °
    depth        # m
    δlat         # °
    δlon         # °
    δdepth       # m
    δx           # m
    δy           # m
    δz           # m
    volume       # m³
    depth_top    # m
    A            # m²
    iswet::Bool
end

"""
    box(g::OceanGrid, i, j, k)

Accesses the individual box of `g::OceanGrid` at index `(i,j,k)`.
Each grid can be looped over, where the `OceanGridBox` are the elements of the grid.
This useful to investigate specific boxes in the grid.
"""
function box(g::OceanGrid, i, j, k)
    return OceanGridBox(
                        CartesianIndex(i,j,k),
                        g.lat[i],
                        g.lon[j],
                        g.depth[k],
                        g.δlat[i],
                        g.δlon[j],
                        g.δdepth[k],
                        g.δx_3D[i,j,k],
                        g.δy_3D[i,j,k],
                        g.δz_3D[i,j,k],
                        g.volume_3D[i,j,k],
                        g.depth_top[k],
                        g.A_2D[i,j],
                        g.wet3D[i,j,k]
                       )
end

"""
    box(g::OceanGrid, I)

Accesses the individual box of `g::OceanGrid` at index `I`.
Each grid can be looped over, where the `OceanGridBox` are the elements of the grid.
This useful to investigate specific boxes in the grid.
"""
function box(g::OceanGrid, I)
    i,j,k = CartesianIndices((g.nlat,g.nlon,g.ndepth))[I].I
    return box(g::OceanGrid, i, j, k)
end

Base.getindex(g::OceanGrid, I) = box(g,I)

"""
    Base.size(g::OceanGrid)

Size of the grid.
"""
Base.size(g::OceanGrid) = g.nlat, g.nlon, g.ndepth

"""
    Base.length(g::OceanGrid)

Size of the grid.
"""
Base.length(g::OceanGrid) = g.nlat * g.nlon * g.ndepth

function Base.show(io::IO, b::OceanGridBox)
    println("$(wet_print(b.iswet)) OceanGridBox at $(b.I):")
    println("  location: $(round(b.lat,digits=1))N, $(round(b.lon,digits=1))E")
    println("  depth: $(round(b.depth,digits=1))")
    println("  size: $(round(b.δx |> u"km",digits=1)) × $(round(b.δy |> u"km",digits=1)) × $(round(b.δz,digits=1)) (δx × δy × δz)")
end

wet_print(iswet) = iswet ? "Wet" : "Dry"

Base.round(q::Quantity; digits=0) = round(q |> ustrip, digits=digits) * unit(q)

area(b::OceanGridBox) = b.A |> u"km^2"
volume(b::OceanGridBox) = b.volume

Base.iterate(g::OceanGrid) = box(g::OceanGrid, 1), 1
function Base.iterate(g::OceanGrid, i)
    if i == g.nboxes
        return nothing
    else
        return box(g::OceanGrid, i+1), i+1
    end
end

Base.isequal(g₁::T, g₂::T) where {T<:OceanGrid} = prod([getfield(g₁,f) == getfield(g₂,f) for f in fieldnames(T)])
Base.:(==)(g₁::T, g₂::T) where {T<:OceanGrid} = isequal(g₁, g₂)

export OceanGrid, OceanCurvilinearGrid, OceanRectilinearGrid, box, OceanGridBox

include("functions.jl")

end # module


text:
@testset "l1 regularizer" begin
    x = [1, 2, 3, 4, 5]
    reg = regularizer_l1(0.5)

    @test g(reg, x) == 7.5

    y, M = prox_αg(reg, x, 4.5)
    @test y == [0.0, 0.0, 0.75, 1.75, 2.75]
    @test M == l1Manifold{5}([0, 0, 1, 1, 1])

    @test ∇M_g(reg, M, y) == [0.0, 0.0, 0.5, 0.5, 0.5]

    @test norm(∇²M_g_ξ(reg, M, y, rand(5))) == 0

    comparison = check_regularizer_gradient_hessian(M, reg)
    for (k, slope) in [(:frgrad, 2), (:frhess, 3)]
        regressiondata = remove_small_functionvals(comparison[k])

        # Some vaues may remain for t close to 1
        if length(regressiondata) > 5
            slope, residual = get_sloperesidual(regressiondata)
            @test isapprox(slope, slope; atol = 1e-1)
            @test isapprox(residual, 0; atol = 5e-2)
        end
    end
end


text:
module GridapGeosciences
  using Gridap
  using FillArrays
  using SparseArrays
  using LinearAlgebra
  using CSV
  using DataFrames
  using WriteVTK
  using Printf
  include("GeoConstantsParameters.jl")
  include("CubedSphereTriangulations.jl")
  include("CubedSphereDiscreteModels.jl")
  include("Operators.jl")
  include("CoordinateTransformations.jl")
  include("DiagnosticTools.jl")
  include("ShallowWaterExplicit.jl")
  include("ShallowWaterRosenbrock.jl")
  include("ShallowWaterIMEX.jl")
  include("ThermalShallowWaterExplicit.jl")
  include("ThermalShallowWaterExplicit_MatAdv.jl")
  include("ShallowWaterThetaMethodFullNewton.jl")
  include("Helpers.jl")
  export rₑ, Ωₑ, g, f
  export CubedSphereDiscreteModel
  export perp,⟂
  export divergence_unit_sphere
  export laplacian_unit_sphere
  export normal_unit_sphere
  export gradient_unit_sphere
  export xyz2θϕr
  export xyz2θϕ
  export θϕ2xyz
  export spherical_to_cartesian_matrix
  export cartesian_to_spherical_matrix
  export Eₖ, Eₚ, Eₜ
  export compute_kin_to_pot!, compute_pot_to_kin!, compute_total_mass!
  export project_shallow_water_initial_conditions
  export dump_diagnostics_shallow_water!
  export shallow_water_explicit_time_stepper
  export shallow_water_rosenbrock_time_stepper
  export shallow_water_imex_time_stepper
  export thermal_shallow_water_explicit_time_stepper
  export thermal_shallow_water_mat_adv_explicit_time_stepper
  export shallow_water_theta_method_full_newton_time_stepper
  export write_to_csv, get_scalar_field_from_csv, append_to_csv, initialize_csv
  export clone_fe_function, setup_mixed_spaces, setup_and_factorize_mass_matrices, new_vtk_step

  # MPI-parallel part
  using PartitionedArrays
  using GridapDistributed
  using GridapP4est
  using GridapPETSc
  using P4est_wrapper
  include("mpi/CubedSphereDistributedDiscreteModels.jl")
end # module

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
export SeisMP,
SeisOMP,
SeisSWCGP
include("SeisMP.jl")
include("SeisOMP.jl")
include("SeisSWCGP.jl")

max_stars_repo_name: Ji-seismic/SeisTest

max_issues_repo_name: Ji-seismic/SeisTest

max_forks_repo_name: Ji-seismic/SeisTest


text:
using HPAT
#using HPAT.API.Kmeans
#HPAT.CaptureAPI.set_debug_level(3)
#HPAT.DomainPass.set_debug_level(3)
#HPAT.DataTablePass.set_debug_level(3)
#HPAT.DistributedPass.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc hpat function q26(category, item_count, num_centroids, iterations, file_name)
    store_sales = DataSource(DataTable{:ss_item_sk=Int64,:ss_customer_sk=Int64}, HDF5, file_name)
    item = DataSource(DataTable{:i_item_sk=Int64,:i_category=Int64,:i_class_id=Int64}, HDF5, file_name)

    sale_items = join(store_sales, item, :ss_item_sk==:i_item_sk, :ss_item_sk)

    sale_items = sale_items[:i_category==category]

    customer_i_class = aggregate(sale_items, :ss_customer_sk, :ss_item_count = length(:ss_item_sk),
                                                               :id1 = sum(:i_class_id==1),
                                                               :id2 = sum(:i_class_id==2),
                                                               :id3 = sum(:i_class_id==3),
                                                               :id4 = sum(:i_class_id==4),
                                                               :id5 = sum(:i_class_id==5),
                                                               :id6 = sum(:i_class_id==6),
                                                               :id7 = sum(:i_class_id==7),
                                                               :id8 = sum(:i_class_id==8),
                                                               :id9 = sum(:i_class_id==9),
                                                               :id10 = sum(:i_class_id==10),
                                                               :id11 = sum(:i_class_id==11),
                                                               :id12 = sum(:i_class_id==12),
                                                               :id13 = sum(:i_class_id==13),
                                                               :id14 = sum(:i_class_id==14),
                                                               :id15 = sum(:i_class_id==15))

    customer_i_class = customer_i_class[:ss_item_count>item_count]

    points = transpose(typed_hcat(Float64,customer_i_class[:ss_item_count], customer_i_class[:id1], customer_i_class[:id2],
    customer_i_class[:id3], customer_i_class[:id4], customer_i_class[:id5], customer_i_class[:id6],
    customer_i_class[:id7], customer_i_class[:id8], customer_i_class[:id9], customer_i_class[:id10],
    customer_i_class[:id11], customer_i_class[:id12], customer_i_class[:id13], customer_i_class[:id14], customer_i_class[:id15]))

    model = Kmeans(points, num_centroids, iterations)
    return model
end

#println(q26(1, 1, 8, 20, "test_q26.hdf5"))
println(q26(90882, 1, 8, 20, "test_q26.hdf5"))


text:
# procga_base.jl
"""
Core routine for ProcGA
    gene must be givin as Array{Int,2} rowwise.
"""

using LinearAlgebra, Statistics
using Random
using DataFrames

# abstract function ==========================
# these must be overriden by user
function listpenalty()
    []
end
# calcutblate all penalty adding its validlength
function penalty()
    return 0
end
# ============================================

function nzminimum(lst::Vector{Int})::Int
    l = lst[lst.>0]
    if length(l) > 0
        minimum(l)
    else
        zero(Int)
    end
end

# taking column
function coltake(lst::Vector{Int},idx)
    lst[idx]
end

function coltake(tbl::DVector{Int},idx)
    [coltake(v,idx) for v in tbl]
end

# return list of sum of each column
# function colsumlist(tbl::DVector{Int})
#     # [sum(coltake(tbl,i)) for i in 1:length(tbl[1])]
#     sum(tbl) # this returns list
# end

function zerotbl(tbl::DVector{Int})
    [ zeros(Int,length(tbl[1])) for x in tbl ]
end

# validation
# get valid length of job time table
function validlength(lst::Vector{Int})
    x = findlast(x->(x>0),lst)
    if x == nothing
        return 0
    else
        return x
    end
end

function validlength(tbl::DVector{Int})
    v = validlength.(tbl)
    maximum(v)
end

function listusing(lst,id)
    [ x == id ? 1 : 0 for x in lst]
end

# ============================== penalty =================================

# job switching between assignment
function listjobswitch(tbl::DVector{Int},id::Int)
    utbl = listusing.(tbl,id)
    p = zerotbl(utbl)

    v1 = coltake(utbl,1)
    for i in 2:length(utbl[1])
        v = coltake(utbl,i)
        if (sum(v)>0) & (sum(v1)>0) & (v != v1)
            cc = findall(x->x>0, v)
            for x in cc
                p[x][i] = 1
            end
        end
        v1 = v
    end
    p
end
function listjobswitch(tbl::DVector{Int})
    # take unique non-zero elements
    xx = Int[]
    [append!(xx,unique(x)) for x in tbl]
    sort!(unique!(xx))
    el = xx[xx.>0]

    p = zerotbl(tbl)
    for e in el
        p .+= listjobswitch(tbl,e)
    end
    p
end

function listcoldup(tbl::DVector{Int})
    n = length(tbl)
    pp = zerotbl(tbl)

    for i in 1:length(tbl[1])
        p = zeros(Int,length(tbl))
        v = coltake(tbl,i)
        v2 = v[v.>0]
        for x in unique(v2)
            ck = v[v.==x]
            if length(ck) > 1 # duplicate
                p .+= [ c == x ? 1 : 0 for c in v ]
            end
        end
        for j in 1:n
            pp[j][i] = p[j]
        end
    end
    pp
end

# add penalty if each column is not same as preceeding one
function listinhomogeneous(tbl::DVector{Int})
    v0 = coltake(tbl,1)
    p = zeros(Int,length(tbl[1]))
    for i in 2:length(tbl[1])
        v1 = coltake(tbl,i)
        if (v1 != v0) & (Set(v1) == Set(v0))
            p[i] = 1
        end
        v0 = v1
    end
    [p for x in tbl] # return as table form
end

function listzero(lst::Vector{Int},n::Int)
    p = zeros(Int,length(lst))
    if n < 0
        l = length(lst)
    else
        l = n
    end
    for i in 1:l
        if lst[i] == 0
            p[i] = 1
        end
    end
    p
end
function listzero(lst::Vector{Int})
    l = validlength(lst)
    listzero(lst,l)
end
function listzero(tbl::DVector{Int}, n::Int)
    [ listzero(lst,n) for lst in tbl]
end
listzero(tbl::DVector{Int}) = listzero.(tbl)

# check overuse of same job overall
function listoveruse(tbl::DVector{Int}, id::Int, lmt::Int)
    vn = listusing.(tbl,id)
    v = sum(vn)
    pl = zeros(Int,length(tbl[1]))

    used = false
    cnt = 0
    for i in 1:validlength(v)
        if v[i] > 0
            if used
                cnt += 1
            else
                used = true
                cnt = 1
            end
            pl[i] = cnt > lmt ? (cnt-lmt) : 0
        else
            if used
                used = false
                cnt = 0
            end
        end
    end
    [pl .* x for x in vn]
end

function listoveruse(tbl::DVector{Int}, lmt::Vector{Int})
    pl = listoveruse(tbl,1,lmt[1])
    for i in 2:length(lmt)
        pl .+= listoveruse(tbl,i,lmt[i])
    end
    pl
end

# check short interval between same work
function listshortinterval(tbl::DVector{Int}, id::Int, lmt::Int)
    vn = listusing.(tbl,id)
    v = sum(vn)
    pl = zeros(Int,length(tbl[1]))

    used = false
    cnt = 0
    for i in 1:validlength(v)
        if v[i] > 0
            if !used
                pl[i] = cnt > 0 ? cnt : 0
            end
            used = true
        else
            if !used
                cnt -=  1
            else
                cnt = lmt-1
            end
            used = false
        end
    end
    [pl .* x for x in vn]
end

# uses defautblt value
function listshortinterval(tbl::DVector{Int}, lmt::Vector{Int})
    pl = listshortinterval(tbl,1,lmt[1])
    for i in 2:length(lmt)
        pl .+= listshortinterval(tbl,i,lmt[i])
    end
    pl
end

# returns 1 if item is in grp
function listgroupable(lst::Vector{Int}, grp::Vector{Int})
    p = [ (x in grp) ? 1 : 0 for x in lst]
end
function listgroupable(tbl::DVector{Int}, grp::Vector{Int})
    p = [ listgroupable(lst,grp) for lst in tbl ]
end
function listgroupable(tbl::DVector{Int}, gtbl::DVector{Int})
    p = zerotbl(tbl)
    for g in gtbl
        p .+= listgroupable(tbl,g)
    end
    p
end

# filter groupable items only
function listgroupableitems(tbl::DVector{Int}, grp::DVector{Int})
    gg = listgroupable(tbl,grp)
    [ t .* g for (t,g) in zip(tbl,gg )]
end

# filter out groupable
function listnongroupableitems(tbl::DVector{Int},grp::DVector{Int})
    gg = listgroupable(tbl,grp)
    gn = listzero(gg,-1)
    [ t .* g for (t,g) in zip(tbl,gn )]
end

# sequencial order
function listbadsequence(lst::Vector{Int})
    p = zeros(Int,length(lst))
    for i in 1:(length(lst)-1)
        m = nzminimum(lst[(i+1):end])
        if (lst[i] > 0) & (lst[i] > m)
            p[i] = m > 0 ? lst[i] - m : 0
        end
    end
    p
end
function listbadsequence(tbl::DVector{Int})
    listbadsequence.(tbl)
end

# penalty for bad grouping.
# if element is appear less than n (for each group) or more than n give penalty as its difference.
function listbadgroup(tbl::DVector{Int},grp::Vector{Int}, n::Int)
    nn = min(n,length(tbl))

    gg = listgroupable(tbl,grp)
    gl = sum(gg) # column sum
    pl = [ x == 0 ? 0 : abs(x-nn) for x in gl ]
    [ g .* pl for g in gg]
end
function listbadgroup(tbl::DVector{Int},grp::DVector{Int}, n::Vector{Int})
    p = zerotbl(tbl)
    for i in 1:length(n)
        p .+= listbadgroup(tbl,grp[i],n[i])
    end
    p
end

# if jobs are located sparce it is not good
function listlaziness(lst::Vector{Int})
    p = zeros(Int,length(lst))

    # simply give 1 to zero element between first and last non-zero elemant
    i1 = findfirst(x->x>0, lst)
    i2 = findlast(x->x>0, lst)
    for i in i1:i2
        p[i] = (lst[i] > 0) ? 0 : 1
    end

    # also give penalty to such that 7,0,7,7
    xl0 = unique(lst)
    xl = xl0[xl0.>0]
    for x in xl
        id = (1:length(lst))[lst.==x]
        sort!(id)
        i0 = id[1]
        for i in 2:length(id)
            i1 = id[i]
            p[i1] = (i1 - i0) > 1 ? (i1-i0-1) : 0
            i1 = i0
        end
    end
    p
end
function listlaziness(tbl::DVector{Int})
    listlaziness.(tbl)
end

# ==================== job table editing =============================

# enlarge job length and rotate last item
function addjobcol!(tbl::DVector{Int}, n::Int)
    t = [append!(x, zeros(Int,n)) for x in tbl]
    t[end] = circshift(t[end],n)
    t
end

# randomly switch column,row of elements based on penalty list
function swapjob!(jtbl::DVector{Int})
    rw = length(jtbl)
    col = validlength(jtbl)
    plst = listpenalty(jtbl)

    c1,r1 = findmax((x->x[2]).(findmax.(plst))) # col, row of maximum

    zeropos = findall.(x->x==0, plst)
    r2 = rand(1:rw)
    c2 = rand(zeropos[r2])

    jtbl[r1][c1],jtbl[r2][c2] = jtbl[r2][c2], jtbl[r1][c1]
    jtbl
end

# swap job list based on penaly list
function swapjobrow!(jlst::Vector{Int},plst::Vector{Int})
    pos = findmax(plst)[2]
    if pos > 0
        pos2 = findall(x->x==0, plst)
        if length(pos2) > 0
            p2 = rand(pos2) # swap with zero penalty element
        else
            p2 = rand(1:length(jlst))
        end
        jlst[pos],jlst[p2] = jlst[p2],jlst[pos]
    end
    jlst
end
function swapjobrow!(tbl::DVector{Int})
    ptbl = listpenalty(tbl)
    swapjobrow!.(tbl,ptbl)
end

# sort job table in order while keeping position of 0 entity
function sortjob!(jlst::Vector{Int})
    st = sort(jlst[jlst.>0]) # sorted value
    pos = range(1,stop=length(jlst))[jlst.>0] # index of non zero

    for k in 1:length(pos)
        jlst[pos[k]] = st[k]
    end
    jlst
end
sortjob!(tbl::DVector{Int}) = sortjob!.(tbl)

# move zero column to end
function sortjobcol!(tbl::DVector{Int})
    col = length(tbl[1])
    vs = sum(tbl)
    ky = sortperm(vs, by = x->(x>0 ? x : Inf))

    for i in 1:length(tbl)
        tbl[i] = tbl[i][ky]
    end
    tbl
    # [vs,ky]
end

# condense all by skipping 0 element
function condensejob!(lst::Vector{Int})
    v = lst[lst.>0]
    vm = length(v)

    for i in 1:length(lst)
        if i > vm
            lst[i] = 0
        else
            lst[i] = v[i]
        end
    end
    lst
end
condensejob!(tbl::DVector{Int}) = condensejob!.(tbl)

# remove trailing zeros
function clipjob(jlst::Vector{Int},n)
    jlst[1:n]
end
function clipjob(jlst::Vector{Int})
    clipjob(jlst,validlength(jlst))
end
clipjob(tbl::DVector{Int},n) = clipjob.(tbl, n)
clipjob(tbl::DVector{Int}) = clipjob.(tbl, validlength(tbl))
function clipjob!(tbl::DVector{Int},n)
    tbl = clipjob(tbl,n)
end

# ============================== crossover and mutation ==============================

# crossover exchange rows of two genes at middle
function crossover(jtbl1::DVector{Int}, jtbl2::DVector{Int})
    c1 = deepcopy(jtbl1)
    c2 = deepcopy(jtbl2)
    if c1 != c2
        rows = length(c1)
        m = rows ÷ 2
        c1 = vcat(c1[1:m],c2[m+1:end])
        c2 = vcat(c2[1:m],c1[m+1:end])
    end
    return(c1,c2)
end

# swap elements between valid index span and all
# obsolete: use swapjob
# function mutatejob!(jlst::Vector{Int})
#     c1 = validlength(jlst)
#     c2 = length(jlst)

#     m1 = rand(1:c1)
#     m2 = rand(1:c2) # swap valid cell with all including zero
#     jlst[m1], jlst[m2] = jlst[m2], jlst[m1]
#     return(jlst)
# end
# mutatejob!(tbl::DVector{Int}) = mutatejob!.(tbl)

# ============================== poputblation ==============================

function populateshuffle(tbl::DVector{Int},n::Int)
    popu = [deepcopy(tbl)]
    for i in 2:n
        jt = deepcopy(tbl)
        shuffle!.(jt)
        push!(popu,jt)
    end
    popu
end

# create initial poputblation of size n from source data
function populatefrom(jlst::DVector{Int},n::Int)
    jt = deepcopy(jlst)
    popu = [jt]
    for i in 2:n
        jt = deepcopy(jt)
        swapjobrow!(jt)
        push!(popu, jt)
    end
    popu
end

# sort poputblation table by its valid length
function sortpopulation!(pptbl::Array{DVector{Int},1})
    sort!(pptbl, by = x->penalty(x))
end

# survival
# inferior genes cannot live long
function survive!(pptbl::Array{DVector{Int},1}, survival = 0.8)
    ll = Int(floor(length(pptbl) * survival))
    splice!(pptbl,ll:length(pptbl))
    pptbl
end

# new generation from parent (uni-sex generation)
function childjob(tbl::DVector{Int}, mutant = 0.05, jobswitch = false)
    # make a copy of selected gene
    c1 = deepcopy(tbl)
    if rand() < mutant
        if jobswitch
            swapjob!(c1)
        else
            swapjobrow!(c1) # change job for eachrow
        end
    end
    c1
end

# evolution

function evoluteone!(pptbl::Array{DVector{Int},1}, survival=0.8, elite=0.2, mutant=0.05, jobswitch = false; stopat = 0)
    s0 = length(pptbl)
    sortpopulation!(pptbl)
    survive!(pptbl,survival)
    s = length(pptbl)
    en = Int(round(elite*s))
    for i in 1:(s0-s)
        ic = rand(1:en)
        c = childjob(pptbl[ic],mutant,jobswitch)
        sortjobcol!(c)
        push!(pptbl,c)
    end
    v = penalty.(pptbl)
    [minimum(v),Int(round(median(v))),maximum(v)]
end

function evolution!(pptbl::Array{DVector{Int},1}, n::Int, nr=10, survival=0.8, elite=0.2, mutant=0.05, jobswitch = false; stopat = 0)
    rep = DVector{Int}()

    for i in 1:n
        dt = evoluteone!(pptbl,survival,elite,mutant, jobswitch)
        push!(rep,dt)

        if i % nr == 0
            println("i:$i => $dt")
        end
        if dt[1] <= stopat
            break
        end
    end
    rep
end

function sortevolution!(pptbl::Array{DVector{Int},1}, n::Int, nr=10, survival=0.8, elite=0.2, mutant=0.05, jobswitch = false; stopat = 0)
    rep = DVector{Int}()

    for i in 1:n
        sortjob!.(pptbl)

        dt = evoluteone!(pptbl,survival,elite,mutant, jobswitch)
        push!(rep,dt)

        if i % nr == 0
            println("i:$i => $dt")
        end
        if dt[1] <= stopat
            break
        end
    end
    rep
end

max_stars_repo_name: Yoshinobu-Ishizaki/ProcGA

max_issues_repo_name: Yoshinobu-Ishizaki/ProcGA

max_forks_repo_name: Yoshinobu-Ishizaki/ProcGA


text:
import Parameters: @with_kw

@with_kw mutable struct parameters
	gamNeg::Int64              = -457
	gL::Int64                  = 10
	gH::Int64                  = 500
	alLow::Float64             = 0.2
	alTot::Float64             = 0.4
	thetaF::Float64            = 1e-3
	al::Float64                = 0.184
	be::Float64                = 0.000402
	B::Float64                 = 0.999
	bRange::Array{Float64,2}   = permutedims(push!(collect(0.1:0.025:0.975),0.999))
	pposL::Float64             = 0
	pposH::Float64             = 0.001
	N::Int64                   = 1000
	n::Int64                   = 500
	Lf::Int64                  = 2*10^5
	rho::Float64               = 0.001
	TE::Float64                = 5.0
	diploid::Bool              = false

	NN::Int64 = 2*N
	nn::Int64 = 2*n
	thetaMidNeutral::Array{Float64,1}   = fill(0.001,nn-1)
	dac::Array{Int64,1} = [1,2,4,5,10,20,50,200,500,700]

end

function rates(;param::parameters,convolutedSamples::binomialDict,gH::Array{Int64,1},gL::Union{Array{Int64,1},Nothing},gamNeg::Array{Int64,1},theta::Union{Float64,Nothing}=nothing,rho::Union{Float64,Nothing}=nothing,shape::Float64=0.184,iterations::Int64,output::String)

	# Iterations = models to solve
	# Factor to modify input Γ(shape) parameter. Flexible Γ distribution over negative alleles
	fac     = rand(-2:0.05:2,iterations)
	afac    = @. param.al*(2^fac)

	# Deleting shape > 1. Negative alpha_x values
	idx = findall(afac .> 1)
	if !isempty(idx)
		afac[idx] = rand(afac[afac .< 1],size(idx,1))
	end

	# Random α values
	nTot    = rand(0.1:0.01:0.9,iterations)

	# Defining αW. It is possible to solve non-accounting for weak fixations
	if isnothing(gL)
		# Setting αW to 0 for all estimations
		nLow    = fill(0.0,iterations)
		# Random strong selection coefficients
		ngl     = rand(repeat([1],iterations),iterations);
	else
		# Setting αW as proportion of α
		lfac    = rand(0.0:0.05:0.9,iterations)
		nLow    = @. nTot * lfac
		# Random weak selection coefficients
		ngl     = rand(repeat(gL,iterations),iterations);
	end

	# Creating N models to iter in threads. Set N models (paramerters) and sampling probabilites (binomialDict)
	nParam  = [param for i in 1:iterations];
	nBinom  = [convolutedSamples for i in 1:iterations];

	# Random strong selection coefficients
	ngh     = rand(repeat(gH,iterations),iterations);
	# Random negative selection coefficients
	ngamNeg = rand(repeat(gamNeg,iterations),iterations);

	# Random θ on coding regions
    obsNeut        = [(sfs[:,3]./sum(sfs[:,3]) for i in 1:iterations]
    n              = 1 ./collect(1:(param.nn-1))
    expNeut        = n ./ sum(n)

	if !isnothing(theta)
		θ = fill(param.thetaMidNeutral .* obsNeut ./ expNeut,iterations)
	else
		θrng = rand(0.0005:0.0005:0.01,iterations)
		θrng = fill.(θrng,param.nn-1)
		θ = map(x -> x .* obsNeut ./ expNeut,θrng)
	end
	# Random ρ on coding regions
	if !isnothing(rho)
		ρ = fill(rho,iterations)
	else
		ρ = rand(0.0005:0.0005:0.05,iterations)
	end

	# Estimations to thread pool.
	# Allocate ouput Array
	#=out    = SharedArray{Float64,3}(size(param.bRange,2),(size(param.dac,1) *2) + 12,iterations)=#
	out    = SharedArray{Float64,3}(size(param.bRange,2),(param.nn * 2 - 2) + 12,iterations)
	@time for i in 1:iterations
		# Each iteration solve 1 model accounting all B value in param.bRange
		@inbounds out[:,:,i] = iterRates(nParam[i], nBinom[i], nTot[i], nLow[i], ngh[i], ngl[i], ngamNeg[i], afac[i], θ[i], ρ[i]);
	end
	#=param = nParam[i]; convolutedSamples=nBinom[i];alTot= nTot[i];alLow= nLow[i];gH= ngh[i];gL= ngl[i];gamNeg= ngamNeg[i];afac= afac[i];θ= θ[i]; ρ= ρ[i];obsNeut = sNeut[i]=#
	# Reducing array
	df = vcat(eachslice(out,dims=3)...);

	# Saving models and rates
	models = DataFrame(df[:,1:8],[:B,:alLow,:alTot,:gamNeg,:gL,:gH,:al,:θ])
	neut   = df[:,9:param.nn-1+8]
	sel    = df[:,(9+param.nn-1):(8+param.nn*2-2)]
	dsdn   = df[:,(end-3):end]

	n = OrderedDict{Int,Array}()
	s = OrderedDict{Int,Array}()
	for i in param.dac
		n[i] = neut[:,i]
		s[i] = sel[:,i]
	end

	#=neut   = df[:,9:(8+size(param.dac,1))]
	sel    = df[:,(9+size(param.dac,1)):(8+size(param.dac,1)*2)]
	dsdn   = df[:,(end-3):end]=#

	# Writting HDF5 file
	JLD2.jldopen(output, "a+") do file
		file[string(param.N)* "/" * string(param.n) * "/models"] = models
		file[string(param.N)* "/" * string(param.n) * "/neut"]   = neut
		file[string(param.N)* "/" * string(param.n) * "/sel"]    = sel
		file[string(param.N)* "/" * string(param.n) * "/dsdn"]   = dsdn
		file[string(param.N)* "/" * string(param.n) * "/dac"]    = param.dac
	end

	return df
end

"""
	iterRates(param::parameters,afac::Float64,bfac::Float64,alTot::Float64,alLow::Float64,divergence::Array,sfs::Array)

Estimating rates given a model for all B range.

# Arguments
 - `param::parameters`
 - `convolutedSamples::binomialDict`
 - `alTot::Float64`
 - `alLow::Float64`
 - `gH::Int64`
 - `gL::Int64`
 - `gamNeg::Int64`
 - `afac::Float64`
 - `ρ::Float64`
 - `θ::Float64`
# Output
 - `Array{Float64,2}`
"""
function iterRates(param::parameters,convolutedSamples::binomialDict,alTot::Float64,alLow::Float64,gH::Int64,gL::Int64,gamNeg::Int64,afac::Float64,θ::Float64,ρ::Float64,obsNeut::Array{Float64})

	# Creating model to solve
	# Γ distribution
	param.al    = afac; param.be = abs(afac/gamNeg); param.gamNeg = gamNeg
	# α, αW
	param.alLow = alLow; param.alTot = alTot;
	# Positive selection coefficients
	param.gH    = gH;param.gL = gL
	# Mutation rate and recomb
	param.thetaMidNeutral = θ; param.rho = ρ
	# Solving θ on non-coding region and probabilites to get α value without BGS
	param.B = 0.999
	setThetaF!(param)
	setPpos!(param)

	# Allocate array to solve the model for all B values
	#=r = zeros(size(param.bRange,2),(size(param.dac,1) * 2) + 12)=#
	r = spzeros(size(param.bRange,2),(param.nn * 2 - 2) + 12)
	for j in eachindex(param.bRange)
		# Set B value
		param.B = param.bRange[j]
		# Solve θ non-coding for the B value.
		setThetaF!(param)
		# Solve model for the B value
		@inbounds r[j,:] = gettingRates(param,convolutedSamples.bn[param.B])
	end
	return r
end

"""
	gettingRates(gammaL,gammaH,pposL,pposH,observedData,nopos)

Estimating analytical rates of fixation and polymorphism to approach α value accouting for background selection, weakly and strong positive selection. Output values will be used to sample from a Poisson distribution the total counts of polymorphism and divergence using observed data.

# Arguments
 - `param::parameters`
 - `cnvBinom::SparseMatrixCSC{Float64,Int64}`
# Returns
 - `Array{Float64,2}` containing solved model, fixation and polymorphic rates
"""
function gettingRates(param::parameters,cnvBinom::SparseMatrixCSC{Float64,Int64})

	################################################
	# Subset rates accounting for positive alleles #
	################################################

	# Fixation
	fN       = param.B*fixNeut(param)
	fNeg     = param.B*fixNegB(param,0.5*param.pposH+0.5*param.pposL)
	fPosL    = fixPosSim(param,param.gL,0.5*param.pposL)
	fPosH    = fixPosSim(param,param.gH,0.5*param.pposH)

	ds       = fN
	dn       = fNeg + fPosL + fPosH

	# Polymorphism
	neut = DiscSFSNeutDown(param,cnvBinom)
	selH = if isinf(exp(param.gH * 2))
		DiscSFSSelPosDownArb(param,param.gH,param.pposH,cnvBinom)
	else
		DiscSFSSelPosDown(param,param.gH,param.pposH,cnvBinom)
	end
	selL = DiscSFSSelPosDown(param,param.gL,param.pposL,cnvBinom)
	selN = DiscSFSSelNegDown(param,param.pposH+param.pposL,cnvBinom)
	# Cumulative rates
	tmp = cumulativeSfs(hcat(neut,selH,selL,selN),false)
	splitColumns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i in 1:size(matrix, 2));
	neut, selH, selL, selN = splitColumns(tmp)
	sel = (selH+selL)+selN

	##########
	# Output #
	##########

	analyticalValues::Array{Float64,2} = vcat(param.B,param.alLow,param.alTot,param.gamNeg,param.gL,param.gH,param.al,param.thetaMidNeutral,neut,sel,ds,dn,fPosL,fPosH)'

	return (analyticalValues)
end

################################
######    Polymorphism    ######
################################
# Expected number of polymorphism above frequency x from the standard diffusion theory
# f(x) = ∫(s) θs* (1/x(1-x)) * ( (ℯ^(Ns)*(1-ℯ^(-4Ns(1-x))) ) / (ℯ^(4Ns-1)))
# Convolote with the binomial to obtain the downsampled sfs
# E[P(x)] = ∑_x=x*→x=1 fB(x*)
############Neutral#############

"""

	DiscSFSNeutDown()

Expected rate of neutral allele frequency reduce by backgrou	nd selection. The spectrum depends on the number of individual []

```math
\\mathbb{E}[Ps_{(x)}] = \\sum{x^{*}=x}{x^{*}=1}f_{B}(x)
```

# Return:
 - `Array{Float64}`: expected rate of neutral alleles frequencies.
"""
function DiscSFSNeutDown(param::parameters,binom::SparseMatrixCSC{Float64,Int64})

	NN2 = convert(Int64,ceil(param.NN*param.B))
	# Allocating variables

	neutralSfs(i::Int64) = 1.0/(i)

	x = collect(0:NN2)
	solvedNeutralSfs = neutralSfs.(x)
	replace!(solvedNeutralSfs, Inf => 0.0)

	# subsetDict = get(param.bn,param.B,1)
	# subsetDict = binom
	out::Array{Float64,1} = param.B.*(param.thetaMidNeutral).*0.25.*(binom*solvedNeutralSfs)
	# out = @view out[2:end-1]
	# out = out[2:end-1]

	return out
end

############Positive############
# Variable gamma in function changed to gammaValue to avoid problem with exported SpecialFunctions.gamma
"""

	DiscSFSSelPosDown(gammaValue,ppos)

Expected rate of positive selected allele frequency reduce by background selection. The spectrum depends on the number of individuals.

# Arguments
 - `gammaValue::Int64`: selection strength.
 - `ppos::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected positive selected alleles frequencies.
"""
function DiscSFSSelPosDown(param::parameters,gammaValue::Int64,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if ppos == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else

		redPlus = phiReduction(param,gammaValue)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		gammaCorrected = gammaValue*param.B

		gammaExp1 = exp(gammaCorrected*2)
		gammaExp2 = exp(gammaCorrected*-2)

		positiveSfs(i::Float64,g1::Float64=gammaExp1,g2::Float64=gammaExp2,ppos::Float64=ppos) = Float64(ppos*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))

		# Original
		# ppos*0.5*(ℯ^(2*gammaCorrected)*(1-ℯ^(-2.0*gammaCorrected*(1.0-i)))/((ℯ^(2*gammaCorrected)-1.0)*i*(1.0-i)))

		# Allocating outputs
		solvedPositiveSfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solvedPositiveSfs, NaN => 0.0)

		# subsetDict = get(param.bn,param.B,1)
		# out               = (param.thetaMidNeutral)*redPlus*0.75*(subsetDict*solvedPositiveSfs)
		out::Array{Float64,1} = (param.thetaMidNeutral).*redPlus.*0.75.*(binom*solvedPositiveSfs)
		# out = out[2:end-1]

	end

	return out
end

function DiscSFSSelPosDownArb(param::parameters,gammaValue::Int64,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})

	if ppos == 0.0
		out = zeros(Float64,param.nn + 1)
		out = out[2:end-1]
	else
		redPlus = phiReduction(param,gammaValue)

		# Solving sfs
		NN2 = convert(Int64,ceil(param.NN*param.B))
		xa1  = collect(0:NN2)
		xa2  = xa1/(NN2)

		# Solving float precision performance using exponential rule. Only one BigFloat estimation.
		gammaCorrected = gammaValue*param.B
		gammaExp1::Arb = exp(Arb(gammaCorrected*2,prec=10))
		gammaExp2::Arb = exp(Arb(gammaCorrected*-2,prec=10))

		positiveSfs(i::Float64,g1::Arb=gammaExp1,g2::Arb=gammaExp2,ppos::Float64=ppos) = Float64(ppos*0.5*(g1*(1- g2^(1.0-i))/((g1-1.0)*i*(1.0-i))))
		# Allocating outputs
		solvedPositiveSfs::Array{Float64,1} = (1.0/(NN2)) * (positiveSfs.(xa2))
		replace!(solvedPositiveSfs, NaN => 0.0)
		out::Array{Float64,1} = (param.thetaMidNeutral).*redPlus.*0.75.*(binom*solvedPositiveSfs)
		# out = out[2:end-1]

	end

	return out
end

######Slightly deleterious######
"""

	DiscSFSSelNegDown(param,ppos)

Expected rate of positive selected allele frequency reduce by background selection. Spectrum drawn on a gamma DFE. It depends on the number of individuals.

# Arguments
 - `ppos::Float64`: positive selected alleles probabilty.
# Return:
 - `Array{Float64}`: expected negative selected alleles frequencies.
"""
function DiscSFSSelNegDown(param::parameters,ppos::Float64,binom::SparseMatrixCSC{Float64,Int64})
	# subsetDict = get(param.bn,param.B,1)
	solvedNegative = DiscSFSSelNeg(param,ppos)
	# out::Array = param.B*(param.thetaMidNeutral)*0.75*(subsetDict*solvedNegative)
	out = param.B.*(param.thetaMidNeutral).*0.75.*(binom*solvedNegative)
	# out = @view out[2:end-1]

	# return out[2:end-1]
	return out
end

function DiscSFSSelNeg(param::parameters,ppos::Float64)

	beta     = param.be/(1.0*param.B)
	NN2      = convert(Int64, ceil(param.NN*param.B))
	xa       = collect(0:NN2)/NN2

	solveZ   = similar(xa)

	z(x::Float64,p::Float64=ppos) = (1.0-p)*(2.0^-param.al)*(beta^param.al)*(-SpecialFunctions.zeta(param.al,x+beta/2.0) + SpecialFunctions.zeta(param.al,(2+beta)/2.0))/((-1.0+x)*x)

	solveZ   = z.(xa)

	if (solveZ[1] == Inf || isnan(solveZ[1]))
		solveZ[1] = 0.0
	end
	if (solveZ[lastindex(solveZ)] == Inf || isnan(solveZ[lastindex(solveZ)]))
		solveZ[lastindex(solveZ)] = 0.0
	end

	return 1.0/(NN2+0.0).*solveZ
end

"""
	cumulativeSfs(sfsTemp)

Changing SFS considering all values above a frequency *x*. The original asymptotic-MK approach takes Pn(x) and Ps(x) as the number of polymorphic sites at frequency *x* rather than above *x*, but this approach scales poorly as sample size increases. We define the polymorphic spectrum as stated above since these quantities trivially have the same asymptote but are less affected by changing sample size.
"""
function cumulativeSfs(sfsTemp::Array,freqs::Bool=true)

	out      = Array{Float64}(undef, size(sfsTemp,1),size(sfsTemp,2))

	if freqs
		idx = 2
	else
		idx = 1
	end

	out[1,idx:end] = sum(sfsTemp[:,idx:end],dims=1)

	@simd for i in 2:(size(sfsTemp)[1])

		#=app = view(out,i-1,:) .- view(sfsTemp,i-1,:)=#
		app = out[i-1,idx:end] .- sfsTemp[i-1,idx:end]

		if sum(app) > 0.0
			out[i,idx:end] = app
		else
			out[i,idx:end] = zeros(length(app))
		end
	end

	if freqs
		out[:,1] = sfsTemp[:,1]
	end

	return out
end

"""
	reduceSfs(sfsTemp,bins)

Function to reduce the SFS into N bins.
"""
function reduceSfs(sfsTemp::Array,bins::Int64)

	freq  = collect(0:(size(sfsTemp,1)-1))/size(sfsTemp,1)
	h1    = fit(Histogram,freq,0:(1/(bins-1)):1)
	xmap1 = StatsBase.binindex.(Ref(h1), freq)

	tmp = hcat(sfsTemp,xmap1)
	out = Array{Float64}(undef,bins-1 ,size(sfsTemp,2))
	vIter =  convert(Array,unique(xmap1)')
	@simd for i = eachindex(vIter)
		@inbounds out[i,2:end] = sum(tmp[tmp[:,end] .== i,2:end-1],dims=1)
	end

	out[:,1] = collect(1:(bins-1)) ./ bins

	return (out)
end


text:
if VERSION < v"0.7.0-DEV"
    pairs(collection) = Base.Generator(=>, keys(collection), values(collection))
end

function simple(dict::Dict)
    string("Dict(", join(map(p->string(repr(p.first),"=>", repr(p.second)), pairs(dict)), ", "), ")")
end


text:
"""
    dump_gate(blk::AbstractBlock) -> Expr

convert a gate to a YaoScript expression for serization.
The fallback is `GateTypeName(fields...)`
"""
function dump_gate end

function dump_gate(blk::ConstantGate)
    Symbol("$(typeof(blk).name)"[1:end-4])
end

function dump_gate(blk::ControlBlock)
    pairs = [:($b => C($c)) for (b, c) in zip(blk.ctrl_locs, blk.ctrl_config)]
    :($(pairs...), $(blk.locs) => $(dump_gate(blk.content)))
end

function dump_gate(blk::ChainBlock)
    Expr(:block, [dump_gate(b) for b in blk]...)
end

function dump_gate(blk::RotationGate)
    :(rot($(dump_gate(blk.block)), $(tokenize_param(blk.theta))))
end

function dump_gate(blk::TimeEvolution)
    :(time($(tokenize_param(blk.dt))) => $(dump_gate(blk.H)))
end

function dump_gate(blk::PutBlock)
    :($(blk.locs) => $(dump_gate(blk.content)))
end

function dump_gate(blk::KronBlock{N}) where {N}
    if any(x -> nqubits(x) != 1, subblocks(blk))
        error("unsupported multi-qubit in kron while dumping to Yao script.")
    end
    if length(occupied_locs(blk)) == N
        :(kron($([dump_gate(blk[i]) for i in 1:N]...)))
    else
        :(($([:($i => $(dump_gate(g))) for (i, g) in blk]...),))
    end
end

function dump_gate(blk::RepeatedBlock)
    :(repeat($(blk.locs...)) => $(dump_gate(blk.content)))
end

function dump_gate(blk::Add)
    :(+($(dump_gate.(subblocks(blk))...)))
end

function dump_gate(blk::Daggered)
    :($(dump_gate(blk.content))')
end

function dump_gate(blk::CachedBlock)
    :(cache($(dump_gate(blk.content))))
end

function dump_gate(blk::Scale)
    :($(factor(blk)) * $(dump_gate(blk.content)))
end

function dump_gate(blk::Measure{N,M}) where {M,N}
    if blk.operator == ComputationalBasis()
        MOP = :(Measure)
    else
        MOP = :(Measure($(dump_gate(blk.operator))))
    end
    locs = blk.locations isa AllLocs ? :ALL : blk.locations
    if blk.postprocess isa NoPostProcess
        :($locs => $MOP)
    elseif blk.postprocess isa ResetTo
        :($locs => $MOP => resetto($(blk.postprocess.x...)))
    elseif blk.postprocess isa RemoveMeasured
        :($locs => $MOP => remove)
    end
end

function dump_gate(blk::Subroutine)
    :(focus($(blk.locs...)) => $(dump_gate(blk.content)))
end

tokenize_param(param::Number) = param

yaotoscript(block::AbstractBlock{N}) where {N} =
    Expr(:block, :(nqubits = $N), dump_gate(block)) |> rmlines
function yaotoscript(block::ChainBlock{N}) where {N}
    ex = dump_gate(block)
    Expr(:let, Expr(:block, :(nqubits = $N), :(version = "0.6")), ex)
end
yaotofile(filename::String, block) = write(filename, string(yaotoscript(block)))

macro dumpload_fallback(blocktype, fname)
    quote
        function YaoBlocks.dump_gate(blk::$(esc(blocktype)))
            vars = [getproperty(blk, x) for x in fieldnames($(esc(blocktype)))]
            Expr(:call, $(QuoteNode(fname)), vars...)
        end
        function YaoBlocks.gate_expr(::Val{$(QuoteNode(fname))}, args, info)
            Expr(:call, $(QuoteNode(fname)), render_arg.(args, Ref(info))...)
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# convert color scheme name to color scheme object
ascolorscheme(name::Symbol) = colorschemes[name]
ascolorscheme(name::AbstractString) = ascolorscheme(Symbol(name))
ascolorscheme(scheme) = scheme

# convert value to colorant, optionally using color scheme object
ascolor(value::Symbol, scheme) = ascolor(string(value), scheme)
ascolor(value::AbstractString, scheme) = parse(Colorant, value)
ascolor(value::CategoricalValue, scheme) = scheme[levelcode(value)]
ascolor(value::Missing, scheme) = parse(Colorant, "rgba(0,0,0,0)")
ascolor(value, scheme) = convert(Colorant, value)

# --------------------------------
# PROCESS COLORS PROVIDED BY USER
# --------------------------------

# STEP 1: convert user input to colors
step1(values, scheme) = ascolor.(values, Ref(scheme))
function step1(numbers::AbstractVector{V}, scheme) where {V<:Union{Number,Missing}}
  # find indices with invalid and valid numbers
  isinvalid(v) = ismissing(v) || isnan(v)
  iinds = findall(isinvalid, numbers)
  vinds = setdiff(1:length(numbers), iinds)

  # invalid numbers are assigned full transparency
  icolors = parse(Colorant, "rgba(0,0,0,0)")

  # valid numbers are assigned colors from scheme
  vnumbers = Number.(numbers[vinds])
  vcolors  = get(scheme, vnumbers, :extrema)

  # build final vector of colors
  colors = Vector{Colorant}(undef, length(numbers))
  colors[iinds] .= icolors
  colors[vinds] .= vcolors

  colors
end

# STEP 2: add transparency to colors
step2(colors, alphas) = coloralpha.(colors, alphas)
step2(colors, alpha::Number) = coloralpha.(colors, Ref(alpha))

function process(values, scheme, alphas)
  colors = step1(values, ascolorscheme(scheme))
  step2(colors, alphas)
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using Random

"""
    RandomBranching()

Branching strategy that picks a fractional variable randomly.
"""
struct RandomBranching <: VariableBranchingRule end

function find_branching_var(rule::RandomBranching, node::Node, pool::NodePool)::Variable
    return shuffle(node.fractional_variables)[1]
end


text:
using Base.Test
using GLPlot
using GLAbstraction

using ModernGL

using GLFW
using Color
using Images
using Reactive
using ImmutableArrays
using FixedPointNumbers

N = 54

window = createdisplay(w=iround(N*11.1), h=N)
camera = OrthographicPixelCamera(window.inputs)

immutable GLGlyph{T <: Real} <: AbstractVector{T}
  glyph::T
end
Base.getindex(collection::GLGlyph, I::Integer) = I == 1 ? collection.glyph : error("Out of bounds")
Base.length{T}(::GLGlyph{T}) 		= 1
Base.length{T}(::Type{GLGlyph{T}})  = 1
Base.eltype{T}(::GLGlyph{T}) 		= T
Base.eltype{T}(::Type{GLGlyph{T}}) 	= T
Base.convert{T}(::Type{GLGlyph{T}}, x::Char) = (int(x) >= 0 && int(x) <= 256) ? GLGlyph(convert(T, x)) : error("This char: ", x, " can't be converted to GLGlyph")

test_images = {
	[float32(33) for i=1:N, j=1:N],
	[GLGlyph{UInt8}(33) for i=1:N, j=1:N],
	[Vector3{Ufixed8}(0, 1, 0) for i=1:N, j=1:N],
	[Vec4(0,1,0,1) for i=1:N, j=1:N],
	[Vec3(0,1,0)   for i=1:N, j=1:N],
	[AlphaColorValue(RGB(0f0,1f0,0f0), 1f0) for i=1:N, j=1:N],
	Images.ColorTypes.AlphaColor{RGB{Float32}, Float32}[Images.ColorTypes.AlphaColor(RGB(0f0,1f0,0f0), 1f0) for i=1:N, j=1:N],
	[AlphaColorValue(Images.ColorTypes.BGR(0f0,1f0, 0f0), 1f0) for i=1:N, j=1:N],
	"test.png",
	"test.tif",
	"test.jpg",
	"test.bmp",
}
scale_matrix 	= scalematrix(Vec3(1, 1, 1))
trans_matrix	= translationmatrix(Vec3(0))

test_renderobjects = map(test_images) do image
	global scale_matrix, trans_matrix
	obj = toopengl(Texture(image), model=scale_matrix*trans_matrix, camera=camera)
	trans_matrix *= translationmatrix(Vec3(N+1, 0, 0))
	obj
end

map(glplot, test_renderobjects)

renderloop(window)