text:
using JGUI
using Base.Test

# Box containers are manipulated like queues

# Adding objects

## push!
w = window(); f = vbox(w); push!(w, f)

## push!
btns = [push!(button(f, string("button $i"))) for i in 1:2]

## append! is more convenient
btns = [button(f, string("button $i")) for i in 3:4]
append!(f, btns)

## prepend!
prepend!(f, [button(f, "button 0")])

## insert! Index is 1:(n+1)
insert!(f, 3, button(f, "button 1.5"))

# removing widgets

## pop! (last)
pop!(f)

## shift! (first
shift!(f)

## splice! (by index)
splice!(f, 2)

## remove by item
item = f[2]
splice!(f, findin(item, f))

# indexing
f[2]


text:
module XEB

export brick_layout, google_layout_53, ustc_layout_60,
    add_noise!, gates, update_noise!, to_ein_code_pauli
export flux
export plot

include("layout.jl")
include("predefines.jl")
include("rewrite.jl")
include("cut.jl")
include("ein_code_pauli.jl")
include("ein_code_xeb.jl")
include("utils.jl")
include("yao_block.jl")
include("plots.jl")
include("plot3d.jl")

end

ext: jl


text:
"""
Definition of `DynamicalSystem` as well as all integrators
used in the ecosystem of DynamicalSystems.jl.

Also contains predefined well-known systems.
"""
module DynamicalSystemsBase

using DelayEmbeddings

include("deprecations.jl")

include("core/dynamicalsystem.jl")
include("core/create_dynamics.jl")
include("core/api_docstrings.jl")
include("core/discrete.jl")
include("core/continuous.jl")

include("predefined/famous_systems.jl") # includes both discrete and continuous

export Systems, reinit!
export SVector, SMatrix, @SVector, @SMatrix

end


text:
# unit test for the options factory
using COSMO, Test, LinearAlgebra, SparseArrays, Pardiso

@testset "Options Factory" begin

P = spzeros(1, 1)
q = [1]
c1 = COSMO.Constraint([1.], [-1.], COSMO.Nonnegatives)

iparms = Dict{Int64, Int64}()
iparms[13] = 100
num_threads = 2

settings = COSMO.Settings(kkt_solver = with_options(COSMO.MKLPardisoKKTSolver, iparm = iparms, num_threads = num_threads, msg_level_on = false))
model = COSMO.Model()
assemble!(model, P, q, c1, settings = settings)
model.ρvec = settings.rho * ones(size(model.p.A, 1))
COSMO._make_kkt_solver!(model)

  @test typeof(model.kkt_solver) <: COSMO.MKLPardisoKKTSolver
  # @test Pardiso.get_nprocs(model.kkt_solver.ps) == num_threads
  @test Pardiso.get_iparm(model.kkt_solver.ps, 13) == iparms[13]

  ## check that this works via the MathOptInterface wrapper as well
  optimizer = COSMO.Optimizer(time_limit = 10., kkt_solver = with_options(COSMO.MKLPardisoKKTSolver, iparm = iparms, num_threads = num_threads, msg_level_on = false))
  @test model.settings.kkt_solver == optimizer.inner.settings.kkt_solver

  # assemble model when only KKT solver type is provided
  settings = COSMO.Settings(kkt_solver = COSMO.QdldlKKTSolver)
  model = COSMO.Model()
  assemble!(model, P, q, c1, settings = settings)
  model.ρvec = settings.rho * ones(size(model.p.A, 1))
  COSMO._make_kkt_solver!(model)
  @test typeof(model.kkt_solver) <: COSMO.QdldlKKTSolver

  # try to assemble with wrong type
  struct MockAccelerator end
  @test_throws MethodError COSMO.Settings(kkt_solver = with_options(MockAccelerator, mem = 5))

end


text:
import Statistics
import Graphs

include("STOCK.jl")

# 1. compute stock price daily changes
X = STOCK.Closes - STOCK.Opens

# 2. compute *squared* correlation matrix
Xcorr = Statistics.cor(X).^2

# 3. compute *maximum* spanning tree
edges = Graphs.kruskal_mst(Graphs.SimpleGraphs.complete_graph(size(Xcorr, 1)), Xcorr, minimize=false)
T = Graphs.SimpleGraphFromIterator(edges)

# 4. write a graphviz file for the graph, see Makefile for generating the pdf figure
open("stocks.graphviz", "w") do file
  write(file, "graph G {\n")

  for (i, name) in enumerate(STOCK.names)
    write(file, "\t\"$(name)\" [style=filled fillcolor=$(STOCK.colors[STOCK.industries[i]])];\n")
  end

  for (i, pa) in enumerate(Graphs.dfs_tree(T, 1).badjlist)
    for p in pa
      write(file, "\t\"$(STOCK.names[i])\" -- \"$(STOCK.names[p])\";\n")
    end
  end
  write(file, "}")
end


text:
# methods which check the subtype relationships between NewTypes

const FV = FreeVariables

isnewsubtype(a::NewType, b::NewType) = issubenv(a, b, FV())
isnewsubtype(a::Type, b::Type) = (a <: b)

"""
    `issubenv(a, b, env)`

Test if type a is subtype of type b with the given variables and their bindings.
"""
issubenv(::NewBottomType, b::NewType, env::FV) = true
issubenv(a::NewType, b::NewUnion, env::FV) = issubenv(a, b.a, env) || issubenv(a, b.b, env)
issubenv(a::NewUnion, b::NewType, env::FV) = issubenv(a.a, b, env) && issubenv(a.b, b, env)
issubenv(a::NewUnion, b::NewUnion, env::FV) = issubenv(a.a, b, env) && issubenv(a.b, b, env)

# is any of the "supertypes" of `a` a direct subtype of `b`?
issubenv(a::NewType, b::NewType, env::FV) = any(isdirect1(st, b, env) for st in linearize(a))

# perform trivial checks before calling specialized functions
function isdirect1(a::NewType, b::NewType, env::FV)
    a === b || b == NewAny || isdirect(a, b, env)
end

# maintain stack of free variables and delegate to body
isdirect(a::NewType, b::NewUnionAll, env::FV) = withvariable(isdirect, a, b, env)

# check name identity then specialize for ordinary data types and tuple types
function isdirect(a::NewType, b::NewDataType, env::FV)
    a = typebody(a)
    a.name == b.name || return false
    if isatuple(b)
        walktuple(a, b, env, iscovar)
    else
        walkordinary(a, b, env, isinvar)
    end
end

# handle the parameter lists of ordinary types
function walkordinary(a::NewDataType, b::NewDataType, env::FV, func::Function)
    pa = a.parameters; na = length(pa)
    pb = b.parameters; nb = length(pb)
    na == nb && all(func(pa[k], pb[k], env) for k = 1:na)
end

# handle the parameter lists of tuples (with varargs)
function walktuple(a::NewDataType, b::NewDataType, env::FV, func::Function)
    pa = a.parameters; na = length(pa)
    pb = b.parameters; nb = length(pb)
    nb == 0 && return na == 0
    tbb = typebody(pb[nb])
    if na == 0
        return nb == 1 && isavararg(tbb) && !isa(tbb.parameters[2], Integer)
    end
    tba = typebody(pa[na])
    if isavararg(tbb)
        vtb = tbb.parameters[1]
        vlb = tbb.parameters[2]
        absorbb = 0
        if isavararg(tba)
            nc = min(na, nb)
            all(func(pa[k], pb[k], env) for k = 1:nc-1) || return false
            vta = tba.parameters[1]
            vla = tba.parameters[2]
            if isa(vla, Integer)
                if isa(vlb, Integer)
                    na + vla == nb + vlb || return false
                    vla == 0 || vlb == 0 || func(vta, vtb, env) || return false
                else
                    na + vla >= nb || return false
                    vla == 0 || func(vta, vtb, env) || return false
                    testbinding!(env, vlb, na - nb + vla)
                end
            else
                if isa(vlb, Integer)
                    return false
                else
                    na >= nb || return false
                    func(vta, vtb, env) || return false
         #TODO: bind type variable vlb to vla + constant!!!
                    testbinding!(env, vlb, vla) # + (nb - na)
                end
            end
            all(func(vta, pb[k], env) for k = nc:nb-1) &&
            all(func(pa[k], vtb, env) for k = nc:na-1)
        else
            na >= nb - 1 || return false
            all(func(pa[k], pb[k], env) for k=1:nb-1) || return false
            if isa(vlb, Integer)
                nb - 1 + vlb == na || return false
            end
            all(func(pa[k], vtb, env) for k = nb:na) &&
            testbinding!(env, vlb, na - nb + 1)
        end
    else
        !isavararg(tba) && na == nb && all(func(pa[k], pb[k], env) for k=1:nb)
    end
end

# parameters checking if type variables may be involved
# parameters can be NewVarType, NewType, or any bits type
# if no type variables are involved, invariant version requires type equality
# while covariant version tests the subtype relation.
# if on of the arguments is a type variable, there is no difference
# between covariant and invariant case.

iscovar(a, b, env::FV) = isinvar(a, b, env)
iscovar(a::NewType, b::NewType, env::FV) = issubenv(a, b, env)

isinvar(a, b, env::FV) = equaltypes(a, b, env)
isinvar(a::NewTypeVar, b, env::FV) = false
#isinvar(a::NewType, b::NewType, env::FV) = equaltypes(a, b, env)
isinvar(a::NewTypeVar, b::NewType, env::FV) = false

# if the second argument is a type variable, the first argument is bound to it
# if it is already bound to another value false is returned
isinvar(a, b::NewTypeVar, env::FV) = testbinding!(env, b, a)
isinvar(a::NewTypeVar, b::NewTypeVar, env::FV) = testbinding!(env, b, a)

# utility functions

function withvariable(f::Function, a, b::NewUnionAll, env::FV)
    push!(env, b.var)
    r = f(a, b.body, env) && checkbounds(binding(env, b.var), b.var, env)
    pop!(env, b.var)
    r
end

checkbounds(a, b, env::FV) = issubenv(a, b, env)
checkbounds(::Void, v::NewTypeVar, env::FV) = throw(UndefVarError(v.name))
checkbounds(a::Any, v::NewTypeVar, env::FV) = a == binding(env, v)
function checkbounds(a::NewType, v::NewTypeVar, env::FV)
    checkbounds(v.lb, a, FV()) &&
    checkbounds(a, v.ub, env)
end

function checkbounds(a::NewTypeVar, b::NewTypeVar, env::FV)
    checkbounds(a.ub, b.ub, env) &&
    checkbounds(b.lb, a.lb, FV())
end

equaltypes(a, b, env::FV) = a == b
equaltypes(a, b::NewTypeVar, env::FV) = isinvar(a, b, env)
equaltypes(a, b::NewUnionAll, env::FV) = withvariable(equaltypes, a, b, env)

equaltypes(a::NewType, b::NewDataType, env::FV) = equaltypes(typebody(a), b, env)

function equaltypes(a::NewDataType, b::NewDataType, env::FV)
    a.name == b.name || return false
    if isatuple(a)
        walktuple(a, b, env, equaltypes)
    else
        walkordinary(a, b, env, equaltypes)
    end
end

typebody(x::Any) = x
typebody(x::NewUnionAll) = typebody(x.body)

isavararg(x) = isaspecial(x, :Vararg)
isatuple(x) = isaspecial(x, :Tuple)
isaspecial(x, ::Symbol) = false
isaspecial(x::NewDataType, s::Symbol) = x.name.name == s
isaspecial(x::NewUnionAll, s::Symbol) = isaspecial(typebody(x), s)


text:
using Images, AxisArrays
using Base: OneTo
using Unitful: Time, Length, s

struct Acquisition <: AbstractAcquisition # {{{
    "Acquired images"
    images::AxisArray{<:AbstractGray,3}
    "Numerical aperature"
    NA::Real
    "Magnification"
    M::Integer
    "Fluorophore marker"
    marker::FluorophoreMarker
    "Pixel size"
    pixelsize::Length
end # }}}

function Acquisition( # {{{
    images::AbstractArray{<:Color,3},
    acquisitiontimes::AbstractVector{<:Time},
    NA,
    M,
    marker::FluorophoreMarker,
    pixelsize::Length,
)
    if !(eltype(images) <: AbstractGray)
        images = Gray.(images)
    end
    images_axis = AxisArray(
        images,
        Axis{:x}(OneTo(size(images, 1))),
        Axis{:y}(OneTo(size(images, 2))),
        Axis{:time}(acquisitiontimes),
    )
    return Acquisition(images_axis, NA, M, marker, pixelsize)
end # }}}

function Acquisition( # {{{
    images::AbstractArray{<:Color,3},
    acquisitiontimestep::Time,
    NA,
    M,
    marker::FluorophoreMarker,
    pixelsize::Length,
)
    return Acquisition(
        images,
        range(0s; length=size(images, 3), step=acquisitiontimestep),
        NA,
        M,
        marker,
        pixelsize,
    )
end # }}}

times(acquisition::AbstractAcquisition) = acquisition.images[Axis{:time}].val


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

using Test
using Random, SimilaritySearch, KCenters, StatsBase

const X = [rand(4) for i in 1:1000]

@testset "Clustering with enet" begin
    for i in 2:5
        p = enet(L2Distance(), X, i^2)
        D = [minimum(p) for p in p.seq]
        @test maximum(D) <= p.dmax
    end
end

