text:
module GnuTLS

using Base.Meta
using BinDeps
@BinDeps.load_dependencies [:gnutls]

import Base: isopen, write, read, readall, readavailable, close, show, nb_available

const gnutls_version = convert(VersionNumber,bytestring(ccall((:gnutls_check_version,gnutls),Ptr{Uint8},(Ptr{Uint8},),C_NULL)))

macro gnutls_since(v,f)
	if isexpr(v,:macrocall) && v.args[1] == symbol("@v_str")
		v = convert(VersionNumber,v.args[2])
	end
	if gnutls_version < v
		msg = """This function is only supported in GnuTLS versions > $v.
				 To force GnuTLS.jl to build a more recent version, you may set
				 ENV[\"GNUTLS_VERSION\"] and run Pkg.fixup()"""
		body = quote
			error($msg)
		end
		if isexpr(f,:function) || isexpr(f,:(=))
			@assert isexpr(f.args[1],:call)
			return esc(Expr(f.head,f.args[1],body))
		else
			return nothing
		end
	end
	esc(f)
end

# GnuTLS initialization

function init()
	ccall((:gnutls_global_set_mem_functions,gnutls),Void,
		(Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
		  cglobal(:jl_gc_counted_malloc),   # Malloc
		  cglobal(:jl_gc_counted_malloc),   # Secure Malloc
		  C_NULL,							# is_secure (may be NULL)
          cglobal(:realloc),  				# Realloc
          cglobal(:jl_gc_counted_free))		# Free
	ccall((:gnutls_global_init,gnutls),Int32,())
end
deinit() = ccall((:gnutls_global_deinit,gnutls),Void,())

# GnuTLS Error handling

include("errormap.jl")

immutable GnuTLSException
	msg::ASCIIString
	code::Int32
end

GnuTLSException(code::Integer) = GnuTLSException("",code)

show(io::IO,err::GnuTLSException) = print(io,"GnuTLS Exception: ",err.msg,error_codes[err.code][1],"(",string(err.code),"): ",error_codes[err.code][2])

gnutls_error(msg::ASCIIString,err::Int32) = err != 0 ? throw(GnuTLSException(err)) : nothing
gnutls_error(err::Int32) = err != 0 ? throw(GnuTLSException(err)) : nothing

# GnuTLS Session support

const GNUTLS_SERVER = 1
const GNUTLS_CLIENT = 1<<1
const GNUTLS_DATAGRAM = 1<<2
const GNUTLS_NONBLOCK = 1<<3
const GNUTLS_NO_EXTENSIONS = 1<<4

type Session <: IO
	handle::Ptr{Void}
	open::Bool
	# For rooting purposes
	read::IO
	write::IO
	function Session(isserver::Bool = false)
		x = Array(Ptr{Void},1)
		gnutls_error(ccall((:gnutls_init,gnutls),Int32,(Ptr{Ptr{Void}},Uint32),x,isserver?GNUTLS_SERVER:GNUTLS_CLIENT))
		ret = new(x[1],false)
		finalizer(ret,free_session)
		ret
	end
end

const GNUTLS_SHUT_RDWR = 0
const GNUTLS_SHUT_WR = 1

free_session(s::Session) = ccall((:gnutls_deinit,gnutls),Void,(Ptr{Void},),s.handle)
isopen(s::Session) = (isopen(s.read) || isopen(s.write)) && s.open

function close(s::Session)
	if !isopen(s.read)
		return
	end
	ret::Int32 = 0
	try # The remote might very well simply shut the stream rather than acknowledge the closure
		ret = ccall((:gnutls_bye,gnutls), Int32, (Ptr{Void},Int32), s.handle, GNUTLS_SHUT_RDWR)
	catch e
		if !isa(e,EOFError)
			rethrow()
		end
	end
	gnutls_error(ret)
	s.open=false
	nothing
end

type CertificateStore
	handle::Ptr{Void}
	function CertificateStore()
		x = Array(Ptr{Void},1)
		gnutls_error(ccall((:gnutls_certificate_allocate_credentials,gnutls),Int32,(Ptr{Ptr{Void}},),x))
		ret = new(x[1])
		finalizer(ret,free_certificate_store)
		ret
	end
end

free_certificate_store(x::CertificateStore) = ccall((:gnutls_certificate_free_credentials,gnutls),Void,(Ptr{Void},),x.handle)
@gnutls_since v"3.0" function set_system_trust!(c::CertificateStore)
	ret = ccall((:gnutls_certificate_set_x509_system_trust,gnutls),Int32,(Ptr{Void},),c.handle)
	if ret == -1250
		return false
	end
	gnutls_error(ret)
	true
end

function load_certificate(c::CertificateStore,certfile::String,keyfile::String,isPEM=false)
	gnutls_error(ccall((:gnutls_certificate_set_x509_key_file,gnutls),Int32,(Ptr{Void},Ptr{Uint8},Ptr{Uint8},Int32),c.handle,certfile,keyfile,isPEM?1:0))
end

@gnutls_since v"3.0" begin
	const system_certificate_store = CertificateStore()
	const has_system_trust = set_system_trust!(system_certificate_store)
end

@gnutls_since v"3.0" system_trust() = system_certificate_store()

function poll_readable{S<:IO}(io::S,ms::Uint32)
	c = Condition()
	@async begin
		sleep(ms/1000)
		notify(c,int32(0))
	end
	@async begin
		err = wait(readnotify(io))
		notify(c,err.code == 0 ? 1 : 0)
	end
	wait(c)::Int32
end

function read_ptr{S<:IO}(io::S,ptr::Ptr{Uint8},size::Csize_t)
	Base.wait_readnb(io,1)
	n = min(nb_available(io.buffer),size)
	read(io,pointer_to_array(ptr,int(n)))
	return signed(n)
end

function associate_stream{S<:IO,T<:IO}(s::Session, read::S, write::T=read)
	s.read = read
	s.write = write
	if write == read
		ccall((:gnutls_transport_set_ptr,gnutls),Void,(Ptr{Void},Any),s.handle,read)
	else
		ccall((:gnutls_transport_set_ptr2,gnutls),Void,(Ptr{Void},Any,Any),s.handle,read,write)
	end
	ccall((:gnutls_transport_set_pull_timeout_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(poll_readable,Int32,(S,Uint32)))
	ccall((:gnutls_transport_set_pull_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(read_ptr,Cssize_t,(S,Ptr{Uint8},Csize_t)))
	ccall((:gnutls_transport_set_push_function,gnutls),Void,(Ptr{Void},Ptr{Void}),s.handle,cfunction(Base.write,Int64,(T,Ptr{Uint8},Csize_t)))
end

handshake!(s::Session) = (gnutls_error(ccall((:gnutls_handshake,gnutls),Int32,(Ptr{Void},),s.handle));s.open = true; nothing)
function set_priority_string!(s::Session,priority::ASCIIString="NORMAL")
	x = Array(Ptr{Uint8},1)
	old_ptr = convert(Ptr{Uint8},priority)
	ret = ccall((:gnutls_priority_set_direct,gnutls),Int32,(Ptr{Void},Ptr{Uint8},Ptr{Ptr{Uint8}}),s.handle,old_ptr,x)
	offset = x[1] - old_ptr
	gnutls_error("At priority string offset $offset: ",ret)
end

const GNUTLS_CRD_CERTIFICATE = 1
const GNUTLS_CRD_ANON 		 = 2
const GNUTLS_CRD_SRP		 = 3
const GNUTLS_CRD_PSK 		 = 4
const GNUTLS_CRD_IA 		 = 5

set_credentials!(s::Session,c::CertificateStore) = gnutls_error(ccall((:gnutls_credentials_set,gnutls),Int32,(Ptr{Void},Int32,Ptr{Void}),s.handle,GNUTLS_CRD_CERTIFICATE,c.handle))

function write(io::Session, data::Ptr{Uint8}, size::Integer)
	ret = ccall((:gnutls_record_send,gnutls), Int, (Ptr{Void},Ptr{Uint8},Csize_t), io.handle, data, size)
	if ret < 0
		gnutls_error(ret)
	end
	ret
end

function write(io::Session, data::Array{Uint8,1})
	ret = ccall((:gnutls_record_send,gnutls), Int, (Ptr{Void},Ptr{Uint8},Csize_t), io.handle, data, length(data))
	if ret < 0
		gnutls_error(ret)
	end
	ret
end

function read(io::Session, data::Array{Uint8,1})
	total = 0
	while total < length(data)
		ret = ccall((:gnutls_record_recv,gnutls), Int, (Ptr{Void},Ptr{Uint8},Csize_t), io.handle, data, length(data))
		if ret < 0
			gnutls_error(ret)
		elseif ret == 0
			throw(EOFError())
		end
		total += ret
	end
	data
end

function readtobuf(io::Session,buf::IOBuffer,nb)
	Base.ensureroom(buf,int(nb))
	ret = ccall((:gnutls_record_recv,gnutls), Int, (Ptr{Void},Ptr{Uint8},Csize_t), io.handle, pointer(buf.data,buf.size+1), nb)
	if ret < 0
		gnutls_error(int32(ret))
	elseif ret == 0
		close(io)
		return false
	end
	buf.size += ret
	return true
end

const TLS_CHUNK_SIZE = 4096
function readall(io::Session)
	buf = IOBuffer(Array(Uint8,TLS_CHUNK_SIZE),true,true,true,true,typemax(Int))
	buf.size = 0
	while readtobuf(io,buf,TLS_CHUNK_SIZE); end
	readall(buf)
end

nb_available(s::Session) = ccall((:gnutls_record_check_pending,gnutls), Csize_t, (Ptr{Void},), s.handle)

function readavailable(io::Session)
	buf = IOBuffer(Array(Uint8,0),true,true,true,true,typemax(Int))
	n = nb_available(io)
	if n == 0
		readtobuf(io,buf,1)
	end
	while (n = nb_available(io))>0
		if !readtobuf(io,buf,n)
			break
		end
	end
	readall(buf)
end

function write{T}(s::Session, a::Array{T})
    if isbits(T)
        write(s,reinterpret(Uint8,a))
    else
        invoke(write, (IO, Array), s, a)
    end
end

read(io::Session, ::Type{Uint8}) = (x=Array(Uint8,1);read(io,x);x[1])

export  handshake!, associate_stream, set_priority_string!, set_credentials!

export SHA1, MD5, RMD160, SHA256, SHA384, SHA512, SHA224, HMACState, initHash, update, digest!, hexdigest!, hash
abstract HashAlgorithm

immutable SHA1 <: HashAlgorithm; end
immutable MD5 <: HashAlgorithm; end
immutable RMD160 <: HashAlgorithm; end
immutable MD2 <: HashAlgorithm; end
immutable SHA256 <: HashAlgorithm; end
immutable SHA384 <: HashAlgorithm; end
immutable SHA512 <: HashAlgorithm; end
immutable SHA224 <: HashAlgorithm; end

# Output size of the algorithm in Bytes
output_size(::Type{SHA1}) 	= 20
output_size(::Type{MD5}) 	= 16
output_size(::Type{RMD160}) = 20
output_size(::Type{MD2}) 	= 16
output_size(::Type{SHA256}) = 32
output_size(::Type{SHA384}) = 48
output_size(::Type{SHA512}) = 64
output_size(::Type{SHA224}) = 28

# Gnutls algorithm id
gnutls_id(::Type{SHA1})	 	= 3
gnutls_id(::Type{MD5}) 		= 2
gnutls_id(::Type{RMD160}) 	= 4
gnutls_id(::Type{MD2}) 		= 5
gnutls_id(::Type{SHA256}) 	= 6
gnutls_id(::Type{SHA384}) 	= 7
gnutls_id(::Type{SHA512}) 	= 8
gnutls_id(::Type{SHA224}) 	= 9

# HMAC
type HMACState{T<:HashAlgorithm}
	handle::Ptr{Void}
end

function deinit_state{T}(state::HMACState{T})
	if state.handle != C_NULL
	    ret = Array(Uint8,output_size(T));
            ccall((:gnutls_hmac_deinit,gnutls),Void,(Ptr{Void},Ptr{Uint8}),
                  state.handle,ret)
	    state.handle = C_NULL
	end
end

function HMACState{T<:HashAlgorithm}(::Type{T},key)
	x = Array(Ptr{Void},1)
	gnutls_error(ccall((:gnutls_hmac_init,gnutls),Int32,(Ptr{Ptr{Void}},Int32,Ptr{Uint8},Csize_t),x,gnutls_id(T),key,sizeof(key)))
	ret = HMACState{T}(x[1])
	finalizer(ret,deinit_state)
	ret
end

function update{T}(state::HMACState{T},data)
	state.handle == C_NULL && error("Cannot update an HMAC that was freed.")
	gnutls_error(ccall((:gnutls_hmac,gnutls),Int32,(Ptr{Void},Ptr{Uint8},Csize_t),state.handle,data,sizeof(data)))
end

function digest!{T}(state::HMACState{T})
	ret = Array(Uint8,output_size(T));
        ccall((:gnutls_hmac_output,gnutls),Void,(Ptr{Void},Ptr{Uint8}),
              state.handle,ret)
    	ret
end

# similar to Python's hmac.HMAC.hexdigest
function hexdigest!(state::HMACState)
    d = digest!(state)
    n = length(d)
    h = Array(Uint8, 2*n)
    for i = 1:n
        x = d[i]
        h[2*i] = Base.digit(x & 0xf)
        x >>= 4
        h[2*i-1] = Base.digit(x & 0xf)
    end
    ASCIIString(h)
end

# Hashing
type HashState{T<:HashAlgorithm}
	handle::Ptr{Void}
end

function initHash{T<:HashAlgorithm}(::Type{T})
	x = Array(Ptr{Void},1)
	gnutls_error(ccall((:gnutls_hash_init,gnutls),Int32,(Ptr{Ptr{Void}},Int32),x,gnutls_id(T)))
	ret = HashState{T}(x[1])
	finalizer(ret,deinit_state)
	ret
end

function update{T}(state::HashState{T},data)
	state.handle == C_NULL && error("Cannot update a Hash that was freed.")
	gnutls_error(ccall((:gnutls_hash,gnutls),Int32,(Ptr{Void},Ptr{Uint8},Csize_t),state.handle,data,sizeof(data)))
end

function digest!{T}(state::HashState{T})
	ret = Array(Uint8,output_size(T));
	ccall((:gnutls_hash_deinit,gnutls),Void,(Ptr{Void},Ptr{Uint8}),state.handle,ret)
	state.handle = C_NULL
	ret
end

hash{T<:HashAlgorithm}(::Type{T},data) = (ret = Array(Uint8,output_size(T));
	ccall((:gnutls_hash_fast,gnutls),Void,(Int32,Ptr{Uint8},Ptr{Uint8},Ptr{Uint8}),gnutls_id(T),data,sizeof(data),ret);ret)

end

GnuTLS.init()


text:
module Molly

using BioStructures
using Colors
using Distributions
using KernelDensity
using ProgressMeter
using Reexport
using Requires

@reexport using StaticArrays

using Base.Threads
using LinearAlgebra

include("types.jl")
include("setup.jl")
include("spatial.jl")
include("forces.jl")
include("simulators.jl")
include("thermostats.jl")
include("neighbours.jl")
include("loggers.jl")
include("analysis.jl")

function __init__()
    @require Makie="ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("makie.jl")
end

end

ext: jl


text:
"""
This part of the scanner constructs the initial dictionary with the contents of the XML.
"""
module Parsing

using LightXML
using ..WaylandCore

# Explicit Base imports for adding methods
import Base.string, Base.show

export wlparse, SProtocol, SInterface, SRequest, SEvent, SArgument, SEnum, SEnumEntry, SDescription

@enum RequestType none destructor
"""
	Scanner types

These corespond to XML descriptions, with comments stating how they are supplied in the XML. Some "required" by the format fields aren't actually necessary to function, so we'll asume they might be missing and just generate a warning if that's the case.
Ordering of member values:
primary:
1. attribute
2. single child node
3. multiple child nodes
secondary:
1. required
2. optional
3. implied
"""
abstract type ScannerStruct end
struct SDescription <: ScannerStruct
	summary::Union{AbstractString,Nothing} # attribute, required
	content::Union{AbstractString,Nothing} # childnode (text content of the element), required
end
struct SEnumEntry <: ScannerStruct
	name::AbstractString # attribute, required. May be redundant because the entries are accessed in the Dict by name anyway.
	value::WlUInt # attribute, required
	since::WlVersion # attribute, implied
	summary::Union{AbstractString,Nothing} # attribute, implied
	description::Union{SDescription,Nothing} # childnode, optional
end
struct SEnum <: ScannerStruct
	name::AbstractString # attribute, required
	since::WlVersion # attribute, implied
	bitfield::Bool # attribute, implied. Signifies if the enum is a plain, C-like enum, or a multi-choice bitfield
	description::Union{SDescription,Nothing} # childnode, optional
	entries::Dict{Symbol,SEnumEntry} # childnodes, required
end
struct SArgument <: ScannerStruct
	name::AbstractString # attribute, required. If the XML is non-compliant, this may be non-existent and still able to work, but the arguments to actual functions are named after it, so making it optional is TODO.
	type::TypeofWlMsgType # attribute, required
	summary::Union{AbstractString,Nothing} # attribute, implied
	interface::Union{AbstractString,Nothing} # attribute, implied. If the argument is of type "object", this describes its interface. XXX: Needs later conversion!
	nullable::Bool # attribute, implied
	enum::Union{AbstractString,Nothing} # attribute, implied. If the argument is actually an enum (has int or uint type and this attribute set), specifies which one is it. XXX: Needs later conversion!
	description::Union{SDescription,Nothing} # childnode, optional
end
struct SEvent <: ScannerStruct
	name::AbstractString # attribute, required
	since::WlVersion # attribute, implied
	description::Union{SDescription,Nothing} # childnode, optional
	args::Union{Vector{SArgument}, Nothing} # childnodes, optional
end
struct SRequest <: ScannerStruct
	name::AbstractString # attribute, required
	since::WlVersion # attribute, implied
	rtype::RequestType # attribute, implied
	description::Union{SDescription,Nothing} # childnode, optional
	args::Union{Vector{SArgument}, Nothing} # childnodes, optional
end
struct SInterface <: ScannerStruct
	name::AbstractString # attribute, required
	version::WlVersion # attribute, required
	description::Union{SDescription,Nothing} # childnode, optional
	requests::Vector{SRequest} # childnodes, optional
	events::Vector{SEvent} # childnodes, optional
	enums::Vector{SEnum} # childnodes, optional
end
struct SProtocol <: ScannerStruct
	name::AbstractString # attribute, required
	description::Union{SDescription,Nothing} # childnode, optional
	interfaces::Vector{SInterface} # childnodes, required (at least one)
end
# scanner structs helpers
"""
    imply_since(s)

Imply the "since" attribute, resolving to 1 if it's not present and valid.
"""
imply_since(::Nothing) = 1
imply_since(s::AbstractString) = (since = tryparse(WlVersion, s)) == nothing ? nothing : since
function parse_wlmsgtype(ts::AbstractString)::TypeofWlMsgType
	if ts == "int"
		WlInt
	elseif ts == "uint"
		return WlUInt
	elseif ts == "fixed"
		WlFixed
	elseif ts == "string"
		WlString
	elseif ts == "object"
		WlObjID
	elseif ts == "new_id"
		WlNewID
	elseif ts == "array"
		WlArray
	elseif ts == "fd"
		WlFD
	else
		throw(ArgumentError("Wayland message argument type \"$ts\" not recognised."))
	end
end
"""
    flatten(s::AbstractString)

Flatten a formatted string into a simple, one-line one. Strip all trailing and leading whitespace and change all interior whitespace into a single space.
"""
flatten(s::AbstractString) = replace(strip(s), r"\s+"=>" ")
"""
    string(desc::SDescription)

Converts an SDescription into a plain string.
"""
function string(desc::SDescription)::String
	if desc.summary != nothing
		if desc.content != nothing
			return "$(desc.summary)\n\n$(desc.content)"
		else
			return desc.summary
		end
	else
		if desc.content != nothing
			return desc.content
		else
			return ""
		end
	end
end
"""
	ScannerStruct(::XMLElement)

Do note that all these constructors trust the element to actually be their XML counterpart, the caller is responsible for supplying correct element nodes.
"""
function SProtocol(element::XMLElement)
	name = attribute(element, "name"; required=true)
	_description = find_element(element, "description")
	description = _description == nothing ? nothing : SDescription(_description)
	interfaces	= Vector{SInterface}()
	for child in child_elements(element)
		if LightXML.name(child) == "interface"
			push!(interfaces, SInterface(child))
		end
	end
	SProtocol(name, description, interfaces)
end
function SInterface(element::XMLElement)
	name = attribute(element, "name"; required=true)
	_version = attribute(element, "version")
	version  = _version == nothing ? 1 : parse(WlVersion, _version)
	_description = find_element(element, "description")
	description  = _description == nothing ? nothing : SDescription(_description)
	enums = Vector{SEnum}()
	requests = Vector{SRequest}()
	events = Vector{SEvent}()
	for child in child_elements(element)
		childname = LightXML.name(child)
		if childname == "enum"
			push!(enums, SEnum(child))
		elseif childname == "request"
			push!(requests, SRequest(child))
		elseif childname == "event"
			push!(events, SEvent(child))
		end
	end
	SInterface(name, version, description, requests, events, enums)
end
function SRequest(element::XMLElement)
	name = attribute(element, "name"; required=true)
	since = imply_since(attribute(element, "since"))
	_rtype = attribute(element, "type")
	if	_rtype == "destructor"
		rtype = destructor
	else
		rtype = none
	end
	_description = find_element(element, "description")
	description = _description == nothing ? nothing : SDescription(_description)
	args = Vector{SArgument}()
	for child in child_elements(element)
		if LightXML.name(child) == "arg"
			push!(args, SArgument(child))
		end
	end
	isempty(args) && (args = nothing)
	SRequest(name, since, rtype, description, args)
end
function SEvent(element::XMLElement)
	name = attribute(element, "name"; required=true)
	since = imply_since(attribute(element, "since"))
	_description = find_element(element, "description")
	description = _description == nothing ? nothing : SDescription(_description)
	args = Vector{SArgument}()
	for child in child_elements(element)
		if LightXML.name(child) == "arg"
			push!(args, SArgument(child))
		end
	end
	isempty(args) && (args = nothing)
	SEvent(name, since, description, args)
end
function SEnum(element::XMLElement)
	name = attribute(element, "name"; required=true)
	since = imply_since(attribute(element, "since"))
	_bitfield = attribute(element, "bitfield")
	if	_bitfield == nothing
		bitfield = false
	else
		bitfield = parse(Bool, _bitfield)
	end
	_description = find_element(element, "description")
	description = _description == nothing ? nothing : SDescription(_description)
	entries = Dict{Symbol,SEnumEntry}()
	for child in child_elements(element)
		if LightXML.name(child) == "entry"
			entry = SEnumEntry(child)
			entries[Symbol(entry.name)] = entry
		end
	end
	SEnum(name, since, bitfield, description, entries)
end
function SEnumEntry(element::XMLElement)
	name = attribute(element, "name"; required=true)
	value = parse(WlUInt, attribute(element, "value"; required=true))
	since = imply_since(attribute(element, "since"))
	_summary = attribute(element, "summary")
	_description = find_element(element, "description")
	if	_description == nothing
		description = nothing
		summary = _summary
	else
		description = SDescription(_description)
		if _summary == nothing
			summary = description.summary
		else
			summary = _summary
		end
	end
	SEnumEntry(name, value, since, summary, description)
end
function SArgument(element::XMLElement)
	name = attribute(element, "name"; required=true)
	type = parse_wlmsgtype(attribute(element, "type"; required=true))
	if type != WlID
		interface = nothing
	else
		interface = attribute(element, "interface")
	end
	_nullable = attribute(element, "nullable")
	if	_nullable == nothing
		nullable = false
	else
		nullable = parse(Bool, _nullable)
	end
	if type <: Integer
		enum = attribute(element, "enum")
	else
		enum = nothing
	end
	_description = find_element(element, "description")
	_summary = attribute(element, "summary")
	if	_description == nothing
		description = nothing
		summary = _summary
	else
		description = SDescription(_description)
		if _summary == nothing
			summary = description.summary
		else
			summary = _summary
		end
	end
	SArgument(name, type, summary, interface, nullable, enum, description)
end
function SDescription(element::XMLElement)
	econtent = flatten(content(element))
	econtent = econtent == "" ? nothing : econtent
	SDescription(attribute(element, "summary"), econtent)
end

# utility functions
stype(::SProtocol) = "Protocol"
stype(::SInterface) = "Interface"
stype(::SRequest) = "Request"
stype(::SEvent) = "Event"
stype(::SEnum) = "Enum"
stype(::SEnumEntry) = "Enum entry"
stype(::SArgument) = "Argument"
stype(::SDescription) = "Description"
"""
    SCollection

Any collection of `ScannerStruct`s.
"""
const SCollection = Union{Set{<: ScannerStruct},Vector{<: ScannerStruct},Dict{<: Any,<: ScannerStruct}}
"""
    getindent(count::Integer, s::AbstractString = "\t")

Get an indent string from the indentation depth (count) and a string representing an indentation level.
"""
getindent(count::Integer, s::AbstractString = "\t") = repeat(s, count)
"""
	indent1

The string representing a single indent level.
"""
indent1 = getindent(1)
"""
	show(io::IO, mime::MIME"text/plain". o::ScannerStruct)

Show a ScannerStruct in a human readable indented format.

Overriden with more specific methods when needed.
"""
function show(io::IO, mime::MIME"text/plain", o::ScannerStruct)
	indent_level = get(io, :indent, 0)
	indent = getindent(indent_level, indent1)
	print(io, "$indent$(stype(o)) $(o.name):\n")
	for fname in fieldnames(typeof(o))
		prop = getproperty(o, fname)
		if prop == nothing
			print("$indent$indent1$fname: none")
		elseif prop isa ScannerStruct
			print("$indent$indent1$fname:\n")
			show(IOContext(io, :indent=>indent_level + 2), mime, prop)
		elseif prop isa SCollection
			print(io, "$indent$indent1$fname:")
			show(IOContext(io, :indent=>indent_level + 2), mime, prop)
		else
			print(io, "$indent$indent1$fname: $(repr(MIME("text/plain"), prop))")
		end
		print(io, '\n')
	end
end
"""
    show(io::IO, mime::MIME"text/plain", o::Set{<: ScannerStruct})

Pretty-print a collection of [`ScannerStruct`](@ref)s.
"""
function show(io::IO, mime::MIME"text/plain", o::SCollection)
	if isempty(o)
		print(io, " none")
	else
		print(io, '\n')
		indent_level = get(io, :indent, 0)
		for member in o
			show(IOContext(io, :indent=>indent_level), mime, member)
		end
	end
end
"""
    show(io::IO, mime::MIME"text/plain", o::Dict{Symbol,<: ScannerStruct})

Pretty-print a [`Symbol`](@ref) [`Dict`](@ref) of [`ScannerStruct`](@ref)s.
"""
function show(io::IO, mime::MIME"text/plain", o::Dict{Symbol,<: ScannerStruct})
	indent_level = get(io, :indent, 0)
	indent = getindent(indent_level, indent1)
	print(io, '\n')
	for (key, val) in o
		print(io, "$indent$key:\n")
		show(IOContext(io, :indent=>indent_level + 1), mime, val)
	end
end
"""
	show(io::IO, ::MIME"text/plain", o::SDescription)

Show an [`SDescription`](@ref) in a human readable indented format, shortening the strings.
"""
function show(io::IO, ::MIME"text/plain", o::SDescription)
	if o.summary == nothing && o.content == nothing
		show(io, o)
	else
		indent_level = get(io, :indent, 0)
		indent = getindent(indent_level, indent1)
		summary = o.summary == nothing ? "none" : length(o.summary) >= 100 ? o.summary[1:100] * "…" : o.summary
		content = o.content == nothing ? "none" : length(o.content) >= 100 ? o.content[1:100] * "…" : o.content
		print(io, "$(indent)Description:\n$(indent)$(indent1)summary: $summary\n$(indent)$(indent1)content: $content")
	end
end

"""
    wlparse(element::XMLElement)

Parse an [`XMLElement`](@ref) into a `Set{SProtocol}`.

Returns a `Set{SProtocol}` so that when parsing multiple protocol files or even a file with multiple protocols they can all be easily merged via `append!`.

Example:
```julia
	import LightXML
	protocols = Set{SProtocol}()
	for file in filenames
		xdoc = LightXML.parse_file(file)
		finalizer(free, xdoc)
		append!(protocols, wlparse(LightXML.root(xdoc)))
	end
```
"""
function wlparse(element::XMLElement)
	protocols = Set{SProtocol}()
	if name(element) == "protocol"
		push!(protocols, SProtocol(element))
	else
		for child in child_elements(element)
			append!(wlparse(child))
		end
	end
	return protocols
end

const DEFAULT_PROTOCOL_PATH = "/usr/share/wayland/wayland.xml"

"""
    wlparse(path::AbstractString)

Parse a file denoted by `path`.

Examples:
```julia
# Parse the base protocol file
wayland_protocol = wlparse("$DEFAULT_PROTOCOL_PATH")
# Parse a collection of files
protocol_files = Set("$DEFAULT_PROTOCOL_PATH", "/usr/share/wayland-protocols/stable/xdg-shell")
protocols = reduce(append!, wlparse.(protocol_files))
```
"""
function wlparse(path::AbstractString)
	xdoc = parse_file(path)
	finalizer(free, xdoc)
	wlparse(root(xdoc))
end

"""
    wlparse()

Parse the default Wayland protocol located at `$DEFAULT_PROTOCOL_PATH`.
"""
function wlparse()
	wlparse(DEFAULT_PROTOCOL_PATH)
end

end #module Parsing


text:
using OrdinaryDiffEq, ParameterizedFunctions, DiffEqBase, Base.Test

# Here's the problem to solve

lorenz = @ode_def_bare Lorenz begin
  dx = σ*(y-x)
  dy = ρ*x-y-x*z
  dz = x*y-β*z
end σ = 10. β = 8./3. ρ => 28.

u0 = [1., 5., 10.]
tspan = (0., 100.)
prob = ODEProblem(lorenz, u0, tspan)
sol = solve(prob,Tsit5())

@test DiffEqBase.interpret_vars([(0,1), (1,3), (4,5)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, (1,3), (4,5)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, 3, 4],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,3), (DiffEqBase.DEFAULT_PLOT_FUNC,0,4)]
@test DiffEqBase.interpret_vars(([1,2,3], [4,5,6]),sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,4), (DiffEqBase.DEFAULT_PLOT_FUNC,2,5), (DiffEqBase.DEFAULT_PLOT_FUNC,3,6)]
@test DiffEqBase.interpret_vars((1, [2,3,4]),sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,1,4)]

