text:
using KitBase, FluxReconstruction, LinearAlgebra, OrdinaryDiffEq, Plots
using ProgressMeter: @showprogress

begin
    set = Setup(
        "gas",
        "cylinder",
        "2d0f",
        "hll",
        "nothing",
        1, # species
        3, # order of accuracy
        "positivity", # limiter
        "euler",
        0.1, # cfl
        1.0, # time
    )

    ps0 = KitBase.PSpace2D(0.0, 1.0, 20, 0.0, 1.0, 20)
    deg = set.interpOrder - 1
    ps = FRPSpace2D(ps0, deg)

    vs = nothing
    gas = Gas(
        1e-6,
        2.0, # Mach
        1.0,
        1.0, # K
        5 / 3,
        0.81,
        1.0,
        0.5,
    )
    ib = nothing

    ks = SolverSet(set, ps0, vs, gas, ib)
end

u0 = zeros(ps.nx, ps.ny, deg + 1, deg + 1)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    u0[i, j, k, l] = max(
        exp(-300 * ((ps.xpg[i, j, k, l, 1] - 0.5)^2 + (ps.xpg[i, j, k, l, 2] - 0.5)^2)),
        1e-2,
    )
end

u0 = zeros(ps.nx, ps.ny, deg + 1, deg + 1, 4)
for i in axes(u0, 1), j in axes(u0, 2), k in axes(u0, 3), l in axes(u0, 4)
    ρ = max(
        exp(-300 * ((ps.xpg[i, j, k, l, 1] - 0.5)^2 + (ps.xpg[i, j, k, l, 2] - 0.5)^2)),
        1e-2,
    )
    prim = [ρ, 0.0, 0.0, 1.0]
    u0[i, j, k, l, :] .= prim_conserve(prim, ks.gas.γ)
end

function dudt!(du, u, p, t)
    du .= 0.0

    J, ll, lr, dhl, dhr, lpdm, γ = p

    nx = size(u, 1) - 2
    ny = size(u, 2) - 2
    nsp = size(u, 3)

    f = zeros(nx, ny, nsp, nsp, 4, 2)
    for i in axes(f, 1), j in axes(f, 2), k = 1:nsp, l = 1:nsp
        fg, gg = euler_flux(u[i, j, k, l, :], γ)
        for m = 1:4
            f[i, j, k, l, m, :] .= inv(J[i, j][k, l]) * [fg[m], gg[m]]
        end
    end

    u_face = zeros(nx, ny, 4, nsp, 4)
    f_face = zeros(nx, ny, 4, nsp, 4, 2)
    for i in axes(u_face, 1), j in axes(u_face, 2), l = 1:nsp, m = 1:4
        u_face[i, j, 1, l, m] = dot(u[i, j, l, :, m], ll)
        u_face[i, j, 2, l, m] = dot(u[i, j, :, l, m], lr)
        u_face[i, j, 3, l, m] = dot(u[i, j, l, :, m], lr)
        u_face[i, j, 4, l, m] = dot(u[i, j, :, l, m], ll)

        for n = 1:2
            f_face[i, j, 1, l, m, n] = dot(f[i, j, l, :, m, n], ll)
            f_face[i, j, 2, l, m, n] = dot(f[i, j, :, l, m, n], lr)
            f_face[i, j, 3, l, m, n] = dot(f[i, j, l, :, m, n], lr)
            f_face[i, j, 4, l, m, n] = dot(f[i, j, :, l, m, n], ll)
        end
    end

    fx_interaction = zeros(nx + 1, ny, nsp, 4)
    for i = 2:nx, j = 1:ny, k = 1:nsp
        fw = @view fx_interaction[i, j, k, :]
        uL = @view u_face[i-1, j, 2, k, :]
        uR = @view u_face[i, j, 4, k, :]
        flux_hll!(fw, uL, uR, γ, 1.0)
        fw .*= 40

        #fx_interaction[i, j, k, :] .=
        #    0.5 .* (f_face[i-1, j, 2, k, :, 1] .+ f_face[i, j, 4, k, :, 1]) .-
        #    ps.dx[i, j] .* (u_face[i, j, 4, k, :] - u_face[i-1, j, 2, k, :])
    end
    fy_interaction = zeros(nx, ny + 1, nsp, 4)
    for i = 1:nx, j = 2:ny, k = 1:nsp
        fw = @view fy_interaction[i, j, k, :]
        uL = local_frame(u_face[i, j-1, 3, k, :], 0.0, 1.0)
        uR = local_frame(u_face[i, j, 1, k, :], 0.0, 1.0)
        flux_hll!(fw, uL, uR, γ, 1.0)
        fw .= global_frame(fw, 0.0, 1.0) .* 40

        #fy_interaction[i, j, k, :] .=
        #    0.5 .* (f_face[i, j-1, 3, k, :, 2] .+ f_face[i, j, 1, k, :, 2]) .-
        #    ps.dy[i, j] .* (u_face[i, j, 1, k, :] - u_face[i, j-1, 3, k, :])
    end

    rhs1 = zeros(nx, ny, nsp, nsp, 4)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
        rhs1[i, j, k, l, m] = dot(f[i, j, :, l, m, 1], lpdm[k, :])
    end
    rhs2 = zeros(nx, ny, nsp, nsp, 4)
    for i = 1:nx, j = 1:ny, k = 1:nsp, l = 1:nsp, m = 1:4
        rhs2[i, j, k, l, m] = dot(f[i, j, k, :, m, 2], lpdm[l, :])
    end

    for i = 2:nx-1, j = 2:ny-1, k = 1:nsp, l = 1:nsp, m = 1:4
        du[i, j, k, l, m] = -(
            rhs1[i, j, k, l, m] +
            rhs2[i, j, k, l, m] +
            (fx_interaction[i, j, l, m] - f_face[i, j, 4, l, m, 1]) * dhl[k] +
            (fx_interaction[i+1, j, l, m] - f_face[i, j, 2, l, m, 1]) * dhr[k] +
            (fy_interaction[i, j, k, m] - f_face[i, j, 1, k, m, 2]) * dhl[l] +
            (fy_interaction[i, j+1, k, m] - f_face[i, j, 3, k, m, 2]) * dhr[l]
        )
    end

    return nothing
end

tspan = (0.0, 1.0)
p = (ps.J, ps.ll, ps.lr, ps.dhl, ps.dhr, ps.dl, ks.gas.γ)
prob = ODEProblem(dudt!, u0, tspan, p)

dt = 0.005
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), save_everystep = false, adaptive = false, dt = dt)

@showprogress for iter = 1:10
    step!(itg)
end

sol = zeros(ps.nx, ps.ny, 4)
for i = 1:ps.nx, j = 1:ps.ny
    sol[i, j, :] .= conserve_prim(itg.u[i, j, 2, 2, :], ks.gas.γ)
    sol[i, j, 4] = 1 / sol[i, j, 4]
end

contourf(ps.x, ps.y, sol[:, :, 1], aspect_ratio = 1, legend = true)


text:
export AbstractSwitching, ControlledSwitching, AutonomousSwitching

"""
    AbstractSwitching

Nature of the switching, e.g. [`AutonomousSwitching`](@ref) or [`ControlledSwitching`](@ref), see Section 1.1.3 of [1]

[1] Liberzon, D.
*Switching in systems and control*.
Springer Science & Business Media, 2012
"""
abstract type AbstractSwitching end

"""
    AutonomousSwitching <: AbstractSwitching

Controlled switching, the switching signal is autonomous.
"""
struct AutonomousSwitching <: AbstractSwitching end

"""
    ControlledSwitching <: AbstractSwitching

Controlled switching, the switching signal is controlled.
"""
struct ControlledSwitching <: AbstractSwitching end

ext: jl


text:
function check_size(x, p, dims)
    dims = (dims==1) ? 2 : 1
    x_s, p_s = size(x, dims), size(p[:, :], 2)
    if x_s != p_s
        throw(DimensionMismatch("second dimension of the target data is $x_s, and first dimension of scaler's parameter $p_s. these dimensions must match(i.e. Number of features of fitting data and number of features of target data must match)."))
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_startup_cost_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation::MorLatRam2013.StartupCosts,
)::Nothing
    eq_startup_choose = _init(model, :eq_startup_choose)
    eq_startup_restrict = _init(model, :eq_startup_restrict)
    S = length(g.startup_categories)
    startup = model[:startup]
    for t in 1:model[:instance].time
        # If unit is switching on, we must choose a startup category
        eq_startup_choose[g.name, t] = @constraint(
            model,
            model[:switch_on][g.name, t] ==
            sum(startup[g.name, t, s] for s in 1:S)
        )

        for s in 1:S
            # If unit has not switched off in the last `delay` time periods, startup category is forbidden.
            # The last startup category is always allowed.
            if s < S
                range_start = t - g.startup_categories[s+1].delay + 1
                range_end = t - g.startup_categories[s].delay
                range = (range_start:range_end)
                initial_sum = (
                    g.initial_status < 0 && (g.initial_status + 1 in range) ? 1.0 : 0.0
                )
                eq_startup_restrict[g.name, t, s] = @constraint(
                    model,
                    startup[g.name, t, s] <=
                    initial_sum + sum(
                        model[:switch_off][g.name, i] for i in range if i >= 1
                    )
                )
            end

            # Objective function terms for start-up costs
            add_to_expression!(
                model[:obj],
                startup[g.name, t, s],
                g.startup_categories[s].cost,
            )
        end
    end
    return
end


text:
function deactivate!(plank, loss)
    @inbounds plank.ac .*= (1.0 .- loss)
end

##### grazing and grazing loss
function grazing!(plank, arch::Architecture, plk, p)
    ##### calculate grazing loss
    calc_loss!(plk.DOC.data, plk.POC.data, plank, plank.ac, plank.xi, plank.yi, plank.zi,
               plank.graz, p.grazFracC, arch)

    ##### deactivate grazed individuals
    deactivate!(plank, plank.graz)

    return nothing
end

##### mortality and mortality loss
function mortality!(plank, arch::Architecture, plk, p)
    ##### calculate mortality loss
    calc_loss!(plk.DOC.data, plk.POC.data, plank, plank.ac, plank.xi, plank.yi, plank.zi,
               plank.mort, p.mortFracC, arch)

    ##### deactivate dead individuals
    deactivate!(plank, plank.mort)

    return nothing
end

@kernel function get_tind_kernel!(idx, con, con_ind, de_ind)
    i = @index(Global, Linear)
    if con[i] == 1.0
        idx[i] = de_ind[con_ind[i]]
    end
end
function get_tind!(idx, con, con_ind, de_ind, arch)
    kernel! = get_tind_kernel!(device(arch), 256, (size(idx,1)))
    event = kernel!(idx, con, con_ind, de_ind)
    wait(device(arch), event)
    return nothing
end

##### copy ready to divide individuals to inactive rows
@kernel function copy_daughter_individuals_kernel!(plank, con, idx)
    i = @index(Global, Linear)
    if (con[i] == 1.0) & (idx[i] ≠ 0)
        # @print("index: $(idx[i]), $i \n")
        @inbounds plank.x[idx[i]]    = plank.x[i]
        @inbounds plank.y[idx[i]]    = plank.y[i]
        @inbounds plank.z[idx[i]]    = plank.z[i]
        @inbounds plank.Sz[idx[i]]   = plank.Sz[i]
        @inbounds plank.Bm[idx[i]]   = plank.Bm[i]
        @inbounds plank.Chl[idx[i]]  = plank.Chl[i]
        @inbounds plank.gen[idx[i]]  = plank.gen[i]
        @inbounds plank.ac[idx[i]]   = plank.ac[i]
        @inbounds plank.dvid[idx[i]] = plank.dvid[i]
        @inbounds plank.graz[idx[i]] = plank.graz[i]
        @inbounds plank.mort[idx[i]] = plank.mort[i]
    end
end
function copy_daughter_individuals!(plank, con, idx::AbstractArray{Int64,1}, arch)
    kernel! = copy_daughter_individuals_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, con, idx)
    wait(device(arch), event)
    return nothing
end

##### cell division
@kernel function divide_to_half_kernel!(plank)
    i = @index(Global)
    @inbounds plank.Sz[i]  *= (2.0 - plank.dvid[i]) / 2
    @inbounds plank.Bm[i]  *= (2.0 - plank.dvid[i]) / 2
    @inbounds plank.Chl[i] *= (2.0 - plank.dvid[i]) / 2
    @inbounds plank.gen[i] += plank.dvid[i]
    @inbounds plank.age[i] *= (1.0 - plank.dvid[i])
    @inbounds plank.iS[i]   = plank.iS[i] * (1.0 - plank.dvid[i]) + plank.Sz[i] * plank.dvid[i]
    @inbounds plank.Th[i]  *= (1.0 - plank.dvid[i])
end
function divide_to_half!(plank, arch)
    kernel! = divide_to_half_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank)
    wait(device(arch), event)
    return nothing
