text:
using Dates
using AxisArrays
using ClimatePlots
using ClimateTools
using PyPlot
using Test

@testset "Mapping" begin

# Shapefile test
filenc = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")

C = load(filenc, "tas")
# msk = inpolygrid(C.longrid, C.latgrid, P)
status, figh = ClimatePlots.mapclimgrid(C); @test status == true; PyPlot.close()

regions = ["World", "WorldAz", "Canada", "Quebec", "Quebec_agricole", "South_Quebec", "Americas", "NorthAmerica", "SouthAmerica", "Greenwich", "Outaouais", "Laurentides", "Estrie", "Arctic", "Antarctic", "Africa", "Europe", "Asia", "West-Asia"]

for iregion in regions
    status, figh = ClimatePlots.mapclimgrid(C, region=iregion);@test status == true; PyPlot.close()
end

# precip
C = load(filenc, "pr", data_units="mm") + 2.0

regions = ["World", "WorldAz", "Canada", "Quebec", "Quebec_agricole", "South_Quebec", "Americas", "NorthAmerica", "SouthAmerica", "Greenwich", "Outaouais", "Laurentides", "Estrie", "Arctic", "Antarctic", "Africa", "Europe", "Asia", "West-Asia"]

for iregion in regions
    status, figh = ClimatePlots.mapclimgrid(C, region=iregion);@test status == true; PyPlot.close()
end

status, figh = ClimatePlots.mapclimgrid(C, start_date=(2000, 05, 16, 12), end_date=(2000, 05, 16, 12));@test status == true; PyPlot.close()

# ua wind
filenc = joinpath(dirname(@__FILE__), "data", "sresa1b_ncar_ccsm3-example.nc")

C = load(filenc, "ua")
status, figh = ClimatePlots.mapclimgrid(C, level = 3);@test status == true; PyPlot.close() # given level
status, figh = ClimatePlots.mapclimgrid(C);@test status == true; PyPlot.close() # feeding a 4D field

# Empty maps
regions = ["World", "WorldAz", "Canada", "Quebec", "Quebec_agricole", "South_Quebec", "Americas", "NorthAmerica", "SouthAmerica", "Greenwich", "Outaouais", "Laurentides", "Estrie", "Arctic", "Antarctic", "Africa", "Europe", "Asia", "West-Asia"]

for iregion in regions
    status, figh = ClimatePlots.mapclimgrid(region=iregion);@test status == true; PyPlot.close()
end

