text:
"""
GuillotineModels.PPG2KP.Heuristic implements the heuristic needed for faithful
reimplementation of the Priced PP-G2KP from \\`F. Furini, E. Malaguti, and D.
Thomopulos, \\`\\`Modeling Two-Dimensional Guillotine Cutting Problems via Integer
Programming,'' INFORMS Journal on Computing, vol. 28, no. 4, pp. 736–751, Oct.
2016, doi: 10.1287/ijoc.2016.0710.' The heuristic itself, however, is older
and described in: \\`M. Dolatabadi, A. Lodi, and M. Monaci, \\`\\`Exact algorithms
for the two-dimensional guillotine knapsack,'' Computers \\& Operations
Research, vol. 39, no. 1, pp. 48–53, Jan. 2012, doi: 10.1016/j.cor.2010.12.018.'

The methods exported by this module are not of general interest, you need to
either want: (1) a fast but shelf-restricted heuristic for the guillotine 2D
knapsack problem; (2) to check if they were implemented correctly;
(3) to use them to make your own implementation of Priced PP-G2KP.
"""
module Heuristic

using Random#: shuffle!
import ...CutPattern
import ...TIMER # Global module timer for use with TimerOutputs.@timeit.
import TimerOutputs.@timeit
using ...Utilities: expand

export first_fit_decr, iterated_greedy
export Levels, fast_first_fit_decr, fast_iterated_greedy

"""
    promising_kfirst(::Type{D}, p, a, bkv, A) :: D

!!! **Internal use.**

Zero if there is no `k` for which `sum(p[1:k]) > bkv && sum(a[1:k]) <= A`;
otherwise, return the smallest `k` value satisfying these conditions.

In other words, without changing the order of the given sequence, return size
of the prefix of smallest size (a piece subset) that has better profit than the
best known value and that it is not immediately obvious that it cannot fit the
given area (does not break the naive area bound). If such prefix does not
exist, then just return zero.
"""
function promising_kfirst(
	::Type{D},
	p :: AbstractVector{P},
	a :: AbstractVector{P},
	bkv :: P,
	A :: P
) :: D where {D, P}
	@assert length(p) == length(a)
	sump, suma = zero(P), zero(P)
	for k in eachindex(p)
		suma += a[k]
		suma > A && return zero(D)
		sump += p[k]
		sump > bkv && return convert(D, k)
	end
	# the line below is only executed if sum(p) <= bkv
	return zero(D)
end

"""
    shelves2cutpattern(shelves :: [[D]], l, w, L, W) :: CutPattern{D,S}

Transform the third return of `iterated_greedy` or `first_fit_decr` into a
CutPattern object.
"""
function shelves2cutpattern(
	shelves :: Vector{Vector{D}}, l :: Vector{S}, w :: Vector{S}, L :: S, W :: S
) :: CutPattern{D, S} where {D, S}
	outer_cp = Vector{CutPattern{D, S}}()
	for l_strip in shelves
		inner_cp = Vector{CutPattern{D, S}}()
		for pii in l_strip
			push!(inner_cp, CutPattern(l[pii], w[pii], pii))
		end
		push!(outer_cp, CutPattern(L, w[first(l_strip)], false, inner_cp))
	end
	return CutPattern(L, W, true, outer_cp)
end

"""
    iterated_greedy(d, p, l, w, L, W, rng, max_iter_since_last_improv = 10^6)

The iterated greedy procedure is a heuristic that creates a random piece subset
that may improve the best known value (BKV) and call the first fit
with-decreasing heuristic over it (if this inner heuristic is able to position
all pieces then the BKV is improved).

The first fit heuristic is deterministic and gives the same solution for a
piece subset (the ordering is irrelevant), so the iterated procedure is adding
diversity to it. The first fit heuristic is also shelf-restricted, so the
solution is always shelf-restricted and may be impossible to reach the
optimality (in the case no optimal solution is shelf-restricted).

The procedure is repeated until `max_iter_since_last_improv` calls to
the first fit width-decreasing were made without improving the BKV.

# Arguments

1. `d::AbstractVector{D}`: The pieces demand.
2. `p::AbstractVector{P}`: The pieces profit.
3. `l::AbstractVector{S}`: The pieces length.
4. `w::AbstractVector{S}`: The pieces width.
5. `L::S`: The original plate length.
6. `W::S`: The original plate width.
7. `rng::AbstractRNG`: The random number generator.
8. `max_iter_since_last_improv::I`: The parameter for the stopping criteria,
  after `max_iter_since_last_improv` iterations since the last change in the
  best known value the algorithm ends.

# Returns

The same of the [`first_fit_decr`](@ref) method (the iterated greedy just
returns the best solution found by it).
"""
@timeit TIMER function iterated_greedy(
	d :: AbstractVector{D},
	p :: AbstractVector{P},
	l :: AbstractVector{S},
	w :: AbstractVector{S},
	L :: S,
	W :: S,
	rng :: AbstractRNG,
	max_iter_since_last_improv :: I = 1000000
) :: Tuple{P, Vector{D}, Vector{Vector{D}}} where {D, S, P, I}
	n = length(d) # number of distinct piece types
	A = convert(P, L) * convert(P, W)
	# assert explanation: all vectors should be the same length.
	@assert [n] == unique!(length.([p, l, w]))
	# loop control variables
	curr_iter = 0 # current iteration
	iter_last_improv = 0 # last iteration in which bkv changed
	# The triple from first_fit_decr (they will be overwritten unless n == 0)
	bkv = zero(P) # best known value
	sel = zeros(D, n) # like 'd' but with the pieces selected for the solution
	pat = Vector{Vector{D}}() # the cut pattern (is always a 2-stage cut)
	# Expand all vector by the demand. This corrects the algorithm for instances
	# in which not all piece types have just a single piece.
	pe, le, we = expand.((d, d, d), (p, l, w))
	# Create a first ordering by efficiency (i.e., cost-benefit, profit-to-area
	# ratio). This way, the same instance with the same RNG will have the same
	# results, independent from the initial ordering of the piece types.
	ie = expand(d, collect(1:n)) # permutation used to keep the original order
	ae = @. convert(P, le) * convert(P, we)
	ee = @. pe / convert(Float64, ae)
	oe = convert.(D, sortperm(ee; rev = true))

	while curr_iter - iter_last_improv <= max_iter_since_last_improv
		curr_iter += 1
		permute!.((pe, le, we, ae, ie), (oe, oe, oe, oe, oe))
		k = promising_kfirst(D, pe, ae, bkv, A)
		if !iszero(k)
			pe_,le_,we_,ae_,ie_ = pe[1:k],le[1:k],we[1:k],ae[1:k],ie[1:k]
			oe_ = sortperm(we_; rev = true)
			permute!.((pe_, le_, we_, ae_, ie_), (oe_, oe_, oe_, oe_, oe_))
			new_v, new_sel, new_pat = first_fit_decr(pe_, le_, we_, ie_, n, L, W)
			@assert sum(new_sel .* p) == new_v
			if new_v > bkv
				bkv, sel, pat = new_v, new_sel, new_pat
			end
		end
		shuffle!(rng, oe)
	end

	#=
	nz_i = collect(i for i in 1:n if !iszero(sel[i]))
	@show nz_i
	nz_sel = collect(sel[i] for i in 1:n if !iszero(sel[i]))
	@show nz_sel
	nz_l = collect(l[i] for i in 1:n if !iszero(sel[i]))
	@show nz_l
	nz_w = collect(w[i] for i in 1:n if !iszero(sel[i]))
	@show nz_w
	=#
	return bkv, sel, pat
end

# Pre-allocated struct for fast_first_fit_decr.
struct Levels{D, S}
	lw :: Vector{S} # width of the level
	ll :: Vector{S} # length of the level
	la :: Vector{D} # the amount of pieces in the level
	qt :: Base.RefValue{D} # quantity of levels (type `Ref` is abstract)
	function Levels{D, S}(max_length :: D) where {D, S}
		return new(
			Vector{S}(undef, max_length),
			Vector{S}(undef, max_length),
			Vector{D}(undef, max_length),
			Ref(zero(D))
		)
	end
end

Base.empty!(c :: Levels{D, S}) where {D, S} = c.qt[] = zero(D)

function _extract_shelves(
	pat :: AbstractArray{D, 2},
	lvls :: Levels{D, S}
) where {D, S}
	shelves = [Vector{D}(undef, lvls.la[i]) for i = 1:lvls.qt[]]
	for i = 1:lvls.qt[]
		shelve = shelves[i]
		for j = 1:lvls.la[i]
			shelve[j] = pat[j, i]
		end
	end
	return shelves
end

# This method changes the first `length(indexes)` positions of `buffer`
# to have the same value as `values[indexes]` would return. Then it returns
# the two first parameters for convenience (the idea here is swap their
# bindings and attribute `buffer` to the old name of `values` and `values` to
# the old name of `buffer`).
function _swap_permute!(buffer, values, indexes, k = length(indexes))
	for (new, old) in zip(1:k, indexes)
		buffer[new] = values[old]
	end
	return buffer, values
end
# The same as above for the inverse of the permutation.
function _swap_invpermute!(buffer, values, indexes, k = length(indexes))
	for (new, old) in zip(1:k, indexes)
		buffer[old] = values[new]
	end
	return buffer, values
end

# Shamelessly stolen from: https://github.com/JuliaLang/julia/blob/
# 81a044f988abd26509d9c1db53f6e9eee2f5da6f/base/sort.jl#L493-L508
# License is MIT: https://julialang.org/license
function _sort_range_rev_by_insertion!(
	v::AbstractVector, lo::Integer, hi::Integer
)
	@inbounds for i = lo+1:hi
		j = i
		x = v[i]
		while j > lo
			if x > v[j-1]
				v[j] = v[j-1]
				j -= 1
				continue
			end
			break
		end
		v[j] = x
	end
	return v
end

@timeit TIMER function fast_iterated_greedy(
	d :: AbstractVector{D},
	p :: AbstractVector{P},
	l :: AbstractVector{S},
	w :: AbstractVector{S},
	L :: S,
	W :: S,
	rng :: AbstractRNG,
	max_iter_since_last_improv :: I = 1000000
) :: Tuple{P, Vector{D}, Vector{Vector{D}}} where {D, S, P, I}
	n = convert(D, length(d)) # number of distinct piece types
	A = convert(P, L) * convert(P, W)
	# assert explanation: all vectors should be the same length.
	@assert [n] == unique!(length.([p, l, w]))
	# loop control variables
	curr_iter = 0 # current iteration
	iter_last_improv = 0 # last iteration in which bkv changed
	# Expand all vector by the demand. This corrects the algorithm for instances
	# in which not all piece types have just a single piece.
	pe = expand(d, p)
	le = expand(d, l)
	we = expand(d, w)
	m = length(pe) # the same as sum(d), the number of pieces (not piece types)
	# The values for first_fit_decr (they will be overwritten unless n == 0)
	bkv = zero(P) # best known value
	sel = zeros(D, n) # like 'd' but with the pieces selected for the solution
	pat = Array{D, 2}(undef, (m, m)) # the cut pattern (is always a 2-stage cut)
	lvls = Levels{D, S}(m)
	# Create a first ordering by efficiency (i.e., cost-benefit, profit-to-area
	# ratio). This way, the same instance with the same RNG will have the same
	# results, independent from the initial ordering of the piece types.
	ie = expand(d, collect(one(D):n)) # permutation used to keep the original order
	ae = @. convert(P, le) * convert(P, we)
	ee = @. pe / convert(Float64, ae)
	oe = convert.(D, sortperm(ee; rev = true))

	# Create unintialized copies of many structures above for using the
	# technique of overwrite-and-swap-bindings, which avoid repeated allocation,
	# or just use as a buffer.
	aux_sel, aux_pat = similar(sel), similar(pat)
	aux_lvls = Levels{D, S}(m)
	aux_pe, aux_le, aux_we, aux_ae, aux_ie = similar.((pe, le, we, ae, ie))
	prefix_oe = Vector{D}(undef, m)
	permsort_buffer = Vector{Tuple{S, D}}(undef, m)
	# I think I do not need UnsafeArrays. The reason for views are:
	# (1) fast_first_fit_decr, that is trivial to receive a `k` value and then
	# just work over that subset; (2) getting the sortperm of a prefix of one
	# of the arrays; (3) permuting the prefix of other arrays inplace considering
	# this prefix.
	while curr_iter - iter_last_improv <= max_iter_since_last_improv
		curr_iter += 1
		pe, aux_pe = _swap_permute!(aux_pe, pe, oe)
		le, aux_le = _swap_permute!(aux_le, le, oe)
		we, aux_we = _swap_permute!(aux_we, we, oe)
		ae, aux_ae = _swap_permute!(aux_ae, ae, oe)
		ie, aux_ie = _swap_permute!(aux_ie, ie, oe)
		k = promising_kfirst(D, pe, ae, bkv, A)
		if !iszero(k)
			# Create a sortperm of the first `k` elements of `we` with as
			# little allocation as possible.
			for i = one(D):k; permsort_buffer[i] = (we[i], i); end
			# If this call to a internal undocumented `Base.Sort.sort!` method ever
			# breaks, just comment it and uncomment the call to InsertionSort right
			# below, it can be a little slower for larger sets, but works fine
			# enough. The code of the method below was not copyed because it makes
			# use of many helper methods (as opposed to InsertionSortAlg).
			Base.Sort.sort!(
				permsort_buffer, 1, k, Base.Sort.QuickSort, Base.Sort.Reverse
			)
			#_sort_range_rev_by_insertion!(permsort_buffer, 1, k)
			for i = one(D):k; prefix_oe[i] = last(permsort_buffer[i]); end
			pe, aux_pe = _swap_permute!(aux_pe, pe, prefix_oe, k)
			le, aux_le = _swap_permute!(aux_le, le, prefix_oe, k)
			we, aux_we = _swap_permute!(aux_we, we, prefix_oe, k)
			ae, aux_ae = _swap_permute!(aux_ae, ae, prefix_oe, k)
			ie, aux_ie = _swap_permute!(aux_ie, ie, prefix_oe, k)
			new_value = fast_first_fit_decr!(
				aux_sel, aux_pat, aux_lvls, pe, le, we, ie, k, L, W
			)
			pe, aux_pe = _swap_invpermute!(aux_pe, pe, prefix_oe, k)
			le, aux_le = _swap_invpermute!(aux_le, le, prefix_oe, k)
			we, aux_we = _swap_invpermute!(aux_we, we, prefix_oe, k)
			ae, aux_ae = _swap_invpermute!(aux_ae, ae, prefix_oe, k)
			ie, aux_ie = _swap_invpermute!(aux_ie, ie, prefix_oe, k)
			#@assert sum(aux_sel .* p) == new_value
			if new_value > bkv
				bkv = new_value
				sel, aux_sel = aux_sel, sel
				pat, aux_pat = aux_pat, pat
				lvls, aux_lvls = aux_lvls, lvls
			end
		end
		shuffle!(rng, oe)
	end # while

	#=
	nz_i = collect(i for i in 1:n if !iszero(sel[i]))
	@show nz_i
	nz_sel = collect(sel[i] for i in 1:n if !iszero(sel[i]))
	@show nz_sel
	nz_l = collect(l[i] for i in 1:n if !iszero(sel[i]))
	@show nz_l
	nz_w = collect(w[i] for i in 1:n if !iszero(sel[i]))
	@show nz_w
	=#
	return bkv, sel, _extract_shelves(pat, lvls)