@testset "Clustering with dnet" begin
    for i in 2:5
        res = dnet(L2Distance(), X, i^2)
        @info mean([minimum(p) for p in res.seq])
    end
end

@testset "Clustering with KCenters" begin
    cfft = KCenters.kcenters(L2Distance(), X, 16)
    cdnet = KCenters.kcenters(L2Distance(), X, 16, initial=:dnet)
    crand = KCenters.kcenters(L2Distance(), X, 16, initial=:rand)
    @show mean(cfft.distances)
    @show mean(cdnet.distances)
    @show mean(crand.distances)
end

@testset "Clustering with KCenters; KnnCentroidSelection" begin
    c1 = KCenters.kcenters(L2Distance(), X, 16; sel=KnnCentroidSelection())
    c2 = KCenters.kcenters(L2Distance(), X, 16; sel=KnnCentroidSelection(sel1=MedoidSelection()))
    c3 = KCenters.kcenters(L2Distance(), X, 16; sel=KnnCentroidSelection(sel1=RandomCenterSelection()))
    d1 = mean(c1.distances)
    d2 = mean(c2.distances)
    d3 = mean(c3.distances)
    @show d1 d2 d3
    @test abs(d1 - d2) < 0.2
    @test abs(d1 - d3) < 0.2
end

@testset "Clustering with KCenters with an approximate index" begin
    cfft = KCenters.kcenters(L2Distance(), X, 16, recall=0.99)
    cdnet = KCenters.kcenters(L2Distance(), X, 16, initial=:dnet, recall=0.99)
    crand = KCenters.kcenters(L2Distance(), X, 16, initial=:rand, recall=0.99)
    d1 = mean(cfft.distances)
    d2 = mean(cdnet.distances)
    d3 = mean(crand.distances)
    @show d1 d2 d3
    @test abs(d1 - d2) < 0.2
    @test abs(d1 - d3) < 0.2
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using FEMBase

mutable struct CraigBampton <: AbstractAnalysis
    r_nodes :: Vector{Int64}
    l_nodes :: Vector{Int64}
    K :: SparseMatrixCSC{Float64}
    M :: SparseMatrixCSC{Float64}
    K_red :: Matrix{Float64}
    M_red :: Matrix{Float64}
    B :: Matrix{Float64}
    dim :: Int64
    nev :: Int64
    time ::Float64
end

function CraigBampton()
    return CraigBampton(
        zeros(Int64, 0),
        zeros(Int64, 0),
        spzeros(0, 0),
        spzeros(0, 0),
        zeros(0, 0),
        zeros(0, 0),
        zeros(0, 0),
        3, 10, 0.0)
end

