text:
# Creation of packed types should be automated with macros

mutable struct PackedVariableNodeData
  vecinitval::Array{Float64,1}
  diminitval::Int
  vecinitstdev::Array{Float64,1}
  diminitdev::Int
  vecval::Array{Float64,1}
  dimval::Int
  vecbw::Array{Float64,1}
  dimbw::Int
  BayesNetOutVertIDs::Array{Int,1}
  dimIDs::Array{Int,1}
  dims::Int
  eliminated::Bool
  BayesNetVertID::Int
  separator::Array{Int,1}
  # groundtruth::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } }
  softtype::String
  initialized::Bool
  PackedVariableNodeData() = new()
  PackedVariableNodeData(x1::Vector{Float64},
                         x2::Int,
                         x3::Vector{Float64},
                         x4::Int,
                         x5::Vector{Float64},
                         x6::Int,
                         x7::Vector{Float64},
                         x8::Int,
                         x9::Vector{Int},
                         x10::Vector{Int},
                         x11::Int,
                         x12::Bool,
                         x13::Int,
                         x14::Vector{Int},
                         x15::String,
                         x16::Bool) = new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16)
end

# Obsolete
# FunctionNodeData() = GenericFunctionNodeData{T, Symbol}()
# typealias PackedFunctionNodeData{T <: PackedInferenceType} GenericFunctionNodeData{T, AbstractString}
# PackedFunctionNodeData{T}() where T = GenericFunctionNodeData{T, AbstractString}()
# PackedFunctionNodeData(x1, x2, x3, x4, x5::S, x6::T) where {T <: PackedInferenceType, S <: AbstractString} = GenericFunctionNodeData{T, AbstractString}(x1, x2, x3, x4, x5, x6)

# TODO stop-gap string storage of Distrubtion types, should be upgraded to more efficient storage
function normalfromstring(str::AS) where {AS <: AbstractString}
  meanstr = match(r"μ=[+-]?([0-9]*[.])?[0-9]+", str).match
  mean = split(meanstr, '=')[2]
  sigmastr = match(r"σ=[+-]?([0-9]*[.])?[0-9]+", str).match
  sigma = split(sigmastr, '=')[2]
  Normal{Float64}(parse(Float64,mean), parse(Float64,sigma))
end

function mvnormalfromstring(str::AS) where {AS <: AbstractString}
  means = split(split(split(str, 'μ')[2],']')[1],'[')[end]
  mean = Float64[]
  for ms in split(means, ',')
    push!(mean, parse(Float64, ms))
  end
  sigs = split(split(split(str, 'Σ')[2],']')[1],'[')[end]
  sig = Float64[]
  for ms in split(sigs, ';')
    for m in split(ms, ' ')
      length(m) > 0 ? push!(sig, parse(Float64, m)) : nothing
    end
  end
  len = length(mean)
  sigm = reshape(sig, len,len)
  MvNormal(mean, sigm)
end

function categoricalfromstring(str::AS)::Distributions.Categorical where {AS <: AbstractString}
  # pstr = match(r"p=\[", str).match
  psubs = split(str, '=')[end]
  psubs = split(psubs, '[')[end]
  psubsub = split(psubs, ']')[1]
  pw = split(psubsub, ',')
  p = parse.(Float64, pw)
  return Categorical(p ./ sum(p))
end

function extractdistribution(str::AS)::Union{Void, Distributions.Distribution} where {AS <: AbstractString}
  # TODO improve use of multidispatch and packing of Distribution types
  if str == ""
    return nothing
  elseif (ismatch(r"Normal", str) && !ismatch(r"FullNormal", str))
    return normalfromstring(str)
  elseif ismatch(r"FullNormal", str)
    return mvnormalfromstring(str)
  elseif ismatch(r"Categorical", str)
    return categoricalfromstring(str)
  elseif ismatch(r"KDE:", str)
    return convert(KernelDensityEstimate.BallTreeDensity, str)
  else
    error("Don't know how to extract distribution from str=$(str)")
  end
end

function convert(::Type{PackedVariableNodeData}, d::VariableNodeData)
  return PackedVariableNodeData(d.initval[:],size(d.initval,1),
                              d.initstdev[:],size(d.initstdev,1),
                              d.val[:],size(d.val,1),
                              d.bw[:], size(d.bw,1),
                              d.BayesNetOutVertIDs,
                              d.dimIDs, d.dims, d.eliminated,
                              d.BayesNetVertID, d.separator,
                              string(d.softtype), d.initialized)
end
function convert(::Type{VariableNodeData}, d::PackedVariableNodeData)

  r1 = d.diminitval
  c1 = r1 > 0 ? floor(Int,length(d.vecinitval)/r1) : 0
  M1 = reshape(d.vecinitval,r1,c1)

  r2 = d.diminitdev
  c2 = r2 > 0 ? floor(Int,length(d.vecinitstdev)/r2) : 0
  M2 = reshape(d.vecinitstdev,r2,c2)

  r3 = d.dimval
  c3 = r3 > 0 ? floor(Int,length(d.vecval)/r3) : 0
  M3 = reshape(d.vecval,r3,c3)

  r4 = d.dimbw
  c4 = r4 > 0 ? floor(Int,length(d.vecbw)/r4) : 0
  M4 = reshape(d.vecbw,r4,c4)

  # TODO -- allow out of module type allocation (future feature, not currently in use)
  st = IncrementalInference.ContinuousMultivariate # eval(parse(d.softtype))

  return VariableNodeData(M1,M2,M3,M4, d.BayesNetOutVertIDs,
    d.dimIDs, d.dims, d.eliminated, d.BayesNetVertID, d.separator,
    nothing, st, d.initialized )
end

function compare(a::VariableNodeData,b::VariableNodeData)
    TP = true
    TP = TP && a.initval == b.initval
    TP = TP && a.initstdev == b.initstdev
    TP = TP && a.val == b.val
    TP = TP && a.bw == b.bw
    TP = TP && a.BayesNetOutVertIDs == b.BayesNetOutVertIDs
    TP = TP && a.dimIDs == b.dimIDs
    TP = TP && a.dims == b.dims
    TP = TP && a.eliminated == b.eliminated
    TP = TP && a.BayesNetVertID == b.BayesNetVertID
    TP = TP && a.separator == b.separator
    return TP
end

function ==(a::VariableNodeData,b::VariableNodeData, nt::Symbol=:var)
  return IncrementalInference.compare(a,b)
end

function packmultihypo(fnc::CommonConvWrapper{T}) where {T<:FunctorInferenceType}
  fnc.hypotheses != nothing ? string(fnc.hypotheses) : ""
end
function parsemultihypostr(str::AS) where {AS <: AbstractString}
  mhcat=nothing
  if length(str) > 0
    mhcat = extractdistribution(str)
  end
  return mhcat
end

# heavy use of multiple dispatch for converting between packed and original data types during DB usage

# TODO simplify and reduce to single pack and unpack converter
## packing converters-----------------------------------------------------------

# function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: InferenceType}
#     # println("convert(::Type{PackedFunctionNodeData{$P}}, d::FunctionNodeData{$T})")
#   error("convert(::Type{PackedFunctionNodeData{P}} : this convert function should probably not be used. Use FunctorInferenceType instead of InferenceType.")
#   mhstr = packmultihypo(d.fnc)
#   return PackedFunctionNodeData{P}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
#           string(d.frommodule), convert(P, d.fnc), mhstr) # TODO -- should use convert(P, d.fnc.usrfnc!) instead
# end
function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: FunctorInferenceType}
  # println("convert(::Type{PackedFunctionNodeData{$P}}, d::FunctionNodeData{$T})")
  mhstr = packmultihypo(d.fnc)
  return PackedFunctionNodeData(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc.usrfnc!), mhstr)
end
function convert(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T}) where {P <: PackedInferenceType, T <: ConvolutionObject}
  # println("convert(::Type{PackedFunctionNodeData{$P}}, d::FunctionNodeData{$T})")
  mhstr = packmultihypo(d.fnc)
  return PackedFunctionNodeData(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc.usrfnc!), mhstr)
end

## unpack converters------------------------------------------------------------

function convert(
            ::Type{IncrementalInference.GenericFunctionNodeData{IncrementalInference.CommonConvWrapper{F},Symbol}},
            d::IncrementalInference.GenericFunctionNodeData{P,String} ) where {F <: FunctorInferenceType, P <: PackedInferenceType}
  #
  # warn("Unpacking Option 2, F=$(F), P=$(P)")
  usrfnc = convert(F, d.fnc)
  # @show d.multihypo
  mhcat = parsemultihypostr(d.multihypo)

  # TODO store threadmodel=MutliThreaded,SingleThreaded in persistence layer
  ccw = prepgenericconvolution(Graphs.ExVertex[], usrfnc, multihypo=mhcat)
  return FunctionNodeData{CommonConvWrapper{typeof(usrfnc)}}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          Symbol(d.frommodule), ccw)
end

# Compare FunctionNodeData
function compare(a::GenericFunctionNodeData{T1,S},b::GenericFunctionNodeData{T2,S}) where {T1, T2, S}
  # TODO -- beef up this comparison to include the gwp
  TP = true
  TP = TP && a.fncargvID == b.fncargvID
  TP = TP && a.eliminated == b.eliminated
  TP = TP && a.potentialused == b.potentialused
  TP = TP && a.edgeIDs == b.edgeIDs
  TP = TP && a.frommodule == b.frommodule
  # TP = TP && typeof(a.fnc) == typeof(b.fnc)
  return TP
end

function convert{PT <: PackedInferenceType, T <:FunctorInferenceType}(::Type{PT}, ::T)
  getfield(T.name.module, Symbol("Packed$(T.name.name)"))
end
function convert{T <: FunctorInferenceType, PT <: PackedInferenceType}(::Type{T}, ::PT)
  getfield(PT.name.module, Symbol(string(PT.name.name)[7:end]))
end

function getmodule(t::T) where T
  T.name.module
end
function getname(t::T) where T
  T.name.name
end
function encodePackedType(topackdata::VariableNodeData)
  # error("IncrementalInference.encodePackedType(::VariableNodeData): Unknown packed type encoding of $(topackdata)")
  convert(IncrementalInference.PackedVariableNodeData, topackdata)
end
function encodePackedType(topackdata::GenericFunctionNodeData{T, Symbol}) where {T <: FunctorInferenceType}
  # println("IncrementalInference.encodePackedType(::GenericFunctionNodeData{T,Symbol}): Unknown packed type encoding T=$(T) of $(topackdata)")
  fnctype = getfnctype(topackdata)
  fnc = getfield(getmodule(fnctype), Symbol("Packed$(getname(fnctype))"))
  convert(PackedFunctionNodeData{fnc}, topackdata)
end
function encodePackedType(topackdata::GenericFunctionNodeData{T, <:AbstractString}) where {T <: PackedInferenceType}
  error("IncrementalInference.encodePackedType(::FunctionNodeData{T, <:AbstractString}): Unknown packed type encoding T=$(T) of $(topackdata)")
  # @show T, typeof(topackdata)
  # warn("Yes, its packed!")
  # fnctype = getfnctype(topackdata)
  # @show fnc = getfield(getmodule(fnctype), Symbol("Packed$(getname(fnctype))"))
  # convert(PackedFunctionNodeData{T}, topackdata)
  topackdata
end

"""
    $(SIGNATURES)

Encode complicated function node type to related 'Packed<type>' format assuming a user supplied convert function .
"""
function convert2packedfunctionnode(fgl::FactorGraph,
      fsym::Symbol,
      api::DataLayerAPI=localapi  )
  #
  fid = fgl.fIDs[fsym]
  fnc = getfnctype(fgl, fid)
  usrtyp = convert(PackedInferenceType, fnc)
  cfnd = convert(PackedFunctionNodeData{usrtyp},getData(fgl, fid, api=api) )
  return cfnd, usrtyp
end

function decodePackedType(packeddata::PackedVariableNodeData, typestring::String)
  # error("IncrementalInference.encodePackedType(::VariableNodeData): Unknown packed type encoding of $(topackdata)")
  convert(IncrementalInference.VariableNodeData, packeddata)
end
function decodePackedType(packeddata::GenericFunctionNodeData{PT,<:AbstractString}, typestring::String) where {PT}
  # warn("decodePackedType($(typeof(packeddata)),$(typestring)) is happening with PT=$(PT) and ")
  functype = getfield(PT.name.module, Symbol(string(PT.name.name)[7:end]))
  fulltype = FunctionNodeData{CommonConvWrapper{functype}}
  convert(fulltype, packeddata)
end

"""
    $(SIGNATURES)

Make a full memory copy of the graph and encode all complicated function node
types with assumed to exist convert to 'Packed<type>' formats. Same converters
as used for database persistence storage with CloudGraphs.jl.
"""
function encodefg(fgl::FactorGraph;
      api::DataLayerAPI=localapi  )
  #
  fgs = deepcopy(fgl)
  fgs.cg = nothing
  fgs.registeredModuleFunctions = nothing
  fgs.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
  @showprogress 1 "Encoding variables..." for (vsym,vid) in fgl.IDs
    cpvert = deepcopy(getVert(fgl, vid, api=api))
    api.addvertex!(fgs, cpvert) #, labels=vnlbls)  # currently losing labels
  end

  @showprogress 1 "Encoding factors..." for (fsym,fid) in fgs.fIDs
    data,ftyp = convert2packedfunctionnode(fgl, fsym)
    # data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    newvert = ExVertex(fid,string(fsym))
    for (key,val) in getVert(fgl,fid,api=api).attributes
      newvert.attributes[key] = val
    end
    setData!(newvert, data)
    api.addvertex!(fgs, newvert)
  end
  fgs.g.inclist = typeof(fgl.g.inclist)()

  # iterated over all edges
  @showprogress 1 "Encoding edges..." for (eid, edges) in fgl.g.inclist
    fgs.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgs.g.vertices[ed.source.index],
          fgs.g.vertices[ed.target.index]  )
      push!(fgs.g.inclist[eid], newed)
    end
  end

  return fgs