# DUMMY MAPS
lon = collect(-180.0:180.0)
lat = collect(-90.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = DateTime(2003,1,1):Day(1):DateTime(2003,01,31)
data = randn(361, 181, 31)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.mapclimgrid(C);@test status == true; PyPlot.close()
status, figh = ClimatePlots.mapclimgrid(C, center_cs=true, caxis=[200, 250]);@test status == true; PyPlot.close()
status, figh = ClimatePlots.mapclimgrid(C, center_cs=true, surface=:pcolormesh);@test status == true; PyPlot.close()

lon = collect(-180.0:180.0)
lat = collect(-90.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = Dates.year(DateTime(2003, 01, 01)):Dates.year(1):Dates.year(DateTime(2004, 01, 01))
data = randn(361, 181, 2)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = ClimatePlots.mapclimgrid(C);@test status == true; PyPlot.close()

lon = collect(-180.0:180.0)
lat = collect(-90.0:90.0)
longrid, latgrid = ndgrid(lon, lat)
timeV = DateTime(2003, 01, 01):Day(1):DateTime(2004, 01, 01)
data = randn(361, 181, 366)
dimension_dict = Dict(["lon" => "lon", "lat" => "lat"])
varattribs = Dict(["standard_name" => "random noise"])
axisdata = AxisArray(data, Axis{:lon}(lon), Axis{:lat}(lat), Axis{:time}(timeV))
C = ClimateTools.ClimGrid(axisdata, variable = "psl", typeofvar="psl", longrid=longrid, latgrid=latgrid, dimension_dict=dimension_dict, varattribs=varattribs)

status, figh = plot(C); @test status == true; PyPlot.close()
status, figh = plot(C, xlimit=[2 5]); @test status==true;PyPlot.close()
status, figh = plot(C, label = "dummy", titlestr="dummytest", gridfig=true); @test status == true; PyPlot.close()

status, figh = hist(C); @test status == true; PyPlot.close()
status, figh = hist(C, ylimit=[2, 5]); @test status == true; PyPlot.close()
status, figh = hist(C, label = "dummy", titlestr="dummytest", gridfig=false); @test status == true; PyPlot.close()

end


text:
export gridBoxCoords, GridBox, gridCoords

"""

    GridBox{NX, NY, NZ} <: SemiMutableParameter{GridBox, Float64}

A `struct` that stores coordinates of grid points in terms of both `Vector`s and
`ParamBox`s.

≡≡≡ Field(s) ≡≡≡

`num::Int`: Total number of the grid points.

`spacing::Float64`: The length between adjacent grid points.

`box::Vector{NTuple{3, ParamBox}}`: The coordinates of grid points in terms of `ParamBox`s.

`coord::Array{Array{Float64, 1}, 1}`: The coordinates of grid points in terms of `Vector`s.

≡≡≡ Initialization Method(s) ≡≡≡

    GridBox(nGrids::NTuple{3, Int}, spacing::Real=10,
            centerCoord::Array{<:Real, 1}=[0.0,0.0,0.0];
            canDiff::Bool=true, index::Int=0) -> GridBox

Construct a general `GridBox` that doesn't have to shape as a cube. `nGrid` is a 3-element
`Tuple` that specifies the number of grids (number of grid points - 1) along 3 dimensions.
`spacing` specifies the length between adjacent grid points. `centerCoord` specifies the
geometry center coordinate of the box. `canDiff` determines whether the `ParamBox` should
be marked as differentiable. `index` defines the index number for the actual parameter:
spacing `L`, with the default value 0 it would be `L₀`.

    GridBox(nGridPerEdge::Int, spacing::Real=10,
            centerCoord::Array{<:Real, 1}=[0.0,0.0,0.0];
            canDiff::Bool=true, index::Int=0) -> GridBox

Method of generating a cubic `GridBox`. `nGridPerEdge` specifies the number of grids
(number of grid points - 1) along each dimension.`spacing` specifies the length between
adjacent grid points. `centerCoord` specifies the geometry center coordinate of the box.
`canDiff` determines whether the `ParamBox` should be marked as differentiable. `index`
defines the index number for the actual parameter: spacing `L`, with the default value 0
it would be `L₀`.
"""
struct GridBox{NX, NY, NZ} <: SemiMutableParameter{GridBox, Float64}
    num::Int
    spacing::Float64
    box::Vector{NTuple{3, ParamBox}}

    function GridBox(nGrids::NTuple{3, Int}, spacing::Real=10,
                     centerCoord::Vector{<:Real}=[0.0,0.0,0.0];
                     canDiff::Bool=true, index::Int=0)
        @assert prod(nGrids .> 0) "The number of gird of each edge should be larger than 0."
        sym = ParamList[:spacing]
        spc = spacing |> Float64
        pbRef = ParamBox(spc; canDiff, index)
        boxes = NTuple{3, ParamBox{Float64}}[]
        n = 0
        supIndex = "ᴳ"*numToSups(nGrids[1])*superscriptSym['-']*numToSups(nGrids[2])*
                   superscriptSym['-']*numToSups(nGrids[3])
        for i=0:nGrids[1], j=0:nGrids[2], k=0:nGrids[3]
            n += 1
            fX0 = L -> centerCoord[1] + (i - 0.5*nGrids[1]) * L
            fY0 = L -> centerCoord[2] + (j - 0.5*nGrids[2]) * L
            fZ0 = L -> centerCoord[3] + (k - 0.5*nGrids[3]) * L
            fXname = (ParamList[:X] |> string) * supIndex * numToSubs(n)
            fYname = (ParamList[:Y] |> string) * supIndex * numToSubs(n)
            fZname = (ParamList[:Z] |> string) * supIndex * numToSubs(n)
            fX = renameFunc(fXname, fX0)
            fY = renameFunc(fYname, fY0)
            fZ = renameFunc(fZname, fZ0)
            X = ParamBox(pbRef.data, ParamList[:X], fX, sym; canDiff, index)
            Y = ParamBox(pbRef.data, ParamList[:Y], fY, sym; canDiff, index)
            Z = ParamBox(pbRef.data, ParamList[:Z], fZ, sym; canDiff, index)
            push!(boxes, (X, Y, Z))
        end
        new{nGrids[1], nGrids[2], nGrids[3]}(prod(nGrids .+ 1), spc, boxes)
    end
end

GridBox(nGridPerEdge::Int, spacing::Real=10, centerCoord::Vector{<:Real}=[0.0,0.0,0.0];
        canDiff::Bool=true, index::Int=0) =
GridBox(fill(nGridPerEdge, 3) |> Tuple, spacing, centerCoord; canDiff, index)

"""

    gridCoords(gb::GridBox) -> Array{Array{Float64, 1}, 1}

Return the grid-point coordinates in `Vector`s given the `GriBox`.
"""
gridCoords(gb::GridBox) = [outValOf.(i) |> collect for i in gb.box]

ext: jl


text:
using Documenter
using Nextion

makedocs(
    format = :html,
    sitename = "Nextion.jl",
    pages = [
        "index.md",
        "wiring.md",
        "hardware.md",
        "constants.md",
        "draw.md",
        "objects.md",
        "controls.md",
        "invisibles.md",
    ]
)

deploydocs(
    repo = "github.com/scls19fr/Nextion.jl.git",
    julia  = "1",
    latest = "master",
    target = "build",
    deps = nothing,  # we use the `format = :html`, without `mkdocs`
    make = nothing,  # we use the `format = :html`, without `mkdocs`
)


text:
"""
Covariance Matrix Adaptation Evolution Strategy Implementation: (μ/μ_{I,W},λ)-CMA-ES

The constructor takes following keyword arguments:

- `μ`/`mu` is the number of parents
- `λ`/`lambda` is the number of offspring
- `c_1` is a learning rate for the rank-one update of the covariance matrix update
- `c_c` is a learning rate for cumulation for the rank-one update of the covariance matrix
- `c_mu` is a learning rate for the rank-``\\mu`` update of the covariance matrix update
- `c_sigma` is a learning rate for the cumulation for the step-size control
- `c_m` is the learning rate for the mean update, ``c_m \\leq 1``
- `σ0`/`sigma0` is the initial step size `σ`
- `weights` are recombination weights, if the weights are set to ``1/\\mu`` then the *intermediate* recombination is activated.
"""
struct CMAES{T} <: AbstractOptimizer
    μ::Int
    λ::Int
    c_1::T
    c_c::T
    c_μ::T
    c_σ::T
    σ₀::T
    cₘ::T
    wᵢ::Vector{T}

    function CMAES(; μ::Int=15, λ::Int=2μ, mu::Int=μ, lambda::Int=0, weights::Vector{T}=zeros(lambda),
                     c_1::Real=NaN, c_c::Real=NaN, c_mu::Real=NaN, c_sigma::Real=NaN,
                     sigma0::Real=1, c_m::Real=1) where {T}
        @assert c_m ≤ 1 "cₘ > 1"
        if lambda == 0
            lambda = μ == mu ? λ : 2*mu
            weights = zeros(lambda)
        end
        @assert length(weights) == lambda "Number of weights must be $lambda"
        new{T}(mu, lambda, c_1, c_c, c_mu, c_sigma, sigma0, c_m, weights)
    end
end
population_size(method::CMAES) = method.μ
default_options(method::CMAES) = (iterations=1500, abstol=1e-15)
summary(m::CMAES) = "($(m.μ),$(m.λ))-CMA-ES"
show(io::IO,m::CMAES) = print(io, summary(m))

mutable struct CMAESState{T,TI} <: AbstractOptimizerState
    N::Int
    μ_eff::T
    c_1::T
    c_c::T
    c_μ::T
    c_σ::T
    fitpop::Vector{T}
    C::Matrix{T}
    s::Vector{T}
    s_σ::Vector{T}
    σ::T
    wᵢ::Vector{T}
    d_σ::T
    parent::TI
    fittest::TI
end
value(s::CMAESState) = first(s.fitpop)
minimizer(s::CMAESState) = s.fittest

"""Initialization of CMA-ES algorithm state"""
function initial_state(method::CMAES, options, objfun, population)
    @unpack μ,λ,c_1,c_c,c_μ,c_σ,σ₀,cₘ,wᵢ = method
    @assert μ < λ "Offspring population must be larger then parent population"

    T = typeof(value(objfun))
    individual = first(population)
    TI = typeof(individual)
    n = length(individual)
    α_cov = 2

    # set parameters
    μ_eff = 1/sum(w->w^2,wᵢ[1:μ])
    if !(1 ≤ μ_eff ≤ μ)
        # default parameters for weighted recombination
        w′ᵢ = T[log((λ+1)/2)-log(i) for i in 1:λ]
        w′⁺ = sum(w′ᵢ[w′ᵢ .≥ 0])
        μ_eff = sum(w′ᵢ[1:μ])^2/sum(w->w^2,w′ᵢ[1:μ])
        μ⁻_eff = sum(w′ᵢ[μ+1:end])^2/sum(w->w^2,w′ᵢ[μ+1:end])
        α⁻ = -sum(w′ᵢ[w′ᵢ .< 0])
        α⁻_eff = 1 + (2*μ⁻_eff)/(μ_eff+2)
        c_1 = isnan(c_1) ? α_cov/((n+1.3)^2+μ_eff) : c_1
        c_μ = isnan(c_μ) ? min(1 - c_1, α_cov*(μ_eff-2+1/μ_eff)/((n+2)^2+α_cov*μ_eff/2)) : c_μ
        c_c = isnan(c_c) ? (4 + μ_eff/n)/(n + 4 + 2*μ_eff/n) : c_c
        c_σ = isnan(c_σ) ? (μ_eff+2)/(n+μ_eff+5) : c_σ
        α⁻_pd =  (1 - c_1 - c_μ)/(n*c_μ)
        wᵢ = [ifelse(w≥ 0, w/w′⁺, min(α⁻,α⁻_eff, α⁻_pd)*w/α⁻) for w in w′ᵢ]
    else
        c_c = isnan(c_c) ? 1/sqrt(n) : c_c
        c_σ = isnan(c_σ) ? 1/sqrt(n) : c_σ
        c_μ = isnan(c_μ) ? min(μ_eff/n^2, 1-c_1) : c_μ
        c_1 = isnan(c_1) ?  min(2/n^2, 1-c_μ) : c_1
    end

    @assert 1 ≤ μ_eff ≤ μ "μ_eff ∉ [1, μ]"
    @assert c_1 + c_μ ≤ 1 "c_1 + c_μ > 1"
    @assert c_σ < 1 "c_σ ≥ 1"
    @assert c_c ≤ 1 "c_c > 1"

    d_σ = 1 + 2 * max(0, sqrt((μ_eff-1)/(n+1))-1) + c_σ

    # setup initial state
    return CMAESState{T,TI}(n, μ_eff, c_1, c_c, c_μ, c_σ,
            fill(convert(T, Inf), μ),
            diagm(0=>ones(T,n)),
            zeros(T, n), zeros(T, n), method.σ₀, wᵢ, d_σ,
            copy(individual), copy(individual) )
end

function update_state!(objfun, constraints, state::CMAESState{T,TI}, population::AbstractVector{TI},
                       method::CMAES, itr) where {T, TI}
    μ, λ, cₘ = method.μ, method.λ, method.cₘ
    N,μ_eff,σ,w,d_σ = state.N, state.μ_eff, state.σ, state.wᵢ, state.d_σ
    c_1,c_c,c_μ,c_σ = state.c_1, state.c_c, state.c_μ, state.c_σ
    weights = view(w, 1:μ)
    E_NormN = sqrt(N)*(1-1/(4*N)+1/(21*N*N))

    parent = reshape(state.parent,N)
    parentshape = size(state.parent)

    z = zeros(T, N, λ)
    # y = zeros(T, N, λ)
    ȳ = zeros(T, N)
    offspring = Array{Vector{T}}(undef, λ)
    fitoff = fill(Inf, λ)

    B, D = try
        C = Symmetric(state.C)
        F = eigen!(Symmetric(state.C))
        F.vectors, Diagonal(sqrt.(max.(0,F.values)))
    catch ex
        @error "Break on eigen decomposition: $ex: $(state.C)"
        return true
    end

    # SqrtC = try
    #     cholesky(Symmetric(state.C)).U
    # catch ex
    #     @error "Break on Cholesky: $ex: $(state.C)"
    #     return true
    # end

    for i in 1:λ
        # offspring are generated by transforming standard normally distributed random vectors using a transformation matrix
        z[:,i] = randn(N)
        # y[:,i] =  B * D * z[:,i]
        # y[:,i] =  SqrtC * z[:,i]
        # offspring[i] = state.parent + σ * y[:,i]
        offspring[i] = apply!(constraints, parent + σ * B * D * z[:,i])
        fitoff[i] = value(objfun, reshape(offspring[i],parentshape...)) # Evaluate fitness
    end
    penalty!(fitoff, constraints, offspring)

    # Select new parent population
    idx = sortperm(fitoff)
    for i in 1:μ
        o = offspring[idx[i]]
        population[i] = reshape(o,parentshape...)
        ȳ .+=(o.-parent).*(w[i]/σ)
        state.fitpop[i] = fitoff[idx[i]]
    end

    # y_λ = view(y,:,idx)
    z_λ = view(z,:,idx)
    z̄ = vec(sum(z_λ[:,1:μ].*weights', dims=2))
    parent += (cₘ*σ).*ȳ  #  forming recombinant perent for next generation
    state.s_σ = (1 - c_σ)*state.s_σ + sqrt(μ_eff*c_σ*(2 - c_σ))*(B*z̄)
    # state.s_σ = (1 - c_σ)*state.s_σ + sqrt(μ_eff*c_σ*(2 - c_σ))*(SqrtC*ȳ)
    state.σ = σ * exp((c_σ/d_σ)*(norm(state.s_σ)/E_NormN - 1))
    h_σ = norm(state.s_σ)/sqrt(1 - (1-c_σ)^(2*(itr+1))) < (1.4+2/(N+1))*E_NormN
    state.s = (1 - c_c)*state.s + (h_σ*sqrt(μ_eff*c_c*(2 - c_c)))*ȳ
    # perform rank-one update
    rank_1 = c_1.*state.s*state.s'
    if !h_σ # correction
        rank_1 += (c_1*c_c*(2 - c_c)).*state.C
    end
    # perform rank-μ update
    rank_μ = c_μ*sum( (w ≥ 0 ? 1 : N/norm(B*zi)^2)*w*(zi*zi') for (w,zi) in zip(w, eachcol(z_λ)) )
    # rank_μ = c_μ*sum( (w ≥ 0 ? 1 : N/norm(SqrtC*yi)^2)*w*(yi*yi') for (w,yi) in zip(w, eachcol(y_λ)) )
    # update covariance
    state.C = (1 - c_1 - c_μ*sum(w)).*state.C + rank_1 + rank_μ

    state.fittest = population[1]
    state.parent = reshape(parent,parentshape...)

    return false
end


text:
module MadNLPSchwarz

import ..MadNLPGraph:
    @kwdef, Logger, @debug, @warn, @error,
    default_linear_solver,SparseMatrixCSC, SubVector, StrideOneVector, get_cscsy_view, nnz,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    AbstractOptions, AbstractLinearSolver, set_options!,
    MonolevelPartition, MonolevelStruc, BilevelPartition, BilevelStruc,
    expand!, get_current_V, get_current_size, get_full_size,
    EmptyLinearSolver, introduce, factorize!, solve!, improve!, is_inertia, inertia,
    set_blas_num_threads, blas_num_threads, @blas_safe_threads, @sprintf

const INPUT_MATRIX_TYPE = :csc

@kwdef mutable struct Options <: AbstractOptions
    schwarz_num_parts_upper::Int = 0
    schwarz_num_parts::Int = 2
    schwarz_subproblem_solver::Module
    schwarz_fully_improve_subproblem_solver::Bool=true
    schwarz_max_expand_factor::Int = 4
    schwarz_custom_partition::Bool = false
    schwarz_part_upper::Vector{Int} = Int[]
    schwarz_part::Vector{Int} = Int[]
end

mutable struct SolverWorker
    struc::Union{MonolevelStruc,BilevelStruc}

    max_size::Float64
    overlap_increase_factor::Float64
    restrictor::Vector{Int}

    csc::SparseMatrixCSC{Float64,Int32}
    csc_view::SubVector{Float64}

    M::AbstractLinearSolver

    p::AbstractVector{Float64}
    q::AbstractVector{Float64}
    q_restricted::AbstractVector{Float64}
end

mutable struct Solver <: AbstractLinearSolver
    partition::Union{MonolevelPartition,BilevelPartition}

    csc::SparseMatrixCSC{Float64}
    inds::Vector{Int}

    p::Vector{Float64}
    w::Vector{Float64}
    sws::Vector{SolverWorker}
    opt::Options
    logger::Logger
end

function Solver(csc::SparseMatrixCSC{Float64};
                option_dict::Dict{Symbol,Any}=Dict(),
                opt=Options(schwarz_subproblem_solver=default_linear_solver()),logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs...)
    if string(opt.schwarz_subproblem_solver) == "MadNLP.Mumps"
        @warn(logger,"When Mumps is used as a subproblem solver, Schwarz is run in serial.")
        @warn(logger,"To use parallelized Schwarz, use Ma27 or Ma57.")
    end

    inds = collect(1:nnz(csc))

    p = zeros(csc.n)
    w = zeros(csc.n)

    if opt.schwarz_num_parts_upper == 0
        partition = MonolevelPartition(
            csc,opt.schwarz_custom_partition ? opt.schwarz_part : Int32[],opt.schwarz_num_parts)
    else
        partition = BilevelPartition(
            csc,
            opt.schwarz_custom_partition ? opt.schwarz_part : Int32[],opt.schwarz_num_parts,
            opt.schwarz_custom_partition ? opt.schwarz_part_upper : Int32[],opt.schwarz_num_parts_upper)
    end

    sws=Vector{SolverWorker}(undef,opt.schwarz_num_parts_upper == 0 ? opt.schwarz_num_parts : opt.schwarz_num_parts_upper)
    copied_option_dict = copy(option_dict)
    @blas_safe_threads for k=1:length(sws)
        sws[k] = SolverWorker(
            partition,csc,inds,p,k,opt.schwarz_max_expand_factor,
            opt.schwarz_subproblem_solver,opt.schwarz_fully_improve_subproblem_solver,
            logger,k==1 ? option_dict : copy(copied_option_dict))
    end

    saturation = maximum(sws[k].csc.n/csc.n*100 for k=1:length(sws))
    @debug(logger,@sprintf("overlap size initialized with %3d%% saturation.\n",saturation))

    return Solver(partition,csc,inds,p,w,sws,opt,logger)
end

function SolverWorker(
    partition,csc::SparseMatrixCSC{Float64},inds::Vector{Int},p::Vector{Float64},
    k::Int,max_expand_factor::Int,SubproblemSolverModule::Module,fully_improve_subproblem_solver::Bool,
    logger::Logger,option_dict::Dict{Symbol,Any})
    struc= partition isa MonolevelPartition ?
        MonolevelStruc(partition,k) : BilevelStruc(partition,k)
    overlap_increase_factor =
        (get_full_size(partition)/get_current_size(struc))^(1/max_expand_factor)
    max_size = get_current_size(struc)*(overlap_increase_factor)
    p = view(p,copy(get_current_V(struc)))
    restrictor = 1:length(get_current_V(struc))
    expand!(struc,partition,max_size)
    csc,csc_view = get_cscsy_view(csc,get_current_V(struc),inds=inds)
    if csc.n == 0
        M = EmptyLinearSolver()
        @warn(logger,"empty subproblem at partition $k")
    else
        M = SubproblemSolverModule.Solver(csc;option_dict = option_dict,logger=logger)
    end
    fully_improve_subproblem_solver && while improve!(M) end # starts with fully improved
    q = Vector{Float64}(undef,length(get_current_V(struc)))
    q_restricted = view(q,restrictor)

    return SolverWorker(struc,max_size,overlap_increase_factor,restrictor,csc,csc_view,M,p,q,q_restricted)
end

function factorize!(M::Solver)
    @blas_safe_threads for k=1:length(M.sws)
        M.sws[k].csc.nzval.=M.sws[k].csc_view
        factorize!(M.sws[k].M)
    end
    return M
end

is_maximal_overlap(M::Solver) = all(sw.max_size>=M.csc.n for sw in M.sws)

function improve!(M::Solver)
    is_maximal_overlap(M) && (@debug(M.logger,"improve quality failed.\n");return false)

    @blas_safe_threads for k=1:length(M.sws)
        M.sws[k].max_size *= M.sws[k].overlap_increase_factor
        expand!(M.sws[k].struc,M.partition,M.sws[k].max_size)
        M.sws[k].csc,M.sws[k].csc_view = get_cscsy_view(M.csc,get_current_V(M.sws[k].struc);inds=M.inds)
        if M.sws[k].csc.n == 0
            M.sws[k].M = EmptyLinearSolver()
        else
            M.sws[k].M = M.opt.schwarz_subproblem_solver.Solver(M.sws[k].csc;opt=M.sws[k].M.opt)
        end
        M.opt.schwarz_fully_improve_subproblem_solver && while improve!(M.sws[k].M) end
        resize!(M.sws[k].q,length(get_current_V(M.sws[k].struc)))
    end

    saturation = maximum(M.sws[k].csc.n/M.csc.n*100 for k=1:length(M.sws));
    saturation == 100. ?
        @warn(M.logger,@sprintf("overlap is maximally saturated")) :
        @debug(M.logger,@sprintf("overlap size increased to %3d%% saturation.\n",saturation))

    return true
end

is_inertia(::Solver)=false
function inertia(M::Solver)
    throw(InertiaException)
end

function solve!(M::Solver,x::AbstractVector{Float64})
    @blas_safe_threads for k=1:length(M.sws)
        M.sws[k].q.=view(x,get_current_V(M.sws[k].struc))
        solve!(M.sws[k].M,M.sws[k].q)
        M.sws[k].p.=M.sws[k].q_restricted
    end
    x.=M.p
end

introduce(M::Solver)="schwarz equipped with "*introduce(M.sws[1].M)

end # module


text:
module NetCDF
include("netcdf_c_wrappers.jl")
import Base.show
export show,NcDim,NcVar,NcFile,new,ncread,ncwrite,nccreate,ncsync,ncinfo,ncclose,ncputatt
#Some constants

jltype2nctype={Int16=>NC_SHORT,
               Int32=>NC_INT,
               Float32=>NC_FLOAT,
               Float64=>NC_DOUBLE}

type NcDim
  ncid::Integer
  dimid::Integer
  varid::Integer
  name::String
  dimlen::Integer
  vals::AbstractArray
  atts::Dict{Any,Any}
end
NcDim(name::String,vals::Union(AbstractArray,Number),atts::Dict{Any,Any})=NcDim(-1,-1,-1,name,length(vals),vals,atts)
NcDim(name::String,vals::Union(AbstractArray,Number))=NcDim(name,vals,{"units"=>"unknown"})

type NcVar
  ncid::Integer
  varid::Integer
  ndim::Integer
  natts::Integer
  nctype::Integer
  name::String
  dimids::Array{}
  dim::Array{NcDim}
  atts::Dict{Any,Any}
end
function NcVar(name::String,dimin,atts::Dict{Any,Any},jltype::Type)
    dim=[dimin]
    return NcVar(-1,-1,length(dim),length(atts),jltype2nctype[jltype],name,Array(Int,length(dim)),dim,atts)
end

type NcFile
  ncid::Integer
  nvar::Integer
  ndim::Integer
  ngatts::Integer
  vars::Dict{String,NcVar}
  dim::Dict{String,NcDim}
  gatts::Dict{Any,Any}
  nunlimdimid::Integer
  name::String
  omode::Uint16
  in_def_mode::Bool
end
NcFile(ncid::Integer,nvar::Integer,ndim::Integer,ngatts::Integer,vars::Dict{String,NcVar},dim::Dict{String,NcDim},gatts::Dict{Any,Any},nunlimdimid::Integer,name::String,omode::Uint16)=NcFile(ncid,nvar,ndim,ngatts,vars,dim,gatts,nunlimdimid,name,omode,false)

include("netcdf_helpers.jl")

global currentNcFiles=Dict{String,NcFile}()

# Read block of data from file
function readvar(nc::NcFile,varname::String,start::Array,count::Array)
  ncid=nc.ncid
  varid=nc.vars[varname].varid
  start=int64(start)-1
  count=int64(count)
  @assert nc.vars[varname].ndim==length(start)
  @assert nc.vars[varname].ndim==length(count)
  NC_VERBOSE ? println(keys(nc.vars)) : nothing
  for i = 1:length(count)
    count[i]= count[i]>0 ? count[i] : nc.vars[varname].dim[i].dimlen
  end
  #Determine size of Array
  p=1
  for i in count
    p=p*i
  end
  count=count[length(count):-1:1]
  start=start[length(start):-1:1]
  NC_VERBOSE ? println("$ncid $varid $p $start $count") : nothing
  if nc.vars[varname].nctype==NC_DOUBLE
    retvalsa=Array(Float64,p)
    _nc_get_vara_double_c(ncid,varid,start,count,retvalsa)
  elseif nc.vars[varname].nctype==NC_FLOAT
    retvalsa=Array(Float32,p)
    _nc_get_vara_float_c(ncid,varid,start,count,retvalsa)
  elseif nc.vars[varname].nctype==NC_INT
    retvalsa=Array(Int32,p)
    _nc_get_vara_int_c(ncid,varid,start,count,retvalsa)
  elseif nc.vars[varname].nctype==NC_SHORT
    retvalsa=Array(Int32,p)
    _nc_get_vara_int_c(ncid,varid,start,count,retvalsa)
  elseif nc.vars[varname].nctype==NC_CHAR
    retvalsa=Array(Uint8,p)
    _nc_get_vara_text_c(ncid,varid,start,count,retvalsa)
  end
  NC_VERBOSE ? println("Successfully read from file ",ncid) : nothing
  if length(count)>1
    return reshape(retvalsa,ntuple(length(count),x->count[length(count)-x+1]))
  else
    return retvalsa
  end
end
function readvar(nc::NcFile,varid::Integer,start,count)
  va=getvarbyid(nc,varid)
  va == nothing ? error("Error: Variable $varid not found in $(nc.name)") : return readvar(nc,va.name,start,count)
end
function readvar(nc::NcFile,varid::NcVar,start,count)
  return readvar(nc,varid.varid,start,count)
end

function putatt(ncid::Integer,varid::Integer,atts::Dict)
  for a in atts
    name=a[1]
    val=a[2]
    _nc_put_att(ncid,varid,name,val)
  end
end

function putatt(nc::NcFile,varname::String,atts::Dict)
  varid = haskey(nc.vars,varname) ? nc.vars[varname].varid : NC_GLOBAL
  chdef=false
  if (!nc.in_def_mode)
    chdef=true
    _nc_redef_c(nc.ncid)
  end
  putatt(nc.ncid,varid,atts)
  chdef ? _nc_enddef_c(nc.ncid) : nothing
end
function ncputatt(nc::String,varname::String,atts::Dict)
  nc= haskey(currentNcFiles,abspath(nc)) ? currentNcFiles[abspath(nc)] : open(nc,NC_WRITE)
  if (nc.omode==NC_NOWRITE)
    close(nc)
    println("reopening file in WRITE mode")
    open(fil,NC_WRITE)
  end
  putatt(nc,varname,atts)
end

function putvar(nc::NcFile,varname::String,start::Array,vals::Array)
  ncid=nc.ncid
  haskey(nc.vars,varname) ? nothing : error("No variable $varname in file $nc.name")
  @assert nc.vars[varname].ndim==length(start)
  coun=size(vals)
  count=Array(Int64,length(coun))
  start=int64(start)-1
  #Determine size of Array
  p=1
  for i in 1:length(coun)
    p=p*coun[i]
    count[i]=coun[i]
  end
  count=count[length(count):-1:1]
  start=start[length(start):-1:1]
  NC_VERBOSE ? println("$ncid $varname $p $count ",nc.vars[varname].nctype) : nothing
  #x=reshape(vals,p)
  x=vals
  varid=nc.vars[varname].varid
  if nc.vars[varname].nctype==NC_DOUBLE
    _nc_put_vara_double_c(ncid,varid,start,count,x)
  elseif nc.vars[varname].nctype==NC_FLOAT
    _nc_put_vara_double_c(ncid,varid,start,count,x)
  elseif nc.vars[varname].nctype==NC_INT
    _nc_put_vara_int_c(ncid,varid,start,count,x)
  elseif nc.vars[varname].nctype==NC_SHORT
    _nc_put_vara_int_c(ncid,varid,start,count,x)
  elseif nc.vars[varname].nctype==NC_CHAR
    _nc_put_vara_text_c(ncid,varid,start,count,x)
  end
  NC_VERBOSE ? println("Successfully wrote to file ",ncid) : nothing
end
function putvar(nc::NcFile,varname::String,vals::Array)
  start=ones(length(size(vals)))
  putvar(nc,varname,start,vals)
end

# Function to synchronize all files with disk
function ncsync()
  for ncf in currentNcFiles
    id=ncf[2].ncid
    _nc_sync_c(int32(id))
  end
end

function sync(nc::NcFile)
  id=nc.ncid
  _nc_sync_c(int32(id))
end

#Function to close netcdf files
function ncclose(fil::String)
  if (haskey(currentNcFiles,abspath(fil)))
    close(currentNcFiles[abspath(fil)])
  end
end
function ncclose()
  for f in keys(currentNcFiles)
    ncclose(f)
  end
end

function create(name::String,varlist::Union(Array{NcVar},NcVar))
  ida=Array(Int32,1)
  vars=Dict{String,NcVar}();
  #Create the file
  _nc_create_c(name,NC_CLOBBER,ida);
  id=ida[1];
  # Unify types
  if (typeof(varlist)==NcVar)
    varlist=[varlist]
  end
  # Collect Dimensions
  dims=Set{NcDim}();
  for v in varlist
    for d in v.dim
      add!(dims,d);
    end
  end
  nunlim=0;
  ndim=int32(length(dims));
  #Create Dimensions in the file
  dim=Dict{String,NcDim}();
  for d in dims
    dima=Array(Int32,1);
    NC_VERBOSE? println("Dimension length ", d.dimlen) : nothing
    _nc_def_dim_c(id,d.name,d.dimlen,dima);
    d.dimid=dima[1];
    dim[d.name]=d;
    #Create dimension variable
    varida=Array(Int32,1)
    dumids=[copy(d.dimid)]
    _nc_def_var_c(id,d.name,NC_DOUBLE,1,dumids,varida)
    putatt(id,d.dimid,d.atts)
    d.varid=varida[1]
    dd=Array(NcDim,1)
    dd[1]=d
    vars[d.name]=NcVar(id,varida[1],1,length(d.atts),NC_DOUBLE,d.name,[d.dimid],dd,d.atts)
  end
  # Create variables in the file
  for v in varlist
    i=1
    for d in v.dim
      v.dimids[i]=d.dimid
      i=i+1
    end
    vara=Array(Int32,1);
    dumids=int32(v.dimids)
    NC_VERBOSE ? println(dumids) : nothing
    _nc_def_var_c(id,v.name,v.nctype,v.ndim,int32(dumids[v.ndim:-1:1]),vara);
    putatt(id,v.varid,v.atts)
    v.varid=vara[1];
    vars[v.name]=v;
  end
  # Leave define mode
  _nc_enddef_c(id)
  #Write dimension variables
  for d in dims
    #Write dimension variable
    y=float64(d.vals)
    diml=d.dimlen
    _nc_put_vara_double_c(id,d.varid,[0],[diml],y)
  end
  #Create the NcFile Object
  nc=NcFile(id,length(vars),ndim,0,vars,dim,Dict{Any,Any}(),0,name,NC_WRITE)
  currentNcFiles[abspath(nc.name)]=nc
end

function vardef(fid::Integer,v::NcVar)
    _nc_redef(ncid)
    i=1
    for d in v.dim
      v.dimids[i]=d.dimid
      i=i+1
    end
    vara=Array(Int32,1);
    dumids=int32(v.dimids)
    NC_VERBOSE? println(dumids) : nothing
    _nc_def_var_c(id,v.name,v.nctype,v.ndim,int32(dumids[v.ndim:-1:1]),vara);
    v.varid=vara[1];
    vars[v.name]=v;
end

function close(nco::NcFile)
  #Close file
  _nc_close_c(nco.ncid)
  delete!(currentNcFiles,abspath(nco.name))
  NC_VERBOSE? println("Successfully closed file ",nco.ncid) : nothing
  return nco.ncid
end

function open(fil::String,omode::Uint16=NC_NOWRITE; readdimvar=false)
  # Open netcdf file
  ncid=_nc_op(fil,omode)
  NC_VERBOSE ? println(ncid) : nothing
  #Get initial information
  (ndim,nvar,ngatt,nunlimdimid)=_ncf_inq(ncid)
  NC_VERBOSE ? println(ndim,nvar,ngatt,nunlimdimid) : nothing
  #Create ncdf object
  ncf=NcFile(ncid,nvar-ndim,ndim,ngatt,Dict{String,NcVar}(),Dict{String,NcDim}(),Dict{Any,Any}(),nunlimdimid,abspath(fil),omode)
  #Read global attributes
  ncf.gatts=_nc_getatts_all(ncid,NC_GLOBAL,ngatt)
  #Read dimensions
  for dimid = 0:ndim-1
    (name,dimlen)=_nc_inq_dim(ncid,dimid)
    ncf.dim[name]=NcDim(ncid,dimid,-1,name,dimlen,[1:dimlen],Dict{Any,Any}())
  end
  #Read variable information
  for varid = 0:nvar-1
    (name,nctype,dimids,natts,vndim,isdimvar)=_ncv_inq(ncf,varid)
    if (isdimvar)
      ncf.dim[name].varid=varid
    end
    atts=_nc_getatts_all(ncid,varid,natts)
    vdim=Array(NcDim,length(dimids))
    i=1;
    for did in dimids
      vdim[i]=ncf.dim[getdimnamebyid(ncf,did)]
      i=i+1
    end
    ncf.vars[name]=NcVar(ncid,varid,vndim,natts,nctype,name,int(dimids[vndim:-1:1]),vdim[vndim:-1:1],atts)
  end
  readdimvar==true ? _readdimvars(ncf) : nothing
  currentNcFiles[abspath(ncf.name)]=ncf
  return ncf
end

# Define some high-level functions
# High-level functions for writing data to files
function ncread(fil::String,vname::String,start::Array,count::Array)
  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  x=readvar(nc,vname,start,count)
  return x
end
function ncread(fil::String,vname::String,ran...)
  s=ones(length(ran))
  c=ones(length(ran))
  for i in 1:length(ran)
    typeof(ran[i])<:Range1 ? nothing : error("Expected range as input for reading netcdf variable")
    s[i]=int(ran[i][1])
    c[i]=int(length(ran[i]))
  end
  return ncread(fil,vname,s,c)
end
function ncread(fil::String,vname::String)
  NC_VERBOSE ? println(vname) : nothing
  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  s=ones(Int,nc.vars[vname].ndim)
  c=s*(-1)
  return ncread(fil,vname,s,c)
end

function ncinfo(fil::String)
  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil)
  return(nc)
end

#High-level functions for writing data to a file
function ncwrite(x,fil::String,vname::String,start::Array)

  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,NC_WRITE)
  if (nc.omode==NC_NOWRITE)
    close(nc)
    println("reopening file in WRITE mode")
    open(fil,NC_WRITE)
  end
  putvar(nc,vname,start,x)
end

function ncwrite(x,fil::String,vname::String)
  nc= haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,NC_WRITE)
  if (nc.omode==NC_NOWRITE)
    close(nc)
    println("reopening file in WRITE mode")
    open(fil,NC_WRITE)
  end
  start=ones(nc.vars[vname].ndim)
  putvar(nc,vname,start,x)
end

#High-level function for creating files and variables
#
# if the file does not exist, it will be created
# if the file already exists, the variable will be added to the file

function nccreate(fil::String,varname::String,atts::Dict,dims...)
  # Checking dims argument for correctness
#  dim=parsedimargs(dims)
  dim = [a for a in dims]
  # to be done
  # open the file
  # create the NcVar object
  v=NcVar(varname,dim,atts,Float64)
  # Test if the file already exists
  if (isfile(fil))
    nc=haskey(currentNcFiles,abspath(fil)) ? currentNcFiles[abspath(fil)] : open(fil,NC_WRITE)
    if (nc.omode==NC_NOWRITE)
      close(nc)
      println("reopening file in WRITE mode")
      open(fil,NC_WRITE)
    end
    haskey(nc.vars,varname) ? error("Variable $varname already exists in file fil") : nothing
    # Check if dimensions exist, if not, create
    i=1
    # Remember if dimension was created
    dcreate=Array(Bool,length(dim))
    for d in dim
      did=_nc_inq_dimid(nc.ncid,d.name)
      if (did==-1)
        dima=Array(Int32,1);
        #_nc_redef_c(nc.ncid)
        if (!nc.in_def_mode)
          _nc_redef_c(nc.ncid)
          nc.in_def_mode=true
        end
        _nc_def_dim_c(nc.ncid,d.name,d.dimlen,dima);
        d.dimid=dima[1];
        v.dimids[i]=d.dimid;
        dcreate[i] = true
      else
        dcreate[i] = false
        v.dimids[i]=nc.dim[d.name].dimid;
      end
      i=i+1
    end
    # Create variable
    vara=Array(Int32,1);
    dumids=int32(v.dimids)
    if (!nc.in_def_mode)
          _nc_redef_c(nc.ncid)
          nc.in_def_mode=true
    end
    _nc_def_var_c(nc.ncid,v.name,v.nctype,v.ndim,int32(dumids[v.ndim:-1:1]),vara);
    v.varid=vara[1];
    nc.vars[v.name]=v;
    if (nc.in_def_mode)
          _nc_enddef_c(nc.ncid)
          nc.in_def_mode=false
    end
    i=1
    for d in dim
      if (dcreate[i])
        ncwrite(d.vals,fil,d.name)
      end
     i=i+1
    end
  else
    println(v.dim)
    nc=create(fil,v)
    for d in dim
      ncwrite(d.vals,fil,d.name)
    end
  end
end

function show(io::IO,nc::NcFile)
  println(io,"")
  println(io,"##### NetCDF File #####")
  println(io,"")
  println(io,nc.name)
  println(io,"")
  println(io,"##### Dimensions #####")
  println(io,"")
  @printf(io,"%15s %8s","Name","Length\n")
  println("-------------------------------")
  for d in nc.dim
    @printf(io,"%15s %8d\n",d[2].name,d[2].dimlen)
  end
  println(io,"")
  println(io,"##### Variables #####")
  println(io,"")
  @printf(io,"%20s    %s","Name","Dimensions\n")
  println("---------------------------------------------------------------")
  for v in nc.vars
    @printf(io,"%20s    ",v[2].name)
    for d in v[2].dim
      @printf(io,"%s, ",d.name)
    end
    @printf(io,"\n")
  end
  println(io,"")
  println(io,"##### Attributes #####")
  println(io,"")
  @printf(io,"%20s %20s %20s\n","Variable","Name","Value")
  println("---------------------------------------------------------------")
  for a in nc.gatts
    an=string(a[1])
    av=string(a[2])
    an=an[1:min(length(an),38)]
    av=av[1:min(length(av),38)]
    @printf(io,"%20s %20s %40s\n","global",an,av)
  end
  for v in nc.vars
    for a in v[2].atts
      an=string(a[1])
      av=string(a[2])
      vn=
      an=an[1:min(length(an),38)]
      av=av[1:min(length(av),38)]
      @printf(io,"%20s %20s %40s\n",v[2].name,an,av)
    end
  end

end

end # Module


text:
#!/usr/bin/env julia

using SimradRaw
using SimradEK60TestData
using Test

# EK60 file

_datagrams = collect(datagrams(EK60_SAMPLE))

@test length(_datagrams) == 2860

# EK80 file

# _datagrams = collect(datagrams(EK80_SAMPLE))

# @test length(_datagrams) == 461


text:
using BinaryBuilder

name = "Linux"
version = v"5.3"
# Collection of sources required to build linux
sources = [
    "https://gitlab.com/virtio-fs/linux.git" =>
    "980bfa3c485ad975383f1f36fef25eafb8eb2ec7",
    "./sources",
]

script = raw"""
# Disable a lot of this stuff, since kernel building is weird
for tool in LD AS AR CC CXX FC RANLIB READELF STRIP OBJDUMP OBJCOPY NM LIPO LIBTOOL; do
    unset ${tool}
    unset "HOST${tool}"
    unset "BUILD_${tool}"
    unset "${tool}_FOR_BUILD"
done

cd ${WORKSPACE}/srcdir/linux/
export PATH=/usr/bin:$PATH

# Use our generated config
mv ../linuxkernel.config arch/x86/configs/binarybuilder_defconfig

# Add necessary packages
apk add gcc openssl-dev libelf-dev musl-dev bc linux-headers

# Build it
make binarybuilder_defconfig
make -j${nproc}

# Copy out the stuff we care about
cp vmlinux $prefix
cp arch/x86/boot/bzImage $prefix
"""

products = [
    ExecutableProduct("vmlinux", :vmlinux, ".")
]

platforms = [
    Platform("x86_64", "linux"; libc="glibc")
]

dependencies = [
]

# Build the given platforms using the given sources
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct Detector
    deadtime::Float64
    resolution::Float64
    jitter::Float64
    efficiency::Float64
    darkcounts::Float64
    function Detector(deadtime::Float64,resolution::Float64,jitter::Float64,efficiency::Float64,darkcounts::Float64)
        @assert 0 ≤ efficiency ≤ 1 "Efficiency must be between 0 and 1"
        new(
            deadtime,
            resolution,
            jitter,
            efficiency,
            darkcounts
        )
    end
end

"""
    Detector(deadtime::Number,resolution::Number,jitter::Number,efficiency::Number,darkcounts::Number)

Returns a Detector object. Also accepts a dictionary with keywords matching the inputs below.

Inputs:
    deadtime   : dead time in nanoseconds
    resolution : time resolution in nanoseconds
    jitter     : timing jitter in nanoseconds
    efficiency : quantum efficiency ∈ [0,1]
    darkcounts : dark count rate in GHz
"""
function Detector(deadtime::Number,resolution::Number,jitter::Number,efficiency::Number,darkcounts::Number)
    return Detector(
        convert(Float64,deadtime),
        convert(Float64,resolution),
        convert(Float64,jitter),
        convert(Float64,efficiency),
        convert(Float64,darkcounts)
    )
end

function Detector(params::Dict)
    return Detector(
            params[:deadtime],
            params[:resolution],
            params[:jitter],
            params[:efficiency],
            params[:darkcounts]
        )
end

export Detector

"""
    sparseReadout(t::Number,source::LightSource,detect::Detector)

Returns a sparse vector containing the nonzero counts received by the detector for the given duration t and LightSource
"""
function sparseReadout(t::Number,source::LightSource,detect::Detector)
    field = eField(source)
    return sparseReadout(t,field,detect)
end

"""
    sparseReadout(t::Number,field::eField,detect::Detector)

Returns a sparse vector containing the nonzero counts received by the detector for the given duration t and EM field
"""
function sparseReadout(t::Number,field::eField,detect::Detector)
    nb = nbar(t,field.source)
    γint = γIntensity(nb,t,detect.resolution,field)
    return sparseReadout(γint,detect)
end

"""
    sparseReadout(γint:γIntensity,detect::Detector)

Returns a sparse vector containing the nonzero counts received by the detector for the given photon intensity time-series
"""
function sparseReadout(γint::γIntensity,detect::Detector)
    out = sparsevec(Integer[],Integer[],length(γint))
    ideadtime = convert(Int,ceil(detect.deadtime/detect.resolution))
    i = 1
    while i ≤ length(γint)
        ct = poissonCount(γint[i])
        if ct != 0
            out[i] = ct
            i+=ideadtime
        end
        i+=1
    end
    return out
end

export sparseReadout

################################################################################
"""
    denseReadout(t::Number,source::LightSource,detect::Detector)

Returns a dense vector containing the counts received by the detector for the given duration t and LightSource
"""
function denseReadout(t::Number,source::LightSource,detect::Detector)
    field = eField(source)
    return denseReadout(t,field,detect)
end

"""
    denseReadout(t::Number,field::eField,detect::Detector)

Returns a dense vector containing the counts received by the detector for the given duration t and EM field
"""
function denseReadout(t::Number,field::eField,detect::Detector)
    nb = nbar(t,field.source)
    γint = γIntensity(nb,t,detect.resolution,field)
    return denseReadout(γint,detect)
end

"""
    denseReadout(γint:γIntensity,detect::Detector)

Returns a dense vector containing the counts received by the detector for the given photon intensity time-series
"""
function denseReadout(γint::γIntensity,detect::Detector)
    rngvec = Poisson.(γint.γvec)
    out = zeros(Int,length(γint))
    ideadtime = convert(Int,ceil(detect.deadtime/detect.resolution))
    i = 1
    while i ≤ length(γint)
        ct = rand(rngvec[i])
        if ct != 0
            out[i] = ct
            i+=ideadtime
        end
        i+=1
    end
    return out
end

export denseReadout


text:
"""
Struct for Carr Madan Method

		bsProcess=CarrMadanMethod(σ::num1) where {num1 <: Number}

Where:\n
		A	    = 	volatility of the process.
		Npow	=	volatility of the process.
"""
mutable struct CarrMadanMethod{num <: Number, num_1 <: Integer} <: AbstractFFTMethod
    A::num
    Npow::num_1
    function CarrMadanMethod(A::num, N::num_1) where {num <: Number, num_1 <: Integer}
        @assert(A > 0.0, "A must be positive")
        @assert(N > 2, "N must be greater than 2")
        @assert(N < 24, "N will cause overflow")
        return new{num, num_1}(A, N)
    end
end
export CarrMadanMethod;

using Interpolations, FFTW;
"""
Pricing European Options through Fast Fourier Transform Method (Carr Madan)

		VectorOfPrice=CarrMadanPricer(mcProcess::FinancialMonteCarlo.BaseProcess,S0::Number,StrikeVec::Array{U},r::Number,T::Number,Npow::Integer,A::Real,d::Number=0.0) where {U <: Number}

Where:\n
		S0 = Spot price.
		StrikeVec = Vector of Strike of the Option to price.
		r= zero rate with tenor T.
		T= tenor of the options.
		Npow= Integer Parameter for the FFT. Represent the log2 of the number of nodes.
		A= Real Parameter of the FFT. Represent a maximum for the axis.

		VectorOfPrice= Price of the European Options with Strike equals to StrikeVec, tenor T and the market prices a risk free rate of r.
"""
function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, StrikeVec::Array{U, 1}, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, T::Number, method::CarrMadanMethod, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode()) where {U <: Number}
    Npow = method.Npow
    N = 2^Npow
    S0 = mcProcess.underlying.S0
    r = FinancialMonteCarlo.integral(zero_rate.r, T) / T
    d = FinancialMonteCarlo.integral(FinancialMonteCarlo.dividend(mcProcess), T) / T
    A = method.A
    #* v-> compute integral as a summation
    dx = A / N
    real_vec = 0:(N-1)
    v = collect(real_vec) * dx            # the final value A is excluded
    v[1] = 1e-312
    cf = FinancialFFT.CharactheristicFunction(mcProcess, T)
    corr = FinancialFFT.CharactheristicExponent(-1im, mcProcess, T)
    CharFunc(v) = cf(v) * exp(-v * 1im * corr)
    integrand_f(v) = exp(1im * (r - d) * v * T) * (CharFunc(v - 1im) - 1) / (1im * v - v^2)
    # Option Price
    weights_ = @. 3 + (-1)^((0:(N-1)) + 1)# Simpson weights
    @views weights_[1] = 1
    @views weights_[end] = 1
    dk = 2 * pi / A
    b = N * dk / 2
    complex_vec_z_k = @. integrand_f(v) * dx * weights_ / 3 * exp(1im * b * v)
    fft!(complex_vec_z_k)
    z_T = @. real_mod(complex_vec_z_k) / pi

    ks = -b .+ dk * real_vec
    spline_cub = CubicSplineInterpolation(ks, z_T)
    k_interp = @. log(StrikeVec / S0)
    C = @. S0 * spline_cub(k_interp) + max(S0 - StrikeVec * exp(-(r - d) * T), 0)
    return C * exp(-d * T)
end

function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, method::AbstractFFTMethod, abstractPayoffs::Array{U}, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode()) where {U <: FinancialMonteCarlo.EuropeanOption}
    TT = unique([opt.T for opt in abstractPayoffs])
    zero_typed = FinancialMonteCarlo.predict_output_type_zero(mcProcess, zero_rate, abstractPayoffs)
    prices = Array{typeof(zero_typed)}(undef, length(abstractPayoffs))

    for T in TT
        index_same_t = findall(op -> (op.T == T), abstractPayoffs)
        payoffs = abstractPayoffs[index_same_t]
        strikes = [opt.K for opt in payoffs]
        prices_call = pricer(mcProcess, strikes, zero_rate, T, method)
        prices_final = [call_to_put(prices_call[i], mcProcess.underlying, zero_rate, payoffs[i]) for i = 1:length(payoffs)]
        prices[index_same_t] = prices_final
    end

    return prices
end

function pricer(mcProcess::FinancialMonteCarlo.BaseProcess, zero_rate::FinancialMonteCarlo.AbstractZeroRateCurve, method::AbstractFFTMethod, abstractPayoff::FinancialMonteCarlo.EuropeanOption, ::FinancialMonteCarlo.BaseMode = FinancialMonteCarlo.SerialMode())
    return first(pricer(mcProcess, [abstractPayoff.K], zero_rate, abstractPayoff.T, method))
end


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.14"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	faces, N = 1:6, 10^6
	mcEstimate = counts(rand(faces,N), faces)/N
end

# ╔═╡ b93d18b8-2161-11eb-299b-0d801d262dd3
begin
	plot(faces, mcEstimate,
		line=:stem, marker=:circle,
		c=:blue, ms=10, msw=0, lw=4, label="MC estimate")
	plot!([i for i in faces], [1/6 for _ in faces],
		line=:stem, marker=:xcross, c=:red,
		ms=6, msw=0, lw=2, label="PMF",
		xlabel="Face number", ylabel="Probability", ylims=(0,0.22))
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.14"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═b93d18b8-2161-11eb-299b-0d801d262dd3
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using CSV,DataFrames
export COVID_data_tool
#COVID_data_tool.diccionario_indicadores
#=arreglo_dif=[]
ind=[]
for k in keys(COVID_data_tool.diccionario_estados)
    push!(ind,k)
end
ind
=#
#=
for k in ind
    @show k
    COVID_data_tool.datos_indicador("Porcentaje de hombres",k)
end
=#

ind = ["Coahuila de Zaragoza","Guerrero","Hidalgo","Nayarit","Querétaro","Quintana Roo"]

for k in ind
    @show k
    COVID_data_tool.datos_indicador("Porcentaje de mujeres",k)
    #push!(arreglo_dif, k)
end

#ind=setdiff(ind,arreglo_dif)


text:
adapters = open("adapters.txt") do f
    sort(parse.(Int, readlines(f)))
end

function find_contiguous_sets(numbers::Array{Int64,1}, n::Int64)::Array{Int64, 1}
    sets = []
    i = 1

    while i < length(numbers)
        if numbers[i] != n
            i += 1
            continue
        end

        j = i
        while j < length(numbers)
            if numbers[j+1] != n
                break
            end
            j += 1
        end

        if j > i
            push!(sets, length(numbers[i:j]))
        end

        i = j + 1
    end

    return sets
end

function permutations(x::Int64)::Int64
    sum = 0

    for i in (x ÷ 3):-1:0
        y = x - 3*i

        for j in (y ÷ 2):-1:0
            z = y - 2*j
            p = factorial(z+j+i) ÷ (factorial(z)*factorial(j)*factorial(i))
            sum += p
            # println("$(x): [$(z), $(j), $(i)]: $(p)")
        end
    end

    return sum
end

pushfirst!(adapters, 0)
diffs = push!(map(i -> adapters[i]-adapters[i-1], 2:length(adapters)), 3)
println(count(i -> i == 1, diffs) * count(i -> i == 3, diffs))
println(prod(map(x -> permutations(x), find_contiguous_sets(diffs, 1))))


text:
using JuMP, BARON, Base.Test

m = Model(solver=BaronSolver())
ub = [9.422, 5.9023, 267.417085245]
@variable(m, 0 ≤ x[i=1:3] ≤ ub[i])

@NLconstraints(m, begin
    250 + 30x[1] -  6x[1]^2 == x[3]
    300 + 20x[2] - 12x[2]^2 == x[3]
    150 + 0.5*(x[1]+x[2])^2 == x[3]
end)

@NLobjective(m, Min, -x[3])

solve(m)

@test isapprox(getvalue(x[1]), 6.2934300, rtol=1e-6)
@test isapprox(getvalue(x[2]), 3.8218391, rtol=1e-6)
@test isapprox(getvalue(x[3]), 201.1593341, rtol=1e-5)
@test isapprox(getobjectivevalue(m), -201.1593341, rtol=1e-6)


text:
export parallel

parallel() = schedule_on(ThreadsScheduler())


text:
"""
Like UnitRange{Int}, but can handle the `last` field, being temporarily
< first (this can happen during compacting)
"""
struct StmtRange <: Base.AbstractUnitRange{Int}
    first::Int
    last::Int
end
Base.first(r::StmtRange) = r.first
Base.last(r::StmtRange) = r.last
Base.start(r::StmtRange) = 0
Base.done(r::StmtRange, state) = r.last - r.first < state
Base.next(r::StmtRange, state) = (r.first + state, state + 1)

struct BasicBlock
    stmts::StmtRange
    preds::Vector{Int}
    succs::Vector{Int}
end
function BasicBlock(stmts::StmtRange)
    BasicBlock(stmts, Int[], Int[])
end
function BasicBlock(old_bb, stmts)
    BasicBlock(stmts, old_bb.preds, old_bb.succs)
end

struct CFG
    blocks::Vector{BasicBlock}
    index::Vector{Int}
end

function Base.show(io::IO, cfg::CFG)
    foreach(pairs(cfg.blocks)) do (idx, block)
        println("$idx\t=>\t", join(block.succs, ", "))
    end
end

struct IRCode
    stmts::Vector{Any}
    types::Vector{Any}
    argtypes::Vector{Any}
    cfg::CFG
    new_nodes::Vector{Tuple{Int, Any, Any}}
    mod::Module
end
IRCode(stmts, cfg, argtypes, mod) = IRCode(stmts, Any[], argtypes, cfg, Tuple{Int, Type, Any}[], mod)

function Base.getindex(x::IRCode, s::SSAValue)
    if s.id <= length(x.stmts)
        return x.stmts[s.id]
    else
        return x.new_nodes[s.id - length(x.stmts)][3]
    end
end

struct OldSSAValue
    id::Int
end

struct NewSSAValue
    id::Int
end

mutable struct UseRefIterator{T}
    stmt::T
end
Base.getindex(it::UseRefIterator) = it.stmt

struct UseRef
    urs::UseRefIterator
    use::Int
end

struct OOBToken
end

struct UndefToken
end

function Base.getindex(x::UseRef)
    stmt = x.urs.stmt
    if isa(stmt, Expr) && is_relevant_expr(stmt)
        x.use > length(stmt.args) && return OOBToken()
        stmt.args[x.use]
    elseif isa(stmt, GotoIfNot)
        x.use == 1 || return OOBToken()
        return stmt.cond
    elseif isa(stmt, ReturnNode) || isa(stmt, PiNode)
        isdefined(stmt, :val) || return OOBToken()
        x.use == 1 || return OOBToken()
        return stmt.val
    elseif isa(stmt, PhiNode)
        x.use > length(stmt.values) && return OOBToken()
        isassigned(stmt.values, x.use) || return UndefToken()
        return stmt.values[x.use]
    else
        return OOBToken()
    end
end

function is_relevant_expr(e::Expr)
    isexpr(e, :call) || isexpr(e, :invoke) ||
    isexpr(e, :new) || isexpr(e, :gc_preserve_begin) || isexpr(e, :gc_preserve_end) ||
    isexpr(e, :foreigncall) || isexpr(e, :isdefined) || isexpr(e, :undefcheck) ||
    isexpr(e, :throw_undef_if_not)
end

function Base.setindex!(x::UseRef, v)
    stmt = x.urs.stmt
    if isa(stmt, Expr) && is_relevant_expr(stmt)
        x.use > length(stmt.args) && throw(BoundsError())
        stmt.args[x.use] = v
    elseif isa(stmt, GotoIfNot)
        x.use == 1 || throw(BoundsError())
        x.urs.stmt = GotoIfNot{Any}(v, stmt.dest)
    elseif isa(stmt, ReturnNode)
        x.use == 1 || throw(BoundsError())
        x.urs.stmt = typeof(stmt)(v)
    elseif isa(stmt, PiNode)
        x.use == 1 || throw(BoundsError())
        x.urs.stmt = typeof(stmt)(v, stmt.typ)
    elseif isa(stmt, PhiNode)
        x.use > length(stmt.values) && throw(BoundsError())
        isassigned(stmt.values, x.use) || throw(BoundsError())
        stmt.values[x.use] = v
    else
        return OOBToken()
    end
end

function userefs(x)
    if (isa(x, Expr) && is_relevant_expr(x)) ||
        isa(x, Union{GotoIfNot, ReturnNode, PiNode, PhiNode})
        UseRefIterator(x)
    else
        ()
    end
end

Base.start(it::UseRefIterator) = 1
function Base.next(it::UseRefIterator, use)
    x = UseRef(it, use)
    v = x[]
    v === UndefToken() && return next(it, use + 1)
    x, use + 1
end
function Base.done(it::UseRefIterator, use)
    x, _ = next(it, use)
    v = x[]
    v === OOBToken() && return true
    false
end

function scan_ssa_use!(used, stmt)
    if isa(stmt, SSAValue)
        push!(used, stmt.id)
    end
    for useref in userefs(stmt)
        val = useref[]
        if isa(val, SSAValue)
            push!(used, val.id)
        end
    end
end

function ssamap(f, stmt)
    urs = userefs(stmt)
    urs === () && return stmt
    for op in urs
        val = op[]
        if isa(val, SSAValue)
            op[] = f(val)
        end
    end
    urs[]
end

function foreachssa(f, stmt)
    for op in userefs(stmt)
        val = op[]
        if isa(val, SSAValue)
            f(val)
        end
    end
end

function print_node(io::IO, idx, stmt, used, maxsize; color = true, print_typ=true)
    if idx in used
        pad = " "^(maxsize-length(string(idx)))
        print(io, "%$idx $pad= ")
    else
        print(io, " "^(maxsize+4))
    end
    if isa(stmt, PhiNode)
        args = map(1:length(stmt.edges)) do i
            e = stmt.edges[i]
            v = !isassigned(stmt.values, i) ? "#undef" :
                sprint() do io′
                    print_ssa(io′, stmt.values[i])
                end
            "$e => $v"
        end
        print(io, "φ ", '(', join(args, ", "), ')')
    elseif isa(stmt, PiNode)
        print(io, "π (")
        print_ssa(io, stmt.val)
        print(io, ", ")
        if color
            printstyled(io, stmt.typ, color=:red)
        else
            print(io, stmt.typ)
        end
        print(io, ")")
    elseif isa(stmt, ReturnNode)
        if !isdefined(stmt, :val)
            print(io, "unreachable")
        else
            print(io, "return ")
            print_ssa(io, stmt.val)
        end
    elseif isa(stmt, GotoIfNot)
        print(io, "goto ", stmt.dest, " if not ")
        print_ssa(io, stmt.cond)
    elseif isexpr(stmt, :call)
        print_ssa(io, stmt.args[1])
        print(io, "(")
        print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args[2:end]), ", "))
        print(io, ")")
        if print_typ && stmt.typ !== Any
            print(io, "::$(stmt.typ)")
        end
    elseif isexpr(stmt, :new)
        print(io, "new(")
        print(io, join(map(arg->sprint(io->print_ssa(io, arg)), stmt.args), ", "))
        print(io, ")")
    else
        print(io, stmt)
    end