end

"""
    first_fit_decr(p, l, w, t, n, L, W) :: (P, [D], [[D]])

The first fit width-decreasing heuristic is a deterministic, fast, and
shelf-restricted heuristic for the knapsack 2D guillotine problem.

This method is mainly called by [`iterated_greedy`](@ref) and because
of this the arguments are adjusted so the method is aware it is receiving a
subset of the complete piece set.

# Arguments

1. `p::AbstractVector{P}`: The pieces profit.
2. `l::AbstractVector{S}`: The pieces length.
3. `w::AbstractVector{S}`: The pieces width.
4. `t::AbstractVector{D}`: The piecex type/index/identifier (index in the
   original piece set).
5. `n::D`: The number of pieces in set (necessary to build the second
   returned vector).
6. `L::S`: The original plate length.
7. `W::S`: The original plate width.

# Returns

1. The best know value, i.e., the value of the returned solution.
2. A vector `v` of the same size as the number of pieces; if `v[i] == x`
   then the piece `i` has `x` copies in the solution returned.
3. A solution with the best known value, as the heuristic only consider
   shelf-restricted configurations the solution may be represented as a
   vector of vectors. Each vector is a length strip (with the width of
   the first piece inside, and the summed length of all pieces inside).
   Inside each vector the pieces are ordered by non-increasing width,
   and the vectors/stripes are sorted by non-increasing width too (
   i.e., the width of their first element).
"""
function first_fit_decr(
	p :: AbstractVector{P},
	l :: AbstractVector{S},
	w :: AbstractVector{S},
	t :: AbstractVector{D},
	n :: D, # number of distinct piece types (n >= maximum(t))
	L :: S, # length of the plate
	W :: S # width of the plate
) :: Tuple{P, Vector{D}, Vector{Vector{D}}} where {D, S, P}
	# assert explanation: all vectors should be the same length.
	@assert isone(length(unique!(length.([p, l, w, t]))))
	# assert explanation: this is first-fit WIDTH-decreasing
	@assert issorted(w; rev = true)
	lw = Vector{S}() # levels width
	ll = Vector{S}() # levels length
	rw = W # remaining width

	bkv = zero(P) # best known value
	sel = zeros(D, n) # if piece i is cut then sel[t[i]] is incremented
	pat = Vector{Vector{D}}() # the pattern (always a 2-stage pattern)
	for i in 1:length(p) # for each piece
		assigned = false
		for j in 1:length(lw) # for each already open level
			# assert reasoning: the levels are created with the width of previous
			# pieces, and the pieces are ordered by decreasing width, therefore
			# no piece may have more width than an already open level.
			@assert lw[j] >= w[i]
			if ll[j] >= l[i] # checks if the level has enough length yet
				ll[j] -= l[i]
				bkv += p[i]
				sel[t[i]] += one(D)
				push!(pat[j], t[i])
				assigned = true
				break
			end
		end
		# If the piece does not fit any already open level...
		if !assigned
			# ... and there is no space to open a new level, then jump to the
			# next piece type.
			rw < w[i] && continue
			# ... but there is yet space to open a new level, create a new level.
			push!(lw, w[i])
			push!(ll, L - l[i])
			rw -= w[i]
			bkv += p[i]
			sel[t[i]] += one(D)
			push!(pat, [t[i]])
		end
	end

	return bkv, sel, pat
end

function fast_first_fit_decr!(
	sel :: AbstractArray{D},
	pat :: AbstractArray{D, 2},
	c :: Levels{D, S},
	p :: AbstractVector{P},
	l :: AbstractVector{S},
	w :: AbstractVector{S},
	t :: AbstractVector{D},
	k :: D, # consider only the first k positions in p, l, w, and, t
	L :: S, # length of the plate
	W :: S # width of the plate
) :: P where {D, S, P}
	rw = W # remaining width
	fill!(sel, zero(eltype(sel)))
	empty!(c)
	# The `pat` parameter is not reset because the code already
	# works on the assumption that it is a dirty buffer.
	bkv = zero(P) # best known value
	for i in one(D):k # for each piece selected
		assigned = false
		for j in one(D):c.qt[] # for each already open level
			# assert reasoning: the levels are created with the width of previous
			# pieces, and the pieces are ordered by decreasing width, therefore
			# no piece may have more width than an already open level.
			#@assert lw[j] >= w[i]
			if c.ll[j] >= l[i] # checks if the level has enough length yet
				c.ll[j] -= l[i]
				bkv += p[i]
				sel[t[i]] += one(D)
				c.la[j] += one(D) # increase the amount of pieces in that level
				pat[c.la[j], j] = t[i]
				assigned = true
				break
			end
		end
		# If the piece does not fit any already open level...
		if !assigned
			# ... and there is no space to open a new level, then jump to the
			# next piece type.
			rw < w[i] && continue
			# ... but there is yet space to open a new level, create a new level.
			c.qt[] += one(D)
			c.lw[c.qt[]] = w[i]
			c.ll[c.qt[]] = L - l[i]
			c.la[c.qt[]] = one(D)
			rw -= w[i]
			bkv += p[i]
			sel[t[i]] += one(D)
			pat[one(D), c.qt[]] = t[i]
		end
	end

	return bkv
end
end # module


text:
#EasyPlot attribute system
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#

#==Main types
===============================================================================#

"Object that supports attributes changes"
abstract type AbstractAttributeReceiver; end

struct AttributeChangeData
	base::Tuple
	named::Base.Iterators.Pairs
end
#_"Constructor", avoiding signature conflicts
_AttributeChangeData(args...; kwargs...) = AttributeChangeData(args, kwargs)
_AttributeChangeData() = _AttributeChangeData(tuple())

struct AttributeChangeSpec{ID}
	data::AttributeChangeData
end
AttributeChangeSpec(id::Symbol, data::AttributeChangeData) =
	AttributeChangeSpec{id}(data)
AttributeChangeSpec(id::Symbol, data) =
	AttributeChangeSpec{id}(_AttributeChangeData(data))

#==Validation
===============================================================================#

#==Helper functions
===============================================================================#

#==Main functions
===============================================================================#

function attributes(;kwargs...)
	#Each key identifies "type of attribute" - as in `AttributeChangeSpec{type}`
	#Each value resolves/gets wrapped to an `AttributeChangeData`
	result = AttributeChangeSpec[]
	for (k, v) in kwargs
		push!(result, AttributeChangeSpec(k, v))
	end
	return result
end

#==apply() interface (implements set interface)
===============================================================================#
#Calls user-defined set functions:
#   _apply(ar::T, ::DS{AT}, args...; kwargs...)
function _apply(ar::AbstractAttributeReceiver, args::Tuple, kwargs::Base.Iterators.Pairs)
	getDS(::AttributeChangeSpec{AT}) where AT = DS(AT)

	for arg in args
		if !isa(arg, Vector{AttributeChangeSpec})
			T = typeof(ar)
			msg = "set(::$T, ...): non-keyword arguments must be created from attributes():\n    $arg"
			throw(ArgumentError(msg))
		end
		for attr in arg
			_apply(ar, getDS(attr), attr.data.base...; attr.data.named...)
		end
	end
	alist = attributes(;kwargs...)
	for attr in alist
		_apply(ar, getDS(attr), attr.data.base...; attr.data.named...)
	end
	return ar
end

#==Register constructors with cons() interface
===============================================================================#
cons(::DS{:a}; kwargs...) = attributes(; kwargs...)
cons(::DS{:attr}; kwargs...) = attributes(; kwargs...)
cons(::DS{:attribute_list}; kwargs...) = attributes(; kwargs...)

#==set() interface
===============================================================================#
"""
    set(args...; kwargs...)

Construct `AttributeChangeData` object if first argument is not special.

# Example
```julia-repl
set(style=:dash, width=3, color=:red)

#Change specs (`AttributeChangeSpec[]`) are created with `attributes`-`set` combinations:
lineattributes = attributes(
    line = set(style=:dash, width=3, color=:red)
)
```
"""
set(args...; kwargs...) = AttributeChangeData(args, kwargs)

"""
    set(ar::AbstractAttributeReceiver, [alist1], [alist2], ...; [akv_expr1], [akv_expr1], ...)

Set attributes to an `::AbstractAttributeReceiver` object using:
 1. `AttributeChangeSpec[]` created by `attributes()` function, and
 2. Attribute key-value expressions, **as used** by `attributes()` function.

Applicable types implementing this feature can be found using:
```julia-repl
julia> subtypes(EasyPlot.AbstractAttributeReceiver)
```

See also: [`attributes`](@ref)
"""
set(ar::AbstractAttributeReceiver, args...; kwargs...) = _apply(ar, args, kwargs)

#==Show functions
===============================================================================#
function Base.show(io::IO, d::AttributeChangeData)
	print(io, "{", d.base, ", ")

	#Show "named" attribute changes:
	rmg = length(d.named)
	for (k,v) in d.named
		print(io, k, "=")
		show(io, v)
		rmg -= 1
		if rmg!=0; print(io, ", "); end
	end

	print(io, "}")
end

function Base.show(io::IO, s::AttributeChangeSpec{T}) where T
	print(io, "Δattr{$T}: ", s.data)
end

#Last line

ext: jl


text:
module DataGeneratorTranslators

export regex_generator, type_generator, xsd_generator, bnf_generator

using ParserCombinator # BNF parsing
using LightXML # XSD parsing

include("parse.jl")
include("transform.jl")
include("build.jl")

# utility function to support direct inclusion of generator resulting from translation
function include_generator(genname::Symbol, translatefn::Function, translateargs...)
	genbuf = IOBuffer()
	translatefn(genbuf, genname, translateargs...)
	genstr = takebuf_string(genbuf)
	include_string(genstr)
	current_module().eval(genname)
end

include(joinpath("regex", "regex.jl"))
#include(joinpath("type", "type.jl"))
include(joinpath("bnf", "bnf.jl"))
include(joinpath("xsd", "xsd.jl"))

end


text:
# References:
# [1] 3-Axis Digital Accelerometer ADXL345 - D07925-0-5/09(0)

module PiADXL345

export adxl_open

using BBSPI
using PiGPIOMEM

"""
GPIOPin with inverted logic (Low = `true`, High = `false`).
"""
struct ActiveLowPin
    pin::GPIOPin
end
Base.setindex!(p::ActiveLowPin, v) = p.pin[] = iszero(v) ? 1 : 0
Base.getindex(p::ActiveLowPin) = iszero(p.pin[])

"""
Bit-delay for bit-bangged SPI.
"""
BBSPI.delay(s::BBSPI.SPISlave) = PiGPIOMEM.spin(50)

