text:
################################################# FILE DESCRIPTION ##########################################################

# This file contains graph compatibility methods for SparseMatrixCSC.

################################################# HELPERS ###################################################################

""" Sort the index entries """
function reorder!(x::SparseMatrixCSC)
   x.nzval[:] = 1 : nnz(x)
   return x
end

################################################# ACCESSORS #################################################################

""" The number of vertices in the adjacency matrix """
nv(x::SparseMatrixCSC{Int}) = x.m

""" The number of edges in the adjacency matrix """
ne(x::SparseMatrixCSC{Int}) = nnz(x)

################################################# GENERATION ################################################################

""" Build an adjacency matrix from an edge iterator """
function SparseMatrixCSC(nv::Int, eit::EdgeIter)
   sparse(eit.us, eit.vs, 1, nv, nv, min)
end

"""
Spawn a random sparse matrix, sort indices and remove self loops

The number of edges in the sparse matrix may not equal the input ne,
and is more likely to be an approximate
"""
function randindxs(nv::Int, ne::Int)
   sv = sprand(Int, nv, nv, ne/(nv*(nv-1)))
   sv = sv - spdiagm(diag(sv), 0)
   reorder!(sv)
   return sv
end

"""
Spawn a random sparse matrix denoting a complete graph without
self loops. The returned matrix resembles a non-sparse matrix,
so it'd be unwise to use this for a large number of vertices
"""
function completeindxs(nv::Int)
   sv = sparse(spzeros(Int, nv, nv) .+ 1 - speye(Int, nv, nv))
   reorder!(sv)
   return sv
end

################################################# PAIR GETINDEX #############################################################

""" Retrieve an edges' index in the ddge dataframe """
Base.getindex(x::SparseMatrixCSC{Int}, e::EdgeID) = x[e.second, e.first]

""" Retrieve the edge dataframe indices for a list of edges """
function Base.getindex(x::SparseMatrixCSC{Int}, es::EdgeList)
   n = length(es)
   vals = sizehint!(Vector{Int}(), n)
   for e in es
      push!(vals, getindex(x, e))
   end
   return vals
end

""" Retrieve the edge dataframe indices for all edges in an edge iterator """
function Base.getindex(x::SparseMatrixCSC{Int}, eit::EdgeIter)
   n = length(eit)
   vals = sizehint!(Vector{Int}(), n)
   for i in 1 : n
      @inbounds u = eit.us[i]
      @inbounds v = eit.vs[i]
      push!(vals, getindex(x, v, u))
   end
   return vals
end

################################################# PAIR SETINDEX #############################################################

""" Change an edge's index in the edge dataframe """
function Base.setindex!(x::SparseMatrixCSC{Int}, val::Int, e::EdgeID)
   x[e.second, e.first] = val
end

""" Change the edge dataframe indices for a list of edges """
function Base.setindex!(x::SparseMatrixCSC{Int}, val, es::EdgeList)
   for e in es
      x[e] = val
   end
end

""" Change the edge dataframe indices for a list of edges """
function Base.setindex!(x::SparseMatrixCSC{Int}, vals::AbstractVector{Int}, es::EdgeList)
   for i in eachindex(vals, es)
      @inbounds e = es[i]
      @inbounds val = vals[i]
      x[e] = val
   end
end

""" Change the edge dataframe indices for all edges in an iterator """
function Base.setindex!(x::SparseMatrixCSC{Int}, val::Int, eit::EdgeIter)
   n = length(eit)
   for i in 1 : n
      @inbounds u = eit.us[i]
      @inbounds v = eit.vs[i]
      x[v,u] = val
   end
end

""" Change the edge dataframe indices for all edges in an iterator """
function Base.setindex!(x::SparseMatrixCSC{Int}, vals::AbstractVector{Int}, eit::EdgeIter)
   n = length(eit)
   for i in 1 : n
      @inbounds u = eit.us[i]
      @inbounds v = eit.vs[i]
      @inbounds val = vals[i]
      x[v,u] = val
   end
end
################################################# ADJACENCY #################################################################

"""
Retrieve a list of vertices connected to vertex v.

This method spwans a new array, so is slow and malloc prone.
"""
function fadj(x::SparseMatrixCSC{Int}, v::VertexID)
   x.rowval[nzrange(x, v)]
end

"""
Retrieve a list of vertices connect to vertex v.

This method copies the adjacencies onto the input array,
and is comparitively faster, and causes no mallocs.
"""
# TODO: Replace this functionaly with fast subarrays, once that gets out.
function fadj!(x::SparseMatrixCSC{Int,Int}, v::VertexID, adj::Vector{Int})
   @inbounds p1 = x.colptr[v]
   @inbounds p2 = x.colptr[v+1]
   resize!(adj, p2 - p1)
   copy!(adj, 1, x.rowval, p1, p2 - p1)
end

""" Compute the outdegree of a vertex """
function outdegree(x::SparseMatrixCSC{Int}, v::VertexID)
   @inbounds p1 = x.colptr[v]
   @inbounds p2 = x.colptr[v+1]
   return p2 - p1
end

""" Compute outdegrees for a list of vertices """
function outdegree(x::SparseMatrixCSC{Int}, vs::VertexList)
   [outdegree(x, v) for v in vs]
end

""" Compute outdegrees for all vertices in the graph """
function outdegree(x::SparseMatrixCSC{Int})
   degs = Vector{Int}(nv(x))
   for v in 1 : nv(x)
      degs[v] = x.colptr[v+1] - x.colptr[v]
   end
   return degs
end

"""
Compute the indegree of a vertex. This method is slow
since reverse adjacencies are not stored
"""
function indegree(x::SparseMatrixCSC{Int}, v::VertexID)
   count = 0
   for i in x.rowval
      if i == v
         count += 1
      end
   end
   return count
end

""" Compute the indegrees for a list of vertices. Note that the list may not be unique! """
function indegree(x::SparseMatrixCSC{Int}, vs::VertexList)
   degs = zeros(Int, nv(x))
   for i in x.rowval
      degs[i] += 1
   end
   return degs[vs]
end

""" Compute the indegrees for all vertices in the graph """
function indegree(x::SparseMatrixCSC{Int})
   degs = zeros(Int, nv(x))
   for i in x.rowval
      degs[i] += 1
   end
   return degs
end

################################################# ADDVERTEX ###############################################################

"""
Grow the input sparsematrix to hold another vertex.
Returns a new matrix, since the SparseMatrixCSC type is immutable.
The input matrix is reused in the construction.
"""
function addvertex!(x::SparseMatrixCSC{Int,Int})
   SparseMatrixCSC{Int,Int}(nv(x)+1, nv(x)+1, push!(x.colptr, x.colptr[end]), x.rowval, x.nzval)
end

################################################# ADDEDGE #################################################################

""" Insert an edge into the sparsematrix """
function addedge!(x::SparseMatrixCSC{Int}, e::EdgeID, erow::Int)
   x[e] = erow
end

"""  Insert a list of edges into the sparsematrix """
function addedge!(x::SparseMatrixCSC{Int}, es::EdgeList, erows::Vector{Int})
   x[es] = erows
end

################################################# RMVERTEX ################################################################

"""
Remove a vertex(s) from the sparsematrx. This method
first computes a list of rows to be removed from the edge dataframe.
It then computes and returns a subset of the input sparsematrix, and reorders
edge indices.
"""
function rmvertex!(x::SparseMatrixCSC{Int,Int}, vs)
   # Reorder entries in the index table
   reorder!(x)

   # Check which entries in the edge table have to be removed
   erows = vcat(x[EdgeIter(x, vs)], x[EdgeIter(x, :, vs)])

   # Delete entries from index table
   vlist = collect(1 : x.m)
   deleteat!(vlist, vs)
   x = x[vlist,vlist]
   reorder!(x)

   return(x, erows)
end

################################################# RMVERTEX ################################################################

"""
Remove an edge from the sparsematrix. This method returns the edge dataframe
index of the removed edge.
"""
function rmedge!(x::SparseMatrixCSC{Int}, e::EdgeID)
   erow = x[e]
   x[e] = 0
   dropzeros!(x)
   x.nzval[:] = 1 : nnz(x)
   return erow
end

"""
Remove edges from the sparsematrix. This method returns the edge dataframe
index of the removed edges.
"""
function rmedge!(x::SparseMatrixCSC{Int}, es::EdgeList)
   erows = x[es]
   x[es] = 0
   dropzeros!(x)
   x.nzval[:] = 1 : nnz(x)
   return erows
end

################################################# SUBGRAPH ################################################################

"""
Compute a subgraph containing the list of input vertices. This method
uses vector getindex and is very fast. It computes the list of rows in
the edge dataframed to be copied into the subgraph, and reorders entries in the
new adjacency matrix
"""
function subgraph(x::SparseMatrixCSC{Int}, vs::VertexList)
   sv = x[vs,vs]
   erows = copy(nonzeros(sv))
   reorder!(sv)
   return(sv, erows)
end

"""
Compute a subgraph from a list of edges. This method builds a new
adjacency matrix from the iterator, reorders the indices and returns a
list of edge dataframe rows to be copied into the subgraph.
"""
function subgraph(x::SparseMatrixCSC{Int}, eit::EdgeIter)
   nv = size(x, 1)
   erows = x[eit]
   sv = sparse(eit.vs, eit.us, collect(1 : length(eit)), nv, nv)
   return(sv, erows)
end

subgraph(x::SparseMatrixCSC{Int}, es::EdgeList) = subgraph(x, EdgeIter(es))

"""
Compute a subgraph from a list of vertices and a list of edges. This method
first computes a subgraph comprising of the input edges, and then proceeds
to preserve only the input list of vertices. Reorders indices and returns a
list of edge table rows to be copied into the subgraph.
"""
function subgraph(x::SparseMatrixCSC{Int}, vs::VertexList, eit::EdgeIter)
   nv = size(x, 1)
   sv = sparse(eit.vs, eit.us, x[eit], nv, nv)[vs,vs]
   erows = sort(nonzeros(sv))
   reorder!(sv)
   return(sv, erows)
end

subgraph(x::SparseMatrixCSC{Int}, vs::VertexList, es::EdgeList) = subgraph(x, vs, EdgeIter(es))


text:
open("data.txt","w") do fp
    for i in 1:10
        write(fp, "$i\n")
    end
end

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module SymbolServer

import Sockets
pipename = length(ARGS) > 1 ? ARGS[2] : nothing
conn = pipename !== nothing ? Sockets.connect(pipename) : nothing

start_time = time_ns()

# Try to lower the priority of this process so that it doesn't block the
# user system.
@static if Sys.iswindows()
    # Get process handle
    p_handle = ccall(:GetCurrentProcess, stdcall, Ptr{Cvoid}, ())

    # Set BELOW_NORMAL_PRIORITY_CLASS
    ret = ccall(:SetPriorityClass, stdcall, Cint, (Ptr{Cvoid}, Culong), p_handle, 0x00004000)
    ret != 1 && @warn "Something went wrong when setting BELOW_NORMAL_PRIORITY_CLASS."
else
    ret = ccall(:nice, Cint, (Cint,), 1)
    # We don't check the return value because it doesn't really matter
end

module LoadingBay
end

!in("@stdlib", LOAD_PATH) && push!(LOAD_PATH, "@stdlib") # Make sure we can load stdlibs
using Pkg, SHA
using Base: UUID

include("faketypes.jl")
include("symbols.jl")
include("utils.jl")
include("serialize.jl")
using .CacheStore

store_path = length(ARGS) > 0 ? ARGS[1] : abspath(joinpath(@__DIR__, "..", "store"))

ctx = try
    Pkg.Types.Context()
catch err
    @info "Package environment can't be read."
    exit()
end
# Add some methods to check whether a package is part of the standard library and so
# won't need recaching.
if isdefined(Pkg.Types, :is_stdlib)
    is_stdlib(uuid::UUID) = Pkg.Types.is_stdlib(uuid)
else
    is_stdlib(uuid::UUID) = uuid in keys(ctx.stdlibs)
end

server = Server(store_path, ctx, Dict{UUID,Package}())

written_caches = String[] # List of caches that have already been written
toplevel_pkgs = deps(project(ctx)) # First get a list of all package UUIds that we want to cache
packages_to_load = []
# Next make sure the cache is up-to-date for all of these

for (pk_name, uuid) in toplevel_pkgs
    isinmanifest(ctx, uuid) || continue
    pe = frommanifest(manifest(ctx), uuid)
    cache_path = joinpath(server.storedir, SymbolServer.get_cache_path(manifest(ctx), uuid)...)

    if isfile(cache_path)
        if is_package_deved(manifest(ctx), uuid)
            try
                cached_version = open(cache_path) do io
                    CacheStore.read(io)
                end
                if sha_pkg(frommanifest(manifest(ctx), uuid)) != cached_version.sha
                    @info "Outdated sha, will recache package $pk_name ($uuid)"
                    push!(packages_to_load, uuid)
                else
                    @info "Package $pk_name ($uuid) is cached."
                end
            catch err
                @info "Couldn't load $pk_name ($uuid) from file, will recache."
            end
        else
            @info "Package $pk_name ($uuid) is cached."
        end
    else
        @info "Will cache package $pk_name ($uuid)"
        push!(packages_to_load, uuid)
    end