end

print_ssa(io::IO, val) = isa(val, SSAValue) ? print(io, "%$(val.id)") : print(io, val)
function Base.show(io::IO, code::IRCode)
    used = Set{Int}()
    println(io, "Code")
    foreach(stmt->scan_ssa_use!(used, stmt), code.stmts)
    foreach(((_a,_b,node),)->scan_ssa_use!(used, node), code.new_nodes)
    if isempty(used)
        maxsize = 0
    else
        maxused = maximum(used)
        maxsize = length(string(maxused))
    end
    cfg = code.cfg
    max_bb_idx_size = length(string(length(cfg.blocks)))
    bb_idx = 1
    perm = sortperm(code.new_nodes, by = x->x[1])
    new_nodes_perm = Iterators.Stateful(perm)
    for (idx, stmt) in enumerate(code.stmts)
        bbrange = cfg.blocks[bb_idx].stmts
        bb_pad = max_bb_idx_size - length(string(bb_idx))
        if idx != last(bbrange)
            if idx == first(bbrange)
                print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
            else
                print(io, "│  "," "^max_bb_idx_size)
            end
        end
        print_sep = false
        if idx == last(bbrange)
            print_sep = true
        end
        floop = true
        while !isempty(new_nodes_perm) && code.new_nodes[Base.peek(new_nodes_perm)][1] == idx
            node_idx = popfirst!(new_nodes_perm)
            _, typ, node = code.new_nodes[node_idx]
            node_idx += length(code.stmts)
            if print_sep
                if floop
                    print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
                else
                    print(io, "│  "," "^max_bb_idx_size)
                end
            end
            print_sep = true
            floop = false
            print_ssa_typ = !isa(node, PiNode) && node_idx in used
            Base.with_output_color(:yellow, io) do io′
                print_node(io′, node_idx, node, used, maxsize; color = false,
                    print_typ=!print_ssa_typ || (isa(node, Expr) && typ != node.typ))
            end
            if print_ssa_typ
                printstyled(io, "::$(typ)", color=:red)
            end
            println(io)
        end
        if print_sep
            if idx == first(bbrange) && floop
                print(io, "$(bb_idx) ","─"^(1+bb_pad)," ")
            else
                print(io, idx == last(bbrange) ? string("└", "─"^(1+max_bb_idx_size), " ") :
                    string("│  ", " "^max_bb_idx_size))
            end
        end
        if idx == last(bbrange)
            bb_idx += 1
        end
        typ = code.types[idx]
        print_ssa_typ = !isa(stmt, PiNode) && idx in used
        print_node(io, idx, stmt, used, maxsize,
            print_typ=!print_ssa_typ || (isa(stmt, Expr) && typ != stmt.typ))
        if print_ssa_typ
            printstyled(io, "::$(typ)", color=:red)
        end
        println(io)
    end