struct ADXL345{T} <: AbstractChannel{Tuple{Float64,Float64,Float64}}
    spi::T
    buf::Vector{UInt8}
    ADXL345(spi::T) where T = new{T}(spi, zeros(UInt8, 10))
end

"""
    adxl_open(;cs=chip select output pin,
              sdo=master input pin,
              sda=master output pin,
              scl=clock output pin)::ADXL345

Open ADXL345 conencted to GPIO pins `cs`, `sdo`, `sda` and `scl`.

Methods must be defined for `setindex!(::PinType, v)` and `getindex(::PinType)`.
See `help?> SPISlave`.
"""
function adxl_open(;cs=nothing, scl=nothing, sda=nothing, sdo=nothing)

    adxl_open(BBSPI.SPISlave(cs=ActiveLowPin(cs),
                             clk=ActiveLowPin(scl),
                             mosi=sda,
                             miso=sdo))
end

function adxl_open(spi)

    spi.chip_select[] = 0
    spi.clock[] = 0

    adxl = ADXL345(spi)

    @assert adxl_is_connected(adxl)

    adxl_enable(adxl)
    adxl_bw_rate_12hz5(adxl)

    adxl
end

adxl_is_connected(adxl) = adxl_read(adxl, 0) == 0xE5 # DEVID [1, p15]

adxl_enable(adxl) = adxl_write(adxl, 0x2D, 0b00001000) # POWER_CTL [1, p16]

adxl_offset(adxl) = [adxl_read(adxl, 0x1E), # OFSX [1, p15]
                     adxl_read(adxl, 0x1F), # OFSY [1, p15]
                     adxl_read(adxl, 0x20)] # OFSZ [1, p15]

function adxl_set_offset(adxl, x, y, z)
    adxl_write(adxl, 0x1E, unsigned(Int8(x))) # OFSX [1, p15]
    adxl_write(adxl, 0x1F, unsigned(Int8(y))) # OFSY [1, p15]
    adxl_write(adxl, 0x20, unsigned(Int8(z))) # OFSZ [1, p15]
end

adxl_bw_rate_12hz5(adxl) = adxl_write(adxl, 0x2C, 0b00000111) # BW_RATE [1, p16]
adxl_bw_rate_25hz(adxl)  = adxl_write(adxl, 0x2C, 0b00001000) # BW_RATE [1, p16]
adxl_bw_rate_50hz(adxl)  = adxl_write(adxl, 0x2C, 0b00001001) # BW_RATE [1, p16]
adxl_bw_rate_100hz(adxl) = adxl_write(adxl, 0x2C, 0b00001010) # BW_RATE [1, p16]

"""
    adxl_read(adxl, address, n=1)

Read `n` bytes from `address`.
"""
function adxl_read(adxl, address, n)
    n += 1
    cout = UInt8[0b11000000 | address] # SPI header [1, p9 Figure 5]
    BBSPI.transfer(adxl.spi, cout, view(adxl.buf, 1:n))
    return view(adxl.buf, 2:n)
end

adxl_read(adxl, address) = adxl_read(adxl, address, 1)[1]

"""
    adxl_write(adxl, address, v)

Write byte `v` to `address`.
"""
function adxl_write(adxl, address, v)
    cout = UInt8[address, v] # SPI header [1, p9 Figure 5]
    BBSPI.transfer(adxl.spi, cout)
    nothing
end

"""
    take!(::ADXL345)

Read [x,y,z] vector from ADXL345.
"""
function Base.take!(adxl::ADXL345)

    v = adxl_read(adxl, 0x32, 6) # DATAX0... Register [1, p18]
    x = signed(v[1] | UInt16(v[2]) << 8)/256 # 0.0039g/LSB [1, Table 1, p3]
    y = signed(v[3] | UInt16(v[4]) << 8)/256
    z = signed(v[5] | UInt16(v[6]) << 8)/256
    (x, y, z)
end

pitch_and_roll(v) = (pitch = atan(-v[1], hypot(v[2], v[3])) * 180 / π,
                     roll  = atan( v[2],             v[3])  * 180 / π)

function adxl_demo()
    adxl = adxl_open(cs=GPIOPin(26; output=true),
                     scl=GPIOPin(21; output=true),
                     sda=GPIOPin(20; output=true),
                     sdo=GPIOPin(16))

    while true
        t = time()
        v = take!(adxl)
        pr = pitch_and_roll(v)
        print(round(pr.pitch; digits=3), " ",
              round(pr.roll; digits=3), "      \r")
        sleep(max(0, 0.08 - (time() - t)))
    end
end

end # module


text:
using Roots

using UUIDs
using Observables
using CSV

using ScottishTaxBenefitModel
using .BCCalcs
using .Definitions
using .ExampleHelpers
using .FRSHouseholdGetter
using .GeneralTaxComponents
using .ModelHousehold
using .Monitor
using .Results
using .Runner
using .RunSettings
using .SimplePovertyCounts: GroupPoverty
using .SingleHouseholdCalculations
using .STBIncomes
using .STBOutput
using .STBParameters
using .Utils

const BASE_UUID = UUID("985c312f-129b-4acd-9e40-cb629d184183")

function initialise_settings()::Settings
    settings = Settings()
	settings.uuid = BASE_UUID
	settings.means_tested_routing = modelled_phase_in
    settings.run_name="run-$(date_string())"
	settings.income_data_source = ds_frs
	settings.dump_frames = false
	settings.do_marginal_rates = false
	settings.poverty_line = 100.0 # doesn't matter so long as +ive
	settings.requested_threads = 4
	settings.dump_frames = true
	return settings
end

const BASE_SETTINGS = initialise_settings()

function load_system()::TaxBenefitSystem
	sys = load_file( joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021_22.jl" ))
	#
	# Note that as of Budget21 removing these doesn't actually happen till May 2022.
	#
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "sys_2021-uplift-removed.jl"))
	# uc taper to 55
	load_file!( sys, joinpath( Definitions.MODEL_PARAMS_DIR, "budget_2021_uc_changes.jl"))
	weeklyise!( sys )
	return sys
end

settings = initialise_settings()

sys = load_system()
chsys = deepcopy( sys )
chsys.scottish_child_payment.amount = 20.0

params = [sys, chsys]

# results = do_one_run( settings, params, obs )

mutable struct Thing
    a :: Float64
end

mutable struct RunParameters{T<:AbstractFloat}
    params :: TaxBenefitSystem{T}
    settings :: Settings
	base_cost :: T
end

function run( x :: Number, things :: RunParameters )
    obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
    nsr = deepcopy( things.params.it.non_savings_rates )
    things.params.it.non_savings_rates .+= x
    results = do_one_run(things.settings, [things.params], obs )
    things.params.it.non_savings_rates = nsr
	summary = summarise_frames(results,settings)
	nc = summary.income_summary[1][1,:net_cost]
	return round( nc - things.base_cost, digits=0 )
    # x^2 + things.params.it.non_savings_rates[1]
end

function baserun_cost()
    obs = Observable(
		Progress(settings.uuid, "",0,0,0,0))
	results = do_one_run(settings, [sys], obs )
	summary = summarise_frames(results,settings)
	return summary.income_summary[1][1,:net_cost]
end

nc = baserun_cost()

things = RunParameters( chsys, settings, nc )

zerorun = ZeroProblem( run, 0.0 )

incch = solve( zerorun, things )

println( "incch=$incch")

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
"""
    DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulation)
"""
function DiscontinuousFESpace(reffes::Vector{<:ReferenceFE}, trian::Triangulation, is_ref=true)

  cell_to_ctype = get_cell_type(trian)
  cell_map = get_cell_map(trian)

  cell_dofs, nfree = compute_discontinuous_cell_dofs(reffes,cell_to_ctype)

  ndirichlet = 0
  dirichlet_dof_tag = Int8[]
  dirichlet_cells = Int[]
  ntags = 0

  cell_shapefuns, cell_dof_basis = compute_cell_space(reffes, cell_to_ctype, cell_map,Val(is_ref))

  UnconstrainedFESpace(
    nfree,
    ndirichlet,
    cell_dofs,
    cell_shapefuns,
    cell_dof_basis,
    dirichlet_dof_tag,
    dirichlet_cells,
    ntags)

end

"""
"""
function compute_discontinuous_cell_dofs(reffes,cell_type)

  ctype_to_nldofs = map(num_dofs,reffes)

  _compute_discontinuous_cell_dofs(cell_type,ctype_to_nldofs)

end

function _compute_discontinuous_cell_dofs(cell_to_ctype,ctype_to_nldofs)

  ncells = length(cell_to_ctype)
  ptrs = zeros(Int32,ncells+1)
  for (cell, ctype) in enumerate(cell_to_ctype)
    nldofs = ctype_to_nldofs[ctype]
    ptrs[cell+1] = nldofs
  end

  length_to_ptrs!(ptrs)

  ndata = ptrs[end]-1
  data = collect(Int,1:ndata)

  (Table(data,ptrs), ndata)

end


text:
__precompile__()

module QuantumWalk

using Expokit
using LightGraphs
using Optim

using LinearAlgebra
using SparseArrays
using Arpack

export
	AbstractCTQW,
	AbstractSzegedy,
	CTQW,
    CTQWDense,
	QSearchState,
	QWDynamics,
	QWEvolution,
	QWModel,
	QWModelCont,
	QWModelDiscr,
	QWSearch,
	QuantumWalk,
	Szegedy,
	check_qwdynamics,
	evolve,
	execute,
	execute_all,
	execute_all_measured,
	execute_single,
	execute_single_measured,
	expected_runtime,
	graph,
	initial_state,
	marked,
	matrix,
	maximize_quantum_search,
	measure,
	model,
	parameters,
	penalty,
	probability,
	runtime,
	sqrtstochastic,
	state

"""
    QuantumWalk

An package for quantum walk and quantum search simulation.

The packaged provides hierarchy of Quantum walk models and their numerical description together
with general function for quantum walk simulation and quantum search evolution. In
particular `CTQW` and `Szegedy` walk models are implemented as an example of continuous
and discrete walk models.

For both continuous and discrete model quantum walk simulator is implemented and
function finding optimal measurement time for quantum search function. For discrete
walk evolution various optimization may be chosen.

Full documentation is available in [GitHub pages](https://quantumwalks.github.io/QuantumWalk.jl/latest/).

"""
QuantumWalk

include("type_hierarchy.jl")
include("dynamics.jl")

include("qwdynamics/qwdynamics_import.jl")
include("qwmodels/qwmodels_import.jl")

end


text:
@testset "Test functionality of System" begin
    sys = create_rts_system()
    summary(devnull, sys)

    generators = collect(get_components(ThermalStandard, sys))
    generator = get_component(ThermalStandard, sys, get_name(generators[1]))
    @test IS.get_uuid(generator) == IS.get_uuid(generators[1])
    @test_throws(IS.ArgumentError, add_component!(sys, generator))

    generators2 = get_components_by_name(ThermalGen, sys, get_name(generators[1]))
    @test length(generators2) == 1
    @test IS.get_uuid(generators2[1]) == IS.get_uuid(generators[1])

    @test isnothing(get_component(ThermalStandard, sys, "not-a-name"))
    @test isempty(get_components_by_name(ThermalGen, sys, "not-a-name"))
    @test_throws(IS.ArgumentError,
                 get_component(ThermalGen, sys, "not-a-name"))
    @test_throws(IS.ArgumentError,
                 get_components_by_name(ThermalStandard, sys, "not-a-name"))

    # Test get_bus* functionality.
    bus_numbers = Vector{Int}()
    for bus in get_components(Bus, sys)
        push!(bus_numbers, bus.number)
        if length(bus_numbers) >= 2
            break
        end
    end

    bus = PowerSystems.get_bus(sys, bus_numbers[1])
    @test bus.number == bus_numbers[1]

    buses = PowerSystems.get_buses(sys, Set(bus_numbers))
    sort!(bus_numbers)
    sort!(buses, by=x -> x.number)
    @test length(bus_numbers) == length(buses)
    for (bus_number, bus) in zip(bus_numbers, buses)
        @test bus_number == bus.number
    end

    initial_times = get_forecast_initial_times(sys)
    @assert length(initial_times) == 1
    initial_time = initial_times[1]

    # Get forecasts with a label and without.
    components = get_components(HydroDispatch, sys)
    forecasts = get_forecasts(Forecast, sys, initial_time, components, "PMax MW")
    @test length(forecasts) > 0

    forecasts = collect(get_forecasts(Forecast, sys, initial_time, components))
    count = length(forecasts)
    @test count > 0

    # Verify that the two accessor functions return the same results.
    all_components = get_components(Component, sys)
    all_forecasts1 = get_forecasts(Forecast, sys, initial_time, all_components)
    all_forecasts2 = get_forecasts(Forecast, sys, initial_time)
    @test length(all_forecasts1) == length(all_forecasts2)

    # Get specific forecasts. They should not match.
    specific_forecasts = get_forecasts(PSY.Deterministic{Bus}, sys, initial_time)
    @test length(specific_forecasts) < length(all_forecasts1)

    @test get_forecasts_horizon(sys) == 24
    @test get_forecasts_initial_time(sys) == Dates.DateTime("2020-01-01T00:00:00")
    @test get_forecasts_interval(sys) == Dates.Hour(1)  # TODO
    @test get_forecasts_resolution(sys) == Dates.Hour(1)  # TODO

    for forecast in collect(get_forecasts(Forecast, sys, initial_time))
        remove_forecast!(sys, forecast)
    end

    @test length(get_forecasts(Forecast, sys, initial_time)) == 0

    # ArgumentError is thrown if the type is concrete and there is no forecast for a
    # component.
    @test_throws(IS.ArgumentError,
                 get_forecasts(Forecast, sys, initial_time, components))

    # But not if the type is abstract.
    new_forecasts = get_forecasts(Forecast, sys, initial_time,
                                  get_components(HydroGen, sys))
    @test length(new_forecasts) == 0

    add_forecasts!(sys, forecasts)

    forecasts = get_forecasts(Forecast, sys, initial_time, components)
    @assert length(forecasts) == count

    @test_throws(IS.ArgumentError,
                 get_forecasts(PSY.Deterministic{Bus}, sys, initial_time, components))

    #Get forecast by type
    res = get_component_forecasts(RenewableDispatch, sys, initial_time)
    for i in res
        @test isa(i,PSY.Deterministic{RenewableDispatch})
    end

    @test_throws(IS.ArgumentError,
                 get_component_forecasts(RenewableGen, sys, initial_time))

    f = forecasts[1]
    forecast = PSY.Deterministic(Bus(nothing), f.label, f.resolution, f.initial_time, f.data)
    @test_throws(IS.ArgumentError, add_forecasts!(sys, [forecast]))

    component = deepcopy(f.component)
    component.internal = IS.InfrastructureSystemsInternal()
    forecast = PSY.Deterministic(component, f.label, f.resolution, f.initial_time, f.data)
    @test_throws(IS.ArgumentError, add_forecasts!(sys, [forecast]))