"""
    get_global_matrices(cb)

Return global matrices K and M. If matrices are already defined to
`cb.properties.K` and `cb.properties.M`, return those instead.
"""
function get_global_matrices(cb)
    if !isempty(cb.properties.K) && !isempty(cb.properties.M)
        return cb.properties.K, cb.properties.M
    end

    t0 = time()
    info("Assembling global matrices")
    assemble!(cb, cb.properties.time; with_mass_matrix=true)

    M = SparseMatrixCOO()
    K = SparseMatrixCOO()
    Kg = SparseMatrixCOO()
    f = SparseMatrixCOO()
    fg = SparseMatrixCOO()

    for problem in get_problems(cb)
        is_field_problem(problem) || continue
        append!(M, problem.assembly.M)
        append!(K, problem.assembly.K)
        append!(Kg, problem.assembly.Kg)
        append!(f, problem.assembly.f)
        append!(fg, problem.assembly.fg)
    end

    dim = size(K, 1)

    M = sparse(M, dim, dim)
    K = sparse(K, dim, dim) + sparse(Kg, dim, dim)
    f = sparse(f, dim, 1) + sparse(fg, dim, 1)

    for problem in get_problems(cb)
        is_boundary_problem(problem) || continue
        eliminate_boundary_conditions!(problem, K, M, f)
    end

    K = 1/2*(K + K')
    M = 1/2*(M + M')

    SparseArrays.droptol!(K, 1.0e-9)
    SparseArrays.droptol!(M, 1.0e-9)
    dt = round(time() - t0)
    info("Assembly of global matrices done in $dt seconds")

    cb.properties.K = K
    cb.properties.M = M
    return K, M
end

function FEMBase.run!(cb::Analysis{CraigBampton})
    dim = cb.properties.dim
    nev = cb.properties.nev
    r_nodes = cb.properties.r_nodes
    l_nodes = cb.properties.l_nodes
    r_dofs = [dim*(i-1)+j for i in r_nodes for j=1:dim]
    l_dofs = [dim*(i-1)+j for i in l_nodes for j=1:dim]
    info("number of retained nodes = ", length(r_nodes))
    info("number of retained dofs = ", length(r_dofs))
    info("number of other nodes = ", length(l_nodes))
    info("number of other dofs = ", length(l_dofs))

    K, M = get_global_matrices(cb)

    K_LL = K[l_dofs, l_dofs]
    K_RR = K[r_dofs, r_dofs]
    K_RL = K[r_dofs, l_dofs]
    K_LR = K[l_dofs, r_dofs]

    M_LL = M[l_dofs, l_dofs]
    M_RR = M[r_dofs, r_dofs]
    M_RL = M[r_dofs, l_dofs]
    M_LR = M[l_dofs, r_dofs]

    info("1. Static condensation")

    t0 = time()
    info("Construct recovery matrix B = K_LL^-1 * K_LR")
    info("Factorize K_LL using Cholesky factorization")
    nz = get_nonzero_rows(K_LL)
    CF = cholfact(K_LL[nz,nz])
    info("Calculate B = K_LL^-1 * K_LR blockwise")
    dim1, dim2 = size(K_LR)
    cb.properties.B = B = zeros(dim1, dim2)
    nz2 = get_nonzero_columns(K_LR)
    for j=1:length(nz2)
        if mod(j,10) == 0
            info("Condensate dof $j / $(length(nz2))...")
        end
        B[nz,nz2[j]] = CF \ full(K_LR[nz,nz2[j]])
    end
    dt = round(time() - t0, 2)
    info("Calculation of recovery matrix B done in $dt seconds")

    t0 = time()
    info("Calculate boundary stiffness matrix K_BB")
    K_BB = K_RR - K_RL*B
    dt = round(time() - t0, 2)
    info("Calculating K_BB done in $dt seconds")

    t0 = time()
    info("Calculate boundary mass matrix M_BB")
    info("Calculate Y = M_RL*B")
    Y = -M_RL*B
    Z = M_LL*B
    info("Calculate M_BB = M_RR + Y + Y' + B'*M_LL*B")
    M_BB = M_RR + Y + Y' + B'*Z
    dt = round(time() - t0, 2)
    info("Calculating M_BB done in $dt seconds")

    info("2. Dynamic condensation")

    info("Calculate eigenvalue problem K_LL + w^2*M_LL = 0")
    t0 = time()
    X_LL = zeros(size(K_LL, 1), nev)
    w_LL, X_LL[nz,:] = eigs(K_LL[nz,nz], M_LL[nz,nz]; which=:SM, nev=nev)
    w_LL = real(w_LL)
    dt = round(time() - t0, 2)
    info("Eigenvalues solved in $dt seconds")

    info("Construct the final system M_red*u'' + K_red*u = f_red")
    t0 = time()
    K_mm = X_LL' * K_LL * X_LL
    M_mm = X_LL' * M_LL * X_LL
    # K_mm and M_mm should be diagonal and contains some roundoff errors ...
    K_mm[abs.(K_mm) .< 1.0e-6] = 0.0
    M_mm[abs.(K_mm) .< 1.0e-6] = 0.0
    K_mB = zeros(nev, length(r_dofs))
    M_mB = X_LL' * (M_LR + Z)
    K_Bm = K_mB'
    M_Bm = M_mB'
    cb.properties.M_red = [M_BB M_Bm; M_mB M_mm]
    cb.properties.K_red = [K_BB K_Bm; K_mB K_mm]
    dt = round(time() - t0, 2)
    info("Reduced system ready in $dt seconds")

    nr_dofs = length(r_dofs)
    nl_dofs = length(l_dofs)
    nt_dofs = nr_dofs + nl_dofs
    p = round(nr_dofs/nt_dofs*100, 2)

    info("--- CALCULATION SUMMARY ---")
    info("Dimension of original system: $nt_dofs dofs")
    info("Dimension of reduced system: boundary $nr_dofs dofs + $nev general dofs")
    info("Dimension of reduced system is $p % of the original system")

    return cb.properties.M, cb.properties.K
end


text:
"""
    sample_vdpoisson(shape::Tuple{Int64},redFac;seed::Int64=1234,kargs...)

generates a Poisson disk sampling pattern for an Vector of size `shape` with a subsampling factor `redFac`.
The arguments are the same as in the 2d case with the exception that shape is of type `Tuple{Int64}`
"""
function sample_vdpoisson(shape::Tuple{Int64},redFac;seed::Int64=1234,kargs...)
  return sample_vdpoisson((shape[1],1),redFac,seed=seed;kargs...)
end

function sample_vdpoisson(shape::Tuple{Int64,Int64,Int64},redFac;seed::Int64=1234,kargs...)
  error("Not implemented")
end

"""
    sample_vdpoisson(shape::Tuple{Int64,Int64},redFac::Float64; seed::Int64=1234,kargs...)

generates a variable density Poisson disk sampling pattern for an Array of size `shape` with a subsampling factor `redFac`.

# Arguments
* `shape::NTuple{2,Int64}` - size of the Array to be sampled
* `redFac::Float64`        - subsampling factor
* (`seed=1234`)            - seed for the random number generator
"""
function sample_vdpoisson(shape::Tuple{Int64,Int64},redFac::Float64; seed::Int64=1234,densityFac::Float64=1/pi,kargs...)
  M,N = shape
  samples = collect(1:M*N)
  chosenSamples = []
  forbiddenSamples =[]

  # set number of samples and mean distance
  numSamples = floor(Int64,M*N/redFac)
  rmean = densityFac*M*N/numSamples

  # srand(seed)
  while length(chosenSamples)<numSamples
    # check if samples are available. Otherwise reduce rmean and recalculate samples
    if isempty(samples)
      rmean  = rmean/2.0
      forbiddenSamples = updateForbiddenSamples(M,N,chosenSamples,rmean)
      samples = setdiff(collect(1:M*N),forbiddenSamples)
      continue
    end

    # pick a sample from the available samples and add them to chosen samples
    chosen = samples[ rand(1:length(samples)) ]
    push!(chosenSamples,chosen)

    # calculate forbidden radius and points within
    fSphere = forbiddenSphere(M,N,chosen,rmean)
    forbiddenSamples = union(forbiddenSamples,fSphere)

    # update available samples
    samples = setdiff(samples,forbiddenSamples)
  end

  return chosenSamples
end

function rmin(kx::Float64, ky::Float64, ai::Float64, a0::Float64, rmean::Float64)
  m = rmean/(ai-3*a0)
  c = rmean*(ai-5*a0)/(2*(ai-3*a0))
  return m*sqrt(kx^2+ky^2)+c
end

function updateForbiddenSamples(M,N,samples,rmean)
  forbiddenSamples = []
  for i=1:length(samples)
    fSphere = forbiddenSphere(M,N,samples[i],rmean)
    forbiddenSamples = union(forbiddenSamples,fSphere)
  end
  return forbiddenSamples
end

function forbiddenSphere(M::Int64,N::Int64,chosen::Int64,rmean::Float64)
  x,y = Tuple(CartesianIndices((M,N))[chosen])
  kx = x-(M+1.)/2.0
  ky = y-(N+1.)/2.0
  rad_min = rmin(kx, ky, 0.5, 0.0, rmean)
  fsphere = []
  for xi = floor(Int64,x-M*rad_min):ceil(Int64,x+M*rad_min)
    # xi must be positive and smaller then M
    (xi<1) && continue
    (xi>M) && continue

    for yi =  floor(Int64,y-N*rad_min):ceil(Int64,y+N*rad_min)
      # yi must be positive and smaller then N
      (yi<1) && continue
      (yi>N) && continue

      if (xi-x)^2+(yi-y)^2 < rad_min^2
        idx = LinearIndices((M,N))[xi,yi]
        push!(fsphere,idx)
      end
    end
  end
  return fsphere
end

function generateMask(M,N,pat)
  mask = zeros(Int64,M,N)
  for i=1:length(pat)
    mask[pat[i]] = 1.0
  end
  return mask
end


text:
# This file was generated, do not modify it. # hide
mach = MLJ.machine(nnregressor, features, targets)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct DictVectorizer <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function DictVectorizer(; kwargs...)
        obj = new(meta(DictVectorizer), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct DictVectorizer
const __meta_DictVectorizer = Ref{ProtoMeta}()
function meta(::Type{DictVectorizer})
    ProtoBuf.metalock() do
        if !isassigned(__meta_DictVectorizer)
            __meta_DictVectorizer[] = target = ProtoMeta(DictVectorizer)
            allflds = Pair{Symbol,Union{Type,String}}[:stringToIndex => StringVector, :int64ToIndex => Int64Vector]
            oneofs = Int[1,1]
            oneof_names = Symbol[Symbol("Map")]
            meta(target, DictVectorizer, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_DictVectorizer[]
    end
end
function Base.getproperty(obj::DictVectorizer, name::Symbol)
    if name === :stringToIndex
        return (obj.__protobuf_jl_internal_values[name])::StringVector
    elseif name === :int64ToIndex
        return (obj.__protobuf_jl_internal_values[name])::Int64Vector
    else
        getfield(obj, name)
    end
end

export DictVectorizer


text:
using ParserCombinator
using REPL: REPL, LineEdit
using ReplMaker

function lisp_reader(s)
    try
        read(String(take!(copy(LineEdit.buffer(s)))))
        true
    catch err
        isa(err, ParserCombinator.ParserException) || rethrow(err)
        false
    end
end

function init_repl(; prompt_text="jλ> ", prompt_color = :red, start_key = ")", sticky=true)
	    ReplMaker.initrepl(lisp_eval_helper,
	                  repl = Base.active_repl,
	                  valid_input_checker = lisp_reader,
	                  prompt_text = prompt_text,
	                  prompt_color = prompt_color,
	                  start_key = start_key,
                          sticky_mode=sticky,
	                  mode_name = "Lisp Mode")
end


text:
function  exasim(pde,mesh)

if isa(pde, Array)
    nmodels = length(pde);
else
    nmodels = 1;
end

res = []
if nmodels==1
    # search compilers and set options
    pde = Gencode.setcompilers(pde);

    # generate input files and store them in datain folder
    pde, mesh, master, dmd = Preprocessing.preprocessing(pde,mesh);

    # generate source codes and store them in app folder
    Gencode.gencode(pde);

    if pde.usecmake==1
        runstr = Main.cmakecompile(pde, nmodels); # use cmake to compile source codes
        compilerstr = "";
    else
        # compile source codes to build an executable file and store it in app folder
        compilerstr = Gencode.compilecode(pde);

        # run executable file to compute solution and store it in dataout folder
        runstr = Gencode.runcode(pde, nmodels);
    end

    # get solution from output files in dataout folder
    sol = Postprocessing.fetchsolution(pde,master,dmd,"dataout");
    if pde.saveResNorm == 1
        fn = "dataout/out_residualnorms0.bin";
        res = reinterpret(Float64,read(fn));
        ne = Int64(round(length(res)/4));
        res = reshape(res,(4,ne));
        res = res';
    end
else
    master = Array{Any, 1}(undef, nmodels);
    dmd = Array{Any, 1}(undef, nmodels);
    sol = Array{Any, 1}(undef, nmodels);
    res = Array{Any, 1}(undef, nmodels);

    # preprocess and generate code for all PDE models
    for m = 1:nmodels
        pde[m],mesh[m],master[m],dmd[m] = generatecode(pde[m],mesh[m]);
    end
    gencodeall(nmodels);

    if pde.usecmake==1
        runstr = Main.cmakecompile(pde[1], nmodels); # use cmake to compile source codes
        compilerstr = "";
    else
        # # compile source codes to build an executable file and store it in app folder
        compilerstr = compilecode(pde[1]);

        # run executable file to compute solution and store it in dataout folder
        runstr = runcode(pde[1],nmodels);
    end

    # get solution from output files in dataout folder
    for m = 1:nmodels
        sol[m] = fetchsolution(pde[m], master[m], dmd[m], "dataout" * string(m));

        if pde[m].saveResNorm
            fn = "dataout/out_residualnorms" * string(m-1) * ".bin";
            tm = reinterpret(Float64,read(fn));
            ne = Int64(round(length(tm)/4));
            tm = reshape(tm,(4,ne));
            res[m] = tm';
        end
    end
end

return sol,pde,mesh,master,dmd,compilerstr,runstr,res

end

max_stars_repo_name: e-dinesh/Exasim

max_issues_repo_name: e-dinesh/Exasim

max_forks_repo_name: e-dinesh/Exasim


text:
using ComputerVision
using Test
using StaticArrays, GeometryTypes

@testset "ComputerVision.jl" begin
    # Write your own tests here.
    @testset "Corresponence Test" begin include("correspondence.jl") end
    @testset "Camera Test" begin include("camera.jl") end
    @testset "Planar Scene" begin include("planar_scene.jl") end
end


text:
"""
包含最基础功能的模块
"""
module Basics
    include("point.jl")
    export Point2D
    export middle_point, absolute_angle

    include("line.jl")
    export Segment

    include("polygon.jl")
    export AbstractPolygon

    include("triangle.jl")
    export Triangle, area
    export RtTriangle, EqTriangle

    include("rectangle.jl")
    export Rectangle, Square, area

    include("hexagon.jl")
    export EqHexagon, area

end


text:
@testset "Places" begin
   Qx, x = FlintQQ["x"]
   K1, a1 = NumberField(x^5 + 2, "a1")
   K2, a2 = NumberField(x^9 + 1*x^4 + 14*x^3 - 14*x^2 + 14*x - 14, "a2")
   K3, a3 = NumberField(x^6 - 2, "a3")

  @testset "Construction" begin

    S = infinite_places(K1)
    @test length(S) == 3

    SR = real_places(K1)
    @test length(SR) == 1
    @test all(isreal, SR)

    SC = complex_places(K1)
    @test length(SC) == 2
    @test all(iscomplex, SC)

    r = conjugates_arb(a1)
    @test overlaps(r[1], SR[1].r)
    @test overlaps(r[2], SC[1].r) && overlaps(r[3], SC[2].r)

    P = infinite_place(K1, 1)
    @test P == S[1]
    @test overlaps(P.r, r[1])

    P = infinite_place(K1, 2)
    @test P == S[2]
    @test overlaps(P.r, r[2])

    @test S[1] == S[1]
    @test S[1] != S[2]

    string(S[1]); # Just check that it does not error

    S = infinite_places(K2)
    @test length(S) == 5

    SR = real_places(K2)
    @test length(SR) == 1
    @test all(isreal, SR)

    SC = complex_places(K2)
    @test length(SC) == 4
    @test all(iscomplex, SC)

    r = conjugates_arb(a2)
    @test overlaps(r[1], SR[1].r)
    @test overlaps(r[2], SC[1].r) && overlaps(r[3], SC[2].r) && overlaps(r[4], SC[3].r) && overlaps(r[5], SC[4].r)

    P = infinite_place(K2, 1)
    @test P == S[1]
    @test overlaps(P.r, r[1])

    P = infinite_place(K2, 2)
    @test P == S[2]
    @test overlaps(P.r, r[2])

    @test S[1] == S[1]
    @test S[1] != S[2]

    string(S[1]); # Just check that it does not error
  end

  @testset "Signs and Positivity" begin
    b = a1
    P = real_places(K1)[1]
    C = complex_places(K1)[1]

    for b in [b, FacElem(b), maximal_order(K1)(b)]
      sb = signs(b)
      @test sb == Dict(P => -1)
      @test sign(b, P) == -1
      @test_throws ErrorException sign(b, C)
      @test signs(b, infinite_places(K1)) == Dict(P => -1)
      @test signs(b, [P]) == Dict(P => -1)
      @test signs(b, complex_places(K1)) == Dict{InfPlc, Int}()
      @test !ispositive(b, P)
      @test_throws ErrorException ispositive(b, C)
      @test !ispositive(b, [P])
      @test !ispositive(b, infinite_places(K1))
      @test ispositive(b, [C])
      @test !istotally_positive(b)

      c = b^10*b^2
      sc = signs(c)
      @test sc == Dict(P => 1)
      @test sign(c, P) == 1
      @test_throws ErrorException sign(c, C)
      @test signs(c, infinite_places(K1)) == Dict(P => 1)
      @test signs(c, [P]) == Dict(P => 1)
      @test signs(c, complex_places(K1)) == Dict{InfPlc, Int}()
      @test ispositive(c, P)
      @test_throws ErrorException ispositive(c, C)
      @test ispositive(c, [P])
      @test ispositive(c, [C])
      @test ispositive(c, infinite_places(K1))
      @test istotally_positive(c)
    end

    b = a3
    P = real_places(K3)[1]
    P2 = real_places(K3)[2]
    C = complex_places(K3)[1]

    for b in [b, FacElem(b), maximal_order(K3)(b)]
      @test signs(b) == Dict(P => -1, P2 => 1)
      @test sign(b, P) == -1
      @test_throws ErrorException sign(b, C)
      @test signs(b, infinite_places(K3)) == Dict(P => -1, P2 => 1)
      @test signs(b, [P]) == Dict(P => -1)
      @test signs(b, complex_places(K3)) == Dict{InfPlc, Int}()
      @test ispositive(b, P2)
      @test !ispositive(b, P)
      @test_throws ErrorException ispositive(b, C)
      @test ispositive(b, [P2])
      @test ispositive(b, [P2, C])
      @test !ispositive(b, [P])
      @test !ispositive(b, [P, C])
      @test !ispositive(b, infinite_places(K3))
      @test !istotally_positive(b)

      c = b^10*b^2
      @test signs(c) == Dict(P => 1, P2 => 1)
      @test sign(c, P) == 1
      @test_throws ErrorException sign(c, C)
      @test signs(c, infinite_places(K3)) == Dict(P => 1, P2 => 1)
      @test signs(c, [P]) == Dict(P => 1)
      @test signs(c, complex_places(K3)) == Dict{InfPlc, Int}()
      @test ispositive(c, P)
      @test ispositive(c, P2)
      @test_throws ErrorException ispositive(c, C)
      @test ispositive(c, [P])
      @test ispositive(c, [P, C])
      @test ispositive(c, infinite_places(K3))
      @test istotally_positive(c)
    end
  end
end


text:
function get_data(dataobject::PartDataType,
                vars::Array{Symbol,1},
                units::Array{Symbol,1},
                direction::Symbol,
                center::Array{<:Any,1},
                mask::MaskType,
                ref_time::Real)

    vars_dict = Dict()
    #vars = unique(vars)
    boxlen = dataobject.boxlen

    if direction == :z
        apos = :x
        bpos = :y
        cpos = :z

        avel = :vx
        bvel = :vy
        cvel = :vz

    elseif direction == :y
        apos = :z
        bpos = :x
        cpos = :y

        avel = :vz
        bvel = :vx
        cvel = :vy
    elseif direction == :x
        apos = :z
        bpos = :y
        cpos = :x

        avel = :vz
        bvel = :vy
        cvel = :vx
    end

    column_names = propertynames(dataobject.data.columns)

    for i in vars

        # quantities that are in the datatable
        if in(i, column_names)
            selected_unit = getunit(dataobject, i, vars, units)

            if i == :x
                #selected_unit = (dataobject, :x, vars, units)
                vars_dict[:x] =  (select(dataobject.data, apos) .-  boxlen * center[1] )  .* selected_unit
            elseif i == :y
                #selected_unit = (dataobject, :y, vars, units)
                vars_dict[:y] =  (select(dataobject.data, bpos) .- boxlen * center[2] )  .* selected_unit
            elseif i == :z
                #selected_unit = (dataobject, :z, vars, units)
                vars_dict[:z] =  (select(dataobject.data, cpos)  .- boxlen * center[3] )  .* selected_unit
            else
                vars_dict[i] =  select(dataobject.data, i) .* selected_unit
            end

        # quantities that are derived from the variables in the data table
        elseif i == :vx2
            selected_unit = getunit(dataobject, :vx2, vars, units)
            vars_dict[:vx2] =  (getvar(dataobject, :vx) .* selected_unit ) .^2

        elseif i == :vy2
            selected_unit = getunit(dataobject, :vy2, vars, units)
            vars_dict[:vy2] =  (getvar(dataobject, :vy) .* selected_unit ) .^2

        elseif i == :v
            selected_unit = getunit(dataobject, :v, vars, units)
            vars_dict[:v] =  sqrt.(select(dataobject.data, :vx).^2 .+
                                   select(dataobject.data, :vy).^2 .+
                                   select(dataobject.data, :vz).^2 ) .* selected_unit
        elseif i == :v2
           selected_unit = getunit(dataobject, :v2, vars, units)
           vars_dict[:v2] = (select(dataobject.data, :vx).^2 .+
                                  select(dataobject.data, :vy).^2 .+
                                  select(dataobject.data, :vz).^2 ) .* selected_unit .^2

       elseif i == :vϕ_cylinder
           radius = getvar(dataobject, :r_cylinder, center=center)
           x = getvar(dataobject, :x, center=center)
           y = getvar(dataobject, :y, center=center)
           vx = getvar(dataobject, :vx)
           vy = getvar(dataobject, :vy)

           selected_unit = getunit(dataobject, :vϕ, vars, units)
           vars_dict[:vϕ_cylinder] =  (x .* vy .- y .* vx) ./ radius .* selected_unit

       elseif i == :vϕ_cylinder2
           radius = getvar(dataobject, :r_cylinder, center=center)
           x = getvar(dataobject, :x, center=center)
           y = getvar(dataobject, :y, center=center)
           vx = getvar(dataobject, :vx)
           vy = getvar(dataobject, :vy)

           selected_unit = getunit(dataobject, :vϕ2, vars, units)
           vars_dict[:vϕ_cylinder2] =  ((x .* vy .- y .* vx) ./ radius .* selected_unit ).^2

       elseif i == :vr_cylinder
           radius = getvar(dataobject, :r_cylinder, center=center)
           x = getvar(dataobject, :x, center=center)
           y = getvar(dataobject, :y, center=center)
           vx = getvar(dataobject, :vx)
           vy = getvar(dataobject, :vy)

           selected_unit = getunit(dataobject, :vr_cylinder, vars, units)
           vars_dict[:vr_cylinder] =  (x .* vx .+ y .* vy) ./ radius .* selected_unit

       elseif i == :vz2
           vz = getvar(dataobject, :vz)

           selected_unit = getunit(dataobject, :vz2, vars, units)
           vars_dict[:vz2] =  (vz .* selected_unit ).^2

       elseif i == :vr_cylinder2
           radius = getvar(dataobject, :r_cylinder, center=center)
           x = getvar(dataobject, :x, center=center)
           y = getvar(dataobject, :y, center=center)
           vx = getvar(dataobject, :vx)
           vy = getvar(dataobject, :vy)
           #end
           selected_unit = getunit(dataobject, :vr_cylinder2, vars, units)
           vars_dict[:vr_cylinder2] =  ((x .* vx .+ y .* vy) ./ radius .* selected_unit ).^2

       elseif i == :r_cylinder
           selected_unit = getunit(dataobject, :r_cylinder, vars, units)

           vars_dict[:r_cylinder] = select( dataobject.data, (apos, bpos)=>p->
                                               selected_unit * sqrt( (p[apos] - center[1] * boxlen )^2 +
                                                                  (p[bpos] - center[2] * boxlen )^2 ) )

       elseif i == :r_sphere
           selected_unit = getunit(dataobject, :r_sphere, vars, units)
           vars_dict[:r_sphere] = select( dataobject.data, (apos, bpos, cpos)=>p->
                                          selected_unit * sqrt( (p[apos] - center[1] * boxlen )^2  +
                                                                  (p[bpos] - center[2] * boxlen )^2 +
                                                                  (p[cpos] - center[3] * boxlen )^2 )  )

        elseif i == :ekin
            selected_unit = getunit(dataobject, :ekin, vars, units)
            vars_dict[:ekin] =   0.5 .* getvar(dataobject, :mass)  .*
                                (select(dataobject.data, :vx).^2 .+
                                select(dataobject.data, :vy).^2 .+
                                select(dataobject.data, :vz).^2 ) .* selected_unit

        elseif i == :age
            selected_unit = getunit(dataobject, :age, vars, units)
            vars_dict[:age] = ( ref_time .- getvar(dataobject, :birth) ) .* selected_unit
        end

    end

    for i in keys(vars_dict)
        vars_dict[i][isnan.(vars_dict[i])] .= 0
    end

    if length(mask) > 1
        for i in keys(vars_dict)
            vars_dict[i]=vars_dict[i][mask]
        end
    end

    if length(vars)==1
            return vars_dict[vars[1]]
    else
            return vars_dict
    end

end


text:
############## param types

abstract type AccParams end

struct AccParamsSingle{V, W} <: AccParams
    fraglength::Int
    incfun::V
    dataentry::Int
    countfun::W
end

struct AccParamsPair{V, W} <: AccParams
    minfragsize::Int
    maxfragsize::Int
    incfun::V
    dataentry::Int
    countfun::W
end

### helper functions

@inline inc_meta_frag!(V, mpl, mpr, w) =  V[intersect(mpl:mpr, 1:size(V, 1))] .+= w
@inline inc_meta_cut!(V, mpl, mpr, w)  =  (inc_meta_index!(V, mpl, w) ; inc_meta_index!(V, mpr, w))
@inline inc_meta_mid!(V, mpl, mpr, w)  =  inc_meta_index!(V, div(mpl + mpr, 2), w)
@inline inc_meta_index!(V, i::Int, w)  =  (1 ≤ i ≤ size(V, 1)) && (V[i] += w)

@inline inc_heat_frag!(V, k, mpl, mpr, w) =  V[intersect(mpl:mpr, 1:size(V, 1)), k] .+= w
@inline inc_heat_cut!(V, k, mpl, mpr, w)  =  (inc_heat_index!(V, k, mpl, w) ; inc_heat_index!(V, k, mpr, w))
@inline inc_heat_left_cut!(V, k, mpl, mpr, w)  =  inc_heat_index!(V, k, mpl, w)
@inline inc_heat_mid!(V, k, mpl, mpr, w)  =  inc_heat_index!(V, k, div(mpl + mpr, 2), w)
@inline @inbounds inc_heat_index!(V, k, i::Int, w)  =  (1 ≤ i ≤ size(V, 1)) && (V[i, k] += w)

inc_meta_mid_width(bw) = (V, mpl, mpr, w) -> begin
    mid = div(mpl + mpr, 2)
    inc_meta_frag!(V, mid - bw, mid + bw, w)
end

inc_heat_mid_width(bw) = (V, k, mpl, mpr, w) -> begin
    mid = div(mpl + mpr, 2)
    bc = intersect((mid - bw):(mid + bw), 1:size(V, 1))
    V[bc, k] .+= w
end

inc_heat_atac_cut_width(bw) = (V, k, mpl, mpr, w, ao=4) -> begin
    inc_heat_frag!(V, k, mpl + ao - bw, mpl + ao + bw, w)
    inc_heat_frag!(V, k, mpr - ao - bw, mpr - ao + bw, w)
end

inc_meta_atac_cut_width(bw) = (V, mpl, mpr, w, ao=4) -> begin
    inc_meta_frag!(V, mpl + ao - bw, mpl + ao + bw, w)
    inc_meta_frag!(V, mpr - ao - bw, mpr - ao + bw, w)
end

@inline inc_heat_atac_cut!(V, k, mpl, mpr, w, ao=4)   =  (inc_heat_index!(V, k, mpl + ao, w) ; inc_heat_index!(V, k, mpr - ao, w))
@inline inc_heat_atac_frag!(V, k, mpl, mpr, w, ao=4)  =  inc_heat_frag!(V, k, mpl+ao, mpr-ao, w);
@inline inc_meta_atac_frag!(V, mpl, mpr, w, ao=4)     =  inc_meta_frag!(V, mpl+ao, mpr-ao, w);

###############################
function readsintersecting(chroms, locations, FM, dataentry=4, fn=identity)
    total = zeros(Int, length(chroms))
    for (k, (c, l)) in enumerate(zip(chroms, locations))

        v = get_frags(c, l, FM)
        for i = 1:size(v, 2)
            total[k] += fn(v[dataentry, i])
        end
    end
    total
end

function fragregion(chrom, location, strand, FM::FragMatrixSingle{T}, inc_fun=inc_meta_frag!, fraglength=120, dataentry=4, fn=identity; pos=true, neg=true) where {T}
    fw = length(location)
    P = zeros(fw)
    V = get_frags(chrom, location, FM)

    for i = 1:size(V, 2)
        rstrand, _ = GenomeFragments.get_strand_chrom_enc(V[3, i])
        !pos && (rstrand == STRAND_POS) && continue
        !neg && (rstrand == STRAND_NEG) && continue
        if rstrand == STRAND_POS
            s = V[1, i] - location.start + 1
            e = s + fraglength - 1
        else
            e = V[2, i]  - location.start + 1
            s = e - fraglength + 1
        end
        if strand == "-"
            s, e = fw - e + 1, fw - s + 1
        end
        w = fn(V[dataentry, i])
        inc_fun(P, s, e,  w)
    end
    P
end

@inline function fragregion(chrom, location, FM::FragMatrixPair{T}, inc_fun=inc_meta_frag!, minfrag=0, maxfrag=1000, dataentry=4, fn=identity) where {T}
    fw = length(location)
    P = zeros(fw)
    V = get_frags(chrom, location, FM)

    for col in eachcol(V)
        fp = col[2] - col[1] + 1
        s = col[1] - location.start + 1
        e = col[2] - location.start + 1
        w = fn(col[dataentry])
        (minfrag ≤ fp ≤ maxfrag) && inc_fun(P, s, e, w)
    end

    P
end

#################### fragheatmap

function fragheatmap(chroms, locations, strands, FM::Vector{FragMatrixSingle{T}}, inc_fun=inc_heat_mid!, fraglength=120, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, length(chroms))
        fragheatmap!(TH, p, chroms, locations, strands, F, inc_fun, fraglength, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
	    elseif norm_type == :frag_region
	        H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        H .*= norm_scale/sum(F.totalecfrags[data_entry - 3] for F in FM)
    else
        H ./= length(FM)
    end
    H
end

function libsize(FM::FragMatrixPair, dataentry, fn=identity, minfrag=0, maxfrag=1000000)

    σ = 0.0
    for i = 1:size(FM.FM, 2)
        fp = FM.FM[2, i] - FM.FM[1, i] + 1
        σ += ifelse(minfrag ≤ fp ≤ maxfrag, fn(FM.FM[dataentry, i]), 0.0)
    end
    σ
end

function libsize(FM::FragMatrixSingle, dataentry, fn=identity)
    σ = 0
    @inbounds for i = 1:size(FM.FM, 2)
        σ += fn(FM.FM[dataentry, i])
    end
    σ
end

function fragheatmap(chroms, locations, strands, FM::Vector{FragMatrixPair{T}}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; norm_scale=1e+6, norm_type=:mpm, show_progress=true) where T

    fw = length(locations[1])
    H = zeros(Float64, fw, length(chroms))
    p = ProgressMeter.Progress(length(chroms)*length(FM))
    for F in FM
        TH = zeros(Float64, fw, length(chroms))
        fragheatmap!(TH, p, chroms, locations, strands, F, inc_fun, minfragsize, maxfragsize, data_entry, fn, show_progress=show_progress)
        if norm_type == :mpm
            H .+= norm_scale.*TH./F.totalecfrags[data_entry - 3]
        elseif norm_type == :mpm_f
            t = sum(fn, view(F.FM, data_entry, :))
            H .+= norm_scale.*TH./t
        elseif norm_type == :frag
            σ = libsize(F, data_entry, fn, minfragsize, maxfragsize)
            H .+= norm_scale.*TH./σ
	elseif norm_type == :frag_region
	    H .+= norm_scale.*TH/sum(TH)
        elseif norm_type == :total
            H += TH
        else
            H += TH
        end
    end
    if norm_type == :total
        H .*= 1e+6/sum(F.totalecfrags[data_entry - 3] for F in FM)
    else
        H ./= length(FM)
    end
    H
end

function fragheatmap!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, inc_fun=inc_heat_mid!, minfragsize=0, maxfragsize=500, data_entry=4, fn=identity; show_progress=true) where {T}

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < minfragsize) || (fp > maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            ### flip missing?

            w = fn(V[data_entry, i])
            inc_fun(H, k, mpl, mpr,  w)
        end
    end
    H
end

function fragheatmap!(H, p, chroms, locations, strands, FM::FragMatrixSingle{T}, inc_fun=inc_heat_mid!, fraglength=120, data_entry=4, fn=identity; show_progress=true) where {T}

    fw = length(locations[1])

    for (k, (c, l, str)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, (l.start - 2*fraglength):(l.stop + 2*fraglength), FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            strand, _ = GenomeFragments.get_strand_chrom_enc(V[3, i])

            if strand == STRAND_POS
                s = V[1, i] - l.start + 1
                e = s + fraglength - 1
            else
                e = V[2, i]  - l.start + 1
                s = e - fraglength + 1
            end
            if str == "-"
                s, e = fw - e + 1, fw - s + 1
            end
            w = fn(V[data_entry, i])
            inc_fun(H, k, s, e,  w)
        end
    end
    H
end

function fragheatmap!(H, p, chroms, locations, strands, FM::FragMatrixPair{T}, P::AccParamsPair{PV, PW}, show_progress=true) where {T, PV, PW}

    fw = length(locations[1])

    for (k, (c, l, s)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, l, FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            fp = V[2, i] - V[1, i] + 1
            ((fp < P.minfragsize) || (fp > P.maxfragsize)) && continue
            if s == "+"
                mpl = V[1, i] - l.start + 1
                mpr = V[2, i] - l.start + 1
                strand = 1
            else
                mpl = l.stop - V[2, i] + 1
                mpr = l.stop - V[1, i] + 1
                strand = -1
            end

            ### flip missing?

            w = P.countfun(V[data_entry, i])
            P.incfun(H, k, mpl, mpr,  w)
        end
    end
    H
end

function fragheatmap!(H, p, chroms, locations, strands, FM::FragMatrixSingle{T}, P::AccParamsSingle{PV, PW}; show_progress=true) where {T, PV, PW}

    fw = length(locations[1])

    for (k, (c, l, str)) in enumerate(zip(chroms, locations, strands))
        show_progress && next!(p)
        V = get_frags(c, (l.start - 2*P.fraglength):(l.stop + 2*P.fraglength), FM)
        isempty(V) && continue
        for i = 1:size(V, 2)
            strand, _ = GenomeFragments.get_strand_chrom_enc(V[3, i])

            if strand == STRAND_POS
                s = V[1, i] - l.start + 1
                e = s + fraglength - 1
            else
                e = V[2, i]  - l.start + 1
                s = e - fraglength + 1
            end
            if str == "-"
                s, e = fw - e + 1, fw - s + 1
            end
            w = P.countfun(V[data_entry, i])
            P.incfun(H, k, s, e,  w)
        end
    end
    H
end


text:
export barrier_stop
barrier_stop(nothing, from, delta) = from + delta

struct CardinalBarriers
    start::Float64
    len::Float64
end

cardinals() = (Vec2(1,0), Vec2(0,1), Vec2(-1,0), Vec2(0,-1))

function barrier_stop(b::CardinalBarriers, from, delta)
    # https://stackoverflow.com/questions/563198/whats-the-most-efficent-way-to-calculate-where-two-line-segments-intersect
    shortest_u = 1.0+2.0*eps()
    q = from
    s = delta
    for dir in cardinals()
        p = b.start*dir
        r = b.len*dir
        rxs = cross(r,s)
        if rxs == 0.0 # parallel (free pass - this won't force a stop)
            continue
        else
            qmp = q-p
            u = cross(qmp, r)/rxs
            t = cross(qmp, s)/rxs
            if 0.0 <= u < shortest_u && 0.0 <= t <= 1.0
                shortest_u = u
            end
        end
    end
    return from + (shortest_u-2.0*eps())*delta
end


text:
function to_standard_form(x::Matrix{Int})::Matrix{Int}
    n,K = size(x)

    mapx = zeros(Int, K)
    new_x = zeros(Int, n, K)

    gmax = 0

    for i=1:n
        r = findfirst(x[i,:] .== 1)
        if mapx[r] == 0
            gmax += 1
            mapx[r] = gmax
        end
        new_x[i,mapx[r]] = 1
    end

    return new_x
end

function randomAssignments(dataset::Dataset; seed::Union{Nothing, Int}=nothing)::Matrix{Int}
    """ Generates random assignments of n nodes to K groups

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.

    Returns
    -------
    z : Array{Int64,2}
        Matrix of assignments of nodes to groups
    """
    n = dataset.n
    K = dataset.n_communities
    if isnothing(seed)
        Random.seed!(seed)
    end

    assign = rand(1:K, n)

    z = zeros(Int, n, K)
    for i=1:n
        z[i,assign[i]] = 1
    end
    @assert sum(z) == n
    return z
end

function maximizationStep(dataset::Dataset, z::Matrix{Int})::Matrix{Float64}
    """ Calculates the optimal probability matrix given the assignments

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    z : Array{Int64,2}
        Matrix of assignments of nodes to groups

    Returns
    -------
    w : Array{Int64,2}
        Matrix of probabilities
    """
    A = dataset.A
    n = dataset.n
    m = dataset.m
    K = dataset.n_communities
    k = dataset.k
    if ~((n,K) == size(z))
        throw(ArgumentError("The dimensions of the assignment matrix do not match the dimensions of the dataset."))
    end

    w = zeros(K,K)
    for g=1:K, h=1:K
        numerator = 0.5*sum(A[i,j]*z[i,g]*z[j,h] for i=1:n, j=1:n)
        denominator = 0.5*sum( ((k[i]*k[j])/(2*m)) * z[i,g]*z[j,h] for i=1:n, j=1:n)
        if numerator == 0
            w[g,h] = 0
        else
            w[g,h] = numerator / denominator
        end
    end
    return w
end

function expectationStep(dataset::Dataset,
                         w::Matrix{Float64},
                         z::Matrix{Int},
                         variant::String,
                         time_limit::Float64)::Matrix{Int}
    if variant in ["ls1","ls2"]
        return expectationStep_LS(dataset, w, z, variant, time_limit)
    elseif variant == "exact"
        return expectationStep_exact(dataset, w, z, time_limit)
    end
end

function expectationStep_LS(dataset::Dataset,
                            w::Matrix{Float64},
                            z::Matrix{Int},
                            variant::String,
                            time_limit::Float64)::Matrix{Int}

    # Check input argument 'e_step'
    if ~(variant in ["ls1","ls2"])
        throw(ArgumentError(string("Invalid variant for E-step: ", variant)))
    end

    n = dataset.n
    K = dataset.n_communities

    start = time(); # Start counting time

    objval = calculateObjective(dataset, w, z)

    improved = true
    while true

        ## Base cases
        # Break if there is no improvement
        if ~improved
            break # Local optimum
        end
        improved = false

        for i in Random.randperm(n)
            for g in Random.randperm(K)
                # Break if time limit is exceeded
                availableTime = time_limit - (time() - start)
                if availableTime <= 0
                    return z
                end

                # Consider z' constructed from z by relocating vertex i to community g
                new_z = copy(z)
                new_z[i,:] = zeros(Int, K)
                new_z[i,g] = 1

                new_w = copy(w)

                # Integrated maximization step
                if variant == "ls2"
                    new_w = maximizationStep(dataset, new_z)
                end

                # Re-calculate objective value
                new_objval = calculateObjective(dataset, new_w, new_z)

                if new_objval < objval
                    w = new_w
                    z = new_z
                    objval = new_objval

                    improved = true
                end
            end
        end
    end

    return z
end

function expectationStep_exact(dataset::Dataset,
                               w::Matrix{Float64},
                               z::Matrix{Int},
                               time_limit::Float64,
                               TOL::Float64=1e-6)::Matrix{Int}

    A = dataset.A
    K = dataset.n_communities
    n = dataset.n
    m = dataset.m
    k = dataset.k

    verbose = true
    symm_break_const = true

    w[w .== 0] .= 1e-12
    W = zeros(n,n,K,K);
    for i=1:n, j=1:n, g=1:K, h=1:K
        if A[i,j] != 0
            W[i,j,g,h] = 0.5*( -A[i,j]*log(w[g,h]) + ((k[i]*k[j])/(2*m))*w[g,h] )
        else
            W[i,j,g,h] = 0.5*((k[i]*k[j])/(2*m))*w[g,h]
        end
    end

    solver = CPLEX.Optimizer
    model = Model(solver)
    set_optimizer_attribute(model, "CPX_PARAM_SCRIND", verbose)
    set_optimizer_attribute(model, "CPX_PARAM_TILIM", time_limit)
    set_optimizer_attribute(model, "CPX_PARAM_EPGAP", TOL)
    set_optimizer_attribute(model, "CPX_PARAM_EPINT", 1e-9) # integrality tolerance
    set_optimizer_attribute(model, "CPX_PARAM_EPRHS", 1e-9) # feasibility tolerance
    set_optimizer_attribute(model, "CPX_PARAM_EPOPT", 1e-9) # optimality tolerance
    set_optimizer_attribute(model, "CPX_PARAM_SCRIND", false)

    # Variables
    @variable(model, z[1:n,1:K], Bin);

    # Objective
    @objective(model, Min, sum(
            (W[i,j,g,h] * z[i,g] * z[j,h])
            for g=1:K, h=1:K, i=1:n, j=1:n));

    # Assignment constraints
    @constraint(model, assign[i = 1:n], sum(z[i,g] for g=1:K) == 1); # assignment constraint

    # Symmetry breaking constraints
    if symm_break_const == true
        @constraint(model, conObj1, z[1,1] == 1); # object 1 must be in cluster 1
        @constraint(model, symBreak[r=2:(K-1), j=r:n], sum(z[i,l] for l=1:(r-1), i=2:(j-1)) - sum(z[j,l] for l=1:r) <= j - 3 )
    end

    # Solve
    start = time()
    optimize!(model)
    status = string(termination_status(model))
    println("status = $status")
    solvetime = time() - start

    obj_lb = objective_bound(model);
    if isnan(obj_lb)
        obj_lb = -Inf
    end
    obj_ub = objective_value(model);
    if isnan(obj_ub)
        obj_ub = Inf
        return z
    end

    ############################
    # Recover variable values
    z = value.(z)
    z = round.(z[:,:])
    z[map(v -> isnan(v) , z)] .= 0 # replace NaN values with 0
    z = Int.(z)

    return z
end

function expectationMaximization(dataset::Dataset,
                                 e_step::String;
                                 time_limit::Float64=600.0,
                                 seed::Union{Nothing, Int}=nothing,
                                 z::Union{Nothing,Matrix{Int}}=nothing)::Tuple{SBM, Matrix{Int}, OptResults}

    # Check input argument 'e_step'
    if ~(e_step in ["ls1","ls2","exact"])
        throw(ArgumentError(string("Invalid argument for E-step: ", e_step)))
    end

    # Start counting time
    start = time();

    # Initialize with random assignments
    if isnothing(z)
        z = randomAssignments(dataset, seed=seed)
    end
    w = maximizationStep(dataset, z)

    # Calculate objective value
    objval = calculateObjective(dataset, w, z)

    status = nothing
    iterations = 0
    improved = true
    while true
        ## Base cases
        # Break if there is no improvement
        if ~improved
            status = "LocalOptimum"
            break
        end
        improved = false

        # Break if time limit is exceeded
        availableTime = time_limit - (time() - start)
        if availableTime <= 0
            status = "TIME_LIMIT"
            break
        end

        # Expectation Step
        z = expectationStep(dataset, w, z, e_step, availableTime)
        z = to_standard_form(z)

        # Maximization Step
        w = maximizationStep(dataset, z)

        # Check if there was an improvement (otherwise it has reached a local optimum)
        new_objval = calculateObjective(dataset, w, z)
        improved = (new_objval < objval) # Check for improvement
        objval = new_objval # Update objective value (it does not deteriorate in the EM)

        # Increment iterations
        iterations += 1
    end
    solvetime = time() - start;

    # SBM
    sbm = SBM(w, "poisson") # assumes a poisson distribution

    # Build OptResults
    obj_lb = -Inf
    obj_ub = objval
    nodecount = nothing
    lazycount = nothing
    num_constraints = nothing
    opt_results = OptResults(obj_lb, obj_ub, status, solvetime, iterations, nodecount, lazycount, num_constraints)
    displayResults(opt_results)
    return sbm, z, opt_results
end

function multiStartEM(dataset::Dataset,
                      e_step::String;
                      time_limit::Float64=600.0,
                      n_trials::Union{Nothing,Int}=nothing)

    # Check input argument 'e_step'
    if ~(e_step in ["ls1","ls2","exact"])
        throw(ArgumentError(string("Invalid argument for E-step: ", e_step)))
    end

    # Start counting time
    start = time();

    best_ub = Inf
    trial = 1
    sbm = nothing
    z = nothing
    status = nothing
    while true
        # Maximum number of trials
        if (~isnothing(n_trials) && trial > n_trials)
            status = "LocalOptimum"
            break
        end

        # Break if time limit is exceeded
        availableTime = time_limit - (time() - start)
        if availableTime <= 0
            status = "TIME_LIMIT"
            break
        end

        sbm_t, z_t, res_t = expectationMaximization(dataset, e_step, time_limit)

        # Check if there was an improvement
        if res_t.UB < best_ub
            best_ub = res_t.UB # Update best solution
            sbm = sbm_t
            z = z_t
        end

        trial += 1
    end

    # Build OptResults
    obj_lb = -Inf
    obj_ub = best_ub
    solvetime = (time() - start)
    nodecount = nothing
    lazycount = nothing
    num_constraints = nothing
    opt_results = OptResults(obj_lb, obj_ub, status, solvetime, trial, nodecount, lazycount, num_constraints)
    displayResults(opt_results)
    return sbm, z, opt_results
end

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
# binary search in interval (low,high)
# assumes that if pred(x)==true, then pred(y)==true ∀y>x
# assumes pred(low)==false and pred(high)==true
function binary_search_interval(pred, low::T, high::T, precision::T) where {T<:AbstractFloat}
    while high-low > precision
        mid = (low+high)/2

        if pred(mid)
            high = mid
        else
            low = mid
        end
    end
    (low+high)/2
end
binary_search_interval(pred, low, high, precision) = binary_search_interval(pred, promote(low, high, precision)...)

# Soft thresholding operator S (See Witten et al.)
function soft_threshold(a,c)
    sign(a) .* max( abs(a)-c, 0 )
end

# computes right singular vector of X
function rightsingularvector(X)
    K = Symmetric(X'X);
    N = size(K,1);
    _,v=eig(K,N:N);
    dropdims(v; dims=2)
end

# computes
#   argmax_u   u'w
#   subject to ||u||_2^2 ≤ 1
# (When used for PCA iterations, w = Xv or w = X'u = (u'X)'.)
function argmax_unpenalized(w)
    # w/norm(w)
    w/max(1e-16,norm(w)) # avoid div by zero
end

# computes
#   argmax_u   uᵀw
#   subject to ||u||_2^2 ≤ 1
#              uᵀO = (0 0 ... 0)          (vector of zeros)
# (When used for PCA iterations, w = Xv or w = X'u = (u'X)' and U is matrix with previous components, i.e. U and V respectively.)
# See Witten et al 2009, section 3.2.
#     uₖ = Cₖ₋₁Cₖ₋₁ᵀXvₖ / ||Cₖ₋₁ᵀXvₖ||₂
#        = (I - Uₖ₋₁Uₖ₋₁ᵀ)Xvₖ / ||Cₖ₋₁ᵀXvₖ||₂
#        = (Xvₖ - Uₖ₋₁Uₖ₋₁ᵀXvₖ) / √( ||Xvₖ||₂² - ||Uₖ₋₁ᵀXvₖ||₂² )
#        = (w - UUᵀw) / √( ||w||₂² - ||Uᵀw||₂² )
#     where C is the orthogonal complement of O (not needed in computations).
function argmax_unpenalized_orthogonal(w, U)
    # Uw = U'w
    # (w-U*Uw) / sqrt(sum(w.^2)-sum(Uw.^2))

    # Simplified
    u = w - U*(U'w)
    u / norm(u)
end

# computes
#   argmax_u   u'a
#   subject to ||u||₂² ≤ 1
#              ||u||₁ ≤ c
# (When used for PCA iterations, a = Xv or a = X'u = (u'X)'.)
function argmax_l1penalized(a,c)
    # test with Δ=0 first (i.e. S(a,0)=a)
    # u = a/norm(a)
    u = a/max(1e-16, norm(a)) # avoid div by zero
    c2 = c*c # we almost only use c squared
    (length(a)<=c2 || norm(u,1)≤c) && return u

    # otherwise chose Δ s.t. ||u||₁=c, where u = S(a,Δ)/||S(a,Δ)||₂

    # we could use a heap for b to avoid sorting all of b (should save a lot for very sparse vectors!)
    heap = binary_maxheap(abs(a)) # work with sorted abs(a) since it doesn't change value of Δ and makes the search simpler

    n1 = 0.0
    n2 = 0.0 # squared
    pk=qk = 0.0
    Δ = 0.0
    b = pop!(heap)

    k = 1
    #while k<length(b)
    while k<length(a)
        # pk += b[k]*b[k]
        # qk += b[k]
        pk += b*b
        qk += b

        bNext = pop!(heap)
        # n1Next = qk - k*b[k+1]
        # n2Next = k*b[k+1]*b[k+1] + pk - 2b[k+1]*qk
        n1Next = qk - k*bNext
        n2Next = k*bNext*bNext + pk - 2bNext*qk

        n1Next*n1Next≥c2*n2Next && break # squared version, ok since everything is positive - now we know bₖ₋₁≥Δ≥bₖ

        n1,n2 = n1Next,n2Next
        b = bNext
        k+=1
    end

    A = (k*k-k*c2)
    B = 2(k*n1 - c2*n1)
    C = n1*n1-c2*n2

    p = B/A
    q = C/A

    #Δ = max(0.0, b[k]-(-p/2+sqrt(max(p*p/4-q,0.0)))) # max with zero only to avoid numerical problems
    Δ = max(0.0, b-(-p/2+sqrt(max(p*p/4-q,0.0)))) # max with zero only to avoid numerical problems

    u=soft_threshold(a,Δ);
    # u/norm(u)
    u/max(1e-16, norm(u)) # avoid div by zero
end

# numerically compute rank 1 approximation
# 2nd component are retrieved using pca_rank(X-u*s*v',...), where (u,s,v) is the result from the first, and so on.
function pcarank1l1(X,c1,c2,tol,maxIters)
    P = size(X,1)
    N = size(X,2)

    v = rightsingularvector(X) # start with optimal v for non-penalized problem.
    u = zeros(P)
    s = 0.0
    converged = false

    # for i=1:maxIters # TODO: convergence criterion
    #     u = c1==Inf ? argmax_unpenalized(X*v) : argmax_l1penalized(X*v,c1)
    #     v = c2==Inf ? argmax_unpenalized(X'u) : argmax_l1penalized(X'u,c2)
    # end
    # s = (u'X*v)[1]

    for i=1:maxIters
        prevU = u
        u = c1==Inf ? argmax_unpenalized(X*v) : argmax_l1penalized(X*v,c1)

        Xu = X'u
        prevV = v
        v = c2==Inf ? argmax_unpenalized(Xu) : argmax_l1penalized(Xu,c2)

        prevS = s
        s = (v'Xu)[1]
        #println("s: ", s)

        #println("i: ", i, ", v: ", norm(v-prevV), ", u: ", norm(u-prevU), ", s:", s-prevS)

        norm(v-prevV)<tol && norm(u-prevU)<tol && (converged=true; break)
        s-prevS < tol && (converged=true; break)

        #norm(v-prevV)<tol && norm(u-prevU)<tol && (converged=true; println(norm(v-prevV), " ", norm(u-prevU)); break)
        #s-prevS < tol && (converged=true; println(s-prevS); break)
    end

    (u,s,v,converged)
end

# numerically compute rank 1 approximation
# 2nd component is retrieved using pcarank1l1_orthogonal(X-u*s*v',c1,V), where (u,s,v) is the result from the first, and so on.
function pcarank1l1_orthogonal(X,c1,V,v0,tol,maxIters)
    P = size(X,1)
    N = size(X,2)

    v = v0

    u = zeros(P)
    s = 0.0
    converged = false

    for i=1:maxIters
        prevU = u
        u = c1==Inf ? argmax_unpenalized(X*v) : argmax_l1penalized(X*v,c1)

        Xu = X'u
        prevV = v
        v = isempty(V) ? argmax_unpenalized(Xu) : argmax_unpenalized_orthogonal(Xu,V)

        prevS = s
        s = (v'Xu)[1]

        norm(v-prevV)<tol && norm(u-prevU)<tol && (converged=true; break)
        s-prevS < tol && (converged=true; break)
    end

    (u,s,v,converged)
end

function pcarank1l1_orthogonal(X,c1,V,tol,maxIters)
    # v0 = rightsingularvector(X) # start with optimal v for non-penalized problem.

    # This enforces orthogonality to V for initial guess.
    v0 = isempty(V) ? rightsingularvector(X) :
                      rightsingularvector( X - (X*V)*V' )

    pcarank1l1_orthogonal(X,c1,V,v0,tol,maxIters)
end

# Sparce PCA with l1-constraints on U-side and orthogonality constraints on V-side.
# If a component with higher optimum than a previous component is found, we backtrack and use it as a starting guess for the previous component.
# This guarantees that S will be sorted and reduces some issues with getting stuck in local optima.
function sparsepcal1(X, k::Integer, c1; backtrack=true, tol=1e-6, maxIters=1000)
    P = size(X,1)
    N = size(X,2)
    @assert N<=P # computations are arranged to exploit N<P (or rather, N<<P)

    # set L1-constraints to Inf if they are not active (better for pcarank1l1)
    c1 = c1<√P ? c1 : Inf

    U = zeros(P,k)
    S = zeros(k)
    V = zeros(N,k)

    XOrig = backtrack ? copy(X) : similar(X,0,0) # only necessary to keep track of original matrix if we are backtracking

    v0 = rightsingularvector(X)

    i=1
    while i<=k
        U[:,i], S[i], V[:,i], converged = pcarank1l1_orthogonal(X,c1,V[:,1:i-1],v0,tol,maxIters)
        converged || @warn("sparsepca component $i did not converge")

        if S[i]<0 || any(abs(V[:,i]'V[:,1:i-1]) .> 1e-6) # In this case, we are not finding any more components. I.e. the rank(XOrig)<i.
            U[:,i], S[i], V[:,i] = .0, .0, .0
            break
        end

        if backtrack && i>1 && S[i]-S[i-1]>tol # Backtrack if we found a component with higher S[i] than the previous one.
            v0 = V[:,i] # Use this optimum as a starting point for the search on the previous matrix. The orthogonality of samples (V) guarantees that the starting point will have the same S value, i.e. we will find a better local optima.
            i = something(findfirst(S[i]-S.>tol),0) - 1 # How far do we need to backtrack?
            X = XOrig - U[:,1:i]*(diagm(S[1:i])*V[:,1:i]') # Restore matrix state s.t. only the i first components have been removed.
        else
            X = X - U[:,i]*(S[i]*V[:,i]') # Remove this component to be able to find next component.
            v0 = rightsingularvector( X - (X*V)*V' ) # Starting guess orthogonal to V.
        end
        i += 1
    end

    # change signs for stability, favouring U's that are positive on average
    sgn = (sum(U;dims=1).>=0)*2.0-1.0
    U .*= sgn
    V .*= sgn

    (U,S,V)
end


text:
module Justly

import AudioSchedules: AudioSchedule
using AudioSchedules:
    add!,
    Cycles,
    fill_all_task_ios,
    Hook,
    Line,
    Map,
    SawTooth,
    Scale,
    triples,
    Weaver,
    write_buffer,
    write_series!
import Base: getproperty, parse, print, Rational, setproperty!, show
using Base.Meta: ParseError
using Base.Threads: nthreads, @spawn
using LightXML:
    add_text, create_root, find_element, new_child, root, set_attribute, XMLDocument
using PortAudio: PortAudioStream
using QML:
    exec,
    ListModel,
    addrole,
    loadqml,
    # to avoid QML bug
    QML,
    qmlfunction,
    setconstructor
using Qt5QuickControls2_jll: Qt5QuickControls2_jll
using SampledSignals: samplerate
using Test: @test
using Unitful: Hz, ms, s, minute, uconvert

using YAML: load_file, write_file
import YAML: _print

# reexport to avoid a QML bug
export QML

const TIME = typeof(0.0s)
const FREQUENCY = typeof(0.0Hz)

function property_model(a_vector, property_names)
    list_model = ListModel(a_vector, false)
    for property_name in property_names
        addrole(
            list_model,
            String(property_name),
            let property_name = property_name
                item -> getproperty(item, property_name)
            end,
            let property_name = property_name, stdout = stdout
                (items, value, index) -> try
                    setproperty!(items[index], property_name, value)
                catch an_error
                    # errors will be ignored, so print them at least
                    showerror(stdout, an_error)
                end
            end,
        )
    end
    setconstructor(list_model, eltype(a_vector))
    list_model
end

"""
    pedal(duration; slope = 1 / 0.1s, peak = 1, overlap = 1/2)

You can use `pedal` to make an envelope with a sustain and ramps at the beginning and end.
`overlap` is the proportion of the ramps that overlap.
"""
function pedal(duration; slope = 1 / 0.1s, peak = 1, overlap = 1 / 2)
    ramp = peak / slope
    ramp_overlap = ramp * overlap
    # if it's too short, there won't be a plateau
    middle = duration - ramp - ramp + ramp_overlap
    if middle <= 0s
        half_duration = (duration + ramp_overlap) / 2
        short_peak = half_duration * slope
        (0, Line => half_duration, short_peak, Line => half_duration, 0)
    else
        (0, Line => ramp, peak, Line => middle, peak, Line => ramp, 0)
    end
end
export pedal

function throw_parse_error(text, description, line_number)
    throw(ParseError(string("Can't parse $text as $description on line $line_number")))
end

include("Interval.jl")
include("Note.jl")
include("Chord.jl")
include("Song.jl")
include("AudioSchedule.jl")
include("MusicXML.jl")

function precompile_song(task_ios, song, buffer)
    for (series, _) in AudioSchedule(song)
        write_series!(task_ios, series, 1, buffer, 0)
    end
end

function get_dummy_envelope(song, frequency)
    ramp = song.ramp
    sample_rate = song.sample_rate
    map(
        (((wave, _, duration),) -> (wave, round(Int, duration * sample_rate))),
        triples(
            sample_rate,
            Map(Scale(song.volume), Map(song.wave, Cycles(frequency))),
            0s,
            0,
            Line => ramp,
            1,
            Line => 0.5s,
            1,
            Line => ramp,
            0,
        ),
    )
end

# cumulative product of previous modulations
function update_key(song, chord_index)
    key = song.initial_key
    for chord in view(song.chords, 1:(chord_index + 1))
        key = key * Rational(chord.modulation.interval)
    end
    key
end

function press!(task_ios, song, presses, releases, buffer)
    for (chord_index, voice_index) in presses
        buffer_at = 0
        if voice_index < 0
            precompile_song(task_ios, song, buffer)
            Base.GC.enable(false)
            for (series, series_total) in AudioSchedule(
                song;
                chords = (@view song.chords[(chord_index + 1):end]),
                initial_key = update_key(song, chord_index - 1),
            )
                if isready(releases)
                    break
                end
                buffer_at = write_series!(task_ios, series, series_total, buffer, buffer_at)
            end
            write_buffer(buffer, buffer_at)
            Base.GC.enable(true)
            take!(releases)
        else
            Base.GC.enable(false)
            # all three will be pairs of iterators and number of frames
            (ramp_up, sustain, ramp_down) = get_dummy_envelope(
                song,
                update_key(song, chord_index) *
                Rational(song.chords[chord_index + 1].notes[voice_index + 1].interval),
            )
            buffer_at = write_series!(task_ios, ramp_up..., buffer, buffer_at)
            buffer_at = write_series!(task_ios, sustain..., buffer, buffer_at)
            buffer_at = write_series!(task_ios, ramp_down..., buffer, buffer_at)
            write_buffer(buffer, buffer_at)
            Base.GC.enable(true)
        end
    end
end

function get_midi_name(midi_code)
    octave, degree = get_note_parts(midi_code)
    string("Initial key: ", NOTE_NAMES[degree + 1], "<sub>", octave, "</sub>")
end

"""
    function edit_song(song_file;
        ramp = 0.1s,
        number_of_tasks = nthreads() - 2,
        test = false,
        keyword_arguments...
    )

Use to edit songs interactively.
The interface might be slow at first while Julia is compiling.

- `song_file` is a YAML string or a vector of [`Chord`](@refs)s. Will be created if it doesn't exist.
- `number_of_tasks` is the number of tasks to use to process data. Defaults to 2 less than the number of threads; we need 1 master thread for QML and 1 master thread for AudioSchedules.
- If `test` is true, will open the editor briefly to test it.
- `keyword_arguments` will be passed to [`read_justly`](@ref).

For more information, see the `README`.

```jldoctest
julia> using Justly

julia> edit_song(joinpath(pkgdir(Justly), "test", "song.justly"); test = true)
```
"""
function edit_song(
    song_file;
    number_of_tasks = nthreads() - 2,
    test = false,
    keyword_arguments...,
)
    if nthreads() < 3
        error("Justly needs at least 3 threads to function")
    end

    song = if isfile(song_file)
        read_justly(song_file; keyword_arguments...)
    else
        Song(; keyword_arguments...)
    end

    qmlfunction("get_midi_name", get_midi_name)

    qmlfunction("get_initial_midi_code", let song = song
        () -> get_initial_midi_code(song)
    end)

    qmlfunction(
        "update_initial_midi_code",
        let song = song
            midi_code -> update_initial_midi_code!(song, midi_code)
        end,
    )

    qmlfunction("get_beats_per_minute", let song = song
        () -> get_beats_per_minute(song)
    end)

    qmlfunction(
        "update_beats_per_minute",
        let song = song
            beats_per_minute -> update_beats_per_minute!(song, beats_per_minute)
        end,
    )

    presses = Channel{Tuple{Int, Int}}(0)
    qmlfunction(
        "press",
        let presses = presses
            (chord_index, voice_index) -> put!(presses, (chord_index, voice_index))
        end,
    )

    releases = Channel{Nothing}(0)
    qmlfunction("release", let releases = releases
        () -> put!(releases, nothing)
    end)

    print_function = let song = song
        io -> print(io, song)
    end

    qmlfunction("to_yaml", let song_file = song_file, print_function = print_function
        () -> open(print_function, song_file, write = true)
    end)

    stream = PortAudioStream(0, 1, writer = Weaver(); warn_xruns = false)
    buffer = stream.sink_messanger.buffer
    task_ios = fill_all_task_ios(buffer; number_of_tasks = number_of_tasks)
    (ramp_up, sustain, ramp_down) = get_dummy_envelope(song, 440.0Hz)
    # each is a pair of waves and number of samples
    # set the first index of the buffer to match each wave
    write_series!(task_ios, ramp_up[1], 1, buffer, 0)
    write_series!(task_ios, sustain[1], 1, buffer, 0)
    write_series!(task_ios, ramp_down[1], 1, buffer, 0)
    precompile_song(task_ios, song, buffer)
    loadqml(
        joinpath(@__DIR__, "Song.qml");
        chords_model = property_model(
            song.chords,
            (:numerator, :denominator, :octave, :beats, :words, :notes_model),
        ),
        test = test,
    )
    press_task = Task(
        let task_ios = task_ios,
            song = song,
            presses = presses,
            releases = releases,
            buffer = buffer

            () -> press!(task_ios, song, presses, releases, buffer)
        end,
    )
    press_task.sticky = false
    schedule(press_task)
    try
        exec()
        if test
            # note: this is 1, 1 in julia
            put!(presses, (0, 0))
            put!(presses, (0, -1))
            put!(releases, nothing)
        end
    catch an_error
        # can't error while QML is running, so just message
        println("Front-end errored:")
        showerror(stdout, an_error)
    finally
        close(presses)
        try
            wait(press_task)
        catch an_error
            # can't error while QML is running, so just message
            println("Back-end errored:")
            showerror(stdout, an_error)
        finally
            close(releases)
            close(stream)
        end
    end
    nothing
end
export edit_song

end

max_stars_repo_name: bramtayl/Justly

max_issues_repo_name: bramtayl/Justly

max_forks_repo_name: bramtayl/Justly


text:
using Plots, Test
using RecipesPipeline

@testset "plot" begin
    pl = plot(1:5)
    pl2 = plot(pl, tex_output_standalone = true)
    @test pl[:tex_output_standalone] == false
    @test pl2[:tex_output_standalone] == true
    plot!(pl, tex_output_standalone = true)
    @test pl[:tex_output_standalone] == true
end

@testset "get_axis_limits" begin
    x = [0.1, 5]
    p1 = plot(x, [5, 0.1], yscale = :log10)
    p2 = plot!(identity)
    @test all(RecipesPipeline.get_axis_limits(p1, :x) .== x)
    @test all(RecipesPipeline.get_axis_limits(p2, :x) .== x)
end

@testset "Slicing" begin

    @test plot(1:5, fillrange = 0)[1][1][:fillrange] == 0
    data4 = rand(4,4)
    for i in axes(data4, 1)
        @test plot(data4, fillrange = 0)[1][i][:fillrange] == 0
        @test plot(data4, fillrange = [1, 2])[1][i][:fillrange] == [1.0, 2.0]
        @test plot(data4, fillrange = [1 2])[1][i][:fillrange] == (iseven(i) ? 2 : 1)
    end
end


text:
prepcomp(;refreshcomp=true) = loadcompdata(refreshcomp=refreshcomp)

function loadcompdata(;companame=COMP_A_NAME, compqname = COMP_Q_NAME, compname = COMP_NAME,
   datapath=DATA_PATH, refreshcomp=true,
   incsvstream::Function = IN_CSV_STREAM,
   injlsstream::Function = IN_JLS_STREAM,
   outjlsstream::Function = OUT_JLS_STREAM)::DataFrame

  local compa::DataFrame
  local compq::DataFrame
  local comp::DataFrame

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcomp || (!isfile("$datapath\\$compname.jls.lz4"))
    compa = incsvstream("$datapath\\$companame.csv.lz4") |> CSV.read |> DataFrame
    compq = incsvstream("$datapath\\$compqname.csv.lz4") |> CSV.read |> DataFrame

    comp = preprocesscomp(compa, compq)
    comp = mergecompccm!(comp)

    outjlsstream("$datapath\\$compname.jls.lz4") do s
      serialize(s, comp) #NOTE: DEBUG ONLY
    end

    #NOTE: FINAL VERSION
    #preprocesscomp(compa, compq)
    #serialize("$datapath\\$compname.jls", comp)
  else
    injlsstream("$datapath\\$compname.jls.lz4") do s
      comp = deserialize(s) #NOTE: DEBUG ONLY
    end
  end

  @info "Compustat and CCM linking data loaded and saved into file $compname.jls.lz4"

  return comp
end
#helper function to clean up names in the quarterly file
function properqname(s::Symbol)
  str::String = string(s)

  if str[end] == 'q' || (length(str) ≥ 3 && str[(end-2):end] == "qtr")
    return s
  end

  s = Symbol(s,:q)
  return s
end

#helper function to convert a fiscal year-quarter to a date
@inline function fyearqtr2date(yq::YearQuarter, fyr::Int)
  ym::YearMonth = YearMonth(yq.y, fyr)
  offsetmonths::Int = -(4 - yq.q)*3
  ym += offsetmonths
  return eom(ym)
end

#this processes both the q and the a comp frames, with as much commonly used code as possible
function preprocesscomp(compa::DataFrame, compq::DataFrame;
  compdateformat::DateFormat = COMP_DATE_FORMAT,
  retainedcolscomp::Vector{Symbol} = RETAINED_COLS_COMP,
  testoutput::Bool = TEST_OUTPUT)

  #do operations that are common to both the annual and the quarterly data
  compa = preprocesscomp!(compa)
  compq = preprocesscomp!(compq)

  compa = preprocesscompa!(compa)
  compq = preprocesscompq!(compq)

  comp::DataFrame = join(
    compa, compq, on=[:gvkey=>:gvkeyq, :fyear=>:fyearq], kind=:inner)

  #NOTE: The below is not actually done in the paper, but I think its important
  #(The paper only uses the annual data)
  comp.iaq = ((iaq,ia)->ismissing(iaq) ? ia : iaq).(comp.iaq, comp.ia)

  #form the initial interval candidates
  #these are accurate if the link doesn't changes (adjust for this case later)
  comp = formannouncementintervals!(comp)
  comp = select!(comp,retainedcolscomp)

  #println(describe(comp))
  if testoutput
    d = Date(2015,12,31)
    testoutput && CSV.write("output\\comp_$d.csv", comp[((bd,ed)->
      (d > bd) && (d ≤ ed)).(comp.begindate, comp.enddate), :])
  end

  return comp
end

#operations that are common to both the annual and the quarterly data
function preprocesscomp!(comp::DataFrame)
  #make all column names lower case
  names!(comp, (s::Symbol->Symbol(lowercase(string(s)))).(names(comp)))

  #format dates
  comp.ddate = (s->parsecomp(Date, s)).(comp.datadate)
  comp.dyear = (d::MDate-> ismissing(d) ? missing : Dates.year(d)).(comp.ddate)

  comp = comp[(!ismissing).(comp.dyear), :] #require the calendar date

  return comp
end

function preprocesscompa!(compa::DataFrame;
  usedcompafields::Vector{Symbol} = USED_COMPA_VALUE_FIELDS)
  compa = compa[(!ismissing).(compa.fyear), :]

  #NOTE: REMINDER- convert to first day of month for date math
  #Date math works better using the first day of the month,
  #even though the last day of the month is technically correct
  compa.fdate = ((y::Int,m::Int)->
    lastdayofmonth(Date(y,m,1))).(compa.fyear, compa.fyr)

  #delete rows with missing data
  compa = compa[(!ismissing).(compa.at), :]

  #now process cases where firms have duplicate date-firm key pairs
  #the date field is used to affirm the latest record
  #(for annual records, that is the annual date,
  #for quarterly its the announcement date)
  compa = dedupcomp!(compa, [:gvkey, :fyear], usedcompafields, :indfmt, :ddate)

  #sort the data
  sort!(compa, [:gvkey, :fyear])

  ####now make the I/A amount, both in a quarterly and an annual form
  compa.ia = Vector{MFloat64}(undef, size(compa,1))

  lagwithin!(compa, :at, :gvkey, :fyear, sorted=true)
  compa.ia .= (compa.at .- compa.Lat) ./ compa.Lat
  compa.ia .= (f::MFloat64 -> (!ismissing(f)) && isfinite(f) ? f : missing).(compa.ia)

  return compa
end

function preprocesscompq!(compq::DataFrame;
  usedcompqfields::Vector{Symbol} = USED_COMPQ_VALUE_FIELDS)
  #use a consistent naming convention
  names!(compq, (properqname).(names(compq)))

  compq = compq[((cq, fq)->(!ismissing(cq)) && (!ismissing(fq))).(
    compq.datacqtr, compq.datafqtr),:]
  compq.fyearqtr =  (s->parsecomp(YearQuarter, s)).(compq.datafqtr)
  compq.cyearqtr =  (s->parsecomp(YearQuarter, s)).(compq.datacqtr)
  compq.adateq = (s->parsecomp(Date, s)).(compq.rdq)
  compq = compq[(!ismissing).(compq.adateq),:]

  #NOTE: REMINDER- convert to first day of month for date math
  #Date math works better using the first day of the month,
  #even though the last day of the month is technically correct
  compq.fdateq = ((yq::YearQuarter,fyr::Int)->
    fyearqtr2date(yq,fyr)).(compq.fyearqtr, compq.fyrq)

  #don't let the fiscal date be too far from the data date (uncomment to use data date)
  compq = compq[(r->(r.adateq - Month(6) ≤ r.fdateq)).(eachrow(compq)),:]

  #under no circumstances allow restatements past the subsequent quarter's announcement date
  validatefadates!(compq; Fgroup=:gvkeyq, Ffdate=:fyearqtr,
    Fadate=:adateq)

  #now process cases where firms have duplicate date-firm key pairs
  #the date field is used to affirm the latest record
  #(for annual records, that is the annual date,
  #for quarterly its the announcement date)
  #Require a 1:1 relationship between announcement dates and year-quarters
  #Mainly want to avoid interm information disclosures
  #NOTE: could revisit in robustness checks
  compq = dedupcomp!(compq, [:gvkeyq, :fyearqtr], usedcompqfields, :indfmtq, :adateq)
  compq = dedupcomp!(compq, [:gvkeyq, :adateq], usedcompqfields, :indfmtq, :fyearqtr)

  #sort the data
  sort!(compq, [:gvkeyq, :fyearqtr])

  ####now make the I/A amount, both in a quarterly and an annual form
  compq.iaq = Vector{MFloat64}(undef, size(compq,1))

  lagwithin!(compq, :atq, :gvkeyq, :fyearqtr, sorted=true, lags=4)
  compq.iaq .= (compq.atq .- compq.L4atq) ./ compq.L4atq
  compq.iaq .= (f::MFloat64 -> (!ismissing(f)) && isfinite(f) ? f : missing).(compq.iaq)

  ###derive ROE
  #first get shareholders equity
  compq.seq = (r::DataFrameRow-> coalesce(r.seqq, r.ceqq+r.pstkq, r.atq - r.ltq)).(eachrow(compq))
  #now add in deferred taxes and inv credits to get book equity, and subtract out preferred stock
  compq.beq = ((r::DataFrameRow)->
    r.seq + coalesce(r.txditcq,0.0) - coalesce(r.pstkrq, r.pstkq,0.0)).(eachrow(compq))
  lagwithin!(compq, :beq, :gvkeyq, :fyearqtr, sorted=true)
  compq.roeq = compq.ibq ./ compq.Lbeq
  compq.roeq .= (f::MFloat64->(!ismissing(f)) && isfinite(f) ? f : missing).(compq.roeq)

  return compq
end

  #under no circumstances allow restatements past the subsequent quarter's announcement date
function validatefadates!(comp::DataFrame; Fgroup::NSymbol=nothing, Ffdate::NSymbol=nothing,
  Fadate::NSymbol=nothing, sort::Bool = true)::Nothing

  (isnothing(Fgroup)) && error("Fgroup missing from validatefadates")
  (isnothing(Ffdate)) && error("Ffdate missing from validatefadates")
  (isnothing(Fadate)) && error("Fadate missing from validatefadates")

  sort!(comp, [Fgroup, Fadate])
  T::Type = eltype(comp[!,Ffdate])
  comp.maxfdate_at_adate = Vector{Union{T,Missing}}(undef, size(comp,1))

  compgrp::GroupedDataFrame = groupby(comp, Fgroup)
  @mpar for i ∈ 1:length(compgrp)
    subdf::SubDataFrame = compgrp[i]
    maxfdate::T = subdf[1,Ffdate]
    for r ∈ eachrow(subdf)
      (r[Ffdate] > maxfdate) &&  (maxfdate = r[Ffdate])
      r.maxfdate_at_adate = maxfdate
    end
  end

  filter!(comp) do r::DataFrameRow
    r.maxfdate_at_adate ≤ r[Ffdate]
  end

  return nothing
end

function dedupcomp!(comp::DataFrame,
  groupbyfields::Vector{Symbol}, usedcompfields::Vector{Symbol},
  formatfield::Symbol, datefield::Symbol)

  #keepflags
  comp.keep = trues(size(comp,1))

  #use this for deduping
  comp.valuesmissings = Vector{MFloat64}(undef, size(comp,1))

  ncols::Int = size(comp,2)
  scomps::GroupedDataFrame = groupby(comp, groupbyfields)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]

    if size(scomp,1) > 1

      #keep the most recent version
      maxdate::Union{Int, YearQuarter, Date} = maximum(scomp[!,datefield])
      mindate::Union{Int, YearQuarter, Date} = minimum(scomp[!,datefield])
      #println("got here")

      #why could this happen? not clear to me, but it does happen a bit
      if maxdate ≠ mindate
        scomp.keep .= false
      end

      #special case where we have two records of two different formats (indl and fs)
      #in this case, keep the format with the least missing values
      if sum(scomp.keep) == 2
        sscomp::SubDataFrame = view(scomp, scomp.keep, :)
        if (length(unique(sscomp.indfmt)) == 2)

          sscomp.valuesmissings[1] = sum((c->ismissing(sscomp[1,c])).(usedcompfields))
          sscomp.valuesmissings[2] = sum((c->ismissing(sscomp[2,c])).(usedcompfields))

          sscomp.keep .= (sscomp.valuesmissings .== minimum(sscomp.valuesmissings))

          #tie breaker is to keep the INDL record (industrial since we are looking at non-financials)
          if sum(sscomp.keep) > 1
            sscomp.keep .= (sscomp[!, formatfield] .== :INDL)
          end
        end
      end

      #if we still have a dup, delete all the records
      if sum(scomp.keep) > 1
        @warn "dedup failed group fields:
          $groupbyfields with values $((f->scomp[1,f]).(groupbyfields)) Records removed."
        scomp.keep .= false
      end
    end
  end

  comp = comp[comp.keep, :] #bandaid until filter is fixed

  select!(comp, Not([:keep, :valuesmissings]))

  return comp
end

#NOTE: messy function
#basically to merge in crsp, the crsp day or month must occur
#in the intersection of the link dates, the period six months following
#the announcement date, and the next announcement date
#this handles everything but the intersection with the link dates
function formannouncementintervals!(comp::DataFrame;
    lastdate::Date = LAST_DATE, time2stale::T = TIME_2_STALE) where T<:DatePeriod

  sort!(comp, [:gvkey, :adateq])

  comp.Nadateq = Vector{MDate}(undef, size(comp,1))

  scomps = groupby(comp, :gvkey)
  @mpar for i ∈ 1:length(scomps)
    scomp::SubDataFrame = scomps[i]
    Nsub::Int = size(scomp,1)

    for i ∈ 1:(Nsub-1)
      scomp[i,:Nadateq] = scomp[i+1,:adateq]

      #check some things that shouldn't happen
      (scomp[i,:adateq] ≥ scomp[i,:Nadateq]) && @warn (
        "announcement dates out of order.\n $(scomp[i:(i+1),:])")

      (scomp[i,:fyearqtr] ≥ scomp[i+1,:fyearqtr]) && @warn (
        "fiscal year dates out of order.\n $(scomp[i:(i+1),:])")
    end

    scomp[Nsub, :Nadateq] = lastdate
  end

  comp.fdateqmax = (lastdayofmonth).(
    (firstdayofmonth).(comp.fdateq) .+ time2stale)

  comp.begindate = comp.adateq
  comp.enddate = (min).(comp.Nadateq,comp.fdateqmax)

  return comp
end

parsecomp(::Type{Date}, s::String;
    compdateformat::DateFormat = COMP_DATE_FORMAT) = Dates.Date(s, compdateformat)

#helper methods and special cases
parsecomp(::Type{<:Any}, v::Missing) = missing
parsecomp(::Type{Date}, i::Int) = parsecomp(Date, "$i")

#parse a value where the type is known: Generic case
function parsecomp(::Type{T}, s::String;
    ##NOTE: uncomment when we have the missings parsedmissings = COMP_PARSED_MISSINGS
    ) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  #NOTE: uncomment when we have the missings
  #(!ismissing(v)) && (v ∈ parsedmissings) && (v=missing)

  return v
end

function parsecomp(::Type{YearQuarter}, s::String)
  local yq::MYearQuarter

  if length(s) ≠ 6
    yq = missing
  else
    y::MInt = parsecomp(Int, "$(s[1:4])")
    q::MInt = parsecomp(Int, "$(s[6])")
    (q==1) || (q==2) || (q==3) || (q==4) || (q=missing)

    yq = (!ismissing(y)) && (!ismissing(q)) ? YearQuarter(y,q) : missing
  end

  return yq
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using DataFrames, CSV, PlotlyJS, RDatasets
using Dash, DashHtmlComponents, DashCoreComponents

iris = dataset("datasets", "iris")

p1 = Plot(iris, x=:SepalLength, y=:SepalWidth, mode="markers", marker_size=8, group=:Species)

app = dash()

app.layout = html_div() do
    html_h4("Iris Sepal Length vs Sepal Width"),
    dcc_graph(
        id = "example-graph-3",
        figure = p1,
    )
end

run_server(app, "0.0.0.0", debug=true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
"""
SquaredMatrixPlusIModule (C) 2020, Martin Gurtner
Efficient implementation of algebraic operations on matrices in the form D*D'+λ*I
"""
module SquaredMatrixPlusIModule

using LinearAlgebra

export SquaredMatrixPlusI

mutable struct SquaredMatrixPlusI{T<:Real} <: AbstractArray{T,2}
    # A = D*D' + λI
    D::AbstractArray{T};
    λ::T;
    tmp::Vector{T};

    function SquaredMatrixPlusI(D::AbstractArray{T}, λ::T;) where {T<:Real}
        tmp_size = length(size(D)) == 1 ? 1 : size(D)[2];
        new{T}(D, λ, zeros(T, tmp_size));
    end
end

function Base.size(sqm::SquaredMatrixPlusI)
    (size(sqm.D)[1], size(sqm.D)[1])
end

function LinearAlgebra.mul!(c::Vector{T}, A::SquaredMatrixPlusI{T}, b::Vector{T}) where T<:Real
    # c = (A.D*A.D' + A.λ*I)*b
    mul!(A.tmp, A.D', b);
    mul!(c, A.D, A.tmp);
    c .+= A.λ.*b;
end

end


text:
const MKL_DOMAIN_PARDISO = @compat Int32(4)

type MKLPardisoSolver <: AbstractPardisoSolver
    pt::Vector{Int}
    iparm::Vector{Int32}
    mtype::MatrixType
    solver::Solver
    phase::Phase
    msglvl::MessageLevel
    maxfct::Int32
    mnum::Int32
    perm::Vector{Int32}
end

function MKLPardisoSolver()
    if !MKL_PARDISO_LOADED
        error("mkl library was not loaded, unable to create solver")
    end

    pt = zeros(Int, 64)
    iparm = zeros(Int32, 64)
    mtype = REAL_NONSYM
    solver = DIRECT_SOLVER
    phase = ANALYSIS_NUM_FACT_SOLVE_REFINE
    msglvl = MESSAGE_LEVEL_OFF
    mnum = 1
    maxfct = 1
    perm = Int32[]

    MKLPardisoSolver(pt, iparm, mtype, solver,
                      phase, msglvl, maxfct, mnum, perm)
end

show(io::IO, ps::MKLPardisoSolver) = print(io, string("$MKLPardisoSolver:\n",
                                  "\tMatrix type: $(MATRIX_STRING[get_matrixtype(ps)])\n",
                                  "\tPhase: $(PHASE_STRING[get_phase(ps)])"))

set_nprocs!(ps::MKLPardisoSolver, n::Integer) = ccall(set_nthreads, Void, (Ptr{Int32},Ptr{Int32}), &(@compat Int32(n)), &MKL_DOMAIN_PARDISO)
get_nprocs(ps::MKLPardisoSolver) = ccall(get_nthreads, Int32, (Ptr{Int32},), &MKL_DOMAIN_PARDISO)

valid_phases(ps::MKLPardisoSolver) = keys(MKL_PHASES)
phases(ps::MKLPardisoSolver) = MKL_PHASES

function get_matrix(ps::MKLPardisoSolver, A, T)
    mtype = get_matrixtype(ps)

    if isposornegdef(mtype)
        T == :C && return conj(tril(A))
        return tril(A)
    end

    if !issymmetric(mtype)
        T == :C && return conj(A)
        return A
    end

    if mtype == COMPLEX_SYM
        T == :C && return conj(tril(A))
        return tril(A)
    end

    error("Unhandled matrix type")
end

function ccall_pardisoinit(ps::MKLPardisoSolver)
    ERR = Ref{Int32}(0)
    ccall(mkl_init, Void,
          (Ptr{Int}, Ptr{Int32}, Ptr{Int32}),
          ps.pt, &ps.mtype, ps.iparm)
    check_error(ps, ERR[])
end

function ccall_pardiso{Tv}(ps::MKLPardisoSolver, N, AA::Vector{Tv}, IA, JA,
                                   NRHS, B::VecOrMat{Tv}, X::VecOrMat{Tv})
    ERR = Ref{Int32}(0)
    ccall(mkl_pardiso_f, Void,
          (Ptr{Int}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Tv}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32}, Ptr{Tv}, Ptr{Tv},
           Ptr{Int32}),
          ps.pt, &ps.maxfct, &ps.mnum, &ps.mtype, &ps.phase,
          &N, AA, IA, JA, ps.perm,
          &NRHS, ps.iparm, &ps.msglvl, B, X,
          ERR)
    check_error(ps, ERR[])
end

function check_error(ps::MKLPardisoSolver, err::Integer)
    err != -1  || throw(PardisoException("Input inconsistent."))
    err != -2  || throw(PardisoException("Not enough memory."))
    err != -3  || throw(PardisoException("Reordering problem."))
    err != -4  || throw(PardisoPosDefException("Zero pivot, numerical fact. or iterative refinement problem."))
    err != -5  || throw(PardisoException("Unclassified (internal) error."))
    err != -6  || throw(PardisoException("Preordering failed (matrix types 11, 13 only)."))
    err != -7  || throw(PardisoException("Diagonal matrix is singular"))
    err != -8  || throw(PardisoException("32-bit integer overflow problem."))
    err != -9  || throw(PardisoException("Not enough memory for OOC."))
    err != -10 || throw(PardisoException("Error opening OOC files."))
    err != -11 || throw(PardisoException("Read/write error with OOC files."))
    err != -12 || throw(PardisoException("pardiso_64 called from 32-bit library"))
    return
end


text:
# This file was generated, do not modify it. # hide
names(iris, Not(:Species))


text:
using DataFrames, CSV, Statistics, StatsBase

df = CSV.read("./preprocessing/input/land-emissions.csv")

# remove last (empty) column
select!(df, Not(:Model))
select!(df, Not(:Region))

df = df[.!ismissing.(df.Scenario),:]

df = stack(df, 4:size(df)[2])

rename!(df, [:year, :value, :scenario, :variable, :unit])
df[df.unit .== "kt N2O/yr",:value] /= 1000

df.source = map(v -> isnothing(match(r"AFOLU", v)) ? "total" : "afolu", df.variable)
# df.gas = ""
df[df.unit .== "kt N2O/yr",:variable] .= "N₂0"
df[df.unit .== "Mt CO2/yr",:variable] .= "CO₂"
df[df.unit .== "Mt CH4/yr",:variable] .= "CH₄"

df[df.unit .== "kt N2O/yr",:unit] .= "Mt/year"
df[df.unit .== "Mt CO2/yr",:unit] .= "Mt/year"
df[df.unit .== "Mt CH4/yr",:unit] .= "Mt/year"

df = unstack(df, :year, :value)

CSV.write("./src/assets/data/land-emissions.csv", df)

max_stars_repo_name: SensesProject/land-use

max_issues_repo_name: SensesProject/land-use

max_forks_repo_name: SensesProject/land-use


text:
"""
    IColourable(nexid)

An interface for colourable widgets.
"""
struct IColourable <: AbstractINextion
    _nid::NexID

    function IColourable(nid::NexID)
        new(nid)
    end
end

function setproperty!(obj::IColourable, property::Symbol, new_val)

    if property == :backcolor
        setnexproperty!(NexID(obj), :bco, new_val)
    elseif property == :forecolor
        setnexproperty!(NexID(obj), :pco, new_val)
    else
        setfield!(obj, property, new_val)
    end

end

function getproperty(obj::IColourable, property::Symbol)

    if property == :backcolor
        getnexproperty(NexID(obj), :bco, Int)  # ToDo convert to Color
    elseif property == :forecolor
        getnexproperty(NexID(obj), :pco, Int)  # ToDo convert to Color
    else
        getfield(obj, property)
    end

end