end

mutable struct IncrementalCompact
    ir::IRCode
    result::Vector{Any}
    result_types::Vector{Any}
    ssa_rename::Vector{Any}
    used_ssas::Vector{Int}
    late_fixup::Vector{Int}
    keep_meta::Bool
    new_nodes_perm::Any
    idx::Int
    result_idx::Int
    function IncrementalCompact(code::IRCode; keep_meta = false)
        new_nodes_perm = Iterators.Stateful(sortperm(code.new_nodes, by=x->x[1]))
        result = Array{Any}(uninitialized, length(code.stmts) + length(code.new_nodes))
        result_types = Array{Any}(uninitialized, length(code.stmts) + length(code.new_nodes))
        ssa_rename = Any[SSAValue(i) for i = 1:(length(code.stmts) + length(code.new_nodes))]
        late_fixup = Vector{Int}()
        used_ssas = fill(0, length(code.stmts) + length(code.new_nodes))
        new(code, result, result_types, ssa_rename, used_ssas, late_fixup, keep_meta, new_nodes_perm, 1, 1)
    end
end

struct TypesView
    compact::IncrementalCompact
end
types(compact::IncrementalCompact) = TypesView(compact)

function Base.getindex(compact::IncrementalCompact, idx)
    if idx < compact.result_idx
        return compact.result[idx]
    else
        return compact.ir.stmts[idx]
    end