end

@testset "Test System iterators" begin
    sys = create_rts_system()

    i = 0
    for component in iterate_components(sys)
        i += 1
    end

    components = get_components(Component, sys)
    @test i == length(components)

    i = 0
    for forecast in iterate_forecasts(sys)
        i += 1
    end

    j = 0
    initial_times = get_forecast_initial_times(sys)
    for initial_time in initial_times
        j += length(get_forecasts(Forecast, sys, initial_time))
    end

    @test i == j
end

@testset "Test remove_component" begin
    sys = create_rts_system()
    generators = get_components(ThermalStandard, sys)
    initial_length = length(generators)
    @assert initial_length > 0
    gen = collect(generators)[1]

    remove_component!(sys, gen)

    @test isnothing(get_component(typeof(gen), sys, get_name(gen)))
    generators = get_components(typeof(gen), sys)
    @test length(generators) == initial_length - 1

    @test_throws(IS.ArgumentError, remove_component!(sys, gen))

    add_component!(sys, gen)
    remove_component!(typeof(gen), sys, get_name(gen))
    @test isnothing(get_component(typeof(gen), sys, get_name(gen)))

    @assert length(get_components(typeof(gen), sys)) > 0
    remove_components!(typeof(gen), sys)
    @test_throws(IS.ArgumentError, remove_components!(typeof(gen), sys))
end

@testset "Test missing Arc bus" begin
    sys = System(100)
    line = Line(nothing)
    @test_throws(IS.ArgumentError, add_component!(sys, line))
end


text:
module TelegramAPI

using HTTP, JSON
include("api.jl")

export get_me, get_updates, send_message

end # module


text:
"""
Simple model with 3 variables, 3 shocks and 3 transition equations.
Like E2, but with more lags/leads.

The NL version is artificially made non-linear by taking one of the
equations (pinf[t]=rhs...) and rewriting it as exp(pinf[t])=exp(rhs...).
The solution of E3nl is identical to the linear E3 and can be used for
testing of the linearization solver.
"""
module  E3nl

using ModelBaseEcon

# start with an empty model
model = Model()
model.options.maxiter = 200

# add parameters
@parameters model begin
    cp = [0.5, 0.02]
    cr = [0.75, 1.5, 0.5]
    cy = [0.5, -0.02]
end

# add variables: a list of symbols
@variables model begin
    pinf
    rate
    ygap
end

# add shocks: a list of symbols
@shocks model begin
    pinf_shk
    rate_shk
    ygap_shk
end

# autoexogenize: define a mapping of variables to shocks
@autoexogenize model begin
    pinf = pinf_shk
    rate = rate_shk
    ygap = ygap_shk
end

# add equations: a sequence of expressions, such that
# use y[t+1] for expectations/leads
# use y[t] for contemporaneous
# use y[t-1] for lags
# each expression must have exactly one "="
@equations model begin
    exp(pinf[t])=exp(cp[1]*pinf[t-1]+0.3*pinf[t+1]+0.05*pinf[t+2]+0.05*pinf[t+3]+cp[2]*ygap[t]+pinf_shk[t])
    rate[t]=cr[1]*rate[t-1]+(1-cr[1])*(cr[2]*pinf[t]+cr[3]*ygap[t])+rate_shk[t]
    ygap[t]=cy[1]/2*ygap[t-2]+cy[1]/2*ygap[t-1]+(.98-cy[1])*ygap[t+1]+cy[2]*(rate[t]-pinf[t+1])+ygap_shk[t]
end

# call initialize! to build internal structures
@initialize model

end


text:
using Documenter, GoogleSheets

makedocs(
    modules = [GoogleSheets],
    sitename="GoogleSheets.jl",
    authors = "Chip Kent",
    format = Documenter.HTML(),
)

deploydocs(
    repo = "github.com/chipkent/GoogleSheets.jl.git",
    devbranch = "main",
    push_preview = true,
)


text:
using HypothesisTests, Test

@testset "Clark-West tests" begin

    e1 = [ 9.78432,  12.73500,   8.67224,   2.62740,   5.60947,
           7.57809,   4.53774,   2.51084,   0.49290,   3.48394,
           9.46152,   9.41220,   2.36289,   0.34495,   3.33599,
           5.31357,   4.28219,  -3.74471,  -5.73575,  -3.71781 ]

    e2 =  [ 2.82053,   4.39754,  -1.78647,  -4.30662,   3.69526,
            3.37259,  -1.09002,  -0.50984,  -0.78973,   3.89077,
            7.52849,   2.82373,  -3.95838,  -0.13606,   4.54444,
            4.18216,   1.67993,  -5.38077,  -0.85686,   2.70479 ]

    atol = 1e-4
    cw_test = ClarkWestTest(e1, e2)
    @test pvalue(cw_test) ≈ 0.0002 atol=atol
    @test pvalue(cw_test, tail=:right) ≈ 0.0001 atol=atol

    cw_test = ClarkWestTest(e1, e2, 3)
    @test pvalue(cw_test) ≈ 0.0157 atol=atol
    @test pvalue(cw_test, tail=:right) ≈ 0.0079 atol=atol

    show(IOBuffer(), cw_test)

    @test_throws DimensionMismatch ClarkWestTest(rand(3), rand(4))

end


text:
using StatsModels:
    FullRank,
    Schema,
    drop_intercept,
    implicit_intercept,
    hasintercept,
    omitsintercept,
    collect_matrix_terms

struct MultiSchema{S}
    base::S
    subs::Dict{Any,S}
end

MultiSchema(s::S) where {S} = MultiSchema(s, Dict{Any,S}())

function StatsModels.apply_schema(
    t::StatsModels.AbstractTerm,
    sch::MultiSchema,
    Ctx::Type
)
    apply_schema(t, sch.base, Ctx)
end

function StatsModels.apply_schema(
    t::StatsModels.TupleTerm,
    sch::MultiSchema,
    Ctx::Type
)
    sum(apply_schema.(t, Ref(sch), Ref(Ctx)))
end

# copied with minimal modifications from StatsModels.jl, in order to wrap the schema
# in MultiSchema.
function StatsModels.apply_schema(t::FormulaTerm, schema::Schema, Mod::Type{<:MixedModel})
    schema = FullRank(schema)

    # Models with the drop_intercept trait do not support intercept terms,
    # usually because they include one implicitly.
    if drop_intercept(Mod)
        if hasintercept(t)
            throw(ArgumentError("Model type $Mod doesn't support intercept " *
                                "specified in formula $t"))
        end
        # start parsing as if we already had the intercept
        push!(schema.already, InterceptTerm{true}())
    elseif implicit_intercept(Mod) && !hasintercept(t) && !omitsintercept(t)
        t = FormulaTerm(t.lhs, InterceptTerm{true}() + t.rhs)
    end

    # only apply rank-promoting logic to RIGHT hand side
    FormulaTerm(apply_schema(t.lhs, schema.schema, Mod),
                collect_matrix_terms(apply_schema(t.rhs, MultiSchema(schema), Mod)))
end


text:
using Documenter, VarBCFiles

makedocs(sitename="VarBCFiles")

deploydocs(
    repo = "github.com/Hirlam/VarBCFiles.jl.git",
    devbranch = "main"
)


text:
## Show the results
function Base.show(io::IO, m::Results)
  compact = get(io, :compact, false)
  if !compact
    println("Fit converged after $(m.iterations) iterations")
    println("Best fit at: -LnL = $(m.min_objective)")
    println("Values of $(m.min_parameters)")
  end
end

## Pretty output
function pretty_results(r::Results)
  df = DataFrame(
                 Name = String[],
                 Fit = Float64[],
                 Interval_Low = Float64[],
                 Interval_High = Float64[],
                )
  for v in r.model.params
    push!( df,
          (
           v.name,
           v.fit,
           v.low,
           v.high,
          )
         )
  end
  return df
end


text:
@testset "Gauging type" begin

    # Test the construction of the couple (level, discharge)
    G = Gauging(0.0,1.0)
    @test G.level ≈ 0.0 atol = sqrt(eps())
    @test G.discharge ≈ 1.0 atol = sqrt(eps())

    # Test the promotion of level and discharge
    G = Gauging(0,1)
    @test typeof(G.level) == Float64
    @test typeof(G.discharge) == Float64

    # Test the promotion of level
    G = Gauging(0,1.1)
    @test typeof(G.level) == Float64
    @test typeof(G.discharge) == Float64

    # Test the empty gauing construction
    G = Gauging()
    @test G.level ≈ 0.0 atol = sqrt(eps())
    @test G.discharge ≈ 0.0 atol = sqrt(eps())

    # Test the broadcast
    G = Gauging.([1,2],[3,4])
    @test G[2].level ≈ 2
    @test G[2].discharge ≈ 4

end

@testset "Gauging methods" begin

    # Test `level()`for a single gauging
    G = Gauging(0.0,1.0)
    @test level(G) ≈ 0.0

    # Test the broadcasting of `level()`
    G = Gauging.([0.0,1.0],[2.0,3.0])
    @test level.(G) ≈ [0.0, 1.0]

    # Test `discharge()`for a single gauging
    G = Gauging(0.0,1.0)
    @test discharge(G) ≈ 1.0

    # Test the broadcasting of `discharge()`
    G = Gauging.([0.0,1.0],[2.0,3.0])
    @test discharge.(G) ≈ [2.0, 3.0]

end

@testset "bootstrap Gauging" begin
    G = Gauging.(1:3, 1:3)

    Gbs = RatingCurves.bootstrap(G)

    @test minimum(G) ∈ Gbs
    @test length(G) == length(Gbs)

end

@testset "crcfit" begin
    G = Gauging.([1/4, 1/3, 1/2, 2,3,4],[1/32, 2/27, 1/4, 8, 18, 32])
    crc = crcfit(G)

    rc₁ = crc.component[1]
    rc₂ = crc.component[2]

    @test crc.threshold[] ≈ 1.0

    @test rc₁.a ≈ 2.0
    @test rc₁.b ≈ 0.0 atol=1e-8
    @test rc₁.c ≈ 3.0

    @test rc₂.a ≈ 2.0
    @test rc₂.b ≈ 0.0 atol=1e-8
    @test rc₂.c ≈ 2.0
end

@testset "crcfit with given k" begin
    G = Gauging.([1/4, 1/3, 1/2, 2,3,4],[1/32, 2/27, 1/4, 8, 18, 32])
    crc = crcfit(G, 1)

    rc₁ = crc.component[1]
    rc₂ = crc.component[2]

    @test crc.threshold[] ≈ 1.0

    @test rc₁.a ≈ 2.0
    @test rc₁.b ≈ 0.0 atol=1e-8
    @test rc₁.c ≈ 3.0

    @test rc₂.a ≈ 2.0
    @test rc₂.b ≈ 0.0 atol=1e-8
    @test rc₂.c ≈ 2.0
end

@testset "getinitialvalues" begin

    G = Gauging.([1,2,3], [4,5,6])
    rc = RatingCurves.getinitialvalues(G)

    @test rc.a ≈ 126.4911 atol = 1e-4
    @test rc.b ≈ .9
    @test rc.c ≈ 3/2

end

@testset "minimum of Gauging" begin
    G = Gauging.(1:3, 1:3)

    @test minimum(G).level ≈ 1.0
    @test minimum(G).discharge ≈ 1.0

end