end

"""
    $(SIGNATURES)

If you get unknown type conversion error when loading a .jld, while using your own
FunctorInferenceTypes, you should:
Copy these functions below, and overload in your package with local extented
FunctorInferenceType definitions.
See RoME/src/fgos.jl for example.
"""
function convertfrompackedfunctionnode(fgl::FactorGraph,
      fsym::Symbol,
      api::DataLayerAPI=localapi  )
  #
  fid = fgl.fIDs[fsym]
  fnc = getData(fgl, fid).fnc #getfnctype(fgl, fid)
  usrtyp = convert(FunctorInferenceType, fnc)
  data = getData(fgl, fid, api=api)
  newtype = FunctionNodeData{CommonConvWrapper{usrtyp}}
  cfnd = convert(newtype, data)
  return cfnd, usrtyp
end

"""
    decodefg(fgs::FactorGraph)

Unpack PackedFunctionNodeData formats back to regular FunctonNodeData.
"""
function decodefg(fgs::FactorGraph; api::DataLayerAPI=localapi)
  fgu = deepcopy(fgs)
  fgu.cg = nothing
  fgu.registeredModuleFunctions = nothing
  fgu.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
  @showprogress 1 "Decoding variables..." for (vsym,vid) in fgs.IDs
    cpvert = deepcopy(getVert(fgs, vid, api=api))
    api.addvertex!(fgu, cpvert) #, labels=vnlbls)  # currently losing labels
  end

  @showprogress 1 "Decoding factors..." for (fsym,fid) in fgu.fIDs
    data,ftyp = convertfrompackedfunctionnode(fgs, fsym)
    # data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    newvert = ExVertex(fid,string(fsym))
    for (key,val) in getVert(fgs,fid,api=api).attributes
      newvert.attributes[key] = val
    end
    setData!(newvert, data)
    api.addvertex!(fgu, newvert)
  end
  fgu.g.inclist = typeof(fgs.g.inclist)()

  # iterated over all edges
  @showprogress 1 "Decoding edges..." for (eid, edges) in fgs.g.inclist
    fgu.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgu.g.vertices[ed.source.index],
          fgu.g.vertices[ed.target.index]  )
      push!(fgu.g.inclist[eid], newed)
    end
  end

  return fgu
end

#


text:
using Pkg
using Statistics

IRIS_DOWNLOAD_PATH = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"

@testset "Create job using remote BanyanDataFrames Github with $pf_dispatch_table pf_dispatch_table" for pf_dispatch_table in
    ["default", "http path"]

    Pkg.activate("./")

    if pf_dispatch_table == "default"
        pf_dispatch_table = ""
    elseif pf_dispatch_table == "http path"
        pf_dispatch_table = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/Banyan/res/pt_lib_info.json"
    end
    files = [IRIS_DOWNLOAD_PATH]

    # Create job
    job_id = create_job(
        cluster_name = ENV["BANYAN_CLUSTER_NAME"],
        nworkers = 2,
        files = files,
        pf_dispatch_table = pf_dispatch_table,
        url = "https://github.com/banyan-team/banyan-julia.git",
        branch = "remove_banyanfile",  # TODO: Change to "v0.1.3",
        directory = "banyan-julia/BanyanDataFrames/test",
        dev_paths = [
            "banyan-julia/Banyan",
            "banyan-julia/BanyanDataFrames",
        ],
        force_clone = true
    )
    @test get_job_id() == job_id
    @test get_job().cluster_name == ENV["BANYAN_CLUSTER_NAME"]
    @test get_job().nworkers == 2
    @test get_cluster_name() == ENV["BANYAN_CLUSTER_NAME"]

    # Describe jobs
    curr_jobs = get_jobs(ENV["BANYAN_CLUSTER_NAME"], status="running")
    @test haskey(curr_jobs, job_id)
    @test curr_jobs[job_id]["status"] == "running"

    # Perform computation 1
    bucket_name = get_cluster_s3_bucket_name(ENV["BANYAN_CLUSTER_NAME"])
    df = read_csv("s3://$(bucket_name)/iris.csv")
    gdf = groupby(df, :species)
    pl_means = combine(gdf, :petal_length => mean)
    res = collect(pl_means)
    @test res[:, :petal_length_mean] == [1.464, 4.26, 5.552]
    @test res[:, :species] == ["setosa", "versicolor", "virginica"]

    # Destroy job
    destroy_job(job_id)
    curr_jobs = get_jobs(ENV["BANYAN_CLUSTER_NAME"], status="running")
    @test !haskey(curr_jobs, job_id)
end

@testset "Create job using local Julia environment" begin
    # Activate an environment
    Pkg.activate("envs/local")

    # Import some packages
    using Banyan
    using BanyanArrays
    using Statistics
    using Distributions

    # Create a job
    job_id = create_job(
        cluster_name = ENV["BANYAN_CLUSTER_NAME"],
        nworkers = 2,
    )
    # Perform some computation
    size = 64_000
    x = BanyanArrays.fill(0.5, size)
    y = BanyanArrays.map(
        t -> cdf(Normal(), t),
        x
    )
    res = collect(sum(y))
    @test round(res) == 44254.0

    # Destroy job
    destroy_job(job_id)
end

# @testset "Create job using remote Github environment"
    # Pkg.activate("./")
# end

ext: jl

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
abstract type Adjoint_rep_Gaugefields{NC,Dim,NumofBasis}  <: AbstractGaugefields{NC,Dim} #Traceless antihermitian matrix
end

include("./4D/Adjoint_rep_gaugefields_4D.jl")

function construct_Adjoint_rep_Gaugefields(U::Array{<:AbstractGaugefields{NC,Dim},1}) where {NC,Dim}
    F = Array{Adjoint_rep_Gaugefields{NC,Dim},1}(undef,Dim)
    for μ=1:Dim
        F[μ] = construct_Adjoint_rep_Gaugefields(U[μ])
    end
    return F
end

function construct_Adjoint_rep_Gaugefields(u::AbstractGaugefields{NC,Dim};verbose_level=2) where {NC,Dim}
    mpi = u.mpi
    if mpi
        error("MPI is not supported")
        if Dim==4
            return construct_Adjoint_rep_Gaugefields_mpi(u)
        elseif Dim==2
            error("Dim = $Dim is not supoorted")
        else
            error("Dim = $Dim is not supoorted")
        end

    else
        if Dim==4
            return construct_Adjoint_rep_Gaugefields_4D_wing(u,verbose_level=u.verbose_print.level)
        elseif Dim==2
            error("Dim = $Dim is not supoorted")
        else
            error("Dim = $Dim is not supoorted")
        end

    end
end


text:
function read_net_file(filename::String)
  M = readdlm(filename)
  ei = M[:,1]
  ej = M[:,2]

  ei = map(x->x[2:end],ei)
  ej = map(x->x[2:end],ej)

  ei = map(x->parse(Int64,x),ei)
  ej = map(x->parse(Int64,x),ej)

  m = maximum(ei)
  n = maximum(ej)
  d = max(m,n)

  return sparse(ei,ej,1,d,d)
end

function read_sim_file(filename::String,m::Int64,n::Int64)
  M = readdlm(filename)
  ei = M[:,1]
  ej = M[:,2]

  ei = map(x->x[2:end],ei)
  ej = map(x->x[2:end],ej)

  ei = map(x->parse(Int64,x),ei)
  ej = map(x->parse(Int64,x),ej)

  w = M[:,3]
  w = convert(Array{Float64},w)

  return sparse(ei,ej,w,m,n)
end

function setup_yeast_networks(file1,file2)
    lb1,A1 = convert_LEDA_to_adjacency(file1)
    lb2,A2 = convert_LEDA_to_adjacency(file2)
    n = size(A2,1)
    ei,ej,ev = findnz(A2)
    eii = n+1-ei # flip to unbias
    ejj = n+1-ej # flip to unbias
    B = sparse(eii,ejj,ev)
    Pref = sparse(1:n,collect(n:-1:1),1,size(A1,1),size(B,1))
    return A1,B,Pref
end

function gs3_evaluate(A,B,mi,mj)
  mi = collect(mi)
  mj = collect(mj)
  num_matches = length(mi)
  NCV = (length(unique(mi)) + length(unique(mj))) / (size(A, 1) + size(B, 1));

  #ii,jj,vv = findnz(A)
  #ei = mi[ii]
  #ej = mi[jj]
  #A_induced = sparse(ei,ej,vv)

  #ii,jj,vv = findnz(B)
  #ei = mj[ii]
  #ej = mj[jj]
  #B_induced = sparse(ei,ej,vv)

  A_induced = A[mi, mi];
  B_induced = B[mj, mj];
  D = A_induced + B_induced;
  D = D - spdiagm(diag(D))
  ii,jj = ind2sub(D,find(D.==2))
  align_graph = sparse(ii, jj, 1, num_matches, num_matches)

  M_align = nnz(align_graph) / 2
  # T_align = full(trace(align_graph^3)/6);

  GS3 = M_align / (nnz(A_induced)/2 + nnz(B_induced)/2 - M_align)
  NCV_GS3 = sqrt(GS3*NCV)
  return NCV_GS3
end

function gs3_evaluate2(A,B,mi,mj)
  mi = collect(mi)
  mj = collect(mj)
  num_matches = length(mi)
  NCV = (length(unique(mi)) + length(unique(mj))) / (size(A, 1) + size(B, 1));

  #ii,jj,vv = findnz(A)
  #ei = mi[ii]
  #ej = mi[jj]
  #A_induced = sparse(ei,ej,vv)

  #ii,jj,vv = findnz(B)
  #ei = mj[ii]
  #ej = mj[jj]
  #B_induced = sparse(ei,ej,vv)

  A_induced = A[mi, mi];
  B_induced = B[mj, mj];
  D = A_induced + B_induced;
  D = D - spdiagm(diag(D))
  ii,jj = ind2sub(D,find(D.==2))
  align_graph = sparse(ii, jj, 1, num_matches, num_matches)

  M_align = nnz(align_graph) / 2
  # T_align = full(trace(align_graph^3)/6);

  GS3 = M_align / (nnz(A_induced)/2 + nnz(B_induced)/2 - M_align)
  NCV_GS3 = 1/(((1/GS3)+(1/NCV))/2)
  return NCV_GS3
end

"""
sortcolsperm(A,true) returns the indices of sorted columns in A in descending order
sortcolsperm(A,false) returns the indices of sorted columns in A in ascending order
```
example:
  julia> W = rand(3,3)
    3×3 Array{Float64,2}:
    0.661943  0.00517749  0.332394
    0.716344  0.61179     0.544258
    0.372336  0.994069    0.297704
  julia> sortcolsperm(W,true)
    3×3 Array{Int64,2}:
    2  3  2
    1  2  1
    3  1  3
  julia> sortcolsperm(W,false)
    3×3 Array{Int64,2}:
    3  1  3
    1  2  1
    2  3  2
```
"""
function sortcolsperm(X::Matrix{T},REV::Bool) where T
    P = Matrix{Int}(undef,size(X,1),size(X,2))

    Threads.@threads for i=1:size(X,2)
        P[:,i] = sortperm(X[:,i]; rev=REV)
    end
    return P
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
function process_simulation_files(path_to_sim_files)

  # Load the data from disk -
  file_extension = ".dat"
  searchdir(path,key) = filter(x->contains(x,key),readdir(path))

  # build src file list -
  list_of_sim_files = searchdir(path_to_sim_files,file_extension);

  # initialize set -
  simulation_data_result_set::Set{Array} = Set{Array}()

  # go thru the src file list, and copy the files to the output path -
  for (sim_file_index,sim_file) in enumerate(list_of_sim_files)

    # path to sim file -
    sim_file_path = "./"*path_to_sim_files*"/"*string(sim_file)

    # load the data array -
    simulation_data_array = readdlm(sim_file_path);

    # add to set -
    push!(simulation_data_result_set,simulation_data_array);
  end

  return simulation_data_result_set
end