end

# Load all packages together
for uuid in packages_to_load
    load_package(ctx, uuid, conn, LoadingBay)
end

# Create image of whole package env. This creates the module structure only.
env_symbols = getenvtree()

# Populate the above with symbols, skipping modules that don't need caching.
# symbols (env_symbols)
visited = Base.IdSet{Module}([Base, Core]) # don't need to cache these each time...
for (pid, m) in Base.loaded_modules
    if pid.uuid !== nothing && is_stdlib(pid.uuid) &&
        isinmanifest(ctx, pid.uuid) &&
        isfile(joinpath(server.storedir, SymbolServer.get_cache_path(manifest(ctx), pid.uuid)...))
        push!(visited, m)
        delete!(env_symbols, Symbol(pid.name))
    end
end

symbols(env_symbols, nothing, getallns(), visited)

# Wrap the `ModuleStore`s as `Package`s.
for (pkg_name, cache) in env_symbols
    !isinmanifest(ctx, String(pkg_name)) && continue
    uuid = packageuuid(ctx, String(pkg_name))
    pe = frommanifest(ctx, uuid)
    server.depot[uuid] = Package(String(pkg_name), cache, uuid, sha_pkg(pe))
end

write_depot(server, server.context, written_caches)

@info "Symbol server indexing took $((time_ns() - start_time) / 1e9) seconds."

end


text:
using SolidStateDetectors
using Plots
using SpecialFunctions

struct CustomChargeDensity{T} <: SolidStateDetectors.AbstractChargeDensity{T}
    A::T
    S::T
    PN::T
    NP::T # N+
    PP::T # P+
end
function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::CartesianPoint{T}) where {T}
    # just a dummy model for a pn-junction with overdoped n+ and p+ layer on both ends
    cdm.A * erf(cdm.S * (pt[1] - cdm.PN)) +
    2 * cdm.A * erf(10 * cdm.S * (pt[1] - cdm.NP)) +
    2 * cdm.A * erf(10 * cdm.S * (pt[1] - cdm.PP))
end
function SolidStateDetectors.get_charge_density(cdm::CustomChargeDensity{T}, pt::CylindricalPoint{T}) where {T}
    SolidStateDetectors.get_charge_density(cdm, CartesianPoint(pt))
end

begin
    T = Float32
    sim = Simulation{T}(SSD_examples[:InfiniteParallelPlateCapacitor])

    ccdm = CustomChargeDensity{Float32}( 1*10^17, -10^4, 0, 0.00035, -0.00035 ) # meter not mm
    sim.detector.semiconductors[1].charge_density_model = ccdm

    xs = -5e-4:1e-6:5e-4
    charge_vec = similar(xs)
    for ix in eachindex(xs)
        pt = CartesianPoint{T}(xs[ix], 0, 0)
        charge_vec[ix] = pt in sim.detector.semiconductors[1] ?
            SolidStateDetectors.get_charge_density(sim.detector.semiconductors[1], pt) :
            0
    end

    calculate_electric_potential!(sim,
        max_refinements = 4,
        init_grid_spacing = T.( (1e-4, 1e-3, 1e-3) ),
        convergence_limit = 1e-3,
        depletion_handling=true
    )

    plot(
        plot(xs, charge_vec, ylabel = "impurity density [1/m^3]", xlabel = "x [m]"),
        plot(sim.ρ, aspect_ratio = :none, y = 0),
        plot(sim.electric_potential, aspect_ratio = :none, y = 0),
        plot(sim.point_types, aspect_ratio = :none, y = 0),
        layout = (4,1), size = (800, 1200)
    )
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
using Images
using PixelArt

img = load("img.jpg")
img_pixel = pixel(img)

save("img_pixel.jpg", img_pixel)


text:
"""
General evolution for spin s field about Kerr
"""

module Evolution

include("Fields.jl")
include("Radial.jl")
include("Sphere.jl")
include("BackgroundNP.jl")
include("GHP.jl")

using  .Fields: Field
import .Radial
import .Sphere
import .BackgroundNP: NP_0
import .GHP: GHP_ops
#import .MetricReconstruction

const half  = 1.0/2.0
const third = 1.0/3.0
const sixth = 1.0/6.0

export Evo_lin_f, Initialize_Evo_lin_f, Evolve_lin_f!

struct Evo_lin_f
   A_pp::Array{Float64,2}
   A_pq::Array{Float64,2}
   B_pp::Array{ComplexF64,2}
   B_pq::Array{ComplexF64,2}
   B_pf::Array{ComplexF64,2}

   S_lapl::Array{Float64,2}
   S_fltr::Array{Float64,2}

   mv::Int64

   function Evo_lin_f(;
         Rvals::Vector{Float64},
         Cvals::Vector{Float64},
         Svals::Vector{Float64},
         bhm::Float64,
         bhs::Float64,
         mv::Int64,
         cl::Float64,
         spin::Int64
      )
      nx = length(Rvals)
      ny = length(Cvals)

      A_pp = zeros(Float64,nx,ny)
      A_pq = zeros(Float64,nx,ny)

      B_pp = zeros(ComplexF64,nx,ny)
      B_pq = zeros(ComplexF64,nx,ny)
      B_pf = zeros(ComplexF64,nx,ny)

      S_lapl = zeros(Float64,ny,ny)
      S_fltr = zeros(Float64,ny,ny)

      pre = zeros(Float64,nx,ny)

      for j=1:ny
         sy = Svals[j]
         for i=1:nx
            R = Rvals[i]

            pre[i,j] = 1.0 / (
               8*bhm*(2*(cl^2)*bhm - (bhs^2)*R)*((cl^2) + 2*bhm*R) / (cl^4)
            -  (bhs*sy)^2
            )
         end
      end

      S_lapl = Sphere.swal_laplacian_matrix(ny,spin,mv)
      S_fltr = Sphere.swal_filter_matrix(   ny,spin,mv)

      for i=1:ny
         for j=1:ny
            S_lapl[i,j] = pre[i,j]*S_lapl[i,j]
            S_fltr[i,j] =          S_fltr[i,j]
         end
      end

      for j=1:ny
         cy = Cvals[j]
         for i=1:nx
            R = Rvals[i]

            A_pp[i,j] = (2 / (cl^4)) * (
               (cl^6)
            +  (cl^2)*(bhs^2 - 8*(bhm^2))*(R^2)
            +  4*(bhs^2)*bhm*(R^3)
            )

            A_pq[i,j] = ((R^2) / (cl^4)) * (
               (cl^4)
            -  2*(cl^2)*bhm*R
            +  (bhs^2)*(R^2)
            )

            B_pp[i,j] = - (
               2*im*bhs*mv*((cl^2) + 4*bhm*R) / (cl^2)
            -  2*(bhs^2)*(cl^2 + 6*bhm*R)*R / (cl^4)
            -  4*bhm*spin
            +  8*(bhm^2)*(2 + spin)*R / (cl^2)
            +  2*im*bhs*spin*cy
            )

            B_pq[i,j] = (2 * R / (cl^4)) * (
               2*(bhs^2)*(R^2)
            +  (1 + spin)*(cl^4)
            -  (im*bhs*mv + (3 + spin)*bhm) * ((cl^2) * R)
            )

            B_pf[i,j] = - (2 * R / (cl^4)) * (
               im*bhs*(cl^2)*mv
            -  (bhs^2)*R
            +  (cl^2)*bhm*(1 + spin)
            )

            A_pp[i,j] *= pre[i,j]
            A_pq[i,j] *= pre[i,j]

            B_pp[i,j] *= pre[i,j]
            B_pq[i,j] *= pre[i,j]
            B_pf[i,j] *= pre[i,j]
         end
      end

      return new(A_pp,A_pq,B_pp,B_pq,B_pf,S_lapl,S_fltr,mv)
   end
end

function Initialize_Evo_lin_f(;
      Rvals::Vector{Float64},
      Cvals::Vector{Float64},
      Svals::Vector{Float64},
      Mvals::Vector{Int64},
      bhm::Float64,
      bhs::Float64,
      cl::Float64,
      spin::Int64)
   return Dict([
      (mv,Evo_lin_f(Rvals=Rvals,Cvals=Cvals,Svals=Svals,bhm=bhm,bhs=bhs,mv=mv,cl=cl,spin=spin))
      for mv in Mvals
     ])
end

"""
set_kp
"""
function set_kp(
      kp::Array{ComplexF64,2},
      f_rd1::Array{ComplexF64,2},
      f_rd2::Array{ComplexF64,2},
      sph_lap::Array{ComplexF64,2},
      p_rd1::Array{ComplexF64,2},
      f::Array{ComplexF64,2},
      p::Array{ComplexF64,2},
      A_pp::Array{Float64,2},
      A_pq::Array{Float64,2},
      B_pp::Array{ComplexF64,2},
      B_pq::Array{ComplexF64,2},
      B_pf::Array{ComplexF64,2}
   )
   nx, ny = size(kp)
   for j=1:ny
      for i=1:nx
         kp[i,j] = (
               A_pp[i,j] * p_rd1[i,j]
             + A_pq[i,j] * f_rd2[i,j]
             + B_pp[i,j] * p[i,j]
             + B_pq[i,j] * f_rd1[i,j]
             + B_pf[i,j] * f[i,j]

             + sph_lap[i,j]
      )
      end
   end
   return nothing
end

"""
Fourth order Runge-Kutta evolution of linear field
"""
function Evolve_lin_f!(
      lin_f,
      lin_p,
      Evo::Evo_lin_f,
      dr::Float64,
      dt::Float64
   )::Nothing
   @assert Evo.mv == lin_f.mv
   @assert Evo.mv == lin_p.mv

   ## simplify names for convenience

   f_n       = lin_f.n
   f_tmp     = lin_f.tmp
   f_np1     = lin_f.np1
   f_k       = lin_f.k
   f_rd1     = lin_f.rad_d1
   f_rd2     = lin_f.rad_d2
   f_sph_lap = lin_f.sph_lap

   p_n       = lin_p.n
   p_tmp     = lin_p.tmp
   p_np1     = lin_p.np1
   p_k       = lin_p.k
   p_rd1     = lin_p.rad_d1

   A_pp = Evo.A_pp
   A_pq = Evo.A_pq
   B_pp = Evo.B_pp
   B_pq = Evo.B_pq
   B_pf = Evo.B_pf

   laplM = Evo.S_lapl
   fltrM = Evo.S_fltr

   nx, ny = size(f_n)

   ## step 1
   Radial.set_d1!(f_rd1, f_n, dr)
   Radial.set_d1!(p_rd1, p_n, dr)
   Radial.set_d2!(f_rd2, f_n, dr)

   Sphere.angular_matrix_mult!(f_sph_lap,f_n,laplM)

   set_kp(p_k,
      f_rd1, f_rd2, f_sph_lap, p_rd1, f_n, p_n,
      A_pp, A_pq,
      B_pp, B_pq, B_pf
   )
   for j=1:ny
      for i=1:nx
         f_k[i,j] = p_n[i,j]

         f_tmp[i,j] = f_n[i,j] + half*dt*f_k[i,j]
         p_tmp[i,j] = p_n[i,j] + half*dt*p_k[i,j]

         f_np1[i,j] = f_n[i,j] + sixth*dt*f_k[i,j]
         p_np1[i,j] = p_n[i,j] + sixth*dt*p_k[i,j]
      end
   end
   ## step 2
   Radial.set_d1!(f_rd1, f_tmp, dr)
   Radial.set_d1!(p_rd1, p_tmp, dr)
   Radial.set_d2!(f_rd2, f_tmp, dr)

   Sphere.angular_matrix_mult!(f_sph_lap,f_tmp,laplM)

   set_kp(p_k,
      f_rd1, f_rd2, f_sph_lap, p_rd1, f_tmp, p_tmp,
      A_pp, A_pq,
      B_pp, B_pq, B_pf
   )
   for j=1:ny
      for i=1:nx
         f_k[i,j] = p_tmp[i,j]

         f_tmp[i,j] = f_n[i,j] + half*dt*f_k[i,j]
         p_tmp[i,j] = p_n[i,j] + half*dt*p_k[i,j]

         f_np1[i,j] += third*dt*f_k[i,j]
         p_np1[i,j] += third*dt*p_k[i,j]
      end
   end
   ## step 3
   Radial.set_d1!(f_rd1, f_tmp, dr)
   Radial.set_d1!(p_rd1, p_tmp, dr)
   Radial.set_d2!(f_rd2, f_tmp, dr)

   Sphere.angular_matrix_mult!(f_sph_lap,f_tmp,laplM)

   set_kp(p_k,
      f_rd1, f_rd2, f_sph_lap, p_rd1, f_tmp, p_tmp,
      A_pp, A_pq,
      B_pp, B_pq, B_pf
   )
   for j=1:ny
      for i=1:nx
         f_k[i,j] = p_tmp[i,j]

         f_tmp[i,j] = f_n[i,j] + dt*f_k[i,j]
         p_tmp[i,j] = p_n[i,j] + dt*p_k[i,j]

         f_np1[i,j] += third*dt*f_k[i,j]
         p_np1[i,j] += third*dt*p_k[i,j]
      end
   end
   ## step 4
   Radial.set_d1!(f_rd1, f_tmp, dr)
   Radial.set_d1!(p_rd1, p_tmp, dr)
   Radial.set_d2!(f_rd2, f_tmp, dr)

   Sphere.angular_matrix_mult!(f_sph_lap,f_tmp,laplM)

   set_kp(p_k,
      f_rd1, f_rd2, f_sph_lap, p_rd1, f_tmp, p_tmp,
      A_pp, A_pq,
      B_pp, B_pq, B_pf
   )
   for j=1:ny
      for i=1:nx
         f_k[i,j] = p_tmp[i,j]

         f_np1[i,j] += sixth*dt*f_k[i,j]
         p_np1[i,j] += sixth*dt*p_k[i,j]
      end
   end
   Radial.filter!(f_np1,f_tmp,0.5)
   Radial.filter!(p_np1,p_tmp,0.5)
   for j=1:ny
      for i=1:nx
         f_tmp[i,j] = f_np1[i,j]
         p_tmp[i,j] = p_np1[i,j]
      end
   end
   Sphere.angular_matrix_mult!(f_np1,f_tmp,fltrM)
   Sphere.angular_matrix_mult!(p_np1,p_tmp,fltrM)

   return nothing