@testset "maximum of Gauging" begin
    G = Gauging.(1:3, 1:3)

    @test maximum(G).level ≈ 3.0
    @test maximum(G).discharge ≈ 3.0

end

@testset "rcfit with given b" begin
    G = Gauging.([2,3,4],[3,12,27])
    rc = rcfit(G, 1)
    @test rc.a ≈ 3.0
    @test rc.b ≈ 1.0
    @test rc.c ≈ 2.0
end

@testset "rcfit" begin
    G = Gauging.([2,3,4],[3,12,27])
    rc = rcfit(G)
    @test rc.a ≈ 3.0
    @test rc.b ≈ 1.0
    @test rc.c ≈ 2.0
end

@testset "rcfit with given b s.t. constraint" begin
    G = Gauging.([2,3,4],[3,12,27])
    rc = rcfit(G, 1, [5, 48])
    @test rc.a ≈ 3.0
    @test rc.b ≈ 1.0
    @test rc.c ≈ 2.0
end

@testset "rcfit s.t. constraint" begin
    G = Gauging.([2,3,4],[3,12,27])
    rc = rcfit(G, [5, 48])
    @test rc.a ≈ 3.0
    @test rc.b ≈ 1.0
    @test rc.c ≈ 2.0
end

@testset "sort Gauging" begin
    G = Gauging.(1:3, 1:3)

    G_ascend = sort(G)
    @test level.(G_ascend) ≈ [1, 2, 3]

    G_descend = sort(G, rev=true)
    @test level.(G_descend) ≈ [3, 2, 1]

end


text:
using BinDeps

@BinDeps.setup

ENV["JULIA_ROOT"] = abspath(JULIA_HOME, "../../")

libportaudio = library_dependency("libportaudio")
libsndfile = library_dependency("libsndfile")

# TODO: add other providers with correct names
provides(AptGet, "portaudio19-dev", libportaudio)
provides(AptGet, "libsndfile1-dev", libsndfile)
provides(Pacman, "portaudio", libportaudio)
provides(Pacman, "libsndfile", libsndfile)

@osx_only begin
    using Homebrew
    provides(Homebrew.HB, "portaudio", libportaudio)
    provides(Homebrew.HB, "libsndfile", libsndfile)
end

@windows_only begin
    using WinRPM
    provides(WinRPM.RPM, "libportaudio2", libportaudio, os = :Windows)
    provides(WinRPM.RPM, "libsndfile1", libsndfile, os = :Windows)
end

@BinDeps.install [:libportaudio => :libportaudio,
                  :libsndfile => :libsndfile]


text:
using Documenter, Qaintessent