@test DiffEqBase.interpret_vars([(:t,:x),(:t,:y)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2)]
@test DiffEqBase.interpret_vars([:x, (0,:x), (:x,:y)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars([:x, :y, :z],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2), (DiffEqBase.DEFAULT_PLOT_FUNC,0,3)]
@test DiffEqBase.interpret_vars(([:x,:x], [:y,:z]),sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3)]
@test DiffEqBase.interpret_vars((:x, [:y,:z]),sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3)]

f = (x,y) -> (x+y,y)
@test DiffEqBase.interpret_vars([(f,0,1), (1,3), (4,5)],sol) == [(f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, (f,1,3), (4,5)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (f,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([(f,:t,:x),(:t,:y)],sol) == [(f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2)]
@test DiffEqBase.interpret_vars([:x, (f,0,:x), (:x,:y)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars([(:x,:y)],sol) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars((f,:x,:y),sol) == [(f,1,2)]


text:
module QJuliaReduce

import QJuliaRegisters
import QJuliaIntrinsics

using MPI

blas_global_reduction = true

@inline function set_blas_global_reduction(flag::Bool)
  global blas_global_reduction = flag
end

@inline function reset_blas_global_reduction()
  global blas_global_reduction = true
end

#SSE
double2   = QJuliaRegisters.double2
float4    = QJuliaRegisters.float4
#AVX/AVX2
double4   = QJuliaRegisters.double4
float8    = QJuliaRegisters.float8
#AVX3
double8   = QJuliaRegisters.double8
float16   = QJuliaRegisters.float16

mm_mul  = QJuliaIntrinsics.mm_mul
mm_add  = QJuliaIntrinsics.mm_add
mm_sub  = QJuliaIntrinsics.mm_sub
mm_mad  = QJuliaIntrinsics.mm_mad

@inline function mm_dot(x::Vector{double4}, y::Vector{double4})
                 local res = double4((0.0, 0.0, 0.0, 0.0))
                 for i in 1:length(x)
                   a   = mm_mul(x[i], y[i])
                   res = mm_add(res,a)
                 end

                 send_buff = [res[1].value, res[2].value, res[3].value, res[4].value]
                 if (blas_global_reduction == true)
                   recv_buff = zeros(Cdouble, 4)
                   MPI.Allreduce!(send_buff, recv_buff, MPI.SUM, MPI.COMM_WORLD)
                 else
                   recv_buff = send_buff
                 end

                 return (recv_buff[1]+recv_buff[2]+recv_buff[3]+recv_buff[4])
end #mm_dot

@inline function gnorm2(x::Vector{T})  where T <: AbstractFloat
	         #set precision for the big float if necessary.
	         if(T == BigFloat); setprecision(T, 128); end
	         local res::T = 0.0

                 for e in x; res += e*e; end

                 if (blas_global_reduction == true)
		   resd = Float64(res)
                   val  = MPI.Allreduce(resd, MPI.SUM, MPI.COMM_WORLD)
                 else
                   val = res
                 end

                 return val
end #gnorm2

@inline function gnorm2(x::Vector{Complex{T}})  where T <: AbstractFloat
                 local res = 0.0
                 for e in x; res += abs2(e); end

                 if (blas_global_reduction == true)
                   val = MPI.Allreduce(res, MPI.SUM, MPI.COMM_WORLD)
                 else
                   val = res
                 end

                 return val
end #gnorm2

@inline function gnorm2(x::AbstractArray)
                 local res = 0.0
                 for i in 1:length(x); res += abs2(x[i]); end

                 if (blas_global_reduction == true)
                   val = MPI.Allreduce(res, MPI.SUM, MPI.COMM_WORLD)
                 else
                   val = res
                 end

                 return res
end #gnorm2

@inline function cDotProductNormX(x::Vector{T}, y::Vector{T})  where T <: AbstractFloat
                 local cres = 0.0+0.0im
                 local rres = 0.0

                 cx = view(reinterpret(Complex{T}, x), :)
                 cy = view(reinterpret(Complex{T}, y), :)

                 for i in 1:length(cx)
                    conjcx    = conj(cx[i])
                    cres += (conjcx * cy[i])
                    rres += abs2(conjcx)
                 end

                 send_buff = [real(cres), imag(cres), rres]

                 if blas_global_reduction == true
                   recv_buff = zeros(Cdouble, 3)
                   MPI.Allreduce!(send_buff, recv_buff, MPI.SUM, MPI.COMM_WORLD)
                 else
                   recv_buff = send_buff
                 end

                 return (Complex{T}(recv_buff[1], recv_buff[2]), recv_buff[3])
end #cDotProductNormX

@inline function reDotProductNormX(x::Vector{T}, y::Vector{T})  where T <: AbstractFloat

	local xdty = 0.0
	local xnrm = 0.0

	for i in 1:length(x)
		xdty += (x[i] * y[i])
		xnrm += abs2(x[i])
	end
	send_buff = [xdty, xnrm]
	if blas_global_reduction == true
		recv_buff = zeros(Cdouble, 2)
		MPI.Allreduce!(send_buff, recv_buff, MPI.SUM, MPI.COMM_WORLD)
	else
		recv_buff = send_buff
	end
	return (recv_buff[1], recv_buff[2])
end #reDotProductNormX

@inline function reDotProduct(x::Vector{T}, y::Vector{T})  where T <: AbstractFloat
	             #set precision for the big float if necessary.
	             if(T == BigFloat); setprecision(T, 128); end
                 local res::T = 0.0

                 for i in 1:length(x); res += x[i]*y[i]; end

				 #println("Info : ", sizeof(res), " res = ", res)
				 #resstrn = string(res)
				 #println("String = ", resstrn)
				 #resd1 = parse(Float64, resstrn[1:16]) # do parsing
				 #resd2 = parse(Int64, resstrn[17:32]) # do parsing
				 resd = Float64(res)

                 return (blas_global_reduction == true ?  MPI.Allreduce(resd, MPI.SUM, MPI.COMM_WORLD) : res)
end #reDotProduct

@inline function cDotProduct(x::Vector{T}, y::Vector{T})  where T <: AbstractFloat
                 local res = 0.0+0.0im

                 cx = view(reinterpret(Complex{T}, x), :)
                 cy = view(reinterpret(Complex{T}, y), :)

                 for i in 1:length(cx); res += conj(cx[i]) * cy[i]; end

                 return (blas_global_reduction == true ?  MPI.Allreduce(res, MPI.SUM, MPI.COMM_WORLD) : res)
end #reDotProduct

end #QJuliaReduce

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
module TemperatureConverter

	export as_celsius

	function as_celsius(temperature, unit)
		if unit == :Celsius
			return temperature
		elseif unit == :Kelvin
			return kelvin_to_celsius(temperature)
		end
	end

	function kelvin_to_celsius(temperature)
		# 'private' function
		return temperature + 273
	end

end


text:
using LatinSquares

if length(ARGS) == 0
	println("Usage: julia run_latin n")
	exit()
end
n = parse(ARGS[1])
println("n = $n")
println("Starting up")

# Here's a place to adjust Gurobi settings
# env = Gurobi.Env()
# setparam!(env,"Threads", 8)
# setparam!(env,"ConcurrentMIP", 2)

tic();
try
	A,B = ortho_latin(n)
	println("A = $A")
	println("B = $B")
end
toc()


text:
using Latexify
using Test
import Latexify: PlainNumberFormatter, PrintfNumberFormatter

@test FancyNumberFormatter() == FancyNumberFormatter(4) == FancyNumberFormatter("%.4g", "\\cdot") == FancyNumberFormatter("%.4g", s"\g<mantissa> \\cdot 10^{\g<sign_exp>\g<mag_exp>}")

x = -23.4728979e7

@test PlainNumberFormatter()(x) == "-2.34728979e8"
@test PrintfNumberFormatter("%.4g")(x) == "-2.347e+08"
@test StyledNumberFormatter()(x) == "-2.347 \\mathrm{e}{8}"
@test FancyNumberFormatter()(x) == "-2.347 \\cdot 10^{8}"
@test FancyNumberFormatter("%.5E", s"\g<mantissa>,\g<before_dp>,\g<sign>,\g<before_dp_nosign>,\g<after_dp>,\g<e_or_E>,\g<raw_exp>,\g<sign_exp>,\g<mag_exp>")(x) == "-2.34729,-2,-,2,34729,E,+08,,8"

xne = -23.4728979e-7

@test FancyNumberFormatter("%.5E", s"\g<mantissa>,\g<before_dp>,\g<sign>,\g<before_dp_nosign>,\g<after_dp>,\g<e_or_E>,\g<raw_exp>,\g<sign_exp>,\g<mag_exp>")(xne) == "-2.34729,-2,-,2,34729,E,-06,-,6"

y = 0xf43

@test StyledNumberFormatter()(y) == FancyNumberFormatter()(y) == "\\mathtt{0x0f43}"


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::InfiniteVariableIndex
    )::InfiniteVariableRef
    return InfiniteVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::VariableData{<:InfiniteVariable}
    )::InfiniteVariableIndex
    return MOIUC.add_item(model.infinite_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(
    model::InfiniteModel,
    ::Type{InfiniteVariable}
    )::MOIUC.CleverDict{InfiniteVariableIndex, VariableData{<:InfiniteVariable}}
    return model.infinite_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    vref::InfiniteVariableRef
    )::MOIUC.CleverDict{InfiniteVariableIndex, VariableData{<:InfiniteVariable}}
    return JuMP.owner_model(vref).infinite_vars
end

# Extend _data_object
function _data_object(vref::InfiniteVariableRef)
    object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
    if isnothing(object)
        error("Invalid infinite variable reference, cannot find ",
              "corresponding variable in the model. This is likely ",
              "caused by using the reference of a deleted variable.")
    end
    return object
end

# Extend _core_variable_object
function _core_variable_object(vref::InfiniteVariableRef)
    return _data_object(vref).variable
end

# Extend _object_numbers
function _object_numbers(vref::InfiniteVariableRef)::Vector{Int}
    return _core_variable_object(vref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(vref::InfiniteVariableRef)::Vector{Int}
    return _core_variable_object(vref).parameter_nums
end

# Define getter function for var.is_vector_start
function _is_vector_start(vref::InfiniteVariableRef)::Bool
    return _core_variable_object(vref).is_vector_start
end

## Set helper methods for adapting data_objects with parametric changes
# No change needed
function _adaptive_data_update(
    vref::InfiniteVariableRef,
    var::V,
    data::VariableData{V}
    )::Nothing where {V <: InfiniteVariable}
    data.variable = var
    return
end

# Reconstruction is necessary
function _adaptive_data_update(
    vref::InfiniteVariableRef,
    var::V1,
    data::VariableData{V2}
    )::Nothing  where {V1, V2}
    new_data = VariableData(var, data.name, data.lower_bound_index,
                            data.upper_bound_index, data.fix_index,
                            data.zero_one_index, data.integrality_index,
                            data.measure_indices, data.constraint_indices,
                            data.in_objective, data.point_var_indices,
                            data.semi_infinite_var_indices,
                            data.derivative_indices, data.deriv_constr_indices)
    _data_dictionary(vref)[JuMP.index(vref)] = new_data
    return
end

# Extend _set_core_variable_object for InfiniteVariableRefs
function _set_core_variable_object(
    vref::InfiniteVariableRef,
    var::InfiniteVariable
    )::Nothing
    _adaptive_data_update(vref, var, _data_object(vref))
    return
end

################################################################################
#                          DEFINTION HELPER METHODS
################################################################################
"""
    Infinite{VT <: VectorTuple} <: InfOptVariableType

A `DataType` to assist in making infinite variables. This can be passed as an
extra argument to `@variable` to make an infinite variable:
```julia
@variable(model, var_expr, Infinite(parameter_refs...), args..., kwargs...)
```
Here `parameter_refs` can be a single parameter reference, a single parameter
array with parameters defined in the same macro call, or multiple arguments where
each argument is either of the first two options listed.

**Fields**
- `parameter_refs::VT`: The infinite parameters the variable will depend on.
"""
struct Infinite{VT <: Collections.VectorTuple} <: InfOptVariableType
    parameter_refs::VT
    function Infinite(vt::VT) where {VT <: Collections.VectorTuple}
        return new{VT}(vt)
    end
end
function Infinite(args...)
    return Infinite(Collections.VectorTuple(args))
end

## Check that each parameter tuple element is formatted correctly
# IndependentParameterRefs
function _check_tuple_element(
    _error::Function,
    prefs::Vector{IndependentParameterRef}
    )::Nothing
    return
end

# DependentParameterRefs
function _check_tuple_element(
    _error::Function,
    prefs::Vector{DependentParameterRef}
    )::Nothing
    if length(prefs) != _num_parameters(first(prefs))
        _error("Infinite parameter tuple elements cannot depend on a subset ",
               "of dependent parameters.")
    end
    return
end

# Fallback
function _check_tuple_element(_error::Function, prefs)
    _error("Cannot have mixed parameter types in a tuple element and can only ",
           "specify infinite parameters.")
end

# Check parameter tuple, ensure all elements contain parameter references
function _check_parameter_tuple(
    _error::Function,
    raw_prefs::Collections.VectorTuple{GeneralVariableRef}
    )::Nothing
    allunique(raw_prefs) || _error("Cannot double specify infinite parameter ",
                                   "references.")
    for i in 1:size(raw_prefs, 1)
        prefs = dispatch_variable_ref.(raw_prefs[i, :])
        _check_tuple_element(_error, prefs)
    end
    return
end
function _check_parameter_tuple(
    _error::Function,
    raw_prefs::Collections.VectorTuple{T}
    ) where {T}
    _error("Expected infinite parameter reference arguments, but got arguments ",
           "of type $T. This may be because `Infinite()` was given (i.e., no ",
           "infinite parameter references were specified as arguments).")
end

## Check and format the variable info considering functional start values
# Just a number given for the start value
function _check_and_format_infinite_info(
    _error::Function,
    info::JuMP.VariableInfo{<:Real, <:Real, <:Real, <:Real},
    prefs::Collections.VectorTuple
    )
    # prepare the start value function and return the info
    start_func = (s::Vector{<:Real}) -> info.start
    return JuMP.VariableInfo(info.has_lb, convert(Float64, info.lower_bound),
                             info.has_ub, convert(Float64, info.upper_bound),
                             info.has_fix, convert(Float64, info.fixed_value),
                             !isnan(info.start), start_func, info.binary,
                             info.integer),
           true
end

# make function for checking start value function properties
function _check_valid_function(
    _error::Function,
    func::Function,
    prefs::Collections.VectorTuple
    )::Nothing
    input_format = typeof(Tuple(Vector{Float64}(undef, length(prefs)), prefs))
    if !hasmethod(func, input_format)
        _error("Specified function `$func` must be able to accept a `Float64` " *
               "support realization of the infinite parameter tuple $(prefs).")
    end
    return
end

# A function is given for the start value generation
function _check_and_format_infinite_info(
    _error::Function,
    info::JuMP.VariableInfo{<:Real, <:Real, <:Real, F},
    prefs::Collections.VectorTuple
    )::Tuple{JuMP.VariableInfo{Float64, Float64, Float64, F}, Bool} where {F <: Function}
    # check the function properties
    _check_valid_function(_error, info.start, prefs)
    # make the info and return
    return JuMP.VariableInfo(info.has_lb, convert(Float64, info.lower_bound),
                             info.has_ub, convert(Float64, info.upper_bound),
                             info.has_fix, convert(Float64, info.fixed_value),
                             info.has_start, info.start, info.binary,
                             info.integer),
           false
end

# Fallback
function _check_and_format_infinite_info(
    _error::Function,
    info::JuMP.VariableInfo,
    prefs::Collections.VectorTuple
    )
    _error("Unrecognized formatting for the variable information.")
end

"""
    JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                        var_type::Infinite)::InfiniteVariable{GeneralVariableRef}

Build and return an infinite variable based on `info` and `var_type`. Errors if
the infinite parameter references included in `var_type` are invalid. See
[`Infinite`](@ref) for more information.

**Example**
```julia-repl
julia> info = VariableInfo(false, 0, false, 0, false, 0, true, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite(t));
```
"""
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    var_type::Infinite;
    extra_kwargs...
    )
    # check for unneeded keywords
    for (kwarg, _) in extra_kwargs
        _error("Keyword argument $kwarg is not for use with infinite variables.")
    end
    # get the parameter_refs
    prefs = var_type.parameter_refs
    # check the VectorTuple for validity and format
    _check_parameter_tuple(_error, prefs)
    # check and format the info (accounting for start value functions)
    new_info, is_vect_func = _check_and_format_infinite_info(_error, info, prefs)
    # get the parameter object numbers
    object_nums = Int[]
    for pref in prefs
        union!(object_nums, _object_number(pref))
    end
    # make the variable and return
    return InfiniteVariable(new_info, prefs,
                            [_parameter_number(pref) for pref in prefs],
                            object_nums, is_vect_func)
end

# check the pref tuple contains only valid parameters
function _check_parameters_valid(
    model::InfiniteModel,
    prefs::Collections.VectorTuple
    )::Nothing
    for pref in prefs
        JuMP.check_belongs_to_model(pref, model)
    end
    return
end

# Used to update parameter-infinite variable mappings
function _update_param_var_mapping(
    vref::InfiniteVariableRef,
    prefs::Collections.VectorTuple
    )::Nothing
    for pref in prefs
        dependency_list = _infinite_variable_dependencies(pref)
        if !(JuMP.index(vref) in dependency_list)
            push!(dependency_list, JuMP.index(vref))
        end
    end
    return
end

"""
    JuMP.add_variable(model::InfiniteModel, var::InfiniteVariable,
                      [name::String = ""])::GeneralVariableRef

Extend the `JuMP.add_variable` function to accomodate infinite variable
types. Adds a variable to an infinite model `model` and returns a
[`GeneralVariableRef`](@ref). Primarily intended to be an internal function of
the constructor macro `@variable`. However, it can be used in combination with
`JuMP.build_variable` to add  infinite variables to an infinite model object.
Errors if invalid parameter reference(s) are included in `var`.

**Example**
```julia-repl
julia> @infinite_parameter(m, t in [0, 10]);

julia> info = VariableInfo(false, 0, false, 0, false, 0, true, 0, false, false);

julia> inf_var = build_variable(error, info, Infinite(t));

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)
```
"""
function JuMP.add_variable(
    model::InfiniteModel,
    v::InfiniteVariable,
    name::String = ""
    )::GeneralVariableRef
    _check_parameters_valid(model, v.parameter_refs)
    data_object = VariableData(v, name)
    vindex = _add_data_object(model, data_object)
    vref = InfiniteVariableRef(model, vindex)
    _update_param_var_mapping(vref, v.parameter_refs)
    gvref = _make_variable_ref(model, vindex)
    _set_info_constraints(v.info, gvref, vref)
    model.name_to_var = nothing
    return gvref
end

################################################################################
#                              RESTRICTION METHODS
################################################################################
## Dispatch functions for functional syntax of making restricted variables
# Point variable
function _restrict_infinite_variable(
    ivref::GeneralVariableRef,
    vt::Collections.VectorTuple{<:Real}
    )::GeneralVariableRef
    info = JuMP.VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN,
                             false, false)
    new_var = JuMP.build_variable(error, info, Point(ivref, vt))
    return JuMP.add_variable(JuMP.owner_model(ivref), new_var,
                             update_info = false)
end

# Infinite variable
function _restrict_infinite_variable(
    ivref::GeneralVariableRef,
    vt::Collections.VectorTuple{<:GeneralVariableRef}
    )::GeneralVariableRef
    if parameter_list(ivref) != vt.values
        error("Unrecognized syntax for infinite variable restriction.")
    end
    @warn("Unnecessary use of functional infinite variable restriction syntax " *
          "that will cause performance degredations. This was probably caused " *
          "by using syntax like `y(t, x)` inside expressions. Instead just " *
          "use the infinite variable reference (e.g. `y`).")
    return ivref
end

# Semi-Infinite variable
function _restrict_infinite_variable(
    ivref::GeneralVariableRef,
    vt::Collections.VectorTuple
    )::GeneralVariableRef
    info = JuMP.VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN,
                             false, false)
    new_var = JuMP.build_variable(error, info, SemiInfinite(ivref, vt))
    return JuMP.add_variable(JuMP.owner_model(ivref), new_var)
end

"""
    restrict(ivref::GeneralVariableRef, supps...)::GeneralVariableRef

Restrict the input domain of an infinite variable/derivative `ivref` in
accordance with the infinite parameters and/or values `supps`. Here `supps` must
match the formatting of `ivref`'s infinite parameters. Here the following
outputs are possible:
- Equivalent to `@variable(model, variable_type = Point(ivref, supps...)` if
  `supps` are a complete support point
- Equivalent to `@variable(model, variable_type = SemiInfinite(ivref, supps...)`
  if `supps` are a partial support point.

Conveniently, we can also invoke this method by calling `ivref(supps...)`.

Errors if ivref is not an infinite variable or derivative or the formatting of
`supps` is incorrect. Will warn if `supps` only contain infinite parameters and
will simply return `ivref`.

**Example**
```julia-repl
julia> restrict(y, 0, x)
y(0, [x[1], x[2]])

julia> restrict(y, 0, [0, 0])
y(0, [0, 0])

julia> y(0, x)
y(0, [x[1], x[2]])

julia> y(0, [0, 0])
y(0, [0, 0])
```
"""
function restrict(ivref::GeneralVariableRef, supps...)::GeneralVariableRef
    idx_type = _index_type(ivref)
    if idx_type != InfiniteVariableIndex && idx_type != DerivativeIndex
        error("The `vref(values..)` restriction syntax is only valid for ",
              "infinite variable and derivative references.")
    end
    return _restrict_infinite_variable(ivref, Collections.VectorTuple(supps))
end

# This enables functional calls (e.g., `y(0, x)` => semi-infinite variable)
function _functional_reference_call(
    ivref::GeneralVariableRef,
    ::Union{Type{InfiniteVariableIndex}, Type{DerivativeIndex}},
    supps...
    )::GeneralVariableRef
    return _restrict_infinite_variable(ivref, Collections.VectorTuple(supps))
end

################################################################################
#                            VARIABLE DEPENDENCIES
################################################################################
# Extend _semi_infinite_variable_dependencies
function _semi_infinite_variable_dependencies(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Vector{SemiInfiniteVariableIndex}
    return _data_object(vref).semi_infinite_var_indices
end

# Extend _point_variable_dependencies
function _point_variable_dependencies(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Vector{PointVariableIndex}
    return _data_object(vref).point_var_indices
end

# Extend _derivative_dependencies
function _derivative_dependencies(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Vector{DerivativeIndex}
    return _data_object(vref).derivative_indices
end

"""
    used_by_semi_infinite_variable(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool

Return a `Bool` indicating if `vref` is used by a semi-infinite variable.

**Example**
```julia-repl
julia> used_by_semi_infinite_variable(vref)
false
```
"""
function used_by_semi_infinite_variable(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Bool
    return !isempty(_semi_infinite_variable_dependencies(vref))
end

"""
    used_by_point_variable(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool

Return a `Bool` indicating if `vref` is used by a point variable.

**Example**
```julia-repl
julia> used_by_point_variable(vref)
false
```
"""
function used_by_point_variable(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Bool
    return !isempty(_point_variable_dependencies(vref))
end

"""
    used_by_derivative(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool

Return a `Bool` indicating if `vref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool
    return !isempty(_derivative_dependencies(vref))
end

"""
    is_used(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```julia-repl
julia> is_used(vref)
false
```
"""
function is_used(vref::Union{InfiniteVariableRef, DerivativeRef})::Bool
    if used_by_measure(vref) || used_by_constraint(vref)
        return true
    end
    for vindex in _point_variable_dependencies(vref)
        if is_used(PointVariableRef(JuMP.owner_model(vref), vindex))
            return true
        end
    end
    for vindex in _semi_infinite_variable_dependencies(vref)
        if is_used(SemiInfiniteVariableRef(JuMP.owner_model(vref), vindex))
            return true
        end
    end
    for dindex in _derivative_dependencies(vref)
        if is_used(DerivativeRef(JuMP.owner_model(vref), dindex))
            return true
        end
    end
    return false
end

################################################################################
#                           PARAMETER REFERENCES
################################################################################
"""
    raw_parameter_refs(vref::InfiniteVariableRef)::VectorTuple

Return the raw [`VectorTuple`](@ref InfiniteOpt.Collections.VectorTuple) of the
parameter references that `vref` depends on. This is primarily an internal method
where [`parameter_refs`](@ref parameter_refs(vref::InfiniteVariableRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(vref::InfiniteVariableRef)
    return _core_variable_object(vref).parameter_refs
end

"""
    parameter_refs(vref::InfiniteVariableRef)::Tuple

Return the parameter references associated with the infinite variable `vref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `vref`.

**Example**
```julia-repl
julia> @variable(model, T, Infinite(t))
T(t)

julia> parameter_refs(T)
(t,)
```
"""
function parameter_refs(vref::InfiniteVariableRef)
    return Tuple(raw_parameter_refs(vref))
end

"""
    parameter_list(vref::InfiniteVariableRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `vref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(vref::InfiniteVariableRef))
is intended as the preferred user function.
"""
function parameter_list(vref::InfiniteVariableRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(vref).values
end

# get parameter list from raw VectorTuple
function parameter_list(
    prefs::Collections.VectorTuple{GeneralVariableRef}
    )::Vector{GeneralVariableRef}
    return prefs.values
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
## format infinite info
# Good to go
function _format_infinite_info(
    info::I
    )::I where {I <: JuMP.VariableInfo{Float64, Float64, Float64, <:Function}}
    return info
end

# Convert as needed
function _format_infinite_info(
    info::JuMP.VariableInfo{V, W, T, F}
    )::JuMP.VariableInfo{Float64, Float64, Float64, F} where {V, W, T, F <: Function}
    return JuMP.VariableInfo(info.has_lb, convert(Float64, info.lower_bound),
                             info.has_ub, convert(Float64, info.upper_bound),
                             info.has_fix, convert(Float64, info.fixed_value),
                             info.has_start, info.start, info.binary,
                             info.integer)
end

# Set info for infinite variables
function _update_variable_info(
    vref::InfiniteVariableRef,
    info::JuMP.VariableInfo
    )::Nothing
    prefs = raw_parameter_refs(vref)
    param_nums = _parameter_numbers(vref)
    obj_nums = _object_numbers(vref)
    is_vect_func = _is_vector_start(vref)
    new_var = InfiniteVariable(_format_infinite_info(info), prefs, param_nums,
                               obj_nums, is_vect_func)
    _set_core_variable_object(vref, new_var)
    return
end

# Specify start_value fallback for infinite variables
function JuMP.start_value(vref::Union{InfiniteVariableRef, DerivativeRef})
    error("`start_value` not defined for infinite variables, consider calling " *
          "`start_value_function` instead.")
end

# Specify set_start_value fallback for infinite variables
function JuMP.set_start_value(
    vref::Union{InfiniteVariableRef, DerivativeRef},
    value::Real
    )
    error("`set_start_value` not defined for infinite variables, consider calling " *
          "`set_start_value_function` instead.")
end

"""
    start_value_function(vref::Union{InfiniteVariableRef, DerivativeRef})::Union{Nothing, Function}

Return the function that is used to generate the start values of `vref` for
particular support values. Returns `nothing` if no start behavior has been
specified.

**Example**
```julia-repl
julia> start_value_function(vref)
my_start_func
```
"""
function start_value_function(
    vref::Union{InfiniteVariableRef, DerivativeRef}
    )::Union{Nothing, Function}
    if _variable_info(vref).has_start
        return _variable_info(vref).start
    else
        return
    end
end

"""
    set_start_value_function(vref::InfiniteVariableRef,
                             start::Union{Real, Function})::Nothing

Set the start value function of `vref`. If `start::Real` then a function is
generated to such that the start value will be `start` for the entire infinite
domain. If `start::Function` then this function should map to a scalar start value
given a support value arguments matching the format of the parameter elements in
`parameter_refs(vref)`.

**Example**
```julia-repl
julia> set_start_value_function(vref, 1) # all start values will be 1

julia> set_start_value_function(vref, my_func) # each value will be made via my_func
```
"""
function set_start_value_function(
    vref::InfiniteVariableRef,
    start::Union{Real, Function}
    )::Nothing
    info = _variable_info(vref)
    set_optimizer_model_ready(JuMP.owner_model(vref), false)
    prefs = raw_parameter_refs(vref)
    temp_info = JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                                  info.upper_bound, info.has_fix, info.fixed_value,
                                  true, start, info.binary, info.integer)
    new_info, is_vect_func = _check_and_format_infinite_info(error, temp_info, prefs)
    obj_nums = _object_numbers(vref)
    param_nums = _parameter_numbers(vref)
    new_var = InfiniteVariable(new_info, prefs, param_nums, obj_nums, is_vect_func)
    _set_core_variable_object(vref, new_var)
    # TODO update point variable start values as appropriate
    return
end

"""
    reset_start_value_function(vref::InfiniteVariableRef)::Nothing

Remove the existing start value function and return to the default. Generally,
this is triggered by deleting an infinite parameter that `vref` depends on.

**Example**
```julia-repl
julia> reset_start_value_function(vref)
```
"""
function reset_start_value_function(vref::InfiniteVariableRef)::Nothing
    info = _variable_info(vref)
    set_optimizer_model_ready(JuMP.owner_model(vref), false)
    start_func = (s::Vector{<:Real}) -> NaN
    new_info = JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                                 info.upper_bound, info.has_fix, info.fixed_value,
                                 false, start_func, info.binary, info.integer)
    prefs = raw_parameter_refs(vref)
    obj_nums = _object_numbers(vref)
    param_nums = _parameter_numbers(vref)
    new_var = InfiniteVariable(new_info, prefs, param_nums, obj_nums, true)
    _set_core_variable_object(vref, new_var)
    # TODO update point variable start values as appropriate
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::InfiniteVariableRef)::Nothing
    # remove variable info constraints associated with vref
    _delete_info_constraints(vref)
    # update parameter mapping
    all_prefs = parameter_list(vref)
    for pref in all_prefs
        filter!(e -> e != JuMP.index(vref), _infinite_variable_dependencies(pref))
    end
    model = JuMP.owner_model(vref)
    # delete associated point variables and mapping
    for index in _point_variable_dependencies(vref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated semi-infinite variables and mapping
    for index in _semi_infinite_variable_dependencies(vref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(vref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    return
end


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Arpack
using StatProfilerHTML
using BenchmarkTools

@testset "full_hbuild" begin
    atoms = []

    clusters = []
if true
    push!(atoms,Atom(1,"H",[0,0,0]))
    push!(atoms,Atom(2,"H",[0,1,0]))
    push!(atoms,Atom(3,"H",[0,0,2]))
    push!(atoms,Atom(4,"H",[0,0,3]))
    push!(atoms,Atom(5,"H",[0,0.4,4]))
    push!(atoms,Atom(6,"H",[0,0,5]))
    #push!(atoms,Atom(7,"H",[0,0,6]))
    #push!(atoms,Atom(8,"H",[0,0,7]))
    #push!(atoms,Atom(9,"H",[0,0,8]))
    #push!(atoms,Atom(10,"H",[0,0,9]))
    #clusters    = [(1:4),(5:6),(7:10)]
    clusters    = [(1:2),(3:4),(5:6),(7:8)]
    clusters    = [(1:4),(5:8)]
    clusters    = [(1:4),(5:6),(7:8)]
    clusters    = [(1:6),(7:10)]
    clusters    = [(1:2),(3:4),(5:6)]
elseif true
    push!(atoms,Atom(1,"H",[-1.30,00,0.00]))
    push!(atoms,Atom(2,"H",[-1.30,00,1.00]))
    push!(atoms,Atom(3,"H",[ 0.00,30,0.00]))
    push!(atoms,Atom(4,"H",[ 0.00,30,1.00]))
    push!(atoms,Atom(5,"H",[ 1.33,60,0.00]))
    push!(atoms,Atom(6,"H",[ 1.30,60,1.00]))
    clusters    = [(1:2),(3:4),(5:6)]
end

#    push!(atoms,Atom(1,"H",[0,0,0]))
#    push!(atoms,Atom(2,"H",[0,0,1]))
#    push!(atoms,Atom(3,"H",[0,0,2]))
#    push!(atoms,Atom(4,"H",[0,0,3]))
#    push!(atoms,Atom(5,"H",[0,0,4]))
#    push!(atoms,Atom(6,"H",[0,0,5]))
#    push!(atoms,Atom(7,"H",[0,0,6]))
#    push!(atoms,Atom(8,"H",[0,0,7]))
    #push!(atoms,Atom(9,"H",[0,0,8]))
    #push!(atoms,Atom(10,"H",[0,0,9]))
    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    na = 3
    nb = 3

    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));
    e_fci, d1_fci, d2_fci = FermiCG.pyscf_fci(ints,na,nb,conv_tol=1e-10,max_cycle=100)
    @printf(" FCI Energy: %12.8f\n", e_fci)

    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Build Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    max_roots = 400
    nroots = 10

    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]

    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots)

    display.(clusters)
    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)

    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    ci_vector = FermiCG.ClusteredState(clusters, nroots = nroots)

    FermiCG.expand_to_full_space!(ci_vector, cluster_bases, na, nb)

    display(ci_vector)
    #display(cluster_bases[2][(2,2)])

    #@time H = FermiCG.build_full_H_serial(ci_vector, cluster_ops, clustered_ham)
    #e,v = Arpack.eigs(H, nev = 1, which=:SR)
    #@printf(" Energy: %18.12f\n",real(e[1]))

    #@profilehtml H = FermiCG.build_full_H(ci_vector, cluster_ops, clustered_ham)
    @time H = FermiCG.build_full_H(ci_vector, cluster_ops, clustered_ham)
    #@btime FermiCG.build_full_H(ci_vector, cluster_ops, clustered_ham)
    e,v = Arpack.eigs(H, nev = nroots, which=:SR)
    for ei in e
        @printf(" Energy: %18.12f\n",real(ei))
    end

    display(size(v))
    display(size(ci_vector))
    FermiCG.set_vector!(ci_vector, v)

    println(" Overlaps")
    for i in 1:nroots
        for j in i:nroots
            d = FermiCG.dot(ci_vector, ci_vector, i, j)
            @printf(" %4i,%4i = %18.12f\n", i, j, d)

        end
    end

    sig1 = FermiCG.open_matvec(ci_vector, cluster_ops, clustered_ham, nbody=3, thresh=-1)
    println(" root1")
    display(ci_vector,root=1)
    println(" root2")
    display(ci_vector,root=2)
    #display(ci_vector,thresh=-1)
    #display(sig1,thresh=-1)

    sig2 = deepcopy(sig1)
    FermiCG.set_vector!(sig2, H*v)
    println(" <i|H|j>")
    for i in 1:nroots
        for j in i:nroots
            d = FermiCG.dot(sig2, ci_vector, i, j)
            @printf(" %4i,%4i = %18.12f\n", i, j, d)

        end
    end

    println(" <i|H|j>")
    for i in 1:nroots
        for j in i:nroots
            d = FermiCG.dot(sig1, ci_vector, i, j)
            @printf(" %4i,%4i = %18.12f\n", i, j, d)

        end
    end

    display(size(sig1))
    display(size(sig2))

    #FermiCG.set_vector!(ci_vector, F.vectors[:,1])
    #display(ci_vector)
    println()

    #@test isapprox(F.values[1], -5.066833300762457, atol=1e-10)

    #maximum(abs.(H-H')) < 1e-14 || error("Hamiltonian not symmetric: ",maximum(abs.(H-H')));
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Solver

A solver for geostatistical problems.
"""
abstract type Solver end

"""
    EstimationSolver

A solver for a geostatistical estimation problem.
"""
abstract type EstimationSolver <: Solver end

"""
    SimulationSolver

A solver for a geostatistical simulation problem.
"""
abstract type SimulationSolver <: Solver end

"""
    LearningSolver

A solver for a geostatistical learning problem.
"""
abstract type LearningSolver <: Solver end

"""
    solve(problem, solver; [options])

Solve the `problem` with the `solver`, optionally
passing `options`.
"""
function solve end

"""
    solve(problem, solver; procs=[myid()])

Solve the simulation `problem` with the simulation `solver`,
optionally using multiple processes `procs`.

### Notes

Default implementation calls `solvesingle` in parallel.
"""
function solve(problem::SimulationProblem, solver::SimulationSolver; procs=[myid()])
  # sanity checks
  @assert targets(solver) ⊆ name.(variables(problem)) "invalid variables in solver"

  # dictionary with variable types
  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

  # optional preprocessing
  preproc = preprocess(problem, solver)

  # pool of worker processes
  pool = CachingPool(procs)

  # list of covariables
  allcovars = covariables(problem, solver)

  # simulation loop
  results = []
  for covars in allcovars
    # simulate covariables
    reals = pmap(pool, 1:nreals(problem)) do _
      solvesingle(problem, covars, solver, preproc)
    end

    # rearrange realizations
    vnames = covars.names
    vtypes = [mactypeof[var] for var in vnames]
    vvects = [Vector{V}[] for V in vtypes]
    rtuple = (; zip(vnames, vvects)...)
    for real in reals
      for var in vnames
        push!(rtuple[var], real[var])
      end
    end

    push!(results, rtuple)
  end

  # merge results into a single dictionary
  pdomain = domain(problem)
  preals  = reduce(merge, results)

  Ensemble(pdomain, preals)
end

"""
    preprocess(problem, solver)

Preprocess the simulation `problem` once before generating each realization
with simulation `solver`.

### Notes

The output of the function is defined by the solver developer.
Default implementation returns nothing.
"""
preprocess(::SimulationProblem, ::SimulationSolver) = nothing

"""
    solvesingle(problem, covariables, solver, preproc)

Solve a single realization of `covariables` in the simulation `problem`
with the simulation `solver`, optionally using preprocessed data in `preproc`.

### Notes

By implementing this function instead of `solve`, the developer is
informing the framework that realizations generated with his/her
solver are indenpendent one from another. GeoStats.jl will trigger
the algorithm in parallel (if enough processes are available).
"""
function solvesingle end

"""
    covariables(var, solver)

Return covariables associated with the variable `var`
in the `solver` and their respective parameters.
"""
function covariables end

"""
    covariables(problem, solver)

Return all covariables in the `solver` based on list of
variables in the `problem`.
"""
function covariables(problem::Problem, solver::Solver)
  pvars = collect(name.(variables(problem)))

  result = []
  while !isempty(pvars)
    # choose a variable from the problem
    var = first(pvars)

    # find covariables of the variable
    covars = covariables(var, solver)

    # save covariables to result
    push!(result, covars)

    # update remaining variables
    pvars = setdiff(pvars, covars.names)
  end

  result
end

"""
    targets(solver)

Return target variables in the `solver`.
"""
function targets end


text:
export NeighborInfo, Neighborhood, PackedNeighborhoods
export NearNeighborDataStructure, MetricNNDS, NullNNDS, MemoizedNNDS
export NeighborhoodCache, StreamingNC, MemoizedNC, OfflineNC
export includes_controls, neighbor_info_type, neighbors

# Utilities
const F_or_B = Union{Val{:F},Val{:B}}
const BoolVal = Union{Val{true},Val{false}}
@inline always_false(::SampleIndex) = false    # used as a filtering predicate
@inline keep_controls_if(nt::NamedTuple, cond::BoolVal) = cond === Val(true) ? nt : (nt..., controls=nothing)
# @inline lattice_neighbors(v::SVector{D}) where {D} = (setindex(v, v[i] + j, i) for i in 1:D for j in (-1, 1))
@inline lattice_neighbors(v::SVector{D}) where {D} = (v + SVector(Tuple(I)) - SVector(ntuple(i -> 2, D)) for
                                                      I in CartesianIndices(ntuple(i -> 3, D)) if
                                                      I != CartesianIndex(ntuple(i -> 2, D)))        # JuliaLang/julia#29440
@inline linear_index(i::Int) = i
@inline linear_index(i::TiledIndex) = i.off_lattice
@inline collect_if_not_Array(x::Array) = x
@inline collect_if_not_Array(x) = collect(x)
@inline includes_controls(x) = includes_controls(typeof(x))
@inline neighbor_info_type(x) = neighbor_info_type(typeof(x))

# Neighborhoods
const NeighborInfo{X<:SampleIndex,D,U} = NamedTuple{(:index,:cost,:controls),Tuple{X,D,U}}
const Neighborhood{X<:SampleIndex,D,U} = AbstractVector{NeighborInfo{X,D,U}}
includes_controls( ::Type{  NeighborInfo{X,D,U}}) where {X,D,U} = U !== Nothing
includes_controls( ::Type{<:Neighborhood{X,D,U}}) where {X,D,U} = U !== Nothing
neighbor_info_type(::Type{<:Neighborhood{X,D,U}}) where {X,D,U} = NeighborInfo{X,D,U}
# indextype(   ::Type{NeighborInfo{X,D,U}}) where {X,D,U} = X
# costtype(    ::Type{NeighborInfo{X,D,U}}) where {X,D,U} = D
# controlstype(::Type{NeighborInfo{X,D,U}}) where {X,D,U} = U
function _neighbor_info_type(::Type{X}, ::Type{NamedTuple{M,Tuple{D,U}}}, include_controls::BoolVal) where {X,M,D,U}
    include_controls === Val(true) ? NeighborInfo{X,D,U} : NeighborInfo{X,D,Nothing}
end
function neighbor_info_type(nodes::SampleSet, bvp::SteeringBVP; include_controls::BoolVal=Val(false))
    X = indextype(nodes)
    BVP = typeof(bvp(nodes.init, nodes.init))
    _neighbor_info_type(X, BVP, include_controls)
end
function shift_lattice_component(n::NeighborInfo{<:TiledIndex}, shift)
    (index=shift_lattice_component(n.index, shift), cost=n.cost, controls=n.controls)
end
function reverse_neighbor_dict(i, neighborhood, ::Type{NeighborInfo{X,D,U}}) where {X,D,U}
    Dict{X,NeighborInfo{X,D,U}}(n.index => (index=i, cost=n.cost, controls=n.controls) for n in neighborhood)
end

## CSC-style Neighborhoods
struct PackedNeighborhoods{X,D,U,N<:Neighborhood{X,D,U}} <: AbstractVector{N}
    ptr::Vector{Int}
    val::Vector{NeighborInfo{X,D,U}}
    neighborhoods::Vector{N}
end
function PackedNeighborhoods(neighborhoods::AbstractVector{<:Neighborhood})
    ptr = pushfirst!(cumsum([length(n) for n in neighborhoods]) + 1, 1)
    val = reduce(vcat, neighborhoods)
    PackedNeighborhoods(ptr, val, [view(val, ptr[i]:ptr[i+1]-1) for i in 1:length(neighborhoods)])
end
Base.length(PN::PackedNeighborhoods) = length(PN.neighborhoods)
Base.getindex(PN::PackedNeighborhoods, i::Int) = PN.neighborhoods[i]
includes_controls( ::Type{<:PackedNeighborhoods{X,D,U}}) where {X,D,U} = U !== Nothing
neighbor_info_type(::Type{<:PackedNeighborhoods{X,D,U}}) where {X,D,U} = NeighborInfo{X,D,U}

# Near Neighborhood Caches
abstract type NeighborhoodCache end
includes_controls(::Type{<:NeighborhoodCache}) = false
## StreamingNC
struct StreamingNC <: NeighborhoodCache end
function neighbors(cache::StreamingNC, nnds, nodes, bvp, v, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false))
    neighbors(nnds, nodes, bvp, v, r, dir=dir, include_controls=include_controls)
end

## MemoizedNC
struct MemoizedNC{N<:Neighborhood} <: NeighborhoodCache
    neighborhoods_F::Vector{N}
    neighborhoods_B::Vector{N}
end
function MemoizedNC(::Type{T}, n::Int=0, issymmetric::Bool=false) where {T}
    if issymmetric
        neighborhoods = Vector{Vector{T}}(undef, n)
        MemoizedNC(neighborhoods, neighborhoods)
    else
        MemoizedNC(Vector{Vector{T}}(undef, n), Vector{Vector{T}}(undef, n))
    end
end
function MemoizedNC(nodes::SampleSet, bvp::SteeringBVP; include_controls::BoolVal=Val(false))
    MemoizedNC(neighbor_info_type(nodes, bvp, include_controls=include_controls), length(nodes), issymmetric(bvp))
end
includes_controls( ::Type{MemoizedNC{N}}) where {N} = includes_controls(N)
neighbor_info_type(::Type{MemoizedNC{N}}) where {N} = neighbor_info_type(N)
function neighbors(cache::MemoizedNC{N}, nnds, nodes, bvp, v, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false)) where {N}
    neighborhoods = dir === Val(:F) ? cache.neighborhoods_F : cache.neighborhoods_B
    i = linear_index(v)
    if i <= 0 || !isassigned(neighborhoods, i)
        neighborhood = neighbors(nnds, nodes, bvp, v, r, dir=dir, include_controls=include_controls)
        i > 0 && (neighborhoods[i] = collect_if_not_Array(neighborhood))
        neighborhood
    else
        neighborhoods[i]
    end
end

## OfflineNC
struct OfflineNC{PN<:PackedNeighborhoods} <: NeighborhoodCache
    neighborhoods_F::PN
    neighborhoods_B::PN
end
function OfflineNC(cache::MemoizedNC)
    if cache.neighborhoods_F === cache.neighborhoods_B
        packed_neighborhoods = PackedNeighborhoods(cache.neighborhoods_F)
        OfflineNC(packed_neighborhoods, packed_neighborhoods)
    else
        OfflineNC(PackedNeighborhoods(cache.neighborhoods_F), PackedNeighborhoods(cache.neighborhoods_B))
    end
end
includes_controls( ::Type{OfflineNC{PN}}) where {PN} = includes_controls(PN)
neighbor_info_type(::Type{OfflineNC{PN}}) where {PN} = neighbor_info_type(PN)
function neighbors(cache::OfflineNC{PN}, nnds, nodes, bvp, v, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false)) where {PN}
    dir === Val(:F) ? cache.neighborhoods_F[linear_index(v)] : cache.neighborhoods_B[linear_index(v)]
end

# Near Neighbor Data Structures
abstract type NearNeighborDataStructure end
abstract type MetricNNDS <: NearNeighborDataStructure end    # `using NearestNeighbors` or `using FLANN`
## NullNNDS (Streaming)
struct NullNNDS <: NearNeighborDataStructure end
function neighbors(nnds::NullNNDS, nodes, bvp, v, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false))
    neighbors(nodes, bvp, v, r, dir=dir, include_controls=include_controls)
end

## MemoizedNNDS
struct MemoizedNNDS{S,D,U,F} <: NearNeighborDataStructure
    bvp_results::Dict{Tuple{S,S},NamedTuple{(:cost,:controls),Tuple{D,U}}}
    memoized_bvp::F
end
function MemoizedNNDS(nodes::ExplicitSampleSet, bvp::SteeringBVP; include_controls::BoolVal=Val(false))
    S = typeof(nodes.init)
    D = typeof(bvp(nodes.init, nodes.init).cost)
    U = include_controls === Val(true) ? typeof(bvp(nodes.init, nodes.init).controls) : Nothing
    bvp_results = Dict{Tuple{S,S},NamedTuple{(:cost,:controls),Tuple{D,U}}}()
    memoized_bvp = (x0, xf, cost_bound) -> begin    # TODO: check if closure performance is disastrous
        get!(bvp_results, (x0, xf)) do              #       also, make subtype of SteeringBVP?
            keep_controls_if(bvp(x0, xf, cost_bound), include_controls)
        end
    end
    MemoizedNNDS(bvp_results, memoized_bvp)
end
includes_controls(::MemoizedNNDS{S,D,U}) where {S,D,U} = U !== Nothing
function neighbors(nnds::MemoizedNNDS, nodes, bvp, v, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false))
    neighbors(nodes, nnds.memoized_bvp, v, r, dir=dir, include_controls=include_controls)
end

# General Methods
function neighbors(nodes::ExplicitSampleSet, bvp, v::Int, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false))
    neighbors(nodes, bvp, nodes[v], r, dir=dir, include_controls=include_controls, omit=isequal(v))
end
function neighbors(nodes::ExplicitSampleSet, bvp, x::State, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false), omit::Function=always_false)
    gen = dir === Val(:F) ? ((index=i, bvp(x, nodes[i], r)...) for i in eachindex(nodes) if !omit(i)) :
                            ((index=i, bvp(nodes[i], x, r)...) for i in eachindex(nodes) if !omit(i))
    (keep_controls_if(n, include_controls) for n in gen if n.cost <= r)
end

function neighbors(nodes::TiledSampleSet, bvp, v::TiledIndex, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false))
    if v.off_lattice <= 0
        neighbors(nodes, bvp, nodes[v], r, dir=dir, include_controls=include_controls)
    else
        if dir === Val(:F)
            x0_xf = i -> (nodes[zero_lattice_component(v)], nodes[shift_lattice_component(i, -v.on_lattice)])
        else
            x0_xf = i -> (nodes[zero_lattice_component(i)], nodes[shift_lattice_component(v, -i.on_lattice)])
        end
        neighbors(nodes, bvp, nodes[v], r, dir=dir, include_controls=include_controls, omit=isequal(v), x0_xf=x0_xf)
    end
end
function neighbors(nodes::TiledSampleSet, bvp, x::State, r;
                   dir::F_or_B=Val(:F), include_controls::BoolVal=Val(false), omit::Function=always_false,
                   x0_xf=(dir === Val(:F) ? i -> (x, nodes[i]) : i -> (nodes[i], x)))
    neighborhood = neighbor_info_type(nodes, bvp, include_controls=include_controls)[]
    close_idx = findclose(nodes, bvp, x, r, dir=dir, omit=omit, x0_xf=x0_xf)
    queue = [close_idx.on_lattice]
    visited = Set(queue)
    while !isempty(queue)
        on_lattice_idx = popfirst!(queue)
        found_a_neighbor = false
        for i in eachindex(nodes.Voff)
            idx = TiledIndex(on_lattice_idx, i)
            n = (index=idx, bvp(x0_xf(idx)..., r)...)
            if n.cost <= r
                found_a_neighbor = true
                !omit(idx) && push!(neighborhood, keep_controls_if(n, include_controls))
            end
        end
        found_a_neighbor && for oli in lattice_neighbors(on_lattice_idx)
            if !(oli in visited)
                push!(queue, oli)
                push!(visited, oli)
            end
        end
    end
    neighborhood
end
function findclose(nodes::TiledSampleSet, bvp, x::State, r; dir::F_or_B=Val(:F), omit::Function=always_false,
                   x0_xf=(dir === Val(:F) ? i -> (x, nodes[i]) : i -> (nodes[i], x)), cutoff::Int=30)
    start_lattice_idx = round.(Int, SVector{size(nodes.A,2)}(Matrix(nodes.A)\(x - nodes.b)))    # TODO: \ for non-square A
    sentinel = fill(typemin(Int), typeof(start_lattice_idx))
    queue = [start_lattice_idx, sentinel]
    visited = Set(queue)
    ct = 0
    while !isempty(queue)
        on_lattice_idx = popfirst!(queue)
        if on_lattice_idx == sentinel
            push!(queue, sentinel)
            ct += 1
            ct > cutoff && error("Either r is too small, or searching near the starting "*
                                 "state isn't working to find a close tiled sample!")
            continue
        end
        for i in eachindex(nodes.Voff)
            idx = TiledIndex(on_lattice_idx, i)
            omit(idx) && continue
            bvp(x0_xf(idx)..., r).cost <= r && return idx
        end
        for oli in lattice_neighbors(on_lattice_idx)
            if !(oli in visited)
                push!(queue, oli)
                push!(visited, oli)
            end
        end
    end
end

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
import Statistics
using Test, Random
using NiLang, NiLang.AD
using Distributions

@testset "statistics" begin
    x = randn(100)
    y = randn(100)
    @test i_mean_sum(0.0, 0.0, x)[1] ≈ Statistics.mean(x)
    info = VarianceInfo(Float64)
    @test almost_same(i_var_mean_sum(info, copy(x))[1], VarianceInfo(Statistics.var(x), Statistics.var(x)*99, Statistics.mean(x), sum(x)))
    @test almost_same((~i_var_mean_sum)(i_var_mean_sum(info, copy(x))...), (info, x))
    @test almost_same(i_cor_cov(0.0, 0.0, copy(x), copy(y)), (Statistics.cor(x,y), Statistics.cov(x,y), x, y))
    @test almost_same((~i_cor_cov)(i_cor_cov(0.0, 0.0, copy(x), copy(y))...), (0.0, 0.0, x, y))
end

@testset "normal log pdf" begin
    out = 0.0
    x = 1.0
    μ = 0.3
    σ = 1.5
    l1 = i_normal_logpdf(out, x, μ, σ)
    distri = Normal(μ, σ)
    l2 = logpdf(distri, x)
    @test l1[1] ≈ l2
    @test check_inv(i_normal_logpdf, (out, x, μ, σ))
end


text:
@testset "PaddedAxis" begin
    x = [:a, :b, :c, :d, :e]
    @test AxisArray(x, circular_pad(first_pad=3, last_pad=3)) == [:c, :d, :e, :a, :b, :c, :d, :e, :a, :b, :c]
    @test AxisArray(x, replicate_pad(first_pad=3, last_pad=3)) == [:a, :a, :a, :a, :b, :c, :d, :e, :e, :e, :e]
    @test AxisArray(x, symmetric_pad(first_pad=3, last_pad=3)) == [:d, :c, :b, :a, :b, :c, :d, :e, :d, :c, :b]
    @test AxisArray(x, reflect_pad(first_pad=3, last_pad=3)) == [:c, :b, :a, :a, :b, :c, :d, :e, :e, :d, :c]
    @test AxisArray(1:2, one_pad(sym_pad=2)) == [1, 1, 1, 2, 1, 1]
    @test AxisArray(1:2, zero_pad(sym_pad=2)) == [0, 0, 1, 2, 0, 0]

    x = reshape(1:6, 3, 2)
    ax = @inferred(AxisArray(x, replicate_pad(sym_pad=2), replicate_pad(sym_pad=2)))
    @test ax == @inferred(replicate_pad(x; sym_pad=2))
    @test eltype(x) <: eltype(ax)
    @test axes(ax) == (-1:5, -1:4)
    cax = collect(ax)
    @test cax == ax
    @test axes(cax, 1) isa OffsetAxis
    @test @inferred(IndexStyle(ax)) isa IndexCartesian

    axis = zero_pad(SimpleAxis(ArrayInterface.StaticInt(2)); sym_pad=ArrayInterface.StaticInt(2))
    @test ArrayInterface.known_first(axis) === -1
    @test ArrayInterface.known_last(axis) === 4
    @test ArrayInterface.known_length(axis) === 6

    axis = zero_pad(1:2, sym_pad=2)
    @test ArrayInterface.known_first(axis) === nothing
    @test ArrayInterface.known_last(axis) === nothing
    @test ArrayInterface.known_length(axis) === nothing

    axis = eachindex(AxisArray(3:4, one_pad(sym_pad=2)))
    x = map(identity, axis)
    @test x == axis
    @test axes(x, 1) == axis
end


text:
"""
    clean([distclean=false])

Clean up build/doc/testing artifacts. Restore to clean checkout state
(distclean)
"""
function clean(; distclean=false, _exit=true)

    _glob(folder, ending) =
        [name for name in readdir(folder; join=true) if (name |> endswith(ending))]
    _glob_star(folder; except=[]) = [
        joinpath(folder, name) for
        name in readdir(folder) if !(name |> startswith(".") || name ∈ except)
    ]
    _exists(name) = isfile(name) || isdir(name)
    _push!(lst, name) = _exists(name) && push!(lst, name)

    ROOT = dirname(@__DIR__)

    ###########################################################################
    CLEAN = String[]
    for folder in ["", "src", "test"]
        append!(CLEAN, _glob(joinpath(ROOT, folder), ".cov"))
    end
    append!(
        CLEAN,
        _glob_star(joinpath(ROOT, "docs", "src", "examples"), except=["index.md"])
    )
    _push!(CLEAN, joinpath(ROOT, "coverage"))
    _push!(CLEAN, joinpath(ROOT, "docs", "build"))
    _push!(CLEAN, joinpath(ROOT, "lcov.info"))
    ###########################################################################

    ###########################################################################
    DISTCLEAN = String[]
    for folder in ["", "docs", "test"]
        _push!(DISTCLEAN, joinpath(joinpath(ROOT, folder), "Manifest.toml"))
    end
    _push!(DISTCLEAN, joinpath(ROOT, "docs", "Project.toml"))
    _push!(DISTCLEAN, joinpath(ROOT, "docs", "src", "examples", ".ipynb_checkpoints"))
    _push!(DISTCLEAN, joinpath(ROOT, ".JuliaFormatter.toml"))
    ###########################################################################

    for name in CLEAN
        @info "rm $name"
        rm(name, force=true, recursive=true)
    end
    if distclean
        for name in DISTCLEAN
            @info "rm $name"
            rm(name, force=true, recursive=true)
        end
        if _exit
            @info "Exiting"
            exit(0)
        end
    end

end

distclean() = clean(distclean=true)


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using LinearAlgebra, Random, Statistics
using StatsBase, Distributions, ArraysOfArrays

@testset "reshaped_dist" begin
    scshape = ScalarShape{Real}()

    a1shape = ArrayShape{Real}(3)
    a2shape = ArrayShape{Real}(2, 3)

    ntshape = NamedTupleShape(
        a = ArrayShape{Real}(2, 3),
        b = ScalarShape{Real}(),
        c = ConstValueShape(4.2),
        x = ConstValueShape([11 21; 12 22]),
        y = ArrayShape{Real}(4)
    )

    names = (:a, :b, :c, :x, :y)

    _mvndist(n::Integer) = MvNormal(fill(2.0, n), Diagonal(fill(3.0, n)))
    _mvnormalrd(shape::AbstractValueShape) = ReshapedDist(_mvndist(totalndof(shape)), shape)
    _mvnormalrd2(shape::AbstractValueShape) = shape(_mvndist(totalndof(shape)))

    @testset "ctors" begin
        @test @inferred(_mvnormalrd(scshape)) isa Distribution{Univariate, Continuous}
        @test @inferred(_mvnormalrd(a1shape)) isa Distribution{Multivariate, Continuous}
        @test @inferred(_mvnormalrd(a2shape)) isa Distribution{Matrixvariate, Continuous}
        @test @inferred(_mvnormalrd(ntshape)) isa Distribution{ValueShapes.NamedTupleVariate{names}, Continuous}

        @test_throws ArgumentError ReshapedDist(MvNormal(Diagonal(fill(3.0, 2))), a2shape)

        @test @inferred(_mvnormalrd2(scshape)) isa ReshapedDist{Univariate, Continuous}
        @test @inferred(_mvnormalrd2(a1shape)) isa MvNormal
        @test @inferred(_mvnormalrd2(a2shape)) isa MatrixReshaped
        @test @inferred(_mvnormalrd2(ntshape)) isa ReshapedDist{ValueShapes.NamedTupleVariate{names}, Continuous}

        @inferred(varshape(_mvnormalrd2(ntshape))) == ntshape
        @inferred(unshaped(_mvnormalrd2(ntshape))) == MvNormal(fill(2.0, totalndof(ntshape)), Diagonal(fill(3.0, totalndof(ntshape))))
    end

    @testset "rand" begin
        @test @inferred(rand(_mvnormalrd(scshape))) isa Real
        @test @inferred(rand(_mvnormalrd(scshape), ())) isa AbstractArray{<:Real,0}
        @test @inferred(rand(_mvnormalrd(scshape), 7)) isa AbstractArray{<:Real,1}
        @test @inferred(rand(_mvnormalrd(scshape), (7,))) isa AbstractArray{<:Real,1}
        @test size(rand(_mvnormalrd(scshape), 7)) == (7,)

        let d = _mvndist(totalndof(a1shape))
            @test @inferred(a1shape(d)) === d
        end
        @test @inferred(rand(_mvnormalrd(a1shape))) isa AbstractVector{<:Real}
        @test @inferred(rand(_mvnormalrd(a1shape), ())) isa AbstractArray{<:AbstractVector{<:Real},0}
        @test @inferred(rand(_mvnormalrd(a1shape), 7)) isa AbstractArray{<:Real,2}
        @test @inferred(rand(_mvnormalrd(a1shape), (7,))) isa AbstractArray{<:AbstractVector{<:Real},1}
        @test size(rand(_mvnormalrd(a1shape), 7)) == (3, 7)
        @test size(rand(_mvnormalrd(a1shape), (7,))) == (7,)

        let d = _mvndist(totalndof(a2shape))
            @test @inferred(a2shape(d)) isa MatrixReshaped
            @test unshaped(a2shape(d)) === d
        end
        @test @inferred(rand(_mvnormalrd(a2shape))) isa AbstractMatrix{<:Real}
        @test @inferred(rand(_mvnormalrd(a2shape), 7)) isa AbstractArray{<:AbstractMatrix{<:Real},1}
        @test size(rand(_mvnormalrd(a2shape), 7)) == (7,)

        @test @inferred(rand(_mvnormalrd(ntshape))) isa NamedTuple{names}
        @test @inferred(rand(_mvnormalrd(ntshape), ())) isa ShapedAsNT{<:NamedTuple{names}}
        @test @inferred(rand(_mvnormalrd(ntshape), 7)) isa ShapedAsNTArray{<:NamedTuple{names},1}
        @test size(rand(_mvnormalrd(ntshape), 7)) == (7,)
        let X = rand(_mvnormalrd(ntshape), 7)
            @test @inferred(rand!(_mvnormalrd(ntshape), view(X, 1))) == view(X, 1)
            @test @inferred(rand!(_mvnormalrd(ntshape), X)) === X
        end
    end

    @testset "stats functions" begin
        @test @inferred(mean(_mvnormalrd(scshape))) ≈ 2
        @test @inferred(mode(_mvnormalrd(scshape))) ≈ 2
        @test @inferred(var(_mvnormalrd(scshape))) ≈ 3
        let rd = _mvnormalrd(scshape), ux = rand(unshaped(rd)), vs = varshape(rd)
            @test @inferred(pdf(rd, vs(ux)[])) == pdf(unshaped(rd), ux)
            @test @inferred(pdf(rd, vs(ux))) == pdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux)[])) == logpdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux))) == logpdf(unshaped(rd), ux)
            @test @inferred(insupport(rd, vs(ux)[])) == true
        end

        @test @inferred(mean(_mvnormalrd(a1shape))) ≈ fill(2, 3)
        @test @inferred(mode(_mvnormalrd(a1shape))) ≈ fill(2, 3)
        @test @inferred(var(_mvnormalrd(a1shape))) ≈ fill(3, 3)
        @test @inferred(cov(_mvnormalrd(a1shape))) ≈ Diagonal(fill(3.0, 3))
        let rd = _mvnormalrd(a1shape), ux = rand(unshaped(rd)), vs = varshape(rd)
            @test @inferred(pdf(rd, vs(ux))) == pdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux))) == logpdf(unshaped(rd), ux)
            @test @inferred(insupport(rd, vs(ux))) == true
        end

        @test @inferred(mean(_mvnormalrd(a2shape))) ≈ fill(2, 2, 3)
        @test @inferred(mode(_mvnormalrd(a2shape))) ≈ fill(2, 2, 3)
        @test @inferred(var(_mvnormalrd(a2shape))) ≈ fill(3, 2, 3)
        let rd = _mvnormalrd(a2shape), ux = rand(unshaped(rd)), vs = varshape(rd)
            @test @inferred(pdf(rd, vs(ux))) == pdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux))) == logpdf(unshaped(rd), ux)
            @test @inferred(insupport(rd, vs(ux))) == true
        end

        @test @inferred(mean(_mvnormalrd(ntshape))) == (a = [2.0 2.0 2.0; 2.0 2.0 2.0], b = 2.0, c = 4.2, x = [11 21; 12 22], y = [2.0, 2.0, 2.0, 2.0])
        @test @inferred(mode(_mvnormalrd(ntshape))) == (a = [2.0 2.0 2.0; 2.0 2.0 2.0], b = 2.0, c = 4.2, x = [11 21; 12 22], y = [2.0, 2.0, 2.0, 2.0])
        @test @inferred(var(_mvnormalrd(ntshape))) == (a = [3.0 3.0 3.0; 3.0 3.0 3.0], b = 3.0, c = 0.0, x = [0 0; 0 0], y = [3.0, 3.0, 3.0, 3.0])
        let rd = _mvnormalrd(ntshape), ux = rand(unshaped(rd)), vs = varshape(rd)
            @test @inferred(pdf(rd, vs(ux)[])) == pdf(unshaped(rd), ux)
            @test @inferred(pdf(rd, vs(ux))) == pdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux)[])) == logpdf(unshaped(rd), ux)
            @test @inferred(logpdf(rd, vs(ux))) == logpdf(unshaped(rd), ux)
            @test @inferred(insupport(rd, vs(ux)[])) == true
            @test @inferred(insupport(rd, vs(ux))) == true
        end

        let rd = ReshapedDist(_mvndist(5), ArrayShape{Real}(5)), X = rand(rd, 10), Xn = nestedview(X)
            @test @inferred(Distributions._pdf(rd, Xn[1])) == pdf(rd, Xn[1])
            @test @inferred(pdf(rd, X)) == pdf.(Ref(rd), Xn)

            @test @inferred(Distributions._logpdf(rd, Xn[1])) == logpdf(rd, Xn[1])
            @test @inferred(logpdf(rd, X)) == logpdf.(Ref(rd), Xn)
        end
    end
end


text:
#------------------------------------------------------------------------------
# hydraulic julia script
# based on mechanical example @ Ramos Jr. & Vieira [ICECE2000]
#
# ricardoaf@lccv.ufal.br (version 21/01/2021)
#------------------------------------------------------------------------------
include("DiscreteModel.jl")

function Hydraulic()

	node, elem, α = init(
		connectivity=[1 2; 2 4; 2 3; 2 5; 3 6; 6 5; 4 6; 6 7; 4 5; 3 5],
		stiffness=[3, 5, 3, 2, 2, 2, 4, 2, 3, 3],
		displacement=[10, 0, 0, 0, 0, 0, 0],
		force=[0, 0, 0, 0, 0, 0, -8],
		fixnode=[1],
		type="hydraulic")

	node, elem = run!(node, elem, α)

	println("     Pressures: ", node.u)
	println("         Flows: ", node.f)
	println("Internal flows: ", elem.f)
end

max_stars_repo_name: ricardoaf/discretemodel

max_issues_repo_name: ricardoaf/discretemodel

max_forks_repo_name: ricardoaf/discretemodel


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

# This file contains the implementation of the Edge DataFrame.

################################################# IMPORT/EXPORT ############################################################

export listeprops, haseprop, geteprop, seteprop!

################################################# MUTATION #################################################################

""" Reorder the edge dataframe to match the order of edges in the index table """
function reorder!(x::AbstractDataFrame, indxs::Vector{Int})
   for eprop in names(x)
      x[eprop] = x[eprop][indxs]
   end
end

""" Add a new row to the edge dataframe """
function addedge!(x::AbstractDataFrame)
   if !isempty(x)
      push!(x, @data fill(NA, ncol(x)))
   end
end

""" Remove a row from the edge dataframe """
function rmedge!(x::AbstractDataFrame, erow::Int)
   if !isempty(x)
      deleterows!(x, erow)
   end
end

""" Remove a row from the edge dataframe """
function rmedge!(x::AbstractDataFrame, erows::AbstractVector{Int})
   erows = unique(sort(erows))
   for prop in names(x)
      deleteat!(x[prop], erows)
   end
end

################################################# EPROPS ACCESSORS ##########################################################

""" List the column headers of the edge dataframe """
listeprops(g::Graph) = names(edata(g))

""" Check if the edge dataframe has the input edge property """
haseprop(g::Graph, eprop::Symbol) = haskey(edata(g), eprop)

################################################# GETEPROP ##################################################################

"""
Retrieve edge properties.

geteprop(g::Graph, e::EdgeID, eprop::Symbol) -> Fetch the value of a property for edge e
"""
geteprop(g::Graph, e::EdgeID, eprop::Symbol) = edata(g)[eprop][indxs(g)[e]]

""" geteprop(g::Graph, es::EdgeList, eprop::Symbol) -> Fetch the value of a property for edge e in es """
geteprop(g::Graph, es::EdgeList, eprop::Symbol) = edata(g)[eprop][indxs(g)[es]]

""" geteprop(g::Graph, ::Colon, eprop::Symbol) -> Fetch the value of a property for all edges """
geteprop(g::Graph, ::Colon, eprop::Symbol) = copy(edata(g)[eprop])

################################################# SETEPROP ##################################################################

"""
Set edge properties.

seteprop!(g::Graph, e::EdgeID, val, eprop::Symbol) -> Set a property for an edge e
"""
function seteprop!(g::Graph, e::EdgeID, val, eprop::Symbol)
   if haseprop(g, eprop)
      edata(g)[indxs(g)[e], eprop] = val
   else
      error("Property $eprop doesn't exist. Please create it on all edges first.")
   end
end

""" seteprop!(g::Graph, es::EdgeList, val(s), eprop::Symbol) -> Set a property for e in es """
function seteprop!(g::Graph, es::EdgeList, val, eprop::Symbol)
   if haseprop(g, eprop)
      edata(g)[eprop][indxs(g)[es]] = val
   else
      error("Property $eprop doesn't exist. Please create it on all edges first.")
   end
end

""" seteprop!(g::Graph, ::Colon, val(s), eprop::Symbol) """
function seteprop!(g::Graph, ::Colon, vals::AbstractVector, eprop::Symbol)
   if length(vals) == ne(g)
      edata(g)[eprop] = vals
   else
      error("Trying to set $(length(vals)) values to $(ne(g)) properties")
   end
end

seteprop!(g::Graph, ::Colon, val, eprop::Symbol) = seteprop!(g, :, fill(val, ne(g)), eprop)

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

""" Retrieve a subset of the edge dataframe for the input edges """
subgraph(x::AbstractDataFrame, es::EdgeList) = edata(g)[indxs(g)[es],:]

"""
Retrieve a subset of the edge dataframe for the input edges, containing a
subset of the edge properties
"""
subgraph(x::AbstractDataFrame, es::EdgeList, eprops::Vector{Symbol}) = edata(g)[indxs(g)[es], eprops]