end

end


text:
@doc """
Return the minium projected distance between two point vectors `p1` and `p2`,
i.e. the minimum absolute value of the differences of the coordinates of `p1` and `p2`.
"""
function projected_distance( p1 :: RVec, p2 :: RVec )
    minimum( abs.( p1 .- p2) )
end

@doc "Return array of projected distances for point `p1` against every point in `P`."
function projected_distance( p1 :: RVec, P :: RVecArr )
    [ projected_distance(p1, p2) for p2 ∈ P ]
end

@doc "Return array of projected distances for point `p1` against every point in `P`, but every value below `threshold` is set to `0.0`."
function projected_distance_thresholded( p1 :: RVec, P :: RVecArr, threshold :: Real = 0.1 )
    pdist_array = projected_distance( p1, P );
    pdist_array[ pdist_array .< threshold ] .= 0
    return pdist_array
end

@doc "Return Euclidean distance from `p1` to every point in `P`."
function distance( p1 :: RVec, P :: RVecArr )
    [ norm(p1 .- p2, 2) for p2 ∈ P ]
end

"Indices of vectors in `set` that violate the box constraints `lb` or `ub`."
function bad_indices( set :: RVecArr, lb :: RVec, ub = RVec )
    [ any( s .< lb ) || any( s .> ub ) for s ∈ set]
end

"Indices of vectors in `set` that conform to the box constraints `lb` or `ub`."
function good_indices( set :: RVecArr, lb :: RVec, ub = RVec )
    .!bad_indices(set, lb, ub)
end

"Delete vectors from array `seeds` that violate the box constraints `lb` and `ub`."
function discard_bad_seeds!( seeds :: RVecArr, lb :: RVec, ub = RVec )
    bad_seeds = bad_indices( seeds, lb, ub )
    deleteat!(seeds, bad_seeds)
end

# Simple scaling from and to the square [0,1]ⁿ for finite box constraints

@memoize function _width( lb :: RVec, ub :: RVec )
    @assert length(ub) == length(lb)
    @assert all( isfinite.(lb) )
    @assert all( isfinite.(ub) )
    return ub .- lb
end

function scale_to_unit_square( p :: RVec, lb :: RVec, ub :: RVec )
    w = _width(lb, ub)
    return ( p .- lb ) ./ w
end

function scale_to_unit_square!( p :: RVec, lb :: RVec, ub :: RVec ) :: Nothing
    w = _width(lb, ub)
    p .-= lb
    p ./= w;
    nothing
end

function scale_to_unit_square( P :: RVecArr, lb :: RVec, ub :: RVec )
    [ scale_to_unit_square( p, lb, ub ) for p ∈ P ]
end

function scale_to_unit_square!( P :: RVecArr, lb :: RVec, ub :: RVec ) :: Nothing
    for p ∈ P
        scale_to_unit_square!( p, lb, ub )
    end
    nothing
end

function unscale_from_unit_square( p :: RVec, lb :: RVec, ub :: RVec )
    w = _width(lb, ub)
    return lb .+ w .* p
end

function unscale_from_unit_square!( p :: RVec, lb :: RVec, ub :: RVec ) :: Nothing
    w = _width(lb, ub)
    p .*= w
    p .+= lb
    nothing
end

function unscale_from_unit_square( P :: RVecArr, lb :: RVec, ub :: RVec )
    [ unscale_from_unit_square(p, lb, ub) for p ∈ P ]
end

function unscale_from_unit_square!( P :: RVecArr, lb :: RVec, ub :: RVec ) :: Nothing
    for p ∈ P
        unscale_from_unit_square!(p, lb, ub)
    end
    nothing
end


text:
module SymbolicMarkov

export Markov, trans!, reward!, initial!, generate, @tr, @reward, @init, @markov, @s
export CTMCModel, ctmc, exrss, exrt, cexrt
export sprob, tprob, ctprob

using SymbolicDiff: AbstractSymbolic, AbstractVectorSymbolic, AbstractMatrixSymbolic, SymbolicEnv, SymbolicCache, SymbolicExpression, @expr, @bind, @vars
import SymbolicDiff: _toexpr, _eval, seval
import NMarkov: gth, stsen, stgs, stsengs, tran, mexp, mexpc
using SparseArrays: SparseMatrixCSC, spzeros
using SparseMatrix: SparseCSC, BlockCOO
import Base
import LinearAlgebra: dot

include("_markov.jl")
include("_ctmc.jl")

include("_ctmc_sprob.jl")
include("_ctmc_tprob.jl")
include("_ctmc_exr.jl")

end


text:
module CodeToString

export code_llvm_to_bytestring, @code_llvm_to_bytestring, code_native_to_bytestring, @code_native_to_bytestring, code_warntype_to_bytestring, @code_warntype_to_bytestring

macro stdout_to_bytestring(ex)
    quote
        old_stdout = STDOUT
        (rd, wr) = redirect_stdout()
        $(esc(ex))
        close(wr)
        data = readavailable(rd)
        close(rd)
        redirect_stdout(old_stdout)
        data
    end
end

function code_llvm_to_bytestring(f, types)
    @stdout_to_bytestring code_llvm(f, types)
end

macro code_llvm_to_bytestring(ex)
    quote
        @stdout_to_bytestring @code_llvm $(esc(ex))
    end
end

function code_native_to_bytestring(f, types)
    @stdout_to_bytestring code_native(f, types)
end

macro code_native_to_bytestring(ex)
    quote
        @stdout_to_bytestring @code_native $(esc(ex))
    end
end

function code_warntype_to_bytestring(f, types)
    @stdout_to_bytestring code_warntype(f, types)
end

macro code_warntype_to_bytestring(ex)
    quote
        @stdout_to_bytestring @code_warntype $(esc(ex))
    end
end

end #module


text:
# This file was generated, do not modify it. # hide
XGBC = @load XGBoostClassifier
dtc = machine(XGBC(), Xtrain, ytrain)
fit!(dtc)
ŷ = MLJ.predict(dtc, Xtrain)
cross_entropy(ŷ, ytrain) |> mean

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using Documenter, JuliaInterpreter

makedocs(
    modules = [JuliaInterpreter],
    clean = false,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "JuliaInterpreter.jl",
    authors = "Keno Fischer, Tim Holy, and others",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "ast.md",
        "internals.md",
        "dev_reference.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaDebug/JuliaInterpreter.jl.git",
)


text:
# https://adventofcode.com/2020/day/2

# For example, suppose you have the following list:

# 1-3 a: abcde
# 1-3 b: cdefg
# 2-9 c: ccccccccc

# Each line gives the password policy and then the password. The password policy
# indicates the lowest and highest number of times a given letter must appear for
#     the password to be valid. For example, 1-3 a means that the password must
#     contain a at least 1 time and at most 3 times.

# In the above example, 2 passwords are valid. The middle password, cdefg, is not;
# it contains no instances of b, but needs at least 1. The first and third passwords
# are valid: they contain one a or nine c, both within the limits of their respective
# policies.

function getpattern(row)
    i, j, char, password = match(r"^(\d+)\-(\d+) ([^:]+): (.*)$", row).captures
    return parse(Int32, i), parse(Int32, j), char[1], password
end

function part1(inputs)
    count = 0
    for row in split(inputs, '\n', keepempty=false)
        min, max, char, password = getpattern(row)

        num_found = sum(collect(password) .== char)
        if min <= num_found <= max
            count += 1
        end
    end
    return count
end

example = "
1-3 a: abcde
1-3 b: cdefg
2-9 c: ccccccccc
"

@assert part1(example) == 2

# Each policy actually describes two positions in the password, where 1 means the first
# character, 2 means the second character, and so on. (Be careful; Toboggan Corporate
# Policies have no concept of "index zero"!) Exactly one of these positions must contain
# the given letter. Other occurrences of the letter are irrelevant for the purposes of
# policy enforcement.

# Given the same example list from above:

#     1-3 a: abcde is valid: position 1 contains a and position 3 does not.
#     1-3 b: cdefg is invalid: neither position 1 nor position 3 contains b.
#     2-9 c: ccccccccc is invalid: both position 2 and position 9 contain c.

# How many passwords are valid according to the new interpretation of the policies?

function part2(inputs)
    count = 0
    for row in split(inputs, '\n', keepempty=false)
        i, j, char, password = getpattern(row)

        if xor(password[i] == char, password[j] == char)
            count += 1
        end
    end
    return count
end

@assert part2(example) == 1