function process_ensemble_data(path_to_sim_files,time_array,species_index)

  # Load the data from disk -
  file_extension = ".dat"
  searchdir(path,key) = filter(x->contains(x,key),readdir(path))

  # build src file list -
  list_of_sim_files = searchdir(path_to_sim_files,file_extension);

  # initialize data array -
  raw_data_array = zeros(length(time_array),1);

  # go thru the src file list, and copy the files to the output path -
  for (sim_file_index,sim_file) in enumerate(list_of_sim_files)

    # path to sim file -
    sim_file_path = "./"*path_to_sim_files*"/"*string(sim_file)

    # load the data array -
    simulation_data_array = readdlm(sim_file_path);

    # Get the time -
    simulation_time_array = simulation_data_array[:,1];
    simulation_state_array = simulation_data_array[:,species_index+1];

    # interpolate -
    interpolated_data_array = np.interp(time_array,simulation_time_array,simulation_state_array);

    @show (sim_file_index,sim_file)
    raw_data_array = hcat(raw_data_array,transpose(interpolated_data_array))
  end

  # Generate the scaled data -
  raw_data_array = raw_data_array[:,2:end]
  (number_of_rows,number_of_cols) = size(raw_data_array);
  scaled_data_array = zeros(number_of_rows,number_of_cols);
  for col_index = 1:number_of_cols

    # grab the first elements -
    scale_element = raw_data_array[1,col_index];
    for row_index = 1:number_of_rows
      scaled_value_element = raw_data_array[row_index,col_index]*(1/scale_element);
      scaled_data_array[row_index,col_index] = scaled_value_element;
    end
  end

  return (raw_data_array,scaled_data_array)
end


text:
"""
Retrieve intraday bar data (size 30 minutes by default) as a time series.
The start date is required, the end date is optional. These arguments may be
dates of the form YYYY-MM-DD or times of the form YYYY-MM-DDTHH:MM:SS.
The event type can be one of TRADE, BID, ASK, BEST_BID, BEST_ASK.

#Examples
```julia
bbars("AAPL US Equity","TRADE","2017-09-27",end_date_time="2017-09-29")
bbars("AAPL US Equity","BID","2017-09-27T06:00:00",end_date_time="2017-09-27T09:00:00",bar_size_minutes=10)
```
"""
function bbars(ticker::String,event::String,start_date_time::String;end_date_time::String="",bar_size_minutes::Int=30,overrides=Dict())
    global session
    global services
    @assert event in ["TRADE", "BID", "ASK", "BEST_BID", "BEST_ASK"]
    service = services["refdata"]
    request = @with_pointer request blpapi_Service_createRequest(service,request,"IntradayBarRequest")
    elements = blpapi_Request_elements(request)
    try
        blpapi_Element_setElementString(elements,"security",C_NULL,ticker)
        blpapi_Element_setElementString(elements,"eventType",C_NULL,event)
        blpapi_Element_setElementInt32(elements,"interval",C_NULL,bar_size_minutes)
        @assert length(start_date_time) >= 10
        blpapi_Element_setElementString(elements,"startDateTime",C_NULL,start_date_time)
        if length(end_date_time) > 0
            @assert length(end_date_time) >= 10
            blpapi_Element_setElementString(elements,"endDateTime",C_NULL,end_date_time)
        end
        if length(overrides) > 0
            overrides_ = @with_pointer overrides_ blpapi_Element_getElement(elements,overrides_,"overrides",C_NULL)
            for (k,v) in overrides
                new_ = @with_pointer new_ blpapi_Element_appendElement(overrides_,new_)
                blpapi_Element_setElementString(new_, "fieldId", C_NULL, string(k))
                blpapi_Element_setElementString(new_, "value", C_NULL, string(v))
            end
        end
        correlation_id = blpapi_CorrelationId(0,0)
        blpapi_Session_sendRequest(session,request,Ref(correlation_id),C_NULL,C_NULL,C_NULL,0)
    finally
        blpapi_Request_destroy(request)
    end
    out = process_tree(retrieve_response())
    if size(out,1) > 0
        TimeSeries.TimeArray(out, timestamp=:time)
    end
end


text:
"""
    side_indices(c::GridContent)::RowCols{Int}

Indices of the rows / cols for each side
"""
function side_indices(c::GridContent)
    return RowCols(
        c.sp.cols.start,
        c.sp.cols.stop,
        c.sp.rows.start,
        c.sp.rows.stop,
    )
end

"""
These functions tell whether an object in a grid touches the left, top, etc. border
of the grid. This means that it is relevant for the grid's own protrusion on that side.
"""
ismostin(sp::GridContent, grid, ::Left) = sp.sp.cols.start == 1
ismostin(sp::GridContent, grid, ::Right) = sp.sp.cols.stop == grid.ncols
ismostin(sp::GridContent, grid, ::Bottom) = sp.sp.rows.stop == grid.nrows
ismostin(sp::GridContent, grid, ::Top) = sp.sp.rows.start == 1

isleftmostin(sp::GridContent, grid) = ismostin(sp, grid, Left())
isrightmostin(sp::GridContent, grid) = ismostin(sp, grid, Right())
isbottommostin(sp::GridContent, grid) = ismostin(sp, grid, Bottom())
istopmostin(sp::GridContent, grid) = ismostin(sp, grid, Top())

# parentlayout(pl::ProtrusionLayout) = pl.parent
#
# function ProtrusionLayout(content)
#     needs_update = Node(false)
#
#     protrusions = protrusionnode(content)
#     csize = computedsizenode(content)
#
#     # the nothing parent here has to be replaced by a grid parent later
#     # when placing this layout in the grid layout
#     pl = ProtrusionLayout(nothing, protrusions, csize, needs_update, content)
#
#     update_func = x -> begin
#         p = parentlayout(pl)
#         if isnothing(p)
#             error("This layout has no parent and can't continue the update signal chain.")
#         end
#         p.needs_update[] = true
#     end
#
#     on(update_func, protrusions)
#     on(update_func, csize)
#
#     pl
# end

# protrusionnode(pl::ProtrusionLayout) = pl.protrusions
# computedsizenode(pl::ProtrusionLayout) = pl.computedsize

# """
# All the protrusion functions calculate how much stuff "sticks out" of a layoutable object.
# This is so that collisions are avoided, while what is actually aligned is the
# "important" edges of the layout objects.
# """
# leftprotrusion(x) = protrusion(x, Left())
# rightprotrusion(x) = protrusion(x, Right())
# bottomprotrusion(x) = protrusion(x, Bottom())
# topprotrusion(x) = protrusion(x, Top())

# protrusion(a::ProtrusionLayout, ::Left) = isnothing(a.protrusions[]) ? 0f0 : a.protrusions[].left
# protrusion(a::ProtrusionLayout, ::Right) = isnothing(a.protrusions[]) ? 0f0 : a.protrusions[].right
# protrusion(a::ProtrusionLayout, ::Top) = isnothing(a.protrusions[]) ? 0f0 : a.protrusions[].top
# protrusion(a::ProtrusionLayout, ::Bottom) = isnothing(a.protrusions[]) ? 0f0 : a.protrusions[].bottom