end

function Base.setindex!(compact::IncrementalCompact, v, idx)
    if idx < compact.result_idx
        # Kill count for current uses
        for ops in userefs(compact.result[idx])
            val = ops[]
            isa(val, SSAValue) && (compact.used_ssas[val.id] -= 1)
        end
        # Add count for new use
        isa(v, SSAValue) && (compact.used_ssas[v.id] += 1)
        return compact.result[idx] = v
    else
        return compact.ir.stmts[idx] = v
    end
end

function Base.getindex(view::TypesView, idx)
    if idx < view.compact.result_idx
        return view.compact.result_types[idx]
    else
        return view.compact.ir.types[idx]
    end
end

function value_typ(ir::IRCode, value)
    isa(value, SSAValue) && return ir.types[value.id]
    isa(value, GlobalRef) && return typeof(getfield(value.mod, value.name))
    isa(value, Argument) && return ir.argtypes[value.n]
    return typeof(value)
end

function value_typ(ir::IncrementalCompact, value)
    isa(value, SSAValue) && return types(ir)[value.id]
    isa(value, GlobalRef) && return typeof(getfield(value.mod, value.name))
    isa(value, Argument) && return ir.ir.argtypes[value.n]
    return typeof(value)
end

Base.start(compact::IncrementalCompact) = (1,1,1)
function Base.done(compact::IncrementalCompact, (idx, _a, _b)::Tuple{Int, Int, Int})
    return idx > length(compact.ir.stmts) && isempty(compact.new_nodes_perm)
end