test = read("data/day-02.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 660
@assert result2 == 530

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
@testset "Nested" begin
  # nested variogram with nugget effect
  γ = NuggetEffect(0.2) + GaussianVariogram(nugget=0.1, sill=0.8, range=50.0)
  @test nugget(γ) ≈ 0.3
  @test sill(γ) ≈ 1.0
  @test range(γ) ≈ 50.0
  @test (@elapsed sill(γ)) < 1e-6
  @test (@elapsed nugget(γ)) < 1e-6
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32
  γ = 2.0*NuggetEffect(0.2)
  @test nugget(γ) ≈ 0.4
  @test sill(γ) ≈ 0.4
  @test range(γ) ≈ 0.0
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32

  # sill is defined for nested models
  γ = GaussianVariogram(sill=1.) + ExponentialVariogram(sill=2.)
  @test sill(γ) == 3.0
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32

  # nugget is defined for nested models
  γ₁ = GaussianVariogram()
  γ₂ = GaussianVariogram() + ExponentialVariogram()
  @test nugget(γ₁) == nugget(γ₂)

  # stationarity of nested models
  γ = GaussianVariogram() + ExponentialVariogram() + SphericalVariogram()
  @test isstationary(γ)
  @test sill(γ) == 3.0
  @test !isstationary(γ + PowerVariogram())
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32

  # result type is defined for nested models
  # see https://github.com/JuliaEarth/GeoStats.jl/issues/121
  γ = GaussianVariogram() + ExponentialVariogram()
  @test Variography.result_type(γ, rand(Point3), rand(Point3)) == Float64
  γ = GaussianVariogram(sill=1.0f0,range=1.0f0,nugget=0.1f0)
  @test Variography.result_type(γ, rand(Point3f), rand(Point3f)) == Float32

  # nested model with matrix coefficients
  C₁ = [1.0 0.5; 0.5 2.0]
  C₂ = [3.0 0.0; 0.0 3.0]
  γ = C₁ * GaussianVariogram(range=1.0) + C₂ * SphericalVariogram(range=2.0)
  @test range(γ) ≈ 2.0
  @test sill(γ)  ≈ C₁ .+ C₂
  @test γ(10.0) ≈ sill(γ)
  @test γ(Point(10.,0.), Point(0.,0.)) ≈ sill(γ)
  @test isstationary(γ)

  # nested model with matrix coefficients
  C = [1. 0.; 0. 1.]
  γ = C*GaussianVariogram() + C*ExponentialVariogram() + C*CubicVariogram()
  @test range(γ) ≈ 1.0
  @test sill(γ) ≈ [3. 0.; 0. 3.]
  @test γ(10.0) ≈ sill(γ)
  @test γ(Point(10.,0.), Point(0.,0.)) ≈ sill(γ)
  @test isstationary(γ)

  # test constructor explicitly
  γ = NestedVariogram((1., 2.), (ExponentialVariogram(), SphericalVariogram()))
  @test sill(γ) == 3.0
  @test range(γ) == 1.0
  @test nugget(γ) == 0.0
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32

  # test individual structures
  γ = SphericalVariogram() + 2ExponentialVariogram() + NuggetEffect(10.0)
  @test structures(γ) == (10.0, (1.0, 2.0), (SphericalVariogram(), ExponentialVariogram()))
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32
  γ = SphericalVariogram(sill=2.0) + ExponentialVariogram(nugget=0.1)
  @test structures(γ) == (0.1, (2.0, 0.9), (SphericalVariogram(), ExponentialVariogram()))
  @test structures(SphericalVariogram()) == (0.0, (1.0,), (SphericalVariogram(),))
  @test (@elapsed sill(γ)) < 1e-5
  @test (@elapsed nugget(γ)) < 1e-5
  @test (@allocated sill(γ)) < 32
  @test (@allocated nugget(γ)) < 32

  # nested model with change of support
  γ = GaussianVariogram() + SphericalVariogram()
  p = Point(1., 2., 3.)
  h = CartesianGrid(10,10,10)[1]
  @test γ(p, h) == γ(h, p)
  @test γ(h, h) < γ(p, h)
end


text:
module CGmethods
    using LinearAlgebra
    import ..Verbose_print:Verbose_level,Verbose_3,Verbose_2,Verbose_1,println_verbose3

    function add!(b,Y,a,X) #b*Y + a*X -> Y
        axpby!(a,X,b,Y) #X*a + Y*b -> Y
    end

    function add!(Y,a,X) #Y + a*X -> Y
        axpby!(a,X,1,Y) #X*a + Y -> Y
    end

    function bicg(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_2()) #Ax=b
        println_verbose3(verbose,"--------------------------------------")
        println_verbose3(verbose,"bicg method")
        res = deepcopy(b)
        temp1 = similar(x)
        p = similar(x)
        q = similar(x)
        s = similar(x)

        mul!(temp1,A,x)
        add!(res,-1,temp1)

        rnorm = real(res⋅res)

        if rnorm < eps
            return
        end
        #println(rnorm)

        mul!(p,A',res)
        c1 = p ⋅ p

        for i=1:maxsteps
            mul!(q,A,p)
            #! ...  c2 = < q | q >
            c2 = q ⋅ q

            alpha = c1 / c2
            #! ...  x   = x   + alpha * p
            add!(x,alpha,p)
            #...  res = res - alpha * q
            add!(res,-alpha,q)
            rnorm = real(res ⋅ res)
            println_verbose3(verbose,"$i-th eps: $rnorm")

            if rnorm < eps
                println_verbose3(verbose,"Converged at $i-th step. eps: $rnorm")
                println_verbose3(verbose,"--------------------------------------")
                return
            end

            mul!(s,A',res)

            c3 = s * s

            beta = c3 / c1
            c1 = c3

            add!(beta,p,1,s) #p = beta*p + s

        end

        error("""
        The BICG is not converged! with maxsteps = $(maxsteps)
        residual is $rnorm
        maxsteps should be larger.""")

    end

    function cg(x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_2()) #Ax=b
        println_verbose3(verbose,"--------------------------------------")
        println_verbose3(verbose,"cg method")
        res = deepcopy(b)
        temp1 = similar(x)

        q = similar(x)

        mul!(temp1,A,x)
        add!(res,-1,temp1)
        p = deepcopy(res)

        rnorm = real(res⋅res)
        #println(rnorm)

        if rnorm < eps
            return
        end

        c1 = p ⋅ p

        for i=1:maxsteps
            mul!(q,A,p)

            c2 = p ⋅ q

            α = c1 / c2
            #! ...  x   = x   + alpha * p
            add!(x,α,p)
            #...  res = res - alpha * q
            add!(res,-α,q)
            c3 = res ⋅ res
            rnorm = real(c3)
            println_verbose3(verbose,"$i-th eps: $rnorm")

            if rnorm < eps
                println_verbose3(verbose,"Converged at $i-th step. eps: $rnorm")
                println_verbose3(verbose,"--------------------------------------")
                return
            end

            β = c3 / c1
            c1 = c3

            add!(β,p,1,res) #p = beta*p + s

        end

        error("""
        The CG is not converged! with maxsteps = $(maxsteps)
        residual is $rnorm
        maxsteps should be larger.""")

    end

    function Base.:*(x::Array{T,1},y::Array{T,1}) where T <: Number
        return x'*y
    end

    function shiftedcg(vec_x,vec_β,x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_2()) #Ax=b
        println_verbose3(verbose,"--------------------------------------")
        println_verbose3(verbose,"shifted cg method")
        N = length(vec_β)
        temp1 = similar(b)
        r = deepcopy(b)
        p = deepcopy(b)
        q = similar(b)

        vec_r = Array{typeof(r),1}(undef,N)
        vec_p = Array{typeof(p),1}(undef,N)
        for j=1:N
            vec_r[j] = deepcopy(b)
            vec_p[j] = deepcopy(b)
        end

        αm = 1.0
        βm = 0.0

        ρm = ones(ComplexF64,N)
        ρ0 = ones(ComplexF64,N)
        ρp = ones(ComplexF64,N)
        residual = 0

        for i=1:maxsteps
            mul!(q,A,p)

            pAp = p ⋅ q
            rr = r*r
            αk = rr / pAp

            #! ...  x   = x   + alpha * p
            add!(x,αk,p)

            #...  r = r - alpha * q
            add!(r,-αk,q)

            βk = (r*r)/ rr
            add!(βk,p,1,r) #p = beta*p + r

            for j=1:N
                ρkj = ρ0[j]
                if abs(ρkj) < eps
                    continue
                end
                ρkmj =ρm[j]
                ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*vec_β[j])+αk*βm*(ρkmj-ρkj))
                αkj = (ρp[j]/ρkj)*αk

                add!(vec_x[j],αkj,vec_p[j])
                βkj = (ρp[j]/ρkj)^2*βk
                add!(βkj,vec_p[j],ρp[j],r)

            end

            ρm[:] = ρ0[:]
            ρ0[:] = ρp[:]
            αm = αk
            βm = βk

            ρMAX = maximum(abs.(ρp))^2
            residual = abs(rr*ρMAX)
            println_verbose3(verbose,"$i-th eps: $residual")

            if abs(residual) < eps
                println_verbose3(verbose,"Converged at $i-th step. eps: $residual")
                println_verbose3(verbose,"--------------------------------------")
                return
            end

        end

        error("""
        The shifted CG is not converged! with maxsteps = $(maxsteps)
        residual is $residual
        maxsteps should be larger.""")

    end

    function reducedshiftedcg(leftvec,vec_β,x,A,b;eps=1e-10,maxsteps = 1000,verbose = Verbose_2()) #Ax=b
        println_verbose3(verbose,"--------------------------------------")
        println_verbose3(verbose,"shifted cg method")
        N = length(vec_β)
        temp1 = similar(b)
        r = deepcopy(b)
        p = deepcopy(b)
        q = similar(b)

        #=
        vec_r = Array{typeof(r),1}(undef,N)
        vec_p = Array{typeof(p),1}(undef,N)
        for j=1:N
            vec_r[j] = deepcopy(b)
            vec_p[j] = deepcopy(b)
        end
        =#

        Σ = leftvec ⋅ b

        θ = zeros(ComplexF64,N)
        Π = ones(ComplexF64,N) .* Σ

        αm = 1.0
        βm = 0.0

        ρm = ones(ComplexF64,N)
        ρ0 = ones(ComplexF64,N)
        ρp = ones(ComplexF64,N)
        residual = 0

        for i=1:maxsteps
            mul!(q,A,p)

            pAp = p ⋅ q
            rr = r*r
            αk = rr / pAp

            #! ...  x   = x   + alpha * p
            add!(x,αk,p)

            #...  r = r - alpha * q
            add!(r,-αk,q)

            βk = (r*r)/ rr
            add!(βk,p,1,r) #p = beta*p + r

            Σ = leftvec ⋅ r

            for j=1:N
                ρkj = ρ0[j]
                if abs(ρkj) < eps
                    continue
                end
                ρkmj =ρm[j]
                ρp[j] = ρkj*ρkmj*αm/(ρkmj*αm*(1.0+αk*vec_β[j])+αk*βm*(ρkmj-ρkj))
                αkj = (ρp[j]/ρkj)*αk
                θ[j] += αkj*Π[j]

                #add!(vec_x[j],αkj,vec_p[j])
                βkj = (ρp[j]/ρkj)^2*βk
                Π[j] = βkj*Π[j] + ρp[j]*Σ

                #mul!(Π,ρp[j],Σ,1,βkj)
                #Π[:] =
                #add!(βkj,vec_p[j],ρp[j],r)

            end

            ρm[:] = ρ0[:]
            ρ0[:] = ρp[:]
            αm = αk
            βm = βk

            ρMAX = maximum(abs.(ρp))^2
            residual = abs(rr*ρMAX)
            println_verbose3(verbose,"$i-th eps: $residual")

            if abs(residual) < eps
                println_verbose3(verbose,"Converged at $i-th step. eps: $residual")
                println_verbose3(verbose,"--------------------------------------")
                return θ
            end

        end

        error("""
        The shifted CG is not converged! with maxsteps = $(maxsteps)
        residual is $residual
        maxsteps should be larger.""")

    end

end


text:
__precompile__()
module SparseInterp
using Base.Threads
include("utils.jl")
include("univariategrids.jl")
include("grids.jl")
include("jlfuncs.jl")
include("adapt.jl")
include("api.jl")

export Linear, Quadratic,
       getW,
       interpolate,
       level,
       grow!,
       @threads,
       NGrid
end


text:
export InputParams

"""
    InputParams

A structure for inputing parameters necessary for processing ambient noise data.

## Fields: InputParams
| **Field** | **Description** |
|:--------------------|:--------------------|
| :pair               | Network.Station.Channel for each station to be correlated. |
| :FFTOUT             | Directory filepath for FFT output. |
| :CORROUT            | Directory filepath for correlation output. |
| :cc_len             | Length of each correlation in seconds. |
| :cc_step            | Spacing between correlation windows in seconds. |
| :fs                 | Sampling frequency in Hz to downsample raw data; set to 0.0. |
| :freqmin            | Minimum frequency for whitening.  |
| :freqmax            | Maximum frequency for whitening. |
| :time_norm          | Apply one-bit whitening with "one_bit". |
| :to_whiten          | Apply whitening to FFTs. |
| :maxlag             | Maximum lag time in seconds to keep from correlations. |
| :smoothing_half_win | Half window in samples for FFT smoothing. |
| :corr_type          | Type of correlation e.g. "cross-correlate", "deconv" or "coherence" |
"""
mutable struct InputParams
           pair::Array{String,1}
           FFTOUT::String
           CORROUT::String
           cc_len::Int
           cc_step::Int
           fs::Float64
           freqmin::Float64
           freqmax::Float64
           time_norm::Union{Bool,String}
           to_whiten::Bool
           maxlag::Float64
           smoothing_half_win::Int
           corr_type::String

     function InputParams(
           pair::Array{String,1},
           FFTOUT::String,
           CORROUT::String,
           cc_len::Int,
           cc_step::Int,
           fs::Float64,
           freqmin::Float64,
           freqmax::Float64,
           time_norm::Union{Bool,String},
           to_whiten::Bool,
           maxlag::Float64,
           smoothing_half_win::Int,
           corr_type::String)

           return new(pair,FFTOUT,CORROUT,cc_len,cc_step,fs,freqmin,freqmax,
                      time_norm, to_whiten, maxlag, smoothing_half_win,
                      corr_type)
           end
end

InputParams(;
            pair::Array{String,1}         = Array{String,1}(undef, 0),
            FFTOUT::String                = "",
            CORROUT::String               = "",
            cc_len::Int                   = zero(Int),
            cc_step::Int                  = zero(Int),
            fs::Float64                   = zero(Float64),
            freqmin::Float64              = zero(Float64),
            freqmax::Float64              = zero(Float64),
            time_norm::Union{Bool,String} = false,
            to_whiten::Bool               = false,
            maxlag::Float64               = zero(Float64),
            smoothing_half_win::Int       = zero(Int),
            corr_type::String             = ""
            ) = InputParams(pair,FFTOUT,CORROUT,cc_len,cc_step,fs,freqmin,freqmax,
                       time_norm, to_whiten, maxlag, smoothing_half_win,
                       corr_type)


text:
export apply_qubit_op, get_transition_amplitude_with_obs

"""
Divide a qubit operator into the hermite and antihermite parts.
"""
divide_real_imag(op::QubitOperator) =
    (op+hermitian_conjugated(op))/2, (op-hermitian_conjugated(op))/2im

"""
Apply a qubit operator op to |state_ket> and fit the result with
circuit * |state_bra>.
The circuit object will be updated on exit.
The squared norm of op * |state_ket> will be returned.
"""
function apply_qubit_op(
    op::QubitOperator,
    state_ket::QuantumState,
    circuit::UCCQuantumCircuit, state0_bra::QuantumState,
    minimizer=mk_scipy_minimize()
    )
    n_qubit = get_n_qubit(state0_bra)
    her, antiher = divide_real_imag(op)

    function cost(thetas::Vector{Float64})
        update_circuit_param!(circuit, thetas)
        re_ = get_transition_amplitude_with_obs(circuit, state0_bra, her, state_ket)
        im_ = get_transition_amplitude_with_obs(circuit, state0_bra, antiher, state_ket)
        -abs2(re_ + im_ * im)
    end

    opt_thetas = minimizer(cost, copy(circuit.thetas))
    -cost(opt_thetas)
end

"""
Compute <state_bra| circuit^+ obs |state_ket>, where obs is a hermite observable.
"""
function get_transition_amplitude_with_obs(
    circuit::UCCQuantumCircuit, state0_bra::QuantumState,
    op::QubitOperator,
    state_ket::QuantumState)
    state_bra = copy(state0_bra)
    update_quantum_state!(circuit, state_bra)
    get_transition_amplitude(op, state_bra, state_ket)
end


text:
const StoreType = Union{<:Tuple, <:Vector}

"""
    LittleDict(keys, vals)<:AbstractDict

A ordered dictionary type for small numbers of keys.
Rather than using `hash` or some other sophisicated measure
to store the vals in a clever arrangement,
it just keeps everything in a pair of lists.

While theoretically this has expected time complexity _O(n)_,
vs the hash-based `OrderDict`/`Dict`'s expected time complexity _O(1)_,
and the search-tree-based `SortedDict`'s expected time complcity _O(log(n))_.
In practice it is really fast, because it is cache & SIMD friendly.

It is reasonable to expect it to outperform an `OrderedDict`,
with up to around 30 elements in general;
or with up to around 50 elements if using a `LittleDict` backed by `Tuples`
(see [`freeze`](@ref))
However, this depends on exactly how long `isequal` and `hash` take,
as well as on how many hash collisions occur etc.
"""
struct LittleDict{K,V,KS<:StoreType,VS<:StoreType} <: AbstractDict{K, V}
    keys::KS
    vals::VS

    function LittleDict{K,V,KS,VS}(keys,vals) where {K,V,KS,VS}
        if length(keys) != length(vals)
            throw(ArgumentError(
                "Number of keys ($(length(keys))) differs from " *
                "number of values ($(length(vals))"
            ))
        end
        K<:eltype(KS) || ArgumentError("Invalid store type $KS, for key type $K")
        V<:eltype(VS) || ArgumentError("Invalid store type $VS, for value type $K")

        return new(keys,vals)
    end
end

function LittleDict{K,V}(ks::KS, vs::VS) where {K,V, KS<:StoreType,VS<:StoreType}
    return LittleDict{K, V, KS, VS}(ks, vs)
end

function LittleDict(ks::KS, vs::VS) where {KS<:StoreType,VS<:StoreType}
    return LittleDict{eltype(KS), eltype(VS)}(ks, vs)
end

# Other iterators should be copyed to a Vector
LittleDict(ks, vs) = LittleDict(collect(ks), collect(vs))

function LittleDict{K,V}(itr) where {K,V}
    ks = K[]
    vs = V[]
    for val in itr
        if !(val isa Union{Tuple{<:Any, <:Any}, Pair})
            throw(ArgumentError(
                "LittleDict(kv): kv needs to be an iterator of tuples or pairs")
            )
        end
        k, v = val
        push!(ks, k)
        push!(vs, v)
    end
    return LittleDict(ks, vs)
end

LittleDict{K,V}(itr...) where {K,V} = LittleDict{K,V}(itr)
LittleDict(itr...) = LittleDict(itr)
LittleDict(itr::T) where T = LittleDict{kvtype(eltype(T))...}(itr)

# Avoid contention between the core constructor, and the list of elements
LittleDict(itr1::Pair, itr2::Pair) = LittleDict(first.([itr1, itr2]), last.([itr1,itr2]))
LittleDict(itr1::Pair) = LittleDict([first(itr1)], [last(itr1)])

kvtype(::Any) = (Any, Any)
kvtype(::Type{Union{}}) = (Any,Any)

kvtype(::Type{Pair{K,V}}) where {K,V} = (K,V)
kvtype(::Type{Pair{<:Any,V}}) where {V} = (Any,V)
kvtype(::Type{Pair{K,<:Any}}) where {K} = (K,Any)

kvtype(::Type{Tuple{K,V}}) where {K,V} = (K,V)
kvtype(::Type{Tuple{<:Any,V}}) where {V} = (Any,V)
kvtype(::Type{Tuple{K,<:Any}}) where {K} = (K,Any)

"""
    freeze(dd::AbstractDict)
Render an dictionary immutable by converting it to a `Tuple` backed
`LittleDict`.
The `Tuple` backed `LittleDict` is faster than the `Vector` backed `LittleDict`,
particularly when the keys are all concretely typed.
"""
function freeze(dd::AbstractDict)
    ks = Tuple(keys(dd))
    vs = Tuple(values(dd))
    return LittleDict(ks, vs)
end

isordered(::Type{<:LittleDict}) = true

# For now these are internal UnionAlls for dispatch purposes
const UnfrozenLittleDict{K,V} = LittleDict{K,V, Vector{K}, Vector{V}}
const FrozenLittleDict{K,V} = LittleDict{K,V, <:Tuple, <:Tuple}

##### Methods that all AbstractDicts should implement

Base.length(dd::LittleDict) = length(dd.keys)

function Base.getkey(dd::LittleDict, key, default)
    if key ∈ dd.keys
        return key
    else
        return default
    end
end

struct NotFoundSentinel end  # Struct to mark not not found
function Base.get(dd::LittleDict, key, default)
    @assert length(dd.keys) == length(dd.vals)
    for ii in 1:length(dd.keys)
        cand = @inbounds dd.keys[ii]
        isequal(cand, key) && return @inbounds(dd.vals[ii])
    end
    return default
end
function get(default::Base.Callable, dd::LittleDict, key)
    got = get(dd, key, NotFoundSentinel())
    if got isa NotFoundSentinel  # not found
        return default()
    else
        return got
    end
end

function Base.iterate(dd::LittleDict, ii=1)
    ii > length(dd.keys) && return nothing
    return (dd.keys[ii] => dd.vals[ii], ii+1)
end

function merge(d1::LittleDict, d2::AbstractDict)
    return merge((x,y)->y, d1, d2)
end

function merge(
    combine::Function,
    d::LittleDict,
    others::AbstractDict...
    )
    K,V = _merge_kvtypes(d, others...)
    dc = LittleDict{K,V}(d)
    for d2 in others
        for (k2,v2) in d2
            got = get(dc, k2, NotFoundSentinel())
            if got isa NotFoundSentinel
                add_new!(dc, k2, v2)
            else
                # GOLDPLATE: ideally we would avoid iterating this twice
                # once for get and once for setindex!
                dc[k2]=combine(got, v2)
            end
        end
    end
    return dc
end

Base.empty(dd::LittleDict{K,V}) where {K,V} = LittleDict{K,V}()

######## Methods that all mutable AbstractDict's should implement

function Base.sizehint!(dd::UnfrozenLittleDict, sz)
    sizehint!(dd.keys, sz)
    sizehint!(dd.vals,sz)
    return dd
end

function add_new!(dd::UnfrozenLittleDict{K, V}, key, value) where {K, V}
    kk = convert(K, key)
    vv = convert(V, value)

    # if we can convert it to the right type, and the dict is unfrozen
    # then neither push can fail, so the dict length with remain in sync
    push!(dd.keys, kk)
    push!(dd.vals, vv)

    return dd
end

function Base.setindex!(dd::LittleDict{K,V, <:Any, <:Vector}, value, key) where {K,V}
    # Note we only care if the Value store is mutable (<:Vector)
    # As we can have immutable keys, if we are setting the value of an existing key

    # Assertion below commented out as by standards of carefully optimised
    # setindex! it has huge code (26%), this does mean that if someone has messed
    # with the fields of the LittleDict directly, then the @inbounds could be invalid
    #@assert length(dd.keys) == length(dd.vals)

    kk = convert(K, key)
    vv = convert(V, value)
    for ii in 1:length(dd.keys)
        cand = @inbounds dd.keys[ii]
        if isequal(cand, kk)
            @inbounds(dd.vals[ii] = vv)
            return dd
        end
    end
    add_new!(dd, key, value)
    return dd
end

function Base.pop!(dd::UnfrozenLittleDict)
    pop!(dd.keys)
    return pop!(dd.vals)
end

function Base.pop!(dd::UnfrozenLittleDict, key)
    @assert length(dd.keys) == length(dd.vals)

    for ii in 1:length(dd.keys)
        cand = @inbounds dd.keys[ii]
        if isequal(cand, key)
            deleteat!(dd.keys, ii)
            val = @inbounds dd.vals[ii]
            deleteat!(dd.vals, ii)
            return val
        end
    end
    # Not found, throw error
    throw(KeyError(key))
end

function Base.delete!(dd::UnfrozenLittleDict, key)
    pop!(dd, key)
    return dd
end

Base.empty!(dd::UnfrozenLittleDict) = (empty!(dd.keys); empty!(dd.vals); dd)

function get!(default::Base.Callable, dd::UnfrozenLittleDict, key)
    got = get(dd, key, NotFoundSentinel())
    if got isa NotFoundSentinel  # not found
        val = default()
        add_new!(dd, key, val)
        return val
    else
        return got
    end
end
get!(dd::UnfrozenLittleDict, key, default) = get!(()->default, dd, key)


text:
qs = 1.0 .+ 2.0.^(-5:2)
times = Array{Float64}(undef,length(qs),4)
means = Array{Float64}(undef,length(qs),4)

using StochasticDiffEq, DiffEqProblemLibrary, Random,
      Plots, ParallelDataTransfer, DiffEqMonteCarlo, Distributed
Random.seed!(99)

using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
full_prob = DiffEqProblemLibrary.SDEProblemLibrary.oval2ModelExample(largeFluctuations=true,useBigs=false)
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additivesystem,
            prob_sde_additive, prob_sde_2Dlinear, prob_sde_linear, prob_sde_wave
prob = remake(full_prob,tspan=(0.0,1.0))

println("Solve once to compile.")
sol = solve(prob,EM(),dt=1/2^(18))
Int(sol.u[end][1]!=NaN)
println("Compilation complete.")
num_runs = 10000

probs = Vector{SDEProblem}(undef,3)
p1 = Vector{Any}(undef,3)
p2 = Vector{Any}(undef,3)
p3 = Vector{Any}(undef,3)
## Problem 1
probs[1] = prob_sde_linear
## Problem 2
probs[2] = prob_sde_wave
## Problem 3
probs[3] = prob_sde_additive

println("Setup Complete")

## Timing Runs

function runAdaptive(i,k)
  sol = solve(prob,SRIW1(),dt=1/2^(8),abstol=2.0^(-15),reltol=2.0^(-10),
              verbose=false,maxIters=Int(1e12),qmax=qs[k])
  Int(any(isnan,sol[end]) || sol.t[end] != 1)
end

#Compile
monte_prob = MonteCarloProblem(probs[1])
test_mc = solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,num_monte=1000,abstol=2.0^(-1),reltol=0)
DiffEqBase.calculate_monte_errors(test_mc);

for k in eachindex(qs)
  global times
  Random.seed!(99)
  adaptiveTime = @elapsed numFails = sum(map((i)->runAdaptive(i,k),1:num_runs))
  println("k was $k. The number of Adaptive Fails is $numFails. Elapsed time was $adaptiveTime")
  times[k,4] = adaptiveTime
end

for k in eachindex(probs)
  global probs, times, means, qs
  println("Problem $k")
  ## Setup
  prob = probs[k]

  for i in eachindex(qs)
    msim = solve(monte_prob,dt=1/2^(4),SRIW1(),adaptive=true,num_monte=num_runs,abstol=2.0^(-13),reltol=0,qmax=qs[i])
    test_msim = DiffEqBase.calculate_monte_errors(msim)
    times[i,k] = test_msim.elapsedTime
    means[i,k] = test_msim.error_means[:final]
    println("for k=$k and i=$i, we get that the error was $(means[i,k]) and it took $(times[i,k]) seconds")
  end
end

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
struct Point
    x::Integer
    y::Integer
    z::Integer
end

function getNeighbours(basePoint::Point, maxSize::Integer)
    neighbours = [
    Point(1,0,0)
    Point(-1,0,0),
    Point(0,1,0),
    Point(0,-1,0),
    Point(0,0,1),
    Point(0,0,-1)
    ]
    # Remove invalid points

    return neighbours
end

function buildArray(arraySize::Point)
    return BitArray(undef, arraySize.x, arraySize.y, arraySize.z)
end

max_stars_repo_name: MattPotticary/floodFill

max_issues_repo_name: MattPotticary/floodFill

max_forks_repo_name: MattPotticary/floodFill


text:
using Test
@time @testset "Vibration in fluids" begin
include("test_basics.jl")
include("test_vib.jl")
end
true

max_stars_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_issues_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids

max_forks_repo_name: PetrKryslUCSD/FinEtoolsVibInFluids


text:
module BakersMap
export bakers_map
using ..ExampleBase: LEDemo, DiscreteExample

"""
Baker's map

* <https://en.wikipedia.org/wiki/Baker%27s_map>
"""
function bakers_map(;
        u0=[0.6, 0.4],
        t_renorm=10,
        t_attr=100000,
        atol=0, rtol=1e-5,
        terminator_options = [:atol => atol, :rtol => rtol],
        kwargs...)
    @inline function phase_dynamics!(u_next, u, p, t)
        if u[1] < 0.5
            u_next[1] = 2 * u[1]
            u_next[2] = u[2] / 2
        else
            u_next[1] = 2 - 2 * u[1]
            u_next[2] = 1 - u[2] / 2
        end
    end
    @inline @views function tangent_dynamics!(u_next, u, p, t)
        phase_dynamics!(u_next[:, 1], u[:, 1], p, t)
        if u[1, 1] < 0.5
            u_next[1, 2:end] .= 2 .* u[1, 2:end]
            u_next[2, 2:end] .= u[2, 2:end] ./ 2
        else
            u_next[1, 2:end] .= - 2 .* u[1, 2:end]
            u_next[2, 2:end] .= - u[2, 2:end] ./ 2
        end
    end
    LEDemo(DiscreteExample(
        "Baker's map",
        phase_dynamics!, u0, t_renorm, nothing,
        tangent_dynamics!,
        t_attr,
        log.([2.0, 0.5]),       # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end


text:
using Documenter, Bio

custom_deps() = run(`pip install --user pygments mkdocs mkdocs-biojulia`)
@osx? makedocs(doctest = false) : makedocs()
deploydocs(
           repo = "github.com/BioJulia/Bio.jl.git",
           julia = "0.4",
           osname = "linux",
           deps = custom_deps
           )


text:
using Base
using GraffSDK

# 1. Get a Synchrony configuration
# Assume that you're running in local directory
cd(joinpath(Pkg.dir("GraffSDK"),"examples"))
config = loadGraffConfig("synchronyConfig.json")
#Update the config (updates internal config too as it's by reference)
robotId = "NewRobotSam4"
config.robotId = "NewRobotSam4"

# 2a. Robot creation and retrieval
# Creating a bot allows us to create session data for that robot.
# and associate it with data sessions
if !isRobotExisting(robotId)
    newRobot = RobotRequest(robotId, "My New Bot", "Cloudy Robot", "Active")
    retRobot = addRobot(newRobot)
end

# Now we can get it as well if we want, just asking for our default robot set in the config
getRobotDetails = getRobot()

# 2b. Get robots
# Just for example's sake - let's retrieve all robots associated with our user
robots = getRobots()
@show robots

# 3. Set some user configuration parameters for the robot
robotConfig = getRobotConfig()
testConfiguration = Dict{String, String}("myUserConfigSetting" => "Testing", "MoreSettings" => "Store your robot settings here!")
updatedConfig = updateRobotConfig(testConfiguration)
refreshedConfig = getRobotConfig()


text:
using Plots
gr(titlefontsize=12,size=(500,300),colorbar=false)

using FourierGPE

# Units: Oscillator
L = (15.,15.,15.)
N = (64,64,64)
sim = Sim(L,N)
@unpack_Sim sim

import FourierGPE.V
V(x,y,z,t) = 0.5*(x^2 + y^2 + z^2)

## set simulation parameters
γ = 0.5
tf = 2/γ
Nt = 200
t = LinRange(0.,tf,Nt)

## useful state functions
ψ0(x,y,z,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,z,0.0)/μ,0.0)+im*0.0)
healing(x,y,z,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,z,μ,g)))

## make initial state
x,y,z = X
ψi = ψ0.(x,y',reshape(z,1,1,length(z)),μ,g)
ψi = randn(N)+im*randn(N)
ϕi = kspace(ψi,sim)

@pack_Sim! sim

## Evolve in k space
sol = runsim(sim)
##

# pull out the ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
zmid = findfirst(z .≈ 0)
showpsi(x,y,ψg[:,:,zmid])

plot(x,y,abs2.(ψg[:,:,zmid]),linetype=:surface)

## Dynamics of Kohn mode
# kick it
simk = Sim(L,N)

ki = .7
ϕk = copy(ϕg)
ψk = xspace(ϕk,sim)
ψk = ψk.*exp.(im*ki*x)
ϕk = kspace(ψk,sim)
ϕi = copy(ϕk)
γ = 0.0
tf = 2*pi
t = LinRange(0.,tf,Nt)

@pack_Sim! simk

@time solk = runsim(simk)

## animate a slice
anim = @animate for i=1:Nt
    ϕ = solk[i]
    zmid = findfirst(z .≈ 0)
    ψ = xspace(ϕ,simk)[:,:,zmid]
    showpsi(x,y,ψ)
end

gif(anim,"./examples/3dKohn.gif",fps=30)
# gif(anim,"./examples/3dKohn.mov",fps=30)

# adimate slice as a surface
anim = @animate for i=1:Nt
    ϕ = solk[i]
    zmid = findfirst(z .≈ 0)
    ψ = xspace(ϕ,simk)[:,:,zmid]
    plot(x,y,abs2.(ψ),linetype=:surface)
end

gif(anim,"./examples/3dKohnSurface.mov",fps=30)

simk

## animate isosurface in Makie
using Makie, AbstractPlotting

function dense(phi)
    ψm = xspace(phi,sim)
    density = abs2.(ψm)
    pmax = maximum(density)
    return density/pmax
end

function densityfilm(sol,Nt)
    saveto="examples/3dKohn.gif"
    scene = Scene()
    tindex = Node(1)
    scene = volume(lift(i -> dense(sol[i]), tindex),
    algorithm = :iso,
    show_axis=false,
    isovalue=3f0(.15))

    R = 70
    eyepos = Vec3f0(R,R,R)
    lookat = Vec3f0(18,18,0)

    record(scene, saveto, 1:Nt) do i
        update_cam!(scene, eyepos, lookat)
        rotate_cam!(scene, 0., -0.25, 0.)
        tindex[] = i
    end
    p = scene[end];
    return
end

p = densityfilm(sol,Nt)

# p.attributes.attribute


text:
function string2cmd(str::String)

ind = findall(" ", str);
n = length(ind);
cmdstr = Array{String,1}(undef,n+1);
cmdstr[1] = str[1:(ind[1][1]-1)];
for i = 2:n
    cmdstr[i] = str[(ind[i-1][1]+1):(ind[i][1]-1)];
end
i = n+1;
cmdstr[i] = str[(ind[i-1][1]+1):end];

return Cmd(cmdstr)

end

function findinstallexec(filename, appname, brew, searchopt)

print("Finding " * appname * "...\n");

dirname = "";
appinstall = 0;

status = Sys.which(filename);
if status != nothing
    dirname = filename;
    print("Exasim found " * dirname * "\n");
    return dirname, appinstall;
end
status = Sys.which("/usr/bin/" * filename);
if status != nothing
    dirname = "/usr/bin/" * filename;
    print("Exasim found " * dirname * "\n");
    return dirname, appinstall;
end
status = Sys.which("/usr/local/bin/" * filename);
if status != nothing
    dirname = "/usr/local/bin/" * filename;
    print("Exasim found " * dirname * "\n");
    return dirname, appinstall;
end
status = Sys.which("/opt/local/bin/" * filename);
if status != nothing
    dirname = "/opt/local/bin/" * filename;
    print("Exasim found " * dirname * "\n");
    return dirname, appinstall;
end

if searchopt==1
    if Sys.isapple()
        a = read(string2cmd("find /Applications -name " * filename * " -type f"), String);
        if length(a)>0
            ii = findlast("/MacOS/" * filename, a);
            dirname = a[1:ii[end]];
            print("Exasim found " * dirname * "\n");
            return dirname, appinstall;
        end
    end
end

print("Exasim could not find " * appname * " on your computer.\n");
if searchopt==10
    dirname = filename;
    print("CUDA Toolkit is not found on your system.\n");
    print("If you have Nividia GPUs on your system, please visit https://docs.nvidia.com/cuda/ to install it.\n");
    return dirname, appinstall;
else
    appinstall = 1;
    if Sys.isapple()
        print("Installing " * appname * " via brew.\n");
        if searchopt==1
            run(string2cmd(brew * " cask install " * appname));
        else
            run(string2cmd(brew * " install " * appname));
        end
    elseif Sys.isunix()
        print("Installing " * appname * " via apt.\n");
        run(string2cmd("sudo apt install " * appname));
    elseif Sys.windows()
        print("Installing " * appname * " via apt.\n");
        run(string2cmd("sudo apt install " * appname));
    end

    return dirname, appinstall;
end

return dirname, appinstall;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using AmplNLWriter, Ipopt

import MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test

const OPTIMIZER = AmplNLWriter.Optimizer(Ipopt.amplexe, ["print_level = 0"])

@test sprint(show, OPTIMIZER) == "An AmplNLWriter model"

const CONFIG = MOIT.TestConfig(
    atol = 1e-4, rtol = 1e-4, optimal_status = MOI.LOCALLY_SOLVED,
    infeas_certificates = false, duals = false
)

@testset "Unit Tests" begin
    MOIT.unittest(OPTIMIZER, CONFIG, [
        "solve_objbound_edge_cases",  # ObjectiveBound not implemented
        "solve_integer_edge_cases",  # Ipopt doesn't handle integer
        "solve_affine_deletion_edge_cases",  # VectorAffineFunction
        "solve_duplicate_terms_vector_affine",  # VectorAffineFunction
        # It seems that the AMPL NL reader declares NL files with no objective
        # and no constraints as corrupt, even if they have variable bounds. Yuk.
        "solve_blank_obj"
    ])
end

@testset "Linear tests" begin
    MOIT.contlineartest(OPTIMIZER, CONFIG, [
        "linear7", "linear15"  # VectorAffineFunction
    ])
end

@testset "Quadratic tests" begin
    MOIT.contquadratictest(OPTIMIZER, CONFIG, [
        "qcp1"  # VectorAffineFunction
    ])
end

@testset "ModelLike tests" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "AmplNLWriter"
    @testset "default_objective_test" begin
         MOIT.default_objective_test(OPTIMIZER)
     end
     @testset "default_status_test" begin
         MOIT.default_status_test(OPTIMIZER)
     end
    @testset "nametest" begin
        MOIT.nametest(OPTIMIZER)
    end
    @testset "validtest" begin
        MOIT.validtest(OPTIMIZER)
    end
    @testset "emptytest" begin
        # Requires VectorOfVariables
        # MOIT.emptytest(OPTIMIZER)
    end
    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(OPTIMIZER)
    end
    @testset "copytest" begin
        # Requires VectorOfVariables
        # MOIT.copytest(OPTIMIZER, AmplNLWriter.Optimizer(Ipopt.amplexe))
    end
end

@testset "MOI NLP tests" begin
    # Requires ExprGraph in MOI tests
    # MOIT.nlptest(OPTIMIZER, CONFIG)
end


text:
#= utils.jl
Philip Barrett, pobarrett@gmail.com
27may2016, Chicago

Defines utilities for Polygon objects
=#

"""
    union( poly1::Polygon, poly2::Polygon )
Returns the convex union of two polygons
"""
function union( poly1::Polygon, poly2::Polygon )
    pts = [ poly1.pts ; poly2.pts ]
        # The matrix of points to work with
    # return grahamScan( pts )
    return chull( Polygon( pts = pts ) )
end

"""
    union( polys::Array{ Polygon, 1 } )
Returns the convex union of an array of polygons
"""
function union( polys::Array{ Polygon, 1 } )
    pts = vcat( [ polys[i].pts for i in 1:length(polys) ] ... )
        # The matrix of all points
    # return grahamScan( pts )
    return chull( Polygon( pts = pts ) )
end

"""
    crop( poly::Polygon, dim::Int, dist, upper=true )
Returns a Polygon cropped in either the x or y dimension, with upper = true
retaining the part of the Polygon above the chop, and if false the part below.
The Polygon should already be oreinted counterclockwise
"""
function crop( poly::Polygon, dim::Int, dist::Float64, upper::Bool=true )
  N = size(poly.pts)[1]

  # The crop doesn't bind
  l = minimum( poly.pts[:,dim] )
  u = maximum( poly.pts[:,dim] )
  if( upper )
      if( dist < l )
          return poly
      end
      if( dist > u )
          return Polygon()
      end
  else
      if( dist < l )
          return Polygon()
      end
      if( dist > u )
          return poly
      end
  end

  # Compute the direction of the crop
  if( dim == 1 && upper )
    dir = [ -1.0, 0.0 ]
    dist = - dist
  elseif( dim == 1 && !upper )
    dir = [ 1.0, 0.0 ]
  elseif( dim==2 && upper )
    dir = [ 0.0, -1.0 ]
    dist = - dist
  else
    dir = [ 0.0, 1.0 ]
  end
  c = dir[1]
  s = dir[2]
      # Sign and cosine of the dir vector

  # Find the points which are to be cropped
  init = 0      # Last point to be retained
  term = 0      # Last point to be cropped
      # Initial and terminal cropping points
  for( i in 1:N )
      if( ( poly.pts[i,:] * dir )[1] < dist )
          if( 0 < init < i-1 && term > 0 )
              break
          end
          init = ( init == 0 ) ? term + 1 : init + 1
      else
          if( 0 < term < i-1 && init > 0 )
              break
          end
          term = ( term == 0 ) ? init + 1 : term + 1
            # Increment the terminal counter if non-zero

      end
  end

  if( init < term )
    newdirs = [ poly.dirs[ 1:init, : ] ;
                dir' ; poly.dirs[ term:end, : ] ]
    newdists = [ poly.dists[ 1:init ]; dist;
                poly.dists[ term:end ] ]
  else
    newdirs = [ dir' ; poly.dirs[ term:init, : ] ]
    newdists = [ dist; poly.dists[ term:init ] ]
  end

  return Polygon( dirs = newdirs, dists = newdists )

end


text:
function update!(app::ApplicationState)
    app.noise = @timeit to "Generate noise" perlin(app.resolution, app.scale)
    app.haschanged = true
end

function texture_data(app::ApplicationState)
    noise_data = remap(app.noise, (0.0, 1.0))
    RGBA{Float16}.(noise_data, noise_data, noise_data, 1.0)
end

function create_staging_buffer!(app::ApplicationState, device::Device)
    # initialize local buffer
    info = BufferCreateInfo(buffer_size(app.noise), BUFFER_USAGE_TRANSFER_DST_BIT | BUFFER_USAGE_TRANSFER_SRC_BIT, SHARING_MODE_EXCLUSIVE, [0])
    local_buffer = unwrap(create_buffer(device, info))
    local_memory = DeviceMemory(local_buffer, texture_data(app))
    local_resource = GPUResource(local_buffer, local_memory, info)
    app.gpu.buffers[:staging] = local_resource
end

function transfer_texture!(gr::GUIRenderer, app::ApplicationState)
    rdr = gr.rdr
    # (re)create staging buffer if necessary
    if !haskey(app.gpu.buffers, :staging) || buffer_size(app.noise) ≠ app.gpu.buffers[:staging].info.size
        create_staging_buffer!(app, device(rdr))
    else
        upload_data(memory(app.gpu.buffers[:staging]), texture_data(app))
    end

    # transfer host-visible memory to device-local image
    image = gr.resources[:perlin][1].image.resource
    cbuffer, _... = unwrap(allocate_command_buffers(device(rdr), CommandBufferAllocateInfo(rdr.gpu.command_pools[:primary], COMMAND_BUFFER_LEVEL_PRIMARY, 1)))
    @record cbuffer begin
        # transition layout to transfer destination
        cmd_pipeline_barrier(
            PIPELINE_STAGE_TOP_OF_PIPE_BIT,
            PIPELINE_STAGE_TRANSFER_BIT,
            [],
            [],
            [
                ImageMemoryBarrier(
                    AccessFlag(0),
                    ACCESS_TRANSFER_WRITE_BIT,
                    IMAGE_LAYOUT_UNDEFINED,
                    IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                    QUEUE_FAMILY_IGNORED,
                    QUEUE_FAMILY_IGNORED,
                    image,
                    ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
                ),
            ],
        )
        cmd_copy_buffer_to_image(
            app.gpu.buffers[:staging].resource,
            image,
            IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
            [BufferImageCopy(0, app.resolution..., ImageSubresourceLayers(IMAGE_ASPECT_COLOR_BIT, 0, 0, 1), Offset3D(0, 0, 0), Extent3D(app.resolution..., 1))],
        )
        # transition to final layout
        cmd_pipeline_barrier(
            PIPELINE_STAGE_TRANSFER_BIT,
            PIPELINE_STAGE_FRAGMENT_SHADER_BIT,
            [],
            [],
            [
                ImageMemoryBarrier(
                    ACCESS_TRANSFER_WRITE_BIT,
                    AccessFlag(0),
                    IMAGE_LAYOUT_TRANSFER_DST_OPTIMAL,
                    IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL,
                    QUEUE_FAMILY_IGNORED,
                    QUEUE_FAMILY_IGNORED,
                    image,
                    ImageSubresourceRange(IMAGE_ASPECT_COLOR_BIT, 0, 1, 0, 1),
                ),
            ],
        )
    end

    transfer = CommandBufferSubmitInfoKHR(cbuffer, 0)
    submit(rdr, [SubmitInfo2KHR([], [transfer], [])])
end


text:
using Interpolations
using QuadGK

function calculate_Vdl_interp(dos_f, Vq_min, Vq_max, C_dl)
    Vq_range = range(Vq_min, Vq_max, step=0.001)
    E_min = dos_f.ranges[1][1]
    E_max = dos_f.ranges[1][end]
    CQ = [compute_cq(E_min, E_max, V, dos_f)*1.602*10 for V in Vq_range]
    Vdl_data = []
    Vappl_data = []
    for i = 1:length(Vq_range)
        V_app_i =  Vq_range[i]*((CQ[i]/C_dl) + 1)
        V_dl_i = V_app_i - Vq_range[i]
        push!(Vdl_data, V_dl_i)
        push!(Vappl_data, V_app_i)
    end

    ## sort data and do V_dl interpolation
    inds = sortperm(Vappl_data)
    v_interp = LinearInterpolation(Vappl_data[inds], Vdl_data[inds])
    return v_interp
end

function QC_integrand(E, eVq, dos_func; kT=.026)
    dos_func(E) * sech((E-eVq)/(2*kT))^2/(4*kT)
end

function compute_cq(E_min, E_max, eVq, dos_func; kT=.026)
    fcn = E->QC_integrand(E, eVq, dos_func; kT=kT)
    quadgk(fcn, E_min, E_max)[1]
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module LinAlg

import Base: \, /, *, ^, +, -, ==
import Base: A_mul_Bt, At_ldiv_Bt, A_rdiv_Bc, At_ldiv_B, Ac_mul_Bc, A_mul_Bc, Ac_mul_B,
    Ac_ldiv_B, Ac_ldiv_Bc, At_mul_Bt, A_rdiv_Bt, At_mul_B
import Base: USE_BLAS64, abs, big, broadcast, ceil, conj, convert, copy, copy!,
    adjoint, eltype, exp, eye, findmax, findmin, fill!, floor, full, getindex,
    hcat, imag, indices, inv, isapprox, isone, IndexStyle, kron, length, log, map,
    ndims, oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round,
    setindex!, show, similar, size, sqrt, transpose, trunc, typed_hcat, vec
using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
    @propagate_inbounds, @pure, reduce, typed_vcat
# We use `_length` because of non-1 indices; releases after julia 0.5
# can go back to `length`. `_length(A)` is equivalent to `length(linearindices(A))`.

export
# Modules
    LAPACK,
    BLAS,

# Types
    RowVector,
    ConjArray,
    ConjVector,
    ConjMatrix,
    SymTridiagonal,
    Tridiagonal,
    Bidiagonal,
    Factorization,
    BunchKaufman,
    Cholesky,
    CholeskyPivoted,
    Eigen,
    GeneralizedEigen,
    GeneralizedSVD,
    GeneralizedSchur,
    Hessenberg,
    LU,
    LDLt,
    QR,
    QRPivoted,
    LQ,
    Schur,
    SVD,
    Hermitian,
    Symmetric,
    LowerTriangular,
    UpperTriangular,
    Diagonal,
    UniformScaling,

# Functions
    axpy!,
    axpby!,
    bkfact,
    bkfact!,
    chol,
    cholfact,
    cholfact!,
    cond,
    condskeel,
    copy!,
    copy_transpose!,
    cross,
    adjoint,
    adjoint!,
    det,
    diag,
    diagind,
    diagm,
    diff,
    dot,
    eig,
    eigfact,
    eigfact!,
    eigmax,
    eigmin,
    eigs,
    eigvals,
    eigvals!,
    eigvecs,
    eye,
    factorize,
    givens,
    gradient,
    hessfact,
    hessfact!,
    isdiag,
    ishermitian,
    isposdef,
    isposdef!,
    issuccess,
    issymmetric,
    istril,
    istriu,
    kron,
    ldltfact!,
    ldltfact,
    linreg,
    logabsdet,
    logdet,
    lu,
    lufact,
    lufact!,
    lyap,
    norm,
    normalize,
    normalize!,
    nullspace,
    ordschur!,
    ordschur,
    peakflops,
    pinv,
    qr,
    qrfact!,
    qrfact,
    lq,
    lqfact!,
    lqfact,
    rank,
    scale!,
    schur,
    schurfact!,
    schurfact,
    svd,
    svdfact!,
    svdfact,
    svds,
    svdvals!,
    svdvals,
    sylvester,
    trace,
    transpose,
    transpose!,
    transpose_type,
    tril,
    triu,
    tril!,
    triu!,
    vecdot,
    vecnorm,

# Operators
    \,
    /,
    A_ldiv_B!,
    A_ldiv_Bc,
    A_ldiv_Bt,
    A_mul_B!,
    A_mul_Bc,
    A_mul_Bc!,
    A_mul_Bt,
    A_mul_Bt!,
    A_rdiv_Bc,
    A_rdiv_Bt,
    Ac_ldiv_B,
    Ac_ldiv_Bc,
    Ac_ldiv_B!,
    Ac_mul_B,
    Ac_mul_B!,
    Ac_mul_Bc,
    Ac_mul_Bc!,
    Ac_rdiv_B,
    Ac_rdiv_Bc,
    At_ldiv_B,
    At_ldiv_Bt,
    At_ldiv_B!,
    At_mul_B,
    At_mul_B!,
    At_mul_Bt,
    At_mul_Bt!,
    At_rdiv_B,
    At_rdiv_Bt,

# Constants
    I

const BlasFloat = Union{Float64,Float32,Complex128,Complex64}
const BlasReal = Union{Float64,Float32}
const BlasComplex = Union{Complex128,Complex64}

if USE_BLAS64
    const BlasInt = Int64
else
    const BlasInt = Int32
end

# Check that stride of matrix/vector is 1
# Writing like this to avoid splatting penalty when called with multiple arguments,
# see PR 16416
@inline chkstride1(A...) = _chkstride1(true, A...)
@noinline _chkstride1(ok::Bool) = ok || error("matrix does not have contiguous columns")
@inline _chkstride1(ok::Bool, A, B...) = _chkstride1(ok & (stride(A, 1) == 1), B...)

"""
    LinAlg.checksquare(A)

Check that a matrix is square, then return its common dimension.
For multiple arguments, return a vector.

# Examples
```jldoctest
julia> A = ones(4,4); B = zeros(5,5);

julia> LinAlg.checksquare(A, B)
2-element Array{Int64,1}:
 4
 5
```
"""
function checksquare(A)
    m,n = size(A)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
    m
end

function checksquare(A...)
    sizes = Int[]
    for a in A
        size(a,1)==size(a,2) || throw(DimensionMismatch("matrix is not square: dimensions are $(size(a))"))
        push!(sizes, size(a,1))
    end
    return sizes
end

function char_uplo(uplo::Symbol)
    if uplo == :U
        'U'
    elseif uplo == :L
        'L'
    else
        throw(ArgumentError("uplo argument must be either :U (upper) or :L (lower)"))
    end
end

copy_oftype(A::AbstractArray{T}, ::Type{T}) where {T} = copy(A)
copy_oftype(A::AbstractArray{T,N}, ::Type{S}) where {T,N,S} = convert(AbstractArray{S,N}, A)

include("conjarray.jl")
include("transpose.jl")
include("rowvector.jl")

include("exceptions.jl")
include("generic.jl")

include("blas.jl")
include("matmul.jl")
include("lapack.jl")

include("dense.jl")
include("tridiag.jl")
include("triangular.jl")

include("factorization.jl")
include("qr.jl")
include("hessenberg.jl")
include("lq.jl")
include("eigen.jl")
include("svd.jl")
include("symmetric.jl")
include("cholesky.jl")
include("lu.jl")
include("bunchkaufman.jl")
include("diagonal.jl")
include("bidiag.jl")
include("uniformscaling.jl")
include("givens.jl")
include("special.jl")
include("bitarray.jl")
include("ldlt.jl")
include("schur.jl")

include("arpack.jl")
include("arnoldi.jl")

function __init__()
    try
        BLAS.check()
        if BLAS.vendor() == :mkl
            ccall((:MKL_Set_Interface_Layer, Base.libblas_name), Void, (Cint,), USE_BLAS64 ? 1 : 0)
        end
    catch ex
        Base.showerror_nostdio(ex,
            "WARNING: Error during initialization of module LinAlg")
    end
end

end # module LinAlg

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
"""
    Pipeline([client::Client=get_global_client(); filter_multi_exec::Bool=false]) -> Pipeline

Creates a Pipeline client instance for executing commands in batch.

# Attributes
- `client::Client`: Reference to the underlying Client connection.
- `resp::Vector{String}`: Batched commands converted to RESP compliant string.
- `filter_multi_exec::Bool`: Set `true` to filter out QUEUED responses in a MULTI/EXEC transaction.
- `multi_exec::Bool`: Used to track and filter MULTI/EXEC transactions.
- `multi_exec_bitmask::Vector{Bool}`: Used to track and filter MULTI/EXEC transactions.

# Examples
```julia-repl
julia> pipe = Pipeline();

julia> set("key", "value"; client=pipe);

julia> get("key"; client=pipe);

julia> execute(pipe)
2-element Array{String,1}:
 "OK"
 "value"
```
"""
mutable struct Pipeline
    client::Client
    resp::Vector{String}
    filter_multi_exec::Bool
    multi_exec::Bool
    multi_exec_bitmask::Vector{Bool}
end
Pipeline(client::Client=get_global_client(); filter_multi_exec::Bool=false) = Pipeline(client, [], filter_multi_exec, true, [])

"""
    add!(pipe::Pipeline, command)

Add a RESP compliant command to a pipeline client.
"""
function add!(pipe::Pipeline, command::AbstractArray)
    push!(pipe.resp, resp(command))

    if pipe.filter_multi_exec
        first = uppercase(command[1])

        if first == "MULTI"
            pipe.multi_exec = false
        elseif first == "EXEC"
            pipe.multi_exec = true
        end

        push!(pipe.multi_exec_bitmask, pipe.multi_exec)
    end
end
function add!(pipe::Pipeline, command::AbstractString)
    add!(pipe, split_on_whitespace(command))
end

"""
    add!(pipe::Pipeline, command)

Flushes the underlying client socket and resets the pipeline in to a clean slate.
"""
function flush!(pipe::Pipeline)
    flush!(pipe.client)
    pipe.resp = []
    pipe.multi_exec = false
    pipe.multi_exec_bitmask = []
end

"""
    pipeline(fn::Function[; clientt=get_global_client(), filter_multi_exec=false])

Execute commands batched in a pipeline client in a do block, optionally filter out MULTI transaction
responses before the EXEC call, e.g. "QUEUED".

# Examples
```julia-repl
julia> pipeline() do pipe
           lpush("example", 1, 2, 3, 4; client=pipe)
           lpop("example"; client=pipe)
           rpop("example"; client=pipe)
           multi_exec(; client=pipe) do
               lpop("example"; client=pipe)
               rpop("example"; client=pipe)
           end
           lpop("example"; client=pipe)
       end
5-element Array{Any,1}:
 4  # Integer response from lpush
 "4"  # String response from lpop
 "1"  # String response from rpop
 ["3", "2"]  # Array of String response from multi_exec do block, with responeses before the exec call filtered out
 nothing  # Nil response from final lpop
```
"""
function Base.pipeline(fn::Function; client::Client=get_global_client(), filter_multi_exec=false)
    pipe = Pipeline(client; filter_multi_exec=filter_multi_exec)
    fn(pipe)
    return execute(pipe)
end


text:
struct GraphBased end
struct HyperNetX end

let div_id = 0
    global get_next_div_id
    function get_next_div_id()
        div_id += 1
    end
end

"""
    function draw(
            h::Hypergraph,
            type::Type{GraphBased};
            element::Union{String, Int}=get_next_div_id(),
            width::Int=500,
            height::Int=500,
            radius::Real=10,
            node_radii::Union{AbstractVector{<:Real}, Nothing}=nothing,
            node_color::String="#999",
            node_colors::Union{AbstractVector{String}, Nothing}=nothing,
            node_stroke::Union{String, Nothing} = nothing,
            node_strokes::Union{AbstractVector{String}, Nothing}=nothing,
            stroke_width::Real=0,
            stroke_widths::Union{AbstractVector{<:Real}, Nothing}=nothing,
            node_opacity::Real=1,
            node_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
            stroke_opacity::Real=1,
            stroke_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
            with_node_labels::Bool=false,
            node_labels::Union{AbstractVector{String}, Nothing}=nothing,
            with_node_metadata_hover::Bool=false,
            with_node_weight::Bool=false,
            he_colors::Union{AbstractVector{String}, Nothing}=nothing,
            with_he_labels::Bool=false,
            he_labels::Union{AbstractVector{String}, Nothing}=nothing,
            with_he_metadata_hover::Bool=false
        )

Draw a hypergraph `h` in a web-based environment (e.g. Jupyter Notebook),
using a js script based on the library (D3)[https://d3js.org/].
Each hyperedge `he` is represented by a *fake* vertex `fv` to which
each vertex `v ∈ he` is connected.

**Arguments**

* `h` : the hypergraph to draw
* `type` : how the hypergraph will be drawn. If `type=GraphBased`, each hyperedge
will be represented as a vertex (see above)
* `width` : width of the figure
* `height` : height of the figure
* `radius` : same default radius for each vertex (represented as a circle)
* `node_radii` : distinct radius values for each vertex
* `node_color` : same default color for each vertex
* `node_colors` : distinct node colors for each vertex
* `node_stroke` : same default stroke for each vertex
* `node_strokes` : distinct node strokes for each vertex
* `stroke_width` : same default stroke-width for each vertex
* `stroke_widths` : distinct stroke-width values for each vertex
* `node_opacity` : same default opacity for each vertex
* `node_opacities` : distinct node-opacity values for each vertex
* `stroke_opacity` :  same default stroke-opacity for each vertex
* `stroke_opacities` : distinct stroke-opacity values for each vertex
* `with_node_labels` : whether displaying node labels
* `node_labels` : node labels to be shown
* `with_node_metadata_hover` : whether displaying node metadata when hovering each vertex
* `with_node_weight` : whether displaying node weights within each hyperedge
* `he_colors` : distinct hyperedge colors for each hyperedge
* `with_he_labels` : whether displaying hyoeredges labels
* `with_he_metadata_hover` : whether displaying hyperedge metadata when hovering each hyperedge

"""
function draw(
        h::Hypergraph,
        type::Type{GraphBased};
        element::Union{String, Int}=get_next_div_id(),
        width::Int=500,
        height::Int=500,
        radius::Real=10,
        node_radii::Union{AbstractVector{<:Real}, Nothing}=nothing,
        node_color::String="#999",
        node_colors::Union{AbstractVector{String}, Nothing}=nothing,
        node_stroke::Union{String, Nothing} = nothing,
        node_strokes::Union{AbstractVector{String}, Nothing}=nothing,
        stroke_width::Real=0,
        stroke_widths::Union{AbstractVector{<:Real}, Nothing}=nothing,
        node_opacity::Real=1,
        node_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
        stroke_opacity::Real=1,
        stroke_opacities::Union{AbstractVector{<:Real}, Nothing}=nothing,
        with_node_labels::Bool=false,
        node_labels::Union{AbstractVector{String}, Nothing}=nothing,
        with_node_metadata_hover::Bool=false,
        with_node_weight::Bool=false,
        he_colors::Union{AbstractVector{String}, Nothing}=nothing,
        with_he_labels::Bool=false,
        he_labels::Union{AbstractVector{String}, Nothing}=nothing,
        with_he_metadata_hover::Bool=false
    )

    w = widget_graph(
        JSON3.write(h.v2he),
        JSON3.write(h.he2v),
        element;
        v_meta=h.v_meta,
        he_meta=h.he_meta,
        width=width,
        height=height,
        radius=radius,
        node_radii=node_radii,
        node_color=node_color,
        node_colors=node_colors,
        node_stroke=node_stroke,
        node_strokes=node_strokes,
        stroke_width=stroke_width,
        stroke_widths=stroke_widths,
        node_opacity=node_opacity,
        node_opacities=node_opacities,
        stroke_opacity=stroke_opacity,
        stroke_opacities=stroke_opacities,
        with_node_labels=with_node_labels,
        node_labels=node_labels,
        with_node_metadata_hover=with_node_metadata_hover,
        with_node_weight=with_node_weight,
        he_colors=he_colors,
        with_he_labels=with_he_labels,
        he_labels=he_labels,
        with_he_metadata_hover=with_he_metadata_hover
        )

    display(w)
end

"""
    draw(
        h::Hypergraph,
        type::Type{HyperNetX};
        width::Int=10,
        height::Int=10,
        node_labels::Union{Dict{Int, String}, Nothing}=nothing,
        edge_labels::Union{Dict{Int, String}, Nothing}=nothing,
        collapse_nodes::Bool=false,
        collapse_edges::Bool=false,
        pos::Union{Dict{Int,Pair{Int,Int}}, Nothing}=nothing,
        with_color::Bool=true,
        with_node_counts::Bool=false,
        with_edge_counts::Bool=false,
        layout::PyObject=nx.spring_layout,
        layout_kwargs::Dict=Dict{String, Any}(),
        ax::Union{PyObject, Nothing}=nothing,
        no_border::Bool=false,
        edges_kwargs::Dict=Dict{String, Any}(),
        nodes_kwargs::Dict=Dict{String, Any}(),
        edge_labels_kwargs::Dict=Dict{String, Any}(),
        node_labels_kwargs::Dict=Dict{String, Any}(),
        with_edge_labels::Bool=true,
        with_node_labels::Bool=true,
        label_alpha::Float64=.35
        )

 Draw a hypergraph `h` as an Euler diagram, using the library [HyperNetX](https://github.com/pnnl/HyperNetX).

 **Arguments**

 * `h` : the hypergraph to draw
 * `type` : how the hypergraph will be drawn. If `type=HyperNetX`, the hypergraph will be represented as a Euler Diagram
 * `width` : width of the figure
 * `height` : height of the figure
 * `node_labels` : node labels to be shown
 * `edge_labels` : edge labels to be shown
 * `collapse_nodes` : draws the hypergraph gotten by identifying nodes contained by the same edges (from HyperNetX)
 * `collapse_edges` : draws the hypergraph gotten by identifying edges containing the same nodes (from HyperNetX)
 * `no_border` : indicates wheter the figure should have a border

For more details about the other parameters, please refer to the library [HyperNetX](https://github.com/pnnl/HyperNetX).
"""
function draw(
        h::Hypergraph,
        type::Type{HyperNetX};
        width::Int=10,
        height::Int=10,
        node_labels::Union{Dict{Int, String}, Nothing}=nothing,
        edge_labels::Union{Dict{Int, String}, Nothing}=nothing,
        collapse_nodes::Bool=false,
        collapse_edges::Bool=false,
        pos::Union{Dict{Int,Pair{Int,Int}}, Nothing}=nothing,
        with_color::Bool=true,
        with_node_counts::Bool=false,
        with_edge_counts::Bool=false,
        layout::PyObject=nx.spring_layout,
        layout_kwargs::Dict=Dict{String, Any}(),
        ax::Union{PyObject, Nothing}=nothing,
        no_border::Bool=false,
        edges_kwargs::Dict=Dict{String, Any}(),
        nodes_kwargs::Dict=Dict{String, Any}(),
        edge_labels_kwargs::Dict=Dict{String, Any}(),
        node_labels_kwargs::Dict=Dict{String, Any}(),
        with_edge_labels::Bool=true,
        with_node_labels::Bool=true,
        label_alpha::Float64=.35
        )
    if (!SimpleHypergraphs.support_hypernetx())
        throw("HyperNetX is not installed in Python used by this Julia. Install HyperNetX and reload SimpleHypergraphs.jl")
    end

    h_hnx = _convert_to_hnx(h;
        node_labels=node_labels,
        edge_labels=edge_labels
        )

    if collapse_nodes
        h_hnx = h_hnx.collapse_nodes()
    end
    if collapse_edges
        h_hnx = h_hnx.collapse_edges()
    end

    if ax === nothing #isnothing(ax)
        fig = plt.figure(figsize=[width,height])
        ax = plt.gca()

        if no_border
            ax.axis("off")
        end
    end

    hnx.draw(h_hnx,
        pos=pos,
        with_color=with_color,
        with_node_counts=with_node_counts,
        with_edge_counts=with_edge_counts,
        layout=layout,
        layout_kwargs=layout_kwargs,
        ax=ax,
        edges_kwargs=edges_kwargs,
        nodes_kwargs=nodes_kwargs,
        edge_labels_kwargs=edge_labels_kwargs,
        node_labels_kwargs=node_labels_kwargs,
        with_edge_labels=with_edge_labels,
        with_node_labels=with_node_labels,
        label_alpha=label_alpha
    )
end

"""
    _convert_to_hnx(
        h::Hypergraph;
        node_labels::Union{Dict{Int, String}, Nothing}=nothing,
        edge_labels::Union{Dict{Int, String}, Nothing}=nothing,
    )

Convert a SimpleHypergraphs hypergraph `h` into a HyperNetX hypergraph `h_hnx`.
`node_labels` represents vertex ids in `h_hnx`. If `node_labels` is none,
h's vertex ids will be used.
`edge_labels` represent edge ids in `h_hnx`. If `edge_labels` is none,
h's edge ids will be used.
"""
function _convert_to_hnx(h::Hypergraph;
        node_labels::Union{Dict{Int, String}, Nothing}=nothing,
        edge_labels::Union{Dict{Int, String}, Nothing}=nothing,
        )
    h_hnx = hnx.Hypergraph()

    for he=1:nhe(h)
        if node_labels === nothing #isnothing(node_labels)
            vertices = collect(keys(getvertices(h, he)))
        else
            vertices = [get(node_labels, v, v) for v in collect(keys(getvertices(h, he)))]
        end

        if edge_labels === nothing #isnothing(edge_labels)
            he_hnx = hnx.Entity(string(he), elements=vertices)
        else
            he_hnx = hnx.Entity(get(edge_labels, he, string(he)), elements=vertices)
        end

        h_hnx.add_edge(he_hnx)
    end

    h_hnx
end


text:
##############################################################################
# Results

abstract type Results end

mutable struct RegPathResults<:Results
    results    # list of PointResults
    imin       # index of lambda that minimizes test loss
end

mutable struct PointResults<:Results
    theta      # optimal theta
    lambda     # lambda used
    trainloss  # number
    testloss   # number
end

mutable struct FoldResults<:Results
    results    # list of results, one per fold
end

mutable struct NoResults<:Results end

##############################################################################
# Data

abstract type Data end

mutable struct FoldedData<:Data
    X
    Y
    nfolds
    foldrows
    nonfoldrows
    results::Results
end

# Note both X and Y are 2 dimensional Arrays
mutable struct SplitData<:Data
    Xtrain
    Ytrain
    Xtest
    Ytest
    trainrows
    testrows
    trainfrac
    results::Results
end

mutable struct NoData<:Data end