function protrusionnode(x::T) where T
    error("protrusionnode() is not defined for type $T, if you want to include
        such an object in a grid layout, you need to create a method for your
        type that returns a Node{RectSides{Float32}}. If your object has no
        protrusions, the values of RectSides can just all be 0f0.")
end

function protrusion(x::T, side::Side) where T
    protrusions = protrusionnode(x)
    @match side begin
        si::Left => protrusions[].left
        si::Right => protrusions[].right
        si::Bottom => protrusions[].bottom
        si::Top => protrusions[].top
        si => error("Can't get a protrusion value for side $(typeof(si)), only
            Left, Right, Bottom, or Top.")
    end
end

function protrusion(sp::GridContent, side::Side)
    prot = @match sp.side begin
        sps::Inner => protrusion(sp.al, side)
        # sps::Outer => BBox(l - pl, r + pr, b - pb, t + pt)
        sps::Union{Left, Right} => @match side begin
            si::typeof(sps) => determinedirsize(sp.al, Col(), sp.side)
            si => 0.0
        end
        sps::Union{Top, Bottom} => @match side begin
            si::typeof(sps) => determinedirsize(sp.al, Row(), sp.side)
            si => 0.0
        end
        sps::TopLeft => @match side begin
            si::Top => determinedirsize(sp.al, Row(), sp.side)
            si::Left => determinedirsize(sp.al, Col(), sp.side)
            si => 0.0
        end
        sps::TopRight => @match side begin
            si::Top => determinedirsize(sp.al, Row(), sp.side)
            si::Right => determinedirsize(sp.al, Col(), sp.side)
            si => 0.0
        end
        sps::BottomLeft => @match side begin
            si::Bottom => determinedirsize(sp.al, Row(), sp.side)
            si::Left => determinedirsize(sp.al, Col(), sp.side)
            si => 0.0
        end
        sps::BottomRight => @match side begin
            si::Bottom => determinedirsize(sp.al, Row(), sp.side)
            si::Right => determinedirsize(sp.al, Col(), sp.side)
            si => 0.0
        end
        sps => error("Invalid side $sps")
    end
    ifnothing(prot, 0.0)
end

function protrusion(gl::GridLayout, side::Side)
    # when we align with the outside there is by definition no protrusion
    if gl.alignmode isa Outside
        return 0.0
    elseif gl.alignmode isa Inside
        prot = 0.0
        for elem in gl.content
            if ismostin(elem, gl, side)
                # take the max protrusion of all elements that are sticking
                # out at this side
                prot = max(protrusion(elem, side), prot)
            end
        end
        return prot
    end
end

# protrusion(s::SolvedProtrusionLayout, ::Left) = left(s.inner) - left(s.outer)
# function protrusion(s::SolvedProtrusionLayout, side::Side)
#     return s.outer[side] - s.inner[side]
# end

function bbox_for_solving_from_side(maxgrid::RowCols, bbox_cell::BBox, idx_rect::RowCols, side::Side)
    pl = maxgrid.lefts[idx_rect.lefts]
    pr = maxgrid.rights[idx_rect.rights]
    pt = maxgrid.tops[idx_rect.tops]
    pb = maxgrid.bottoms[idx_rect.bottoms]

    l = left(bbox_cell)
    r = right(bbox_cell)
    b = bottom(bbox_cell)
    t = top(bbox_cell)

    @match side begin
        s::Inner => bbox_cell
        s::Outer => BBox(l - pl, r + pr, b - pb, t + pt)
        s::Left => BBox(l - pl, l, b, t)
        s::Top => BBox(l, r, t, t + pt)
        s::Right => BBox(r, r + pr, b, t)
        s::Bottom => BBox(l, r, b - pb, b)
        s::TopLeft => BBox(l - pl, l, t, t + pt)
        s::TopRight => BBox(r, r + pr, t, t + pt)
        s::BottomRight => BBox(r, r + pr, b - pb, b)
        s::BottomLeft => BBox(l - pl, l, b - pb, b)
        s => error("Invalid side $s")
    end
end

startside(c::Col) = Left()
stopside(c::Col) = Right()
startside(r::Row) = Top()
stopside(r::Row) = Bottom()

span(sp::GridContent, dir::Col) = sp.sp.cols
span(sp::GridContent, dir::Row) = sp.sp.rows

"""
Determine the size of a protrusion layout along a dimension. This size is dependent
on the `Side` at which the layout is placed in its parent grid. An `Inside` side
means that the protrusion layout reports its width but not its protrusions. `Left`
means that the layout reports only its full width but not its height, because
an element placed in the left protrusion loses its ability to influence height.
"""
function determinedirsize(content, gdir::GridDir, side::Side)
    computedsize = computedsizenode(content)
    if gdir isa Row
        @match side begin
            # TODO: is computedsize the correct thing to return? or plus protrusions depending on the side
            si::Union{Inner, Top, Bottom, TopLeft, TopRight, BottomLeft, BottomRight} =>
                ifnothing(computedsize[][2], nothing)
            si::Union{Left, Right} => nothing
            si => error("$side not implemented")
        end
    else
        @match side begin
            si::Union{Inner, Left, Right, TopLeft, TopRight, BottomLeft, BottomRight} =>
                ifnothing(computedsize[][1], nothing)
            si::Union{Top, Bottom} => nothing
            si => error("$side not implemented")
        end
    end
end

# function solve(ua::ProtrusionLayout, innerbbox)
#     ib = innerbbox
#     bbox = BBox(
#         left(ib) - protrusion(ua, Left()),
#         right(ib) + protrusion(ua, Right()),
#         bottom(ib) - protrusion(ua, Bottom()),
#         top(ib) + protrusion(ua, Top()))
#     SolvedProtrusionLayout(innerbbox, ua.content)
# end

function to_ranges(g::GridLayout, rows::Indexables, cols::Indexables)
    if rows isa Int
        rows = rows:rows
    elseif rows isa Colon
        rows = 1:g.nrows
    end
    if cols isa Int
        cols = cols:cols
    elseif cols isa Colon
        cols = 1:g.ncols
    end
    rows, cols
end

function adjust_rows_cols!(g::GridLayout, rows, cols)
    rows, cols = to_ranges(g, rows, cols)

    if rows.start < 1
        n = 1 - rows.start
        prependrows!(g, n)
        # adjust rows for the newly prepended ones
        rows = rows .+ n
    end
    if rows.stop > g.nrows
        n = rows.stop - g.nrows
        appendrows!(g, n)
    end
    if cols.start < 1
        n = 1 - cols.start
        prependcols!(g, n)
        # adjust cols for the newly prepended ones
        cols = cols .+ n
    end
    if cols.stop > g.ncols
        n = cols.stop - g.ncols
        appendcols!(g, n)
    end

    rows, cols
end


text:
# Autogenerated wrapper script for Bzip2_jll for x86_64-unknown-freebsd
export libbzip2

JLLWrappers.@generate_wrapper_header("Bzip2")
JLLWrappers.@declare_library_product(libbzip2, "libbz2.so.1.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libbzip2,
        "lib/libbz2.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#-----------------------------------------------------------------------------------------------------
# Spin Operators
#-----------------------------------------------------------------------------------------------------
# Dictionary
function spin_dict(D::Integer)
    J = (D-1)/2
    coeff = [sqrt(i*(D-i)) for i = 1:D-1]
    sp = sparse(1:D-1, 2:D, coeff, D, D)
    sm = sp'
    sx = (sp+sm)/2
    sy = (sp-sm)/2
    sz = sparse(1:D, 1:D, J:-1:-J)
    s0 = sparse(1.0I, D, D)
    Dict([('+', sp), ('-', sm), ('x', sx), ('y', sy), ('z', sz), ('1', s0)])
end
# Atomic spin matrix
function spin_atom(s::String, dic::Dict)
    n = length(s)
    ny = sum(si == 'y' for si in s)
    temp = n == 1 ? dic[s[1]] : kron([dic[si] for si in s]...)
    P = mod(ny, 2) == 0 ? (-1)^(ny÷2) : (1im)^ny
    P * temp
end
# General spin matrix
function spin(spins; D::Integer=2)
    dic = spin_dict(D)
    res = sum(ci * spin_atom(si, dic) for (ci, si) in spins)
    Array(res)
end

#---------------------------------------------------------------------------------------------------
# Other properties
#---------------------------------------------------------------------------------------------------
export entropy
entropy(S::AbstractVector; cutoff::Real=1e-14) = -sum(i*log(i) for i in S if i > cutoff)

export inner_product
function inner_product(T1::Array, T2::Array)
    α, β = size(T1, 3), size(T2, 3)
    K = Kraus(T1, conj(T2), :r)
    ρ = rand(ComplexF64, α, β) |> vec
    e, v = eigsolve(K, ρ)
    abs(e[1])
end

inner_product(T::Array) = inner_product(T, T)


text:
#####
##### Utilities for launching kernels
#####

using CUDA
using KernelAbstractions
using Oceananigans.Architectures
using Oceananigans.Grids

const MAX_THREADS_PER_BLOCK = 256

flatten_reduced_dimensions(worksize, dims) = Tuple(i ∈ dims ? 1 : worksize[i] for i = 1:3)

"""
    work_layout(grid, dims; include_right_boundaries=false, location=nothing)

Returns the `workgroup` and `worksize` for launching a kernel over `dims`
on `grid`. The `workgroup` is a tuple specifying the threads per block in each
dimension. The `worksize` specifies the range of the loop in each dimension.

Specifying `include_right_boundaries=true` will ensure the work layout includes the
right face end points along bounded dimensions. This requires the field `location`
to be specified.

For more information, see: https://github.com/CliMA/Oceananigans.jl/pull/308
"""
function work_layout(grid, dims; include_right_boundaries=false, location=nothing, reduced_dimensions=())

    Nx, Ny, Nz = size(grid)

    workgroup = Nx == 1 && Ny == 1 ?

                    # One-dimensional column models:
                    (1, 1) :

                Nx == 1 ?

                    # Two-dimensional y-z slice models:
                    (1, min(MAX_THREADS_PER_BLOCK, Ny)) :

                Ny == 1 ?

                    # Two-dimensional x-z slice models:
                    (1, min(MAX_THREADS_PER_BLOCK, Nx)) :

                    # Three-dimensional models
                    (Int(√MAX_THREADS_PER_BLOCK), Int(√MAX_THREADS_PER_BLOCK))

    Nx′, Ny′, Nz′ = include_right_boundaries ? size(location, grid) : (Nx, Ny, Nz)

    Nx′, Ny′, Nz′ = flatten_reduced_dimensions((Nx′, Ny′, Nz′), reduced_dimensions)

    worksize = dims == :xyz ? (Nx′, Ny′, Nz′) :
               dims == :xy  ? (Nx′, Ny′) :
               dims == :xz  ? (Nx′, Nz′) :
               dims == :yz  ? (Ny′, Nz′) : throw(ArgumentError("Unsupported launch configuration: $dims"))

    return workgroup, worksize
end

"""
    launch!(arch, grid, layout, kernel!, args...; dependencies=nothing, kwargs...)

Launches `kernel!`, with arguments `args` and keyword arguments `kwargs`,
over the `dims` of `grid` on the architecture `arch`.

Returns an `event` token associated with the `kernel!` launch.

The keyword argument `dependencies` is an `Event` or `MultiEvent` specifying prior kernels
that must complete before `kernel!` is launched.
"""
function launch!(arch, grid, dims, kernel!, args...;
                 dependencies = nothing,
                 include_right_boundaries = false,
                 reduced_dimensions = (),
                 location = nothing,
                 kwargs...)

    workgroup, worksize = work_layout(grid, dims,
                                      include_right_boundaries = include_right_boundaries,
                                            reduced_dimensions = reduced_dimensions,
                                                      location = location)

    loop! = kernel!(Architectures.device(arch), workgroup, worksize)

    @debug "Launching kernel $kernel! with worksize $worksize"

    event = loop!(args...; dependencies=dependencies, kwargs...)

    return event
end


text:
# BCF Writer
# ==========

"""
    BCFWriter(output::IO, header::VCFHeader)

Create a data writer of the BCF file format.

# Arguments
* `output`: data sink
* `header`: VCF header object
"""
type BCFWriter{T<:IO} <: Bio.IO.AbstractWriter
    stream::BGZFStream{T}
end

function BCFWriter(output::IO, header::VCFHeader)
    stream = BGZFStream(output, "w")
    write(stream, b"BCF\x02\x02")
    buf = IOBuffer()
    len = write(buf, header)
    if len > typemax(Int32)
        error("too long header")
    end
    write(stream, htol(Int32(len)))
    data = takebuf_array(buf)
    @assert length(data) == len
    write(stream, data)
    return BCFWriter(stream)
end

function Bio.IO.stream(writer::BCFWriter)
    return writer.stream
end

function Base.write(writer::BCFWriter, record::BCFRecord)
    n = 0
    n += write(writer.stream, htol(record.sharedlen))
    n += write(writer.stream, htol(record.indivlen))
    n += write(writer.stream, record.data)
    return n
end


text:
module Config

export SimulationConfig, constant_scale_factor
export TimeStepOptions

struct SimulationConfig
    "function defining a(t)"
    a
    phase_grad_limit::Float64
    density_threshold::Float64
    time_step_options
end

function constant_scale_factor(t)
    1.0
end

function SimulationConfig(
                          ;
                          a=constant_scale_factor,
                          phase_grad_limit=π/4,
                          density_threshold=1e-6,
                          time_step_options=TimeStepOptions(),
                         )
    SimulationConfig(constant_scale_factor, phase_grad_limit, density_threshold, time_step_options)
end

"""
    TimeStepOptions(update_period=10, multiplier=1.0)

struct containing options controlling the size and calculation of time steps.

update_period::Int64 controls how many steps are taken before the timestep is
updated.

multiplier::Float64 multiplies the calculated maximum time step by a constant

See also: [`SimulationConfig`](@

# Examples
```jldoctest
julia> using UltraDark

julia> TimeStepOptions()
TimeStepOptions(10, 1.0)

```
"""
struct TimeStepOptions
    update_period::Int64
    multiplier::Float64
end

function TimeStepOptions(
                         ;
                         update_period=10,
                         multiplier=1.0,
                        )
    TimeStepOptions(update_period, multiplier)
end

end # module


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using IterativeSolvers
using LinearAlgebra
using Random
using SparseArrays

@testset "ILAN" begin

    @testset "Compute as many eigenpairs as possible (neigs=Inf)" begin
        n=100
        Random.seed!(1) # reset the random seed
        K = [1:n;2:n;1:n-1]; J=[1:n;1:n-1;2:n]
        A1 = sparse(K, J, rand(3*n-2)); A1 = A1+A1';
        A2 = sparse(K, J, rand(3*n-2)); A2 = A2+A2';
        A3 = sparse(K, J, rand(3*n-2)); A3 = A3+A3';
        nep=DEP([A1,A2,A3],[0,1,0.8])
        v0=rand(n)
        λ,W=ilan(nep,σ=0,γ=1;neigs=Inf,logger=0,maxit=30,tol=eps()*100,check_error_every=Inf,v=v0,errmeasure=ResidualErrmeasure(nep),inner_solver_method=NEPSolver.IARInnerSolver(tol=1e2,maxit=50))
        verify_lambdas(3, nep, λ, W, eps()*100)
    end

    @testset "Compute eigenpairs via Ritz extraction" begin
        nep=nep_gallery("dep_symm_double",10);
        v0=ones(size(nep,1));
        λ,W=ilan(nep;v=v0,tol=1e-5,neigs=3);
        verify_lambdas(3, nep, λ, W, 1e-5)
    end

    @testset "Errors thrown" begin
        n=100
        Random.seed!(1) # reset the random seed
        K = [1:n;2:n;1:n-1]; J=[1:n;1:n-1;2:n]
        A1 = sparse(K, J, rand(3*n-2)); A1 = A1+A1';
        A2 = sparse(K, J, rand(3*n-2)); A2 = A2+A2';
        A3 = sparse(K, J, rand(3*n-2)); A3 = A3+A3';
        nep=DEP([A1,A2,A3],[0,1,0.8])
        v0=rand(n)
        @test_throws NEPCore.NoConvergenceException λ,W,V2,H2,ω2,HH2=ilan(nep,σ=0,γ=1;neigs=2,logger=0,maxit=3,tol=eps()*100,check_error_every=Inf,v=v0,errmeasure=ResidualErrmeasure(nep))
    end

    @testset "Different format" begin
        n=100
        Random.seed!(1) # reset the random seed
        K = [1:n;2:n;1:n-1]; J=[1:n;1:n-1;2:n]
        A1 = sparse(K, J, rand(3*n-2)); A1 = A1+A1';
        A2 = sparse(K, J, rand(3*n-2)); A2 = A2+A2';
        nep1=DEP([A1,A2],[0,1])
        nep2=SPMF_NEP([one(A1),A1,A2],[S -> -S,S -> one(S),S -> exp(-S)]);

        v0=rand(n)
        λ,W,err,V,H,ω,HH=ilan(nep1,σ=0,γ=1;neigs=Inf,logger=0,maxit=10,tol=eps()*100,check_error_every=Inf,v=v0)

        λ2,W2,err2,V2,H2,ω2,HH2,=ilan(nep2,σ=0,γ=1;neigs=Inf,logger=0,maxit=10,tol=eps()*100,check_error_every=Inf,v=v0)
        #@test opnorm(W - W2) < 1e-6
        @test norm(λ - λ)<1e-6
        @test norm(W - W2)<1e-6
        @test norm(V - V2)<1e-6
        @test norm(H - H2)<1e-6
        @test norm(ω - ω2)<1e-6
        @test norm(HH - HH2)<1e-6
        @test opnorm(V'*V - I) < 1e-6
        @test opnorm(V2'*V2 - I) < 1e-6

    end
end


text:
# convert CorrData to SeisData
using SeisIO, SeisNoise, Dates, Logging
export corr2seis
""" Example of CorrData

CorrData with 1 Corrs
      NAME: "YW.FACR..BHZ.YW.FACR..BHZ"
        ID: "2010-01-01"
       LOC: 0.0 N, 0.0 E, 0.0 m
      COMP: "ZZ"
   ROTATED: false
 CORR_TYPE: "cross-correlation"
        FS: 5.0
      GAIN: 1.0
   FREQMIN: 0.1
   FREQMAX: 2.4
    CC_LEN: 3600
   CC_STEP: 1800
  WHITENED: false
 TIME_NORM: false
      RESP: c = 0.0, 0 zeros, 0 poles
      MISC: 2 entries
     NOTES: 3 entries
    MAXLAG: 100.0
         T: 2010-01-01T00:00:00.000            …
      CORR: 1001×1 Array{Float32,2}
"""

""" Example of CorrData.misc
Dict{String,Any} with 27 entries:
  "e"        => 86399.8
  "dist"     => 0.0
  "scale"    => 1.0
  "location" => Dict("YW.FACO..BHZ"=>0.0 N, 0.0 E, 0.0 m,"TA.G05D..BHZ"=>0.0 N, 0.0 E, 0.0 m)
  "iztype"   => 9
  "khole"    => "       \xb0"
  "lcalda"   => 0
  "nzhour"   => 0
  "nvhdr"    => 6
  "nzjday"   => 40
  "kcmpnm"   => "BHZ\0\0\0\0\0"
  "depmin"   => -7.54567e-8
  "iftype"   => 1
  "depmax"   => 9.53069e-8
  "nzmin"    => 0
  "b"        => 0.0
  "delta"    => 0.2
  "npts"     => 432000
  "kstnm"    => "G05D\0\0\0\0"
  "nzyear"   => 2010
  "nzmsec"   => 0
  "leven"    => 1
  "lpspol"   => 1
  "depmen"   => -7.11344e-16
  "nzsec"    => 0
  "lovrok"   => 1
  "knetwk"   => "TA\0\0\0\0\0\0"
"""

""" Example of SeisData

SeisChannel with 432000 samples
    ID: YW.FACO..BHZ
  NAME:
   LOC: 0.0 N, 0.0 E, 0.0 m
    FS: 5.0
  GAIN: 1.0
  RESP: c = 0.0, 0 zeros, 0 poles
 UNITS:
   SRC: /Users/xiaotaoyang/Work/Cascadia/Raw/Event_2010_001/YW.FACO.BHZ.sac
  MISC: 0 entries
 NOTES: 1 entries
     T: 2010-01-01T00:00:00.010 (0 gaps)
     X: +1.595e-18
        +1.421e-14
            ...
        +6.933e-11
        (nx = 432000)

"""
# richmisc: true if include location
function corr2seis(C::CorrData)
    #get the size of CorrData, to determine how many channels needed.

    #create an empty SeisData object with the right size.
    S = SeisData(size(C.corr,2))

    #loop through the CorrData (corr) and pass the header values to SeisData object
    # in loop, save customzed headers to misc dictionary.
    # save freqmin and freqmax twice:
    # 1. use real field names. 2. use USER0 and USER1, to be saved to SAC for later use.
    stemp = split(C.name,".")
    evname = join([stemp[1],stemp[2],stemp[3],stemp[4]],".")  #use the first entry in location as the event name
    stname = join([stemp[5],stemp[6],stemp[7],stemp[8]],".")
    evla = 0.0
    evlo = 0.0
    evdp = 0.0
    stla = 0.0
    stlo = 0.0
    stel = 0.0
    if in("location",collect(keys(C.misc)))
        # println(C.name,": ",collect(keys(C.misc["location"])))
        # evname = collect(keys(C.misc["location"]))[1]  #use the first entry in location as the event name
        # stname = collect(keys(C.misc["location"]))[2]
        evla = C.misc["location"][evname].lat
        evlo = C.misc["location"][evname].lon
        evdp = C.misc["location"][evname].el
        stla = C.misc["location"][stname].lat
        stlo = C.misc["location"][stname].lon
        stel = C.misc["location"][stname].el
    else
        @warn "$(C.name) : 'location' field not defined in misc. Use default locations (0.0)."
    end

    for i = 1:S.n
        SC = SeisChannel();
        # for tv in C.corr[:,i]
        #     push!(SC.x,tv)
        # end
        SC.x = C.corr[:,i]
        SC.id = C.id
        SC.name = C.name
        setfield!(SC.loc,:lat,stla)
        setfield!(SC.loc,:lon,stlo)
        setfield!(SC.loc,:dep,stel)
        # = C.loc
        SC.fs = C.fs
        SC.gain = C.gain
        SC.resp = C.resp
        SC.src = C.name
        SC.t = [1 Int64(C.t[i]*1e6);size(C.corr,1) 0]
        SC.notes = C.notes
        # TO-DO get misc values
        SC.misc = C.misc

        #update some values in misc for SAC headers.
        SC.misc["b"] = -1.0*C.maxlag
        SC.misc["e"] = C.maxlag
        SC.misc["delta"] = 1.0/C.fs
        SC.misc["npts"] = size(C.corr,1)
        SC.misc["depmin"] = minimum(C.corr[:,i])
        SC.misc["depmax"] = maximum(C.corr[:,i])
        SC.misc["depmen"] = sum(C.corr[:,i])/size(C.corr,1)
        SC.misc["kevnm"] = evname
        SC.misc["evla"] = evla
        SC.misc["evlo"] = evlo
        SC.misc["evdp"] = evdp
        SC.misc["stla"] = stla
        SC.misc["stlo"] = stlo
        SC.misc["stel"] = stel
        SC.misc["nzyear"] = Dates.year(u2d(C.t[i]))
        SC.misc["nzjday"] = Dates.date2epochdays(Date(u2d(C.t[i]))) - Dates.date2epochdays(Date(string(Dates.year(u2d(C.t[i])))*"-01-01")) + 1
        SC.misc["nzhour"] = Dates.hour(u2d(C.t[i]))
        SC.misc["nzmin"] = Dates.minute(u2d(C.t[i]))
        SC.misc["nzsec"] = Dates.second(u2d(C.t[i]))
        SC.misc["nzmsec"] = Dates.millisecond(u2d(C.t[i]))
        SC.misc["comp"] = C.comp
        SC.misc["rotated"] = C.rotated
        SC.misc["corr_type"] =  C.corr_type
        # SC.misc["kuser0"] =  C.corr_type
        SC.misc["freqmin"] = C.freqmin
        SC.misc["freqmax"] = C.freqmax
        # SC.misc["user0"] = C.freqmin
        # SC.misc["user1"] = C.freqmax
        SC.misc["cc_len"] = C.cc_len
        SC.misc["cc_step"] = C.cc_step
        # SC.misc["user2"] = C.cc_len
        # SC.misc["user3"] = C.cc_step
        # SC.misc["whitened"] = C.whitened
        # SC.misc["time_norm"] = C.time_norm
        # if C.whitened
        #     SC.misc["kuser1"] =  "NOWHITND"
        # else
        #     SC.misc["kuser1"] =  "WHITENED"
        # end
        # if C.time_norm
        #     SC.misc["kuser2"] =  "TIMENORM"
        # else
        #     SC.misc["kuser2"] =  "NOTIMENM"
        # end
        if C.whitened
            SC.misc["whitened"] =  "WHITENED"
        else
            SC.misc["whitened"] =  "NOWHITND"
        end
        if typeof(C.time_norm) == Bool
            if C.time_norm
                SC.misc["time_norm"] =  "TIMENORM"
            else
                SC.misc["time_norm"] =  "NOTIMENM"
            end
        else
            SC.misc["time_norm"] =  C.time_norm
        end
        S[i] = SC
    end

    return(S)
end


text:
# MetidaNLopt
# Copyright © 2019-2020 Vladimir Arnautov aka PharmCat <mail@pharmcat.net>
using Metida
using MetidaNLopt
using  Test, CSV, DataFrames, StatsModels, CategoricalArrays

df0 = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","df0.csv"); types = [String, String, String, String, Float64, Float64]) |> DataFrame
transform!(df0, :subject => categorical, renamecols=false)
transform!(df0, :period => categorical, renamecols=false)
transform!(df0, :sequence => categorical, renamecols=false)
transform!(df0, :formulation=> categorical, renamecols=false)

ftdf         = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","1fptime.csv"); types = [String, String, Float64, Float64]) |> DataFrame

ftdf2        = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","1freparma.csv"); types = [String, String, Float64, Float64]) |> DataFrame

ftdf3        = CSV.File(joinpath(dirname(pathof(Metida)),"..","test","csv","2f2rand.csv"); types =
[String,  Float64, Float64, String, String, String, String, String]) |> DataFrame

@testset "  Basic test                                               " begin
    lmm = LMM(@formula(var~sequence+period+formulation), df0;
    random = VarEffect(@covstr(formulation|subject), CSH),
    repeated = VarEffect(@covstr(formulation|subject), DIAG),
    )
    fit!(lmm; solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    #Metida.fit_nlopt!(lmm; solver = :nlopt, rholinkf = :sigm, f_tol=0.0, x_tol = 0.0, f_rtol =0.0, x_rtol =1e-18)
    #Metida.m2logreml(lmm) ≈ 10.065239006121315
    #10.065239006121315
    #10.065456008797781
    @test Metida.m2logreml(lmm) ≈ 10.065238620486195 atol=1E-6

    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
    repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CS),
    )
    Metida.fit!(lmm, solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm) ≈ 10.3039977509049 atol=1E-6

    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
    repeated = Metida.VarEffect(Metida.@covstr(period|subject), Metida.CSH),
    )
    Metida.fit!(lmm; solver = :nlopt, f_tol=1e-16, x_tol=1e-16,
    )
    @test Metida.m2logreml(lmm) ≈ 8.740095378772942 atol=1E-8

    lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
    random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
    )
    Metida.fit!(lmm; solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm) ≈ 1300.1807598168923 atol=1E-6

    lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;
    repeated = Metida.VarEffect(Metida.@covstr(time|subject&factor), Metida.ARMA),
    )
    Metida.fit!(lmm; solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm) ≈ 715.4528559688382 atol = 1E-6

    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = [Metida.VarEffect(Metida.@covstr(r1|subject), Metida.CS), Metida.VarEffect(Metida.@covstr(r2|subject), Metida.CS)],
    )
    Metida.fit!(lmm; solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm)  ≈ 710.4250214813896 atol=1E-8

    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = Metida.VarEffect(Metida.@covstr(r1|subject), Metida.AR),
    repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.DIAG),
    )
    Metida.fit!(lmm; solver = :nlopt, f_tol=1e-10, x_tol=1e-10)
    @test Metida.m2logreml(lmm)  ≈ 698.8792511057682 atol=1E-8

    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3; contrasts=Dict(:factor => DummyCoding(; base=1.0)),
    random = Metida.VarEffect(Metida.@covstr(p|r1&r2), Metida.ARMA),
    )
    Metida.fit!(lmm, solver = :nlopt)
    @test Metida.m2logreml(lmm)  ≈ 913.9176298311813 atol=1E-8

    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(1 + r1 + r2|subject), Metida.TOEPHP(3); coding = Dict(:r1 => DummyCoding(), :r2 => DummyCoding())),
    )
    Metida.fit!(lmm, solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm)  ≈ 705.9946274598822 atol=1E-8

    lmm = Metida.LMM(@formula(response ~ 1 + factor), ftdf3;
    random = Metida.VarEffect(Metida.@covstr(r2|subject), Metida.DIAG),
    repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.TOEPP(3)),
    )
    Metida.fit!(lmm, solver = :nlopt, f_tol=1e-16, x_tol=1e-16)
    @test Metida.m2logreml(lmm)  ≈ 773.9575538254085 atol=1E-8