function process_node!(result, result_idx, ssa_rename, late_fixup, used_ssas, stmt, idx, processed_idx, keep_meta)
    ssa_rename[idx] = SSAValue(result_idx)
    if stmt === nothing
        ssa_rename[idx] = stmt
    elseif !keep_meta && (isexpr(stmt, :meta) || isa(stmt, LineNumberNode))
        # eliminate this node
    elseif isa(stmt, GotoNode)
        result[result_idx] = stmt
        result_idx += 1
    elseif isexpr(stmt, :call) || isexpr(stmt, :invoke) || isa(stmt, ReturnNode) || isexpr(stmt, :gc_preserve_begin) ||
           isexpr(stmt, :gc_preserve_end) || isexpr(stmt, :foreigncall)
        result[result_idx] = renumber_ssa!(stmt, ssa_rename, true, used_ssas)
        result_idx += 1
    elseif isa(stmt, PhiNode)
        values = Vector{Any}(uninitialized, length(stmt.values))
        for i = 1:length(stmt.values)
            isassigned(stmt.values, i) || continue
            val = stmt.values[i]
            if isa(val, SSAValue)
                if val.id > processed_idx
                    push!(late_fixup, result_idx)
                    val = OldSSAValue(val.id)
                else
                    val = renumber_ssa!(val, ssa_rename, true, used_ssas)
                end
            end
            values[i] = val
        end
        result[result_idx] = PhiNode(stmt.edges, values)
        result_idx += 1
    elseif isa(stmt, SSAValue) || (!isa(stmt, Expr) && !isa(stmt, PhiNode) && !isa(stmt, PiNode) && !isa(stmt, GotoIfNot))
        # Constant or identity assign, replace uses of this
        # ssa value with its result
        stmt = isa(stmt, SSAValue) ? ssa_rename[stmt.id] : stmt
        ssa_rename[idx] = stmt
    else
        result[result_idx] = renumber_ssa!(stmt, ssa_rename, true, used_ssas)
        result_idx += 1
    end
    return result_idx
end
function process_node!(compact::IncrementalCompact, result_idx, stmt, idx, processed_idx)
    process_node!(compact.result, result_idx, compact.ssa_rename,
        compact.late_fixup, compact.used_ssas, stmt, idx, processed_idx, compact.keep_meta)
end

function Base.next(compact::IncrementalCompact, (idx, active_bb, old_result_idx)::Tuple{Int, Int, Int})
    if length(compact.result) < old_result_idx
        resize!(compact.result, old_result_idx)
        resize!(compact.result_types, old_result_idx)
    end
    bb = compact.ir.cfg.blocks[active_bb]
    if !isempty(compact.new_nodes_perm) && compact.ir.new_nodes[Base.peek(compact.new_nodes_perm)][1] == idx
        new_idx = popfirst!(compact.new_nodes_perm)
        _, typ, new_node = compact.ir.new_nodes[new_idx]
        new_idx += length(compact.ir.stmts)
        compact.result_types[old_result_idx] = typ
        result_idx = process_node!(compact, old_result_idx, new_node, new_idx, idx)
        (old_result_idx == result_idx) && return next(compact, (idx, result_idx))
        compact.result_idx = result_idx
        return (old_result_idx, compact.result[old_result_idx]), (compact.idx, active_bb, compact.result_idx)
    end
    # This will get overwritten in future iterations if
    # result_idx is not, incremented, but that's ok and expected
    compact.result_types[old_result_idx] = compact.ir.types[idx]
    result_idx = process_node!(compact, old_result_idx, compact.ir.stmts[idx], idx, idx)
    if idx == last(bb.stmts)
        # If this was the last statement in the BB and we decided to skip it, insert a
        # dummy `nothing` node, to prevent changing the structure of the CFG
        if result_idx == first(bb.stmts)
            compact.result[old_result_idx] = nothing
            result_idx = old_result_idx + 1
        end
        compact.ir.cfg.blocks[active_bb] = BasicBlock(bb, StmtRange(first(bb.stmts), result_idx-1))
        active_bb += 1
        if active_bb <= length(compact.ir.cfg.blocks)
            new_bb = compact.ir.cfg.blocks[active_bb]
            compact.ir.cfg.blocks[active_bb] = BasicBlock(new_bb,
                StmtRange(result_idx, last(new_bb.stmts)))
        end
    end
    (old_result_idx == result_idx) && return next(compact, (idx + 1, active_bb, result_idx))
    compact.idx = idx + 1
    compact.result_idx = result_idx
    if !isassigned(compact.result, old_result_idx)
        ccall(:jl_, Cvoid, (Any,), compact.ir.stmts)
        ccall(:jl_, Cvoid, (Any,), (compact.ir.stmts[idx], old_result_idx, result_idx, idx))
        ccall(:jl_, Cvoid, (Any,), compact.result)
        @assert false
    end
    return (old_result_idx, compact.result[old_result_idx]), (compact.idx, active_bb, compact.result_idx)
end

function maybe_erase_unused!(extra_worklist, compact, idx)
   if stmt_effect_free(compact.result[idx], compact.ir, compact.ir.mod)
        for ops in userefs(compact.result[idx])
            val = ops[]
            if isa(val, SSAValue)
                if compact.used_ssas[val.id] == 1
                    if val.id < idx
                        push!(extra_worklist, val.id)
                    end
                end
                compact.used_ssas[val.id] -= 1
            end
        end
        compact.result[idx] = nothing
    end
end

function finish(compact::IncrementalCompact)
    for idx in compact.late_fixup
        stmt = compact.result[idx]
        @assert isa(stmt, PhiNode)
        values = Vector{Any}(uninitialized, length(stmt.values))
        for i = 1:length(stmt.values)
            isassigned(stmt.values, i) || continue
            val = stmt.values[i]
            if isa(val, OldSSAValue)
                val = compact.ssa_rename[val.id]
                if isa(val, SSAValue)
                    compact.used_ssas[val.id] += 1
                end
            end
            values[i] = val
        end
        compact.result[idx] = PhiNode(stmt.edges, values)
    end
    # Record this somewhere?
    result_idx = compact.result_idx
    resize!(compact.result, result_idx-1)
    resize!(compact.result_types, result_idx-1)
    bb = compact.ir.cfg.blocks[end]
    compact.ir.cfg.blocks[end] = BasicBlock(bb,
                StmtRange(first(bb.stmts), result_idx-1))
    # Perform simple DCE for unused values
    extra_worklist = Int[]
    for (idx, nused) in enumerate(compact.used_ssas)
        idx >= result_idx && break
        nused == 0 || continue
        maybe_erase_unused!(extra_worklist, compact, idx)
    end
    while !isempty(extra_worklist)
        maybe_erase_unused!(extra_worklist, compact, pop!(extra_worklist))
    end
    cfg = CFG(compact.ir.cfg.blocks, Int[first(bb.stmts) for bb in compact.ir.cfg.blocks[2:end]])
    IRCode(compact.result, compact.result_types, compact.ir.argtypes, cfg, Any[], compact.ir.mod)
end

function compact!(code::IRCode)
    compact = IncrementalCompact(code)
    # Just run through the iterator without any processing
    foreach(_->nothing, compact)
    return finish(compact)
end

struct BBIdxStmt
    ir::IRCode
end

bbidxstmt(ir) = BBIdxStmt(ir)

Base.start(x::BBIdxStmt) = (1,1)
Base.done(x::BBIdxStmt, (idx, bb)) = idx > length(x.ir.stmts)
function Base.next(x::BBIdxStmt, (idx, bb))
    active_bb = x.ir.cfg.blocks[bb]
    next_bb = bb
    if idx == last(active_bb.stmts)
        next_bb += 1
    end
    return (bb, idx, x.ir.stmts[idx]), (idx + 1, next_bb)
end


text:
@with_kw struct PlaneWaveSource{GeometryType,E <: AbstractExcitation} <: AbstractSource
    geometry::GeometryType
    excitation::E
    polarization::Vec
    propagation_direction::Vec
end

excitation(planewave_source::PlaneWaveSource) = planewave_source.excitation
propagation_direction(planewave_source::PlaneWaveSource) = planewave_source.propagation_direction

k⃗₀(planewave_source::PlaneWaveSource) =  k₀(excitation(planewave_source))  .* propagation_direction(planewave_source)


text:
using Knet
include("../helper.jl")
include("accuracy.jl")
include("loss.jl")
include("mpii.batcher.jl")
include("../models.jl")
include("../training.jl")
include("../plots.jl")
include("../modular.resnet.jl")
include("batch_loading_script.jl")

# example usage
# continue_training("29-12-20", "29-12-20-cont", [(0.002, 10), (0.001, 12)])

function init_training(; save_tag, learning_rate_by_epochs, should_use_pmap = true, check_loss = false)
    dtrn, dval, dtst, trn_data_items, val_data_items, tst_data_items = load_and_save_batches(save_tag; should_use_pmap = should_use_pmap)

    init_training_with_resnet152(
        dtrn,
        dval,
        val_data_items;
        save_tag = save_tag,
        learning_rate_by_epochs = learning_rate_by_epochs,
        check_loss = check_loss
    )
end

function init_training_with_resnet152(
    dtrn,
    dval,
    val_data_items;
    save_tag,
    learning_rate_by_epochs,
    check_loss = true
)
    CUDA.reclaim()
    CUDA.memory_status()
    GC.gc(true)
    println("cleared gpu memory")
    # Uses ResNet 152
    deeper_cut_with_loc_ref =
        generate_deeper_cut(; should_use_resnet50 = false, is_loc_ref_enabled = true)
    println("prepared deepercut")

    CUDA.reclaim()
    CUDA.memory_status()
    GC.gc(true)
    println("cleared gpu memory")

    # [(0.005, 1), (0.02, 17), (0.002, 12), (0.001, 12)]
    println("started training")
    deeper_cut_with_loc_ref_results = trainresults(
        save_tag,
        deeper_cut_with_loc_ref,
        dtrn,
        dval,
        learning_rate_by_epochs;
        accuracy_func = nothing,
        error_func = nothing,
        secondary_accuracy_func = modelized_PCKh_sigm,
        data_items_trn = nothing,
        data_items_val = val_data_items,
        check_loss = check_loss
    )
end

# learning_rate_by_epochs: [(0.005, 1), (0.02, 17), (0.002, 12), (0.001, 12)]
function continue_training(save_tag, new_save_tag, learning_rate_by_epochs)
    println("start loading dtrn")
    dtrn = Knet.load("$(pre_full_path)$(save_tag)-dtrn.jld2", "$(save_tag)-dtrn")

    println("start loading dval")
    dval = Knet.load("$(pre_full_path)$(save_tag)-dval.jld2", "$(save_tag)-dval")

    println("start loading val data items")
    val_data_items =
        Knet.load("$(pre_full_path)$(save_tag)-val_data_items.jld2", "$(save_tag)-val_data_items")

    println("start loading model")
    model = Knet.load("$(pre_full_path)$(save_tag)-training_model.jld2", "$(save_tag)-training_model")

    CUDA.reclaim()
    CUDA.memory_status()
    GC.gc(true)
    println("cleared gpu memory")

    trainresults(
        new_save_tag,
        model,
        dtrn,
        dval,
        learning_rate_by_epochs;
        accuracy_func = nothing,
        error_func = nothing,
        secondary_accuracy_func = modelized_PCKh_sigm,
        data_items_trn = nothing,
        data_items_val = val_data_items,
    )

end


text:
"""
	characteristicMatrix( FV::Cells )::ChainOp

Binary matrix representing by rows the `p`-cells of a cellular complex.
The input parameter must be of `Cells` type. Return a sparse binary matrix,
providing the basis of a ``Chain`` space of given dimension. Notice that the
number of columns is equal to the number of vertices (0-cells).

# Example

```julia
V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);

julia> full(characteristicMatrix(FV))
6×8 Array{Int8,2}:
	1  1  1  1  0  0  0  0
	0  0  0  0  1  1  1  1
	1  1  0  0  1  1  0  0
	0  0  1  1  0  0  1  1
	1  0  1  0  1  0  1  0
	0  1  0  1  0  1  0  1

julia> full(characteristicMatrix(CV))
1×8 Array{Int8,2}:
	1  1  1  1  1  1  1  1

julia> full(characteristicMatrix(EV))
12×8 Array{Int8,2}:
	1  1  0  0  0  0  0  0
	0  0  1  1  0  0  0  0
	0  0  0  0  1  1  0  0
	0  0  0  0  0  0  1  1
	1  0  1  0  0  0  0  0
	0  1  0  1  0  0  0  0
	0  0  0  0  1  0  1  0
	0  0  0  0  0  1  0  1
	1  0  0  0  1  0  0  0
	0  1  0  0  0  1  0  0
	0  0  1  0  0  0  1  0
	0  0  0  1  0  0  0  1
```
"""
function characteristicMatrix( FV::Cells )::ChainOp
	I,J,V = Int64[],Int64[],Int8[]
	for f=1:length(FV)
		for k in FV[f]
		push!(I,f)
		push!(J,k)
		push!(V,1)
		end
	end
	M_2 = sparse(I,J,V)
	return M_2
end

"""
	boundary_1( EV::Cells )::ChainOp

Computation of sparse signed boundary operator ``C_1 -> C_0``.

# Example
```julia
julia> V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);

julia> EV
12-element Array{Array{Int64,1},1}:
	[1, 2]
	[3, 4]
	...
	[2, 6]
	[3, 7]
	[4, 8]

julia> boundary_1( EV::Cells )
8×12 SparseMatrixCSC{Int8,Int64} with 24 stored entries:
	[1 ,  1]  =  -1
	[2 ,  1]  =  1
	[3 ,  2]  =  -1
	...       ...
	[7 , 11]  =  1
	[4 , 12]  =  -1
	[8 , 12]  =  1

julia> full(boundary_1(EV::Cells))
8×12 Array{Int8,2}:
	-1   0   0   0  -1   0   0   0  -1   0   0   0
	1   0   0   0   0  -1   0   0   0  -1   0   0
	0  -1   0   0   1   0   0   0   0   0  -1   0
	0   1   0   0   0   1   0   0   0   0   0  -1
	0   0  -1   0   0   0  -1   0   1   0   0   0
	0   0   1   0   0   0   0  -1   0   1   0   0
	0   0   0  -1   0   0   1   0   0   0   1   0
	0   0   0   1   0   0   0   1   0   0   0   1
```
"""
function boundary_1( EV::Cells )::ChainOp
	sp_boundary_1 = characteristicMatrix(EV)'
	for e = 1:length(EV)
		sp_boundary_1[EV[e][1],e] = -1
	end
	return sp_boundary_1