makedocs(
    sitename="Qaintessent.jl Documentation",
    pages = [
        "Home" => "index.md",
        "Section" => [
            "gates.md",
            "circuitgates.md",
            "circuit.md",
            "densitymatrices.md",
            "gradients.md",
            "view.md",
            "qasm.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/Qaintum/Qaintessent.jl.git",
    push_preview = true
)


text:
using PrimeFactorsKata
using Documenter

makedocs(;
    modules=[PrimeFactorsKata],
    authors="David Wilde",
    repo="https://github.com/davidchris/PrimeFactorsKata.jl/blob/{commit}{path}#L{line}",
    sitename="PrimeFactorsKata.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://davidchris.github.io/PrimeFactorsKata.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/davidchris/PrimeFactorsKata.jl",
)


text:
abstract type HamiltonianModel end

@doc raw"
    hamiltonian(model<:HamiltonianModel)

return the hamiltonian of the `model` as a two-site tensor operator.
"
function hamiltonian end

struct Ising <: HamiltonianModel end

@doc raw"
    TFIsing(hx::Real)

return a struct representing the transverse field ising model with magnetisation `hx`.
"
struct TFIsing{T<:Real} <: HamiltonianModel
    hx::T
end

"""
    hamiltonian(model::TFIsing)

return the transverse field ising hamiltonian for the provided `model` as a
two-site operator.
"""
function hamiltonian(model::TFIsing)
    hx = model.hx
    -2 * ein"ij,kl -> ijkl"(σz,σz) -
        hx/2 * ein"ij,kl -> ijkl"(σx, id2) -
        hx/2 * ein"ij,kl -> ijkl"(id2, σx)
end

@doc raw"
    Heisenberg(Jz::T,Jx::T,Jy::T) where {T<:Real}

return a struct representing the heisenberg model with magnetisation fields
`Jz`, `Jx` and `Jy`..
"
struct Heisenberg{T<:Real} <: HamiltonianModel
    Jz::T
    Jx::T
    Jy::T
end
Heisenberg() = Heisenberg(1.0,1.0,1.0)

"""
    hamiltonian(model::Heisenberg)

return the heisenberg hamiltonian for the `model` as a two-site operator.
"""
function hamiltonian(model::Heisenberg)
    h = model.Jz * ein"ij,kl -> ijkl"(σz,σz) -
        model.Jx * ein"ij,kl -> ijkl"(σx, σx) -
        model.Jy * ein"ij,kl -> ijkl"(σy, σy)
    h = ein"ijcd,kc,ld -> ijkl"(h,σx,σx')
    real(h ./ 2)
end


text:
using DelaySSAToolkit
using Catalyst, DiffEqJump
using Statistics
using Test

reltol = 0.1
@parameters t
@variables X(t)
burst_sup = 30
a, b = [0.0282, 3.46]
rxs = [Reaction(a*b^i/(1+b)^(i+1),nothing,[X],nothing,[i]) for i in 1:burst_sup]
rxs = vcat(rxs)
@named rs = ReactionSystem(rxs,t,[X],[])

# convert the ReactionSystem to a JumpSystem
jumpsys = convert(JumpSystem, rs, combinatoric_ratelaws=false)

u0 = [0]
de_chan0 = [[]]
tf = 200.
tspan = (0,tf)
dprob = DiscreteProblem(jumpsys,u0,tspan)
τ = 130.

delay_trigger = Dict([Pair(i, [1=>fill(τ, i)]) for i in 1:burst_sup])
delay_complete = Dict(1=>[1=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)

algos = [DelayDirect(), DelayMNRM(), DelayRejection(), DelayDirectCR()]
timestamps = [10, 20, 50, 200]

bursty_mean(t) = a*b*min(t,τ)
bursty_var(t) = 2*a*b^2*min(t,τ) + a*b*min(t,τ)

# Check mean and variance
samplesize = Int64(5e4)
for alg in algos
    jprob = DelayJumpProblem(jumpsys, dprob, alg, delayjumpset, de_chan0, save_positions=(false, false))
    ensprob = EnsembleProblem(jprob)
    @time ens = solve(ensprob, SSAStepper(), EnsembleSerial(), trajectories=samplesize, saveat = timestamps)
    for idx in eachindex(timestamps)
    sol_end = reduce(vcat, [ens[i].u[idx+1] for i in 1:samplesize])
        t_ = timestamps[idx]
        @test mean(sol_end) ≈ bursty_mean(t_) atol = reltol*bursty_mean(t_)
        @test var(sol_end) ≈ bursty_var(t_) atol = reltol*bursty_var(t_)
    end
end


text:
## Run Metropolis-Hastings using a Bayesian logit model with a Normal prior
## N(0, priorVar*I) on the Swiss banknote data

using Distributions, GeometricMCMC

include("swiss.jl")
include("logitNormalPrior.jl")

# Create Model instance
model =
  Model(nPars, data, logPrior, logLikelihood, gradLogPosterior, randPrior);

# Create instance mhOpts of Metropolis-Hastings options
mhOpts = MhOpts(55000, 5000, 0.1);

# Run Metropolis-Hastings simulation
mcmc, z = mh(model, mhOpts);

# Compute ZV-MH mean estimators based on linear polynomial
linearZvMcmc, linearCoef = linearZv(mcmc, z);

# Compute ZV-MH mean estimators based on quadratic polynomial
quadraticZvMcmc, quadraticCoef = quadraticZv(mcmc, z);


text:
using Revise
using Jtac

G = MetaTac
bp = 0.02

chain = @chain G Conv(256) Pool() Dropout(0.2) Deconv(512) Conv(512) Dropout(0.2) Conv(256)
#chain = @chain G Dense(1000) Dense(1000) Dense(1000) Dense(1000) Dense(Jtac.policy_length(G)+1)

#model = ShallowConv(G, 500)
model = NeuralModel(G, chain)

amodel = Async(model, max_batchsize = 100)

games = [G() for i in 1:100]

println("CPU performance for 1000 games")
@time model(games);
@time map(model, games);
@time asyncmap(amodel, games);

@time record_self(MCTSPlayer(model, power = 100), branching = bp)
#@time record_selfplay(amodel, 20, branch_prob = bp)

gmodel = model |>  to_gpu
gamodel = Async(gmodel, max_batchsize = 100)

games = [MetaTac() for i in 1:1000]
#games = map(1:100) do _
#  game = G()
#  for i in 1:rand(1:20)
#    if !Jtac.is_over(game)
#      Jtac.random_turn!(game)
#    end
#  end
#  game
#end

println()
println("GPU performance for 1000 games")
@time gmodel(games)
@time map(gmodel, games)
@time asyncmap(gamodel, games)

@time record_self(MCTSPlayer(gmodel, power = 100), branching = bp)
@time record_self(MCTSPlayer(gamodel, power = 100), 20, branching = bp)

max_stars_repo_name: roSievers/jtac


text:
## delete this on NCI
# using Distributed
# addprocs(2)
## end delete
## MPI Init
using MPIClusterManagers, Distributed
import MPI
MPI.Init()
rank = MPI.Comm_rank(MPI.COMM_WORLD)
size = MPI.Comm_size(MPI.COMM_WORLD)
if rank == 0
    @info "size is $size"
end
manager = MPIClusterManagers.start_main_loop(MPI_TRANSPORT_ALL)
@info "there are $(nworkers()) workers"
@everywhere @info gethostname()
include("01_read_data_aseggdf.jl")
include("02_set_options.jl")
## MPI checks
# split into sequential iterations of parallel soundings
nsoundings = length(soundings)
ncores = nworkers()
nsequentialiters = ceil(Int, nsoundings/ncores)
@info "will require $nsequentialiters iterations of $ncores soundings"
## set up transD_GP
@everywhere using Distributed
@everywhere using HiQGA.transD_GP
## do the parallel soundings
@info "starting"
transD_GP.SkyTEM1DInversion.loopacrosssoundings(soundings, σstart, σ0,
                            nsequentialiters   = nsequentialiters,
                            zfixed             = zfixed          ,
                            ρfixed             = ρfixed          ,
                            zstart             = zstart          ,
                            extendfrac         = extendfrac      ,
                            dz                 = dz              ,
                            ρbg                = ρbg             ,
                            nlayers            = nlayers         ,
                            ntimesperdecade    = ntimesperdecade ,
                            nfreqsperdecade    = nfreqsperdecade ,
                            modelprimary       = modelprimary    ,
                            regtype            = regtype         ,
                            nstepsmax          = nstepsmax       ,
                            ntries             = ntries          ,
                            target             = target          ,
                            lo                 = lo              ,
                            hi                 = hi              ,
                            λ²min              = λ²min           ,
                            λ²max              = λ²max           ,
                            λ²frac             = λ²frac          ,
                            β²                 = β²              ,
                            ntestdivsλ²        = ntestdivsλ²     ,
                            αmin               = αmin            ,
                            αmax               = αmax            ,
                            αfrac              = αfrac           ,
                            ntestdivsα         = ntestdivsα      ,
                            regularizeupdate   = regularizeupdate,
                            knownvalue         = knownvalue      ,
                            firstvalue         = firstvalue      ,
                            κ                  = κ               ,
                            breakonknown       = breakonknown    ,
                            dobo               = dobo,
                            zipsaveprefix      = zipsaveprefix)

MPIClusterManagers.stop_main_loop(manager)
rmprocs(workers())
exit()

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
#
# This example implements a DFA (deterministiv finite automaton)
#
#   A = ({s,q,p,r},{'a','b','c'},δ,s,{r})
#
#   with the following transition table:

#         | input symbols
#  states | 'a' | 'b' | 'c'
#  -------+-----+-----+-----
#   -> s  |  s  |  q  |  s
#      q  |  p  |  r  |  r
#      p  |  p  |  q  |  p
#    * r  |  r  |  r  |  r
#
using YAActL

s(::Val{'b'}) = become(q)
s(x) = nothing
q(::Val{'a'}) = become(p)
q(::Val{'b'}) = become(r)
q(::Val{'c'}) = become(r)
q(x) = nothing
p(::Val{'b'}) = become(q)
p(x) = nothing
r(x) = nothing

mydfa = Actor(s)

function check(lk::Link, str::String)
    become!(lk, s)          # start behavior
    foreach(str) do c
        cast!(lk, Val(c))   # send all chars to the actor
    end
    query!(lk, :bhv) == r   # check the actors behavior
end

check(mydfa, "baab")        # false
check(mydfa, "baabc")       # true
check(mydfa, "babaccabb")   # true


text:
module NilNumber

using NNlib
import Random

"""
    Nil <: Number

Nil is a singleton type with a single instance `nil`.
Unlike `Nothing` and `Missing` it subtypes `Number`.
"""
struct Nil <: Number end

const nil = Nil()

Nil(::T) where T<:Number = nil
(::Type{T})(::Nil) where T<:Number = nil
Base.convert(::Type{Nil}, ::Number) = nil

Base.float(::Type{Nil}) = Nil

for f in [:copy, :zero, :one, :oneunit,
          :+, :-, :abs, :abs2, :inv,
          :exp, :log, :log1p, :log2, :log10,
          :sqrt, :tanh, :conj]
  @eval Base.$f(::Nil) = nil
end

for f in [:+, :-, :*, :/, :^, :mod, :div, :rem]
  @eval Base.$f(::Nil, ::Nil) = nil
end

Base.isless(::Nil, ::Nil) = true
Base.isless(::Nil, ::Number) = true
Base.isless(::Number, ::Nil) = true

Base.isnan(::Nil) = false

Base.typemin(::Type{Nil}) = nil
Base.typemax(::Type{Nil}) = nil

Base.promote_rule(x::Type{Nil}, y::Type{<:Number}) = Nil

Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{Nil}) = nil

end  # module

using .NilNumber: Nil, nil

"""
    outputsize(m, inputsize::Tuple; padbatch=false)

Calculate the output size of model `m` given the input size.
Obeys `outputsize(m, size(x)) == size(m(x))` for valid input `x`.
Keyword `padbatch=true` is equivalent to using `(inputsize..., 1)`, and
returns the final size including this extra batch dimension.

This should be faster than calling `size(m(x))`. It uses a trivial number type,
and thus should work out of the box for custom layers.

If `m` is a `Tuple` or `Vector`, its elements are applied in sequence, like `Chain(m...)`.

# Examples
```jldoctest
julia> using Flux: outputsize

julia> outputsize(Dense(10, 4), (10,); padbatch=true)
(4, 1)

julia> m = Chain(Conv((3, 3), 3 => 16), Conv((3, 3), 16 => 32));

julia> m(randn(Float32, 10, 10, 3, 64)) |> size
(6, 6, 32, 64)

julia> outputsize(m, (10, 10, 3); padbatch=true)
(6, 6, 32, 1)

julia> outputsize(m, (10, 10, 3, 64))
(6, 6, 32, 64)

julia> try outputsize(m, (10, 10, 7, 64)) catch e println(e) end
DimensionMismatch("Input channels must match! (7 vs. 3)")

julia> outputsize([Dense(10, 4), Dense(4, 2)], (10, 1))
(2, 1)

julia> using LinearAlgebra: norm

julia> f(x) = x ./ norm.(eachcol(x));

julia> outputsize(f, (10, 1)) # manually specify batch size as 1
(10, 1)

julia> outputsize(f, (10,); padbatch=true) # no need to mention batch size
(10, 1)
```
"""
function outputsize(m, inputsize::Tuple; padbatch=false)
  inputsize = padbatch ? (inputsize..., 1) : inputsize

  return size(m(fill(nil, inputsize)))
end

## make tuples and vectors be like Chains

outputsize(m::Tuple, inputsize::Tuple; padbatch=false) = outputsize(Chain(m...), inputsize; padbatch=padbatch)
outputsize(m::AbstractVector, inputsize::Tuple; padbatch=false) = outputsize(Chain(m...), inputsize; padbatch=padbatch)

## bypass statistics in normalization layers

for layer in (:LayerNorm, :BatchNorm, :InstanceNorm, :GroupNorm)
  @eval (l::$layer)(x::AbstractArray{Nil}) = x
end

## fixes for layers that don't work out of the box

for (fn, Dims) in ((:conv, DenseConvDims), (:depthwiseconv, DepthwiseConvDims))
  @eval begin
    function NNlib.$fn(a::AbstractArray{Nil}, b::AbstractArray{Nil}, dims::$Dims)
      fill(nil, NNlib.output_size(dims)..., NNlib.channels_out(dims), size(a)[end])
    end

    function NNlib.$fn(a::AbstractArray{<:Real}, b::AbstractArray{Nil}, dims::$Dims)
      NNlib.$fn(fill(nil, size(a)), b, dims)
    end

    function NNlib.$fn(a::AbstractArray{Nil}, b::AbstractArray{<:Real}, dims::$Dims)
      NNlib.$fn(a, fill(nil, size(b)), dims)
    end
  end
end


text:
using PiGPIO

red_pin = 23
green_pin = 24
blue_pin = 25
p=Pi()

using Colors

dc = distinguishable_colors(10)

for r in dc
    set_PWM_dutycycle(p, red_pin, round(Int, r.r*255))
    set_PWM_dutycycle(p, green_pin, round(Int,r.g*255))
    set_PWM_dutycycle(p, blue_pin, round(Int,r.b*255))
    sleep(0.5)
end

set_PWM_dutycycle(p, red_pin, 0)
set_PWM_dutycycle(p, green_pin, 0)
set_PWM_dutycycle(p, blue_pin, 0)

using SIUnits
using SIUnits.ShortUnits

function normalize0(c::XYZ)
    d=convert(xyY, c)
    xyY(d.x, d.y, 1.0)
end

const hc_k  = 0.0143877696*K*m
const twohc²= 1.19104287e-16*Watt*m^2
planck{S1<:Real, S2<:Real}(λ::quantity(S1,Meter); T::quantity(S2,Kelvin)=5778.0K) =
    λ≤0m ? zero(λ)*Watt*m^-4 : twohc²*λ^-5.0/(exp(hc_k/(λ*T))-1)

Base.convert{S<:Real}(::Type{xyY}, T::quantity(S, K)) =
  mapreduce(λ->planck(λ*nm,T=T)*m^3/Watt*colormatch(λ), +, 380:780) |>
  normalize0

blackbodies = xyY[convert(xyY, T) for T in 100K:200K:10000K]

for b in blackbodies
     r=convert(RGB, b)
     set_PWM_dutycycle(p, red_pin, round(Int, r.r*255))
     set_PWM_dutycycle(p, green_pin, round(Int,r.g*255))
     set_PWM_dutycycle(p, blue_pin, round(Int,r.b*255))
     sleep(0.5)
end


text:
export ChannelCreate,
    ChannelUpdate,
    ChannelDelete,
    ChannelPinsUpdate

"""
Sent when a new [`DiscordChannel`](@ref) is created.
"""
struct ChannelCreate <: AbstractEvent
    channel::DiscordChannel
end
@boilerplate ChannelCreate :docs :mock
ChannelCreate(; kwargs...) = ChannelCreate(DiscordChannel(; kwargs...))
ChannelCreate(d::Dict{Symbol, Any}) = ChannelCreate(; d...)

"""
Sent when a [`DiscordChannel`](@ref) is updated.
"""
struct ChannelUpdate <: AbstractEvent
    channel::DiscordChannel
end
@boilerplate ChannelUpdate :docs :mock
ChannelUpdate(; kwargs...) = ChannelUpdate(DiscordChannel(; kwargs...))
ChannelUpdate(d::Dict{Symbol, Any}) = ChannelUpdate(; d...)

"""
Sent when a [`DiscordChannel`](@ref) is deleted.
"""
struct ChannelDelete <: AbstractEvent
    channel::DiscordChannel
end
@boilerplate ChannelDelete :docs :mock
ChannelDelete(; kwargs...) = ChannelDelete(DiscordChannel(; kwargs...))
ChannelDelete(d::Dict{Symbol, Any}) = ChannelDelete(; d...)

"""
Sent when a [`DiscordChannel`](@ref)'s pins are updated.
"""
struct ChannelPinsUpdate <: AbstractEvent
    channel_id::Snowflake
    last_pin_timestamp::Nullable{DateTime}
end
@boilerplate ChannelPinsUpdate :constructors :docs :mock


text:
function print_load_info(cur_cabinet::FileCabinet)

  println("\nCurrent load order:\n")

  for cur_file in cur_cabinet.load_order
    println(" + ", cur_file)
  end

  if !isempty(cur_cabinet.export_list)

    println("\nExported the following undefined variables:\n")

    for cur_export in cur_cabinet.export_list
      println(" + ", cur_export)
    end

  end

end


text:
; equals_test
import 'lib_src/print.jl'

class Animal {
    field name, age
    def new(name, age) {
        self.name = name
        self.age = age
    }

    def eq(other) {
        return self.name == other.name & self.age == other.age
    }
}

class Dog : Animal {
    field breed
    def new(name, age, breed) : (name, age) {
        self.breed = breed
    }
}

class Cat : Animal {
    field size
    def new(name, age, size) : (name, age) {
        self.size = size
    }
}

def test_is(obj, cls) {
    print_strln(if obj is cls 'Success' else 'Failure')
}

def test_isnt(obj, cls) {
    print_strln(if obj isnt cls 'Success' else 'Failure')
}

def main {
    dog1 = Dog.new('dog1', 4, 'Daschund')
    cat1 = Cat.new('cat1', 1, 3)

    test_is(dog1, Dog)
    test_is(dog1, Animal)

    test_is(cat1, Cat)
    test_is(cat1, Animal)

    test_isnt(dog1, Cat)
    test_isnt(cat1, Dog)

    println(hash_uint32_t_(5))
    println(hash_uint32_t_(7.2))
    println(hash_uint32_t_('abcdef'))
    println(hash_uint32_t_([1, 2, 3, 4, 5, 6]))
    println(hash_uint32_t_([1, 2, 3, 4, 5]))
    println(hash_uint32_t_([1, 2, 3, 4]))
    println(hash_uint32_t_([1, 2, 3]))
    println(hash_uint32_t_([1, 2]))
    println(hash_uint32_t_([1]))
    println(hash_uint32_t_([]))

}

max_stars_repo_name: jeffreymanzione/jeff-vm-lang

max_issues_repo_name: jeffreymanzione/jeff-vm-lang

max_forks_repo_name: jeffreymanzione/jeff-vm-lang


text:
"""
TestSolution

"""
mutable struct TestSolution{T,N,hasinterp,tType,uType,iType} <: AbstractTimeseriesSolution{T,N,uType}
  t::tType
  u::uType
  interp::iType
  dense::Bool
  retcode::Symbol
end
(T::TestSolution)(t) = T.interp(t)
function TestSolution(t,u)
  T = eltype(eltype(u))
  N = length((size(u[1])..., length(u)))
  TestSolution{T,N,false,typeof(t),typeof(u),Nothing}(t,u,nothing,false,:Success)
end
function TestSolution(t,u,interp)
  T = eltype(eltype(u))
  N = length((size(u[1])..., length(u)))
  TestSolution{T,N,true,typeof(t),typeof(u),typeof(interp)}(t,u,interp,true,:Success)
end
TestSolution(interp::DESolution) = TestSolution{Nothing,0,true,Nothing,Nothing,typeof(interp)}(nothing,nothing,interp,true,:Success)
hasinterp(::TestSolution{T,N,hi,tType,uType,iType}) where {T,N,hi,tType,uType,iType} = hi
"""
`appxtrue(sol::AbstractODESolution,sol2::TestSolution)`

Uses the interpolant from the higher order solution sol2 to approximate
errors for sol. If sol2 has no interpolant, only the final error is
calculated.
"""
function appxtrue(sol::AbstractODESolution,sol2::TestSolution)
  if sol2.u == nothing && hasinterp(sol2)
    _sol = TestSolution(sol.t,sol2(sol.t),sol2)
  else
    _sol = sol2
  end

  errors = Dict(:final=>recursive_mean(abs.(sol[end]-_sol[end])))
  if _sol.dense
    timeseries_analytic = _sol(sol.t)
    errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol-timeseries_analytic))
    errors[:l2] = sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol-timeseries_analytic)))
    densetimes = collect(range(sol.t[1],stop=sol.t[end],length=100))
    interp_u = sol(densetimes)
    interp_analytic = _sol(densetimes)
    interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),
                         :L2=>sqrt(recursive_mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
    errors = merge(errors,interp_errors)
  else
    timeseries_analytic = sol2.u
    if sol.t == sol2.t
      errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol-timeseries_analytic))
      errors[:l2] = sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol-timeseries_analytic)))
    end
  end
  DiffEqBase.build_solution(sol,timeseries_analytic,errors)
end

"""
`appxtrue(sol::AbstractODESolution,sol2::AbstractODESolution)`

Uses the interpolant from the higher order solution sol2 to approximate
errors for sol. If sol2 has no interpolant, only the final error is
calculated.
"""
function appxtrue(sol::AbstractODESolution,sol2::AbstractODESolution;timeseries_errors=sol2.dense,dense_errors=sol2.dense)
  errors = Dict(:final=>recursive_mean(abs.(sol[end]-sol2[end])))
  if sol2.dense
    timeseries_analytic = sol2(sol.t)
    errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol-timeseries_analytic))
    errors[:l2] = sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol-timeseries_analytic)))
    if dense_errors
      densetimes = collect(range(sol.t[1],stop=sol.t[end],length=100))
      interp_u = sol(densetimes)
      interp_analytic = sol2(densetimes)
      interp_errors = Dict(:L∞=>maximum(vecvecapply((x)->abs.(x),interp_u-interp_analytic)),
                           :L2=>sqrt(recursive_mean(vecvecapply((x)->float(x).^2,interp_u-interp_analytic))))
      errors = merge(errors,interp_errors)
    end
  else
    timeseries_analytic = sol2.u
    if timeseries_errors && sol.t == sol2.t
      errors[:l∞] = maximum(vecvecapply((x)->abs.(x),sol-timeseries_analytic))
      errors[:l2] = sqrt(recursive_mean(vecvecapply((x)->float(x).^2,sol-timeseries_analytic)))
    end
  end
  DiffEqBase.build_solution(sol,timeseries_analytic,errors)