end
function divide!(plank, deactive_ind, arch::Architecture)
    con_ind = cumsum(plank.dvid)
    get_tind!(plank.idx, plank.dvid, Int.(con_ind), deactive_ind, arch)
    copy_daughter_individuals!(plank, plank.dvid, Int.(plank.idx), arch)
    divide_to_half!(plank, arch)
    # println("time step ends")
    return nothing
end


text:
test("Add path to offset") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = ClipperOffset()

    return add_path!(c, path, JoinTypeRound, EndTypeClosedPolygon)
end

test("Add paths to offset") do
    path1 = Vector{IntPoint}()

    push!(path1, IntPoint(0, 0))
    push!(path1, IntPoint(0, 1))

    path2 = Vector{IntPoint}()

    push!(path2, IntPoint(5, 0))
    push!(path2, IntPoint(2, 1))

    paths = Vector{IntPoint}[path1, path2]

    c = ClipperOffset()

    return add_paths!(c, paths, JoinTypeRound, EndTypeClosedPolygon)
end

test("Clear") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = ClipperOffset()

    add_path!(c, path, JoinTypeRound, EndTypeClosedPolygon)

    return Clipper.clear!(c)
end

test("Offset") do
    path = Vector{IntPoint}()

    push!(path, IntPoint(0, 0))
    push!(path, IntPoint(0, 1))

    c = ClipperOffset()

    add_path!(c, path, JoinTypeRound, EndTypeOpenRound)

    poly = execute(c, 1.0)

    @test poly == [[IntPoint(1, 2), IntPoint(-1, 2), IntPoint(-1, -1), IntPoint(1, -1)]]
end


text:
"""
    Cartpole{T}

Canoncial nonlinear system, with a prismatic joint followed by a revolute joint, or a pendulum
attached to a cart, that is only actuated by sliding the cart. It has 4 states and 1 control.

# Constructors
    Cartpole(; kwargs...)

with keyword arguments
* `mc` - mass of the cart, in kg (default = 1)
* `mp` - mass of the pendulum, in kg (default = 0.2)
* `l` - length of the pendulum, in m (default = 0.5)
* `g` - gravity, in m/s² (default = 9.81)
"""
@autodiff struct Cartpole{T} <: ContinuousDynamics
    mc::T
    mp::T
    l::T
    g::T
    function Cartpole(mc, mp, l, g)
        T = eltype(promote(mc, mp, l, g))
        new{T}(mc, mp, l, g)
    end
end

Cartpole(; mc=1.0, mp=0.2, l=0.5, g=9.81) = Cartpole(mc, mp, l, g)

function dynamics(model::Cartpole, x, u)
    mc = model.mc  # mass of the cart in kg (10)
    mp = model.mp   # mass of the pole (point mass at the end) in kg
    l = model.l   # length of the pole in m
    g = model.g  # gravity m/s^2

    q = x[ @SVector [1,2] ]
    qd = x[ @SVector [3,4] ]

    s = sin(q[2])
    c = cos(q[2])

    H = @SMatrix [mc+mp mp*l*c; mp*l*c mp*l^2]
    C = @SMatrix [0 -mp*qd[2]*l*s; 0 0]
    G = @SVector [0, mp*g*l*s]
    B = @SVector [1, 0]

    qdd = -H\(C*qd + G - B*u[1])
    return [qd; qdd]
end

function dynamics!(model::Cartpole, xdot, x, u)
    xdot .= dynamics(model, x, u)
end

RobotDynamics.state_dim(::Cartpole) = 4
RobotDynamics.control_dim(::Cartpole) = 1

# Base.position(::Cartpole, x::StaticVector) = SA[0,x[1],0]
# RobotDynamics.orientation(::Cartpole, x::StaticVector) = UnitQuaternion(expm(SA[1,0,0]*x[2]))


text:
# test of optBL with the example in STICR panmixia test
# Claudia January 2015

include("types.jl")
include("functions.jl")

df = readtable("/Users/Clauberry/Documents/phylo/software/CFimplementation/sticr/quartetCF.csv")
d2 = readTableCF(df)

net = readTopologyUpdate("/Users/Clauberry/Documents/phylo/software/CFimplementation/sticr/tree.tre")
# fixit: cannot read _ inside the taxon name, need to change readsubtree
printEdges(net)

net.ht
realht = [0.1,0.2,0.1,1.0]

@time fmin,xmin=optBL(net,d2)


text:
################################################################################
# this implements the right pseudoinverse
# is defined if A has linearly independent columns
# ⁻¹, ⁺ syntax
struct PseudoInverse{T, M} <: Factorization{T}
    parent::M
end
PseudoInverse(A) = PseudoInverse{eltype(A), typeof(A)}(A)