end


text:
export marginalrule

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, m_Λ::PointMass) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_μ), mean(m_Λ)), m_out), μ = m_μ, Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, m_Λ::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, MvNormalMeanPrecision(mean(m_out), mean(m_Λ))), Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, m_Λ::PointMass) = begin
    xi_y, W_y  = weightedmean_precision(m_out)
    xi_m, W_m  = weightedmean_precision(m_μ)

    W_bar = mean(m_Λ)

    Λ  = [ W_y + W_bar -W_bar; -W_bar W_m + W_bar ]
    ξ  = [ xi_y; xi_m ]

    return (out_μ = MvNormalWeightedMeanPrecision(ξ, Λ), Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, q_Λ::Any) = begin
    xi_y, W_y  = weightedmean_precision(m_out)
    xi_m, W_m  = weightedmean_precision(m_μ)

    W_bar = mean(q_Λ)

    Λ  = [ W_y + W_bar -W_bar; -W_bar W_m + W_bar ]
    ξ  = [ xi_y; xi_m ]

    return MvNormalWeightedMeanPrecision(ξ, Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, q_Λ::Any) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_out), mean(q_Λ)), m_μ))
end

@marginalrule MvNormalMeanPrecision(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, q_Λ::Any) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanPrecision(mean(m_μ), mean(q_Λ)), m_out), μ = m_μ)
end

## ProdFinal / BIFM related

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::ProdFinal, m_μ::PointMass, m_Λ::PointMass) = begin
    return (out = getdist(m_out), μ = m_μ, Λ = m_Λ)
end

@marginalrule MvNormalMeanPrecision(:out_μ_Λ) (m_out::PointMass, m_μ::ProdFinal, m_Λ::PointMass) = begin
    return (out = m_out, μ = getdist(m_μ), Λ = m_Λ)
end


text:
"""
    AbstractNumbers

An abstract type to represent the number system on which a manifold is built.

This provides concrete number types for dispatch. The two most common number types are
the fields [`RealNumbers`](@ref) (`ℝ` for short) and [`ComplexNumbers`](@ref) (`ℂ`).
"""
abstract type AbstractNumbers end

"""
    ℝ = RealNumbers()

The field of real numbers.
"""
struct RealNumbers <: AbstractNumbers end

"""
    ℂ = ComplexNumbers()

The field of complex numbers.
"""
struct ComplexNumbers <: AbstractNumbers end

"""
    ℍ = QuaternionNumbers()

The division algebra of quaternions.
"""
struct QuaternionNumbers <: AbstractNumbers end

const ℝ = RealNumbers()
const ℂ = ComplexNumbers()
const ℍ = QuaternionNumbers()

Base.show(io::IO, ::RealNumbers) = print(io, "ℝ")
Base.show(io::IO, ::ComplexNumbers) = print(io, "ℂ")
Base.show(io::IO, ::QuaternionNumbers) = print(io, "ℍ")

^(𝔽::AbstractNumbers, n) = Euclidean(n...; field = 𝔽)

@doc doc"""
    real_dimension(𝔽::AbstractNumbers)

Return the real dimension $\dim_ℝ 𝔽$ of the [`AbstractNumbers`] system `𝔽`.
The real dimension is the dimension of a real vector space with which a number in `𝔽` can be
identified.
For example, [`ComplexNumbers`](@ref) have a real dimension of 2, and
[`QuaternionNumbers`](@ref) have a real dimension of 4.
"""
function real_dimension(𝔽::AbstractNumbers)
    error("real_dimension not defined for number system $(𝔽)")
end
real_dimension(::RealNumbers) = 1
real_dimension(::ComplexNumbers) = 2
real_dimension(::QuaternionNumbers) = 4