end

function appxtrue(sim::EnsembleSolution,appx_setup;kwargs...)
  _new_sols = Vector{DESolution}(length(sim.u))
  for i in eachindex(sim)
    prob = sim[i].prob
    prob2 = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,noise=NoiseWrapper(sim[i].W))
    true_sol = solve(prob2,appx_setup[:alg];appx_setup...)
    _new_sols[i] = appxtrue(sim[i],true_sol)
  end
  new_sols = convert(Vector{typeof(_new_sols[1])},_new_sols)
  calculate_ensemble_errors(new_sols;converged=sim.converged,elapsedTime=sim.elapsedTime,kwargs...)
end


text:
#%%
using Revise
using Glimpse
const Gl = Glimpse
using Glimpse: Spatial, Line, BufferColor, LineGeometry, VectorGeometry
using NearestNeighbors
using StaticArrays
@component @with_kw struct WallPlane <: ComponentData
    w1::Vec3{Float64} = Vec3(1.0, 0.0, 0.0)
    w2::Vec3{Float64} = Vec3(0.0, 1.0, 0.0)
    normal::Vec3{Float64} = normalize(cross(w1, w2))
end
function cube_planes(dio, origin, right)
    r = right - origin
    widths = (Vec3(abs(r[1]), 0.0, 0.0), Vec3(0.0, abs(r[2]), 0.0),
              Vec3(0.0, 0.0, abs(r[3])))
    heights = (Vec3(0.0, abs(r[2]), 0.0), Vec3(0.0, 0.0, abs(r[3])),
               Vec3(abs(r[1]), 0.0, 0.0))

    for (w, h) in zip(widths, heights)
        Entity(dio, WallPlane(; w1 = w, w2 = h), Spatial(; position = origin))
        Entity(dio, WallPlane(; w1 = -w, w2 = -h), Spatial(; position = right))
    end
    return Entity(dio,
                  Gl.assemble_wire_box(; left = Vec3f0(origin), right = Vec3f0(right))...)
end

struct WallBouncer <: System end
Overseer.requested_components(::WallBouncer) = (Spatial, WallPlane)

function Gl.update(::WallBouncer, m::Diorama)
    spat = m[Spatial]
    wp   = m[WallPlane]
    dt   = m[Gl.TimingData][1].dtime
    for we in @entities_in(spat && wp)
        n = wp[we].normal
        w1 = wp[we].w1
        w2 = wp[we].w2
        origin = spat[we].position
        for e in @entities_in(m[Spatial] && !wp)
            Threads.@spawn begin
                e_spat = spat[e]
                if !isapprox(dot(n, e_spat.velocity), 0)
                    p0 = e_spat.position
                    pr = dt * e_spat.velocity
                    t = dot(-p0 + origin, n) / dot(pr, n)
                    if 0 <= t <= 2.0
                        int_p = p0 + t * pr
                        if 0 <= dot(int_p - origin, normalize(w1)) <= norm(w1) &&
                           0 <= dot(int_p - origin, normalize(w2)) <= norm(w2)
                            spat[e] = Spatial(e_spat.position,
                                              e_spat.velocity -
                                              2n * dot(n, e_spat.velocity))
                        end
                    end
                end
            end
        end
    end
    return push!(m[Gl.UpdatedComponents][1], Spatial)
end

struct VelocityDrawer <: System end
Overseer.requested_components(::VelocityDrawer) = (Spatial, Gl.UniformColor)

function Gl.update(::VelocityDrawer, m::Diorama)
    spat = m[Spatial]
    geom = m[VectorGeometry]
    ucolor = m[Gl.UniformColor]
    r = [0; range(0, 0.7; length = 3)]
    for e in @entities_in(spat && geom)
        e = Entity(it)
        if norm(s.velocity) != 0.0
            geom[e] = LineGeometry([Point3f0(s.velocity * i) for i in r])
        end
    end
end

@component struct Boid end

struct Boids <: System end
Overseer.requested_components(::Boids) = (Boid,)
#%%
function Gl.update(::Boids, m::Diorama)
    spat = m[Spatial]
    dt = m[Gl.TimingData][1].dtime
    geom = m[Gl.PolygonGeometry]
    boid = m[Boid]
    it = @entities_in(spat && boid)
    points = map(e -> spat[e].position, it)
    tree = KDTree(points; leafsize = 10)
    wp = m[WallPlane]
    for e in @entities_in(spat && boid)
        Threads.@spawn begin
            s1      = spat[e]
            prev_v  = norm(s1.velocity)
            added_v = zero(Vec3f0)
            avg_pos = zero(Point3f0)
            avg_v   = zero(Vec3f0)
            ids     = inrange(tree, s1.position, 7, false)
            tot     = length(ids)
            for id in ids
                s2 = spat[it[id]]
                r = s1.position - s2.position
                if norm(r) < 1
                    added_v += r
                end
                avg_pos += s2.position
                avg_v   += s2.velocity
            end
            if tot != 0
                added_v += (avg_pos / tot - s1.position) / 200 +
                           (avg_v / tot - s1.velocity) / 20
            end
            spat[e] = Spatial(s1.position, normalize(s1.velocity + added_v) * prev_v)
        end
    end
end
dio = Gl.Diorama(Stage(:boids, [Boids(), WallBouncer(), VelocityDrawer()]);
                 background = RGBAf0(0.0, 0.0, 0.0, 1.0));
cube_planes(dio, Vec3(-60.0), Vec3(60.0))
expose(dio);
for i in -1000:1000
    t = Entity(dio,
               Gl.assemble_sphere(Point3f0(10 * (0.5 - rand(Point3f0)));
                                  velocity = 8 * normalize(0.5f0 - rand(Vec3f0)),
                                  radius = 0.5f0)..., Gl.LineOptions(),
               Gl.Spring(; k = 0.00001), Boid())
end
empty!(dio)
#%%
Gl.glfw_destroy_current_context()

using BenchmarkTools
@btime Gl.update(Boids(), $(dio.manager))
using StaticArrays
Point3f0 <: SVector
supertype(Point3f0)


text:
#####
##### Statistics helper functions
#####

"""
Returns the log probability of the data under a uniform multinomial distribution.
"""
function ll_uniform_multinomial(data)
    k = length(data)
    n = sum(data)
    n == 0 && return 0.
    logpdf(Multinomial(n, ones(k)./k), data)
end
"""
Returns the maximum likelihood estimate of the Multinomial rates and the log
probability of the counts under the Multinomial with these rates
"""
function MLE(counts)
    Σ = sum(counts)
    Σ == 0 && return zeros(length(counts)), 0.
    mle = counts ./ Σ
    mle, logpdf(Multinomial(Σ, mle), counts)
end
"""
Returns the normalization constant of a Dirichlet distribution with parameters
counts .+ 1
"""
function logmultinomialnormalization(counts)
    N = sum(counts)
    loggamma(N + 1) - +(map(x -> loggamma(x + 1), counts)...)
end
"""
Returns the log probability of the data under a Dirichlet Multinomial with
parameters `alpha` (`= ones(length(counts))` by default).
"""
function logpdf_dirichletmultinomial(counts, alpha = ones(length(counts)))
    n = sum(counts)
    n == 0 && return 0.
    logpdf(DirichletMultinomial(n, alpha), counts)
end

function _f(y, j)
    j == 1 && return prod(y)
    d = length(y)
    j == d + 1 && return 1 - y[1]
    prod(y[i] for i in 1:d-j+1) * (1 - y[d-j+2])
end
function get_constraints(args...)
    y -> begin
        (&)([a[2](_f(y, a[1]), _f(y, a[3])) for a in args]...)
    end
end
function _integrand(x, constraints)
    n = cumsum(x)
    K = length(n)
    y -> constraints(y) ? prod(y[i]^(n[K - i] - 1)*(1 - y[i])^(x[K - i + 1] - 1) for i in eachindex(y)) : 0.
end
function logpdf_dirichlet_int(x; constraints = y -> true)
    f = _integrand(x, constraints)
    K = length(x)
    res = hcubature(f, SVector{K-1}(zeros(K-1)), SVector{K-1}(ones(K-1)), maxevals = 10^5)
    res[2]/res[1] > .05 && @warn "Relative error in numerical integration = $(res[2]/res[1])"
    log(res[1])
end
logpdf_dirichletmultinomial_int(counts, c; alpha = ones(length(counts))) =
    logpdf_dirichletmultinomial_int(counts, get_constraints(c...), alpha = alpha)
function logpdf_dirichletmultinomial_int(counts, constraints::Function = y -> true;
                                         alpha = ones(length(counts)))
    x = counts .+ alpha
    logpdf_dirichlet_int(x, constraints = constraints) + loggamma(sum(counts) + 1) - sum(loggamma.(counts .+ 1)) - logpdf_dirichlet_int(alpha, constraints = constraints)
end

"""
Akaike information criterion for parameters `mle` with log likelihood `ll`
"""
AIC(n, ll) = 2n - 2ll

#####
##### helper functions
#####

"""
    partialsumsperbird(data, compartement_indep, food_indep, N)

Sums those counts that are assumed to have the same rate and returns additionally
the log probability that the summed counts are produced by a uniform Multinomial,
(see Methods, Eq. 2).

## Example
In the following we compute the sums and the log probability of uniform
distribution under the assumption that 30 decisions are taken in total and the
rates are independent of the compartment, i.e. they are the same for each
compartment, but they are different for each type of food. Data for different
compartments is in the rows, for different foodtypes in the columns.

```julia-repl
julia> partialsumsperbird([1 2;
                    3 4;
                    5 6], true, false, 30)
([9, 9, 12], -7.311519521215035)
```
"""
function partialsumsperbird(data, compartement_indep, food_indep, N)
    if compartement_indep && food_indep
        v = [sum(data)]
        ll = ll_uniform_multinomial(data[:])
    elseif compartement_indep
        v = sum(data, dims = 1)[:]
        ll = +([ll_uniform_multinomial(data[:, i]) for i in 1:length(v)]...)
    elseif food_indep
        v = sum(data, dims = 2)[:]
        ll = +([ll_uniform_multinomial(data[i, :]) for i in 1:length(v)]...)
    else
        v = data[:]
        ll = 0
    end
    if N > 0
        n = N - sum(data)
        @assert n >= 0
        push!(v, n), ll
    else
        v, ll
    end
end

"""
    partialsums(data, bird_indep, compartement_indep, food_indep, N)

Performs the partial sums per bird and, if the rates are assumed to be
independent of the bird identity, performs the partial sums over birds and adds
the correction factors (see Methods, Eq. 3).
"""
function partialsums(data, bird_indep, compartement_indep, food_indep, N)
    countsll = [partialsumsperbird(data[b], compartement_indep, food_indep, N)
                for b in 1:length(data)]
    counts = [c[1] for c in countsll]
    ll = +([c[2] for c in countsll]...)
    if bird_indep
        sum(counts),
        ll + +(map(logmultinomialnormalization, counts)...) -
            logmultinomialnormalization(sum(counts))
    else
        counts, ll
    end
end

"""
    logcategoricalposterior(data, bird_indep, compartement_indep, food_indep, N = 40)

Computes the log probability of the posterior given assumptions about the
independence of the different parameters. The strategy is to sum all counts that
are assumed to have the same rate and add the log probability that the individual
counts are generated with a uniform Multinomial.
"""
function logcategoricalposterior(data; bird_indep, compartement_indep, food_indep, N = 40, group = nothing, constraints = nothing, hyp = constraints)
    counts, ll = partialsums(data, bird_indep, compartement_indep, food_indep, N)
    if hyp == :cch
        if group === nothing
            constraints = ((1, <, 2), (3, <, 2), (4, >, 5), (6, >, 5))
        elseif group == :FN
            constraints = ((1, <, 2), (2, >, 3))
        elseif group == :NF
            constraints = ((1, >, 2), (2, <, 3))
        end
    elseif hyp == :fph2
        if group === nothing
            constraints = ((1, >, 2), (3, >, 2), (4, <, 5), (6, <, 5))
        elseif group == :FN
            constraints = ((1, >, 2), (2, <, 3))
        elseif group == :NF
            constraints = ((2, >, 1), (2, >, 3))
        end
    elseif hyp == :fph1
        if group === nothing
            constraints = ((1, >, 4), (1, >, 2), (1, >, 5), (1, >, 3), (1, >, 6))
        elseif group == :FN
            constraints = ((1, >, 2), (1, >, 3))
        elseif group == :NF
            constraints = ((2, >, 1), (2, >, 3))
        end
    elseif constraints !== nothing
        constraints = constraints
    else
        constraints = y -> true
    end
    if bird_indep
        if hyp === nothing
            logpdf_dirichletmultinomial(counts) + ll
        else
            logpdf_dirichletmultinomial_int(counts, constraints) + ll
        end
    else
        if hyp === nothing
            +(map(logpdf_dirichletmultinomial, counts)...) + ll
        else
            +(map(c -> logpdf_dirichletmultinomial_int(c, constraints), counts)...) + ll
        end
    end