end

"""
	coboundary_0(EV::LinearAlgebraicRepresentation.Cells)

Return the `coboundary_0` signed operator `C_0` -> `C_1`.
"""
coboundary_0(EV::Cells) = boundary_1(EV::Cells)'

"""
	u_coboundary_1( FV::Cells, EV::::Cells)::ChainOp

Compute the sparse *unsigned* coboundary_1 operator ``C_1 -> C_2``.
Notice that the output matrix is `m x n`, where `m` is the number of faces, and `n`
is the number of edges.

# Example

```julia
julia> V,(VV,EV,FV,CV) = LinearAlgebraicRepresentation.cuboid([1.,1.,1.], true);

julia> u_coboundary_1(FV,EV)
6×12 SparseMatrixCSC{Int8,Int64} with 24 stored entries:
	[1 ,  1]  =  1
	[3 ,  1]  =  1
	[1 ,  2]  =  1
	[4 ,  2]  =  1
	...		...
	[4 , 11]  =  1
	[5 , 11]  =  1
	[4 , 12]  =  1
	[6 , 12]  =  1

julia> full(u_coboundary_1(FV,EV))
6×12 Array{Int8,2}:
	1  1  0  0  1  1  0  0  0  0  0  0
	0  0  1  1  0  0  1  1  0  0  0  0
	1  0  1  0  0  0  0  0  1  1  0  0
	0  1  0  1  0  0  0  0  0  0  1  1
	0  0  0  0  1  0  1  0  1  0  1  0
	0  0  0  0  0  1  0  1  0  1  0  1

julia> unsigned_boundary_2 = u_coboundary_1(FV,EV)';
```
"""
function u_coboundary_1( FV::LinearAlgebraicRepresentation.Cells, EV::LinearAlgebraicRepresentation.Cells)::LinearAlgebraicRepresentation.ChainOp
	cscFV = LinearAlgebraicRepresentation.characteristicMatrix(FV)
	cscEV = LinearAlgebraicRepresentation.characteristicMatrix(EV)
	temp = cscFV * cscEV'
	I,J,V = Int64[],Int64[],Int8[]
	for j=1:size(temp,2)
		for i=1:size(temp,1)
		if temp[i,j] == 2
			push!(I,i)
			push!(J,j)
			push!(V,1)
		end
		end
	end
	sp_u_coboundary_1 = sparse(I,J,V)
	return sp_u_coboundary_1
end

"""
	u_boundary_2(FV::LinearAlgebraicRepresentation.Cells, EV::LinearAlgebraicRepresentation.Cells)::LinearAlgebraicRepresentation.ChainOp

Return the unsigned `boundary_2` operator `C_2` -> `C_1`.
"""
u_boundary_2(EV, FV) = (LinearAlgebraicRepresentation.u_coboundary_1(FV, EV))'

"""
Local utility function. Storage of information need to build face cycles.
"""
function columninfo(infos,EV,next,col)
	infos[1,col] = 1
	infos[2,col] = next
	infos[3,col] = EV[next][1]
	infos[4,col] = EV[next][2]
	vpivot = infos[4,col]
end

"""
	coboundary_1( FV::Cells, EV::Cells)::ChainOp

Compute the sparse *signed* coboundary_1 operator ``C_1 -> C_2``.
The sparse matrix generated by `coboundary_1` contains by row a representation
of faces as oriented cycles of edges. The orientation of cycles is arbitrary.
```
julia> coboundary_1( FV,EV )
6×12 SparseMatrixCSC{Int8,Int64} with 24 stored entries:
	[1 ,  1]  =  -1
	[3 ,  1]  =  -1
	[1 ,  2]  =  1
	[4 ,  2]  =  -1
	...		  ...
	[4 , 11]  =  1
	[5 , 11]  =  -1
	[4 , 12]  =  -1
	[6 , 12]  =  -1

julia> full(coboundary_1( FV,EV ))
6×12 Array{Int8,2}:
	-1   1   0  0   1  -1  0   0  0   0   0   0
	0   0  -1  1   0   0  1  -1  0   0   0   0
	-1   0   1  0   0   0  0   0  1  -1   0   0
	0  -1   0  1   0   0  0   0  0   0   1  -1
	0   0   0  0  -1   0  1   0  1   0  -1   0
	0   0   0  0   0  -1  0   1  0   1   0  -1

julia> boundary_2(FV,EV) = coboundary_1(FV,EV)'
12×6 Array{Int8,2}:
	-1   0  -1   0   0   0
	1   0   0  -1   0   0
	0  -1   1   0   0   0
	...			...
	0   0  -1   0   0   1
	0   0   0   1  -1   0
	0   0   0  -1   0  -1
```
"""
function coboundary_1( FV::Cells, EV::Cells)::ChainOp
	sp_u_coboundary_1 = u_coboundary_1(FV,EV)
	larEV = characteristicMatrix(EV)
	# unsigned incidence relation
	FE = [findn(sp_u_coboundary_1[f,:]) for f=1:size(sp_u_coboundary_1,1) ]
	I,J,V = Int64[],Int64[],Int8[]
	vedges = [findn(larEV[:,v]) for v=1:size(larEV,2)]

	# Loop on faces
	for f=1:length(FE)
		fedges = Set(FE[f])
		next = pop!(fedges)
		col = 1
		infos = zeros(Int64,(4,length(FE[f])))
		vpivot = infos[4,col]
		vpivot = columninfo(infos,EV,next,col)
		while fedges != Set()
			nextedge = intersect(fedges, Set(vedges[vpivot]))
			fedges = setdiff(fedges,nextedge)
			next = pop!(nextedge)
			col += 1
			vpivot = columninfo(infos,EV,next,col)
			if vpivot == infos[4,col-1]
				infos[3,col],infos[4,col] = infos[4,col],infos[3,col]
				infos[1,col] = -1
				vpivot = infos[4,col]
			end
		end
		for j=1:size(infos,2)
			push!(I, f)
			push!(J, infos[2,j])
			push!(V, infos[1,j])
		end
	end

	sp_coboundary_1 = sparse(I,J,V)
	return sp_coboundary_1
end

