text:
# !INPUT/OUTPUT PARAMETERS:
#   sol    : speed of light in atomic units (in,real)
#   ptnucl : .true. if the nucleus is a point particle (in,logical)
#   zn     : nuclear charge (in,real)
#   nst    : number of states to solve for (in,integer)
#   n      : priciple quantum number of each state (in,integer(nst))
#   l      : quantum number l of each state (in,integer(nst))
#   k      : quantum number k (l or l+1) of each state (in,integer(nst))
#   occ    : occupancy of each state (inout,real(nst))
#   xctype : exchange-correlation type (in,integer(3))
#   xcgrad : 1 for GGA functional, 0 otherwise (in,integer)
#   nr     : number of radial mesh points (in,integer)
#   r      : radial mesh (in,real(nr))
#   eval   : eigenvalue without rest-mass energy for each state (out,real(nst))
#   rho    : charge density (out,real(nr))
#   vr     : self-constistent potential (out,real(nr))
#   rwf    : major and minor components of radial wavefunctions for each state
#            (out,real(nr,2,nst))
function solve_atom!(
    sol, ptnucl, zn, nst, n, l, k, occ,
    xctype, xcgrad, nr, r, evals, rho, vr, rwf
)

# real(8), intent(in) :: sol
# logical, intent(in) :: ptnucl
# real(8), intent(in) :: zn
# integer, intent(in) :: nst
# integer, intent(in) :: n(nst),l(nst),k(nst)
# real(8), intent(inout) :: occ(nst)
# integer, intent(in) :: xctype(3),xcgrad
# integer, intent(in) :: nr
# real(8), intent(in) :: r(nr)
# real(8), intent(out) :: eval(nst)
# real(8), intent(out) :: rho(nr),vr(nr)
# real(8), intent(out) :: rwf(nr,2,nst)

    #local variables
    maxscl = 200

    # potential convergence tolerance
    SMALL = 1.0e-6

    @assert nst > 0

    # allocate local arrays
    vn = zeros(Float64,nr)
    vh = zeros(Float64,nr)
    ex = zeros(Float64,nr)
    ec = zeros(Float64,nr)
    vx = zeros(Float64,nr)
    vc = zeros(Float64,nr)
    vrp = zeros(Float64,nr)

    ri = zeros(Float64,nr)
    wpr = zeros(Float64,4,nr)
    fr1 = zeros(Float64,nr)
    fr2 = zeros(Float64,nr)
    gr1 = zeros(Float64,nr)
    gr2 = zeros(Float64,nr)

    #if (xcgrad.eq.1) then
    #  allocate(grho(nr),g2rho(nr),g3rho(nr))
    #end if

    # find total electronic charge
    ze = 0.0
    for ist in 1:nst
        ze = ze + occ[ist]
    end

    # set up nuclear potential
    potnucl!(ptnucl, nr, r, zn, vn)

    for ir in 1:nr
        ri[ir] = 1.0/r[ir]
        # initialise the Kohn-Sham potential to the nuclear potential
        vr[ir] = vn[ir]
    end

    # determine the weights for radial integration
    wsplintp!(nr, r, wpr)

    # initialise mixing parameter
    betamix = 0.5

    # initialise eigenvalues to relativistic values (minus the rest mass energy)
    for ist in 1:nst
        t1 = sqrt( k[ist]^2 - (zn/sol)^2 )
        t1 = ( n[ist] - abs(k[ist]) + t1)^2
        t1 = 1.0 + (zn/sol)^2/t1
        evals[ist] = sol^2/sqrt(t1) - sol^2
    end

    dvp = 0.0
    # start of self-consistent loop
    is_converged = false
    for iscl in 1:maxscl

        # solve the Dirac equation for each state
        for ist in 1:nst
            @views evals[ist] = rdirac!( sol, n[ist], l[ist], k[ist], nr, r, vr, evals[ist],
                rwf[:,1,ist], rwf[:,2,ist] )
        end

        # compute the charge density
        for ir in 1:nr
            ss = 0.0
            for ist in 1:nst
              ss = ss + occ[ist]*( rwf[ir,1,ist]^2 + rwf[ir,2,ist]^2 )
            end
            fr1[ir] = ss
            fr2[ir] = ss*ri[ir]
            rho[ir] = ss*ri[ir]^2 / (4*pi)
        end
        splintwp!(nr, wpr, fr1, gr1)
        splintwp!(nr, wpr, fr2, gr2)

        # find the Hartree potential
        t1 = gr2[nr]
        for ir in 1:nr
           vh[ir] = gr1[ir]*ri[ir] + t1 - gr2[ir]
        end

        # normalise charge density and potential
        t1 = ze/gr1[nr]
        rho[:] .= t1*rho[:]
        vh[:] .= t1*vh[:]

        # compute the exchange-correlation energy and potential
        if xcgrad == 1
            # GGA functional
            # |grad rho|
            fderiv!(1, nr, r, rho, grho)
            # grad^2 rho
            fderiv!(2, nr, r, rho, g2rho)
            for ir in 1:nr
                g2rho[ir] = g2rho[ir] + 2.0*ri[ir]*grho[ir]
            end
            # approximate (grad rho).(grad |grad rho|)
            for ir in 1:nr
                g3rho[ir] = grho[ir]*g2rho[ir]
            end
            #xcifc(xctype,n=nr,rho=rho,grho=grho,g2rho=g2rho,g3rho=g3rho,ex=ex, ec=ec,vx=vx,vc=vc)
        else
            # LDA functional
            #call xcifc(xctype,n=nr,rho=rho,ex=ex,ec=ec,vx=vx,vc=vc)
            for ir in 1:nr
                if rho[ir] <= 10e-12
                    #println("small rho = ", rho[ir])
                    ex[ir] = vx[ir] = 0.0
                    ec[ir] = vc[ir] = 0.0
                else
                    ex[ir], vx[ir] = XC_x_slater( rho[ir] )
                    ec[ir], vc[ir] = XC_c_vwn( rho[ir] )
                end
            end
        end

        #println("sum(vx) = ", sum(vx))
        #println("sum(vc) = ", sum(vc))
        #exit()

        # self-consistent potential
        @views vr[1:nr] = vh[1:nr] + vx[1:nr] + vc[1:nr]

        # determine change in potential
        ss = 0.0
        for ir in 1:nr
          ss = ss + ( vr[ir] - vrp[ir] )^2
        end
        dv = sqrt(ss)/nr

        if iscl > 2
            # reduce beta if change in potential is diverging
            if dv > dvp
                betamix = betamix*0.8
            end
            betamix = max(betamix,0.01)
        end

        dvp = dv

        for ir in 1:nr
            # mix old and new potentials
            vr[ir] = (1.0 - betamix)*vrp[ir] + betamix*vr[ir]
            vrp[ir] = vr[ir]
            # add nuclear potential
            vr[ir] = vr[ir] + vn[ir]
        end

        @printf("iscl = %5d dv = %18.10f\n", iscl, dv)

        # check for convergence
        if ( (iscl > 2) && (dv < SMALL) )
            println(".... Converged ....")
            is_converged = true
            break
        end
    end

    if !is_converged
        println("Warning(atom): maximum iterations exceeded")
    end

    return

end


text:
## Load packages
using Random
using MAT
using StatsBase
using Revise
using Pkg; Pkg.activate(".")
using HypergraphModularity
using SparseArrays
using Statistics
include("src/AON_hyperlouvain.jl")
include("src/synthetic_hypergraphs.jl")
include("src/louvain_utils.jl")
include("src/graph_louvain.jl")
## Run experiment

kmin = 2
kmax = 4
davg = 10

lower = 4
v = range(lower,stop=6,length=10)
Nvals = round.(Int64,[10^i for i in v])

s = 1
Kvals = round.(Int64,sqrt.(Nvals))

# Louvain parameters
maxits = 100
randflag = false
verbose = false
clusterpenalty = 0
r_sizes = ones(kmax-kmin+1) # hyperedges of different sizes, equally likely
N = maximum(Nvals)

aris = zeros(length(Nvals),s)
nmis = zeros(length(Nvals),s)
runs = zeros(length(Nvals),s)
cnum = zeros(length(Nvals),s)

aris_refine = zeros(length(Nvals),s)
nmis_refine = zeros(length(Nvals),s)
runs_refine = zeros(length(Nvals),s)
cnum_refine = zeros(length(Nvals),s)

aris_dyadic = zeros(length(Nvals),s)
nmis_dyadic = zeros(length(Nvals),s)
runs_dyadic = zeros(length(Nvals),s)
cnum_dyadic = zeros(length(Nvals),s)

aris_urefine = zeros(length(Nvals),s)
nmis_urefine = zeros(length(Nvals),s)
runs_urefine = zeros(length(Nvals),s)
cnum_urefine = zeros(length(Nvals),s)

aris_udyadic = zeros(length(Nvals),s)
nmis_udyadic = zeros(length(Nvals),s)
runs_udyadic = zeros(length(Nvals),s)
cnum_udyadic = zeros(length(Nvals),s)

println("Louvain Runtime Performance")
print(rpad("Method", 20))
print(rpad("n", 10))
print(rpad("ARI", 15))
print(rpad("NMI", 15))
print(rpad("#Clusters", 15))
println(rpad("Runtime", 15))
println(rpad("",  80, "-"))
# Kmax = maximum(Kvals)

clusterpenalty = 0
for ni = 1:length(Nvals)
    if ni == 1
        println("")
    end
    global n
    n = round(Int64,Nvals[ni])
    m = davg*n
    edgeweights = ones(m)

    if option == 1
        tag = "sqrtnclusters"
        K = Kvals[ni]
    elseif option == 2
        tag = "200clusters"
        K = 200
    else
        tag = "size200clusters"
        K = round(Int64,n/200)
    end

    cluster_sizes=ones(K)
    cluster_prefs=ones(K)

    pvals = [.6, 1/n^3, 1/n^4]

    # Draw s samples and cluster them with the parameters
    for sample = 1:s
        He2n, e2n, elen, deg, ground_truth = GenerateHypergraphAll(n,m,K,pvals,kmin,kmax,cluster_sizes,r_sizes,cluster_prefs)
        n2e = incidence2elist(SparseArrays.sparse(He2n'))
        H = Elist_to_Hypergraph(e2n)
        #cut_weights, vol_weights, learntime, α = learn_alpha_wrapper(H,ground_truth,kmax,n)
        cut_weights, vol_weights, omega, EdgesAndCuts = learn_omega_aon(e2n,ground_truth,kmax,deg,n)
        γ̂ = computeDyadicResolutionParameter(H, ground_truth)
        gu = computeDyadicResolutionParameter(H, ground_truth;weighted = false)

        # Run with hypergraph version
        Zwarm = collect(1:n)  # no warm start
        tic = time()
        Zs = SuperNode_PPLouvain(n2e,e2n,edgeweights,deg,elen,cut_weights,vol_weights,kmax,randflag,maxits,verbose,Zwarm,clusterpenalty)
        toc = time()-tic
        Z = Zs[:,end]
        clusts = maximum(Z)
        ARI = ari(ground_truth,Z)
        NMI = nmi(ground_truth,Z)

        aris[ni,sample] = ARI
        nmis[ni,sample] = NMI
        runs[ni,sample] = toc
        cnum[ni,sample] = clusts

        # Run dyadic version, weighted
        tic =  time()
        Z_dyadic = CliqueExpansionModularity(H,γ̂;maxits = maxits)
        toc = time()-tic
        clusts = maximum(Z_dyadic)
        ARI = ari(ground_truth,Z_dyadic)
        NMI = nmi(ground_truth,Z_dyadic)
        aris_dyadic[ni,sample] = ARI
        nmis_dyadic[ni,sample] = NMI
        runs_dyadic[ni,sample] = toc
        cnum_dyadic[ni,sample] = clusts

        # Run dyadic version, unweighted
        tic =  time()
        Z_udyadic = CliqueExpansionModularity(H,gu;maxits = maxits,weighted = false)
        toc = time()-tic
        clusts = maximum(Z_udyadic)
        ARI = ari(ground_truth,Z_udyadic)
        NMI = nmi(ground_truth,Z_udyadic)
        aris_udyadic[ni,sample] = ARI
        nmis_udyadic[ni,sample] = NMI
        runs_udyadic[ni,sample] = toc
        cnum_udyadic[ni,sample] = clusts

        # Refine weighted dyadic output
        Zwarm = Z_dyadic
        tic = time()
        Zs = SuperNode_PPLouvain(n2e,e2n,edgeweights,deg,elen,cut_weights,vol_weights,kmax,randflag,maxits,verbose,Zwarm,clusterpenalty)
        toc = time()-tic
        Z = Zs[:,end]
        clusts = maximum(Z)
        ARI = ari(ground_truth,Z)
        NMI = nmi(ground_truth,Z)

        aris_refine[ni,sample] = ARI
        nmis_refine[ni,sample] = NMI
        runs_refine[ni,sample] = toc
        cnum_refine[ni,sample] = clusts

        # Refine unweighted dyadic output
        Zwarm = Z_udyadic
        tic = time()
        Zs = SuperNode_PPLouvain(n2e,e2n,edgeweights,deg,elen,cut_weights,vol_weights,kmax,randflag,maxits,verbose,Zwarm,clusterpenalty)
        toc = time()-tic
        Z = Zs[:,end]
        clusts = maximum(Z)
        ARI = ari(ground_truth,Z)
        NMI = nmi(ground_truth,Z)

        aris_urefine[ni,sample] = ARI
        nmis_urefine[ni,sample] = NMI
        runs_urefine[ni,sample] = toc
        cnum_urefine[ni,sample] = clusts
    end

    runmed = median(runs_udyadic[ni,:])
    arimed = median(aris_udyadic[ni,:])
    nmimed = median(nmis_udyadic[ni,:])
    cnummed = median(cnum_udyadic[ni,:])
    print(rpad("UDyadic", 20))
    print(rpad("$n", 10))
    print(rpad("$(round(arimed,digits = 5))", 15))
    print(rpad("$(round(nmimed,digits = 5))", 15))
    print(rpad("$(round(cnummed))", 15))
    println(rpad("$(round(runmed; digits=3))", 10))

    runmed = median(runs_dyadic[ni,:])
    arimed = median(aris_dyadic[ni,:])
    nmimed = median(nmis_dyadic[ni,:])
    cnummed = median(cnum_dyadic[ni,:])
    print(rpad("WDyadic", 20))
    print(rpad("$n", 10))
    print(rpad("$(round(arimed,digits = 5))", 15))
    print(rpad("$(round(nmimed,digits = 5))", 15))
    print(rpad("$(round(cnummed))", 15))
    println(rpad("$(round(runmed; digits=3))", 10))

    print(rpad("AONLouv", 20))
    runmed = median(runs[ni,:])
    arimed = median(aris[ni,:])
    nmimed = median(nmis[ni,:])
    cnummed = median(cnum[ni,:])
    print(rpad("$n", 10))
    print(rpad("$(round(arimed,digits = 5))", 15))
    print(rpad("$(round(nmimed,digits = 5))", 15))
    print(rpad("$(round(cnummed))", 15))
    print(rpad("$(round(runmed; digits=3))", 10))
    println(rpad("$(K)", 10))
    # println(rpad("",  80, "-"))

    print(rpad("U-Ref", 20))
    runmed = median(runs_urefine[ni,:])
    arimed = median(aris_urefine[ni,:])
    nmimed = median(nmis_urefine[ni,:])
    cnummed = median(cnum_urefine[ni,:])
    print(rpad("$n", 10))
    print(rpad("$(round(arimed,digits = 5))", 15))
    print(rpad("$(round(nmimed,digits = 5))", 15))
    print(rpad("$(round(cnummed))", 15))
    print(rpad("$(round(runmed; digits=3))", 10))
    println(rpad("$(K)", 10))

    print(rpad("W-Ref", 20))
    runmed = median(runs_refine[ni,:])
    arimed = median(aris_refine[ni,:])
    nmimed = median(nmis_refine[ni,:])
    cnummed = median(cnum_refine[ni,:])
    print(rpad("$n", 10))
    print(rpad("$(round(arimed,digits = 5))", 15))
    print(rpad("$(round(nmimed,digits = 5))", 15))
    print(rpad("$(round(cnummed))", 15))
    print(rpad("$(round(runmed; digits=3))", 10))
    println(rpad("$(K)", 10))
    println(rpad("",  90, "-"))

    matwrite("../../fig/performance_throughput/Toughest_regime_experiment_$tag.mat",
    Dict("aris"=>aris,"nmis"=>nmis,"runs"=>runs,"cnum"=>cnum,"pvals"=>pvals,"Kvals"=>Kvals,"Nvals"=>Nvals,
    "aris_dyadic"=>aris_dyadic,"nmis_dyadic"=>nmis_dyadic,"runs_dyadic"=>runs_dyadic,"cnum_dyadic"=>cnum_dyadic,
    "aris_refine"=>aris_refine,"nmis_refine"=>nmis_refine,"runs_refine"=>runs_refine,"cnum_refine"=>cnum_refine,
    "aris_udyadic"=>aris_udyadic,"nmis_udyadic"=>nmis_udyadic,"runs_udyadic"=>runs_udyadic,"cnum_udyadic"=>cnum_udyadic,
    "aris_urefine"=>aris_urefine,"nmis_urefine"=>nmis_urefine,"runs_urefine"=>runs_urefine,"cnum_urefine"=>cnum_urefine
    ))

end

ext: jl

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
## -------------------------------------------------------------------
# proj 2D
let
    me_method = :ME_MAX_POL
    fba_method = :FBA_Z_FIX_MIN_COST

    for (iJR, Data) in [
            # (ChK.iJR904, ChK.KayserData),
            (ChN.iJR904, ChN.NanchenData),
            # (ChF.iJR904, ChF.FolsomData),
            (ChH.iJR904, ChH.HeerdenData),
        ]

        src = nameof(Data)

        DAT = ChE.load_DAT(src)
        EXCH_FLX_IDERS = DAT[:EXCH_FLX_IDERS]
        EXPS = DAT[:EXPS]

        ps_pool = Dict()
        for exp in EXPS

            for ider in EXCH_FLX_IDERS

                # 2D Projection
                p = plot(;title = string(src, " exp:", exp),
                    xlabel = "D", ylabel = ider,
                    legend = :right
                )
                proj = DAT[me_method, :proj, ider, exp]
                ChP.plot_projection2D!(p, proj; l = 50)

                # bounds
                lb, ub = DAT[:bounds, ider, exp]
                hline!(p, [lb]; lw = 3,
                    label = "fva lb",
                    color = :blue, ls = :solid
                )
                hline!(p, [ub]; lw = 3,
                    label = "fva ub",
                    color = :red, ls = :solid
                )

                # EXPERIMENTAL FLXS
                exp_biom = DAT[:exp, :flx, "D", exp]
                exp_exch = DAT[:exp, :flx, ider, exp]
                scatter!(p, [exp_biom], [exp_exch];
                    m = 8, color = :black, label = "exp"
                )

                # MAXENT FLXS
                me_biom = DAT[me_method, :flx, "D", exp]
                me_biom_err = DAT[me_method, :err, "D", exp]
                me_exch = DAT[me_method, :flx, ider, exp]
                me_exch_err = DAT[me_method, :err, ider, exp]
                scatter!(p, [me_biom], [me_exch];
                    xerr = [me_biom_err], yerr = [me_exch_err],
                    m = 8, color = :red, label = "maxent"
                )

                # FBA FLXS
                fba_biom = DAT[fba_method, :flx, "D", exp]
                fba_exch = DAT[fba_method, :flx, ider, exp]
                scatter!(p, [fba_biom], [fba_exch];
                    m = 8, color = :blue, label = "fba"
                )

                ps_pool[(exp, ider)] = deepcopy(p)
            end
        end

        # collect
        for exp in EXPS
            ps = Plots.Plot[ps_pool[(exp, ider)] for ider in EXCH_FLX_IDERS]
            mysavefig(ps, "polytope"; exp, src, me_method, fba_method)
        end

        for ider in EXCH_FLX_IDERS
            ps = Plots.Plot[ps_pool[(exp, ider)] for exp in EXPS]
            mysavefig(ps, "polytope"; ider, src, me_method, fba_method)
        end
    end # for (iJR, Data)
end

## -------------------------------------------------------------------
# TODO: update ChemostatPlots
function plot_projection2D!(p, proj; l = 10, spkwargs...)

    def_spkwargs = (;label = "", alpha = 0.5)

    flx1s = sort(collect(keys(proj)))
    rflx1s = reverse(flx1s)
    flx2lbs = [first(proj[flx1]) for flx1 in flx1s]
    flx2ubs = [last(proj[flx1]) for flx1 in rflx1s]

    s = Shape([flx1s; rflx1s], [flx2lbs; flx2ubs])
    plot!(p, [s]; label = "", def_spkwargs..., spkwargs...)

    return p
end

## -------------------------------------------------------------------
# blue polytope
let
    me_method = :ME_MAX_POL
    # fba_method = :FBA_Z_FIX_MIN_COST

    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
        # (ChH.iJR904, ChH.HeerdenData),
    ]
        src = nameof(Data)

        DAT = ChE.load_DAT(src)
        EXCH_FLX_IDERS = DAT[:EXCH_FLX_IDERS]
        EXPS = DAT[:EXPS]

        for ider in EXCH_FLX_IDERS
            # 2D Projection
            exp = EXPS[1]
            fontsize = 13
            p = plot(;
                # title = string(src, " exp:", exp),
                xlabel = "D", ylabel = ider,
                legend = :right,
                dpi = 1000,
                thickness_scaling = 1.3,
                xguidefontsize = fontsize, yguidefontsize = fontsize
            )
            proj = DAT[me_method, :proj, ider, exp]
            # ChP.plot_projection2D!(p, proj)
            plot_projection2D!(p, proj)

            # EXPERIMENTAL FLXS
            for exp in EXPS
                exp_biom = DAT[:exp, :flx, "D", exp]
                exp_exch = DAT[:exp, :flx, ider, exp]
                marker = (8, source_markers[Data])
                scatter!(p, [exp_biom], [exp_exch];
                    marker, color = :white, label = "",
                )
            end

            # plot!(p; ylim = [-10, 0.0])
            mysavefig(p, "proj2D"; src, ider)
        end
    end
end


text:
include("indexing.jl")

include("windows.jl")


text:
function hessian(f::Function, x::Array{T, 1}) where T <: AbstractFloat

  m = length(x)

	e = eps(T)^(1/4)*maximum(abs,[x;one(T)])
	dh = Matrix{Float64}(I, m, m)*e
  hess = Array{T}(undef,m,m)

  for i = 1:m
		hess[i,i] = (-f(x+2*dh[:,i])+16*f(x+dh[:,i])-30*f(x)+16*f(x-dh[:,i])-f(x-2*dh[:,i]))[1,1]/(12*e^2)
  end

    for i = 1:m-1
      for j = i+1:m
			  hess[i,j] = (f(x+dh[:,i]+dh[:,j])+f(x-dh[:,i]-dh[:,j])-f(x-dh[:,i]+dh[:,j])-f(x+dh[:,i]-dh[:,j]))[1,1]/(4*e^2)
        hess[j,i] = hess[i,j]
      end
    end

    return hess

end


text:
module Models

using LinearAlgebra
using Random
using Logging
using Printf

import Base: eltype, size, length, *
import LinearAlgebra: mul!, ldiv!, transpose!

using ..IterativeSolvers: IterativeSolver, GMRES, ConjugateGradient, BiCGStab, solve!
using ..Checkerboard: checkerboard_mul!, checkerboard_transpose_mul!
using ..Utilities: get_index

export mulM!, mulMᵀ!, mulMᵀM!, muldMdx!, construct_M, write_M_matrix!

"""
Abstract type to represent continuous real or complex numbers.
"""
Continuous = Union{AbstractFloat,Complex{<:AbstractFloat}}

"""
Abstract type to represent type of bond/hopping in lattice.
"""
abstract type AbstractBond end

"""
Represent a type of bare hopping/bond in lattice.
"""
struct Bond{T<:Continuous} <: AbstractBond

    "Average hopping energy."
    t::T

    "Standard Deviation of hopping energy."
    σt::T

    "Starting site/orbital in unit cell."
    o₁::Int

    "Ending site/orbital in unit cell."
    o₂::Int

    "Displacement in unit cells."
    v::Vector{Int}

    function Bond(t::T,σt::T,o₁::Int,o₂::Int,v::AbstractVector{Int}) where {T<:Continuous}

        @assert length(v)==3
        v′ = zeros(T,3)
        copyto!(v′,v)
        return new{T}(t,σt,o₁,o₂,v′)
    end
end

"""
Abstract type to represent models.
    T1: data type for parameter values that are always real
    T2: data type for matrix elements of M
    T3: data type describing which iterative solver is being used
"""
abstract type AbstractModel{T1<:AbstractFloat,T2<:Continuous,T3<:IterativeSolver,T4<:AbstractRNG} end

# include code for models
include("HolsteinModels.jl")
include("SSHModels.jl")

"""
Iteratively solve the linear system M⋅x=b ==> x=M⁻¹⋅b or MᵀM⋅x=b ==> x=[MᵀM]⁻¹⋅b.
"""
function ldiv!(x::AbstractVector, model::AbstractModel{T1,T2,T3}, b::AbstractVector, P;
               maxiter::Int=0)::Tuple{Int,T1,Int} where {T1,T2,T3}

    # set default maxiter
    if maxiter == 0
        maxiter = model.solver.maxiter
    end

    # rever to solving linear system w/o preconditioner is preconditioner is the identity operator
    if P == I

        # solve linear system
        iters, residual_error, flag = ldiv!(x, model, b, maxiter=maxiter)

    else

        # solve linear system
        iters = solve!(x, model, b, model.solver, P, maxiter=maxiter)

        # calculate residual error
        v    = model.v‴
        mul!(v,model,x)
        @. v = v - b
        residual_error = norm(v)/norm(b)

        # detect large residual error
        if residual_error > sqrt(model.solver.tol)

            # maxed out iterations
            if iters == maxiter

                flag = 1
                @info("Hit Max Iters, Residual Error = $residual_error, Iterations = $iters, W/ Preconditioner")
                logger = global_logger()
                flush(logger.stream)
                fill!(x,0)

            # false convergence
            else

                flag = 2
                @info("Large Residual Error = $residual_error, Iterations = $iters, W/ Preconditioner")
                logger = global_logger()
                flush(logger.stream)
                fill!(x,0)
            end

        # converged normally
        else

            flag = 0

        end

        # if failed to converge attempt without preconditioner
        if flag > 0

            # solve linear system
            iters, residual_error, flag = ldiv!(x, model, b, maxiter=10*maxiter)
        end
    end

    return iters, residual_error, flag
end

function ldiv!(x::AbstractVector, model::AbstractModel{T1,T2,T3}, b::AbstractVector;
               maxiter::Int=0)::Tuple{Int,T1,Int} where {T1,T2,T3}

    # set default maxiter
    if maxiter == 0
        maxiter = model.solver.maxiter
    end

    # solve linear system
    iters = solve!(x, model, b, model.solver, maxiter=maxiter)

    # calculate residual error
    v    = model.v‴
    mul!(v,model,x)
    @. v = v - b
    residual_error  = norm(v)/norm(b)

    # detect large residual error
    if residual_error > sqrt(model.solver.tol) # detect large residual error

        # maxed out iterations
        if iters == model.solver.maxiter

            flag = 1
            @info("Hit Max Iters, Residual Error = $residual_error, Iterations = $iters, W/O Preconditioner")
            logger = global_logger()
            flush(logger.stream)
            fill!(x,0)

        # false convergence
        else

            flag = 2
            @info("Large Residual Error = $residual_error, Iterations = $iters, W/O Preconditioner")
            logger = global_logger()
            flush(logger.stream)
            fill!(x,0)
        end

    # converged normally
    else

        flag = 0

    end

    return iters, residual_error, flag
end

"""
Default multiplication routine for AbstractModel type.
"""
function mul!(y::AbstractVector{T2},model::AbstractModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    if model.mul_by_M
        if model.transposed
            mulMᵀ!(y,model,v)
        else
            mulM!(y,model,v)
        end
    else
        if model.transposed
            mulMMᵀ!(y,model,v)
        else
            mulMᵀM!(y,model,v)
        end
    end

    return nothing
end

"""
Perform the multiplication y = MᵀM⋅v
"""
function mulMᵀM!(y::AbstractVector{T2},model::AbstractModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    # y' = M⋅v
    mulM!(model.v′, model, v)

    # y  = Mᵀ⋅y' = MᵀM⋅v
    mulMᵀ!(y, model, model.v′)

    return nothing
end

"""
Perform the multiplication y = MMᵀ⋅v
"""
function mulMMᵀ!(y::AbstractVector{T2},model::AbstractModel{T1,T2},v::AbstractVector{T2}) where {T1,T2}

    # y' = Mᵀ⋅v
    mulMᵀ!(model.v′, model, v)

    # y  = M⋅y' = MMᵀ⋅v
    mulM!(y, model, model.v′)

    return nothing
end

"""
Tranpose M ⇆ Mᵀ (MᵀM ⇆ MMᵀ) with regards to application of mul! routine.
"""
function transpose!(model::AbstractModel)

    model.transposed = !model.transposed
    return nothing
end

"""
Data type of matrix elements of M matrix.
"""
function eltype(model::AbstractModel{T1,T2,T3})::DataType where {T1<:AbstractFloat,T2<:Number,T3<:IterativeSolver}

    return T2
end

"""
Dimension of M matrix.
"""
function length(model::AbstractModel)::Int

    return model.Ndim
end

"""
Dimension of M matrix.
"""
function size(model::AbstractModel)::Tuple{Int,Int}

    return model.Ndim, model.Ndim
end

"""
Dimension of M matrix.
"""
function size(model::AbstractModel,dim::Int)::Int

    return model.Ndim
end

"""
Assign data folder.
"""
function assign_datafolder!(model::AbstractModel,datafolder::String)

    model.datafolder = datafolder
    return nothing
end

"""
Returns the M matrix as a sparse matrix.
"""
function construct_M(model::AbstractModel{T1,T2}, threshold::T1=1e-14) where {T1,T2}

    # to contain M[row,col]=val info for constructing sparse matrix
    rows = Int[]
    cols = Int[]
    vals = T2[]

    # size of model model
    NLτ = size(model,1)

    # represents unit vector
    unitvector = zeros(T2,NLτ)

    # stores columns vector
    colvector = zeros(T2,NLτ)

    # iterating over rows
    for col in 1:NLτ

        # constructing unitvector
        unitvector[mod1(col-1,NLτ)] = 0.0
        unitvector[col]             = 1.0

        # multiply unit vector by M matrix
        mulM!(colvector,model,unitvector)

        # iterate of column vecto
        for row in 1:NLτ

            # if nonzero
            if abs(colvector[row])>threshold

                # save matrix element
                append!(rows,row)
                append!(cols,col)
                append!(vals,colvector[row])
            end
        end
    end

    return rows,cols,vals
end

"""
Write M matrix to file.
"""
function write_M_matrix!(model::AbstractModel{T1,T2}, filename::String, threshold::T1=1e-10) where {T1,T2}

    # construct M matrix
    rows, cols, vals = construct_M(model,threshold)

    # open file
    open(filename,"w") do file

        # write header to file
        write(file,"col row real imag\n")

        # iterate over non-zero matrix elements
        for i in 1:length(rows)

            # write matrix element to file
            @printf file "%d %d %.10f %.10f\n" cols[i] rows[i] real(vals[i]) imag(vals[i])
        end
    end

    return nothing
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics


text:
# Partial correlation

"""
    partialcor(x, y, Z)

Compute the partial correlation of the vectors `x` and `y` given `Z`, which can be
a vector or matrix.
"""
function partialcor(x::AbstractVector, y::AbstractVector, Z::AbstractVecOrMat)
    length(x) == length(y) == size(Z, 1) ||
        throw(DimensionMismatch("Inputs must have the same number of observations"))
    length(x) > 0 || throw(ArgumentError("Inputs must be non-empty"))
    return Statistics.clampcor(_partialcor(x, mean(x), y, mean(y), Z))
end

function _partialcor(x::AbstractVector, μx, y::AbstractVector, μy, Z::AbstractMatrix)
    p = size(Z, 2)
    p == 1 && return _partialcor(x, μx, y, μy, vec(Z))
    z₀   = view(Z, :, 1)
    Zmz₀ = view(Z, :, 2:p)
    μz₀ = mean(z₀)
    rxz = _partialcor(x,  μx,  z₀, μz₀, Zmz₀)
    rzy = _partialcor(z₀, μz₀, y,  μy,  Zmz₀)
    rxy = _partialcor(x,  μx,  y,  μy,  Zmz₀)::typeof(rxz)
    return (rxy - rxz * rzy) / (sqrt(1 - rxz^2) * sqrt(1 - rzy^2))
end

function _partialcor(x::AbstractVector, μx, y::AbstractVector, μy, z::AbstractVector)
    μz = mean(z)

    # Initialize all of the accumulators to 0 of the appropriate types
    Σxx = abs2(zero(eltype(x)) - zero(μx))
    Σyy = abs2(zero(eltype(y)) - zero(μy))
    Σzz = abs2(zero(eltype(z)) - zero(μz))
    Σxy = zero(Σxx * Σyy)
    Σxz = zero(Σxx * Σzz)
    Σzy = zero(Σzz * Σyy)

    # We only want to make one pass over all of the arrays
    @inbounds begin
        @simd for i in eachindex(x, y, z)
            xi = x[i] - μx
            yi = y[i] - μy
            zi = z[i] - μz

            Σxx += abs2(xi)
            Σyy += abs2(yi)
            Σzz += abs2(zi)

            Σxy += xi * yi
            Σxz += xi * zi
            Σzy += zi * yi
        end
    end

    # Individual pairwise correlations
    rxy = Σxy / sqrt(Σxx * Σyy)
    rxz = Σxz / sqrt(Σxx * Σzz)
    rzy = Σzy / sqrt(Σzz * Σyy)

    return (rxy - rxz * rzy) / (sqrt(1 - rxz^2) * sqrt(1 - rzy^2))
end


text:
"""
    generate_data(params::Dict)

Generates DFT/reference data based on a particular potential.
"""
function generate_data(model::String, params::Dict)
    # Load the potential model (E.g. GaN)
    p_model = Symbol(params[model]["potential"])
    p = @eval $p_model($(params[model]))
    # Load parameters
    fit_forces = params["global"]["fit_forces"]
    a = params["global"]["no_train_atomic_conf"]
    b = params["global"]["no_atomic_conf"]
    positions_per_conf = params["global"]["positions_per_conf"]
    rcutfac = params["global"]["rcutfac"]
    # Calculate DFT training data
    training_data = generate_data_aux(p, positions_per_conf, 1, a, rcutfac, fit_forces)
    # Calculate DFT validation data
    validation_data = generate_data_aux(p, positions_per_conf, a+1, b, rcutfac, fit_forces)

    return training_data, validation_data
end

"""
    generate_data_aux(p::Potential, positions_per_conf::Vector,
                      a::Int64, b::Int64, rcutfac::Float64, fit_forces::Bool)

Auxiliar function. See `generate_data`.
"""
function generate_data_aux(p::Potential, positions_per_conf::Vector,
                           a::Int64, b::Int64, rcutfac::Float64, fit_forces::Bool)
    potentials  = [potential_energy(p, positions_per_conf[j], rcutfac) for j = a:b]
    lin_forces = fit_forces ? linearize([forces(p, positions_per_conf[j], rcutfac)
                                         for j = a:b]) : Vector{Float64}()
    return [potentials; lin_forces]
end

"""
    linearize(data::Vector{Force})

Linearizes a vector of forces.
E.g. [Force(1.0,2.0,3.0), Force(4.0,5.0,6.0)] => [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
"""
function linearize(data::Vector{Vector{Force}})
    lin_data = Vector{Float64}()
    for i = 1:length(data), j = 1:length(data[i]), k = 1:3
        push!(lin_data, data[i][j][k])
    end
    return lin_data
end

function get_dft_data(params::Dict)
    path = params["global"]["path"]
    fit_forces = params["global"]["fit_forces"]
    n = params["global"]["no_train_atomic_conf"]
    m = params["global"]["no_atoms_per_conf"]

    dirs = readdir(joinpath(path, "DATA"))
    dft_energies = Vector{Float64}()
    for d in dirs
        open(joinpath(path, "DATA", d, "ENERGY") , "r") do file
            e = parse(Float64, readline(file))
            push!(dft_energies, e)
        end
    end
    training_data = dft_energies[1:n]
    validation_data = dft_energies[n+1:end]

    if fit_forces
        dft_forces = Vector{Float64}()
        for d in dirs
            open(joinpath(path, "DATA", d, "FORCES") , "r") do file
                while !eof(f)
                    f = parse.(Float64, split(readline(file)))
                    push!(dft_forces, f[1])
                    push!(dft_forces, f[2])
                    push!(dft_forces, f[3])
                end
            end
        end
        training_data = [training_data; dft_forces[1:3*n*m]]
        validation_data = [training_data; dft_forces[3*n*m+1:end]]
    end
    return training_data, validation_data
end


text:
__precompile__()
module NeMo

export initnemo
export delconf, stdoutlog, synapseswo
export delnet, addizneuron, addsynapse
export ncount, walltime, simtime, timereset
export newsim, delsim, simstep
export getmembranes

const libnemo = "libnemo"

## Configuration
"""
    newconf()

Create a configuration, return its pointer. Not exported (use `initnemo()` instead)
"""
newconf() = ccall((:nemo_new_configuration, libnemo), Ptr{UInt8}, ())

## FIXME: prevent segfault on multiple calls
"""
    delconf(conf)

Delete a configuration. `conf` is the pointer to the config
Segfault if called on unexisting configuration
"""
function delconf(conf)
    ccall((:nemo_delete_configuration, libnemo),
                  Void, (Ptr{UInt8},), conf)
end

"Enables NeMosim log"
stdoutlog(conf) = ccall((:nemo_log_stdout, libnemo), Ptr{UInt8}, (Ptr{UInt8},), conf) ## Returns null?
#TODO: add cpu/gpu conf

"Set the synapses Write-Only (disable in-simulation read). Faster net creation"
synapseswo(conf) = ccall((:nemo_set_write_only_synapses, libnemo), Ptr{UInt8}, (Ptr{UInt8},), conf)

## Construction
"Create a network (returns pointer)"
newnet() = ccall((:nemo_new_network, libnemo), Ptr{UInt8}, ())
## FIXME: prevent segfault on multiple calls
"Delete a network. Segfault on non-existing"
delnet(net) = ccall((:nemo_delete_network, libnemo), Void, (Ptr{UInt8},), (net));

## TODO: write
## function addntype(net; name = "Izhikevich")
##     ntype = Ref{Cuint}(0)
##     ccall((:nemo_add_neuron_type, libnemo), Ptr{UInt8}, (Ptr{UInt8}, Cstring, Ref{Cuint}), net, name, ntype)
##     ntype
## end

"Add an Izhikevich neuron to the net
args: net, index, iz-params ..."
function addizneuron(net, idx, a, b, c, d, u, v, sigma)
    res = ccall((:nemo_add_neuron_iz, libnemo), Ptr{UInt8}, (Ptr{UInt8}, UInt32, Float32, Float32, Float32, Float32, Float32, Float32, Float32), net, idx, a, b, c, d, u, v, sigma)
    if res != C_NULL
        error("You have added 2 neurons with the same id! Exiting")
    end
    res
end

function addsynapse(net, source, target, delay, weight, isplastic::Bool)
    sid = Ref{UInt8}(0)
    ccall((:nemo_add_synapse, libnemo), Ptr{UInt8}, (Ptr{UInt8}, UInt32, UInt32, UInt32, Float32, UInt8, Ptr{UInt8}), net, source, target, delay, weight, isplastic, sid)
    Int(sid[])
end

function ncount(net)
    count = Ref{Cuint}(0)
    ccall((:nemo_neuron_count, libnemo), Ptr{UInt8}, (Ptr{UInt8}, Ref{Cuint}), net, count);
    Int(count[])
end

## Simulation

newsim(net, conf) = ccall((:nemo_new_simulation, libnemo), Ptr{UInt8}, (Ptr{UInt8}, Ptr{UInt8}), net, conf)

delsim(sim) = ccall((:nemo_delete_simulation, libnemo), Void, (Ptr{UInt8},), sim)

function simstep(sim, neurons, fstim::Vector{UInt32}, istim::Vector{UInt32}, Iistim::Vector{Float32})
    nfired = Ref{Cuint}(0)
    fired::Ref{Vector{Cuint}}
    fired = zeros(neurons)
    istim::Array{UInt32}
    fstim::Array{UInt32}
    Iistim::Array{Float32}
    res = ccall((:nemo_step, libnemo)
                , Ptr{UInt8}
                , (Ptr{UInt8}
                   , Ptr{Cuint}, Csize_t
                   , Ptr{Cuint}, Ptr{Cfloat}, Csize_t
                   , Ptr{Vector{Cuint}}, Ptr{Cuint})
                , sim
                , convert(Ptr{Cuint}, pointer(fstim)), length(fstim)
                , convert(Ptr{Cuint}, pointer(istim)), convert(Ptr{Cuint}, pointer(Iistim)), length(istim)
                , fired, nfired)
    res != C_NULL ? error(errorname()) : fired[][1:nfired[]]
end

const fstim = Array{UInt32}[] ## Probably never used, speedup
const lfstim = length(fstim)
const pointfstim = convert(Ptr{Cuint}, pointer(fstim))
"Overloaded simstep with static fstims. ~3% faster"
function simstep(sim, neurons, istim::Vector{UInt32}, Iistim::Vector{Float32})
    nfired = Ref{Cuint}(0)
    fired::Ref{Vector{Cuint}} = zeros(neurons)
    res = ccall((:nemo_step, libnemo)
                , Ptr{UInt8}
                , (Ptr{UInt8}
                   , Ptr{Cuint}, Csize_t
                   , Ptr{Cuint}, Ptr{Cfloat}, Csize_t
                   , Ptr{Vector{Cuint}}, Ptr{Cuint})
                , sim
                , pointfstim, lfstim
                , convert(Ptr{Cuint}, pointer(istim)), convert(Ptr{Cuint}, pointer(Iistim)), length(istim)
                , fired, nfired)
    res != C_NULL ? error(errorname()) : fired[][1:nfired[]]
end

## Querying the network
"Returns the memrane value of a single neuron (selected by its index)"
function getmembrane(sim::Ptr{UInt8}, idx::UInt32)
    v = Ref{Float32}(0)
    ccall((:nemo_get_membrane_potential, libnemo), Ptr{UInt8},
          (Ptr{UInt8}, Cuint, Ref{Float32}), sim, idx, v)
    v[]
end
"Returns a list of tuples with (idx, potential)
for all membrane valus of the neuron in a given indexes array.
If you don't want automatic indexing, use getmembranes2"
getmembranes(sim::Ptr{UInt8}, idxs::Array{UInt32}) =
    [ (id, getmembrane(sim, id)) for id in idxs]
"Same as getmembranes, but returns just a list
(then ie. enumerate() it)"
getmembranes2(sim::Ptr{UInt8}, idxs::Array{Int}) =
    [ getmembrane(sim, id) for id in idxs]
## How to use vectorize?
## @vectorize_2arg Real getmembrane;

## Macro?
"Wall time elapsed from start/reset"
function walltime(sim)
    time = Ref{UInt}(0)
    ccall((:nemo_elapsed_wallclock, libnemo), Ptr{UInt8},
          (Ptr{UInt8}, Ref{UInt}), sim, time)
    UInt(time[])
end

"Simulation time elapsed"
function simtime(sim)
    time = Ref{UInt}(0)
    ccall((:nemo_elapsed_simulation, libnemo), Ptr{UInt8},
          (Ptr{UInt8}, Ref{UInt}), sim, time)
    UInt(time[])
end
"Reset the time count"
timereset() = ccall((:nemo_reset_timer, libnemo), Ptr{UInt8}, ())

## Error handling
errorname() = bytestring(ccall((:nemo_strerror, libnemo), Ptr{UInt8}, ()))

#######################################################################
"Public. Initialize nemo, returns network, config"
initnemo() = newnet(), newconf()

"Private: add 1 neuron to net"
function neuronadd(net, idx, vals)
    idx == 0 ? error("Theres an error in your config file: Neuron index can't be 0") : "" ## Actually can, but forbidding simplifies julia numbering
    if ismatch(r"^I[zZ]", vals["type"])
        params = vals["params"]
        addizneuron(net, idx, params["a"], params["b"], params["c"], params["d"], params["u"], params["v"], params["s"])
    else error("Only IZ neurons are supported now!") ## FIXME
    end

end

## TODO: move to spikesIO?
"Add a dict of neurons to net. Returns neurons added #"
function neuronsadd(net, nlist)
    for (idx, neuron) in enumerate(nlist)
        neuronadd(net, idx, neuron)
    end
    length(nlist)
end

"Private: add 1 synapse per dest to the net"#source, dest, delay, weight, plastic
function synapseadd(net, s)
    for dest in s["to"]
        addsynapse(net,
                   s["from"],
                   dest,
                   s["values"]["latency"],
                   s["values"]["intensity"],
                   s["values"]["learning"])
    end
end

## TODO: move to spikesIO?
"Add a dict of synapses (improted from JSON) to net"
synapsesadd(net, sdict) = [ synapseadd(net, s) for s in sdict ]

end ## end module

## TODO: move the following to spikesIO
## TODO: refactor synapse definiton
## "1": {
##       "from": 1,
##       "to": [2],
##       "latency": 4,
##       "strength": -20,
##       "learning": false
##     },


text:
# Part 1 of final project: Diffusion equation
if !@isdefined USE_GPU
    const USE_GPU = length(ARGS) > 0 && ARGS[1] == "gpu"
end
using CUDA
import MPI
using ImplicitGlobalGrid
using ParallelStencil
using ParallelStencil.FiniteDifferences3D

## Trick so we can include this file with cpu or gpu
if abspath(PROGRAM_FILE) == @__FILE__
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 3)
        println("Init parallel stencil for GPU.")
    else
        @init_parallel_stencil(Threads, Float64, 3)
        println("Init parallel stencil for CPU.")
    end
end

include("part1_array_programming.jl")
include("part1_kernel_programming.jl")

function main()
    local nx_, ny_, nz_ = 32, 32, 32
    local algorithm = ""
    local bench = "0"
    try
        algorithm = ARGS[2]
    catch
        println("Usage: mpiexecjl -np <np> julia --project scripts/part1.jl [cpu/gpu] [array/kernel] [nx ny nz]")
        exit()
    end

    if length(ARGS) > 2
        try
            nx_, ny_, nz_ = parse.(Int, ARGS[3:5])
        catch
            println("Usage: mpiexecjl -np <np> julia --project scripts/part1.jl [cpu/gpu] [array/kernel] [nx ny nz]")
            exit()
        end
    end
    if length(ARGS) > 5
        bench = ARGS[6]
    end
    bench = (bench == "1")

    if algorithm == "array"
        _, _, bench_res = diffusion_3D_array_programming(nx = nx_, ny = ny_, nz = nz_, verbose=!bench)
    elseif algorithm == "kernel"
        _, _, bench_res = diffusion_3D_kernel_programming(nx = nx_, ny = ny_, nz = nz_, verbose=!bench)
    else
        println("Usage: mpiexecjl -np <np> julia --project scripts/part1.jl [cpu/gpu] [array/kernel] [nx ny nz]")
        exit()
    end
    if bench
        @show bench_res
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
import LinearAlgebra

import RungeKutta.Tableaus: get_lobatto_nodes, get_lobatto_weights, get_gauss_nodes

@doc raw"""
The projective Lobatto-GLRK coefficients are implicitly given by
```math
\sum \limits_{j=1}^{s} a_{ij} c_{j}^{k-1} = \frac{\bar{c}_i^k}{k}  \qquad i = 1 , \, ... , \, \sigma , \; k = 1 , \, ... , \, s ,
```
where $c$ are Gauß-Legendre nodes with $s$ stages and $\bar{c}$ are Gauß-Lobatto nodes with $\sigma$ stages.
"""
function get_lobatto_glrk_coefficients(s, σ=s+1, T=Float64)
    if σ == 1
        @error "Lobatto III coefficients for one stage are not defined."
    end

    c = get_gauss_nodes(s)
    b̄ = get_lobatto_weights(σ)
    c̄ = get_lobatto_nodes(σ)
    M = [ c[j]^(k-1) for k in 1:s, j in 1:s ]

    row(i) = begin
        r = [ c̄[i]^k / k for k in 1:s ]
        M \ r
    end

    ā = vcat([row(i)' for i in 1:σ]...)

    CoefficientsIRK{T}(:LobattoIIIGLRK, s^2, s, σ, ā, b̄, c̄)
end

function get_lobatto_ω_matrix(s)
    as = TableauLobattoIIIA(s).a[2:s,1:s]
    es = zeros(s)
    es[s] = 1

    Q = vcat( hcat(as, zeros(s-1)), hcat(zeros(s)', 1) )
    L = vcat( as, es' )
    ω = inv(L) * Q

    return ω
end


text:
using SHA, Compat

function isdir(filename)
	s = stat(filename)
    s.mode & 0x4000 > 0
end

struct FileLists
	all
	files
	dirs
end

function walk(dirs::Array; kargs...)
	r = walk(first(dirs))
	for i = 2:length(dirs)
		a = walk(dirs[i])
		append!(r.all, a.all)
		append!(r.files, a.files)
		append!(r.dirs, a.dirs)
	end
	r
end

function walk(dir = "", lists = FileLists(Any[],Any[],Any[]))
	items = Base.map(x->joinpath(dir,x), isempty(dir) ? readdir() : readdir(dir))
	dirs = filter(isdir, items)
	files = filter(x->!isdir(x), items)
    append!(lists.all, items)
    append!(lists.files, files)
    append!(lists.dirs, dirs)
	Base.map(x->walk(x, lists), dirs)
	lists
end

watchfile(filename, f) = @async (watch_file(filename); f(filename); watchfile(filename,f))

function watchfiles(f, filenames, watchers)
	for filename in filenames
        if isfile(filename)
            t = mtime(filename)
            watchers[filename] = t
            watchfile(filename, f)
        end
	end
end

firstn(a, n) = a[1:min(n, length(a))]

function parseargs(ARGS = ARGS)
	# -w=dir1,dir2   ... directories to watch, default is all dirs
	# -f=jl,txt      ... filetype to watch, default "jl"
    # --now           ... already execute for the first time on startup, then watch
    # --run           ... everything after this is the command

    splitind = findall("--run" .== ARGS)[1]
    parseuntil = splitind > 0 ? splitind : length(ARGS)
    toparse = ARGS[1:parseuntil]

    isw(a) = startswith(a, "-w=")
    isf(a) = startswith(a, "-f=")
    isnow(a) = startswith(a, "--now")
    isrun(a) = startswith(a, "--run")
    isarg(a) = isw(a) || isf(a) || isnow(a) || isrun(a)
    w = filter(isw, toparse)
    f = filter(isf, toparse)
    now = filter(isnow, toparse)

    for x in ARGS[1:parseuntil]
        if !isarg(x)
            warn("Watcher: Ignoring unknown argument \"$x\"")
        end
    end
	if splitind > 0
        cmd = ARGS[splitind+1:end]
    else
		cmd = ["julia", "test/runtests.jl"]
	end

	dirs = filter(x -> isdir(x) && x[1]!='.', readdir())
	if length(w) > 0
		dirs = split(w[1][4:end],",")
	end
	filetypes = ["jl"]
	if length(f) > 0
		filetypes = map(x-> join([".",x]), split(f[1][4:end],","))
	end
    # @show ARGS
	# @show (dirs, filetypes, cmd, length(now)>0)
	(dirs, filetypes, cmd, length(now)>0)
end

function filehash(filename)
	return hash(read(filename))
end

function runcmd(cmd::Vector{String}, processes::Vector, watchers=Dict(), filename=nothing)
    if filename != nothing
        # t = mtime(filename)
        h = filehash(filename)
        if watchers[filename] == h
            return
        end
        watchers[filename] = h
    end

    @async try
        while !isempty(processes)
            p = pop!(processes)
            kill(p,9)
        end
        println("\n\n-------------------------\n\n")
        process = open(`$cmd`)
        push!(processes, process)
        while in(process, processes) && !eof(process.out)
            x = read(process.out, String)
            in(process, processes) && print(x)
            flush(stdout)
        end
    catch e
        if isa(e, InterruptException)
            kill(process)
        end
    end
end


text:
# Some tests
using Statistics

A = Array{Float32}(1:1e-6:2-1e-6)
Ao = set_one(A,8)
As = shave(A,8)
Ar = round(A,8)
Ag = groom(A,8)

# check bits that are always zero
zero_check = UInt32(0)
for x in As
    global zero_check |= reinterpret(UInt32,x)
end
bitstring(zero_check)

# bias?
mean(As)-mean(A)
mean(Ao)-mean(A)
mean(Ar)-mean(A)
mean(Ag)-mean(A)

# mean error
mean(As-A)
mean(Ao-A)
mean(Ar-A)
mean(Ag-A)

# mean absolute error
mean(abs.(As-A))
mean(abs.(Ao-A))
mean(abs.(Ar-A))
mean(abs.(Ag-A))

# max absolute error
maximum(abs.(As-A))
maximum(abs.(Ao-A))
maximum(abs.(Ar-A))
maximum(abs.(Ag-A))

# max decimal error
maximum(abs.(log10.(As./A)))
maximum(abs.(log10.(Ao./A)))
maximum(abs.(log10.(Ar./A)))
maximum(abs.(log10.(Ag./A)))


text:
@symbol_func function K_WL(cur_reactor::AbstractReactor)
  cur_K = K_W(cur_reactor)

  cur_K *= K_n(cur_reactor) ^ 2

  cur_K /= cur_reactor.max_P_W

  cur_K ^= (1/3)

  cur_K
end


text:
"""
Convenience function for setting an ICD (Installable Client Driver) used by Vulkan.
Only SwiftShader is currently supported. To add another ICD, you must do it by hand. You can achieve that by setting the environment variable "VK_ICD_FILENAMES" to point to your own ICD JSON file, as described in https://github.com/KhronosGroup/Vulkan-Loader/blob/master/loader/LoaderAndLayerInterface.md#icd-discovery.

Available drivers:
- SwiftShader: a CPU implementation of Vulkan. Requires `SwiftShader_jll` to be imported in `mod`.
"""
function set_driver(mod::Module, backend::Symbol)
    @match backend begin
        :SwiftShader => begin
                            libvulkan = mod.SwiftShader_jll.libvulkan
                            ENV["JULIA_VULKAN_LIBNAME"] = basename(libvulkan)
                            libdir = dirname(libvulkan)
                            sep = Sys.iswindows() ? ';' : ':'
                            old_vk_icd_filenames = get(ENV, "VK_ICD_FILENAMES", [])
                            ENV["VK_ICD_FILENAMES"] = join([old_vk_icd_filenames; joinpath(libdir, "vk_swiftshader_icd.json")], sep)
                        end
        _ => error("Backend $backend not available. Only 'SwiftShader' is currently supported.")
    end
end

"""
Call `set_driver` passing the current module as the first argument.
"""
macro set_driver(backend) :(set_driver($__module__, $(esc(backend)))) end


text:
include("../molecule_tools/read_xyz.jl")
include("correlation_functions.jl")
using .Correlation_Functions

header, labels, geoms = read_xyz("data/w6_dynamics_velocity_production_full.xyz")

series = form_atomic_time_series(geoms)

signal = Correlation_Functions.autocorrelation_function(series, 5000)
Correlation_Functions.plot_correlation_function(signal, 2.5)
freqs, vdos = Correlation_Functions.vdos(signal, 2.5)
Correlation_Functions.plot_vdos(freqs, vdos)

max_stars_repo_name: heindelj/pyMD

max_issues_repo_name: heindelj/pyMD

max_forks_repo_name: heindelj/pyMD


text:
using Test

include("darts.jl")

@testset "Missed target" begin
    @test score(-9, 9) == 0
end

@testset "On the outer circle" begin
    @test score(0, 10) == 1
end

@testset "On the middle circle" begin
    @test score(-5, 0) == 5
end

@testset "On the inner circle" begin
    @test score(0, -1) == 10
end

@testset "Exactly on centre" begin
    @test score(0, 0) == 10
end

@testset "Near the centre" begin
    @test score(-0.1, -0.1) == 10
end

@testset "Just within the inner circle" begin
    @test score(0.7, 0.7) == 10
end

@testset "Just outside the inner circle" begin
    @test score(0.8, -0.8) == 5
end

@testset "Just within the middle circle" begin
    @test score(-3.5, 3.5) == 5
end

@testset "Just outside the middle circle" begin
    @test score(-3.6, -3.6) == 1
end

@testset "Just within the outer circle" begin
    @test score(-7.0, 7.0) == 1
end

@testset "Just outside the outer circle" begin
    @test score(7.1, -7.1) == 0
end

@testset "Asymmetric position between the inner and middle circles" begin
    @test score(0.5, -4) == 5
end


text:
using ManifoldLearning

include("nearestneighbors.jl")

k=13
X, L = ManifoldLearning.swiss_roll()

# Use default distance matrix based method to find nearest neighbors
M1 = fit(Isomap, X)
Y1 = transform(M1)

# Use NearestNeighbors package to find nearest neighbors
M2 = fit(Isomap, X, nntype=KDTree)
Y2 = transform(M2)

# Use FLANN package to find nearest neighbors
M3 = fit(Isomap, X, nntype=FLANNTree)
Y3 = transform(M3)

using Plots
plot(
    plot(X[1,:], X[2,:], X[3,:], zcolor=L, m=2, t=:scatter3d, leg=false, title="Swiss Roll"),
    plot(Y1[1,:], Y1[2,:], c=L, m=2, t=:scatter, title="Distance Matrix"),
    plot(Y2[1,:], Y2[2,:], c=L, m=2, t=:scatter, title="NearestNeighbors"),
    plot(Y3[1,:], Y3[2,:], c=L, m=2, t=:scatter, title="FLANN")
, leg=false)


text:
"""
> Yura Malitsky. Projected reflected gradient methods for monotone variational inequalities.
> SIAM Journal on Optimization, 25(1):502–520, 2015.
"""
function ProjectedReflectedGradientVI(stp::AbstractStopping; rho0::Float64 = 0.5)
  xk = stp.current_state.x
  xkp = similar(xk)
  yk = copy(xk)
  rho = rho0
  Fx = similar(xk)

  OK = update_and_start!(stp)
  while !OK
    abcresidual!(stp.pb, xk, rho, yk, Fx)
    project!(stp.pb, Fx, xk) # possible failure here
    yk .= 2 .* xkp .- xk

    if norm(xk - xkp, Inf) < stp.meta.atol * rho
      stp.meta.optimal = true
    end
    OK = update_and_stop!(stp, x = xk)
  end

  return xk
end

export ProjectedReflectedGradientVI

function ProjectedReflectedGradientVI(
  model::AbstractVIModel,
  x0::AbstractVector;
  rho0::Float64 = 0.5,
  kwargs...,
)
  stp = GenericStopping(model, x0; kwargs...)
  return ProjectedReflectedGradientVI(stp, rho0 = rho0)
end

function abcresidual!(model, xk, rho, yk, Fx) # xk - rho * F(yk)
  residual!(model, yk, Fx)
  Fx .*= -rho
  Fx .+= xk
  return Fx
end


text:
export fwi_objective, lsrtm_objective, fwi_objective!, lsrtm_objective!

function multi_src_fg(model_full::Model, source::judiVector, dObs::judiVector, dm, options::JUDIOptions, nlind::Bool, lin::Bool)
# Setup time-domain linear or nonlinear foward and adjoint modeling and interface to OPESCI/devito

    # assert this is for single source LSRTM
    @assert source.nsrc == 1 "Multiple sources are used in a single-source fwi_objective"
    @assert dObs.nsrc == 1 "Multiple-source data is used in a single-source fwi_objective"

    # Load full geometry for out-of-core geometry containers
    dObs.geometry = Geometry(dObs.geometry)
    source.geometry = Geometry(source.geometry)

    # Limit model to area with sources/receivers
    if options.limit_m == true
        model = deepcopy(model_full)
        model, dm = limit_model_to_receiver_area(source.geometry, dObs.geometry, model, options.buffer_size; pert=dm)
    else
        model = model_full
    end

    # Set up Python model
    modelPy = devito_model(model, options, dm)
    dtComp = convert(Float32, modelPy."critical_dt")

    # Extrapolate input data to computational grid
    qIn = time_resample(make_input(source), source.geometry, dtComp)[1]
    dObserved = time_resample(make_input(dObs), dObs.geometry, dtComp)[1]

    # Set up coordinates
    src_coords = setup_grid(source.geometry, model.n)  # shifts source coordinates by origin
    rec_coords = setup_grid(dObs.geometry, model.n)    # shifts rec coordinates by origin

    length(options.frequencies) == 0 ? freqs = nothing : freqs = options.frequencies
    argout1, argout2 = pycall(ac."J_adjoint", Tuple{Float32, PyArray}, modelPy,
                  src_coords, qIn, rec_coords, dObserved, t_sub=options.subsampling_factor,
                  space_order=options.space_order, checkpointing=options.optimal_checkpointing,
                  freq_list=freqs, isic=options.isic, is_residual=false, born_fwd=lin, nlind=nlind,
                  dft_sub=options.dft_subsampling_factor[1], f0=options.f0, return_obj=true)

    argout2 = remove_padding(argout2, modelPy.padsizes; true_adjoint=options.sum_padding)
    if options.limit_m==true
        argout2 = extend_gradient(model_full, model, argout2)
    end
    return Ref{Float32}(argout1), PhysicalParameter(argout2, model_full.d, model_full.o)
end

# Find number of experiments
"""
    get_nexp(x)

Get number of experiments given a JUDI type. By default we have only one experiment unless we input
a Vector of judiType such as [model, model] to compute gradient for different cases at once.
"""
get_nexp(x) = 1
for T in [judiVector, Model, judiWeights, judiWavefield, PhysicalParameter, Vector{Float32}]
    @eval get_nexp(v::Vector{<:$T}) = length(v)
    @eval get_nexp(v::Tuple{N, <:$T}) where N = length(v)
end

# Filter arguments for given task
"""
    get_exp(x, i)

Filter input `x`` for experiment number `i`. Returns `x` is a constant not depending on experiment.
"""
get_exp(x, i) = x
get_exp(x::Tuple{}, i::Any) = x[i]
for T in [judiVector, Model, judiWeights, judiWavefield, Array{Float32}, PhysicalParameter]
    @eval get_exp(v::Vector{<:$T}, i) = v[i]
    @eval get_exp(v::NTuple{N, <:$T}, i) where N = v[i]
end

function check_args(args...)
    n = [get_nexp(a) for a in args]
    nexp = maximum(n)
    check = all(ni -> (ni==nexp || ni==1), n)
    check || throw(ArgumentError("Incompatible number of experiements"))
    return nexp
end

# Type of accepted input
Dtypes = Union{<:judiVector, NTuple{N, <:judiVector} where N, Vector{<:judiVector}}
MTypes = Union{Model, NTuple{N, Model} where N, Vector{Model}}
dmTypes = Union{dmType, NTuple{N, dmType} where N, Vector{dmType}}

"""
    fwi_objective(model, source, dobs; options=Options())

    Evaluate the full-waveform-inversion (reduced state) objective function. Returns a tuple with function value and vectorized \\
gradient. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======
    function_value, gradient = fwi_objective(model, source, dobs)
"""
function fwi_objective(model::MTypes, q::Dtypes, dobs::Dtypes; options=Options(), kw...)
    n_exp = check_args(model, q, dobs)
    G = n_exp == 1 ? similar(model.m) : [similar(get_exp(model, i).m) for i=1:n_exp]
    f = fwi_objective!(G, model, q, dobs; options=options, kw...)
    f, G
end

"""
    lsrtm_objective(model, source, dobs, dm; options=Options(), nlind=false)

Evaluate the least-square migration objective function. Returns a tuple with function value and \\
gradient. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======
    function_value, gradient = lsrtm_objective(model, source, dobs, dm)
"""
function lsrtm_objective(model::MTypes, q::Dtypes, dobs::Dtypes, dm::dmTypes; options=Options(), nlind=false, kw...)
    n_exp = check_args(model, q, dobs, dm)
    G = n_exp == 1 ? similar(model.m) : [similar(get_exp(model, i).m) for i=1:n_exp]
    f = lsrtm_objective!(G, model, q, dobs, dm; options=options, nlind=nlind, kw...)
    f, G
end

"""
    fwi_objective!(G, model, source, dobs; options=Options())

    Evaluate the full-waveform-inversion (reduced state) objective function. Returns a the function value and assigns in-place \\
the gradient to G. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======
    function_value = fwi_objective!(gradient, model, source, dobs)
"""
function fwi_objective!(G, model::MTypes, q::Dtypes, dobs::Dtypes; options=Options(), kw...)
    n_exp = check_args(G, model, dobs, q)
    return multi_exp_fg!(Val(n_exp), G, model, q, dobs, nothing; options=options, nlind=false, lin=false, kw...)
end

"""
    lsrtm_objective!(G, model, source, dobs, dm; options=Options(), nlind=false)

    Evaluate the least-square migration (data-space) objective function. Returns the function value and assigns in-place \\
the gradient to G. `model` is a `Model` structure with the current velocity model and `source` and `dobs` are the wavelets and \\
observed data of type `judiVector`.

Example
=======
    function_value = lsrtm_objective!(gradient, model, source, dobs, dm; options=Options(), nlind=false)
"""
function lsrtm_objective!(G, model::MTypes, q::Dtypes, dobs::Dtypes, dm::dmTypes; options=Options(), nlind=false, kw...)
    n_exp = check_args(G, model, q, dobs, dm)
    return multi_exp_fg!(Val(n_exp), G, model, q, dobs, dm; options=options, nlind=nlind, lin=true, kw...)
end

multi_exp_fg!(n::Val{1}, ar...; kw...) = multi_src_fg!(ar...; kw...)

function multi_exp_fg!(n::Val{N}, ar...; kw...) where N
    f = zeros(Float32, N)
    @sync for i=1:N
        ai = (get_exp(a, i) for a in ar)
        @async f[i] = multi_src_fg!(ai...; kw...)
    end
    sum(f)
end


text:
function ACQR_bench(bench,nt,seed=1)
	#
	# benchmark models to test the functions
	#

	# random seed for the simulations
	srand(seed)

	if bench==1
		#
		A = 0.6
		C = 1
		Q = 2
		R = 4

		v = sqrt(R)*randn(nt)
		w = sqrt(Q)*randn(nt)

		x = zeros(nt)
		for t in 1:nt-1
			x[t+1] = A*x[t] + w[t]
		end
		y = C*x + v

		salida = y,A,C,Q,R

	elseif bench==2
		#
		# Based on:
		# "Subspace Identification for Linear Systems"
		# P. Van Overschee and B. De Moor, page 81
		#
		A = [0.6 0.6 0.;-0.6 0.6 0.;0. 0. 0.4]
		C = [0.7831 0.5351 0.9701]
		R = 6.0
		Q = [3 0.5 0.2;0.5 2 0.3;0.2 0.3 1]
		U = chol(Q)

		v = sqrt(R)*randn(1,nt)
		ns = size(A,1)
		w = U'*randn(ns,nt)

 		x = zeros(ns,nt)
		x[:,1] = w[:,1]
		for t in 1:nt-1
			x[:,t+1] = A*x[:,t] + w[:,t]
		end
		y = C*x + v

		salida = y,A,C,Q,R

	end

	return salida

end


text:
using JuMP
using AxisArrays
using MathOptInterface
using AxisArrays
using LinearAlgebra
using MutableArithmetics
using Printf
using Pkg
using JSON
using Logging

@info pwd()
mi = JSON.parsefile(joinpath(ENV["PROJECT_PATH"], "julia_envs", "Xpress", "modelinputs.json"))

Pkg.build("PyCall")
Pkg.build("Xpress")

# somehow changing working dir during the image build?
include(joinpath(ENV["PROJECT_PATH"], "reo", "src", "reopt_xpress_model.jl"))
include(joinpath(ENV["PROJECT_PATH"], "reo", "src", "reopt.jl"))

model = reopt_model(300)
results = reopt(model, mi)


text:
function singular_values_experiments(testing=true)

    if testing
        #  -- test_params --  #
        n_sizes = [100]
        trials = 2
        result_path = all_results_path*"/rank1_singular_values/test/"
    else
        n_sizes = [100, 500, 1000, 2000, 5000, 10000]
        trials = 50
        result_path = all_results_path*"rank1_singular_values/"
    end

    if !isdir(result_path)
        mkpath(result_path)
    end

    alphas = [1.0]
    betas = [0.0]

    # -- ER params -- #
    p_remove = [.05]
    # -- Duplication params -- #
    edge_inclusion_p = [.5]
    step_percentage = [.25]

    X = distributed_random_trials(trials,ErdosRenyiNoise(),true;method=TAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), p_remove,n_sizes,profile=true,no_matching=true,
                                    postProcessing=noPostProcessing(),tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"TAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(p_remove)_noiseModel:ER_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)

    X = distributed_random_trials(trials,DuplicationNoise(),true;method=TAME_M(),graph=RandomGeometric(),iter=15,
                                    degreedist=LogNormal(log(5),1), edge_inclusion_p, step_percentage,n_sizes,
                                    profile=true,no_matching=true, postProcessing=noPostProcessing(),
                                    tol=1e-12,alphas,betas)
    filename=filter(x-> !isspace(x),"TAME_noMatch_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_p:$(edge_inclusion_p)_noiseModel:Duplication_sp:$(step_percentage)_trials:$(trials)_MaxRankResults.json")
    save_results(result_path*filename,X)

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
function DFS!(G, p, k)
  #inicijalizacija cvorova za obilazak
  for u = 1:length(G.V)
      G.V[u].color = 'W'
  end
  #poziv za pretragu puta od pocetnog do krajnjeg cvora
  return DFS_Visit!(G, p, k)
end

function DFS_Visit!(G, p, k)
   #trazimo put od cvora p do cvora k
   path = Int[]
   #ako smo stigli do krajnjeg cvora,
   #vracamo putanju sa tim cvorom kao rezultat
   if p == k
      path = [p]
      return path
   end
   #trazimo putanju preko susednih cvorova
   G.V[p].color = 'G';
   for v in findall(G.AdjMatrix[p,:] .==1)
      if G.V[v].color == 'W' && G.V[v].type  == 'P'
         path = DFS_Visit!(G, v, k);
         # ako je nadjen put od suseda do krajnjeg cvora
         # taj put prosirujemo sa trenutnim cvorom i prekidamo pretragu
         if length(path) != 0
             return [p; path]
         end
      end
   end
   G.V[p].color = 'B'
   return path
end


text:
################################################################################

"""

    extractFFT(edfDf::DataFrames.DataFrame, params::Dict)

# Description
Use `extractFFT` on EDF file per channel from shell arguments. Returns a dictionary with channel names as keys.

See also: [`extractSignalBin`](@ref)
"""
function extractFFT(edfDf::DataFrames.DataFrame, params::Dict)
  if haskey(params, "window-size") && haskey(params, "bin-overlap")
    return extractFFT(edfDf, binSize = params["window-size"], binOverlap = params["bin-overlap"])
  else
    @error "Variables are not defined in dictionary"
  end
end

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

"""

    extractFFT(channel::Array;
    binSize::T, binOverlap::T) where T <: Number

# Description
Apply fast fourier transform (FFT) to channel.

See also: [`extractSignalBin`](@ref)
"""
function extractFFT(channel::Array; binSize::T, binOverlap::T) where T <: Number
  # define variables
  stepSize = floor(Int32, binSize / binOverlap)
  signalSteps = 1:stepSize:length(channel)
  freqAr = Array{Float64}(
    undef,
    length(signalSteps),
    binSize
  )

  # iterate over signal bins
  for ι ∈ eachindex(signalSteps)
    signalBoundry = signalSteps[ι] + binSize - 1

    # extract channel
    if signalBoundry <= length(channel)
      channelExtract = [
        channel[signalSteps[ι]:signalBoundry];
        zeros(binSize)
      ]

    # adjust last bin
    elseif signalBoundry > length(channel)
      channelExtract = [
        channel[signalSteps[ι]:end];
        (signalBoundry - length(channel) |> abs |> zeros);
        zeros(binSize)
      ]
    end

    # calculate fourier transform
    fftChannel = fft(channelExtract)
    realFft = abs.(fftChannel)
    freqAr[ι, :] = realFft[1:binSize]

  end
  return freqAr
end

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

"""

    extractFFT(edfDf::DataFrames.DataFrame;
    binSize::T, binOverlap::T) where T <: Number

# Description
Use `extractFFT` on EDF file per channel. Returns a dictionary with channel names as keys.

See also: [`extractSignalBin`](@ref)
"""
function extractFFT(edfDf::DataFrames.DataFrame; binSize::T, binOverlap::T) where T <: Number
  @info "Extracting channels frecuencies..."
  channelDc = Dict()
  freqAr = begin
    stepSize = floor(Int32, binSize / binOverlap)
    signalSteps = 1:stepSize:size(edfDf, 1)
    Array{Float64}(
      undef,
      1,
      binSize,
      length(signalSteps)
    )
  end

  # iterate on dataframe channels
  for (ψ, ε) ∈ enumerate(names(edfDf))
    α = extractFFT(edfDf[:, ψ], binSize = binSize, binOverlap = binOverlap)
    for β ∈ 1:size(α, 1)
      freqAr[:, :, β] = α[β, :]
    end
    channelDc[ε] = copy(freqAr)
  end
  return channelDc
end

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

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
# simplified form of https://github.com/JuliaData/DataFrames.jl/blob/master/src/other/index.jl
abstract type AbstractIndex end

struct Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{AbstractString, Int}      # name => names array position
    names::Vector{AbstractString}
end

function Index(names::Array{T, 1}) where T <: AbstractString
    @assert allunique(names) "names must be unique check for $names"
    lookup = Dict{AbstractString, Int}(zip(names, 1:length(names)))
    Index(lookup, names)
end

Index() = Index(Dict{AbstractString, Int}(), String[])
Base.length(x::Index) = length(x.names)
Base.names(x::Index) = copy(x.names)
_names(x::Index) = x.names
Base.copy(x::Index) = Index(copy(x.lookup), copy(x.names))
Base.isequal(x::AbstractIndex, y::AbstractIndex) = _names(x) == _names(y) # it is enough to check names
Base.:(==)(x::AbstractIndex, y::AbstractIndex) = isequal(x, y)

Base.haskey(x::Index, key::AbstractString) = haskey(x.lookup, key)
Base.haskey(x::Index, key::Integer) = 1 <= key <= length(x.names)
Base.haskey(x::Index, key::Bool) =
    throw(ArgumentError("invalid key: $key of type Bool"))
Base.keys(x::Index) = names(x)

@inline Base.getindex(x::AbstractIndex, idx::Bool) = throw(ArgumentError("invalid index: $idx of type Bool"))

@inline function Base.getindex(x::AbstractIndex, idx::Integer)
    if !(1 <= idx <= length(x))
        throw(BoundsError("attempt to access a Index with $(length(x)) columns at index $idx"))
    end
    Int(idx)
end

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{Int})
    isempty(idx) && return idx
    minidx, maxidx = extrema(idx)
    if minidx < 1
        throw(BoundsError("attempt to access a Index with $(length(x)) columns at index $minidx"))
    end
    if maxidx > length(x)
        throw(BoundsError("attempt to access a Index with $(length(x)) columns at index $maxidx"))
    end
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    idx
end

@inline function Base.getindex(x::AbstractIndex, idx::AbstractRange{Int})
    isempty(idx) && return idx
    minidx, maxidx = extrema(idx)
    if minidx < 1
        throw(BoundsError("attempt to access a Index with $(length(x)) columns at index $minidx"))
    end
    if maxidx > length(x)
        throw(BoundsError("attempt to access a Index with $(length(x)) columns at index $maxidx"))
    end
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    idx
end

@inline Base.getindex(x::AbstractIndex, idx::AbstractRange{<:Integer}) = getindex(x, collect(Int, idx))
@inline Base.getindex(x::AbstractIndex, ::Colon) = Base.OneTo(length(x))

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{<:Integer})
    if any(v -> v isa Bool, idx)
        throw(ArgumentError("Bool values except for AbstractVector{Bool} are not allowed for column indexing"))
    end
    getindex(x, Vector{Int}(idx))
end

@inline Base.getindex(x::AbstractIndex, idx::AbstractRange{Bool}) = getindex(x, collect(idx))

@inline function Base.getindex(x::AbstractIndex, idx::AbstractVector{Bool})
    length(x) == length(idx) || throw(BoundsError(x, idx))
    findall(idx)
end

# catch all method handling cases when type of idx is not narrowest possible, Any in particular
@inline function Base.getindex(x::AbstractIndex, idxs::AbstractVector)
    isempty(idxs) && return Int[] # special case of empty idxs
    if idxs[1] isa Real
        if !all(v -> v isa Integer && !(v isa Bool), idxs)
            throw(ArgumentError("Only Integer values allowed when indexing by vector of numbers"))
        end
        return getindex(x, convert(Vector{Int}, idxs))
    end
    idxs[1] isa AbstractString && return getindex(x, convert(Vector{AbstractString}, idxs))
    throw(ArgumentError("idxs[1] has type $(typeof(idxs[1])); "*
                        "Only Integer or String values allowed when indexing by vector"))
end

@inline function Base.getindex(x::AbstractIndex, rx::Regex)
    getindex(x, filter(name -> occursin(rx, String(name)), _names(x)))
end

"""
    fuzzymatch(l::Dict, idx::AbstractString)
# Fuzzy matching rules:
# 1. ignore case
# 2. maximum Levenshtein distance is 2
# 3. always show matches with 0 difference (wrong case)
# 4. on top of 3. do not show more than 8 matches in total
# Returns candidates ordered by (distance, name) pair
"""
function fuzzymatch(l::Dict{AbstractString, Int}, idx::AbstractString)
        idxs = uppercase(idx)
        dist = [(REPL.levenshtein(uppercase(x), idxs), x) for x in keys(l)]
        sort!(dist)
        c = [count(x -> x[1] <= i, dist) for i in 0:2]
        maxd = max(0, searchsortedlast(c, 8) - 1)
        [s for (d, s) in dist if d <= maxd]
end

@inline function lookupname(l::Dict{AbstractString, Int}, idx::AbstractString)
    i = get(l, idx, nothing)
    if i === nothing
        candidates = fuzzymatch(l, idx)
        if isempty(candidates)
            throw(ArgumentError("column name :$idx not found in the data frame"))
        end
        candidatesstr = join(string.(':', candidates), ", ", " and ")
        throw(ArgumentError("column name :$idx not found in the data frame; " *
                            "existing most similar names are: $candidatesstr"))
    end
    i
end

@inline Base.getindex(x::Index, idx::AbstractString) = lookupname(x.lookup, idx)
@inline function Base.getindex(x::Index, idx::AbstractVector{AbstractString})
    allunique(idx) || throw(ArgumentError("Elements of $idx must be unique"))
    [lookupname(x.lookup, i) for i in idx]
end


text:
# Adjoint test for convolution functions provided in NNlib
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using NNlib, Test, LinearAlgebra

###################################################################################################
# Adjoint test 2D

d = 4
n_in = 10
n_out = 20

X = glorot_uniform(32, 32, n_in, 1)
Y = glorot_uniform(8, 8, n_out, 1)
W = glorot_uniform(d, d, n_in, n_out)

# Convolution using NNlib functions
cdims = DenseConvDims(X, W; stride=(d,d), padding=(1,1))

Y_ = conv(X, W, cdims)
X_ = ∇conv_data(Y, W, cdims)    # corresponds to transpose operation

a = vec(Y)'*vec(Y_)
b = vec(X)'*vec(X_)

@test isapprox(a/b - 1f0, 0f0, atol=1f-4)

###################################################################################################
# Adjoint test 3D

d = 4
n_in = 10
n_out = 20

X = glorot_uniform(32, 32, 32, n_in, 1)
Y = glorot_uniform(8, 8, 8, n_out, 1)
W = glorot_uniform(d, d, d, n_in, n_out)

# Convolution using NNlib functions
cdims = DenseConvDims(X, W; stride=(d,d,d), padding=(1,1,1))

Y_ = conv(X, W, cdims)
X_ = ∇conv_data(Y, W, cdims)    # corresponds to transpose operation

a = vec(Y)'*vec(Y_)
b = vec(X)'*vec(X_)

@test isapprox(a/b - 1f0, 0f0, atol=1f-4)

###################################################################################################
# Gradient test 2D

# Input
X = glorot_uniform(64, 64, n_in, 1)
W = glorot_uniform(d, d, n_in, n_out)
W0 = glorot_uniform(d, d, n_in, n_out)
dW = W - W0
b = glorot_uniform(n_out)
b0 = glorot_uniform(n_out)
db = b - b0

# Operators
d = 4
k = 20

cdims = DenseConvDims(X, W; stride=(d,d), padding=(1,1))

function objective(W, b, X, Y)
    ΔY = (conv(X, W, cdims) .+ reshape(b, 1, 1, :, 1)) - Y
    f = .5f0*norm(ΔY)^2f0
    gW = ∇conv_filter(X, ΔY, cdims)
    gb = sum(ΔY, dims=[1,2,4])[1, 1, :, 1]
    return f, gW, gb
end

Y = conv(X, W, cdims) .+ reshape(b, 1, 1, :, 1)

# Gradient test weights
f0, gW, gb = objective(W0, b, X, Y)
h = 0.01f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = objective(W0 + h*dW, b, X, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dW, gW))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test bias
f0, gW, gb = objective(W, b0, X, Y)
h = 0.01f0
maxiter = 6
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = objective(W, b0 + h*db, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(db, gb))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###################################################################################################
# Gradient test 3D

# Input
X = glorot_uniform(64, 64, 64, n_in, 1)
W = glorot_uniform(d, d, d, n_in, n_out)
W0 = glorot_uniform(d, d, d, n_in, n_out)
dW = W - W0
b = glorot_uniform(n_out)
b0 = glorot_uniform(n_out)
db = b - b0

# Operators
d = 4
k = 20

cdims = DenseConvDims(X, W; stride=(d,d,d), padding=(1,1,1))

function objective(W, b, X, Y)
    ΔY = (conv(X, W, cdims) .+ reshape(b, 1, 1, 1, :, 1)) - Y
    f = .5f0*norm(ΔY)^2f0
    gW = ∇conv_filter(X, ΔY, cdims)
    gb = sum(ΔY, dims=[1,2,3,5])[1, 1, 1, :, 1]
    return f, gW, gb
end

Y = conv(X, W, cdims) .+ reshape(b, 1, 1, 1, :, 1)

# Gradient test weights
f0, gW, gb = objective(W0, b, X, Y)
h = 0.01f0
maxiter = 6
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = objective(W0 + h*dW, b, X, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dW, gW))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test bias
f0, gW, gb = objective(W, b0, X, Y)
h = 0.01f0
maxiter = 6
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)

print("\nGradient test convolutions\n")
for j=1:maxiter
    f = objective(W, b0 + h*db, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(db, gb))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
@testset "Conversion from and to SBML model" begin
    json_model = model_paths["iJO1366.json"]

    jm = load_json_model(json_model)
    sm = convert(StandardModel, jm)
    jm2 = convert(JSONModel, sm)

    @test Set(reactions(jm)) == Set(reactions(sm))
    @test Set(reactions(jm)) == Set(reactions(jm2))
end

@testset "JSONModel generic interface" begin
    model = load_model(model_paths["e_coli_core.json"])

    @test reaction_stoichiometry(model, "EX_ac_e") == Dict("ac_e" => -1)
end


text:
export
    astrology,
    biorhythm

"""
    biorhythm(jd_origin, jd_current)

Some pseudoscience. Don't know where these come from, but they're quite funny.

Return some biorhythm stuff

Parameters:
    jd_origin
    jd_current

Return:
    four biorhythms
"""
function biorhythm(jd_origin, jd_current)
    pi2 = pi * 2

    t = jd_current - jd_origin
    physical     = round(sin(pi2 * t/23), digits=2)
    emotional    = round(sin(pi2 * t/28), digits=2)
    intellectual = round(sin(pi2 * t/33), digits=2)
    intuitive    = round(sin(pi2 * t/38), digits=2)

    return (physical, emotional, intellectual, intuitive)
end

"""
    astrology(body, jd)

Return star sign for solar system body for a particular date

Parameters:
    body : "Sun" or planet, one of "Mercury", "Venus", "Earth", "Mars",
        "Jupiter", "Saturn", "Uranus", "Neptune"

Return:
    starsign: string of star sign
"""
function astrology(body, jd)
    global days_per_second
    signs = ["Aries", "Taurus", "Gemini", "Cancer", "Leo", "Virgo", "Libra", "Scorpio",
    "Sagittarius", "Capricorn", "Aquarius", "Pisces"]

    if body == "Moon"
        return "not available"
    end
    if body == "Sun"
    	(L, B, R) = sun_dimension3(jd)
    	long = L
    else
        deltaPsi = nut_in_lon(jd)
        epsilon  = true_obliquity(jd)
        ra, decl = geocentric_planet(jd, body, deltaPsi, epsilon, days_per_second)
        long, lat = equ_to_ecl(ra, decl, epsilon)
    end
    i = 1 + floor(Integer, rad2deg(long) / 30)
    return signs[i]
end


text:
module TestDistributedReduce

include("preamble.jl")
using BangBang
using Distributed

if get(ENV, "CI", "false") == "true"
    addprocs(3)
end
@info "Testing with:" nworkers()

@testset begin
    fname = gensym(:attach_pid)
    @everywhere $fname(x) = [(x, getpid())]
    fun = getproperty(Main, fname)

    pids = fetch.(map(id -> remotecall(getpid, id), workers()))
    xs = 1:10
    ys = dreduce(append!!, Map(fun), xs; init=Union{}[])
    @test first.(ys) == xs
    @test Set(last.(ys)) == Set(pids)

    ys = dreduce(append!!, Map(fun), withprogress(xs; interval=0); init=Union{}[])
    @test first.(ys) == xs
    @test Set(last.(ys)) == Set(pids)
end

end  # module


text:
module Dashboards
import HTTP, JSON2
using MacroTools
include("ComponentPackages.jl")
include("ComponentMetas.jl")
include("Components.jl")

import .ComponentPackages
using .ComponentMetas
using .Components

export Dash, Component, @use, <|, @props_str, @callid_str, CallbackId, callback!, make_handler

ComponentPackages.@reg_components()

@doc """
    module Dashboards

Julia backend for [Plotly Dash](https://github.com/plotly/dash)

# Examples
```julia
import HTTP
using Dashboards
app = Dash("Test", external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"]) do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        dcc_graph(id="graph",
            figure = (
                data = [(x = [1,2,3], y = [3,2,8], type="bar")],
                layout = Dict(:title => "Graph")
            )
        )

    end
end
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children") do type, value
    "You've entered: '\$(value)' into a '\$(type)' input control"
end
callback!(app, callid"graphTitle.value => graph.figure") do value
    (
        data = [
            (x = [1,2,3], y = abs.(randn(3)), type="bar"),
            (x = [1,2,3], y = abs.(randn(3)), type="scatter", mode = "lines+markers", line = (width = 4,))
        ],
        layout = (title = value,)
    )
end
handle = make_handler(app, debug = true)
HTTP.serve(handle, HTTP.Sockets.localhost, 8080)
```

# Available components
$(ComponentPackages.component_doc_list())
""" Dashboards

const IdProp = Tuple{Symbol, Symbol}

struct CallbackId
    state ::Vector{IdProp}
    input ::Vector{IdProp}
    output ::Vector{IdProp}
end

Base.convert(::Type{Vector{IdProp}}, v::IdProp) = [v]

CallbackId(;input ::Union{Vector{IdProp}, IdProp},  output ::Union{Vector{IdProp}, IdProp}, state ::Union{Vector{IdProp}, IdProp} = Vector{IdProp}()) = CallbackId(state, input, output)

struct Callback
    func ::Function
    id ::CallbackId
end

"""
    struct Dash <: Any

Representation of Dash application
"""
struct Dash
    name ::String
    layout ::Component
    callbacks ::Dict{Symbol, Callback}
    external_stylesheets ::Vector{String}
    url_base_pathname ::String
    assets_folder ::String
    function Dash(name::String, layout::Component; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname="/", assets_folder::String = "assets")
        new(name, layout, Dict{Symbol, Callback}(), external_stylesheets, url_base_pathname, assets_folder)
    end
end
"""
    Dash(layout_maker::Function, name::String; external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname::String="/")::Dash

Construct a Dash app using callback for layout creation

# Arguments
- `layout_maker::Function` - function for layout creation. Must has signature ()::Component
- `name::String` - Dashboard name
- `external_stylesheets::Vector{String} = Vector{String}()` - vector of external css urls
- `url_base_pathname::String="/"` - base url path for dashboard, default "/"
- `assets_folder::String` - a path, relative to the current working directory,
for extra files to be used in the browser. Default `"assets"`

# Examples
```jldoctest
julia> app = Dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
```
"""
function Dash(layout_maker::Function, name::String;  external_stylesheets ::Vector{String} = Vector{String}(), url_base_pathname="/", assets_folder::String = "assets")
    Dash(name, layout_maker(), external_stylesheets=external_stylesheets, url_base_pathname=url_base_pathname, assets_folder = assets_folder)
end

function parse_props(s)
    function make_prop(part)
        m = match(r"^(?<id>[A-Za-z]+[\w\-\:\.]*)\.(?<prop>[A-Za-z]+[\w\-\:\.]*)$", strip(part))
        if isnothing(m)
            error("expected <id>.<property>[,<id>.<property>...] in $(part)")
        end
        return (Symbol(m[:id]), Symbol(m[:prop]))
    end

    props_parts = split(s, ",", keepempty = false)

    return map(props_parts) do part
        return make_prop(part)
    end
end

macro props_str(s)
    return parse_props(s)
end
"""
    @callid_str"

Macro for crating Dash CallbackId.
Parse string in form "[{State1[, ...]}] Input1[, ...] => Output1[, ...]"

#Examples
```julia
    id1 = callid"{inputDiv.children} input.value => output1.value, output2.value"
```
"""
macro callid_str(s)
    rex = r"(\{(?<state>.*)\})?(?<input>.*)=>(?<output>.*)"ms
    m = match(rex, s)
    if isnothing(m)
        error("expected {state} input => output")
    end
    input = parse_props(strip(m[:input]))
    output = parse_props(strip(m[:output]))
    state = isnothing(m[:state]) ? Vector{IdProp}() : parse_props(strip(m[:state]))
    return CallbackId(state, input, output)
end

idprop_string(idprop::IdProp) = "$(idprop[1]).$(idprop[2])"

function output_string(id::CallbackId)
    if length(id.output) == 1
        return idprop_string(id.output[1])
    end
    return ".." *
    join(map(idprop_string, id.output), "...") *
    ".."
end
"""
    callback!(func::Function, app::Dash, id::CallbackId)

Create a callback that updates the output by calling function `func`.

#Examples
```julia
app = Dash("Test") do
    html_div() do
        dcc_input(id="graphTitle", value="Let's Dance!", type = "text"),
        html_div(id="outputID"),
        html_div(id="outputID2")

    end
end
callback!(app, CallbackId(
    state = [(:graphTitle, :type)],
    input = [(:graphTitle, :value)],
    output = [(:outputID, :children), (:outputID2, :children)]
    )
    ) do stateType, inputValue
    return (stateType * "..." * inputValue, inputValue)
end
```
You can use macro `callid` string macro for make CallbackId :
```julia
callback!(app, callid"{graphTitle.type} graphTitle.value => outputID.children, outputID2.children") do stateType, inputValue

    return (stateType * "..." * inputValue, inputValue)
end
```
"""
function callback!(func::Function, app::Dash, id::CallbackId)

    for out in id.output
        if any(x->out in x.id.output, values(app.callbacks))
            error("output \"$(out)\" already registered")
        end
    end

    function check_arr(ids)
        for id in ids
            if !Components.is_valid_idprop(app.layout, id)
                error("Layout hasn't component with id `$(id[1])`` or component with id `$(id[1])` hasn't property `$(id[2])``")
            end
        end
    end

    check_arr(id.state)
    check_arr(id.input)
    check_arr(id.output)

    out_symbol = Symbol(output_string(id))

    push!(app.callbacks, out_symbol => Callback(func, id))
end

function index_page(app::Dash; debug = false)
    metas = ["""<meta http-equiv="X-UA-Compatible" content="IE=edge">""", """<meta charset="UTF-8">"""]
    title = app.name
    css = join(map(app.external_stylesheets) do s
        """<link rel="stylesheet" href="$(s)">"""
    end, "\n"
    )

    app_entry = """

                Loading...

    """
    config = (
        url_base_pathname = "$(app.url_base_pathname)",
        requests_pathname_prefix = "$(app.url_base_pathname)",
        ui = true,
        props_check = debug,
        show_undo_redo = false
    )

    scripts = ComponentPackages.components_js_include(app.url_base_pathname, debug = debug)

    """<!DOCTYPE html>

            $(join(metas, ""))
            <title>$(title)</title>
            $(css)

            $(app_entry)

            $scripts

    """
end

function dependencies_json(app::Dash)
    id_prop_named(p::IdProp) = (id = p[1], property = p[2])
    result = map(values(app.callbacks)) do dep
        (inputs = id_prop_named.(dep.id.input),
        state = id_prop_named.(dep.id.state),
        output = output_string(dep.id)
        )
    end
    return JSON2.write(result)
end

function process_callback(app::Dash, body::String)
    params = JSON2.read(body)
    output = Symbol(params[:output])
    if !haskey(app.callbacks, output)
        return []
    end
    args = []
    if haskey(params, :state)
        append!(args, map(x->haskey(x, :value) ? x[:value] : nothing, params.state))
    end
    if haskey(params, :inputs)
        append!(args, map(x->haskey(x, :value) ? x[:value] : nothing, params.inputs))
    end
    res = app.callbacks[output].func(args...)
    if length(app.callbacks[output].id.output) == 1
        return Dict(
            :response => Dict(
                :props => Dict(
                    Symbol(app.callbacks[output].id.output[1][2]) => res
                )
            )
        )
    end
    response = Dict{Symbol, Any}()
    for (ind, out) in enumerate(app.callbacks[output].id.output)
        push!(response,
        Symbol(out[1]) => Dict(
            Symbol(out[2]) => res[ind]
        )
        )
    end
    return Dict(:response=>response, :multi=>true)

end

function process_assets(app::Dash, path)
    assets_path = "$(app.url_base_pathname)assets/"
    filename = joinpath(app.assets_folder, replace(path, assets_path=>""))
    try
        return HTTP.Response(200, [], body = read(filename))
    catch
        return HTTP.Response(404)
    end
end

"""
    make_handler(app::Dash; debug = false)

Make handler for routing Dash application in HTTP package

#Arguments
- `app::Dash` - Dash application
- `debug::Bool = false` - Enable/disable all the dev tools

#Examples
```jldoctest
julia> app = Dash("Test") do
    html_div() do
        html_h1("Test Dashboard")
    end
end
julia> handler = make_handler(app)
julia> HTTP.serve(handler, HTTP.Sockets.localhost, 8080)
```

"""
function make_handler(app::Dash; debug::Bool = false)
    function (req::HTTP.Request)
        uri = HTTP.URI(req.target)
        ComponentPackages.@register_js_sources(uri.path, app.url_base_pathname)
        if uri.path == "$(app.url_base_pathname)"
            return HTTP.Response(200, index_page(app, debug = debug))
        end
        if uri.path == "$(app.url_base_pathname)_dash-layout"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = JSON2.write(app.layout))
        end
        if uri.path == "$(app.url_base_pathname)_dash-dependencies"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = dependencies_json(app))
        end
        if startswith(uri.path, "$(app.url_base_pathname)assets/")
            return process_assets(app, uri.path)
        end
        if uri.path == "$(app.url_base_pathname)_dash-update-component" && req.method == "POST"
            return HTTP.Response(200, ["Content-Type" => "application/json"],
                body = JSON2.write(
                    process_callback(app, String(req.body))
                )
             )
        end
        return HTTP.Response(404)
    end
end

end # module

max_stars_repo_name: EricForgy/Dashboards

max_issues_repo_name: EricForgy/Dashboards

max_forks_repo_name: EricForgy/Dashboards


text:
__precompile__()
module MobiusUtils

using Dierckx

export bisection, falseposition, interpolate

"""
Simple bisection method to solve nonlinear equation
f(x)=0for f: ℝ → ℝ  on the interval (a,b). The result
is the solution c such that f(c-ϵ/2)×f(c+ϵ/2) < 0, i.e.
is the centre of the interval of the length ϵ over
which f(x) changes the sign.
"""
function bisection(f,a,b,ϵ)
    fa = f(a)
    fb = f(b)
    while abs(a-b)≥ ϵ
        c = (a+b)/2
        fc = f(c)
        (a,fa) = fc*fb ≤ 0? (c,fc) : (a,fa)
        (b,fb) = fa*fc ≤ 0? (c,fc) : (b,fb)
    end
    (a+b)/2
end

"""
Algorithm of false position to solve nonlinear
equation f(x)=0 for f: ℝ → ℝ on the interval x∈(a,b)
with the correction for functions that change from
"almost flat" to "rapid growth" type (see Hamming for
    more details).
"""
function falseposition(f,a,b,ϵ)
    fa = f(a)
    fb = f(b)
    keptabefore = false
    kepta = false
    while abs(a-b)≥ ϵ
        c = a - fa*(b-a)/(fb-fa)
        fc = f(c)
        keptabefore = kepta
        if fc*fb ≤ 0
            a,fa = c,fc
            kepta = false
        end
        if fa*fc ≤ 0
            b,fb = c, fc
            kepta = true
        end
        if keptabefore && kepta
            fa = fa/2
        elseif (!keptabefore) && (!kepta)
            fb = fb/2
        end
    end
    a - fa*(b-a)/(fb-fa)
end

"""
Interpolates the data between the points (x,y). Produces
the function of x.

Parameters:
- **order**: order of interpolation(1-5), defualts to 1
- **outside**: what to do when argument is outside of
  interpolation region: "extrapolate" (default), "error",
  "zero", "nearest"
- **spline**: if `false` (default) returns a function, otherwise
  returns the spline object of Dierckx (can be used as a function
  but is not always recognised as such).
"""
function interpolate(x, y, order = 1, outside = "extrapolate", spline = false)
    spl = Spline1D(x, y, k = order, bc = outside)
    spine? spl : z -> spl(z)
end

end # module


text:
using Plots

t0 = -300:.01:0.01
r = .03
f(r,t0) = log.(2 .- (1+r).^(t0))./log(1+r)
plot(2019 .+ t0, [f(r,t0) f(.07,t0)], labels=["3%", "7%"],xaxis="year X", yaxis="time in the future we will have as much growth as has happened since year X")