text:
module photovoltaic
struct pv_module
    module_type # 模块类型
    stc_watts
    stc_eff
    ff
    aspect_ratio
    width
    length
    area
    vmp
    ndiode
    gamma
    ar_glass
    bifaciality
    function pv_module(;module_type::Int=0,bifaciality::Float64=0.001
    )
        new(
        tilt,
        array_type,
        module_type,
        losses,
        azimuth,
        gcr,
        radius,
        name,
        location,
        existing_kw,
        min_kw,
        max_kw,
        installed_cost_per_kw,
        om_cost_per_kw,
        degradation_pct,
        macrs_option_years,
        macrs_bonus_pct,
        macrs_itc_reduction,
        kw_per_square_foot,
        acres_per_kw,
        inv_eff,
        dc_ac_ratio,
        prod_factor_series_kw,
        federal_itc_pct,
        federal_rebate_per_kw,
        state_ibi_pct,
        state_ibi_max,
        state_rebate_per_kw,
        state_rebate_max,
        utility_ibi_pct,
        utility_ibi_max,
        utility_rebate_per_kw,
        utility_rebate_max,
        production_incentive_per_kwh,
        production_incentive_max_benefit,
        production_incentive_years,
        production_incentive_max_kw
    )
    end
end

struct pv_array
    array_type

    dc_nameplate
    dc_ac_ratio
    ac_nameplate
    xfmr_rating
    ac_plant_max
    inv_eff_percent
    dc_loss_percent
    tilt
    azimuth
    rotlim
    function pv_array(;)
        new(

        )
    end
end

function prodfactor(pv::PV,latitude::Real, longitude::Real; timeframe="hourly")

end

end

max_stars_repo_name: mingtaoli/hresOpt

max_issues_repo_name: mingtaoli/hresOpt

max_forks_repo_name: mingtaoli/hresOpt


text:
module Wire
    mutable struct Sender
        acc::Int
    end
    mutable struct Receiver
        acc::Int
    end

    # test No Overdrafts
    # @test "No overdrafts" begin
    #     all(p.val >= 0 for p in people)
    # end

    function withdraw!(x::Sender, amount)
       x.acc -= amount
    end

    function deposit!(x::Receiver, amount)
        x.acc += amount
    end

    function transact!(x::Sender, y::Receiver, amount)
        @sync begin
            @async withdraw!(x, amount)
            @async deposit!(y, amount)
        end
    end

    function main(account = 5, amount = 3)
        alice = Sender(account)
        bob = Receiver(account)

        transact!(alice, bob, amount)

        return all(p.acc >= 0 for p in (alice, bob))
    end

end

using .Wire
Wire.main()


text:
struct MvNormal
    mu::Vector{Float64}
    Sigma::Matrix{Float64}
end

struct Cartesian
    a::Vector{Float64}
    b::Vector{Float64}
    orders::Vector{Int64}
end

@testset "minilanguage" begin

import Dolang: Language, construct_object, add_language_elements!, eval_node
import Dolang: FromGreek, ToGreek
import Dolang: LangInvalidCall, LangUnknownTag, ParsingError, yaml_node_from_file
# the two objects replace  imports from Dolo

# we create a mini-language with two keywords:
# UNormal(μ=0.0, σ=1.0) # keyword object with greek arguments
# Cartesian # positional

UNormal(;μ::Float64=0.0,σ::Float64=1.0) = MvNormal([μ], σ*ones(1,1))
minilang = Language(Dict())
add_language_elements!(minilang, Dict(
    "!Cartesian"=>Cartesian,
    "!UNormal"=>UNormal
))

obj = construct_object(minilang, "!UNormal", Dict(:σ=>0.1, :μ=> 0.001))
@assert typeof(obj) <: MvNormal

obj = construct_object(minilang, "!Cartesian", ([0.0], [1.0], [20]))
@assert typeof(obj) <: Cartesian

# UNormal doesn't take positional arguments:
@test_throws LangInvalidCall construct_object(minilang, "!UNormal", (0.1, 0.001))
@test_throws LangInvalidCall construct_object(minilang, "!Cartesian", ([0.0], [1.0]) )
@test_throws LangUnknownTag construct_object(minilang, "!NonExisting", ([0.0], [1.0], [20]) )

##

calibration = Dict(:a=>0.2, :b=>0.5)

# use test file

data = yaml_node_from_file("minilanguage.yaml", minilang)

# evaluate terminal nodes directly
res = eval_node( data["equations"][1], calibration )
@assert res == 0.7

# if expression is not recognized or calibration not given, expression is returned as such
res = eval_node( data["equations"][2], calibration )
@assert res=="a + b + c"

# one can evaluate a yaml list directly
res = eval_node( data["equations"], calibration)
@assert res[1] == 0.7
@assert res[2] == "a + b + c"
@assert res[3] == 8
@assert res[4] == 9.0 # we return same types as parsed by yaml

# this raises a parsing exception
@test_throws ParsingError{LangUnknownTag} eval_node(data["exogenous"], calibration)
@test_throws ParsingError{LangUnknownTag} eval_node(data["shocks"], calibration)

# evaluate typed section
gg = eval_node(data["grid"], calibration, minilang)
@test typeof(gg) <: Cartesian
gg = eval_node(data["exogenous"], calibration, minilang)
@test typeof(gg) <: MvNormal

# test greek tolerance

gg = eval_node(data["calibration"], calibration, minilang)
@test tuple(sort([keys(gg)...])... ) == (:mu, :β, :λ, :σ)

gg = eval_node(data["calibration"], calibration, minilang, FromGreek())
@test tuple(sort([keys(gg)...])... ) == ( :beta, :lambda, :mu, :sigma)

@test_throws ParsingError{LangInvalidCall} eval_node(data["exogenous2"], calibration, minilang)
@test_throws ParsingError{LangInvalidCall} eval_node(data["exogenous2"], calibration, minilang, FromGreek())

# that one works because object !
gg = eval_node(data["exogenous2"], calibration, minilang, ToGreek())
typeof(gg) <: MvNormal

end

max_issues_repo_name: EconForge/Dolang


text:
using SafeTestsets
const LONGER_TESTS = false

if haskey(ENV,"GROUP")
    group = ENV["GROUP"]
else
    group = "All"
end

is_APPVEYOR = ( Sys.iswindows() && haskey(ENV,"APPVEYOR") )

#Start Test Script

@time begin
if group == "All" || group == "Interface"
  @time @safetestset "Discrete Algorithm Tests" begin include("discrete_algorithm_test.jl") end
  @time @safetestset "Tstops Tests" begin include("ode/ode_tstops_tests.jl") end
  @time @safetestset "Backwards Tests" begin include("ode/ode_backwards_test.jl") end
  @time @safetestset "Initdt Tests" begin include("ode/ode_initdt_tests.jl") end
  @time @safetestset "Mass Matrix Tests" begin include("mass_matrix_tests.jl") end
  @time @safetestset "Differentiation Trait Tests" begin include("differentiation_traits_tests.jl") end
  @time @safetestset "Inf Tests" begin include("inf_handling.jl") end
  @time @safetestset "saveat Tests" begin include("ode/ode_saveat_tests.jl") end
  @time @safetestset "save_idxs Tests" begin include("ode/ode_saveidxs_tests.jl") end
  @time @safetestset "Static Array Tests" begin include("static_array_tests.jl") end
  @time @safetestset "Data Array Tests" begin include("data_array_test.jl") end
  @time @safetestset "u_modifed Tests" begin include("umodified_test.jl") end
  @time @safetestset "Composite Algorithm Tests" begin include("composite_algorithm_test.jl") end
  @time @safetestset "Complex Tests" begin include("complex_tests.jl") end
  @time @safetestset "Stiffness Detection Tests" begin include("stiffness_detection_test.jl") end
  @time @safetestset "Composite Interpolation Tests" begin include("composite_interpolation.jl") end
  @time @safetestset "Export tests" begin include("export_tests.jl") end
  @time @safetestset "Derivative Utilities Tests" begin include("utility_tests.jl") end
  @time @safetestset "Discrete Callback Dual Tests" begin include("discrete_callback_dual_test.jl") end
end