Base.size(P::PseudoInverse) = size(P.parent')
Base.size(P::PseudoInverse, k::Integer) = size(P.parent', k::Integer)
function Base.AbstractMatrix(P::PseudoInverse)
	A = P.parent
    A isa Number ? fill(inv(A), 1, 1) : (A \ I(size(A, 1)))
end
function Base.Matrix(Inv::PseudoInverse)
	M = AbstractMatrix(Inv)
	M isa Matrix ? M : Matrix(M) # since it could be e.g. a Diagonal
end
# Base.Matrix(P::PseudoInverse) = AbstractMatrix(P)
LinearAlgebra.factorize(P::PseudoInverse) = P # same reasoning as for Inverse

# smart constructor
# calls regular inverse if matrix is square
function pseudoinverse end
const pinverse = pseudoinverse
function pseudoinverse(A::AbstractMatOrFac)
    size(A, 1) == size(A, 2) ? inverse(A) : PseudoInverse(A)
end
pseudoinverse(A::Union{Number, UniformScaling}) = inv(A)
pseudoinverse(P::PseudoInverse) = P.parent
pseudoinverse(P::Inverse) = inverse(P) # in case P.parent is scalar


text:
macro nogc(ex)
  quote
    local val
    try
      gc_enable(false)
      val = $(esc(ex))
    finally
      gc_enable(true)
    end
    val
  end
end

max_stars_repo_name: xrounder/optunity

max_issues_repo_name: xrounder/optunity

max_forks_repo_name: xrounder/optunity


text:
import Pkg
Pkg.add("Plots")
Pkg.add("FFTW")
Pkg.add("DSP")
Pkg.add("WAV")
Pkg.add("LibSndFile")


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

# implement methods for spatial statistics
import Statistics: mean, var, quantile

"""
    SpatialStatistic

A spatial statistic defined over a spatial domain.
"""
struct SpatialStatistic{T,N,D<:AbstractDomain{T,N}} <: AbstractData{T,N}
  data::Dict{Symbol,<:AbstractArray}
  domain::D
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, statistic::SpatialStatistic)
  N = ndims(domain(statistic))
  print(io, "$(N)D SpatialStatistic")
end

function Base.show(io::IO, ::MIME"text/plain", statistic::SpatialStatistic)
  println(io, statistic)
  println(io, "  domain: ", domain(statistic))
  print(  io, "  variables: ", join(keys(variables(statistic)), ", ", " and "))
end

include("statistics/data.jl")
include("statistics/solutions.jl")


text:
module Ls

import LsoBase
import Obj
import Obj.Objective

"""
    LineSearch(name, strategy)

    Creates an LineSearch object for the given strategy method.
    The strategy method should have the following signature:

    (w, s, obj, b, fw, gw) -> (α, iter)

    Where:
       w: The current point
       s: The current step direction
     obj: The objective function
       b: Index batch to evaluate objective function on. May be ignored by
          the implementation.
      fw: f(w) (present for performance reasons)
      gw: g(w) (present for performance reasons)
       α: The resulting step size
    iter: The number of LS iterations required
"""
type LineSearch
        name::ASCIIString
    strategy::Function
end

"""
    ls(linesearch, w, s [, obj, b, fw, gw])

    Performs a Line Search with the given strategy for current point w and step
    direction s. When the LS strategy depends on an objective function value, that
    function should be specified using the obj parameter.

    Optionally, b can specify a batch of indices to evaluate the objective function
    on. For performance reasons, f(w) and g(w) can be provided as fw and gw.

    This helper function should ease the access to objects of type Ls.LineSearch.
"""
function ls(linesearch::LineSearch, w::Array{Float64,1}, s::Array{Float64,1},
            obj::Union{Objective,Void}=nothing, b::Array{Int,1}=Int[],
            fw::Float64=NaN, gw::Array{Float64,1}=Float64[])

    # sanitize arguments
    if obj != nothing
      if isnan(fw)
          if length(b) == 0
              fw = Obj.f(obj, w)
          else
              fw = Obj.f(obj, w, b)
          end
      end
      if length(gw) == 0
          if length(b) == 0
              gw = Obj.g(obj, w)
          else
              gw = Obj.g(obj, w, b)
          end
      end
    end

    # run ls
    return linesearch.strategy(w, s, obj, b, fw, gw)

end

# include ./opt/*
LsoBase.includedir(dirname(@__FILE__)*"/ls")

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
using Plots
using HDF5

dir = "StochasticNeighborhoodEmbeddings/examples/"
f = h5open(dir * "MNIST_TSNE_full.h5", "r")
labels = read(f["labels"])
Y = read(f["Y"])
close(f)

YC = [Y[:, ==(i).(labels)] for i in 1:10]
scatter()
for (i, yc) in enumerate(YC)
    scatter!(yc[1, :], yc[2, :], label = "$(i-1)")
end
gui()


text:
include("UWavelets.jl")
include("VWavelets.jl")
include("preconditioners.jl")
include("update_preconditioners.jl")
include("update_wavelets.jl")


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

function updateStationStats!(station::Station; numIdleAmbsChange::Int = 0, time::Float = nullTime)
	if numIdleAmbsChange != 0
		@assert(time != nullTime)
		@assert(time >= station.currentNumIdleAmbsSetTime)

		# stats
		station.numIdleAmbsTotalDuration[station.currentNumIdleAmbs] += time - station.currentNumIdleAmbsSetTime
		station.currentNumIdleAmbs += numIdleAmbsChange
		station.currentNumIdleAmbsSetTime = time
	end
end


text:
module Analyze
    #using Plots
    using Requires
    using Distributions

    function analyze(dirname)
        return Analyze_data(dirname)
    end

    mutable struct Analyze_data
        plaquette::Array{Float64,1}
        polyakov::Array{ComplexF64,1}
        trjs::Array{Int64,1}
        numdata::Int64
        dirname::String

        function Analyze_data(dirname)
            plaquette = Float64[]
            polyakov = ComplexF64[]
            trjs = Int64[]
            filename_plaquette = pwd()*"/"*dirname*"/Plaquette.txt"
            filename_polyakov = pwd()*"/"*dirname*"/Polyakov_loop.txt"
            if isfile(filename_plaquette)
                open(filename_plaquette) do f
                    for i in eachline(f)
                        if occursin("plaq",i)
                            tmp = split(i)
                            itrj = parse(Int,tmp[1])
                            plaq = parse(Float64,tmp[2])
                            append!(trjs,itrj)
                            append!(plaquette,plaq)
                        end

                    end
                end
            else
                println("$filename_plaquette does not exist")
            end

            if isfile(filename_polyakov)
                open(filename_polyakov) do f
                    for i in eachline(f)
                        if occursin("poly",i)
                            tmp = split(i)
                            itrj = parse(Int,tmp[1])
                            poly = parse(Float64,tmp[2])+im*parse(Float64,tmp[3])
                            append!(polyakov,poly)
                        end

                    end
                end
            else
                println("$filename_polyakov does not exist")
            end

            return new(plaquette,polyakov,trjs,length(trjs),dirname)
        end
    end

    function get_plaquette_average(dirname::String)
        return mean( get_plaquette(dirname) )
    end

    function get_polyakov_average(dirname::String)
        return mean( get_polyakov(dirname) )
    end

    function get_plaquette(dirname::String)
        a = analyze(dirname)
        return a.plaquette
    end

    function get_polyakov(dirname::String)
        a = analyze(dirname)
        return a.polyakov
    end

    function get_trjs(dirname::String)
        a = analyze(dirname)
        return a.trjs
    end

    function get_plaquette(a::Analyze_data)
        return a.plaquette
    end

    function get_polyakov(a::Analyze_data)
        return a.polyakov
    end

    function get_trjs(a::Analyze_data)
        return a.trjs
    end

    function get_dirname(a::Analyze_data)
        return a.dirname
    end

    function __init__()
        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
            using Plots
            function plot_plaquette(dirname::String)
                a = analyze(dirname)
                plot_plaquette(a)
            end

            function plot_plaquette(a::Analyze_data)
                plt1 = histogram([0],label=nothing) #plot1
                plt2 = plot([], [],label=nothing) #plot2
                ylabel!("Plaquette")
                xlabel!("MC time")
                plot(plt1, plt2,  layout = 2)
                hist_plaq=[]

                Nsteps = a.numdata

                function plot_refresh!(plt1,plt2,hist_plaq,plaq,itrj)
                    bins=round(Int,log(itrj)*5+1)
                    append!(hist_plaq, plaq)

                    plt1 = histogram(hist_plaq,bins=bins,label=nothing) #plot1
                    xlabel!("Plaquette")

                    push!(plt2, 1, itrj, plaq)

                    plot(plt1, plt2, layout = 2)

                end

                for itrj=1:Nsteps
                    plaq = a.plaquette[itrj]

                    plot_refresh!(plt1,plt2,hist_plaq,plaq,itrj)

                    flush(stdout)
                end

                savefig(pwd()*"/"*a.dirname*"/Plaquette.pdf")
            end

            function plot_polyakov(dirname::String)
                a = analyze(dirname)
                plot_polyakov(a)
            end

            function plot_polyakov(a::Analyze_data)
                plt1 = histogram([0],label=nothing) #plot1
                plt2 = plot([], [],label=nothing) #plot2
                ylabel!("Polyakov loop")
                xlabel!("MC time")
                plot(plt1, plt2,  layout = 2)
                hist_poly=[]

                Nsteps = a.numdata

                function plot_refresh!(plt1,plt2,hist_poly,poly,itrj)
                    bins=round(Int,log(itrj)*5+1)
                    append!(hist_poly,poly)

                    plt1 = histogram(hist_poly,bins=bins,label=nothing) #plot1
                    xlabel!("Polyakov loop")

                    push!(plt2, 1, itrj, poly)

                    plot(plt1, plt2, layout = 2)

                end

                for itrj=1:Nsteps
                    poly = a.polyakov[itrj]
                    plot_refresh!(plt1,plt2,hist_poly,abs(poly),itrj)

                    flush(stdout)
                end
                savefig(pwd()*"/"*a.dirname*"/Polyakov.pdf")
            end

            function plot_plaq_and_poly(dirname::String)
                a = analyze(dirname)
                plot_plaq_and_poly(a)
            end

            function plot_plaq_and_poly(a::Analyze_data)
                plt1 = histogram([0],label=nothing) #plot1
                plt2 = plot([], [],label=nothing) #plot2
                ylabel!("Plaquette")
                xlabel!("MC time")
                plt3 = histogram([0],label=nothing) #plot3
                plt4 = plot([], [],label=nothing) #plot4
                ylabel!("Polyakov loop")
                xlabel!("MC time")
                plot(plt1, plt2, plt3, plt4, layout = 4)

                hist_plaq=[]
                hist_poly=[]

                function plot_refresh!(plt1,plt2,plt3,plt4,hist_plaq,hist_poly, plaq, poly,itrj)
                    bins=round(Int,log(itrj)*5+1)
                    append!(hist_plaq, plaq)
                    append!(hist_poly, poly)
                    #
                    plt1 = histogram(hist_plaq,bins=bins,label=nothing) #plot1
                    xlabel!("Plaquette")
                    plt3 = histogram(hist_poly,bins=bins,label=nothing) #plot3
                    xlabel!("Polyakov loop")
                    #
                    push!(plt2, 1, itrj, plaq)
                    push!(plt4, 1, itrj, poly)
                    #
                    plot(plt1, plt2, plt3, plt4, layout = 4)
                    #gui()
                end

                Nsteps = a.numdata

                for itrj=1:Nsteps
                    poly = a.polyakov[itrj]
                    plaq = a.plaquette[itrj]
                    plot_refresh!(plt1,plt2,plt3,plt4,hist_plaq,hist_poly, plaq, abs(poly),itrj)

                    flush(stdout)
                end
                savefig(pwd()*"/"*a.dirname*"/Plaquette_Polyakov.pdf")
            end
        end
    end

end


text:
module RenderingGeometry
    using StaticArrays
    using Rematch
    include("utils.jl")
    include("component_geometry/lib.jl")
    include("ray.jl")
    include("bounds/bounds.jl")
    include("bounds/bounds_behavior.jl")
end


text:
using DWARF

tag_filter(which_tag, arr) = filter(x->DWARF.tag(x) == which_tag, arr)

function investigate2(oh; dumpsource = false, dumpdies = false, dumploctables = false)
    sects = ELF.Sections(oh)

    # Relocate the object file
    new_buf = copy(oh.io)
    seekstart(new_buf)
    ELF.relocate!(new_buf,oh)
    seekstart(new_buf)
    oh = readmeta(new_buf)

    dbgs = debugsections(oh);
    tree = read(oh, DWARF.DIETree)
    strtab = ELF.load_strtab(dbgs.debug_str)

    dumpdies && show(STDOUT, tree; strtab = strtab)

    subprograms = collect(tag_filter(DWARF.DW_TAG_subprogram, tree.children))
    linkage_names = map(subprograms) do SP
        linkage = bytestring(first(tag_filter(DWARF.DW_AT_linkage_name, DWARF.attributes(SP))),strtab)
    end
    names = map(subprograms) do SP
        linkage = bytestring(first(tag_filter(DWARF.DW_AT_name, DWARF.attributes(SP))),strtab)
    end
    llvmf = map(FindFunctionNamed,linkage_names)
    juliaf = map(x->x == C_NULL ? nothing : getLam(x),llvmf)
    if dumpsource
        for i = 1:length(llvmf)
            print(Base.uncompressed_ast(juliaf[i]).args[3])
            icxx" $(llvmf[i])->getParent()->dump(); "
        end
    end
    for i = 1:length(juliaf)
        println("Debug Info for ",names[i])
        if juliaf[i] == nothing
            println("AST not found")
            continue
        end
        ast = Base.uncompressed_ast(juliaf[i])
        println("Parameters:")
        params = ast.args[1]
        for var in params
            println(" - ",var.args[1])
        end
        lvars = ast.args[2][1]
        if !isempty(lvars)
            println("Local Variables:")
            for var in lvars
                println(" - ",var[1])
            end
        end
    end
    if dumploctables
        seek(oh, sectionoffset(dbgs.debug_loc))
        while position(oh) < sectionoffset(dbgs.debug_loc)+sectionsize(dbgs.debug_loc)
            show(read(oh, DWARF.LocationList{UInt64}))
        end
    end
    text = first(filter(x->sectionname(x)==".text",ELF.Sections(oh)))
    disassemble2(sectionaddress(text),sectionsize(text); rangelists=rangelists(oh))
    (llvmf, juliaf)
end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
# This file was generated, do not modify it. # hide
data = dataset("ISLR", "OJ")

X = select(data, [:PriceCH, :PriceMM, :DiscCH, :DiscMM, :SalePriceMM,
                  :SalePriceCH, :PriceDiff, :PctDiscMM, :PctDiscCH]);


text:
@with_gap @testset "Lattices for symmetric group" begin
  # Look at the hook partition [2, 1^(n-1)] of S_{n + 1}.
  # Will have dimension n

  # n = 4
  n = 4
  g1 = matrix(QQ, 4, 4, [ 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 1, -1, 1, -1 ])
  g2 = matrix(QQ, 4, 4, [ -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, -1, 1, -1, 1 ])
  A = matrix_algebra(QQ, 4)
  O = Order(A, [g1, g2])
  L = Hecke.natural_lattice(O)
  for p in [2, 3, 5, 7]
    S = Hecke.sublattice_classes(L, p)
    @test length(S) == valuation(4 + 1, p) + 1
  end

  @test_throws ErrorException isirreducible(L.V)

  n = 5
  g1 = matrix(QQ, 5, 5, [ 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, -1, 1, -1, 1])
  g2 = matrix(QQ, 5, 5, [ -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 1, -1, 1, -1, 1 ])
  A = matrix_algebra(QQ, 5)
  O = Order(A, [g1, g2])
  L = Hecke.natural_lattice(O)
  for p in [2, 3, 5, 7]
    S = Hecke.sublattice_classes(L, p)
    @test length(S) == valuation(5 + 1, p) + 1
  end
end

@with_gap @testset "Lattices for matrix order" begin
  # provided by Gabi Nebe
  A = matrix(QQ, 3, 3, [1,3^4,0,0,1,0,0,0,1]);
  B = matrix(QQ, 3, 3, [1,0,0,27,1,0,0,0,1]);
  C = matrix(QQ, 3, 3, [125,3^4,0,3*18,125-90,0,0,0,1]);
  D = matrix(QQ, 3, 3, [1,0,0,0,1,0,0,0,-1]);
  E = matrix(QQ, 3, 3, [1,0,3^4,0,1,0,0,0,1]);
  F = matrix(QQ, 3, 3, [1,0,0,0,1,3^5,0,0,1]);
  G = matrix(QQ, 3, 3, [1,0,0,0,1,0,3^3,0,1]);
  H = matrix(QQ, 3, 3, [1,0,0,0,1,0,0,9,1]);
  l = [A, B, C, D, E, F, G, H]
  A = matrix_algebra(QQ, 3)
  O = Order(A, l)
  L = Hecke.natural_lattice(O)
  S = Hecke.sublattices(L, 3)
  @time length(S) == 216
end


text:
@testset "Neighborhoods" begin
  @testset "BallNeighborhood" begin
    # Euclidean metric
    n = BallNeighborhood{1}(0.5)
    @test isneighbor(n, [0.], [0.])
    @test !isneighbor(n, [0.], [1.])

    n = BallNeighborhood{2}(1.0)
    @test isneighbor(n, [0.,0.], [0.,0.])
    @test isneighbor(n, [0.,0.], [1.,0.])
    @test isneighbor(n, [0.,0.], [0.,1.])

    # Chebyshev metric
    n = BallNeighborhood{1}(0.5, Chebyshev())
    @test isneighbor(n, [0.], [0.])
    @test !isneighbor(n, [0.], [1.])

    for r in [1.,2.,3.,4.,5.]
      n = BallNeighborhood{2}(r, Chebyshev())
      for i in 0.0:1.0:r, j in 0.0:1.0:r
        @test isneighbor(n, [0.,0.], [i,j])
      end
    end

    n = BallNeighborhood{2}(1.0)
    @test sprint(show, n) == "2D BallNeighborhood(1.0, Euclidean(0.0))"
    @test sprint(show, MIME"text/plain"(), n) == "2D BallNeighborhood\n  radius: 1.0\n  metric: Euclidean(0.0)"
  end

  @testset "CylinderNeighborhood" begin
    g = RegularGrid{Float64}(3,3,3)
    n = CylinderNeighborhood{3}(1., 1.)
    for i in [1,2,4,10,11,13]
      @test isneighbor(n, coordinates(g, 1), coordinates(g, i))
    end
    n = CylinderNeighborhood{3}(2., 1.)
    for i in [1,2,3,4,5,7,10,11,12,13,14,16]
      @test isneighbor(n, coordinates(g, 1), coordinates(g, i))
    end
    n = CylinderNeighborhood{3}(1., 2.)
    for i in [1,2,4,10,11,13,19,20,22]
      @test isneighbor(n, coordinates(g, 1), coordinates(g, i))
    end

    g = RegularGrid{Float64}(100,100,100)
    n = CylinderNeighborhood{3}(1., 1.)
    for i in [1,2,101,10001,10002,10101]
      @test isneighbor(n, coordinates(g, 1), coordinates(g, i))
    end

    n = CylinderNeighborhood{3}(1., 1.)
    @test sprint(show, n) == "3D CylinderNeighborhood(1.0, 1.0)"
    @test sprint(show, MIME"text/plain"(), n) == "3D CylinderNeighborhood\n  radius: 1.0\n  height: 1.0"
  end
end


text:
function wrap_fun!(ctx::AbstractContext, cursor::CLFunctionDecl)
    func_type = type(cursor)
    func_name = isempty(ctx.force_name) ? Symbol(spelling(cursor)) : ctx.force_name
    ret_type = :Cint
    args = function_args(cursor)
    arg_types = [argtype(func_type, i) for i in 0:length(args)-1]
    arg_reps = clang2julia.(arg_types)

    for (i, arg) in enumerate(arg_reps)
        if Meta.isexpr(arg, :curly) && first(arg.args) == :NTuple
            arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
        end
        if Meta.isexpr(arg, :curly) && first(arg.args) == :Ptr
            arg_reps[i] = Expr(:curly, :Ref, last(arg.args))
        end
    end

    arg_count = 0
    arg_names = map(args) do x
        n = name_safe(name(x))
        s = !isempty(n) ? n : "arg"*string(arg_count+=1)
        Symbol(s)
    end

    signature = efunsig(Symbol(spelling(cursor)[5:end]), arg_names, arg_reps)
    body = eccall(func_name, Symbol(ctx.libname), ret_type, arg_names, arg_reps)
    push!(ctx.api_buffer, Expr(:function, signature, Expr(:block, body)))

    return ctx
end

function eccall(func_name::Symbol, libname::Symbol, rtype, args, types)
    :(ccall(($(QuoteNode(func_name)), $(QuoteNode(libname))),
        $rtype,
        $(Expr(:tuple, types...)),
        $(args...))
     )
end

function is_ptr_type_expr(@nospecialize t)
    (t === :Cstring || t === :Cwstring) && return true
    isa(t, Expr) || return false
    t = t::Expr
    t.head === :curly && t.args[1] === :Ptr
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = [is_ptr_type_expr(t) ? a : Expr(:(::), a, t) for (a,t) in zip(args,types)]
    Expr(:call, name, x...)
end

function get_comment(cursor::CLCursor)
    stringptr = clang_Cursor_getRawCommentText(cursor.cursor)
    if stringptr.data == C_NULL
        return nothing
    else
        cstring = clang_getCString(stringptr)
        return unsafe_string(cstring)
    end
end

function sort_cursors(cursors::Vector{CLMacroDefinition})
    cursor_names = spelling.(cursors)
    cursor_dict = Dict(zip(cursor_names, cursors))
    sort!(unique!(cursor_names))
    sorted_cursors = [cursor_dict[cursor_names[j]] for j in 1:length(cursor_names)]
    return sorted_cursors
end


text:
CategoricalPool{T, R, V}(ordered::Bool=false) where {T, R, V} =
    CategoricalPool{T, R, V}(T[], ordered)
CategoricalPool{T, R}(ordered::Bool=false) where {T, R} =
    CategoricalPool{T, R}(T[], ordered)
CategoricalPool{T}(ordered::Bool=false) where {T} =
    CategoricalPool{T, DefaultRefType}(T[], ordered)

CategoricalPool{T, R}(levels::Vector, ordered::Bool=false) where {T, R} =
    CategoricalPool{T, R, CategoricalValue{T, R}}(convert(Vector{T}, levels), ordered)
CategoricalPool(levels::Vector{T}, ordered::Bool=false) where {T} =
    CategoricalPool{T, DefaultRefType}(convert(Vector{T}, levels), ordered)

CategoricalPool(invindex::Dict{T, R}, ordered::Bool=false) where {T, R <: Integer} =
    CategoricalPool{T, R, CategoricalValue{T, R}}(invindex, ordered)

Base.convert(::Type{T}, pool::T) where {T <: CategoricalPool} = pool

Base.convert(::Type{CategoricalPool{S}}, pool::CategoricalPool{T, R}) where {S, T, R <: Integer} =
    convert(CategoricalPool{S, R}, pool)

function Base.convert(::Type{CategoricalPool{T, R}}, pool::CategoricalPool) where {T, R <: Integer}
    if length(levels(pool)) > typemax(R)
        throw(LevelsException{T, R}(levels(pool)[typemax(R)+1:end]))
    end

    levelsT = convert(Vector{T}, pool.levels)
    invindexT = convert(Dict{T, R}, pool.invindex)
    return CategoricalPool{T, R, CategoricalValue{T, R}}(levelsT, invindexT, pool.ordered)
end

Base.copy(pool::CategoricalPool{T, R, V}) where {T, R, V} =
    CategoricalPool{T, R, V}(copy(pool.levels), copy(pool.invindex), pool.ordered)

function Base.show(io::IO, pool::CategoricalPool{T, R}) where {T, R}
    @printf(io, "%s{%s,%s}([%s])", typeof(pool).name, T, R,
            join(map(repr, levels(pool)), ", "))

    pool.ordered && print(io, " with ordered levels")
end

Base.length(pool::CategoricalPool) = length(pool.levels)

Base.getindex(pool::CategoricalPool, i::Integer) = pool.valindex[i]
Base.get(pool::CategoricalPool, level::Any) = pool.invindex[level]
Base.get(pool::CategoricalPool, level::Any, default::Any) = get(pool.invindex, level, default)

"""
add the returned value to pool.invindex, this function doesn't do this itself to
avoid doing a dict lookup twice
"""
@inline function push_level!(pool::CategoricalPool{T, R}, level) where {T, R}
    x = convert(T, level)
    n = length(pool)
    if n >= typemax(R)
        throw(LevelsException{T, R}([level]))
    end

    i = R(n + 1)
    push!(pool.levels, x)
    push!(pool.valindex, CategoricalValue(i, pool))
    i
end

function mergelevels(ordered, levels...)
    T = Base.promote_eltype(levels...)
    res = Vector{T}(undef, 0)

    nonempty_lv = findfirst(!isempty, levels)
    if nonempty_lv === nothing
        # no levels
        return res, ordered
    elseif all(l -> isempty(l) || l == levels[nonempty_lv], levels)
        # Fast path if all non-empty levels are equal
        append!(res, levels[nonempty_lv])
        return res, ordered
    end

    for l in levels
        levelsmap = indexin(l, res)

        i = length(res)+1
        for j = length(l):-1:1
            @static if VERSION >= v"0.7.0-DEV.3627"
                if levelsmap[j] === nothing
                    insert!(res, i, l[j])
                else
                    i = levelsmap[j]
                end
            else
                if levelsmap[j] == 0
                    insert!(res, i, l[j])
                else
                    i = levelsmap[j]
                end
            end
        end
    end

    # Check that result is ordered
    if ordered
        levelsmaps = [indexin(res, l) for l in levels]

        # Check that each original order is preserved
        for m in levelsmaps
            issorted(Iterators.filter(x -> x != nothing, m)) || return res, false
        end

        # Check that all order relations between pairs of subsequent elements
        # are defined in at least one set of original levels
        pairs = fill(false, length(res)-1)
        for m in levelsmaps
            @inbounds for i in eachindex(pairs)
                pairs[i] |= (m[i] != nothing) & (m[i+1] != nothing)
            end
            all(pairs) && return res, true
        end
    end

    res, false
end

@inline function Base.get!(pool::CategoricalPool, level::Any)
    get!(pool.invindex, level) do
        if isordered(pool)
            throw(OrderedLevelsException(level, pool.levels))
        end

        push_level!(pool, level)
    end
end

@inline function Base.get!(pool::CategoricalPool, level::CategoricalValue)
    pool === level.pool && return level.level
    # TODO: use a global table to cache subset relations for all pairs of pools
    if level.pool ⊈ pool
        if isordered(pool)
            throw(OrderedLevelsException(level, pool.levels))
        end
        newlevs, ordered = mergelevels(isordered(pool), pool.levels, level.pool.levels)
        # Exception: empty pool marked as ordered if new value is ordered
        if length(pool) == 0 && isordered(level.pool)
            ordered!(pool, true)
        end
        levels!(pool, newlevs)
    end
    get!(pool, get(level))
end

@inline function Base.push!(pool::CategoricalPool, level)
    get!(pool.invindex, level) do
        push_level!(pool, level)
    end
    return pool
end

# TODO: optimize for multiple additions
function Base.append!(pool::CategoricalPool, levels)
    for level in levels
        push!(pool, level)
    end
    return pool
end

# Do not override Base.merge as for internal use we need to use the type and orderedness
# of the first pool rather than promoting both pools
function merge_pools(a::CategoricalPool{T, R}, b::CategoricalPool) where {T, R}
    if length(a) == 0 && length(b) == 0
        newlevs = T[]
        ordered = isordered(a)
    elseif length(a) == 0
        newlevs = Vector{T}(levels(b))
        ordered = isordered(b)
    elseif length(b) == 0
        newlevs = copy(levels(a))
        ordered = isordered(a)
    else
        nl, ordered = mergelevels(isordered(a), a.levels, b.levels)
        newlevs = convert(Vector{T}, nl)
    end
    CategoricalPool{T, R}(newlevs, ordered)
end

Base.issubset(a::CategoricalPool, b::CategoricalPool) = issubset(a.levels, keys(b.invindex))

# Contrary to the CategoricalArray one, this method only allows adding new levels at the end
# so that existing CategoricalValue objects still point to the same value
function levels!(pool::CategoricalPool{S, R}, newlevels::Vector) where {S, R}
    levs = convert(Vector{S}, newlevels)
    if !allunique(levs)
        throw(ArgumentError(string("duplicated levels found in levs: ",
                                   join(unique(filter(x->sum(levs.==x)>1, levs)), ", "))))
    elseif length(levs) < length(pool) || view(levs, 1:length(pool)) != pool.levels
        throw(ArgumentError("removing or reordering levels of existing CategoricalPool is not allowed"))
    end

    n = length(levs)

    if n > typemax(R)
        throw(LevelsException{S, R}(setdiff(levs, levels(pool))[typemax(R)-length(levels(pool))+1:end]))
    end

    empty!(pool.invindex)
    resize!(pool.levels, n)
    resize!(pool.valindex, n)
    for i in 1:n
        v = levs[i]
        pool.levels[i] = v
        pool.invindex[v] = i
        pool.valindex[i] = CategoricalValue(i, pool)
    end

    return pool
end

DataAPI.levels(pool::CategoricalPool) = pool.levels

isordered(pool::CategoricalPool) = pool.ordered
ordered!(pool::CategoricalPool, ordered) = (pool.ordered = ordered; pool)

# LevelsException
function Base.showerror(io::IO, err::LevelsException{T, R}) where {T, R}
    levs = join(repr.(err.levels), ", ", " and ")
    print(io, "cannot store level(s) $levs since reference type $R can only hold $(typemax(R)) levels. Use the decompress function to make room for more levels.")
end

# OrderedLevelsException
function Base.showerror(io::IO, err::OrderedLevelsException)
    print(io, "cannot add new level $(err.newlevel) since ordered pools cannot be extended implicitly. Use the levels! function to set new levels, or the ordered! function to mark the pool as unordered.")
end


text:
###
# This file contains BLAS/BandedMatrix overrides for operators
# that depend on the structure of BandedMatrix
####

#####
# Conversions
#####

function BandedMatrix(S::SubOperator{T,ConcreteConversion{Chebyshev{DD,RR},Ultraspherical{Int,DD,RR},T},
                              Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,DD,RR}
    # we can assume order is 1
    ret = BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    @assert -bandwidth(ret,1) ≤ dg ≤ bandwidth(ret,2)-2

    ret[band(dg)] .= 0.5
    ret[band(dg+1)] .= 0.0
    ret[band(dg+2)] .= -0.5

    # correct first entry
    if 1 in kr && 1 in jr
        ret[1,1] = 1.0
    end

    ret
end

function BandedMatrix(V::SubOperator{T,ConcreteConversion{Ultraspherical{LT,DD,RR},Ultraspherical{LT,DD,RR},T},
                                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,LT,DD,RR}

    n,m = size(V)
    V_l, V_u = bandwidths(V)
    ret = BandedMatrix{eltype(V)}(undef, (n,m), (V_l,V_u))
    kr,jr = parentindices(V)

    (isempty(kr) || isempty(jr)) && return ret

    dg = diagindshift(V)

    λ = order(rangespace(parent(V)))
    c = λ-one(T)

    # need to drop columns

    1-n ≤ dg ≤ m-1 && (ret[band(dg)] .= c./(jr[max(0,dg)+1:min(n+dg,m)] .- 2 .+ λ))
    1-n ≤ dg+1 ≤ m-1 && (ret[band(dg+1)] .= 0)
    1-n ≤ dg+2 ≤ m-1 && (ret[band(dg+2)] .= c./(2 .- λ .- jr[max(0,dg+2)+1:min(n+dg+2,m)]))

    ret
end

#####
# Derivatives
#####

function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Chebyshev{DD,RR},K,T},
                                                     Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR}

    n,m = size(S)
    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    D = parent(S)
    k = D.order
    d = domain(D)

    C=convert(T,pochhammer(one(T),k-1)/2*(4/(complexlength(d)))^k)

    # need to drop columns

    if 1-n ≤ dg+k ≤ m-1
        ret[band(dg+k)] .= C.*(jr[max(0,dg+k)+1:min(n+dg+k,m)] .- one(T))
    end

    ret
end

function BandedMatrix(S::SubOperator{T,ConcreteDerivative{Ultraspherical{LT,DD,RR},K,T},
                                                  Tuple{UnitRange{Int},UnitRange{Int}}}) where {T,K,DD,RR,LT}
    n,m = size(S)
    ret = BandedMatrix{eltype(S)}(undef, (n,m), bandwidths(S))
    kr,jr = parentindices(S)
    dg = diagindshift(S)

    D = parent(S)
    k = D.order
    λ = order(domainspace(D))
    d = domain(D)

    C = convert(T,pochhammer(one(T)*λ,k)*(4/(complexlength(d)))^k)
    ret[band(dg+k)] .= C

    ret
end


text:
using ArgParse

include("functions_bootstrap.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--infile", "-i"
            help = "the alignment to bootstrap, in fasta format"
            required = true
        "--out-prefix", "-p"
            help = "prefix for output filenames"
            default = "bootstrap_"
        "--n-bootstraps", "-n"
            help = "how many bootstraps to carry out"
            arg_type = Int
            default = 1
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    bootstrap(parsed_args["infile"], parsed_args["n-bootstraps"], parsed_args["out-prefix"])
end

main()

max_stars_repo_name: cov-ert/julialign

max_issues_repo_name: cov-ert/julialign

max_forks_repo_name: cov-ert/julialign


text:
    export dummyQZ

    """
    `Basics.read()`  ... reads in data from different files and sources.

    + `("CSF list: Grasp92", "cslFilename")`
        ... reads in the CSF list as given by the (existing) .csl file cslFilename; an basis::Basis is returned with
            basis.isDefined = true and with a proper sequence of orbitals, but with basis.orbitals = Orbital[].

    + `("orbital list: Grasp92", "orbFilename")`
        ... read in the orbitals from a (formatted) .rwf file orbFilename; a list of orbitals::Array{OrbitalsR,1}
            is returned with all subfields specified.
    """
    function Basics.read(sa::String, filename::String)
        if      sa == "CSF list: Grasp92"                          wa = readCslFileGrasp92(filename)
        elseif  sa == "orbital list: Grasp92"                      wa = readOrbitalFileGrasp92(filename)
        elseif  sa == "energies & mixing coefficients: Grasp92"    wa = readMixingFileGrasp92(filename)
        else    error("Unsupported keystring = $sa ")
        end

        return( wa )
    end

    """
    Basics.readCslFileGrasp92(filename::String)`
        ... reads in the CSF list from a Grasp92 .csl file; a basis::Basis is returned.
    """
    function Basics.readCslFileGrasp92(filename::String)
        coreSubshells = Subshell[];    peelSubshells =  Subshell[]

        f  = open(filename)
        sa = readline(f);   sa[1:14] != "Core subshells"   &&   error("Not a Grasp92 .cls file.")
        sa = readline(f);   if  length(sa) > 3    go = true    else  go = false    end
        i  = -1
        while go
            i = i + 1;    sh = Basics.subshellGrasp( strip(sa[5i+1:5i+5]) );    push!(coreSubshells, sh)
            if  5i + 6 >  length(sa)    break    end
        end
        sa = readline(f);   println("sa = $sa")
        sa = readline(f);   if  length(sa) > 3    go = true    else  go = false    end
        i  = -1
        while go
            i = i + 1;    sh = Basics.subshellGrasp( strip(sa[5i+1:5i+5]) );    push!(peelSubshells, sh)
            if  5i + 6 >  length(sa)    break    end
        end
        # Prepare a list of all subshells to define the common basis
        subshells = deepcopy(coreSubshells);    append!(subshells, peelSubshells)

        # Now, read the CSF in turn until the EOF
        readline(f)
        NoCSF = 0;    csfs = CsfR[]
        Defaults.setDefaults("relativistic subshell list", subshells)

        while true
            sa = readline(f);      if length(sa) == 0   break    end
            sb = readline(f);      sc = readline(f)
            NoCSF = NoCSF + 1
            push!(csfs, ManyElectron.CsfRGrasp92(subshells, coreSubshells, sa, sb, sc) )
        end

        println("  ... $NoCSF CSF read in from Grasp92 file  $filename)")
        NoElectrons = sum( csfs[1].occupation )

        Basis(true, NoElectrons, subshells, csfs, coreSubshells, Dict{Subshell,Radial.Orbital}() )
    end

    """
    `Basics.readOrbitalFileGrasp92(filename::String, grid.Radial.Grid)`
        ... reads in the orbitals list from a Grasp92 .rwf or .out file; a dictionary
            orbitals::Dict{Subshell,Radial.Orbital} is returned.
    """
    function Basics.readOrbitalFileGrasp92(filename::String, grid::Radial.Grid)
        # using FortranFiles
        f = FortranFile(filename)

        orbitals = Dict{Subshell,Radial.Orbital}();    first = true

        if (String(Base.read(f, FString{6})) != "G92RWF")
        close(f)
        error("File \"", filename, "\" is not a proper Grasp92 RWF file")
        end

        Defaults.setDefaults("standard grid", grid);

        while true
        try
            n, kappa, energy, mtp = Base.read(f, Int32, Int32, Float64, Int32)
            pz, P, Q = Base.read(f, Float64, (Float64, mtp), (Float64, mtp))
            ra = Base.read(f, (Float64, mtp))
                # Now place the data into the right fields of Orbital
                subshell = Subshell( Int64(n), Int64(kappa) )
                if  energy < 0   isBound = true    else    isBound = false   end
                useStandardGrid = true
                println("Basics.readOrbitalFileGrasp92-aa: warning ... Radial functions still need to be interpolated upon the given grid.")
                Px = P
                Qx = Q

                orbitals = Base.merge( orbitals, Dict( subshell => Orbital(subshell, isBound, useStandardGrid, energy, Px, Qx, Radial.Grid() ) ))
        catch ex
            if     ex isa EOFError  break
            else   throw(ex)
            end
        end
        end
        close(f)

        return( orbitals )
    end

    """
    `Basics.readMixFileRelci(filename::String)`
        ... reads in the mixing coefficients from a RELCI .mix file; a multiplet::Multiplet is returned.
    """
    function Basics.readMixFileRelci(filename::String, basis::Basis)
        levels = Level[];    name =  "from Relci"

        f  = open(filename)
        sa = readline(f);   sa[1:14] != "G92MIX (format"   &&   error("Not a formatted relci.mix file.")
        sa = readline(f);   NoElectrons = Base.parse(Int64, strip(sa[1:6]) );    NoCsf = Base.parse(Int64, strip(sa[7:12]) )
                            NoSubshells = Base.parse(Int64, strip(sa[13:18]) )
        sa = readline(f);   NoLevels    = Base.parse(Int64, strip(sa[1:6]) )
        #
        sa = readline(f)
        sa = readline(f);   list2J = AngularJ64[AngularJ64(0)  for  i = 1:NoLevels];    listp = Basics.Parity[Basics.plus  for  i = 1:NoLevels]
        for  i = 1:NoLevels
            jj = Base.parse(Int64, strip(sa[(i-1)*6+1:(i-1)*6+5]) );   list2J[i] = AngularJ64(jj//2)
            pp = string( sa[(i-1)*6+6] );                              listp[i]  = Basics.Parity( pp )
        end
        #
        sa   = readline(f);   listEnergy = Float64[0.  for  i = 1:NoLevels]
        enav = Base.parse(Float64, strip(sa[1:26]) )
        for  i = 1:NoLevels
            en = Base.parse(Float64, strip(sa[26i+1:26i+26]) );   listEnergy[i] = en
        end
        #
        eigenvectors = zeros(NoCsf,NoLevels)
        for  j = 1:NoCsf
            sa   = readline(f)
            for  i = 1:NoLevels
            cm = Base.parse(Float64, strip(sa[16*(i-1)+1:16(i-1)+16]) );   eigenvectors[j,i] = cm
            end
        end
        #
        for  j = 1:NoLevels
            mc = Float64[0.  for  i = 1:NoCsf]
            for  i = 1:NoCsf    mc[i] = eigenvectors[i,j]    end
            Jx = list2J[j];   Mx = AngularM64( Jx.num//Jx.den )
            push!( levels, Level(Jx, Mx, listp[j], j, enav+listEnergy[j], 0., true, basis, mc) )
        end

        multiplet = Multiplet(name, levels)
        return( multiplet )
    end

    """
    `Basics.readMixingFileGrasp92(filename::String, basis::Basis)`
        ... reads in the mixing coefficients from a Grasp92 .mix file by using the given basis; A multiplet::Multiplet
            is returned.
    """
    function Basics.readMixingFileGrasp92(filename::String, basis::Basis)
        name = "from Grasp"
        f = FortranFile(filename)

        if (String(Base.read(f, FString{6})) != "G92MIX")
            close(f)
        error("File \"", filename, "\" does not seem to be a Grasp92 MIX file")
        end

        no_electrons, no_csf, no_subshells = read(f, Int32, Int32, Int32)
        no_asf = Base.read(f, Int32)
        level_numbers = Base.read(f, (Int32, no_asf))
        ijip = Base.read(f, (Int32, 2, no_asf))

        mjs = ijip[1,:] # 2j+1
        mparities = ijip[2,:] # 1=='+', -1=='-'

        mixfile_to_j(i::Int32) = J((i-1)//2)
        mixfile_to_p(i::Int32) = i==1 ?  Basics.plus : Basics.minus

        js = mixfile_to_j.(mjs)
        parities = mixfile_to_p.(mparities)

        avg_energy, eval = read(f, Float64, (Float64, no_asf))
        energies = avg_energy + eval
        eigenvectors = read(f, (Float64, no_csf, no_asf))

        close(f)
        multiplet = Multiplet(name, Level[] )
        return( multiplet )
    end

    """
    `Basics.recast()`
        ... recasts some data from one number/representation into another one; cf. the supported keystrings and return values.

    + `("rate: radiative, to decay width", value::Float64)`
        ... to recast a given radiative rate (Einstein A in atomic units) into a decay withs, taking the selected energy unit
            into account. A Float64 is returned.

    + `("rate: radiative, to Einstein A", value::Float64)`
        ... to recast a given spontaneous radiative rate (= Einstein A-coefficient), taking the selected unit into account.
            A Float64 is returned.

    + `("rate: radiative, to Einstein B", value::Float64)`
        ... to recast a given radiative rate (Einstein A in atomic units) into an Einstein B-coefficient, taking the selected
            unit into account. A Float64 is returned.

    + `("rate: radiative, to g_f", value::Float64)`
        ... to recast a given radiative rate (Einstein A in atomic units) into an oscillator strength g_f;
            a Float64 is returned.
    """
    function Basics.recast(sa::String, line::Union{Einstein.Line,PhotoEmission.Line}, wa::Float64)

        if       sa == "rate: radiative, to decay width"
            width = Defaults.convertUnits("energy: from atomic", wa)
            return( width )

        elseif   sa == "rate: radiative, to Einstein A"
            einsteinA = Defaults.convertUnits("rate: from atomic", wa)
            return( einsteinA )

        elseif   sa == "rate: radiative, to Einstein B"
            einsteinB = (Basics.twice(line.initialLevel.J) + 1) / (Basics.twice(line.finalLevel.J) + 1) * pi^2 *
                        Defaults.getDefaults("speed of light: c")^3 / line.omega^3  * wa
            einsteinB = Defaults.convertUnits("Einstein B: from atomic", einsteinB)
            return( einsteinB )

        elseif   sa == "rate: radiative, to g_f"
            gf = (Basics.twice(line.initialLevel.J) + 1) / (Basics.twice(line.finalLevel.J) + 1) / 2. *
                Defaults.getDefaults("speed of light: c")^3 / line.omega^2 * wa
            return( gf )

        else     error("Unsupported keystring = $sa")
        end
    end

    """
    `Basics.selectLevel(level::Level, levelSelection::LevelSelection)`
        ... returns true::Bool if the levelSelection is inactive or if the level has been selected due to its
            indices or symmetries; in all other case, false::Bool is returned.
    """
    function Basics.selectLevel(level::Level, levelSelection::LevelSelection)
        if  levelSelection.active
            # Test for level index
            ##x @show level.index, levelSelection.indices
            if  level.index  in  levelSelection.indices                                 return( true )   end
            if  LevelSymmetry(level.J, level.parity)  in  levelSelection.symmetries     return( true )   end
        else                                                                            return( true )
        end
        return( false )
    end

    """
    `Basics.selectLevelPair(iLevel::Level, fLevel::Level, lineSelection::LineSelection)`
        ... returns true::Bool if the lineSelection is inactive or if the pair (iLevel, fLevel) has been selected due to its
            indices or symmetries; in all other case, false::Bool is returned.
    """
    function Basics.selectLevelPair(iLevel::Level, fLevel::Level, lineSelection::LineSelection)
        if     lineSelection.active
            # Test for level indexPairs
            for ip in  lineSelection.indexPairs
                if      ip[1] == 0  &&  ip[2] == fLevel.index    return( true )
                elseif  ip[2] == 0  &&  ip[1] == iLevel.index    return( true )
                elseif  ip == (iLevel.index, fLevel.index)       return( true )
                end
            end
            # Test for level symmetries
            for sp in  lineSelection.symmetryPairs
                if      ip == (LevelSymmetry(iLevel.J, iLevel.parity),  LevelSymmetry(fLevel.J, fLevel.parity))   return( true )
                end
            end
        else                                                     return( true )
        end
        return( false )
    end

    """
    `Basics.selectLevelTriple(iLevel::Level, nLevel::Level, fLevel::Level, pathwaySelection::PathwaySelection)`
        ... returns true::Bool if the pathwaySelection is inactive or if the triple (iLevel, nLevel, fLevel) has been selected
            due to its indices or symmetries; in all other case, false::Bool is returned.
    """
    function Basics.selectLevelTriple(iLevel::Level, nLevel::Level, fLevel::Level, pathwaySelection::PathwaySelection)
        if     pathwaySelection.active
            # Test for level indexTriples
            for ip in  pathwaySelection.indexTriples
                if      ip[1] == 0  &&  ip[2] == 0             &&  ip[3] == 0               return( true )
                elseif  ip[1] == 0  &&  ip[2] == 0             &&  ip[3] == fLevel.index    return( true )
                elseif  ip[1] == 0  &&  ip[2] == nLevel.index  &&  ip[3] == 0               return( true )
                elseif  ip[2] == 0  &&  ip[3] == 0             &&  ip[1] == iLevel.index    return( true )
                elseif  ip[1] == 0  &&  ip[2] == nLevel.index  &&  ip[3] == fLevel.index    return( true )
                elseif  ip[2] == 0  &&  ip[1] == iLevel.index  &&  ip[3] == fLevel.index    return( true )
                elseif  ip[3] == 0  &&  ip[1] == iLevel.index  &&  ip[2] == nLevel.index    return( true )
                elseif  ip == (iLevel.index, nLevel.index, fLevel.index)                    return( true )
                end
            end
            # Test for level symmetries
            for sp in  pathwaySelection.symmetryTriples
                if      ip == (LevelSymmetry(iLevel.J, iLevel.parity),  LevelSymmetry(nLevel.J, nLevel.parity),
                               LevelSymmetry(fLevel.J, fLevel.parity))                      return( true )
                end
            end
        else                                                                                return( true )
        end
        return( false )
    end

    """
    `Basics.selectSymmetry(sym::LevelSymmetry, levelSelection::LevelSelection)`
        ... returns true::Bool if the levelSelection is inactive or if the symmetry sym has been selected;
            in all other case, false::Bool is returned.
    """
    function Basics.selectSymmetry(sym::LevelSymmetry, levelSelection::LevelSelection)
        if  levelSelection.active
            if      length(levelSelection.symmetries) == 0    return( true )
            elseif  sym in levelSelection.symmetries          return( true )
            end
        else                                                  return( true )
        end
        return( false )
    end

    """
    `Basics.shiftTotalEnergies(multiplet::Multiplet, energyShift::Float64)`
        ... to shift the energies of all levels in the multiplet by energyShift [a.u.]; a (new) multiplet::Multiplet is returned.
    """
    function Basics.shiftTotalEnergies(multiplet::Multiplet, energyShift::Float64)
        newLevels = Level[]
        for lev in multiplet.levels
            push!(newLevels, Level(lev.J, lev.M, lev.parity, lev.index, lev.energy + energyShift, lev.relativeOcc, lev.hasStateRep, lev.basis, lev.mc) )
        end

        newMultiplet = Multiplet(multiplet.name, newLevels)

        return( newMultiplet )
    end

    """
    `Basics.sortByEnergy(multiplet::Multiplet)`
        ... to sort all levels in the multiplet into a sequence of increasing energy; a (new) multiplet::Multiplet is returned.
    """
    function Basics.sortByEnergy(multiplet::Multiplet)
        sortedLevels = Base.sort( multiplet.levels , lt=Base.isless)
        newLevels = Level[];   index = 0
        for lev in sortedLevels
            index = index + 1
            push!(newLevels, Level(lev.J, lev.M, lev.parity, index, lev.energy, lev.relativeOcc, lev.hasStateRep, lev.basis, lev.mc) )
        end

        newMultiplet = Multiplet(multiplet.name, newLevels)

        return( newMultiplet )
    end

    """
    `Basics.sortByEnergy(multiplet::Hfs.IJF_Multiplet)`
        ... to sort all hyperfine levels in the multiplet into a sequence of increasing energy; a (new) multiplet::Hfs.IJF_Multiplet
            is returned.
    """
    function Basics.sortByEnergy(multiplet::Hfs.IJF_Multiplet)
        sortedLevels = Base.sort( multiplet.levelFs , lt=Base.isless)
        newLevels = Hfs.IJF_Level[];   index = 0
        for lev in sortedLevels
            index = index + 1
            push!(newLevels, Hfs.IJF_Level(lev.I, lev.F, lev.M, lev.parity, lev.energy, lev.hasStateRep, lev.basis, lev.mc) )
        end

        newMultiplet = Hfs.IJF_Multiplet(multiplet.name, newLevels)

        return( newMultiplet )
    end

    """
    `Basics.tabulate(stream::IO, sa::String, multiplet::Multiplet, levelNos::Array{Int64,1})`
        ... tabulates the energies from the multiplet with level numbers in levelNos due to different criteria.

    + `(stream::IO, "multiplet: energies", multiplet::Multiplet)`
        ... to tabulate the energies of all levels of the given multiplet into a neat format; nothing is returned.
    + `(stream::IO, "multiplet: energy relative to immediately lower level", multiplet::Multiplet, levelNos::Array{Int64,1})`
        ... to tabulate the energy splitting between neighboured levels of all levels of the given multiplet into a neat
            format; nothing is returned.
    + `(stream::IO, "multiplet: energy of each level relative to lowest level", multiplet::Multiplet, levelNos::Array{Int64,1}`
        ... to tabulate the energy splitting of all levels with regard to the lowest level of the given multiplet into
            a neat format; nothing is returned.
    """
    function Basics.tabulate(stream::IO, sa::String, multiplet::Multiplet, levelNos::Array{Int64,1})
        if        sa == "multiplet: energies"
            println(stream, "\n  Eigenenergies:")
            sb = "  Level  J Parity          Hartrees       " * "             eV                   " *  TableStrings.inUnits("energy")
            println(stream, "\n", sb, "\n")
            for  i = 1:length(multiplet.levels)
                if  i in levelNos
                    lev = multiplet.levels[i]
                    en  = lev.energy;    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
                    sc  = " " * TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
                    @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
                end
            end

        elseif    sa == "multiplet: energy relative to immediately lower level"
            println(stream, "\n  Energy of each level relative to immediately lower level:")
            sb = "  Level  J Parity          Hartrees       " * "             eV                   " * TableStrings.inUnits("energy")
            println(stream, "\n", sb, "\n")
            for  i = 2:length(multiplet.levels)
                if  i in levelNos
                    lev = multiplet.levels[i]
                    en    = lev.energy - multiplet.levels[i-1].energy;
                    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
                    sc  = " " * TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
                    @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
                end
            end

        elseif    sa == "multiplet: energy of each level relative to lowest level"
            println(stream, "\n  Energy of each level relative to lowest level:")
            sb = "  Level  J Parity          Hartrees       " * "             eV                   " * TableStrings.inUnits("energy")
            println(stream, "\n", sb, "\n")
            for  i = 2:length(multiplet.levels)
                if  i in levelNos
                    lev = multiplet.levels[i]
                    en    = lev.energy - multiplet.levels[1].energy;
                    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
                    sc  = " " * TableStrings.level(i) * "    " * string(LevelSymmetry(lev.J, lev.parity)) * "    "
                    @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
                end
            end
        else
            error("Unsupported keystring.")
        end

        return( nothing )
    end

    """
    `Basics.tabulate(sa::String, multiplet::Hfs.IJF_Multiplet; stream::IO=stdout)`
        ... tabulates the energies from the multiplet due to different criteria.

    + `("multiplet: energies", multiplet::Hfs.IJF_Multiplet; stream::IO=stdout)`
        ... to tabulate the energies of all hyperfine levels of the given multiplet into a neat format; nothing is returned.
    + `("multiplet: energy of each level relative to lowest level", multiplet::Hfs.IJF_Multiplet; stream::IO=stdout)`
        ... to tabulate the energy splitting of all levels with regard to the lowest level of the given multiplet into
            a neat format; nothing is returned.
    """
    function Basics.tabulate(sa::String, multiplet::Hfs.IJF_Multiplet; stream::IO=stdout)
        if        sa == "multiplet: energies"
            println(stream, "\n  Eigenenergies for nuclear spin I = $(multiplet.levelFs[1].I):")
            sb = "  Level  F Parity          Hartrees       " * "             eV                   " *  TableStrings.inUnits("energy")
            println(stream, "\n", sb, "\n")
            for  i = 1:length(multiplet.levelFs)
                lev = multiplet.levelFs[i]
                en  = lev.energy;    en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
                sc  = " " * TableStrings.level(i) * "    " * string(LevelSymmetry(lev.F, lev.parity)) * "    "
                @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
            end

        elseif    sa == "multiplet: energy of each level relative to lowest level"
            println(stream, "\n  Energy of each level relative to lowest level for nuclear spin I = $(multiplet.levelFs[1].I):")
            sb = "  Level  F Parity          Hartrees       " * "             eV                   " * TableStrings.inUnits("energy")
            println(stream, "\n", sb, "\n")
            for  i = 2:length(multiplet.levelFs)
                lev = multiplet.levelFs[i]
                en    = lev.energy - multiplet.levelFs[1].energy;
                en_eV = Defaults.convertUnits("energy: from atomic to eV", en);    en_requested = Defaults.convertUnits("energy: from atomic", en)
                sc  = " " * TableStrings.level(i) * "    " * string(LevelSymmetry(lev.F, lev.parity))  * "    "
                @printf(stream, "%s %.15e %s %.15e %s %.15e %s", sc, en, "  ", en_eV, "  ", en_requested, "\n")
            end
        else
            error("Unsupported keystring.")
        end

        return( nothing )
    end

    """
    `Basics.tabulateKappaSymmetryEnergiesDirac(kappa::Int64, evalues::Array{Float64,1}, ns::Int64, nuclearModel::Nuclear.Model)`
        ... tabulates the eigenenergies for a given symmetry block kappa together with the corresponding Dirac energies for a
            point-like nucleus. The index ns tells the number of 'negative-continnum' energies in the given evalues. nothing is
            returned.
    """
    function Basics.tabulateKappaSymmetryEnergiesDirac(kappa::Int64, evalues::Array{Float64,1}, ns::Int64, nuclearModel::Nuclear.Model)
        Z = nuclearModel.Z;    nx = 77
        # Determine the allowed principal quantum numbers n
        l = Basics.subshell_l(Subshell(101,kappa))
        println("  ", TableStrings.hLine(nx))
        sa = "  "
        sa = sa * TableStrings.center( 7, "Index";            na=2)
        sa = sa * TableStrings.center(10, "Subshell";         na=3)
        sa = sa * TableStrings.center(17, "Energies [a.u.]";  na=2)
        sa = sa * TableStrings.center(17, "Dirac-E  [a.u.]";  na=2)
        sa = sa * TableStrings.center(17, "Delta-E / |E|";    na=2)
        println(sa)
        println("  ", TableStrings.hLine(nx))
        for  i = ns+1:ns+38
            sa = " " * TableStrings.center( 6, TableStrings.level(i-ns); na=2)
            sa = sa *  TableStrings.flushright(10, string(Subshell(i-ns+l, kappa)); na=6)
            en = Basics.computeDiracEnergy(Subshell(i-ns+l, kappa), Z)
            if  evalues[i] >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", evalues[i])                       * "    "
            if  en         >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", en)                               * "    "
            if  evalues[i]-en >= 0.   sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", (evalues[i]-en)/abs(evalues[i])) * "    "
            println(sa)
        end
        println("      :       :    ")
        for  i = length(evalues)-5:length(evalues)
            sa = " " * TableStrings.center( 6, TableStrings.level(i-ns); na=2)
            sa = sa *  TableStrings.flushright(10, string(Subshell(i-ns+l, kappa)); na=6)
            en = Basics.computeDiracEnergy(Subshell(i-ns+l, kappa), Z)
            if  evalues[i] >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", evalues[i])                       * "    "
            if  en         >= 0.      sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", en)                               * "    "
            if  evalues[i]-en >= 0.   sa = sa * "+"   end;    sa = sa * @sprintf("%.8e", (evalues[i]-en)/abs(evalues[i])) * "    "
            println(sa)
        end
        println("  ", TableStrings.hLine(nx))

        return( nothing )
    end

    """
    `Basics.tools(dict::Dict)`
        ... select different tools from a menu if proper results::dict are given in terms of a dictionary as obtained
            usually from an Atomic.Computation.
    """
    function Basics.tools(dict::Dict)
        println("Tools.select(): A menu if a dict comes in; $dict ")
    end

    """
    `Basics.yesno(question::String, sa::String)`
        ... Returns true if the answer 'yes' or 'y' is found, and false otherwise; sa = {"Y", "N"} determines how the zero-String
            "" is interpreted. The given question is repeated until a proper answer is obtained.
    """
    function yesno(question::String, sa::String)
        while  true
            print(question * "  ");    reply = strip( readline(STDIN) )
            if      sa == "Y"
                if      reply in ["", "Y", "Yes", "y", "yes"]   return( true )
                elseif  reply in [    "N", "No",  "n", "no" ]   return( false )
                        println("... answer not recognized ... redo:")
                end
            elseif  sa == "N"
                if      reply in [    "Y", "Yes", "y", "yes"]   return( true )
                elseif  reply in ["", "N", "No",  "n", "no" ]   return( false )
                        println("... answer not recognized ... redo:")
                end
            end
        end
    end


text:
export get_lbi_timescale, get_lbi!

"""
	get_lbi_timescale(lineage, segment)

Return the timescale used for computing LBI values, in branch length units.
 """
function get_lbi_timescale(lineage, segment)
	if typeof(lbi_integration_time) == Year
		ct = lbi_integration_time.value
	elseif typeof(lbi_integration_time) == Month
		ct = lbi_integration_time.value / 12
	elseif typeof(lbi_integration_time) == Day
		ct = lbi_integration_time.value / 365
	else
		ct = lbi_integration_time
	end
	τ = substitution_rate[lineage, segment] * ct
	return τ
end

"""
	get_lbi!(t::Tree{LBIData}, fp::FluPop, strainnames::Array{String}, datemin, datemax, τ, lbi_field)

Compute lbi of `strainsnames` according to tree `t`. Leaves of tree posterior to `datemax` or anterior to `datemin` are set to dead. Tree is modified in the process (life-state and lbi value of nodes).
Store the resulting lbi value in `fp.strains[lbi_field]`
"""
function get_lbi!(t::Tree{LBIData}, fp::FluPop, strainnames::Array{<:AbstractString}, datemin, datemax, τ, lbi_field, verbose=false)
	# Setting life-state of tree nodes
	for node in values(t.leaves)
		if haskey(fp.strains, node.label)
			node.data.alive = datemin <= fp.strains[node.label][:date] < datemax
		else
			node.data.alive = false
		end
	end
	TreeTools.set_live_nodes!(t, set_leaves=false)
	# Computing LBI
	TreeTools.lbi!(t, τ, normalize=true) && verbose && begin
		println("LBI computation failed")
	    println("-- ", length(strainnames), " strains")
	    println("-- Date range $datemin - $datemax")
	end
	for (i,s) in enumerate(strainnames)
		if haskey(t.lleaves, s)
			fp.strains[s].data[lbi_field] = t.lleaves[s].data.lbi
		else
			fp.strains[s].data[lbi_field] = missing
		end
	end
	nothing
end

"""
	function get_lbi!(fp::FluPop, t::Tree{LBIData}, datestyle=:all_anterior;
				lineage="h3n2",
				segment="ha",
				τ=get_tau(lineage,segment),
				lbi_field=:lbi)

Compute lbi of strains in `fp` based on tree `t` and date binning of `fp.datebin`. Tree is modified in the process (life-state of nodes).
The `datestyle` arguments modifies the life-state of nodes when LBI is computed. If `:all_anterior`, all nodes anterior to the current datebin are set to live. If `:datebin`, only nodes in the current datebin are set to live.
If a strain does not exist in the tree, its lbi will be set to `missing`.

"""
function get_lbi!(fp::FluPop, t::Tree{LBIData}, datestyle=:all_anterior;
				lineage="h3n2",
				segment="ha",
				τ=get_lbi_timescale(lineage,segment),
				lbi_field=:lbi,
				verbose=false)
	# Find `datemin` as a function of `datestyle`
	if datestyle == :all_anterior
		datemin = let val = findmin(collect(keys(fp.datebin)))[1][1]
			(x) -> val
		end
	elseif datestyle == :datebin
		datemin = (x) -> minimum(x)
	else
		error("Accepted values for `datestyle`: `:all_anterior` or `:datebin`")
	end
	# Compute LBI for all datebins
	for (datebin,strains) in sort(OrderedDict(fp.datebin), by=x->x[1])
		if length(strains)>1
			labels = [x[:strain] for x in strains]
			get_lbi!(t, fp, labels, datemin(datebin), max(datebin...), τ, lbi_field, verbose)
		else
			for s in strains
				s.data[lbi_field] = missing
			end
		end
	end
	# Setting lbi to missing for strains not in a datebin
	for s in values(fp.strains)
		if !haskey(s.data, lbi_field)
			s.data[lbi_field] = missing
		end
	end
end

"""
	local_lbi_maximas(t::Tree{LBIData}, fp::FluPop, τ)

The idea is to find local maximas of the LBI for each date bin of `fp`. In order to determine if `n` is a local maximum, one needs to know the LBI of its ancestor `a`. Thus, the lbi of `a` has to be recomputed for every datebin, when new leaves come live and old ones die.
## Notes
- `fp` should be date-binned already.
- Strains in `fp` and `t` should match **exactly**.
"""
function local_lbi_maximas!(t::Tree{LBIData}, fp::FluPop, τ=get_lbi_timescale("h3n2","ha"))
	out = Dict()
	# Setting nodes alive in date order
	for (db, strains) in sort(OrderedDict(fp.datebin), by=x->x[1])
		for n in values(t.lnodes)
			n.data.alive = false
		end
		for s in strains
			t.lleaves[s[:strain]].data.alive = true
		end
		TreeTools.set_live_nodes!(t, set_leaves=false)
		TreeTools.lbi!(t, τ, normalize=true)
		# Local maximas for the current tree
		llm = local_lbi_maximas(t)
		out[db] = [(label=llm[i], lbi=t.lnodes[llm[i]].data.lbi) for i in 1:length(llm)]
	end
	return out
end
"""
	local_lbi_maximas(t::Tree{LBIData})

Find all nodes that are a local maxima of the LBI. Does not compute the LBI.
"""
function local_lbi_maximas(t::Tree{LBIData})
	out = String[]
	for (l,n) in t.lnodes
		if n.isroot || n.data.lbi > n.anc.data.lbi
			if mapreduce(x-> n.data.lbi > x.data.lbi, *, n.child, init=true)
				push!(out, n.label)
			end
		end
	end
	return out
end

# if false

# """
# 	local_lbi_maximas(t::Tree{LBIData}, sp::StrainPop, τ)

# The idea is to find local maximas of the LBI for each date bin of `sp`. In order to determine if `n` is a local maximum, one needs to know the LBI of its ancestor `a`. Thus, the lbi of `a` has to be recomputed for every datebin, when new leaves come live and old ones die.
# ## Notes
# - `sp` should be date-binned already.
# - Strains in `sp` and `t` should match **exactly**.
# """
# function local_lbi_maximas!(t::Tree{LBIData}, sp::StrainPop, τ)
# 	out = Dict()
# 	# Setting nodes alive in date order
# 	for (db, strains) in sort(sp.datebin, by=x->x[1])
# 		for n in values(t.lnodes)
# 			n.data.alive = false
# 		end
# 		for s in strains
# 			t.lleaves[s.strain].data.alive = true
# 		end
# 		TreeTools.set_live_nodes!(t, set_leaves=false)
# 		TreeTools.lbi!(t, τ, normalize=true)
# 		# Local maximas for the current tree
# 		llm = local_lbi_maximas(t)
# 		out[db] = [(label=llm[i], lbi=t.lnodes[llm[i]].data.lbi) for i in 1:length(llm)]
# 	end
# 	return out
# end
# """
# 	local_lbi_maximas(t::Tree{LBIData})

# Find all nodes that are a local maxima of the LBI. Does not recompute LBI.
# """
# function local_lbi_maximas(t::Tree{LBIData})
# 	out = String[]
# 	for (l,n) in t.lnodes
# 		if n.isroot || n.data.lbi > n.anc.data.lbi
# 			if mapreduce(x-> n.data.lbi > x.data.lbi, *, n.child, init=true)
# 				push!(out, n.label)
# 			end
# 		end
# 	end
# 	return out
# end

# end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
#=
Popis

Parametre

Navratova hodnota

Poznamky

=#

function rescale_graph(betweenness, n, normalized, k)

  if normalized == true

    if n <= 2
      proportion = nothing
    else
      proportion = 1.0 / ((n - 1) * (n - 2))
    end

  else
    proportion = 1.0 / 2.0
  end

  if proportion != nothing

    if k != nothing
      proportion = proportion * n / k
    end

        for v in keys(betweenness)
            betweenness[v] *= proportion
        end

  end

  return betweenness

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
"""
$(TYPEDEF)

A special basis over the states with parameters , independent variable  and possible exogenous controls.
It extends an `AbstractBasis`, which also stores information about the lifted dynamics, specified by a sufficient
matrix factorization, an output mapping and internal variables to update the equations. It can be called with the typical SciML signature, meaning out of place with `f(u,p,t)`
or in place with `f(du, u, p, t)`. If control inputs are present, it is assumed that no control corresponds to
zero for all inputs. The corresponding function calls are `f(u,p,t,inputs)` and `f(du,u,p,t,inputs)` and need to
be specified fully.

If `linear_independent` is set to `true`, a linear independent basis is created from all atom functions in `f`.

If `simplify_eqs` is set to `true`, `simplify` is called on `f`.

Additional keyworded arguments include `name`, which can be used to name the basis, and
`observed` for defining observeables.

# Fields
$(FIELDS)

## Note

The keyword argument `eval_expression` controls the function creation
behavior. `eval_expression=true` means that `eval` is used, so normal
world-age behavior applies (i.e. the functions cannot be called from
the function that generates them). If `eval_expression=false`,
then construction via GeneralizedGenerated.jl is utilized to allow for
same world-age evaluation. However, this can cause Julia to segfault
on sufficiently large basis functions. By default eval_expression=false.
"""
mutable struct Koopman{O,M,G,T} <: AbstractKoopman
    """The equations of the basis"""
    eqs::Vector{Equation}
    """Dependent (state) variables"""
    states::Vector
    """Control variables"""
    ctrls::Vector
    """Parameters"""
    ps::Vector
    """Observed"""
    observed::Vector
    """Independent variable"""
    iv::Num
    """Internal function representation of the basis"""
    f::Function
    """Associated lifting of the operator"""
    lift::Function
    """Name of the basis"""
    name::Symbol
    """Internal systems"""
    systems::Vector{Basis}
    """Discrete or time continuous"""
    is_discrete::Bool
    """The operator/generator of the dynamics"""
    K::O
    """Mapping back onto the observed states"""
    C::M
    """Internal matrix `Q` used for updating"""
    Q::G
    """Internal matrix `P` used for updating"""
    P::T
end

function Koopman(eqs::AbstractVector{Equation}, states::AbstractVector;
    K::O = diagm(ones(eltype(states), length(eqs))),
    C::AbstractMatrix = diagm(ones(eltype(K), length(eqs))),
    Q::AbstractMatrix = zeros(eltype(states), 0,0),
    P::AbstractMatrix = zeros(eltype(states), 0,0),
    lift::Function = (X, args...)->identity(X),
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], observed::AbstractVector = [],
    name = gensym(:Koopman), is_discrete::Bool = true,
    digits::Int = 10, # Round all coefficients
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...) where O <: Union{AbstractMatrix, Eigen, Factorization}

    iv === nothing && (iv = Variable(:t))
    iv = value(iv)
    eqs = scalarize(eqs)
    states, controls, parameters, observed = value.(scalarize(states)), value.(scalarize(controls)), value.(scalarize(parameters)), value.(scalarize(observed))

    eqs = [eq for eq in eqs if ~isequal(Num(eq),zero(Num))]

    lhs = Num[x.lhs for x in eqs]
    eqs_ = Num[x.rhs for x in eqs]

    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs_, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs_) : eqs_
    end

    unique!(eqs_, !simplify)

    f = DataDrivenDiffEq._build_ddd_function(eqs_, states, parameters, iv, controls, eval_expression)

    eqs = [lhs[i] ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Koopman{typeof(K), typeof(C), typeof(Q), typeof(P)}(eqs,
    states, controls, parameters, observed , iv, f, lift, name, Basis[],
    is_discrete, K, C, Q, P)
end

function Koopman(eqs::AbstractVector{Num}, states::AbstractVector;
    K::O = diagm(ones(eltype(states), length(eqs))),
    C::AbstractMatrix = diagm(ones(eltype(K), length(eqs))),
    Q::AbstractMatrix = zeros(eltype(states), 0,0),
    P::AbstractMatrix = zeros(eltype(states), 0,0),
    lift::Function = (X, args...)->identity(X),
    parameters::AbstractVector = [], iv = nothing,
    controls::AbstractVector = [], observed::AbstractVector = [],
    name = gensym(:Koopman), is_discrete::Bool = true,
    digits::Int = 10, # Round all coefficients
    simplify = false, linear_independent = false,
    eval_expression = false,
    kwargs...) where O <: Union{AbstractMatrix, Eigen, Factorization}

    iv === nothing && (iv = Variable(:t))
    iv = value(iv)
    eqs = scalarize(eqs)
    states, controls, parameters, observed = value.(scalarize(states)), value.(scalarize(controls)), value.(scalarize(parameters)), value.(scalarize(observed))

    eqs_ = [eq for eq in eqs if ~isequal(Num(eq),zero(Num))]

    if linear_independent
        eqs_ = create_linear_independent_eqs(eqs_, simplify)
    else
        eqs_ = simplify ? ModelingToolkit.simplify.(eqs_) : eqs_
    end

    unique!(eqs_, !simplify)

    f = DataDrivenDiffEq._build_ddd_function(eqs_, states, parameters, iv, controls, eval_expression)
    D = Differential(iv)

    eqs = [D(states[i]) ~ eq for (i,eq) ∈ enumerate(eqs_)]

    return Koopman{typeof(K), typeof(C), typeof(Q), typeof(P)}(eqs,
    states, controls, parameters, observed, iv, f, lift, name, Basis[],
    is_discrete, K, C, Q, P)
end

# We assume that we only have real valued observed
Base.Matrix(k::AbstractKoopman) = real.(Matrix(_get_K(k)))

# Get the lifting function
lifting(k::AbstractKoopman) = getfield(k, :lift)

# Get K
_get_K(k::AbstractKoopman) = getfield(k, :K)

implicit_variables(k::AbstractKoopman) = []

"""
$(SIGNATURES)

Returns `true` if the `AbstractKoopmanOperator` `k` is discrete in time.
"""
is_discrete(k::AbstractKoopman) = getfield(k, :is_discrete)

"""
$(SIGNATURES)

Returns `true` if the `AbstractKoopmanOperator` `k` is continuous in time.
"""
is_continuous(k::AbstractKoopman) = !is_discrete(k)

"""
$(SIGNATURES)

Return the eigendecomposition of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigen(k::AbstractKoopman) = begin
    K = _get_K(k)
    isa(K, Eigen) && return K
    eigen(K)
end

"""
$(SIGNATURES)

Return the eigenvalues of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigvals(k::AbstractKoopman) = eigvals(_get_K(k))

"""
$(SIGNATURES)

Return the eigenvectors of the `AbstractKoopmanOperator`.
"""
LinearAlgebra.eigvecs(k::AbstractKoopman) = eigvecs(_get_K(k))

"""
$(SIGNATURES)

Return the eigenvectors of a continuous `AbstractKoopmanOperator`.
"""
modes(k::AbstractKoopman) = is_continuous(k) ? eigvecs(k) : throw(AssertionError("Koopman is discrete."))

"""
$(SIGNATURES)

Return the eigenvalues of a continuous `AbstractKoopmanOperator`.
"""
frequencies(k::AbstractKoopman) = is_continuous(k) ? eigvals(k) : throw(AssertionError("Koopman is discrete."))

"""
$(SIGNATURES)

Return the approximation of the discrete Koopman operator stored in `k`.
"""
operator(k::AbstractKoopman) = is_discrete(k) ? _get_K(k) : throw(AssertionError("Koopman is continouos."))

"""
$(SIGNATURES)

Return the approximation of the continuous Koopman generator stored in `k`.
"""
generator(k::AbstractKoopman) = is_continuous(k) ? _get_K(k) : throw(AssertionError("Koopman is discrete."))

"""
$(SIGNATURES)

Return the array `C`, mapping the Koopman space back onto the state space.
"""
outputmap(k::AbstractKoopman) = Symbolics.unwrap(k.C)

"""
$(SIGNATURES)

Returns `true` if the `AbstractKoopmanOperator` is updatable.
"""
updatable(k::AbstractKoopman) = !isempty(k.Q) && !isempty(k.P)

"""
$(SIGNATURES)

Returns `true` if either:

+ the Koopman operator has just eigenvalues with magnitude less than one or
+ the Koopman generator has just eigenvalues with a negative real part
"""
is_stable(k::AbstractKoopman) = begin
    K = _get_K(k)
    is_discrete(k) && all(real.(eigvals(k)) .< real.(one(eltype(K))))
    all(real.(eigvals(k)) .< real.(zero(eltype(K))))
end

# TODO This does not work, since we are using the reduced basis instead of the
# original, lifted dynamics...
"""
$(SIGNATURES)

Update the Koopman `k` given new data `X` and `Y`. The operator is updated in place if
the L2 error of the prediction exceeds the `threshold`.

`p` and `t` are the parameters of the basis and the vector of timepoints, if necessary.
"""
function update!(k::AbstractKoopman,
    X::AbstractArray, Y::AbstractArray;
    p::AbstractArray = [], t::AbstractVector = [],
    U::AbstractArray = [],
    threshold::T = eps()) where {T <: Real}
    @assert updatable(k) "Linear Koopman is not updatable."

    Ψ₀ = k(X, p, t, U)
    Ψ₁ = k(Y, p, t, U)

    ϵ = norm(Ψ₁-Matrix(k)*Ψ₀, 2)

    if ϵ < threshold
        return
    end

    k.Q += Ψ₁*Ψ₀'
    k.P += Ψ₀*Ψ₀'
    k.operator .= k.Q / k.P

    if norm(Y - outputmap(k)*Matrix(k)*Ψ₀) < threshold
        return
    end

    # TODO Make this a proper rank 1 update
    k.output .= X / Ψ₀

    return
end


text:
# A Rudimentary clock built with papers and scissors
mutable struct Clock
    local_origin::DateTime # Time of clock creation
    _last::Millisecond # last pass of _tick
    delta::Millisecond # Cached version of Millisecond(now()-Clock.local_origin)
    tick::Float64 # Frequency of the clock in hertz
    _tick::Float64 # 1000/tick
    on_fire::Function # Handle for when the clock fires
end

function Clock(f::Function; tick::Float64=100.0) # default to 100hz
    if tick>1000.0 error("Tick must be below 1000hz") end
    c = Clock(now(), Millisecond(0), Millisecond(0), tick, 1000.0/tick, f)
    @spawn begin
        clock_running = true
        count = 0
        @sync while clock_running
            c.delta = Millisecond(now()-c.local_origin)
            if c.delta.value>c._tick+c._last.value
                count+=1
                if f(count, Millisecond(round(Int, datetime2unix(c.local_origin))), c.delta) == -1
                    clock_running = false
                end
                c._last = c.delta
            end
        end
    end
end


text:
function write_system(
  w::SystemWrite,
  K::SparseMatrixCSC{T, Int},
  rhs::Vector{T},
  step::Symbol,
  iter::Int,
) where {T <: Real}
  if rem(iter + 1 - w.kfirst, w.kgap) == 0
    if step != :cc
      K_str = string(w.name, "K_iter", iter + 1, ".mtx")
      MatrixMarket.mmwrite(K_str, K)
    end
    rhs_str = string(w.name, "rhs_iter", iter + 1, "_", step, ".rhs")
    open(rhs_str, "w") do io
      writedlm(io, rhs)
    end
  end
end


text:
abstract type AbstractRenewableFormulation end

abstract type AbstractRenewableDispatchForm <: AbstractRenewableFormulation end

struct RenewableCurtail <: AbstractRenewableDispatchForm end

struct RenewableCapacityCurve <: AbstractRenewableDispatchForm end

include("renewable_generation/renewable_variables.jl")
include("renewable_generation/output_constraints.jl")
include("renewable_generation/renewablegen_cost.jl")


text:
module MnistTest

using Learn
import MNIST
using MLDataUtils
using StatsBase

ENV["GKS_WSTYPE"] = "x11"
using StatPlots; gr(leg=false, linealpha=0.5)

# ----------------------------------------------------------------------------
# TODO: add this to MLPlots??

# a helper class to track many variables at once over time
type TracePlot{I,T}
    indices::I
    plt::Plot{T}
end
function TracePlot(n::Int = 1; maxn::Int = 500, kw...)
    indices = if n > maxn
        # sample maxn
        shuffle(1:n)[1:maxn]
    else
        1:n
    end
    plt = plot(length(indices); kw...)
    TracePlot(indices, plt)
end
function add_data(tp::TracePlot, x::Number, y::AbstractVector)
    for (i,idx) in enumerate(tp.indices)
        push!(tp.plt.series_list[i], x, y[idx])
    end
end
add_data(tp::TracePlot, x::Number, y::Number) = add_data(tp, x, [y])

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

function my_test_loss(obj, testdata, totcount = 500)
    totloss = 0.0
    totcorrect = 0
    for (x,y) in eachobs(rand(eachobs(testdata), totcount))
        totloss += transform!(obj,y,x)

        # logistic version:
        # ŷ = output_value(obj.transformation)[1]
        # correct = (ŷ > 0.5 && y > 0.5) || (ŷ <= 0.5 && y < 0.5)

        # softmax version:
        ŷ = output_value(obj.transformation)
        chosen_idx = indmax(ŷ)
        correct = y[chosen_idx] > 0

        totcorrect += correct
    end
    totloss, totcorrect/totcount
end

function doit()

    # our data:
    x_train, y_train = MNIST.traindata()
    x_test, y_test = MNIST.testdata()

    # normalize the input data given μ/σ for the input training data
    μ, σ = rescale!(x_train)
    rescale!(x_test, μ, σ)

    # note: needs my update to MLDataUtils
    y_train, y_test = map(to_one_hot, (y_train, y_test))

    # to_isone(y::AbstractVector) = (z = Array(eltype(y), 1, length(y)); map!(yi->float(yi==1.0), z, y))
    # y_train, y_test = map(to_isone, (y_train, y_test))

    # # keep only 0's and 1's
    # train = filterobs(i -> y_train[i] < 1.5, x_train, y_train)
    # test = filterobs(i -> y_test[i] < 1.5, x_test, y_test)

    train = (x_train, y_train)
    test = (x_test, y_test)

    # At this point we have train and test where each column of x is length-784 corresponding to pixel intensities,
    # and each column of y is length-10 corresponding to output class.

    nin, nh, nout = 784, [50,50], 10

    # build our objective... a neural net with nh hidden nodes,
    # tanh activation on the hidden layer, and logistic output
    t = nnet(nin, nout, nh, :softplus, :softmax)
    # penalty = L2Penalty(1e-3)
    penalty = ElasticNetPenalty(1e-5)
    obj = objective(t, penalty)
    @show obj

    # parameter plots
    pidx = 1:2:length(t)
    pvalplts = [TracePlot(length(params(t[i])), title="$(t[i])") for i=pidx]
    ylabel!(pvalplts[1].plt, "Param Vals")
    pgradplts = [TracePlot(length(params(t[i]))) for i=pidx]
    ylabel!(pgradplts[1].plt, "Param Grads")

    # nnet plots of values and gradients
    valinplts = [TracePlot(input_length(t[i]), title="input", yguide="Layer Value") for i=1:1]
    valoutplts = [TracePlot(output_length(t[i]), title="$(t[i])", titlepos=:left) for i=1:length(t)]
    gradinplts = [TracePlot(input_length(t[i]), yguide="Layer Grad") for i=1:1]
    gradoutplts = [TracePlot(output_length(t[i])) for i=1:length(t)]

    # loss/accuracy plots
    lossplt = TracePlot(title="Test Loss", ylim=(0,Inf))
    accuracyplt = TracePlot(title="Accuracy", ylim=(0.6,1))

    doanim = false
    # anim = Animation()

    # early_stopping = ConvergenceFunction((obj,i) -> begin
    #     false
    # end)

    tracer = IterFunction((obj, i) -> begin
        n = 100
        mod1(i,n)==n || return false

        # add param data
        for (j,k) in enumerate(pidx)
            add_data(pvalplts[j], i, params(t[k]))
            add_data(pgradplts[j], i, grad(t[k]))
        end

        # add input/output data
        for j=1:length(t)
            if j==1
                add_data(valinplts[j], i, input_value(t[j]))
                add_data(gradinplts[j], i, input_grad(t[j]))
            end
            add_data(valoutplts[j], i, output_value(t[j]))
            add_data(gradoutplts[j], i, output_grad(t[j]))
        end

        # sample points from the test set and compute/save the loss
        @show i
        if mod1(i,500)==500
            totloss, accuracy = my_test_loss(obj, test, 500)
            @show totloss, accuracy
            add_data(lossplt, i, totloss)
            add_data(accuracyplt, i, accuracy)
        end

        # build a heatmap of the total outgoing weight from each pixel
        pixel_importance = reshape(sum(t[1].params.views[1],1), 28, 28)
        hmplt = heatmap(pixel_importance, ratio=1)

        # build a nested-grid layout for all the trace plots
        getplt(p) = p.plt
        plot(
            map(getplt, vcat(
                    pvalplts, pgradplts,
                    valinplts, valoutplts,
                    gradinplts, gradoutplts,
                    lossplt, accuracyplt
                ))...,
            hmplt,
            size = (1400,1000),
            layout=@layout([
                grid(2,length(pvalplts))
                grid(2,length(valoutplts)+1)
                grid(1,3){0.2h}
            ])
        )

        if doanim
            lastframe = 5000
            if i < lastframe
                frame(anim)
            elseif i == lastframe
                gif(anim, fps=10)
            end
        end

        gui()
        # if i>0
        #     @profile gui()
        # else
        #     gui()
        # end
    end)

    # trace once before we start learning to see initial values
    tracer.f(obj, 0)

    # create a gradient descent learner and learn over infinite minibatches
    learner = make_learner(
        # GradientLearner(5e-3, RMSProp(0.9)),
        GradientLearner(5e-2, Adadelta()),
        # GradientLearner(1e-1, SGD(0.3)),
        # early_stopping,
        tracer,
        maxiter = 50000
    )
    learn!(obj, learner, infinite_batches(train, size=5))

    obj, learner
end

obj, learner = doit()

end #module

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization