text:
import PySCF: pyscf, pyscf_atom_from_xyz, get_4idx
import PySCF: diis

abstract type HartreeFockModel end

struct RestrictedHartreeFock <: HartreeFockModel
   Emol
   e
   C
   #S
   #D
end

function buildFock(hcore,D,nao,eri)
   F = copy(hcore)
   for i in 1:nao, j in 1:nao, # add Gμν*P
       k in 1:nao, l in 1:nao
      ijkl = get_4idx(i,j,k,l)
      ikjl = get_4idx(i,k,j,l)
      F[i,j] += D[k,l]*(2eri[ijkl]-eri[ikjl])
   end
   F
end

# Löwdin symmetric Orthogonalization (with A = X^(-1/2))
function lowdinOrtho(A,B)
   A' * B * A
end

#########################################################################
# Project #3: The Hartree-Fock self-consistent field (SCF) procedure.
#########################################################################
function scf_rhf(mol; verbose=true)
   nelec = mol.nelectron
   nocc  = mol.nelectron ÷ 2
   nao   = mol.nao_nr() |> n -> convert(Int,n) # why not converted ??!
   enuc  = mol.energy_nuc()
   s     = mol.intor("cint1e_ovlp_sph")
   hcore = mol.intor("cint1e_nuc_sph") + mol.intor("cint1e_kin_sph")
   eri   = mol.intor("cint2e_sph", aosym="s8")

   ########################################################
   #4: BUILD THE ORTHOGONALIZATION MATRIX
   ########################################################
   # NOTE: The AO basis is non-orthogonal

   # Symmetric orthogonalization matrix
   s_minushalf = s^(-0.5) # Matrix power (via squaring)

   ########################################################
   #5: BUILD THE INITIAL (GUESS) DENSITY
   ########################################################

   P = zeros(nao,nao) # D = 2P, P often called density matrix

   ########################################################
   #6: COMPUTE THE INITIAL SCF ENERGY
   ########################################################

   Eelec = 0.0

   #*******************************************************
   #*******************************************************
   #
   #             SCF ITERATIVE PROCESS
   #
   #*******************************************************
   #*******************************************************

   iteri     = 0    # extant iteration
   itermax   = 100  # max number of iterations
   δE        = 1e-8 # energy threshold for convergence
   δP        = 1e-8 # P-matrix threshold for convergence
   ΔE        = 1.0; @assert ΔE > δE # initial E difference
   ΔP        = 1.0; @assert ΔP > δP # initial P difference
   Eelec     = 0.0
   oldEelec  = 0.0
   Emol      = 0.0  # total molecular energy
   e         = zeros(nao)
   C         = zeros(nao,nao)
   DIISstack = [] # for DIIS
   Fstack    = [] # for DIIS
   if verbose
      println(" Iter        E(elec)        E(tot)        ΔE(elec)         ΔP")
      println("-----   --------------  -------------  ------------  -------------")
   end

   while abs(ΔE) > δE && iteri < itermax && ΔP > δP

      iteri += 1

      ########################################################
      #7: COMPUTE THE NEW FOCK MATRIX
      ########################################################

      F = buildFock(hcore,P,nao,eri)

      F, Fstack, DIISstack = diis(F, Fstack, DIISstack, s_minushalf, s, P, iteri)

      ########################################################
      #8: BUILD THE NEW DENSITY MATRIX
      ########################################################

      # Transformed Fock matrix
      Fp = lowdinOrtho(s_minushalf,F)
      # MO coeffs as linear combination of orthonormal MO basis functions
      e, Cp = eigen(Hermitian(Fp))
      # New MO Coefficients as linear combination of AO basis functions
      C = s_minushalf * Cp
      # Renew P-matrix.
      P, oldP = zeros(nao,nao), copy(P)

      mocc = C[:,1:nocc]
      P  = mocc * mocc'
      ΔP = norm(P - oldP) # Frobenius norm

      ########################################################
      #9: COMPUTE THE NEW SCF ENERGY
      ########################################################

      Eelec, oldEelec = tr(P * (hcore + F)), Eelec # eq. (238) Janos

      ΔE   = Eelec - oldEelec
      Emol = Eelec + enuc

      if verbose
         printfmt(" {1:3d}  {2:14.8f}  {3:14.8f}  {4:14.8f} {5:14.8f}\n",
                    iteri,   Eelec,       Emol,      ΔE   ,    ΔP)
      end
      ########################################################
      #10: TEST FOR CONVERGENCE: WHILE LOOP (go up)
      ########################################################
   end
   @assert iteri < itermax "NOT CONVERGED!!"

#   covalent_bond_order(P,s,mol)
#
#   rdm1 = fock_dirac_denmat(C,s,2P)
#   nelec_rdm1 = tr(rdm1)
#   @assert isapprox(nelec, nelec_rdm1, atol=1e-8)
#   printfmt("Population: {:.3f} elec \n", nelec_rdm1)
#
#   rdm2 = second_denmat(rdm1)
#   # TODO: energy from density matrices eqs. (224,227) Janos

   RestrictedHartreeFock(Emol,e,C)
end


text:
# QuObj Tests
using Schrodinger
using Test, LinearAlgebra, SparseArrays, Statistics

println("Testing Quantum Info Functions...")

@testset "Quantum Information" begin
# Build a few different variable for testing
N = 13
α = complex(1.1,-0.33)

g = basis(2,0)
e1 = basis(2,1)
coherentα = coherent(N,α)
coherentαhalf = coherent(N,0.5α)
numop = numberop(N)

N! = normalize!

@testset "Expectation Value" begin
    # Test calculation of expectation valu
    @test expect(σz,g) == 1
    @test expect(e1,σz) == -1
    @test expect(σx,N!(g+e1)) ≈ 1
    @test expect(N!(g-e1),σx) ≈ -1
    @test expect(N!(g+1im*e1),σy) ≈ 1
    @test expect(σy,N!(g-1im*e1)) ≈ -1
    @test expect(coherentα,destroy(N)) ≈ α
    @test expect(destroy(N),coherent(N,0.5α,true)) ≈ 0.5α
    @test [expect(basis(N,i),numop) for i=0:N-1] == collect(0:N-1)
    @test expect(N!(basis(N,2)+2*basis(N,3)),numop) == (2+2^2*3)/5
    @test expect(N!(Operator(basis(N,2))+2*Operator(basis(N,3))),numop) == (1*2+2*3)/3
    @test expect(numop,thermal(N,0.21)) ≈ 0.21
    @test expect(projectorop(N,[2,4]),maxmixed(N)) == 2/N
end

ONEQUBIT = [g,e1,N!(g+e1),N!(g-e1),N!(g+1im*e1),N!(g-1im*e1)]
include(joinpath(@__DIR__,"twoqubit.jl"))

@testset "Fidelity" begin # Test calculation of various fidelities
    # fidelity
    @test fidelity(g,g) == 1
    @test fidelity(g,e1) == 0
    @test fidelity(g,N!(g+e1)) == 1/√(2)
    @test fidelity(coherentα,basis(N,3)) ≈ abs(coherentα[4])
    @test fidelity(Operator(coherentα),basis(N,3)) ≈ abs(coherentα[4])
    @test fidelity(basis(N,3),Operator(coherentαhalf)) ≈ abs(coherentαhalf[4])
    @test fidelity(Operator(coherentα),Operator(basis(N,3))) ≈ abs(coherentα[4])
    @test fidelity(Operator(basis(N,2)),Operator(coherent(N,1))) ≈ abs(coherent(N,1)[3])
    for i = 1:20
        P = 10*rand()*randn(N,N) + 10*rand()*im*randn(N,N); P = normalize!(Operator(P'*P))
        Q = 10*rand()*randn(N,N) + 10*rand()*im*randn(N,N); Q = normalize!(Operator(Q'*Q))
        @test fidelity(P,Q) ≈ (rtQ=√(Q); real(trace(√(rtQ*P*rtQ))))
    end
    # fidelity2
    @test fidelity2(g,g) == 1
    @test fidelity2(g,e1) == 0
    @test fidelity2(g,N!(g+e1)) ≈ 1/2
    @test fidelity2(coherentαhalf,basis(N,3)) ≈ abs2(coherentαhalf[4])
    @test fidelity2(Operator(coherentαhalf),basis(N,3)) ≈ abs2(coherentαhalf[4])
    @test fidelity2(basis(N,3),Operator(coherentα)) ≈ abs2(coherentα[4])
    @test fidelity2(Operator(coherentαhalf),Operator(basis(N,3))) ≈ abs2(coherentαhalf[4])
    @test fidelity2(Operator(coherent(N,α,true)),Operator(basis(N,3))) ≈ abs2(coherent(N,α,true)[4])
    @test fidelity2(Operator(basis(N,2)),Operator(coherent(N,1))) ≈ abs2(coherent(N,1)[3])
    # entanglement and average gate fidelity
    U1, V1 = rand_unitary(5), rand_unitary(5)
    @test entanglement_fidelity(U1,V1) ≈ (ε=qeye(5)⊗(U1'V1);ϕ=maxentangled(2,5);abs2(ϕ⋅(ε*ϕ)))
    U2, V2 = rand_unitary(2), rand_unitary(2)
    @test gate_fidelity(U2,V2) ≈ mean(fidelity2(ψ,U2'*(V2*ψ)) for ψ ∈ ONEQUBIT)
    @test gate_fidelity(U2,V2) ≈ (2*entanglement_fidelity(U2,V2)+1)/(2+1)
    U3, V3 = rand_unitary(4,(2,2)), rand_unitary(4,(2,2))
    @test gate_fidelity(U3,V3) ≈ mean(fidelity2(ψ,U3'*(V3*ψ)) for ψ ∈ TWOQUBIT)
    @test gate_fidelity(U3,V3) ≈ (4*entanglement_fidelity(U3,V3)+1)/(4+1)
end

@testset "Purity" begin
    U = rand_unitary(N)
    ρ = N!(Operator(coherentαhalf)+numop)
    @test purity(g) == purity(e1) == 1
    @test purity(coherentα) == purity(coherentαhalf) == 1
    @test purity(ρ) ≈ inner(ρ,ρ)
    @test purity(U*ρ*U') ≈ purity(ρ)
    @test purity(maxmixed(4)) == 1/4
end

@testset "Level Probabilities" begin
    # Test calculation of level probabilities
    @test levelprobs(g) == [1,0]
    @test levelprobs(e1) == [0,1]
    @test levelprobs(N!(e1+g)) ≈ [0.5,0.5]
    @test levelprobs(e1⊗g) == [0,0,1,0]
    @test levelprobs(N!(e1⊗g+g⊗e1)) ≈ [0,0.5,0.5,0]
    @test levelprobs(N!(e1⊗g+g⊗g),1:2) ≈ [[0.5,0.5],[1,0]]
end

function inner_cs_test(pe=0.1,ϕ=1.23)
    # σx on qubit 1 of 2-qubit system
    # no phase difference
    V1 = ((1-pe)*σx + pe*qeye(2)) .⊗ (Operator(basis(2,0)), Operator(basis(2,1)))
    # with phase difference
    V2 = (V1[1], cis(ϕ)*V1[2])
    # Ideal gate
    U = σx⊗Operator(basis(2,0)) + σx⊗Operator(basis(2,1))
    # inner products
    w = sqrt(sum(abs2,diag(sum(V2)'U)))/2
    x = abs(inner(sum(data.(V1)),data(U)))/4
    y = sum(abs,inner.(data.(V2),(data(U),)))/4
    z = abs(inner(sum(data.(V2)),data(U)))/4
    w,x,y,z
end

function anharm_NOT_RF(steps)
    dragpulse = Schrodinger.dragpulse_cos
    α = -0.4 * 2π # anharmonicity (GHz)
    Hd = α*Operator(basis(3,2)) # drift Hamiltonian
    σX, σY = 1/2*(create(3)+destroy(3)), 1im/2*(create(3)-destroy(3))
    λ = √2
    T = 5 # (ns)
    tspan = (-T/2,T/2)
    paramsx, paramsy = [0,T,0,α,λ,0,π], [0,T,0,α,λ,π/2,π]
    V = Operator(SchrodingerProp(Hd,((σX,dragpulse,paramsx),(σY,dragpulse,paramsy)),tspan,steps))
    #U = Operator([0 1 0; 1 0 0; 0 0 1]) # 3lvl NOT gate
    ## Fidelity
    #ϕ = N!(basis(3,0)⊗basis(3,0) + basis(3,1)⊗basis(3,1))
    #ϕ' * (qeye(3)⊗U'V) * Operator(ϕ) * (qeye(3)⊗V'U) * ϕ
    # Create objective function type
    #O = CoherentSubspaces(Ut,1:2,Hd,Hc,T,n) # care only about computational subspace
    #grape(O,ui)
    return V
end

end

ext: jl


text:
# All functionality corresponding to the FNN-statistic as Loss function Γ.

## Constructors:

"""
    FNN_statistic <: AbstractLoss

    Constructor for the FNN-statistic loss function (false nearest neighbor) based
    on Hegger & Kantz [^Hegger1999].

    ## Fieldnames
    * `threshold::Float`: A threshold for the tolerable cumulative FNN decrease
      for the current embedding. When the fraction of FNNs fall below this threshold
      in an embedding cycle the embedding stops.
    * `r::Float = 2`: The FNN-distance-expansion threshold (typically set to 2).
    * `samplesize::Real = 1.`: the fraction of all phase space points
      to be considered in the computation of the L-statistic(s).

    ## Defaults
    * When calling `FNN_statistic()`, a FNN_statistic-object is created, which uses no
      threshold and uses the FNN-inter threshold `r=2`.
    * When calling `FNN_statistic(threshold)`, a FNN_statistic-object is created, which uses
      the given `threshold` and uses the FNN-inter threshold `r=2`.

    [^Hegger1999]: Hegger & Kantz, [Improved false nearest neighbor method to detect determinism in time series data. Physical Review E 60, 4970](https://doi.org/10.1103/PhysRevE.60.4970).
"""
struct FNN_statistic <: AbstractLoss
    threshold::AbstractFloat
    r::AbstractFloat
    samplesize::Real
    # Constraints and Defaults
    FNN_statistic(x=0.,y=2.,z=1.) = begin
        @assert x >= 0 "Threshold for FNN-statistic must be ≥ 0"
        @assert y > 0 "FNN-distance-expansion threshold must be >0"
        @assert 0 < z ≤ 1. "The samplesize must be in the interval (0 1]."
        typeof(x) <: Int ? new(convert(AbstractFloat, x),y,z) : new(x,y,z)
        typeof(y) <: Int ? new(x,convert(AbstractFloat,y),z) : new(x,y,z)
    end
end

## Functions:

"""
    Return the loss based on the FNN-statistic `FNN` and indices `max_idx`  for all local maxima in dps
"""
function compute_loss(Γ::FNN_statistic, Λ::AbstractDelayPreselection, dps::Vector{P}, Y_act::Dataset{D, T}, Ys, τs, w::Int, ts::Int, τ_vals, ts_vals; metric=Euclidean(), kwargs...) where {P, D, T}

    r = Γ.r
    s = Ys[:,ts]
    samplesize = Γ.samplesize

    max_idx = get_max_idx(Λ, dps, τ_vals, ts_vals, ts) # get the candidate delays
    isempty(max_idx) && return Float64[], Int64[], []

    FNN_trials = zeros(Float64, length(max_idx))

    # compute nearest-neighbor-distances for actual trajectory
    Y_act2 = Y_act[1:end-τs[maximum(max_idx)-1],:]
    Y_act2 = DelayEmbeddings.standardize(Y_act2)
    vtree = KDTree(Y_act2, metric)

    NN = length(Y_act2)
    if samplesize==1
        ns = 1:NN
        Nfp = length(ns)
    else
        Nfp = Int(floor(samplesize*NN)) # number of considered fiducial points
        ns = sample(vec(1:NN), Nfp, replace = false)  # indices of fiducial points
    end

    _, NNdist_old = DelayEmbeddings.all_neighbors(vtree, Y_act2[ns], ns, 1, w)

    for (i,τ_idx) in enumerate(max_idx)
        # create candidate phase space vector for this peak/τ-value
        Y_trial = DelayEmbeddings.hcat_lagged_values(Y_act,s,τs[τ_idx-1])
        Y_trial = DelayEmbeddings.standardize(Y_trial)
        vtree = KDTree(Y_trial, metric)
        _, NNdist_new = DelayEmbeddings.all_neighbors(vtree, Y_trial[ns], ns, 1, w)
        # compute FNN-statistic
        FNN_trials[i] = fnn_embedding_cycle(NNdist_old, NNdist_new[1:length(NNdist_old)], r)
    end
    return FNN_trials, max_idx, [nothing for i in max_idx]
end

"""
    fnn_embedding_cycle(NNdist, NNdistnew, r=2) -> FNNs

    Compute the amount of false nearest neighbors `FNNs` [^Hegger1999], when adding another component
    to a given (vector-) time series. This new component is the `τ`-lagged version
    of a univariate time series. `NNdist` is storing the distances of the nearest
    neighbor for all considered fiducial points and `NNdistnew` is storing the
    distances of the nearest neighbor for each fiducial point in one embedding
    dimension higher using a given `τ`. The obligatory threshold `r` is by default
    set to 2.

"""
function fnn_embedding_cycle(NNdist, NNdistnew, r::Real=2)
    @assert length(NNdist) == length(NNdistnew) "Both input vectors need to store the same number of distances."
    N = length(NNdist)
    fnns = 0
    fnns2= 0
    inverse_r = 1/r
    @inbounds for i = 1:N
        if NNdistnew[i][1]/NNdist[i][1] > r && NNdist[i][1] < inverse_r
            fnns +=1
        end
        if NNdist[i][1] < inverse_r
            fnns2 +=1
        end
    end
    if fnns==0
        return 1
    else
        return fnns/fnns2
    end
end

max_stars_repo_path: src/Cost functions/FNN_statistic.jl

max_issues_repo_path: src/Cost functions/FNN_statistic.jl

max_forks_repo_path: src/Cost functions/FNN_statistic.jl


text:
# Place here configuration options that will be set for all environments
ENV["GENIE_ENV"] = "prod"

max_stars_repo_name: StatisticalMice/performance-tests

max_issues_repo_name: StatisticalMice/performance-tests

max_forks_repo_name: StatisticalMice/performance-tests


text:
"""
    ReliableInformationDistribution

Distribution based on Schneidman et al 2011, where only "reliable" patterns are
used to influence the coefficients in the max-ent distribution.

"""
mutable struct ReliableInformationDistribution
    X::BitMatrix
    I::Vector{Int}
    metadata::Dict{Any,Any}
    cache::Dict{Any,Any}

    function ReliableInformationDistribution(X::Union{Matrix{Bool},BitMatrix}, I=1:size(X,1); kwargs)
        # MAYBEDO implement this if it seems necessary
    end
end

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
using LibGit2, InteractiveUtils

mutable struct Server
    storedir::String
    context::Pkg.Types.Context
    depot::Dict
end

abstract type SymStore end
struct ModuleStore <: SymStore
    name::VarRef
    vals::Dict{Symbol,Any}
    doc::String
    exported::Bool
    exportednames::Vector{Symbol}
    used_modules::Vector{Symbol}
end

ModuleStore(m) = ModuleStore(VarRef(m), Dict{Symbol,Any}(), "", true, unsorted_names(m), Symbol[])
Base.getindex(m::ModuleStore, k) = m.vals[k]
Base.setindex!(m::ModuleStore, v, k) = (m.vals[k] = v)
Base.haskey(m::ModuleStore, k) = haskey(m.vals, k)

const EnvStore = Dict{Symbol,ModuleStore}

struct Package
    name::String
    val::ModuleStore
    uuid::Base.UUID
    sha::Union{Vector{UInt8},Nothing}
end
Package(name::String, val::ModuleStore, uuid::String, sha) = Package(name, val, Base.UUID(uuid), sha)

struct MethodStore
    name::Symbol
    mod::Symbol
    file::String
    line::Int32
    sig::Vector{Pair{Any,Any}}
    kws::Vector{Symbol}
    rt::Any
end

struct DataTypeStore <: SymStore
    name::FakeTypeName
    super::FakeTypeName
    parameters::Vector{Any}
    types::Vector{Any}
    fieldnames::Vector{Any}
    methods::Vector{MethodStore}
    doc::String
    exported::Bool
end

function DataTypeStore(@nospecialize(t), parent_mod, exported)
    ur_t = Base.unwrap_unionall(t)
    parameters = if isdefined(ur_t, :parameters)
        map(ur_t.parameters) do p
            _parameter(p)
        end
    else
        []
    end
    types = if isdefined(ur_t, :types)
        map(ur_t.types) do p
            FakeTypeName(p)
        end
    else
        []
    end
    DataTypeStore(FakeTypeName(ur_t), FakeTypeName(ur_t.super), parameters, types, ur_t.isconcretetype && fieldcount(ur_t) > 0 ? collect(fieldnames(ur_t)) : Symbol[], MethodStore[], _doc(t), exported)
end

struct FunctionStore <: SymStore
    name::VarRef
    methods::Vector{MethodStore}
    doc::String
    extends::VarRef
    exported::Bool
end

function FunctionStore(@nospecialize(f), parent_mod, exported)
    if f isa Core.IntrinsicFunction
        FunctionStore(VarRef(VarRef(Core.Intrinsics), nameof(f)), MethodStore[], _doc(f), VarRef(VarRef(parentmodule(f)), nameof(f)), exported)
    else
        FunctionStore(VarRef(VarRef(parent_mod), nameof(f)), MethodStore[], _doc(f), VarRef(VarRef(parentmodule(f)), nameof(f)), exported)
    end
end

struct GenericStore <: SymStore
    name::VarRef
    typ::Any
    doc::String
    exported::Bool
end

# adapted from https://github.com/timholy/CodeTracking.jl/blob/afc73a957f5034cc7f02e084a91283c47882f92b/src/utils.jl#L87-L122

"""
    path = maybe_fix_path(path)

Return a normalized, absolute path for a source file `path`.
"""
function maybe_fix_path(file)
    if !isabspath(file)
        # This may be a Base or Core method
        newfile = Base.find_source_file(file)
        if isa(newfile, AbstractString)
            file = normpath(newfile)
        end
    end
    return maybe_fixup_stdlib_path(file)
end

safe_isfile(x) = try isfile(x); catch; false end
const BUILDBOT_STDLIB_PATH = dirname(abspath(joinpath(String((@which versioninfo()).file), "..", "..", "..")))
replace_buildbot_stdlibpath(str::String) = replace(str, BUILDBOT_STDLIB_PATH => Sys.STDLIB)
"""
    path = maybe_fixup_stdlib_path(path::String)

Return `path` corrected for julia issue [#26314](https://github.com/JuliaLang/julia/issues/26314) if applicable.
Otherwise, return the input `path` unchanged.

Due to the issue mentioned above, location info for methods defined one of Julia's standard libraries
are, for non source Julia builds, given as absolute paths on the worker that built the `julia` executable.
This function corrects such a path to instead refer to the local path on the users drive.
"""
function maybe_fixup_stdlib_path(path)
    if !safe_isfile(path)
        maybe_stdlib_path = replace_buildbot_stdlibpath(path)
        safe_isfile(maybe_stdlib_path) && return maybe_stdlib_path
    end
    return path
end

const _global_method_cache = IdDict{Any,Vector{Any}}()
function methodinfo(@nospecialize(f); types = Tuple, world = typemax(UInt))
    key = (f, types, world)
    cached = get(_global_method_cache, key, nothing)
    if cached === nothing
        cached = Base._methods(f, types, -1, world)
        _global_method_cache[key] = cached
    end
    return cached
end

function methodlist(@nospecialize(f))
    ms = methodinfo(f)
    Method[x[3]::Method for x in ms]
end

function sparam_syms(meth::Method)
    s = Symbol[]
    sig = meth.sig
    while sig isa UnionAll
        push!(s, Symbol(sig.var.name))
        sig = sig.body
    end
    return s
end

function cache_methods(@nospecialize(f), name, env, get_return_type)
    if isa(f, Core.Builtin)
        return MethodStore[]
    end
    types = Tuple
    world = typemax(UInt)
    ms = Tuple{Module,MethodStore}[]
    methods0 = try
        methodinfo(f; types = types, world = world)
    catch err
        return ms
    end
    ind_of_method_w_kws = Int[] # stores the index of methods with kws.
    i = 1
    for m in methods0
        # Get inferred method return type
        if get_return_type
            sparams = Core.svec(sparam_syms(m[3])...)
            rt = try
                @static if isdefined(Core.Compiler, :NativeInterpreter)
                Core.Compiler.typeinf_type(Core.Compiler.NativeInterpreter(), m[3], m[3].sig, sparams)
            else
                Core.Compiler.typeinf_type(m[3], m[3].sig, sparams, Core.Compiler.Params(world))
            end
            catch e
                Any
            end
        else
            rt = Any
        end
        file = maybe_fix_path(String(m[3].file))
        MS = MethodStore(m[3].name, nameof(m[3].module), file, m[3].line, [], Symbol[], FakeTypeName(rt))
        # Get signature
        sig = Base.unwrap_unionall(m[1])
        argnames = getargnames(m[3])
        for i = 2:m[3].nargs
            push!(MS.sig, argnames[i] => FakeTypeName(sig.parameters[i]))
        end
        kws = getkws(m[3])
        if !isempty(kws)
            push!(ind_of_method_w_kws, i)
        end
        for kw in kws
            push!(MS.kws, kw)
        end
        push!(ms, (m[3].module, MS))
        i += 1
    end
    # Go back and add kws to methods defined in the same place as others with kws.
    for i in ind_of_method_w_kws
        for j = 1:length(ms) # only need to go up to `i`?
            if ms[j][2].file == ms[i][2].file && ms[j][2].line == ms[i][2].line && isempty(ms[j][2].kws)
                for kw in ms[i][2].kws
                    push!(ms[j][2].kws, kw)
                end
            end
        end
    end

    func_vr = VarRef(VarRef(parentmodule(f)), name)
    for i = 1:length(ms)
        mvr = VarRef(ms[i][1])
        modstore = _lookup(mvr, env)
        if modstore !== nothing
            if !haskey(modstore, name)
                modstore[name] = FunctionStore(VarRef(mvr, name), MethodStore[ms[i][2]], "", func_vr, false)
            elseif !(modstore[name] isa DataTypeStore || modstore[name] isa FunctionStore)
                modstore[name] = FunctionStore(VarRef(mvr, name), MethodStore[ms[i][2]], "", func_vr, false)
            else
                push!(modstore[name].methods, ms[i][2])
            end
        else
        end
    end
    return ms
end

getargnames(m::Method) = Base.method_argnames(m)
@static if length(first(methods(Base.kwarg_decl)).sig.parameters) == 2
    getkws = Base.kwarg_decl
else
    function getkws(m::Method)
        sig = Base.unwrap_unionall(m.sig)
        length(sig.parameters) == 0 && return []
        sig.parameters[1] isa Union && return []
        !isdefined(Base.unwrap_unionall(sig.parameters[1]), :name) && return []
        fname = Base.unwrap_unionall(sig.parameters[1]).name
        if isdefined(fname.mt, :kwsorter)
            Base.kwarg_decl(m, typeof(fname.mt.kwsorter))
        else
            []
        end
    end
end

function apply_to_everything(f, m = nothing, visited = Base.IdSet{Module}())
    if m isa Module
        push!(visited, m)
        for s in unsorted_names(m, all = true, imported = true)
            (!isdefined(m, s) || s == nameof(m)) && continue
            x = getfield(m, s)
            f(x)
            if x isa Module && !in(x, visited)
                apply_to_everything(f, x, visited)
            end
        end
    else
        for m in Base.loaded_modules_array()
            in(m, visited) || apply_to_everything(f, m, visited)
        end
    end
end

function oneverything(f, m = nothing, visited = Base.IdSet{Module}())
    if m isa Module
        push!(visited, m)
        state = nothing
        for s in unsorted_names(m, all = true, imported = true)
            !isdefined(m, s) && continue
            x = getfield(m, s)
            state = f(m, s, x, state)
            if x isa Module && !in(x, visited)
                oneverything(f, x, visited)
            end
        end
    else
        for m in Base.loaded_modules_array()
            in(m, visited) || oneverything(f, m, visited)
        end
    end
end

const _global_symbol_cache_by_mod = IdDict{Module,Base.IdSet{Symbol}}()
function build_namecache(m, s, @nospecialize(x), state::Union{Base.IdSet{Symbol},Nothing} = nothing)
    if state === nothing
        state = get(_global_symbol_cache_by_mod, m, nothing)
        if state === nothing
            state = _global_symbol_cache_by_mod[m] = Base.IdSet{Symbol}()
        end
    end
    push!(state, s)
end

function getnames(m::Module)
    cache = get(_global_symbol_cache_by_mod, m, nothing)
    if cache === nothing
        oneverything(build_namecache, m)
        cache = _global_symbol_cache_by_mod[m]
    end
    return cache
end

function allmodulenames()
    symbols = Base.IdSet{Symbol}()
    oneverything((m, s, x, state) -> (x isa Module && push!(symbols, s); return state))
    return symbols
end

function allthingswithmethods()
    symbols = Base.IdSet{Any}()
    oneverything(function (m, s, x, state)
        if !Base.isvarargtype(x) && !isempty(methodlist(x))
            push!(symbols, x)
        end
        return state
    end)
    return symbols
end

function allmethods()
    ms = Method[]
    oneverything(function (m, s, x, state)
        if !Base.isvarargtype(x) && !isempty(methodlist(x))
            append!(ms, methodlist(x))
        end
        return state
    end)
    return ms
end

usedby(outer, inner) = outer !== inner && isdefined(outer, nameof(inner)) && getproperty(outer, nameof(inner)) === inner && all(isdefined(outer, name) || !isdefined(inner, name) for name in unsorted_names(inner))
istoplevelmodule(m) = parentmodule(m) === m || parentmodule(m) === Main

function getmoduletree(m::Module, amn, visited = Base.IdSet{Module}())
    push!(visited, m)
    cache = ModuleStore(m)
    for s in unsorted_names(m, all = true, imported = true)
        !isdefined(m, s) && continue
        x = getfield(m, s)
        if x isa Module
            if istoplevelmodule(x)
                cache[s] = VarRef(x)
            elseif m === parentmodule(x)
                cache[s] = getmoduletree(x, amn, visited)
            else
                cache[s] = VarRef(x)
            end
        end
    end
    for n in amn
        if n !== nameof(m) && isdefined(m, n)
            x = getfield(m, n)
            if x isa Module
                if !haskey(cache, n)
                    cache[n] = VarRef(x)
                end
                if x !== Main && usedby(m, x)
                    push!(cache.used_modules, n)
                end
            end
        end
    end
    cache
end

function getenvtree(names = nothing)
    amn = allmodulenames()
    EnvStore(nameof(m) => getmoduletree(m, amn) for m in Base.loaded_modules_array() if names === nothing || nameof(m) in names)
end

# faster and more correct split_module_names
all_names(m) = all_names(m, x -> isdefined(m, x))
function all_names(m, pred, symbols = Set(Symbol[]), seen = Set(Module[]))
    push!(seen, m)
    ns = unsorted_names(m; all = true, imported = false)
    for n in ns
        isdefined(m, n) || continue
        Base.isdeprecated(m, n) && continue
        val = getfield(m, n)
        if val isa Module && !(val in seen)
            all_names(val, pred, symbols, seen)
        end
        if pred(n)
            push!(symbols, n)
        end
    end
    symbols
end

function symbols(env::EnvStore, m::Union{Module,Nothing} = nothing, allnames::Base.IdSet{Symbol} = getallns(), visited = Base.IdSet{Module}();  get_return_type = false)
    if m isa Module
        cache = _lookup(VarRef(m), env, true)
        cache === nothing && return
        push!(visited, m)
        ns = all_names(m)
        for s in ns
            !isdefined(m, s) && continue
            x = getfield(m, s)
            if Base.unwrap_unionall(x) isa DataType # Unions aren't handled here.
                if parentmodule((x)) === m
                    cache[s] = DataTypeStore(x, m, s in getnames(m))
                    cache_methods(x, s, env, get_return_type)
                elseif nameof(x) !== s
                    # This needs some finessing.
                    cache[s] = DataTypeStore(x, m, s in getnames(m))
                    cache_methods(x, s, env, get_return_type)
                else
                    # These are imported variables that are reexported.
                    cache[s] = VarRef(VarRef(parentmodule(x)), nameof(x))
                end
            elseif x isa Function
                if parentmodule(x) === m || (x isa Core.IntrinsicFunction && m === Core.Intrinsics)
                    cache[s] = FunctionStore(x, m, s in getnames(m))
                    cache_methods(x, s, env, get_return_type)
                elseif !haskey(cache, s)
                    # This will be replaced at a later point by a FunctionStore if methods for `x` are defined within `m`.
                    if x isa Core.IntrinsicFunction
                        cache[s] = VarRef(VarRef(Core.Intrinsics), nameof(x))
                    else
                        cache[s] = VarRef(VarRef(parentmodule(x)), nameof(x))
                    end
                elseif !((cache[s] isa FunctionStore || cache[s] isa DataTypeStore) && !isempty(cache[s].methods))
                    # These are imported variables that are reexported.
                    # We don't want to remove Func/DT stores that have methods (these will be specific to the module)
                    if x isa Core.IntrinsicFunction
                        cache[s] = VarRef(VarRef(Core.Intrinsics), nameof(x))
                    else
                        cache[s] = VarRef(VarRef(parentmodule(x)), nameof(x))
                    end
                end
            elseif x isa Module
                if x === m
                    cache[s] = VarRef(x)
                elseif parentmodule(x) === m
                    symbols(env, x, allnames, visited)
                else
                    cache[s] = VarRef(x)
                end
            else
                cache[s] = GenericStore(VarRef(VarRef(m), s), FakeTypeName(typeof(x)), _doc(x), s in getnames(m))
            end
        end
    else
        for m in Base.loaded_modules_array()
            in(m, visited) || symbols(env, m, allnames, visited)
        end
    end
end

function load_core(; get_return_type = false)
    c = Pkg.Types.Context()
    cache = getenvtree([:Core,:Base])
    symbols(cache)
    cache[:Main] = ModuleStore(VarRef(nothing, :Main), Dict(), "", true, [], [])

    # This is wrong. As per the docs the Base.include each module should have it's own
    # version.
    push!(cache[:Base].exportednames, :include)

    # Add special cases for built-ins
    let f = cache[:Base][:include]
        cache[:Base][:include] = FunctionStore(f.name, cache[:Base][:MainInclude][:include].methods, f.doc, f.extends, true)
    end

    cache[:Base][Symbol("@.")] = cache[:Base][Symbol("@__dot__")]
    cache[:Core][:Main] = GenericStore(VarRef(nothing, :Main), FakeTypeName(Module), _doc(Main), true)
    # Add built-ins
    builtins = Symbol[nameof(getfield(Core, n).instance) for n in unsorted_names(Core, all = true) if isdefined(Core, n) && getfield(Core, n) isa DataType && isdefined(getfield(Core, n), :instance) && getfield(Core, n).instance isa Core.Builtin]
    cnames = unsorted_names(Core)
    for f in builtins
        if !haskey(cache[:Core], f)
            cache[:Core][f] = FunctionStore(getfield(Core, Symbol(f)), Core, Symbol(f) in cnames)
        end
    end
    haskey(cache[:Core], :_typevar) && push!(cache[:Core][:_typevar].methods, MethodStore(:_typevar, :Core, "built-in", 0, [:n => FakeTypeName(Symbol), :lb => FakeTypeName(Any), :ub => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:_apply].methods, MethodStore(:_apply, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core].vals, :_apply_iterate) && push!(cache[:Core][:_apply_iterate].methods, MethodStore(:_apply_iterate, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    if isdefined(Core, :_call_latest)
        push!(cache[:Core][:_call_latest].methods, MethodStore(:_call_latest, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
        push!(cache[:Core][:_call_in_world].methods, MethodStore(:_call_in_world, :Core, "built-in", 0, [:world => FakeTypeName(UInt), :f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    else
        if isdefined(Core, :_apply_in_world)
            push!(cache[:Core][:_apply_in_world].methods, MethodStore(:_apply_in_world, :Core, "built-in", 0, [:world => FakeTypeName(UInt), :f => FakeTypeName(Function), :args => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
        end
        push!(cache[:Core][:_apply_latest].methods, MethodStore(:_apply_latest, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    end
    push!(cache[:Core][:_apply_pure].methods, MethodStore(:_apply_pure, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:_expr].methods, MethodStore(:_expr, :Core, "built-in", 0, [:head => FakeTypeName(Symbol), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Expr)))
    haskey(cache[:Core].vals, :_typevar) && push!(cache[:Core][:_typevar].methods, MethodStore(:_typevar, :Core, "built-in", 0, [:name => FakeTypeName(Symbol), :lb => FakeTypeName(Any), :ub => FakeTypeName(Any)], Symbol[], FakeTypeName(TypeVar)))
    push!(cache[:Core][:applicable].methods, MethodStore(:applicable, :Core, "built-in", 0, [:f => FakeTypeName(Function), :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Bool)))
    push!(cache[:Core][:apply_type].methods, MethodStore(:apply_type, :Core, "built-in", 0, [:T => FakeTypeName(UnionAll), :types => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(UnionAll)))
    push!(cache[:Core][:arrayref].methods, MethodStore(:arrayref, :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any), :c => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:arrayset].methods, MethodStore(:arrayset, :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any), :c => FakeTypeName(Any), :d => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:arraysize].methods, MethodStore(:arraysize, :Core, "built-in", 0, [:a => FakeTypeName(Array), :i => FakeTypeName(Int)], Symbol[], FakeTypeName(Int)))
    haskey(cache[:Core], :const_arrayref) && push!(cache[:Core][:const_arrayref].methods, MethodStore(:const_arrayref, :Core, "built-in", 0, [:args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:fieldtype].methods, MethodStore(:fieldtype, :Core, "built-in", 0, [:t => FakeTypeName(DataType), :field => FakeTypeName(Symbol)], Symbol[], FakeTypeName(Type{T} where T)))
    push!(cache[:Core][:getfield].methods, MethodStore(:setfield, :Core, "built-in", 0, [:object => FakeTypeName(Any), :item => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:ifelse].methods, MethodStore(:ifelse, :Core, "built-in", 0, [:condition => FakeTypeName(Bool), :x => FakeTypeName(Any), :y => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:invoke].methods, MethodStore(:invoke, :Core, "built-in", 0, [:f => FakeTypeName(Function), :x => FakeTypeName(Any), :argtypes => FakeTypeName(Type{T} where T) , :args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:isa].methods, MethodStore(:isa, :Core, "built-in", 0, [:a => FakeTypeName(Any), :T => FakeTypeName(Type{T} where T)], Symbol[], FakeTypeName(Bool)))
    push!(cache[:Core][:isdefined].methods, MethodStore(:getproperty, :Core, "built-in", 0, [:value => FakeTypeName(Any), :field => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:nfields].methods, MethodStore(:nfields, :Core, "built-in", 0, [:x => FakeTypeName(Any)], Symbol[], FakeTypeName(Int)))
    push!(cache[:Core][:setfield!].methods, MethodStore(:setfield!, :Core, "built-in", 0, [:value => FakeTypeName(Any), :name => FakeTypeName(Symbol), :x => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:sizeof].methods, MethodStore(:sizeof, :Core, "built-in", 0, [:obj => FakeTypeName(Any)], Symbol[], FakeTypeName(Int)))
    push!(cache[:Core][:svec].methods, MethodStore(:svec, :Core, "built-in", 0, [:args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:throw].methods, MethodStore(:throw, :Core, "built-in", 0, [:e => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:tuple].methods, MethodStore(:tuple, :Core, "built-in", 0, [:args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:typeassert].methods, MethodStore(:typeassert, :Core, "built-in", 0, [:x => FakeTypeName(Any), :T => FakeTypeName(Type{T} where T)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:typeof].methods, MethodStore(:typeof, :Core, "built-in", 0, [:x => FakeTypeName(Any)], Symbol[], FakeTypeName(Type{T} where T)))

    push!(cache[:Core][:getproperty].methods, MethodStore(:getproperty, :Core, "built-in", 0, [:value => FakeTypeName(Any), :name => FakeTypeName(Symbol)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:setproperty!].methods, MethodStore(:setproperty!, :Core, "built-in", 0, [:value => FakeTypeName(Any), :name => FakeTypeName(Symbol), :x => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:setproperty!].methods, MethodStore(:setproperty!, :Core, "built-in", 0, [:value => FakeTypeName(Any), :name => FakeTypeName(Symbol), :x => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_abstracttype) && push!(cache[:Core][:_abstracttype].methods, MethodStore(:_abstracttype, :Core, "built-in", 0, [:m => FakeTypeName(Module), :x => FakeTypeName(Symbol), :p => FakeTypeName(Core.SimpleVector)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_primitivetype) && push!(cache[:Core][:_primitivetype].methods, MethodStore(:_primitivetype, :Core, "built-in", 0, [:m => FakeTypeName(Module), :x => FakeTypeName(Symbol), :p => FakeTypeName(Core.SimpleVector), :n => FakeTypeName(Core.Int)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_equiv_typedef) && push!(cache[:Core][:_equiv_typedef].methods, MethodStore(:_equiv_typedef, :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_setsuper!) && push!(cache[:Core][:_setsuper!].methods, MethodStore(:_setsuper!, :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_structtype) && push!(cache[:Core][:_structtype].methods, MethodStore(:_structtype, :Core, "built-in", 0, [:m => FakeTypeName(Module), :x => FakeTypeName(Symbol), :p => FakeTypeName(Core.SimpleVector), :fields => FakeTypeName(Core.SimpleVector), :mut => FakeTypeName(Bool), :z => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    haskey(cache[:Core], :_typebody) && push!(cache[:Core][:_typebody!].methods, MethodStore(:_typebody!, :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:(===)].methods, MethodStore(:(===), :Core, "built-in", 0, [:a => FakeTypeName(Any), :b => FakeTypeName(Any)], Symbol[], FakeTypeName(Any)))
    push!(cache[:Core][:(<:)].methods, MethodStore(:(<:), :Core, "built-in", 0, [:a => FakeTypeName(Type{T} where T), :b => FakeTypeName(Type{T} where T)], Symbol[], FakeTypeName(Any)))

    for bi in builtins
        if haskey(cache[:Core], bi) && isempty(cache[:Core][bi].methods)
            # Add at least one arbitrary method for anything left over
            push!(cache[:Core][bi].methods, MethodStore(bi, :none, "built-in", 0, [:x => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)))
        end
    end

    cache[:Core][:ccall] = FunctionStore(VarRef(VarRef(Core), :ccall),
        MethodStore[
            MethodStore(:ccall, :Core, "built-in", 0, [:args => FakeTypeName(Vararg{Any})], Symbol[], FakeTypeName(Any)) # General method - should be fixed
        ],
        "`ccall((function_name, library), returntype, (argtype1, ...), argvalue1, ...)`\n`ccall(function_name, returntype, (argtype1, ...), argvalue1, ...)`\n`ccall(function_pointer, returntype, (argtype1, ...), argvalue1, ...)`\n\nCall a function in a C-exported shared library, specified by the tuple (`function_name`, `library`), where each component is either a string or symbol. Instead of specifying a library, one\ncan also use a `function_name` symbol or string, which is resolved in the current process. Alternatively, `ccall` may also be used to call a function pointer `function_pointer`, such as one\nreturned by `dlsym`.\n\nNote that the argument type tuple must be a literal tuple, and not a tuple-valued variable or expression.\n\nEach `argvalue` to the `ccall` will be converted to the corresponding `argtype`, by automatic insertion of calls to `unsafe_convert(argtype, cconvert(argtype, argvalue))`. (See also the documentation for `unsafe_convert` and `cconvert` for further details.) In most cases, this simply results in a call to `convert(argtype, argvalue)`.",
        VarRef(VarRef(Core), :ccall),
        true)
    push!(cache[:Core].exportednames, :ccall)
    cache[:Core][Symbol("@__doc__")] = FunctionStore(VarRef(VarRef(Core), Symbol("@__doc__")), [], "", VarRef(VarRef(Core), Symbol("@__doc__")), true)
    cache_methods(getfield(Core, Symbol("@__doc__")), Symbol("@__doc__"), cache, false)
    # Accounts for the dd situation where Base.rand only has methods from Random which doesn't appear to be explicitly used.
    # append!(cache[:Base][:rand].methods, cache_methods(Base.rand, cache))
    for m in cache_methods(Base.rand, :rand, cache, get_return_type)
        push!(cache[:Base][:rand].methods, m[2])
    end

    # Intrinsics
    cache[:Core][:add_int] = VarRef(VarRef(VarRef(nothing, :Core), :Intrinsics), :add_int)
    cache[:Core][:sle_int] = VarRef(VarRef(VarRef(nothing, :Core), :Intrinsics), :sle_int)
    return cache
end

function collect_extended_methods(depot::EnvStore, extendeds = Dict{VarRef,Vector{VarRef}}())
    for m in depot
        collect_extended_methods(m[2], extendeds, m[2].name)
    end
    extendeds
end

function collect_extended_methods(mod::ModuleStore, extendeds, mname)
    for (n, v) in mod.vals
        if (v isa FunctionStore) && v.extends != v.name
            haskey(extendeds, v.extends) ? push!(extendeds[v.extends], mname) : (extendeds[v.extends] = VarRef[v.extends.parent, mname])
        elseif v isa ModuleStore
            collect_extended_methods(v, extendeds, v.name)
        end
    end
end

getallns() = let allns = Base.IdSet{Symbol}(); oneverything((m, s, x, state) -> push!(allns, s)); allns end

"""
    split_module_names(m::Module, allns)

Return two lists of names accessible from calling getfield(m, somename)`. The first
contains those symbols returned by `Base.names(m, all = true)`. The second contains
all others, including imported symbols and those introduced by the `using` of modules.
"""
function split_module_names(m::Module, allns)
    internal_names = getnames(m)
    availablenames = Set{Symbol}([s for s in allns if isdefined(m, s)])
    # usinged_names = Set{Symbol}()

    for n in availablenames
        if (n in internal_names)
            pop!(availablenames, n)
        end
    end
    allms = get_all_modules()
    for u in get_used_modules(m, allms)
        for n in unsorted_names(u)
            if n in availablenames
                pop!(availablenames, n)
                # push!(usinged_names, pop!(availablenames, n))
            end
        end
    end
    internal_names, availablenames
end

get_all_modules() = let allms = Base.IdSet{Module}(); apply_to_everything(x -> if x isa Module push!(allms, x) end); allms end
get_used_modules(M, allms = get_all_modules()) = [m for m in allms if usedby(M, m)]


text:
module RPC
typealias SlotNum Int
typealias Hostname AbstractString
typealias Port Int64
# Base Message Type, not too useful
abstract BaseMessage
# Command type; cause some action to happen on a server
abstract Command <: BaseMessage
# Resonse as a result of causing a command to happen
abstract Response <: BaseMessage
# Commands that do not need replication
abstract NonCriticalCommand <: Command
# Commands that take a user defined function and run them, without consensus
abstract ExternalNonCriticalCommand <: NonCriticalCommand
abstract InternalNonCriticalCommand <: NonCriticalCommand
# Critical commands need fault tolerance, not instance specific
abstract CriticalCommand <: Command
# Internal critical commands for actions such as membership changes
abstract InternalCriticalCommand <: CriticalCommand
# External critical commands for the user subclass
abstract ExternalCriticalCommand <: CriticalCommand
# Non Critical for consistency, but critical for progress/consensus forming messages.
abstract ConsensusCommand <: Command
abstract ConsensusResponse <: Response
#####################################
##    Components of the log        ##
#####################################
@enum Fate Decided Pending Forgotten
type Slot
    num::SlotNum
    n_p::Int64
    n_a::Int64
    v_a::Nullable{CriticalCommand}
    fate::Fate
end

#####################################
#   Core Command Section            #
#####################################

type PrepareArgs <: ConsensusCommand
    slot_num::SlotNum
    n::Int64
end

type PrepareReply <: ConsensusResponse
    slot_num::SlotNum
    Ok::Bool
    n::Int64
    n_a::Int64
    v_a::Nullable{CriticalCommand}
end
type AcceptArgs <: ConsensusCommand
    slot_num::SlotNum
    n::Int64
    v::CriticalCommand
end
type AcceptReply <: ConsensusResponse
    slot_num::SlotNum
    Ok::Bool
    n::Int64
    v::CriticalCommand
end
type DecideArgs <: ConsensusCommand
    slot_num::SlotNum
    v::CriticalCommand
end
type DecideReply <: ConsensusResponse
    slot_num::SlotNum
    Ok::Bool
end

#########################
## Other msg Types     ##
#########################
type Ping <: InternalNonCriticalCommand
    msg::ASCIIString
end
type PingReply <: Response
end

# TODO
type DeadLeader <: InternalCriticalCommand
    id::Int64
end

type AddPeer <: InternalCriticalCommand
    new::Tuple{Hostname, Port}
end
type DeletePeer <: InternalCriticalCommand
    dead_peer::Tuple{Hostname, Port}
end
type SwapPeer <: InternalCriticalCommand
    old::DeletePeer
    new::AddPeer
end
type NoOp <: InternalCriticalCommand
end
type ExtraInfo
    completed::Array{Int64, 1}
end

immutable WithExtraInfo <: BaseMessage
    command::BaseMessage
    info::ExtraInfo
end
end

max_stars_repo_name: JeremyRubin/Rock


text:
using DynamicalSystemsBase

ti = time()

# Systems:
include("dynsys_types.jl")
include("dynsys_tangent.jl")
include("dynsys_inference.jl")
include("continuous_systems.jl")
include("discrete_systems.jl")
include("integrators_with_callbacks.jl")
include("norm_tests.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits=3), " seconds or ", round(ti/60, digits=3), " minutes")


text:
module C
const NAMES = [
    "InhCDK"
    "CycBCdk1T"
    "CycACdk2T"
    "PP1T"
    "kapp1"
    "kapp1a"
    "kipp1"
    "kipp1C"
    "kPP1Gw"
    "ENSAtot"
    "B55tot"
    "SubT"
    "kass"
    "kdis"
    "kcatB55"
    "kGwENSA"
    "kppxGwl"
    "kcBc1Sub"
    "kcBc1G"
    "Gwtot"
    "kB55G"
    "kB55Sub"
    "kcAc2G"
    "k25S"
    "k25F"
    "kweeS"
    "kweeF"
    "kcBc1W1"
    "kcBc125"
    "kppxY15"
    "kcAc2W1"
    "kcAc225"
    "kB55W1"
    "kB5525"
    "Kd"
]

#name2idx
for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
struct SynthesisPlan{T, P1, P2}
    planθ::P1
    planφ::P2
    C::ColumnPermutation
    temp::Vector{T}
end

function plan_synthesis(A::Matrix{T}) where T<:fftwNumber
    m, n = size(A)
    x = FFTW.FakeArray(T, m)
    y = FFTW.FakeArray(T, n)
    planθ = FFTW.plan_r2r!(x, FFTW.REDFT01), FFTW.plan_r2r!(x, FFTW.RODFT01)
    planφ = FFTW.plan_r2r!(y, FFTW.HC2R)
    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
    SynthesisPlan(planθ, planφ, C, zeros(T, n))
end

function plan_synthesis2(A::Matrix{T}) where T<:fftwNumber
    m, n = size(A)
    x = FFTW.FakeArray(T, m)
    y = FFTW.FakeArray(T, n)
    planθ = FFTW.plan_r2r!(x, FFTW.REDFT00), FFTW.plan_r2r!(FFTW.FakeArray(T, m-2), FFTW.RODFT00)
    planφ = FFTW.plan_r2r!(y, FFTW.HC2R)
    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
    SynthesisPlan(planθ, planφ, C, zeros(T, n))
end

struct AnalysisPlan{T, P1, P2}
    planθ::P1
    planφ::P2
    C::ColumnPermutation
    temp::Vector{T}
end

function plan_analysis(A::Matrix{T}) where T<:fftwNumber
    m, n = size(A)
    x = FFTW.FakeArray(T, m)
    y = FFTW.FakeArray(T, n)
    planθ = FFTW.plan_r2r!(x, FFTW.REDFT10), FFTW.plan_r2r!(x, FFTW.RODFT10)
    planφ = FFTW.plan_r2r!(y, FFTW.R2HC)
    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
    AnalysisPlan(planθ, planφ, C, zeros(T, n))
end

function plan_analysis2(A::Matrix{T}) where T<:fftwNumber
    m, n = size(A)
    x = FFTW.FakeArray(T, m)
    y = FFTW.FakeArray(T, n)
    planθ = FFTW.plan_r2r!(x, FFTW.REDFT00), FFTW.plan_r2r!(FFTW.FakeArray(T, m-2), FFTW.RODFT00)
    planφ = FFTW.plan_r2r!(y, FFTW.R2HC)
    C = ColumnPermutation(vcat(1:2:n, 2:2:n))
    AnalysisPlan(planθ, planφ, C, zeros(T, n))
end

function LAmul!(Y::Matrix{T}, P::SynthesisPlan{T, Tuple{r2rFFTWPlan{T,(FFTW.REDFT01,),true,1}, r2rFFTWPlan{T,(FFTW.RODFT01,),true,1}}}, X::Matrix{T}) where T
    M, N = size(X)

    # Column synthesis
    PCe = P.planθ[1]
    PCo = P.planθ[2]

    X[1] *= two(T)
    mul_col_J!(Y, PCe, X, 1)
    X[1] *= half(T)

    for J = 2:4:N
        mul_col_J!(Y, PCo, X, J)
        J < N && mul_col_J!(Y, PCo, X, J+1)
    end
    for J = 4:4:N
        X[1,J] *= two(T)
        J < N && (X[1,J+1] *= two(T))
        mul_col_J!(Y, PCe, X, J)
        J < N && mul_col_J!(Y, PCe, X, J+1)
        X[1,J] *= half(T)
        J < N && (X[1,J+1] *= half(T))
    end
    lmul!(half(T), Y)

    # Row synthesis
    lmul!(inv(sqrt(π)), Y)
    invsqrttwo = inv(sqrt(2))
    @inbounds for i = 1:M Y[i] *= invsqrttwo end

    temp = P.temp
    planφ = P.planφ
    C = P.C
    for I = 1:M
        copy_row_I!(temp, Y, I)
        row_synthesis!(planφ, C, temp)
        copy_row_I!(Y, temp, I)
    end
    Y
end

function LAmul!(Y::Matrix{T}, P::SynthesisPlan{T, Tuple{r2rFFTWPlan{T,(FFTW.REDFT00,),true,1}, r2rFFTWPlan{T,(FFTW.RODFT00,),true,1}}}, X::Matrix{T}) where T
    M, N = size(X)

    # Column synthesis
    PCe = P.planθ[1]
    PCo = P.planθ[2]

    X[1] *= two(T)
    X[M,1] *= two(T)
    mul_col_J!(Y, PCe, X, 1)
    X[1] *= half(T)
    X[M,1] *= half(T)

    for J = 2:4:N
        mul_col_J!(Y, PCo, X, J, false)
        J < N && mul_col_J!(Y, PCo, X, J+1, false)
    end
    for J = 4:4:N
        X[1,J] *= two(T)
        X[M,J] *= two(T)
        J < N && (X[1,J+1] *= two(T); X[M,J+1] *= two(T))
        mul_col_J!(Y, PCe, X, J)
        J < N && mul_col_J!(Y, PCe, X, J+1)
        X[1,J] *= half(T)
        X[M,J] *= half(T)
        J < N && (X[1,J+1] *= half(T); X[M,J+1] *= half(T))
    end
    lmul!(half(T), Y)

    # Row synthesis
    lmul!(inv(sqrt(π)), Y)
    invsqrttwo = inv(sqrt(2))
    @inbounds for i = 1:M Y[i] *= invsqrttwo end

    temp = P.temp
    planφ = P.planφ
    C = P.C
    for I = 1:M
        copy_row_I!(temp, Y, I)
        row_synthesis!(planφ, C, temp)
        copy_row_I!(Y, temp, I)
    end
    Y
end

function LAmul!(Y::Matrix{T}, P::AnalysisPlan{T, Tuple{r2rFFTWPlan{T,(FFTW.REDFT10,),true,1}, r2rFFTWPlan{T,(FFTW.RODFT10,),true,1}}}, X::Matrix{T}) where T
    M, N = size(X)

    # Row analysis
    temp = P.temp
    planφ = P.planφ
    C = P.C
    for I = 1:M
        copy_row_I!(temp, X, I)
        row_analysis!(planφ, C, temp)
        copy_row_I!(Y, temp, I)
    end

    # Column analysis
    PCe = P.planθ[1]
    PCo = P.planθ[2]

    mul_col_J!(Y, PCe, Y, 1)
    Y[1] *= half(T)
    for J = 2:4:N
        mul_col_J!(Y, PCo, Y, J)
        J < N && mul_col_J!(Y, PCo, Y, J+1)
    end
    for J = 4:4:N
        mul_col_J!(Y, PCe, Y, J)
        J < N && mul_col_J!(Y, PCe, Y, J+1)
        Y[1,J] *= half(T)
        J < N && (Y[1,J+1] *= half(T))
    end
    lmul!(sqrt(π)*inv(T(M)), Y)
    sqrttwo = sqrt(2)
    @inbounds for i = 1:M Y[i] *= sqrttwo end

    Y
end

function LAmul!(Y::Matrix{T}, P::AnalysisPlan{T, Tuple{r2rFFTWPlan{T,(FFTW.REDFT00,),true,1}, r2rFFTWPlan{T,(FFTW.RODFT00,),true,1}}}, X::Matrix{T}) where T
    M, N = size(X)

    # Row analysis
    temp = P.temp
    planφ = P.planφ
    C = P.C
    for I = 1:M
        copy_row_I!(temp, X, I)
        row_analysis!(planφ, C, temp)
        copy_row_I!(Y, temp, I)
    end

    # Column analysis
    PCe = P.planθ[1]
    PCo = P.planθ[2]

    mul_col_J!(Y, PCe, Y, 1)
    Y[1] *= half(T)
    Y[M, 1] *= half(T)
    for J = 2:4:N
        mul_col_J!(Y, PCo, Y, J, true)
        J < N && mul_col_J!(Y, PCo, Y, J+1, true)
        Y[M-1,J] = zero(T)
        J < N && (Y[M-1,J+1] = zero(T))
    end
    for J = 4:4:N
        mul_col_J!(Y, PCe, Y, J)
        J < N && mul_col_J!(Y, PCe, Y, J+1)
        Y[1,J] *= half(T)
        Y[M,J] *= half(T)
        J < N && (Y[1,J+1] *= half(T); Y[M,J+1] *= half(T))
    end
    lmul!(sqrt(π)*inv(T(M-1)), Y)
    sqrttwo = sqrt(2)
    @inbounds for i = 1:M Y[i] *= sqrttwo end

    Y
end

function row_analysis!(P, C, vals::Vector{T}) where T
    n = length(vals)
    cfs = lmul!(two(T)/n,P*vals)
    cfs[1] *= half(T)
    if iseven(n)
        cfs[n÷2+1] *= half(T)
    end

    negateeven!(reverseeven!(lmul!(C, cfs)))
end

function row_synthesis!(P, C, cfs::Vector{T}) where T
    n = length(cfs)
    lmul!(C', reverseeven!(negateeven!(cfs)))
    if iseven(n)
        cfs[n÷2+1] *= two(T)
    end
    cfs[1] *= two(T)
    P*lmul!(half(T), cfs)
end

function copy_row_I!(temp::Vector, Y::Matrix, I::Int)
    M, N = size(Y)
    @inbounds @simd for j = 1:N
        temp[j] = Y[I+M*(j-1)]
    end
    temp
end

function copy_row_I!(Y::Matrix, temp::Vector, I::Int)
    M, N = size(Y)
    @inbounds @simd for j = 1:N
        Y[I+M*(j-1)] = temp[j]
    end
    Y
end

function reverseeven!(x::Vector)
    n = length(x)
    if iseven(n)
        @inbounds @simd for k=2:2:n÷2
            x[k], x[n+2-k] = x[n+2-k], x[k]
        end
    else
        @inbounds @simd for k=2:2:n÷2
            x[k], x[n+1-k] = x[n+1-k], x[k]
        end
    end
    x
end

function negateeven!(x::Vector)
    @inbounds @simd for k = 2:2:length(x)
        x[k] *= -1
    end
    x
end

function mul_col_J!(Y::Matrix{T}, P::r2rFFTWPlan{T}, X::Matrix{T}, J::Int) where T
    unsafe_execute_col_J!(P, X, Y, J)
    return Y
end

function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int) where T<:fftwDouble
    M = size(X, 1)
    ccall((:fftw_execute_r2r, libfftw), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
end

function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int) where T<:fftwSingle
    M = size(X, 1)
    ccall((:fftwf_execute_r2r, libfftwf), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
end

function mul_col_J!(Y::Matrix{T}, P::r2rFFTWPlan{T}, X::Matrix{T}, J::Int, TF::Bool) where T
    unsafe_execute_col_J!(P, X, Y, J, TF)
    return Y
end

function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int, TF::Bool) where T<:fftwDouble
    M = size(X, 1)
    if TF
        ccall((:fftw_execute_r2r, libfftw), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+2), pointer(Y, M*(J-1)+1))
    else
        ccall((:fftw_execute_r2r, libfftw), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+2))
    end
end

function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int, TF::Bool) where T<:fftwSingle
    M = size(X, 1)
    if TF
        ccall((:fftwf_execute_r2r, libfftwf), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+2), pointer(Y, M*(J-1)+1))
    else
        ccall((:fftwf_execute_r2r, libfftwf), Nothing, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+2))
    end
end


text:
using MLFS
using PrettyPrint
using Setfield
using MLFS.HM
using MLFS.TypedIO
using MLStyle
using JSON

stmts = @surf_toplevel begin
    op_Getter :: implicit[Fn[field[(i64, :a, i64)], i64]]
    op_Getter = @julia (identity)

    f :: implicit[field[(i64, :a, i64)]]
    f = @julia (x -> x)
    x :: implicit[i64]
    x = 1

    g :: (i64 -> i64 -> _)
    g = x -> y -> (y.a).?(aa)
end

mr = Surf.ModuleRecord(:a, [], stmts)
open("a.mlfs", "w") do f
    write(f, JSON.json(toVec(Surf.ModuleRecord, mr)))
end

smlfsCompile(["a.mlfs"], String[], "a")

ir2jl("./", "a_gen_asm.jl")


text:
using Distributions, StatsBase, StatsPlots
using Plots
using Polynomials
using Random
using DelimitedFiles
using KernelDensity
using LaTeXStrings
include("ToolBox.jl")
include("predictionGSBR.jl")
include("./data/datasets.jl")

# generate some data from chaotic logistic map
nf = 210
ntrain = 200
θ = [1.,0.,-1.71]
x₀ = .5
data = noisypolynomial(x₀, θ, noisemix2; n=nf, seed=11) # already "seeded"
plot(data)
### set up sampler parameters for the GSBR model
# degree of modeling polynomial, max number of iterations, burnin period,
#         sampler seet, prediction horizon
xtrain = data[2:ntrain]
ytest = data[ntrain+1:end]
deg,maxiter,burnin,seed,T = 5, 40000, 5000, 1, 10
@time thetas, zp, sx0, pred, clusters = predictionGSBR(x,deg,maxiter,burnin,seed,T,"/sims/logistic")

# estimates of coefficients

thetahat = mean(thetas, dims=1)
histogram(sx0)

zpl, zpr = -0.5, 0.5;
zp = zp[zp.>zpl];
zp = zp[zp.<zpr];

noiseplot = plot(kde(zp), color=:red, lw=1.5, label=L"\mathrm{est }\hat{f}(z)", grid=:false) #histogram(zp, bins=200, normalize=:pdf, color=:lavender, label="predictive");
x_range = range(zpl, stop=zpr, length=120);
ddnoise = noisemixdensity2.(x_range);
plot!(noiseplot, x_range, ddnoise, color=:black, lw=1.5, label=L"\mathrm{true } f(z)", ylim=(0,30), ylabel="pdf")

#std prediction plots
stdplt = scatter(data, colour = :blue, label = "Data", ylim=(-1.5, 2), grid=:false)
plot!(stdplt, [ntrain], seriestype =:vline, colour = :green, linestyle =:dash, label = "Training Data End")

tsteps = 1:nf
fit = zeros(length(xtrain))
fit[1] = polyMap(thetahat, xtrain[1])
for i in 2:length(fit)
    fit[i] = polyMap(thetahat, fit[i-1])
end
plot!(stdplt, tsteps[2:ntrain], fit, alpha=0.4, colour =:blue, label = "fitted model")

allpredictions = hcat(pred...)
thinnedpredictions = allpredictions[1:10:end,:]

plot!(stdplt, tsteps[ntrain+1:end], mean(thinnedpredictions, dims=1)', color=:purple, ribbon=std(thinnedpredictions,dims=1)', alpha=0.4, label="prediction")

bestplt = scatter(tsteps[ntrain+1:end], ytest, colour = :blue, label = "Data", ylim=(-1.5, 2), grid=:false)
plot!(bestplt, tsteps[ntrain+1:end], mean(thinnedpredictions, dims=1)', color=:purple, ribbon=std(thinnedpredictions,dims=1)', alpha=0.4, label="prediction")

mean((mean(thinnedpredictions, dims=1)' .- ytest).^2)

max_stars_repo_name: cmerkatas/GSBR

max_issues_repo_name: cmerkatas/GSBR

max_forks_repo_name: cmerkatas/GSBR


text:
using CxxWrap
using ngraph_jll

@info "Building Lib"

# Paths for linking
cxxhome = CxxWrap.prefix_path()
juliahome = dirname(Base.Sys.BINDIR)
ngraphhome = dirname(dirname(ngraph_jll.libngraph_path))

# Use Clang since it seems to get along better with Julia
cxx = "clang++"

cxxflags = [
    "-g",
    "-O3",
    "-Wall",
    "-fPIC",
    "-std=c++17",
    "-DPCM_SILENT",
    "-DJULIA_ENABLE_THREADING",
    "-Dexcept_EXPORTS",
    # Surpress some warnings from Cxx
    "-Wno-unused-variable",
    "-Wno-unused-lambda-capture",
]

includes = [
    "-I$(joinpath(cxxhome, "include"))",
    "-I$(joinpath(juliahome, "include", "julia"))",
    "-I$(joinpath(ngraphhome, "include"))",
]

loadflags = [
    # Linking flags for Julia
    "-L$(joinpath(juliahome, "lib"))",
    "-Wl,--export-dynamic",
    "-Wl,-rpath,$(joinpath(juliahome, "lib"))",
    "-ljulia",
    # Linking Flags for CxxWrap
    "-L$(joinpath(cxxhome, "lib"))",
    "-Wl,-rpath,$(joinpath(cxxhome, "lib"))",
    "-lcxxwrap_julia",
    # Linking Flags for nGraph
    "-L$(joinpath(ngraphhome, "lib"))",
    "-Wl,-rpath,$(joinpath(ngraphhome, "lib"))",
    "-lngraph",
]

src = joinpath(@__DIR__, "ngraph-julia.cpp")
so = joinpath(@__DIR__, "libngraph-julia.so")

cmd = `$cxx $cxxflags $includes -shared $src -lpthread -o $so $loadflags`
@show cmd
run(cmd)


text:
##
using ENN.Neurons, ENN.TimePetriNets, GLMakie, GraphMakie, DataStructures

function combine(all_neurons, all_connections)
    neurons = merge(all_neurons...)
    connections = all_connections[1]
    for cons in all_connections[2:end]
        for (key, con) in pairs(cons)
            if key in keys(connections)
                append!(connections[key], con)
            else
                connections[key] = con
            end
        end
    end
    (; neurons, connections)
end

function prune!(neurons, connections)
    all_targets = Tuple{Symbol,Symbol}[]
    for connection in values(connections)
        append!(all_targets, connection)
    end
    unique!(all_targets)

    for (key,neuron) in pairs(neurons)
        for segment in neuron.all_segments
            to_delete = [i for (i,inp) in enumerate(segment.exc_inputs) if (key, inp) ∉ all_targets]
            deleteat!(segment.exc_inputs, to_delete)

            to_delete = [i for (i,inp) in enumerate(segment.inh_inputs) if (key, inp) ∉ all_targets]
            deleteat!(segment.inh_inputs, to_delete)
        end
    end
end

function make_up_counter(name, size, init_name=Symbol("$(name)_init"), inc_name=Symbol("$(name)_inc"))
    neurons = OrderedDict((Symbol("$(name)_$(i)") => Neuron{Int}((((), [:input1, :input2], [:inh], :dendrite),), [:inc], [], :soma; axon_delay=2) for i in 1:size)...)
    connections = OrderedDict(
        (Symbol("$(name)_$(i)") => [(Symbol("$(name)_$(i+1)"), :input1)] for i in 1:size-1)...,
        init_name => [(Symbol("$(name)_1"), :input1)],
        inc_name => [(key, :inc) for key in keys(neurons)]
    )
    (; neurons, connections)
end

function make_loop(name, loop_label, loop_position, loop_origin, condition_name=Symbol("$(name)_condition"))
    neurons = OrderedDict(Symbol("$(name)_loop_$(loop_label)") => Neuron{Int}((((), [:condition], [], :dendrite),), [:trigger], [], :soma; spike_duration=1))
    connections = OrderedDict(
        Symbol("$(name)_loop_$(loop_label)") => [(Symbol("$(name)_$(loop_origin)"), :input2), (Symbol("$(name)_$(loop_position+1)"), :inh)],
        Symbol("$(name)_$(loop_position)") => [(Symbol("$(name)_loop_$(loop_label)"),:trigger)],
        condition_name => [(Symbol("$(name)_loop_$(loop_label)"), :condition)]
    )
    (; neurons, connections)
end

function make_up_down_counter(name, size, init_name=Symbol("$(name)_init"), inc_name=Symbol("$(name)_inc"), dec_name=Symbol("$(name)_dec"), refresh_name=Symbol("$(name)_refresh"))
    neurons = OrderedDict(
        Symbol("$(name)_$(:interneuron)")=>Neuron{Int}((),[:input_dec, :input_inc, :input_refresh, :input_init],[],:soma; axon_delay=1),
        (Symbol("$(name)_$(i)") => Neuron{Int}(
        (
            # left dendrite branch
            (
                (
                # left top branch
                    ((), [:input_121], [:input_122], :dendrite_12),
                ),
                [:input_111], [:input_112],
                :dendrite_11
            ),
            # center dendrite branch
            (
                (
                # left top branch
                    ((), [:input_321], [:input_322], :dendrite_32),
                ),
                [:input_311], [:input_312],
                :dendrite_31
            ),
            # right dendrite branch
            (
                (
                # right top branch
                    ((), i==1 ? [:input_221,:input_223] : [:input_221], [:input_222], :dendrite_22),
                ),
                [:input_211], [:input_212],
                :dendrite_21
            )
            # soma
        ), i==1 ? [:input_01, :input_02, :input_03, :input_04] : [:input_01, :input_02, :input_03], [], :soma; axon_delay=10)
                           for i in 1:size)...
    )

    connections = OrderedDict(
        # from i to i+1, left top dendrite branch

        (Symbol("$(name)_$(i)") => [
            (Symbol("$(name)_$(i+1)"), :input_121),
            (Symbol("$(name)_$(i)"), :input_321),
            (Symbol("$(name)_$(i-1)"), :input_221)
        ] for i in 2:size-1)...,
        dec_name => [
            [(Symbol("$(name)_$(i)"), :input_02) for i in 1:size-1];
            [(Symbol("$(name)_$(i)"), :input_211) for i in 1:size-1];
             (Symbol("$(name)_$(:interneuron)"), :input_dec)
        ],
        inc_name => [
            [(Symbol("$(name)_$(i)"), :input_01) for i in 2:size];
            [(Symbol("$(name)_$(i)"), :input_111) for i in 2:size];
             (Symbol("$(name)_1"), :input_01);
             (Symbol("$(name)_$(:interneuron)"), :input_inc)
        ],
        refresh_name => [
            [(Symbol("$(name)_$(i)"), :input_03) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_311) for i in 1:size];
             (Symbol("$(name)_$(:interneuron)"), :input_refresh)
        ],
        Symbol("$(name)_$(:interneuron)") => [
            [(Symbol("$(name)_$(i)"), :input_222) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_212) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_122) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_112) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_322) for i in 1:size];
            [(Symbol("$(name)_$(i)"), :input_312) for i in 1:size]
        ],
        init_name => [
            (Symbol("$(name)_1"), :input_04),
            (Symbol("$(name)_1"), :input_111),
            (Symbol("$(name)_1"), :input_121),
            (Symbol("$(name)_$(:interneuron)"), :input_init)
        ],
    )
    if size>1
        push!(connections,
            Symbol("$(name)_1") => [
                (Symbol("$(name)_2"), :input_121),
                (Symbol("$(name)_1"), :input_321),
                (Symbol("$(name)_1"), :input_223)
            ],
            Symbol("$(name)_$(size)") => [
                (Symbol("$(name)_$(size)"), :input_321),
                (Symbol("$(name)_$(size-1)"), :input_221)
            ]
        )
    else
        push!(connections,
            Symbol("$(name)_1") => [
                (Symbol("$(name)_1"), :input_321),
                (Symbol("$(name)_1"), :input_223)
            ]
        )
    end
    (; neurons, connections)
end

function make_relay_neuron(name, sources)
    neurons = OrderedDict(name => Neuron{Int}((), [Symbol("input_$(src)") for src in sources], [], :soma))

    connections = OrderedDict(
        (src => [(name, Symbol("input_$(src)"))] for src in sources)...
    )
    (; neurons, connections)
end

function make_input_trigger(netplot, input_name)
    net_tpn = netplot[:tpn][]
    state = netplot[:state]
    t_id = net_tpn.T_index[Symbol("input_$(input_name)_start")]

    function trigger()
        fire!(net_tpn, state[], t_id)
        fire_necessary!(net_tpn, state[])
        notify(state)
    end
end

function make_time_advance(netplot, trigger_period)
    net_tpn = netplot[:tpn][]
    state = netplot[:state]

    function advance_trigger()
        sim_time = 0
        # execute all transitions until next trigger time
        while true
            τ = max_delay(net_tpn, state[])
            if sim_time + τ > trigger_period
                break
            end
            sim_time += τ
            elapse!(net_tpn, state[], τ)
            fire_necessary!(net_tpn, state[])
        end

        # advance to next trigger time
        elapse!(net_tpn, state[], trigger_period - sim_time)
        notify(state)
    end
end

##

prog_counter=make_up_counter(:program_counter, 5, :init);
prog_loop = make_loop(:program_counter, "loop_1", 3, 1, :condition);
counter1 = make_up_down_counter(:counter1, 3, :init);
relay_inc1 = make_relay_neuron(:counter1_inc, [:program_counter_2, :program_counter_3])
relay_dec1 = make_relay_neuron(:counter1_dec, Symbol[:program_counter_4,:program_counter_5])
relay_refresh1 = make_relay_neuron(:counter1_refresh, Symbol[:program_counter_1])
prog = combine([prog_counter.neurons, prog_loop.neurons, relay_inc1.neurons, relay_dec1.neurons, relay_refresh1.neurons, counter1.neurons], [prog_counter.connections, counter1.connections, prog_loop.connections, relay_inc1.connections, relay_dec1.connections, relay_refresh1.connections]);
prune!(prog.neurons, prog.connections)

net = NeuralNetwork(
    prog.neurons,
    [:init,:program_counter_inc, :condition],
    Symbol[],
    prog.connections
);

f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
res=netplot!(ax, net)
display(f)
trigger_init = make_input_trigger(res, :init)
trigger_condition = make_input_trigger(res, :condition)
trigger_inc = make_input_trigger(res, :program_counter_inc)
advance_cycle = make_time_advance(res, 75)

## Run through without looping
trigger_init()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()

## Run through with looping
trigger_init()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
trigger_condition()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()

##
net = NeuralNetwork(
    res.neurons,
    [:up1_init,:up1_inc,:up1_dec],
    Symbol[],
    res.connections
);

##

net = NeuralNetwork(
    program_counter_neurons,
    [:init,:trigger],
    Symbol[],
    counter_ff_connections
)
##
function parse_program(expr::Expr, counters...)
    neurons = []
    connections = []
    groups = Vector{Symbol}[]

    counter_sizes = OrderedDict(counters...)

    counter_ops = OrderedDict()
    for key in keys(counter_sizes)
        counter_ops[key] = (inc=Int[],dec=Int[])
    end
    # go through program and create loop neurons as we go
    i=0
    for instruction in expr.args
        if (instruction isa Expr) && (instruction.head == :call)
            i+=1
            if(length(instruction.args)==3 && instruction.args[1] == :JZ)
                (counter, loopback) = instruction.args[2:3]
                (n,c) = make_loop(:program_counter,Symbol("inst_$(i)_loop"),i+1,loopback,Symbol("counter_$(counter)_1"))
                if !haskey(counter_ops, counter)
                    counter_ops[counter] = (inc=Int[],dec=Int[])
                end
                push!(neurons,n)
                push!(connections,c)
            elseif (length(instruction.args)==2 && instruction.args[1] == :INC)
                counter = instruction.args[2]
                if !haskey(counter_ops, counter)
                    counter_ops[counter] = (inc=Int[],dec=Int[])
                end
                push!(counter_ops[counter].inc, i)
            elseif (length(instruction.args)==2 && instruction.args[1] == :DEC)
                counter = instruction.args[2]
                if !haskey(counter_ops, counter)
                    counter_ops[counter] = (inc=Int[],dec=Int[])
                end
                push!(counter_ops[counter].dec, i)
            else
                error("Cannot parse instruction '$(instruction)'.")
            end
        end
    end
    num_instructions = i+2 # +1 for the start instruction, +1 to have a buffer if the program ends in a JZ
    # make program counter
    (n,c) = make_up_counter(:program_counter, num_instructions, :init, :step)
    pushfirst!(neurons, n)
    pushfirst!(connections, c)

    prog = combine(neurons, connections)
    push!(groups, collect(keys(prog.neurons)))

    # make up-down counters & relay neurons
    for (counter,ops) in pairs(counter_ops)
        sz = get(counter_sizes, counter, length(ops.inc))+1
        cnt = make_up_down_counter(Symbol("counter_$(counter)"), sz, :init)
        rly_inc = make_relay_neuron(Symbol("counter_$(counter)_inc"),[Symbol("program_counter_$(i)") for i in ops.inc])
        rly_dec = make_relay_neuron(Symbol("counter_$(counter)_dec"),[Symbol("program_counter_$(i)") for i in ops.dec])
        rly_refresh = make_relay_neuron(Symbol("counter_$(counter)_refresh"),[Symbol("program_counter_$(i)") for i in 1:num_instructions if ((i ∉ ops.inc) && (i∉ ops.dec))])

        push!(neurons, cnt.neurons, rly_inc.neurons, rly_dec.neurons, rly_refresh.neurons)
        push!(connections, cnt.connections, rly_inc.connections, rly_dec.connections, rly_refresh.connections)
        push!(groups, collect(keys(cnt.neurons)))
        append!(groups[1],collect(keys(rly_inc.neurons)),collect(keys(rly_dec.neurons)),collect(keys(rly_refresh.neurons)))
    end

    prog = combine(neurons, connections)
    prune!(prog.neurons, prog.connections)

    net=NeuralNetwork(
            (prog.neurons),
            [:init,:step],
            Symbol[],
            (prog.connections)
    )
    return net,groups
end

macro program(expr, counters...)
    net,groups = parse_program(expr, eval.(counters)...)
    return :($net,$groups)
end

net,groups = @program begin
    INC(1)
    INC(1)
    INC(1)
    INC(1)
    INC(1)
    INC(2)
    INC(2)
    INC(2)
    INC(2)
    INC(2)
    DEC(2)
    DEC(2)
    DEC(2)
    DEC(2)
    DEC(2)
    DEC(1)
    DEC(1)
    DEC(1)
    DEC(1)
    DEC(1)
end 1=>5 2=>5

net,groups = @program begin
    INC(1)
    INC(1)
    INC(1)
    INC(1)
    INC(1)
    INC(2)
    INC(2)
    INC(2)
    JZ(2,13)
    DEC(2)
    DEC(1)
    JZ(3,9)
    DEC(3)
end 1=>5 2=>5 3=>1

begin
f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
res=netplot!(ax, net, groups)
hidedecorations!(ax)
display(f)
trigger_init = make_input_trigger(res, :init)
trigger_inc = make_input_trigger(res, :step)
advance_cycle = make_time_advance(res, 75)
advance_minimal = make_time_advance(res, 1)
end

## Run through without looping
begin
sleep(2)
trigger_init()
for i in 1:75
    sleep(0.001)
    advance_minimal()
end
for j in 1:15
    trigger_inc()
    for i in 1:75
        sleep(0.001)
        advance_minimal()
    end
end
end
advance_cycle()
trigger_inc()
## Run through with looping
trigger_init()
advance_cycle()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
trigger_condition()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()
trigger_inc()
advance_cycle()

##
f = Figure()
ax = Axis(f[1, 1], aspect=DataAspect())
res=netplot!(ax, net)
net_tpn = res[:tpn][]
state = res[:state]
t_init_id = res[:tpn][].T_index[:input_init_start]
t_trigger_id = res[:tpn][].T_index[:input_trigger_start]
display(f)

trigger_period = 50
last_trigger = 0
sim_time = 0
##
function trigger_init()
    global sim_time = 0
    global last_trigger = 0
    fire!(net_tpn, state[], t_init_id)
    notify(state)
    fire_necessary!(net_tpn, state[])
    notify(state)
    println("Time: $(sim_time)")
end

function advance_trigger()
    next_trigger = last_trigger + trigger_period
    # execute all transitions until next trigger time
    while true
        τ = max_delay(net_tpn, state[])
        if sim_time + τ > next_trigger
            break
        end
        global sim_time += τ
        elapse!(net_tpn, state[], τ)
        fire_necessary!(net_tpn, state[])
        notify(state)
    end

    # advance to next trigger time
    elapse!(net_tpn, state[], next_trigger - sim_time)
    fire!(net_tpn, state[], t_trigger_id)
    notify(state)

    global sim_time = next_trigger
    global last_trigger = next_trigger
    # execute all transitions until next trigger time
    while true
        τ = max_delay(net_tpn, state[])
        if sim_time + τ > next_trigger
            break
        end
        global sim_time += τ
        elapse!(net_tpn, state[], τ)
        fire_necessary!(net_tpn, state[])
        notify(state)
    end

    println("Time: $(sim_time)")
end

##
trigger_init()
advance_trigger()
##
τ=max_delay(net_tpn,state[])
elapse!(net_tpn, state[], τ)
fire_necessary!(net_tpn, state[])
fire!(net_tpn, state[], t_trigger_id)
notify(state)


text:
module QuantumRelay

using SymPy
using Iterators
using PolyLib
using JuMP
using Clp
using LLLplus
using GSL
using PyPlot
using Lora
using Distributions

export
	qrelay_op,
	op_mat,
	scan_maker,
	QRelaySampler

include("util.jl")
include("operators.jl")
include("scan.jl")
include("distributions.jl")

type QRelaySampler
	prob::Function
	psetproposal::Function

	function QRelaySampler{T<:Int}(mat::Array{T, 2}, coef, omega, pdet0)
		U, S, V = PolyLib.smith_normal_form(mat)
		Ui = PolyLib.inverse(U)
		Vi = PolyLib.inverse(V)
		s = diag(S)
		r = countnz(s)
		s0 = s[1:r]
		@assert s0 == ones(r)
		ui1 = Ui[1:r, :]
		ui2 = Ui[r+1:end, :]
		vi1 = Vi[:, 1:r]
		vi2 = Vi[:, r+1:end]
		vi2 = lll(vi2)[1]
		T0 = vi1*ui1
		ui2oc = orthocomp(ui2)
		setc, scan = scan_maker(vi2)

		function prob(na)
		    @assert countnz(ui2*na) == 0
		    b = T0*na
		    setc(-b)
		    total = 0.0
		    for x in Task(scan)
		        nab = vi2*x + b
		        total += prod([c.^complex(n)/factorial(n) for (c, n) in zip(coef, nab)])
		    end
		    return abs(total*omega)^2
		end

		function prob(q, na, mask)
		    q0 = round(Int, q.>0)
		    m0 = round(Int, mask)
		    return prod((q0 + (1-2q0).*pdet0(na)).^m0)
		end

		psetproposal(x::Vector) = QuantumRelay.OrthoNNDist(x, ui2oc)

		new(prob, psetproposal)
	end
end

end


text:
type Elem <: XML
    name::AbstractString
    attrs::Dict{AbstractString, AbstractString}
    elems::Vector
    Elem() = new("", Dict(), [])
end

type Document <: XML
    doctype::Nullable{AbstractString}
    elems::Vector
    Document() = new(Nullable(), [])
end

function Document(l)
    doc = Document()
    stack = XML[doc]
    for e in l
        if isa(e, AbstractString)
            push!(stack[end].elems, e)
        elseif isa(e, Tag)
            t = Elem()
            t.name = e.name
            if e.t == 2
                #end tag
                t = pop!(stack)
                # maybe malformed documents (end tag that hasn't started)
                # should fail more grazefully than no field error
                while t.name != e.name
                    t = pop!(stack)
                end
                continue
            elseif e.t == 1
                #start tag
                push!(stack[end].elems, t)
                push!(stack, t)
            elseif e.t == 3
                # nothing
                push!(stack[end].elems, t)
            else
                unreachable
            end
            for i in e.attrs
                if length(i.value) ≥ 2 && first(i.value) == last(i.value) == '"'
                    t.attrs[i.name] = i.value[2:end-1]
                else
                    t.attrs[i.name] = i.value
                end
            end
        else
            unreachable
        end

    end
    doc
end

function Base.show(io::IO, d::Document)
    if !isnull(d.doctype)
        println(io, "<?", get(d.doctype), "?>")
    end
    for e in d.elems
        print(io, e)
    end
end

function Base.show(io::IO, d::Elem)
    print(io, "<", d.name)
    for (k,v) in d.attrs
        print(io, " ", k, "=\"", v, "\"")
    end
    print(io, ">")

    for e in d.elems
        print(io, e)
    end
    print(io, "</", d.name, ">")
end


text:
ModelType = MLJFlux.NeuralNetworkRegressor

@testset "equality" begin
    model = MLJFlux.ImageClassifier()
    clone = deepcopy(model)
    @test model == clone
    clone.optimiser.eta *= 10
    @test model != clone

    clone = deepcopy(model)
    clone.builder.dropout *= 0.5
    @test clone != model
end

@testset "clean!" begin
    model = @test_logs (:warn, r"`lambda") begin
        ModelType(lambda = -1)
    end
    @test model.lambda == 0

    model = @test_logs (:warn, r"`alpha") begin
        ModelType(alpha = -1)
    end
    @test model.alpha == 0

    model = @test_logs (:warn, r"`epochs") begin
        ModelType(epochs = -1)
    end
    @test model.epochs == 10

    model = @test_logs (:warn, r"`batch_size") begin
        ModelType(batch_size = -1)
    end
    @test model.batch_size == 1

    if MLJFlux.gpu_isdead()
        model  = @test_logs (:warn, r"`acceleration") begin
            ModelType(acceleration=CUDALibs())
        end
        @test model.acceleration == CUDALibs()
    end
end

@testset "iteration api" begin
    model = MLJFlux.NeuralNetworkRegressor(epochs=10)
    @test MLJBase.supports_training_losses(model)
    @test MLJBase.iteration_parameter(model) == fieldnames(typeof(model))[4]

    # integration test:
    X, y = MLJBase.make_regression(10)
    mach = MLJBase.machine(model, X, y)
    MLJBase.fit!(mach, verbosity=0)
    losses = MLJBase.training_losses(mach)
    @test losses == MLJBase.report(mach).training_losses[2:end]
    @test length(losses) == 10
end