if group == "All" || group == "Integrators"
  @time @safetestset "Reinit Tests" begin include("reinit_test.jl") end
  @time @safetestset "Events Tests" begin include("ode/ode_event_tests.jl") end
  @time @safetestset "Cache Tests" begin include("ode/ode_cache_tests.jl") end
  @time @safetestset "Iterator Tests" begin include("iterator_tests.jl") end
  @time @safetestset "Integrator Interface Tests" begin include("integrator_interface_tests.jl") end
  @time @safetestset "Add Steps Tests" begin include("ode/ode_add_steps_tests.jl") end
  @time @safetestset "Error Check Tests" begin include("check_error.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "Regression" )
  @time @safetestset "Linear Tests" begin include("ode/ode_twodimlinear_tests.jl") end
  @time @safetestset "Dense Tests" begin include("ode/ode_dense_tests.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "AlgConvergence_I" )
  # ~ 250 s
  @time @safetestset "Partitioned Methods Tests" begin include("partitioned_methods_tests.jl") end
  # ~ 400 s
  @time @safetestset "Convergence Tests" begin include("ode/ode_convergence_tests.jl") end
  # ~ 2 s
  @time @safetestset "Adams Variable Coefficients Tests" begin include("ode/adams_tests.jl") end
  # ~ 50 s
  @time @safetestset "Nordsieck Tests" begin include("ode/nordsieck_tests.jl") end
  @time @safetestset "Linear Methods Tests" begin include("linear_method_tests.jl") end
  # ~ 170 s
  @time @safetestset "SSPRK Tests" begin include("ode/ode_ssprk_tests.jl") end
  # ~ 25 s
  @time @safetestset "OwrenZen Tests" begin include("owrenzen_tests.jl") end
  @time @safetestset "Runge-Kutta-Chebyshev Tests" begin include("ode/rkc_tests.jl") end
end

if !is_APPVEYOR && ( group == "All" || group == "AlgConvergence_II" )
  # ~ 110 s
  @time @safetestset "Split Methods Tests" begin include("split_methods_tests.jl") end
  # ~ 550 s
  @time @safetestset "Rosenbrock Tests" begin include("ode/ode_rosenbrock_tests.jl") end
  # ~ 40 s
  @time @safetestset "Linear-Nonlinear Methods Tests" begin include("linear_nonlinear_convergence_tests.jl") end
  # ~ 140 s
  @time @safetestset "Linear-Nonlinear Krylov Methods Tests" begin include("linear_nonlinear_krylov_tests.jl") end
end
end # @time


text:
@testset "876.middle-of-the-linked-list.jl" begin
    head = ListNode{Int}([1, 2, 3, 4, 5])
    expected = head.next.next
    @test middle_node_by_double_pointers(head) === expected

    head = ListNode{Int}([1, 2, 3, 4, 5, 6])
    expected = head.next.next.next
    @test middle_node_by_double_pointers(head) === expected
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using OrdinaryDiffEq, DynamicalSystems
using LinearAlgebra: mul!
import FFTW

# In the following `u` is the field in real space, while `y` is the Fourier
# transform of the field, which lives in inverse space.

b1 = 30 # total domain size
T = 500.0
dt = 0.25 # timestep for integration
saveat = 0:dt:T
dt = 0.25
dx = 0.2 # spatial discretization
xs = range(0, b1, step = dx) # space
u0 = @. cos(xs) + 0.1*sin(xs/8) + 0.01*cos((2π/b)*xs)
ks = Vector(FFTW.rfftfreq(length(u0))/dx) # conjugate space (wavenumbers)

forward_plan = FFTW.plan_rfft(u0)
y0 = forward_plan * u0
inverse_plan = FFTW.plan_irfft(y0, length(u0))
ik2 = -im .* ks ./ 2
k²_k⁴ = @. ks^2 - ks^4

ydummy = copy(y0)
udummy = copy(u0)
ksparams = (forward_plan, inverse_plan, udummy, ydummy, k²_k⁴, ik2)

function kse_spectral!(dy, y, p, t)
    forward_plan, inverse_plan, udummy, ydummy, k²_k⁴, ik2 = p
    y² = begin # nonlinear term
        mul!(udummy, inverse_plan, y) # create current u in x-space
        udummy .*= udummy # square current u
        mul!(ydummy, forward_plan, udummy) # transform to k-space
    end
    # KS equation in spectral space
    @. dy = y²*ik2 + k²_k⁴*y
    return nothing
end

prob = ODEProblem(kse_spectral!, y0, (0.0, T), ksparams)
@time sol = solve(prob, Tsit5(); saveat)

u = [inverse_plan*y for y in sol.u]
U = hcat(u...)

# %% Make figure
fig, axs = subplots(1,2)
pm = axs[1].pcolormesh(saveat, xs, U, cmap = "PRGn", vmin = -3, vmax = 3)
cb = colorbar(pm, ax = axs[1])
axs[1].set_xticks([0, saveat[end]])
axs[1].set_yticks([0, b1])
axs[1].set_xlabel("\$t\$"; labelpad = -20)
axs[1].set_ylabel("\$x\$"; labelpad = -20)

cb.set_ticks([-3, 3])
cb.set_label("\$u\$", labelpad = -40)

# load exponents
λs = []

N = 4000    # steps to run the QR normalization
Δt = 1.0    # timestep between normalizations
dx = 0.2    # spatial discretization
bs = 20:20:60

for b in bs
    try
        N = 4000
        config = @strdict N dx Δt b
        file = wload(datadir("ksiva", savename("ksiva_spectrum", config, "jld2")))
        push!(λs, file["λs"])
    catch
        N = 1000
        config = @strdict N dx Δt b
        file = wload(datadir("ksiva", savename("ksiva_spectrum", config, "jld2")))
        push!(λs, file["λs"])
    end
end

N = 4000    # steps to run the QR normalization
for dx in (0.15, 0.1)
    b = 20
    config = @strdict N dx Δt b
    file = wload(datadir("ksiva", savename("ksiva_spectrum", config, "jld2")))
    push!(λs, file["λs"])
end

axs[2].axhline(0; color = "C3", lw = 1)

for i in 1:length(λs)
    l = λs[i]
    if i ≤ length(bs)
        axs[2].plot(1:length(l), cumsum(l); color = "C$(i-1)",
        label = "\$b=$(bs[i])\$")

        Δ = kaplanyorke_dim(sort(l; rev = true))
        axs[2].plot([Δ, Δ], [-25, 0]; color = "C$(i-1)", marker = "o", ms = 12,
        mec = "C$(i-1)", mfc = :white, lw = 1, zorder = 99)
    else
        axs[2].plot(1:length(l), cumsum(l); color = "C0",
        ls = i == length(bs) + 1 ? "--" : ":", zorder = 1)
    end
end

axs[2].set_ylim(-5, 5)
axs[2].set_xticks(0:40:120)
axs[2].set_xlim(0, 120)
axs[2].set_xlabel("\$n\$"; labelpad = -20)
axs[2].set_ylabel("\$\\Lambda_n\$", labelpad = -10)
axs[2].set_ylim(-5, 5)
axs[2].legend(ncol=1, loc = "upper right", fontsize = 24, handlelength = 1)
add_identifiers!(fig, axs)
fig.tight_layout(pad=0.3)
wsave(plotsdir("11", "ksiva"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# -----------------------------------------------------------------------------
##### CUBOID #####-------------------------------------------------------------
function subregioncuboid(dataobject::HydroDataType;
    xrange::Array{<:Any,1}=[missing, missing],
    yrange::Array{<:Any,1}=[missing, missing],
    zrange::Array{<:Any,1}=[missing, missing],
    center::Array{<:Any,1}=[0., 0., 0.],
    range_unit::Symbol=:standard,
    cell::Bool=true,
    inverse::Bool=false,
    verbose::Bool=verbose_mode)

    printtime("", verbose)

    boxlen = dataobject.boxlen
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    # convert given ranges and print overview on screen
    ranges = prepranges(dataobject.info,range_unit, verbose, xrange, yrange, zrange, center)

    xmin, xmax, ymin, ymax, zmin, zmax = ranges

    #if !(xrange == [dataobject.ranges[1], dataobject.ranges[2]] &&
    #   yrange == [dataobject.ranges[3], dataobject.ranges[4]] &&
    #   zrange == [dataobject.ranges[5], dataobject.ranges[6]])
    if !(xrange[1] === missing && xrange[2] === missing &&
         yrange[1] === missing && yrange[2] === missing &&
         zrange[1] === missing && zrange[2] === missing)

        if inverse == false
            if isamr
                if cell == true
                    sub_data = filter(p->   p.cx >=floor(Int, 2^p.level * xmin) &&
                                            p.cx <=ceil(Int,  2^p.level * xmax) &&
                                            p.cy >=floor(Int, 2^p.level * ymin) &&
                                            p.cy <=ceil(Int,  2^p.level * ymax) &&
                                            p.cz >=floor(Int, 2^p.level * zmin) &&
                                            p.cz <=ceil(Int,  2^p.level * zmax), dataobject.data)
                else
                    sub_data = filter(p->   p.cx >=( 2^p.level * xmin) &&
                                            p.cx <=( 2^p.level * xmax) &&
                                            p.cy >=( 2^p.level * ymin) &&
                                            p.cy <=( 2^p.level * ymax) &&
                                            p.cz >=( 2^p.level * zmin) &&
                                            p.cz <=( 2^p.level * zmax), dataobject.data)
                end
            else # for uniform grid
                if cell == true
                    sub_data = filter(p->   p.cx >=floor(Int, 2^lmax * xmin) &&
                                            p.cx <=ceil(Int,  2^lmax * xmax) &&
                                            p.cy >=floor(Int, 2^lmax * ymin) &&
                                            p.cy <=ceil(Int,  2^lmax * ymax) &&
                                            p.cz >=floor(Int, 2^lmax * zmin) &&
                                            p.cz <=ceil(Int,  2^lmax * zmax), dataobject.data)
                else
                    sub_data = filter(p->   p.cx >=( 2^lmax * xmin) &&
                                            p.cx <=( 2^lmax * xmax) &&
                                            p.cy >=( 2^lmax * ymin) &&
                                            p.cy <=( 2^lmax * ymax) &&
                                            p.cz >=( 2^lmax * zmin) &&
                                            p.cz <=( 2^lmax * zmax), dataobject.data)
                end

            end
        elseif inverse == true
            if cell == true
                if isamr
                    sub_data = filter(p->   p.cx <floor(Int, 2^p.level * xmin) ||
                                            p.cx >ceil(Int,  2^p.level * xmax) ||
                                            p.cy <floor(Int, 2^p.level * ymin) ||
                                            p.cy >ceil(Int,  2^p.level * ymax) ||
                                            p.cz <floor(Int, 2^p.level * zmin) ||
                                            p.cz >ceil(Int,  2^p.level * zmax) , dataobject.data)
                else # for uniform grid
                    sub_data = filter(p->   p.cx <floor(Int, 2^lmax * xmin) ||
                                            p.cx >ceil(Int,  2^lmax * xmax) ||
                                            p.cy <floor(Int, 2^lmax * ymin) ||
                                            p.cy >ceil(Int,  2^lmax * ymax) ||
                                            p.cz <floor(Int, 2^lmax * zmin) ||
                                            p.cz >ceil(Int,  2^lmax * zmax) , dataobject.data)
                end
            else
                if isamr
                    sub_data = filter(p->   p.cx <( 2^p.level * xmin) ||
                                            p.cx >( 2^p.level * xmax) ||
                                            p.cy <( 2^p.level * ymin) ||
                                            p.cy >( 2^p.level * ymax) ||
                                            p.cz <( 2^p.level * zmin) ||
                                            p.cz >( 2^p.level * zmax), dataobject.data)
                else # for uniform grid
                    sub_data = filter(p->   p.cx <( 2^lmax * xmin) ||
                                            p.cx >( 2^lmax * xmax) ||
                                            p.cy <( 2^lmax * ymin) ||
                                            p.cy >( 2^lmax * ymax) ||
                                            p.cz <( 2^lmax * zmin) ||
                                            p.cz >( 2^lmax * zmax), dataobject.data)
                end
            end

                ranges = dataobject.ranges

        end

        printtablememory(sub_data, verbose)

        hydrodata = HydroDataType()
        hydrodata.data = sub_data
        hydrodata.info = dataobject.info
        hydrodata.lmin = dataobject.lmin
        hydrodata.lmax = dataobject.lmax
        hydrodata.boxlen = dataobject.boxlen
        hydrodata.ranges = ranges
        hydrodata.selected_hydrovars = dataobject.selected_hydrovars
        hydrodata.used_descriptors = dataobject.used_descriptors
        hydrodata.smallr = dataobject.smallr
        hydrodata.smallc = dataobject.smallc
        hydrodata.scale = dataobject.scale
        return hydrodata

    else
        return dataobject
        #println("[Mera]: Nothing to do! Given ranges match data ranges!")
    end

end

# -----------------------------------------------------------------------------
##### CYLINDER #####-----------------------------------------------------------

function cell_shift(level, value_shift, cell)

    if cell == false
        return 2^level * value_shift
    else
        return (ceil(Int, 2^level * value_shift))
    end
end

function get_radius_cylinder(x,y, level, cx_shift, cy_shift)
    xcenter = 2^level * cx_shift
    ycenter = 2^level * cy_shift
    x = (x - xcenter)
    y = (y - ycenter)
    return sqrt( x^2 + y^2)
end

function get_height_cylinder(z, level, cz_shift)
    center = 2^level * cz_shift
    z = (z - center)
    return abs(z)
end

function subregioncylinder(dataobject::HydroDataType;
                            radius::Real=0.,
                            height::Real=0.,
                            center::Array{<:Any,1}=[0., 0., 0.],
                            range_unit::Symbol=:standard,
                            direction::Symbol=:z,
                            cell::Bool=true,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    if radius == 0. || height == 0. || in(0., center)
        error("[Mera]: given radius, height or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_shift, height_shift = prepranges(dataobject.info, center, radius, height, range_unit, verbose)

    if inverse == false
        if isamr
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                <= (cell_shift(p.level, radius_shift,cell))  &&

                                get_height_cylinder(p.cz, p.level, cz_shift)
                                <= (cell_shift(p.level, height_shift,cell)),

                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                <= (cell_shift(lmax, radius_shift,cell))  &&

                                get_height_cylinder(p.cz, lmax, cz_shift)
                                <= (cell_shift(lmax, height_shift,cell)),

                                dataobject.data)
        end

    elseif inverse == true
        if isamr
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                > (cell_shift(p.level, radius_shift,cell))  ||

                                get_height_cylinder(p.cz, p.level, cz_shift)
                                > (cell_shift(p.level, height_shift,cell)),

                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                > (cell_shift(lmax, radius_shift,cell))  ||

                                get_height_cylinder(p.cz, lmax, cz_shift)
                                > (cell_shift(lmax, height_shift,cell)),

                                dataobject.data)
        end

        ranges = dataobject.ranges

    end
        printtablememory(sub_data, verbose)

        hydrodata = HydroDataType()
        hydrodata.data = sub_data
        hydrodata.info = dataobject.info
        hydrodata.lmin = dataobject.lmin
        hydrodata.lmax = dataobject.lmax
        hydrodata.boxlen = dataobject.boxlen
        hydrodata.ranges = ranges
        hydrodata.selected_hydrovars = dataobject.selected_hydrovars
        hydrodata.used_descriptors = dataobject.used_descriptors
        hydrodata.smallr = dataobject.smallr
        hydrodata.smallc = dataobject.smallc
        hydrodata.scale = dataobject.scale
        return hydrodata

end

# -----------------------------------------------------------------------------
##### SPHERE #####-------------------------------------------------------------

function get_radius_sphere(x,y,z, level, cx_shift, cy_shift, cz_shift)
    xcenter = 2^level * cx_shift
    ycenter = 2^level * cy_shift
    zcenter = 2^level * cz_shift
    x = (x - xcenter)
    y = (y - ycenter)
    z = (z - zcenter)
    return sqrt( x^2 + y^2 + z^2)
end

function subregionsphere(dataobject::HydroDataType;
                            radius::Real=0.,
                            center::Array{<:Any,1}=[0., 0., 0.],
                            range_unit::Symbol=:standard,
                            cell::Bool=true,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    if radius == 0. || in(0., center)
        error("[Mera]: given radius or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    # convert given ranges and print overview on screen
    height = 0.
    ranges, cx_shift, cy_shift, cz_shift, radius_shift = prepranges(dataobject.info, center, radius, height, range_unit, verbose)

    if inverse == false
        if isamr
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)

                                <= cell_shift(p.level, radius_shift, cell),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)

                                <= cell_shift(lmax, radius_shift, cell),
                                dataobject.data)
        end
    elseif inverse == true
        if isamr
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)

                                > cell_shift(p.level, radius_shift, cell),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)

                                > cell_shift(lmax, radius_shift, cell),
                                dataobject.data)
        end

        ranges = dataobject.ranges

    end

    printtablememory(sub_data, verbose)

    hydrodata = HydroDataType()
    hydrodata.data = sub_data
    hydrodata.info = dataobject.info
    hydrodata.lmin = dataobject.lmin
    hydrodata.lmax = dataobject.lmax
    hydrodata.boxlen = dataobject.boxlen
    hydrodata.ranges = ranges
    hydrodata.selected_hydrovars = dataobject.selected_hydrovars
    hydrodata.used_descriptors = dataobject.used_descriptors
    hydrodata.smallr = dataobject.smallr
    hydrodata.smallc = dataobject.smallc
    hydrodata.scale = dataobject.scale
    return hydrodata

end


text:
@testset "search tests" begin

    @testset "basic usage tests" begin
        for GraphT in [HeapKNNGraph, LockHeapKNNGraph]
            data = [rand(5) for _ in 1:20]
            queries = [rand(5) for _ in 1:5]
            graph = nndescent(GraphT, data, 4, Euclidean(); max_iters=1)
            # vector of vector usage
            @inferred search(graph, queries, 3)

            data = rand(5, 20)
            queries = rand(5, 5)
            # matrix usage
            @inferred search(graph, queries, 3)
        end
        @test true
    end

    Random.seed!(0)
    data = [[0., 0., 0.],
            [0., 0., 1.],
            [0., 1., 0.],
            [0., 1., 1.],
            [1., 0., 0.],
            [1., 0., 1.],
            [1., 1., 0.],
            [1., 1., 1.]]
    queries = [[0., 0., 0.4],
               [0., 1., 0.4]]
    true_inds = [1 3;
                 2 4]
    true_dists = [.4 .4;
                  .6 .6]

    for GraphT in [HeapKNNGraph, LockHeapKNNGraph]
        graph = nndescent(GraphT, data, 3, Euclidean())
        inds, dists = search(graph, queries, 2)

        @test all(inds .== true_inds)
        @test all(dists .== true_dists)
    end

end


text:
const cube_edges = (0x00, 0x01, 0x00, 0x02, 0x00, 0x04, 0x01, 0x03, 0x01, 0x05, 0x02, 0x03,
                    0x02, 0x06, 0x03, 0x07, 0x04, 0x05, 0x04, 0x06, 0x05, 0x07, 0x06, 0x07)

const sn_edge_table = (0x0007, 0x0019, 0x001e, 0x0062, 0x0065, 0x007b, 0x007c, 0x00a8,
                       0x00af, 0x00b1, 0x00b6, 0x00ca, 0x00cd, 0x00d3, 0x00d4, 0x0304,
                       0x0303, 0x031d, 0x031a, 0x0366, 0x0361, 0x037f, 0x0378, 0x03ac,
                       0x03ab, 0x03b5, 0x03b2, 0x03ce, 0x03c9, 0x03d7, 0x03d0, 0x0510,
                       0x0517, 0x0509, 0x050e, 0x0572, 0x0575, 0x056b, 0x056c, 0x05b8,
                       0x05bf, 0x05a1, 0x05a6, 0x05da, 0x05dd, 0x05c3, 0x05c4, 0x0614,
                       0x0613, 0x060d, 0x060a, 0x0676, 0x0671, 0x066f, 0x0668, 0x06bc,
                       0x06bb, 0x06a5, 0x06a2, 0x06de, 0x06d9, 0x06c7, 0x06c0, 0x0a40,
                       0x0a47, 0x0a59, 0x0a5e, 0x0a22, 0x0a25, 0x0a3b, 0x0a3c, 0x0ae8,
                       0x0aef, 0x0af1, 0x0af6, 0x0a8a, 0x0a8d, 0x0a93, 0x0a94, 0x0944,
                       0x0943, 0x095d, 0x095a, 0x0926, 0x0921, 0x093f, 0x0938, 0x09ec,
                       0x09eb, 0x09f5, 0x09f2, 0x098e, 0x0989, 0x0997, 0x0990, 0x0f50,
                       0x0f57, 0x0f49, 0x0f4e, 0x0f32, 0x0f35, 0x0f2b, 0x0f2c, 0x0ff8,
                       0x0fff, 0x0fe1, 0x0fe6, 0x0f9a, 0x0f9d, 0x0f83, 0x0f84, 0x0c54,
                       0x0c53, 0x0c4d, 0x0c4a, 0x0c36, 0x0c31, 0x0c2f, 0x0c28, 0x0cfc,
                       0x0cfb, 0x0ce5, 0x0ce2, 0x0c9e, 0x0c99, 0x0c87, 0x0c80, 0x0c80,
                       0x0c87, 0x0c99, 0x0c9e, 0x0ce2, 0x0ce5, 0x0cfb, 0x0cfc, 0x0c28,
                       0x0c2f, 0x0c31, 0x0c36, 0x0c4a, 0x0c4d, 0x0c53, 0x0c54, 0x0f84,
                       0x0f83, 0x0f9d, 0x0f9a, 0x0fe6, 0x0fe1, 0x0fff, 0x0ff8, 0x0f2c,
                       0x0f2b, 0x0f35, 0x0f32, 0x0f4e, 0x0f49, 0x0f57, 0x0f50, 0x0990,
                       0x0997, 0x0989, 0x098e, 0x09f2, 0x09f5, 0x09eb, 0x09ec, 0x0938,
                       0x093f, 0x0921, 0x0926, 0x095a, 0x095d, 0x0943, 0x0944, 0x0a94,
                       0x0a93, 0x0a8d, 0x0a8a, 0x0af6, 0x0af1, 0x0aef, 0x0ae8, 0x0a3c,
                       0x0a3b, 0x0a25, 0x0a22, 0x0a5e, 0x0a59, 0x0a47, 0x0a40, 0x06c0,
                       0x06c7, 0x06d9, 0x06de, 0x06a2, 0x06a5, 0x06bb, 0x06bc, 0x0668,
                       0x066f, 0x0671, 0x0676, 0x060a, 0x060d, 0x0613, 0x0614, 0x05c4,
                       0x05c3, 0x05dd, 0x05da, 0x05a6, 0x05a1, 0x05bf, 0x05b8, 0x056c,
                       0x056b, 0x0575, 0x0572, 0x050e, 0x0509, 0x0517, 0x0510, 0x03d0,
                       0x03d7, 0x03c9, 0x03ce, 0x03b2, 0x03b5, 0x03ab, 0x03ac, 0x0378,
                       0x037f, 0x0361, 0x0366, 0x031a, 0x031d, 0x0303, 0x0304, 0x00d4,
                       0x00d3, 0x00cd, 0x00ca, 0x00b6, 0x00b1, 0x00af, 0x00a8, 0x007c,
                       0x007b, 0x0065, 0x0062, 0x001e, 0x0019, 0x0007)


text:
module MeshCatMechanisms

export MechanismVisualizer,
       animate,
       MeshCatSink,
       setelement!,
       manipulate!

# Re-export from MeshCat.jl
export render,
       Triad,
       setanimation!

# Re-export from MechanismGeometries.jl
export VisualElement,
       Skeleton,
       URDFVisuals,
       visual_elements

# Re-export from RigidBodyDynamics.jl
export set_configuration!,
       Point3D

using MeshCat
using CoordinateTransformations
using RigidBodyDynamics
const rbd = RigidBodyDynamics
using RigidBodyDynamics.Graphs: ancestors, edge_to_parent, source, vertices, root
using Interpolations: interpolate, Gridded, Linear
using LoopThrottle: @throttle
using MechanismGeometries: visual_elements, VisualElement, Skeleton, URDFVisuals, AbstractGeometrySource, MeshFile
using GeometryTypes: HyperSphere, Point

include("visualizer.jl")
include("animate.jl")
include("ode_callback.jl")
include("manipulate.jl")
using .Manipulate

end # module


text:
function calc_pi(n::Integer)
   x = rand(n)
   y = rand(n)
   num_in_unit_circle = sum(x.*x.+y.*y .< 1)
   num_in_unit_circle * 4/n
end

n = 10_000;
calc_pi(1)  # Force compile
@time pi_estimate = calc_pi(n)
println("# After ", n, " itterations, estimated pi to be...")
println(pi_estimate)


text:
using AutoViz
using Distributions
using Interact
using AutomotiveDrivingModels
using Dates
using Parameters

using NNMPC

# include("/home/sbae/NNMPC.jl/behaviors/baffling_drivers.jl")
# include("/home/sbae/NNMPC.jl/behaviors/baffling_lane_changer.jl")
# include("/home/sbae/NNMPC.jl/behaviors/baffling_longitudinal_tracker.jl")
# include("/home/sbae/NNMPC.jl/behaviors/nn_mpc_driver.jl")

@time begin
roadway = gen_stadium_roadway(3)
timestep = 0.2
scene = Scene()

# ==========================
# Generate Vehicles
# ==========================
n_vehs_top = 15; origin_top = 3.0;
n_vehs_middle = 16; origin_middle = 0.0;
n_vehs_bottom = 15; origin_bottom = 3.0;
n_vehs_tot = n_vehs_top + n_vehs_middle + n_vehs_bottom
dist_betw_cars = 6.0
global num_vehs = 0
# top lane
for i in 1: n_vehs_top
    global num_vehs
    num_vehs += 1
    if i == 1
        curr_speed = 5.0
    else
        curr_speed = 10.0
    end
    push!(scene, Vehicle(VehicleState(VecSE2(origin_top+dist_betw_cars*(i-1),0.0,0.0),roadway,curr_speed),
                         VehicleDef(),
                         num_vehs));
end
# middle lane
for i in 1: n_vehs_middle
    global num_vehs
    num_vehs += 1
    if i == 1
        curr_speed = 5.0
    elseif i == 2
        curr_speed = 10.0
    else
        curr_speed = 10.0
    end
    push!(scene, Vehicle(VehicleState(VecSE2(origin_middle+dist_betw_cars*(i-1),-DEFAULT_LANE_WIDTH,0.0),roadway,curr_speed),
                 VehicleDef(),
                 num_vehs));
end
# bottom lane
for i in 1: n_vehs_bottom
    global num_vehs
    num_vehs += 1
    if i == 1
        curr_speed = 5.0
    else
        curr_speed = 10.0
    end
    push!(scene, Vehicle(VehicleState(VecSE2(origin_bottom+dist_betw_cars*(i-1),-2*DEFAULT_LANE_WIDTH,0.0),roadway,curr_speed),
                 VehicleDef(),
                 num_vehs));
end

# ==========================
# Assign Driver models
# ==========================
nticks = 10
pred_model_type = "perfect" # {sgan, perfect}
ind_ego = n_vehs_top+1
models = Dict{Int, DriverModel}()
for j in 1:num_vehs
    if j == ind_ego
        models[j] = NNMPCdriver(timestep,
                                        n_ticks = nticks+1,
                                        N_sim = 200.0,
                                        T = 1.0, # receding horizon
                                        λ_div = 10000.0,
                                        λ_v = 1000.0,
                                        λ_δ = 500.0,
                                        λ_Δδ = 10.0,
                                        λ_a = 100.0,
                                        δ_max= 0.15,
                                        δ_min = -0.15,
                                        a_min = -9.0,
                                        Δa_max = 10.0,
                                        Δδ_max = 0.2,
                                        width = 2.0,
                                        height = 4.2,
                                        T_obs = 1.6,
                                        thred_safety = 0.3,
                                        pred_model_type = pred_model_type,
                                        isDebugMode = false,
                                        lane_change_action = LaneChangeChoice(DIR_LEFT),
                                        veh_shape_type = "circle"
                                        )
    else
        models[j] = IntelligentDriverModel()
    end
    set_desired_speed!(models[j], 10.0)
    # if j == 4 || j == 8 || j == 12
    #     set_desired_speed!(models[j], 3.0)
    # end
end
set_other_models!(models[ind_ego],models)
set_desired_speed!(models[ind_ego], 10.0)

# ==========================
# Set vehicle colors
# ==========================
# Set colors
saturation=0.85
value=0.85
car_colors = Dict{Int,Colorant}()
n = length(scene)
for (i,veh) in enumerate(scene)
    car_colors[veh.id] = convert(RGB, HSV(0, saturation, value))
end
car_colors[ind_ego] = convert(RGB, HSV(100, saturation, value))
# car_colors = get_pastel_car_colors(scene)
cam = FitToContentCamera()

# ==========================
# Run simulation
# ==========================
rec = SceneRecord(nticks+1, timestep)
simulate!(rec, scene, roadway, models, nticks)

#
# # ===========================
# # Store the variables
# # ===========================
# using JLD
# JLD.save("../sim_records/$(Dates.DateTime(Dates.now()))_N$(trunc(Int,models[ind_ego].N_sim))_tick$(nticks)_vehs$(n_vehs_tot)_Nhorz$(trunc(Int,round(models[ind_ego].T/timestep)))_pred$(pred_model_type).jld",
#     "rec", rec, "roadway", roadway, "cam", cam,
#     "colors", car_colors, "ind", ind_ego)
#
# ===========================
# Rendering
# ===========================
render(rec[0], roadway, cam=cam, car_colors=car_colors)
@manipulate for frame_index in 1:nframes(rec)
    render(rec[frame_index-nframes(rec)],roadway,cam=cam,car_colors=car_colors)
end

end


text:
"""
    xlogx(x)

Return `x * log(x)` for `x ≥ 0`, handling `x == 0` by taking the limit from above, to get zero.
"""
function xlogx(x)
  result = x * log(x)
  ifelse(iszero(x), zero(result), result)
end

"""
    xlogy(x, y)

Return `x * log(y)` for `y > 0`, and zero when `x == 0`.
"""
function xlogy(x, y)
  result = x * log(y)
  ifelse(iszero(x), zero(result), result)
end

@adjoint function broadcasted(::typeof(xlogy), x::Zygote.Numeric, y::Zygote.Numeric)
  res = xlogy.(x, y)
  res, Δ -> (nothing, Zygote.unbroadcast(x, xlogy.(Δ, y)), Zygote.unbroadcast(y, Δ .* x ./ y))
end

# This can be made an error in Flux v0.13, for now just a warning
function _check_sizes(ŷ::AbstractArray, y::AbstractArray)
  for d in 1:max(ndims(ŷ), ndims(y))
    if size(ŷ,d) != size(y,d)
      @warn "Size mismatch in loss function! In future this will be an error. In Flux <= 0.12 broadcasting accepts this, but may not give sensible results" summary(ŷ) summary(y) maxlog=3 _id=hash(size(y))
    end
  end
end
_check_sizes(ŷ, y) = nothing  # pass-through, for constant label e.g. y = 1

Zygote.@nograd _check_sizes


text:
using DevilSpawn
using Test

@testset "DevilSpawn.jl" begin
    spawned = [fetch(@spawn_background Threads.threadid()) for _ in 1:100]
    if Threads.nthreads() > 1
        @test !in(1, spawned)
    else
        @test all(==(1), spawned)
    end
end