"""
	chaincomplex( W::Points, EW::Cells )::Tuple{Array{Cells,1},Array{ChainOp,1}}

Chain 2-complex construction from basis of 1-cells.

From the minimal input, construct the whole
two-dimensional chain complex, i.e. the bases for linear spaces C_1 and
C_2 of 1-chains and  2-chains, and the signed coboundary operators from
C_0 to C_1 and from C_1 to C_2.

# Example
```julia
julia> W =
	[0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  2.0  2.0  2.0  2.0  3.0  3.0  3.0  3.0
	0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0]
# output
	2×16 Array{Float64,2}: ...

julia> EW =
[[1, 2],[2, 3],[3, 4],[5, 6],[6, 7],[7, 8],[9, 10],[10, 11],[11, 12],[13, 14],
	[14, 15],[15, 16],[1, 5],[2, 6],[3, 7],[4, 8],[5, 9],[6, 10],[7, 11],[8, 12],
	[9, 13],[10, 14],[11, 15],[12, 16]]
# output
24-element Array{Array{Int64,1},1}: ...

julia> V,bases,coboundaries = chaincomplex(W,EW)

julia> bases[1]	# edges
24-element Array{Array{Int64,1},1}: ...

julia> bases[2] # faces -- previously unknown !!
9-element Array{Array{Int64,1},1}: ...

julia> coboundaries[1] # coboundary_1
24×16 SparseMatrixCSC{Int8,Int64} with 48 stored entries: ...

julia> full(coboundaries[2]) # coboundary_1: faces as oriented 1-cycles of edges
9×24 Array{Int8,2}:
	-1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0  0
	0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0
	0  0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0
	0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0
	0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0
	0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0
	0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0
	0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0  0
	0  0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1
```
"""
function chaincomplex( W, EW )
	V = W'
	EV = LinearAlgebraicRepresentation.boundary_1(EW)'
	V,cscEV,cscFE = LinearAlgebraicRepresentation.planar_arrangement(V,EV)
	ne,nv = size(cscEV)
	nf = size(cscFE,1)
	EV = [findn(cscEV[e,:]) for e=1:ne]
	FV = [collect(Set(vcat([EV[e] for e in findn(cscFE[f,:])]...)))  for f=1:nf]
	function ord(cells)
		return [sort(cell) for cell in cells]
	end
	temp = copy(cscEV')
	for k=1:size(temp,2)
		h = findn(temp[:,k])[1]
		temp[h,k] = -1
	end
	cscEV = temp'
	bases, coboundaries = (ord(EV),ord(FV)), (cscEV,cscFE)
	return V',bases,coboundaries
end

"""
	chaincomplex( W::Points, FW::Cells, EW::Cells )
		::Tuple{ Array{Cells,1}, Array{ChainOp,1} }

Chain 3-complex construction from bases of 2- and 1-cells.

From the minimal input, construct the whole
two-dimensional chain complex, i.e. the bases for linear spaces C_1 and
C_2 of 1-chains and  2-chains, and the signed coboundary operators from
C_0 to C_1  and from C_1 to C_2.

# Example
```julia
julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )

julia> cube_2 = LinearAlgebraicRepresentation.Struct([LinearAlgebraicRepresentation.t(0,0,0.5), LinearAlgebraicRepresentation.r(0,0,pi/3), cube_1])

julia> V,FV,EV = LinearAlgebraicRepresentation.struct2lar(LinearAlgebraicRepresentation.Struct([ cube_1, cube_2 ]))

julia> V,bases,coboundaries = LinearAlgebraicRepresentation.chaincomplex(V,FV,EV)

julia> (EV, FV, CV), (cscEV, cscFE, cscCF) = bases,coboundaries

julia> FV # bases[2]
18-element Array{Array{Int64,1},1}:
	[1, 3, 4, 6]
	[2, 3, 5, 6]
	[7, 8, 9, 10]
	[1, 2, 3, 7, 8]
	[4, 6, 9, 10, 11, 12]
	[5, 6, 11, 12]
	[1, 4, 7, 9]
	[2, 5, 11, 13]
	[2, 8, 10, 11, 13]
	[2, 3, 14, 15, 16]
	[11, 12, 13, 17]
	[11, 12, 13, 18, 19, 20]
	[2, 3, 13, 17]
	[2, 13, 14, 18]
	[15, 16, 19, 20]
	[3, 6, 12, 15, 19]
	[3, 6, 12, 17]
	[14, 16, 18, 20]

julia> CV # bases[3]
3-element Array{Array{Int64,1},1}:
	[2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20]
	[2, 3, 5, 6, 11, 12, 13, 17]
	[1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]

julia> cscEV # coboundaries[1]
34×20 SparseMatrixCSC{Int8,Int64} with 68 stored entries: ...

julia> cscFE # coboundaries[2]
18×34 SparseMatrixCSC{Int8,Int64} with 80 stored entries: ...

julia> cscCF # coboundaries[3]
4×18 SparseMatrixCSC{Int8,Int64} with 36 stored entries: ...
```
"""
function chaincomplex(W,FW,EW)
	V = W'
	EV = LinearAlgebraicRepresentation.build_copEV(EW)
	FE = LinearAlgebraicRepresentation.coboundary_1(FW,EW)
	V,cscEV,cscFE,cscCF = LinearAlgebraicRepresentation.spatial_arrangement(V,EV,FE)
	ne,nv = size(cscEV)
	nf = size(cscFE,1)
	nc = size(cscCF,1)
	EV = [findn(cscEV[e,:]) for e=1:ne]
	FV = [collect(Set(vcat([EV[e] for e in findn(cscFE[f,:])]...)))  for f=1:nf]
	CV = [collect(Set(vcat([FV[f] for f in findn(cscCF[c,:])]...)))  for c=2:nc]
	function ord(cells)
		return [sort(cell) for cell in cells]
	end
	temp = copy(cscEV')
	for k=1:size(temp,2)
		h = findn(temp[:,k])[1]
		temp[h,k] = -1
	end
	cscEV = temp'
	bases, coboundaries = (ord(EV),ord(FV),ord(CV)), (cscEV,cscFE,cscCF)
	return V',bases,coboundaries
end

# Collect LAR models in a single LAR model
function collection2model(collection)
	W,FW,EW = collection[1]
	shiftV = size(W,2)
	for k=2:length(collection)
		V,FV,EV = collection[k]
		W = [W V]
		FW = [FW; FV + shiftV]
		EW = [EW; EV + shiftV]
		shiftV = size(W,2)
	end
	return W,FW,EW
end

# 	"""
# 		facetriangulation(V::Points, FV::Cells, EV::Cells, cscFE::ChainOp, cscCF::ChainOp)

# 	Triangulation of a single facet of a 3-complex.

# 	# Example
# 	```julia
# 	julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
# 	[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
# 	[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )

# 	julia> cube_2 = LinearAlgebraicRepresentation.Struct([LinearAlgebraicRepresentation.t(0,0,0.5), LinearAlgebraicRepresentation.r(0,0,pi/3), cube_1])

# 	julia> W,FW,EW = LinearAlgebraicRepresentation.struct2lar(LinearAlgebraicRepresentation.Struct([ cube_1, cube_2 ]))

# 	julia> V,(EV,FV,EV),(cscEV,cscFE,cscCF) = LinearAlgebraicRepresentation.chaincomplex(W,FW,EW)
# 	```
# 	"""
#    function facetriangulation(V,FV,EV,cscFE,cscCF)
#       function facetrias(f)
#          vs = [V[:,v] for v in FV[f]]
#          vs_indices = [v for v in FV[f]]
#          vdict = Dict([(i,index) for (i,index) in enumerate(vs_indices)])
#          dictv = Dict([(index,i) for (i,index) in enumerate(vs_indices)])
#          es = findn(cscFE[f,:])

#          vts = [v-vs[1] for v in vs]

#          v1 = vts[2]
#          v2 = vts[3]
#          v3 = cross(v1,v2)
#          err, i = 1e-8, 1
#          while norm(v3) < err
#             v2 = vts[3+i]
#             i += 1
#             v3 = cross(v1,v2)
#          end

#          M = [v1 v2 v3]

#          vs_2D = hcat([(inv(M)*v)[1:2] for v in vts]...)'
#          pointdict = Dict([(vs_2D[k,:],k) for k=1:size(vs_2D,1)])
#          edges = hcat([[dictv[v] for v in EV[e]]  for e in es]...)'

#          trias = Triangle.constrained_triangulation_vertices(
#             vs_2D, collect(1:length(vs)), edges)

#          triangles = [[pointdict[t[1,:]],pointdict[t[2,:]],pointdict[t[3,:]]]
#             for t in trias]
#          mktriangles = [[vdict[t[1]],vdict[t[2]],vdict[t[3]]] for t in triangles]
#          return mktriangles
#       end
#       return facetrias
#    end

#    # Triangulation of the 2-skeleton
# 	"""

# 	"""
#    function triangulate(cf,V,FV,EV,cscFE,cscCF)
#       mktriangles = facetriangulation(V,FV,EV,cscFE,cscCF)
#       TV = Array{Int64,1}[]
#       for (f,sign) in zip(cf[1],cf[2])
#          triangles = mktriangles(f)
#          if sign == 1
#             append!(TV,triangles )
#          elseif sign == -1
#             append!(TV,[[t[2],t[1],t[3]] for t in triangles] )
#          end
#       end
#       return TV
#    end

#    # Map 3-cells to local bases
# 	"""

# 	"""
#    function map_3cells_to_localbases(V,CV,FV,EV,cscCF,cscFE)
#       local3cells = []
#       for c=1:length(CV)
#          cf = findnz(cscCF[c+1,:])
#          tv = triangulate(cf,V,FV,EV,cscFE,cscCF)
#          vs = sort(collect(Set(hcat(tv...))))
#          vsdict = Dict([(v,k) for (k,v) in enumerate(vs)])
#          tvs = [[vsdict[t[1]],vsdict[t[2]],vsdict[t[3]]] for t in tv]
#          v = hcat([V[:,w] for w in vs]...)
#          cell = [v,tvs]
#          append!(local3cells,[cell])
#       end
#       return local3cells
#    end


text:
using PGFPlotsX, LaTeXStrings

include("plot_params.jl")

function add_identif_point!(ps, pb, stats, inds, vals; color="black")
    minval = pb.objopt

    k1, k2 = collect(keys(stats.hist[:iter_support_x]))
    logstep = k2-k1

    reg = typeof(pb.pb).parameters[1]
    opt_support = get_support(reg, pb.pb.x0)
    indices_identif = SortedSet(keys(filter(kv->collect(kv.second) == collect(opt_support), stats.hist[:iter_support_x])))

    if length(indices_identif) > 0
        itmax = maximum(keys(stats.hist[:iter_support_x]))
        itidentif = itmax + logstep
        while itidentif-logstep in indices_identif
            itidentif -= logstep
        end

        if itidentif <= itmax
            res = findall(x->x==itidentif, inds)# find it st inds[it] = itidentif, res = vals[it]
            indidentif = res[1]

            push!(ps, PlotInc(
                PGFPlotsX.Options(
                    "forget plot" => nothing,
                    "draw" => "none",
                    "mark" => "oplus",
                    "color" => color
                ),
                Coordinates([ (inds[indidentif], vals[indidentif]) ]),
            ))
        end
    end

    return
end

function build_save_suboptimality(pb, algo_to_stats, FIGS_FOLDER)
    minval = pb.objopt

    ps = []
    algo_id = 1
    for (algo, stats) in algo_to_stats
        inds = collect(keys(stats.hist[:iter_f]))
        vals = collect(abs.(values(stats.hist[:iter_f]) .- minval))
        coords = [ (inds[i], vals[i]) for i in 1:length(vals)]

        push!(ps, PlotInc(
            PGFPlotsX.Options(
                get_suboptimality_algoparams(algo.name)...
            ),
            Coordinates(coords)
        ))
        push!(ps, LegendEntry(get_legend(algo.name)))

        ## Add identification point
        add_identif_point!(ps, pb, stats, inds, vals, color=get_alg_color(algo.name))

        algo_id+=1
    end

    fig = @pgf Axis(
        {
            # height = "12cm",
            # width = "12cm",
            xlabel = "number of proximal gradient steps",
            legend_cell_align = "left",
            xmin = 0,
            ylabel = L"F(x_k)-F^\star",
            ymode  = "log",
            legend_style = "{font=\\footnotesize}",
            # xmode  = "log",
            # title = "Problem $(pb.name) -- Suboptimality",
            legend_pos = "north east",
        },
        ps...
    )
    pgfsave(joinpath(FIGS_FOLDER, "$(pb.name)_subopt.pdf"), fig)
    pgfsave(joinpath(FIGS_FOLDER, "$(pb.name)_subopt.tikz"), fig)
end

###############################################################################
#### Iterates figure functions

function add_contour!(ps, f, x, y)
    push!(ps, PlotInc(
            PGFPlotsX.Options(
                "forget plot" => nothing,
                "no marks" => nothing,
                "ultra thin" => nothing
            ),
            Table(contours(x, y, f.(x, y'), 10)),
    ))
    return
end

function add_traj!(ps, coords, algoname; params=Dict())
    if length(coords)>0
        push!(ps, PlotInc(
            PGFPlotsX.Options(params...),
            Coordinates(coords),
        ))
        push!(ps, LegendEntry(get_legend(algoname)))
    end
    return
end

function add_point!(ps, xopt)
    coords = [ (xopt[1], xopt[2]) ]

    push!(ps, PlotInc(
        PGFPlotsX.Options(
            "forget plot" => nothing,
            "only marks" => nothing,
            "mark" => "star",
            "thick" => nothing,
            "color" => "black"
            ),
        Coordinates(coords),
    ))

    return
end

function pgf_build_iteratesfig(ps, xmin, xmax, ymin, ymax)
    return @pgf Axis(
        {
            contour_prepared,
            # height = "12cm",
            # width = "12cm",
            # xlabel = "x",
            # ylabel = "y",
            xmin = xmin,
            xmax = xmax,
            ymin = ymin,
            ymax = ymax,
            legend_pos = "north east",
            legend_cell_align = "left",
            legend_style = "font=\\footnotesize",
            # title = "Problem $(pb.name) -- Iterates postion",
        },
        ps...,
    )
end

function add_manifold!(ps, coords)
    push!(ps, PlotInc(
        PGFPlotsX.Options(
            "forget plot" => nothing,
            "no marks" => nothing,
            "smooth" => nothing,
            "thick" => nothing,
            "solid" => nothing,
            "black!50!white" => nothing,
            # "mark size" => "1pt"
        ),
        Coordinates(coords),
    ))
    return
end


text:
__precompile__()

module Processing
	include("./Screen.jl")
	include("./Led.jl")
	include("./Types.jl")

	using BenchmarkTools
	using .Screen
	using .Led
	using .Types

	struct ProcessingConfig
		screenConfig::ScreenConfig
		ledConfig::LedConfig
	end

	function init_processing(monitorNr::Int, checkHeight::Int, sectorCount::Vector{Int}, serial::String)
		scrConfig = config_screen(monitorNr, checkHeight, sectorCount)
		ledConfig = init_serial(serial, sectorCount)
		clear_frame(ledConfig)

		return ProcessingConfig(scrConfig, ledConfig)
	end

	function update_frame(config::ProcessingConfig, prevFrame = Vector{RGB}())::Vector{RGB}
		frameData = process_screen_gpu(config.screenConfig)
		send_frame(config.ledConfig, frameData, prevFrame)
	end

	function start_processing(config::ProcessingConfig, maxFps::Int)
		prevFrame = Vector{RGB}()
		frameTimes = Vector{Float32}()
		maxFrametime = 1 / maxFps
		avgTime = 0

		while true
			start = time()

			prevFrame = update_frame(config)

			push!(frameTimes, time() - start)

			if (length(frameTimes) >= 10)
				deleteat!(frameTimes, 1)
			end

			avgTime = sum(frameTimes) / size(frameTimes, 1)

			if (avgTime < maxFrametime)
				sleep(maxFrametime - avgTime)
			end
		end
	end

	export init_processing, start_processing
end


text:
export idrs, idrs!

using Random

"""
    idrs(A, b; s = 8, kwargs...) -> x, [history]

Same as [`idrs!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
idrs(A, b; kwargs...) = idrs!(zerox(A,b), A, b; kwargs...)

"""
    idrs!(x, A, b; s = 8, kwargs...) -> x, [history]

Solve the problem ``Ax = b`` approximately with IDR(s), where `s` is the dimension of the
shadow space.

# Arguments

- `x`: Initial guess, will be updated in-place;
- `A`: linear operator;
- `b`: right-hand side.

## Keywords

- `s::Integer = 8`: dimension of the shadow space;
- `abstol::Real = zero(real(eltype(b)))`,
  `reltol::Real = sqrt(eps(real(eltype(b))))`: absolute and relative
  tolerance for the stopping condition
  `|r_k| / |r_0| ≤ max(reltol * resnorm, abstol)`, where `r_k = A * x_k - b`
  is the residual in the `k`th iteration;
- `maxiter::Int = size(A, 2)`: maximum number of iterations;
- `log::Bool`: keep track of the residual norm in each iteration;
- `verbose::Bool`: print convergence information during the iterations.

# Return values

**if `log` is `false`**

- `x`: approximate solution.

**if `log` is `true`**

- `x`: approximate solution;
- `history`: convergence history.
"""
function idrs!(x, A, b;
               s = 8,
               abstol::Real = zero(real(eltype(b))),
               reltol::Real = sqrt(eps(real(eltype(b)))),
               maxiter=size(A, 2),
               log::Bool=false,
               kwargs...)
    history = ConvergenceHistory(partial=!log)
    history[:abstol] = abstol
    history[:reltol] = reltol
    log && reserve!(history, :resnorm, maxiter)
    idrs_method!(history, x, A, b, s, abstol, reltol, maxiter; kwargs...)
    log && shrink!(history)
    log ? (x, history) : x
end

#########################
# Method Implementation #
#########################

@inline function omega(t, s)
    angle = sqrt(2.)/2
    ns = norm(s)
    nt = norm(t)
    ts = dot(t,s)
    rho = abs(ts/(nt*ns))
    om = ts/(nt*nt)
    if rho < angle
        om = om*convert(typeof(om),angle)/rho
    end
    om
end

function idrs_method!(log::ConvergenceHistory, X, A, C::T,
    s::Number, abstol::Real, reltol::Real, maxiter::Number; smoothing::Bool=false, verbose::Bool=false
    ) where {T}

    verbose && @printf("=== idrs ===\n%4s\t%7s\n","iter","resnorm")
    R = C - A*X
    normR = norm(R)
    iter = 1
    tol = max(reltol * normR, abstol)

    if smoothing
        X_s = copy(X)
        R_s = copy(R)
        T_s = zero(R)
    end

    if normR <= tol # Initial guess is a good enough solution
        setconv(log, 0<=normR<tol)
        return X
    end

    Z = zero(C)

    P = T[rand!(copy(C)) for k in 1:s]
    U = T[copy(Z) for k in 1:s]
    G = T[copy(Z) for k in 1:s]
    Q = copy(Z)
    V = copy(Z)

    M = Matrix{eltype(C)}(I,s,s)
    f = zeros(eltype(C),s)
    c = zeros(eltype(C),s)

    om::eltype(C) = 1
    while normR > tol && iter ≤ maxiter
        for i in 1:s
            f[i] = dot(P[i], R)
        end
        for k in 1:s
            nextiter!(log,mvps=1)

            # Solve small system and make v orthogonal to P

            c = LowerTriangular(M[k:s,k:s])\f[k:s]
            V .= c[1] .* G[k]
            Q .= c[1] .* U[k]

            for i = k+1:s
                V .+= c[i-k+1] .* G[i]
                Q .+= c[i-k+1] .* U[i]
            end

            # Compute new U[:,k] and G[:,k], G[:,k] is in space G_j
            V .= R .- V

            U[k] .= Q .+ om .* V
            mul!(G[k], A, U[k])

            # Bi-orthogonalise the new basis vectors

            for i in 1:k-1
                alpha = dot(P[i],G[k])/M[i,i]
                G[k] .-= alpha .* G[i]
                U[k] .-= alpha .* U[i]
            end

            # New column of M = P'*G  (first k-1 entries are zero)

            for i in k:s
                M[i,k] = dot(P[i],G[k])
            end

            #  Make r orthogonal to q_i, i = 1..k

            beta = f[k]/M[k,k]
            R .-= beta .* G[k]
            X .+= beta .* U[k]

            normR = norm(R)
            if smoothing
                T_s .= R_s .- R

                gamma = dot(R_s, T_s)/dot(T_s, T_s)

                R_s .-= gamma .* T_s
                X_s .-= gamma .* (X_s .- X)

                normR = norm(R_s)
            end
            push!(log, :resnorm, normR)
            verbose && @printf("%3d\t%1.2e\n",iter,normR)
            if normR < tol || iter == maxiter
                setconv(log, 0<=normR<tol)
                return X
            end
            if k < s
                f[k+1:s] .-=  beta*M[k+1:s,k]
            end
            iter += 1
        end

        # Now we have sufficient vectors in G_j to compute residual in G_j+1
        # Note: r is already perpendicular to P so v = r
        copyto!(V, R)
        mul!(Q, A, V)
        om = omega(Q, R)
        R .-= om .* Q
        X .+= om .* V

        normR = norm(R)
        if smoothing
            T_s .= R_s .- R

            gamma = dot(R_s, T_s)/dot(T_s, T_s)

            R_s .-= gamma .* T_s
            X_s .-= gamma .* (X_s .- X)

            normR = norm(R_s)
        end
        iter += 1
        nextiter!(log, mvps=1)
        push!(log, :resnorm, normR)
    end
    if smoothing
        copyto!(X, X_s)
    end
    verbose && @printf("\n")
    setconv(log, 0<=normR<tol)
    X
end