end

"""
    compute_posterior(data, models; choices = [0], group = nothing)

This is the main entry point to compute the posteriors.
"""
function compute_posterior(data, models; choices = [0], group = nothing, unnormalizedlog = false)
    logposteriors = zeros(length(models))
    for k in choices
        logposteriors .+= [logcategoricalposterior(data; N = k,
                                                   group = group, c...) for c in models]
    end
    logposteriors ./= length(choices)                # uniform prior over choices
    unnormalizedlog && return logposteriors
    logposteriors .-= minimum(logposteriors)         # shift for numerical stability
    posteriors = exp.(logposteriors)/sum(exp.(logposteriors))  # exp and normalize
end

"""
    AIC(data, bird_indep, compartement_indep, food_indep, N = 40)

Computes the AIC given assumptions about the independence of the different
parameters.
"""
function AIC(data, bird_indep, compartement_indep, food_indep, N = 40)
    counts, ll = partialsums(data, bird_indep, compartement_indep, food_indep, N)
    if bird_indep
        mle, llmle = MLE(counts)
        AIC(length(mle), llmle + ll), mle
    else
        n = 0
        llmle = 0
        totalmle = []
        for c in counts
            mle, llmle_tmp = MLE(c)
            append!(totalmle, mle)
            n += length(mle)
            llmle += llmle_tmp
        end
        AIC(n, llmle + ll), totalmle
    end
end

"""
Returns an array where the numbers in row 1 correspond to the cached items in
the the first compartment etc. and numbers in column 1 corresponds to the
cached items of the first type the bird was prefed with etc.
"""
function sort_semantically(data)
    comp = copy(data.comp_order)
    prefed = copy(data.prefed_order)
    res = copy(data.cached)
    while comp[1] != :A
        comp = circshift(comp, 1)
        if size(res, 2) > 1
            res = circshift(res, (-1, 0))
        else
            res = circshift(res, -1)
        end
    end
    while prefed[1] != :P && size(res, 2) > 1
        prefed = circshift(prefed, 1)
        res = circshift(res, (0, 1))
    end
    (comp_order = comp, prefed_order = prefed, cached = res)
end

function load_data(file_name)
    data = YAML.load_file(file_name)
    result = Dict()
    for (key, content) in data
        cached = content["cached"]
        # Data is in row first form in the yaml file.
        length(cached) > 3 && (cached = Array(reshape(cached, :, 3)'))
        result[Symbol(key)] = (prefed_order = Symbol.(content["prefed_order"]),
                               comp_order = Symbol.(content["comp_order"]),
                               cached = cached)
    end
    data_sorted = Dict(name => sort_semantically(content) for (name, content) in result)
    (data = result,
     data_sorted = data_sorted,
     cached = [b.cached for b in values(data_sorted)])
end


text:
using Distributed
addprocs(2)

p1 = Vector{Any}(undef,3)
p2 = Vector{Any}(undef,3)
p3 = Vector{Any}(undef,3)

@everywhere begin
  using DiffEqMonteCarlo, StochasticDiffEq, DiffEqProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
  using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
  import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additive,
              prob_sde_linear, prob_sde_wave
end

using DiffEqMonteCarlo, StochasticDiffEq, DiffEqProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_additive,
            prob_sde_linear, prob_sde_wave

probs = Matrix{SDEProblem}(undef,3,3)
## Problem 1
prob = prob_sde_linear
probs[1,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
probs[1,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
probs[1,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
## Problem 2
prob = prob_sde_wave
probs[2,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
probs[2,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
probs[2,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
## Problem 3
prob = prob_sde_additive
probs[3,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
probs[3,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
probs[3,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))

fullMeans = Vector{Array}(undef,3)
fullMedians = Vector{Array}(undef,3)
fullElapsed = Vector{Array}(undef,3)
fullTols = Vector{Array}(undef,3)
offset = 0

Ns = [17,23,
17]

for k in 1:size(probs,1)
  global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
  println("Problem $k")
  ## Setup
  N = Ns[k]

  msims = Vector{Any}(undef,N)
  elapsed = Array{Float64}(undef,N,3)
  medians = Array{Float64}(undef,N,3)
  means   = Array{Float64}(undef,N,3)
  tols    = Array{Float64}(undef,N,3)

  #Compile
  prob = probs[k,1]
  ParallelDataTransfer.sendto(workers(), prob=prob)
  monte_prob = MonteCarloProblem(prob)
  solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,num_monte=1000,abstol=2.0^(-1),reltol=0)

  println("RSwM1")
  for i=1+offset:N+offset
    tols[i-offset,1] = 2.0^(-i-1)
    msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
                                            num_monte=1000,abstol=2.0^(-i-1),
                                            reltol=0,force_dtmin=true))
    elapsed[i-offset,1] = msims[i-offset].elapsedTime
    medians[i-offset,1] = msims[i-offset].error_medians[:final]
    means[i-offset,1]   = msims[i-offset].error_means[:final]
  end

  println("RSwM2")
  prob = probs[k,2]

  ParallelDataTransfer.sendto(workers(), prob=prob)
  monte_prob = MonteCarloProblem(prob)
  solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,num_monte=1000,abstol=2.0^(-1),reltol=0)

  for i=1+offset:N+offset
    tols[i-offset,2] = 2.0^(-i-1)
    msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
                                            num_monte=1000,abstol=2.0^(-i-1),
                                            reltol=0,force_dtmin=true))
    elapsed[i-offset,2] = msims[i-offset].elapsedTime
    medians[i-offset,2] = msims[i-offset].error_medians[:final]
    means[i-offset,2]   = msims[i-offset].error_means[:final]
  end

  println("RSwM3")
  prob = probs[k,3]
  ParallelDataTransfer.sendto(workers(), prob=prob)
  monte_prob = MonteCarloProblem(prob)
  solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,num_monte=1000,abstol=2.0^(-1),reltol=0)

  for i=1+offset:N+offset
    tols[i-offset,3] = 2.0^(-i-1)
        msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
                                    adaptive=true,num_monte=1000,abstol=2.0^(-i-1),
                                    reltol=0,force_dtmin=true))
    elapsed[i-offset,3] = msims[i-offset].elapsedTime
    medians[i-offset,3] = msims[i-offset].error_medians[:final]
    means[i-offset,3]   = msims[i-offset].error_means[:final]
  end

  fullMeans[k] = means
  fullMedians[k] =medians
  fullElapsed[k] = elapsed
  fullTols[k] = tols
end

gr(fmt=:svg)
lw=3
leg=String["RSwM1","RSwM2","RSwM3"]

titleFontSize = 16
guideFontSize = 14
legendFontSize= 14
tickFontSize  = 12

for k in 1:size(probs,1)
  global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
  p1[k] = Plots.plot(fullTols[k],fullMeans[k],xscale=:log10,yscale=:log10,  xguide="Absolute Tolerance",yguide="Mean Final Error",title="Example $k"  ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
  p2[k] = Plots.plot(fullTols[k],fullMedians[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Median Final Error",title="Example $k",linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
  p3[k] = Plots.plot(fullTols[k],fullElapsed[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Elapsed Time",title="Example $k"      ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
end

Plots.plot!(p1[1])
Plots.plot(p1[1],p1[2],p1[3],layout=(3,1),size=(1000,800))

#savefig("meanvstol.png")
#savefig("meanvstol.pdf")

plot(p3[1],p3[2],p3[3],layout=(3,1),size=(1000,800))
#savefig("timevstol.png")
#savefig("timevstol.pdf")

plot(p1[1],p3[1],p1[2],p3[2],p1[3],p3[3],layout=(3,2),size=(1000,800))

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


text:
using Pkg.TOML

@testset "TOML.jl" begin
	params = TOML.parsefile("TOML.toml")

	# test getTable function
	@test MPISimulations.getTable(params,"Test.Table") == Dict{AbstractString,Any}()
	@test_throws MissingException MPISimulations.getTable(params,"missingTable")
	@test ismissing(MPISimulations.getTable(params,"missingTable",tableOptional=true))

	#test getValue function
	@test MPISimulations.getValue(params,"Test","key") == 42
	@test_throws MissingException  MPISimulations.getValue(params,"Test","missingKey")
	@test_throws MissingException  MPISimulations.getValue(params,"missingTable","missingKey")
	@test ismissing(MPISimulations.getValue(params,"Test","missingKey",valueOptional=true))
	@test ismissing(MPISimulations.getValue(params,"missingTable","missingKey",valueOptional=true))

	# test getType function
	@test MPISimulations.getType(params,"Test") == ComplexF64
	@test ismissing(MPISimulations.getType(params,"Test.Table"))

	# test initialize function
	@test MPISimulations.initialize(Complex,"ComplexNumber",params,Float64) == ComplexF64(1.0,0.5)
	@test MPISimulations.initialize(Complex,"ComplexNumber",params,Float32) == ComplexF32(1.0,0.5)

	# test getTableName function
	@test MPISimulations.getTableName(Complex) == "Number.Complex"
end


text:
function save_result(res_dir::String, data;
        overwrite::Bool = false,
        comment = "",
        sources::Vector{String} = Vector{String}(),
        verbose = true)

    # Tree
    tree = Core.ResDirTree(res_dir);
    if !overwrite && isdir(tree)
        error("$(Core.resdir(tree)) already exist, to overwrite you must explicitly use the keyword!!!");
    end

    try

        # Dirs
        if !isdir(tree)
            mkpath(Core.resdir(tree));
            verbose && info("$(Core.resdir(tree)) created!!!")
        else
            verbose && info("$(Core.resdir(tree)) overwritten!!!")
        end
        mkpath(Core.sourcedir(tree));

        # Description
        desc = Core.build_description(res_dir, comment, sources);
        Core.save_desc(tree, desc);
        Core.add_to_log(tree, desc);

        # sources
        for src in sources
            if isdir(src)
                MyTools.copy_folder_gitless(src, Core.sourcedir(tree),
                    follow_symlinks = true, overwrite = overwrite);
                if verbose
                    info("copying: $(src) to $(Core.sourcedir(tree))");
                end
            elseif isfile(src)
                dest = joinpath(Core.sourcedir(tree), basename(src));
                cp(src, dest, follow_symlinks = true,
                    remove_destination = overwrite)
                verbose && info("copying: $(src) to $(Core.sourcedir(tree))");

            else
                ms = "Warn: copying: $(src) to $(Core.sourcedir(tree)) fails!!!"
                verbose && warn(ms);
                Core.add_to_log(tree, ms);
            end
        end

        # Data
        Core.save_data(tree, data);
    catch err
        if !overwrite
            rm(Core.resdir(tree), force = true, recursive = true);
        end
        throw(err);
    end

end

max_stars_repo_name: josePereiro/MyResultsManager

max_issues_repo_name: josePereiro/MyResultsManager

max_forks_repo_name: josePereiro/MyResultsManager


text:
# find midpoint of a range without causing an overflow to Inf
function robustmidpoint(a::Float64,b::Float64)
	l = min(a,b)
	u = max(a,b)
	if sign(l) == sign(u)
		m = l + (u-l)/2
	else
		m = (u+l)/2
	end
	m
end

# used mainly to adjust parameters of Geometric distribution away from boundary values
function adjusttoopeninterval(p, lower, upper)
	if p == lower
		p += 0.00001
	end
	if p == upper
		p -= 0.00001
	end
	p
end

#
# dummy DerivationState, Choice Context etc. to facilitate testing of internal choice model functions
#
import DataGenerators: getcurrentrecursiondepth, getimmediateancestorcpseqnumber, godelnumber

# for all apart from the ConditionalSampler, the choice context is not required
type DummyDerivationState <: DataGenerators.DerivationState
	currentrecursiondepth::Int
	function DummyDerivationState()
		new(1)
	end
end

# required for sampler choice model
getcurrentrecursiondepth(st::DummyDerivationState) = st.currentrecursiondepth
getimmediateancestorcpseqnumber(st::DummyDerivationState) = 0

dummyChoiceContext() = DataGenerators.ChoiceContext(DummyDerivationState(), :rule, 0, Int, 0, 0)

# choice model that only returns the lower bound Godel number (used to test generator)
type MinimumValueChoiceModel <: DataGenerators.ChoiceModel; end
function godelnumber(cm::MinimumValueChoiceModel, cc::DataGenerators.ChoiceContext)
	return cc.lowerbound, Dict()
end

function setminimumvaluechoicemodel!(g::DataGenerators.Generator)
	DataGenerators.setchoicemodel!(g, MinimumValueChoiceModel())
end