text:
##########################################
# Some uncertain datasets
##########################################
UV = UncertainValueDataset(example_uvals)

##########################################
# Apply functions to datasets `n` times
##########################################
n = 3

@test resample(median, UV, n) isa Vector{T} where T <: Real
@test resample(cor, UV, UV, n) isa Vector{T} where T <: Real


text:
"""
    grammify(x)

Turns a univariate 1-D array into a Grammian Angular Field sporting a 2-D matrix shape.

# Examples
```jldoctest
julia> grammify(1:4)
4×4 Array{Float64,2}:
  1.0        0.333333  -0.333333  -1.0
  0.333333  -0.777778  -1.0       -0.333333
 -0.333333  -1.0       -0.777778   0.333333
 -1.0       -0.333333   0.333333   1.0
```
"""
function grammify(x::AbstractVector)
    θ = acos.(standardize(x, -1, 1))
    # θ*θ'
    [cos(θ[i]+θ[j]) for i in 1:length(θ), j in 1:length(θ)]
end

"""
    plotgrammian(g)

Plot a grammified univariate array.

# Examples
```jldoctest
julia> GrammianAngularFields.plotgrammian(grammify(1:10))
data: [
  "heatmap with fields type and z"
]

layout: "layout with field margin"
```
"""
plotgrammian(g) = PlotlyJS.plot(PlotlyJS.heatmap(z=g))

"""
    plotgrammian2(x)

# Examples
```jldoctest
julia> GrammianAngularFields.plotgrammian2(sin.(1:30))
data: [
  "scatter with fields type, x, xaxis, y, and yaxis",
  "heatmap with fields type, xaxis, yaxis, and z"
]

layout: "layout with fields margin, xaxis1, xaxis2, yaxis1, and yaxis2"
```
"""
function plotgrammian2(x::AbstractVector)
    p1 = PlotlyJS.plot(PlotlyJS.scatter(x=1:length(x), y=x))
    p2 = PlotlyJS.plot(PlotlyJS.heatmap(z=grammify(x)))
    p = [p1 p2]
    p
end

"""
    testme(x)

Test any kind of functionality which may change at any time.
This is indeed the very definition of a random function.

# Examples
```jldoctest
julia> foo = 42
42
```
"""
function testme(x::AbstractVector)
    p1 = PlotlyJS.plot(PlotlyJS.scatter(x=1:length(x), y=x))
    p2 = PlotlyJS.plot(PlotlyJS.heatmap(z=grammify(x)))
    p = [p1 p2]
    p
end

ext: jl


text:
#========================================================================================#
"""
	Casinos

Module Casinos: A repository for fast access to arrays random numbers.

Author: Niall Palfreyman, 9/02/22
"""
module Casinos

# Externally callable methods of Casinos
export Casino, draw, shuffle!

using Random

#-----------------------------------------------------------------------------------------
# Module types:

"""
	Casino

A Casino can return arrays of random numbers in the range [0,1), up to a maximum number of rows
(maxrows), and a maximum number of columns (maxcols). It also contains a vault of prepared
random numbers from which it draws the arrays.
"""
struct Casino
	maxrows::Int					# Maximum number of drawable rows
	maxcols::Int					# Maximum number of drawable columns
	randomness::Int					# How randomised will our withdrawals be?
	vault::Matrix					# Repository of random numbers in [0,1)
	bernoulli						# Dictionary of Bernoulli outcomes

	"The one and only Casino constructor"
	function Casino(maxrows::Int,maxcols::Int,randomness::Int=5)
    	new(
			maxrows, maxcols, randomness,
			rand((maxrows+1)*randomness,(maxcols+1)*randomness),
			Dict{Float64,BitMatrix}()
		)
	end
end

#-----------------------------------------------------------------------------------------
# Module methods:

"""
	draw( casino, nrows, ncols)

Draw the required number of rows and columns from the casino vault, first ensuring that the
vault is large enough to support the withdrawal.
"""
function draw( casino::Casino, nrows::Int, ncols::Int)
	if nrows > casino.maxrows || ncols > casino.maxcols
		# Repository is too small - throw exception:
		error( "Requested withdrawal is too large")
	end

	# Access the vault:
	access( casino.vault, nrows, ncols)
end

#-----------------------------------------------------------------------------------------
"""
	draw( casino, nrows, ncols, bernoulli)

Draw (nrows x ncols) Boolean coin-toss values from the casino, using a bernoulli-biased coin.
"""
function draw( casino::Casino, nrows::Int, ncols::Int, bernoulli::Float64)
	if nrows > casino.maxrows || ncols > casino.maxcols
		# Vault is too small - throw exception:
		error( "Requested withdrawal is too large")
	end

	# Ensure bernoulli is a valid probability between 0 and 1:
	bernoulli = max( 0, min( 1, bernoulli))

	# Ensure Bernoulli entry exists in the Dictionary:
	if !haskey( casino.bernoulli, bernoulli)
		casino.bernoulli[bernoulli] = (casino.vault .< bernoulli)
	end

	# Access the bernoulli-biased matrix entry:
	access( casino.bernoulli[bernoulli], nrows, ncols)
end

#-----------------------------------------------------------------------------------------
"""
	shuffle!( casino)

Reassign random values in the vault.
"""
function shuffle!( casino::Casino)
    rand!( casino.vault)
end

#-----------------------------------------------------------------------------------------
"""
    access( matrix, nrows, ncols)

Access the required number of rows and columns from a matrix, on the
assumption that the matrix is big enough.
"""
function access( matrix, nrows::Int, ncols::Int)
	reprows, repcols = size(matrix)
	offset_r = rand( 1 : (reprows-nrows))
	stride_r = (nrows <= 1) ? 1 :
					rand( 1 : (reprows-offset_r) ÷ (nrows-1))
	offset_c = rand( 1 : (repcols-ncols))
	stride_c = (ncols <= 1) ? 1 :
					rand( 1 : (repcols-offset_c) ÷ (ncols-1))

	# Return a randomly chosen table of slices from the matrix:
	matrix[
		(offset_r : stride_r : (offset_r + (nrows-1)*stride_r)),
		(offset_c : stride_c : (offset_c + (ncols-1)*stride_c))
	]
end

#-----------------------------------------------------------------------------------------
"""
	unittest( casino)

Unit-test the Casinos module.
"""
function unittest()
	println("\n============ Unit test Casinos: ===============")
	println("Casino vault of randomness 2 for matrix withdrawals up to size (2x3):")
	casino = Casino(2,3,2)
	display( casino.vault)
	println()

	println("Draw several (2x3) matrices from the casino:")
	display( draw( casino,2,3)); println()
	display( draw( casino,2,3)); println()
	display( draw( casino,2,3)); println()

	println("Finally, reshuffle the casino and redisplay its vault:")
	shuffle!(casino)
	display(casino.vault)
end

end		# ... of module Casinos


text:
"""
This scheme adds support for multiple-objective based fitness approaches,
through the addition of domination calculation and pareto rank-based fitness.
"""
abstract ParetoFitnessScheme{T} <: FitnessScheme

type ParetoFitness{T}
  scores::Vector{T}
  rank::Integer

  ParetoFitness(s::Vector{T}) = new(s)
end

uses{T}(::ParetoFitnessScheme{T}) = ParetoFitness{T}
maximise(::ParetoFitnessScheme) = false
assign{T}(::ParetoFitnessScheme{T}, v::Vector{T}) = ParetoFitness{T}(v)

function compare(::ParetoFitnessScheme, x::ParetoFitness, y::ParetoFitness)
  if x.rank != y.rank
    return y.rank > x.rank ? -1 : 1
  end
  return 0
end

function dominates{T}(x::Vector{T}, y::Vector{T}, maximise::Vector{Bool})
  dom = false
  for (i, xi) in enumerate(x)
    if xi != y[i]
      if (xi > y[i]) == maximise[i]
        dom = true
      else
        return false
      end
    end
  end
  return dom
end


text:
module MaxwellGuidedMode

using Reexport
@reexport using MaxwellBase
using AbbreviatedTypes
using SimpleConstants
using VoxelwiseConstantMapping
using Arpack
using Dierckx
using LinearAlgebra: normalize!

include("model/model.jl")
include("scan_param/scan_param.jl")

end


text:
using Formatting

function prob(n, l)
    x = BigFloat(l * log(10))
    a = 1
    s = 1

    for k in 1:n-1
        a *= x/k
        s += a
    end

    s /= exp(x)
    return s
end

function main(n)
    l = Int(n * 0.4)
    h = Int(n * 0.5)
    e = 0.01

    while h - l > e
        if prob(n, (h + l)/2) > 0.25
            l = (h + l)/2
        else
            h = (h + l)/2
        end
    end

    println(format((h + l)/2, precision = 2))
end

main(10000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
# iterable which contains the options associated with Deflated Continuation
@with_kw struct DefContIterable{Tit, TperturbSolution, TacceptSolution, TupdateDeflationOp}
	it::Tit									# replicate PALC iterator
	maxBranches::Int64						# maximum number of (active) branches to be computed
	seekEveryStep::Int64 = 1				# whether to seek new (deflated) solution at every step
	perturbSolution::TperturbSolution 		# perturb function
	acceptSolution::TacceptSolution	  		# accept (solution) function
	updateDeflationOp::TupdateDeflationOp	# function to update the deflation operator
end

# state specific to Deflated Continuation, it is updated during the continuation process
mutable struct DCState{T, Tstate}
	tmp::T
	contState::Tstate
	isactive::Bool
	DCState(sol::T) where T = new{T, Nothing}(copy(sol), nothing, true)
	DCState(sol::T, state::ContState) where {T} = new{T, typeof(state)}(copy(sol), state, true)
end
# whether the branch is active
isActive(dc::DCState) = dc.isactive
# getters
getx(dc::DCState) = getx(dc.contState)
getp(dc::DCState) = getp(dc.contState)

function updatebranch!(iter::DefContIterable, dcstate::DCState, contResult::ContResult, defOp::DeflationOperator; current_param, step)
	isActive(dcstate) == false &&  return false, 0
	state = dcstate.contState 	# continuation state
	it = iter.it 				# continuation iterator
	@unpack step, ds, theta = state
	@unpack verbosity = it
	state.z_pred.p = current_param

	getPredictor!(state, it)
	sol1, fval, converged, itnewton = newton(it.F, it.J, getx(state), setParam(it,current_param), it.contParams.newtonOptions, defOp; normN = it.normC, callback = it.callbackN, iterationC = step, z0 = state.z_old)
	if converged
		# record previous parameter (cheap) and update current solution
		copyto!(state.z_old.u, sol1); state.z_old.p = current_param
		state.z_pred.p = current_param

		# Get tangent, it only mutates tau
		getTangent!(state.tau, state.z_pred, state.z_old, it, ds, theta, it.tangentAlgo, verbosity)

		# call user function to deal with DeflationOperator, allows to tackle symmetries
		iter.updateDeflationOp(defOp, sol1, current_param)

		# compute stability and bifurcation points
		computeEigenElements(it.contParams) && computeEigenvalues!(it, state)
		if it.contParams.detectBifurcation > 1 && detectBifucation(state)
			# we double-ckeck that the previous line, which mutated `state`, did not remove the bifurcation point
			if detectBifucation(state)
				_T  = eltype(it)
				_, bifpt = getBifurcationType(it.contParams, state, it.normC, it.printSolution, it.verbosity, :guess, getinterval(current_param, current_param-ds))
				if bifpt.type != :none; push!(contResult.specialpoint, bifpt); end
			end
		end
		state.step += 1
		save!(contResult, iter.it, state)
	else
		dcstate.isactive = false
		# save the last solution
		push!(contResult.sol, (x = copy(getx(state)), p = getp(state), step = state.step))
	end
	return converged, itnewton
end

"""
$(SIGNATURES)

This function computes the set of curves of solutions `γ(s) = (x(s), p(s))` to the equation `F(x,p)=0` based on the algorithm of **deflated continuation** as described in Farrell, Patrick E., Casper H. L. Beentjes, and Ásgeir Birkisson. “The Computation of Disconnected Bifurcation Diagrams.” ArXiv:1603.00809 [Math], March 2, 2016. http://arxiv.org/abs/1603.00809.

Depending on the options in `contParams`, it can locate the bifurcation points on each branch. Note that you can specify different predictors using `tangentAlgo`.

# Arguments:
- `F` is a function with input arguments `(x, p)`, where `p` is the set of parameters passed to `F`, and returning a vector `r` that represents the functional. For type stability, the types of `x` and `r` should match. In particular, it is not **inplace**,
- `J` is the jacobian of `F` at `(x, p)`. It can assume three forms.
    1. Either `J` is a function and `J(x, p)` returns a `::AbstractMatrix`. In this case, the default arguments of `contParams::ContinuationPar` will make `continuation` work.
    2. Or `J` is a function and `J(x, p)` returns a function taking one argument `dx` and returning `dr` of the same type as `dx`. In our notation, `dr = J * dx`. In this case, the default parameters of `contParams::ContinuationPar` will not work and you have to use a Matrix Free linear solver, for example `GMRESIterativeSolvers`,
    3. Or `J` is a function and `J(x, p)` returns a variable `j` which can assume any type. Then, you must implement a linear solver `ls` as a composite type, subtype of `AbstractLinearSolver` which is called like `ls(j, rhs)` and which returns the solution of the jacobian linear system. See for example `examples/SH2d-fronts-cuda.jl`. This linear solver is passed to `NewtonPar(linsolver = ls)` which itself passed to `ContinuationPar`. Similarly, you have to implement an eigensolver `eig` as a composite type, subtype of `AbstractEigenSolver`.
- `par` initial set of parameters,
- `lens::Lens` specifies which parameter axis among `par` is used for continuation. For example, if `par = (α = 1.0, β = 1)`, we can perform continuation w.r.t. `α` by using `lens = (@lens _.α)`. If you have an array `par = [ 1.0, 2.0]` and want to perform continuation w.r.t. the first variable, you can use `lens = (@lens _[1])`. For more information, we refer to `SetField.jl`,
- `contParams` parameters for continuation. See [`ContinuationPar`](@ref) for more information about the options,
- `defOp::DeflationOperator` a Deflation Operator (see [`DeflationOperator`](@ref)) which contains the set of solution guesses for the parameter `par`.

# Optional Arguments:
- `seekEveryStep::Int = 1` we look for additional solution, using deflated newton, every `seekEveryStep` step,
- `maxBranches::Int = 100` maximum number of branches considered,
- `maxIterDefOp::Int` maximum number of deflated Newton iterations
- `showplot = false` whether to plot the solution while computing,
- `printSolution = (x, p) -> norm(x)` function used to plot in the continuation curve. It is also used in the way results are saved. It could be `norm` or `(x, p) -> x[1]`. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...),
- `plotSolution = (x, p; kwargs...) -> nothing` function implementing the plot of the solution,
- `callbackN` callback for newton iterations. see docs for `newton`. Can be used to change preconditioners or affect the newton iterations. In the deflation part of the algorithm, when seeking for new branches, the callback is passed the keyword argument `fromDeflatedNewton = true` to tell the user can it is not in the continuation part (regular newton) of the algorithm,
- `tangentAlgo = NaturalPred()` controls the algorithm used to predict the tangents along the curve of solutions or the corrector. Can be `NaturalPred`, `SecantPred` or `BorderedPred`,
- `verbosity::Int` controls the amount of information printed during the continuation process. Must belong to `{0,⋯,5}`,
- `normN = norm` norm used in the different Newton solves,
- `dotPALC = (x, y) -> dot(x, y) / length(x)`, dot product used to define the weighted dot product (resp. norm) ``\\|(x, p)\\|^2_\\theta`` in the constraint ``N(x, p)`` (see below). This argument can be used to remove the factor `1/length(x)` for example in problems where the dimension of the state space changes (mesh adaptation, ...),
- `perturbSolution = (x, p, id) -> x` perturbation applied to the solution when trying to find new solutions using Deflated Newton. You can use for example `(x, p, id) -> x .+ (1 .+ 0.001 * rand(size(x)...))`

# Outputs:
- `contres::Vector{ContResult}` composite type which contains the computed branches. See [`ContResult`](@ref) for more information,
- the solutions at the last parameter value,
- current parameter value.
"""
function continuation(F, J, par, lens::Lens, contParams::ContinuationPar, defOp::DeflationOperator;
			verbosity::Int = 2,
			maxBranches::Int = 100,
			seekEveryStep::Int = 1,
			maxIterDefOp::Int = 5contParams.newtonOptions.maxIter,
			showplot::Bool = true,
			tangentAlgo = SecantPred(),
			linearAlgo = BorderingBLS(contParams.newtonOptions.linsolver),
			dotPALC = (x,y) -> dot(x,y) / length(x),
			printSolution = (x, p) -> norm(x),
			plotSolution = (x, p ;kwargs...) -> plot!(x; kwargs...),
			perturbSolution = (x, p, id) -> x,
			callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> true,
			acceptSolution = (x, p) -> true,
			updateDeflationOp = (defOp, x, p) -> push!(defOp, x),
			normN = norm) where vectype

	# allow to remove the corner case and associated specific return variables, type stable
	@assert length(defOp) > 0 "You must provide at least one guess"

	# function to get new solutions based on Deflated Newton
	function getNewSolution(_st::DCState, _p::Real, _idb)
		newton(F, J, perturbSolution(getx(_st), _p, _idb), set(par, lens, _p), setproperties(optnewton; maxIter = maxIterDefOp), deflationOp; normN = normN, callback = callbackN, fromDeflatedNewton = true)
	end

	# we make a copy of the deflation operator
	deflationOp = DeflationOperator(defOp.power, defOp.dot, defOp.α, deepcopy(defOp.roots))

	verbosity > 0 && printstyled(color=:magenta, "#"^51*"\n")
	verbosity > 0 && printstyled(color=:magenta, "--> There are $(length(deflationOp)) branches\n")

	# underlying continuation iterator
	contIt = ContIterable(F, J, defOp[1], par, lens, ContinuationPar(contParams, saveSolEveryStep = contParams.saveSolEveryStep == 0 ? Int(1e14) : contParams.saveSolEveryStep), linearAlgo; tangentAlgo = tangentAlgo, plot=showplot, plotSolution = plotSolution, printSolution = printSolution, normC = normN, dotPALC = dotPALC, finaliseSolution = finaliseDefault, callbackN = callbackN, verbosity = verbosity-2, filename = nothing)

	# we "hack" the saveSolEveryStep option because we always want to record the first point on each branch
	iter = DefContIterable(contIt, maxBranches, seekEveryStep, perturbSolution, acceptSolution, updateDeflationOp)

	# allocate states, one per branch
	states = [DCState(sol, iterate(@set contIt.x0 = sol)[1]) for sol in deflationOp.roots]

	# allocate branches to hold the result
	branches = [ContResult(contIt, state.contState) for state in states]

	S1 = empty(copy(deflationOp.roots))
	current_param = get(par, lens)

	# we extract the newton options
	optnewton = contParams.newtonOptions

	ii = 0
	while ((contParams.pMin < current_param < contParams.pMax) || ii == 0) &&
		 		ii < contParams.maxSteps
		# we update the parameter value
		current_param += contParams.ds
		current_param = clampPred(current_param, contIt)

		verbosity > 0 && println("#"^51)
		verbosity > 0 && println("--> step = $ii has $(length(branches)) branche(s), p = $current_param")

		# we empty the set of known solutions
		empty!(deflationOp.roots)

		# update the known branches
		for (idb, state) in enumerate(states)
			# this computes the solution for the new parameter value current_param
			# it also updates deflationOp
			flag, itnewton = updatebranch!(iter, state, branches[idb], deflationOp;
					current_param = current_param, step = ii)
			(verbosity >= 2 && isActive(state)) && println("----> Continuation for branch $idb in $itnewton iterations")
			verbosity >= 1 && ~flag && itnewton>0 && printstyled(color=:red, "--> Fold for branch $idb ?\n")
		end

		verbosity>1 && printstyled(color = :magenta,"--> looking for new branches\n")
		# number of branches
		nc = length(states)
		# number of active branches
		nactive = mapreduce(x -> x.isactive, +, states)
		if showplot && mod(ii, contParams.plotEveryStep) == 0
			plot(branches..., label = "", title  = "$nc branches, actives = $(nactive), step = $ii")
			for br in branches
				length(br) > 1 && plot!([br.branch[end-1:end].param], [getproperty(br.branch,1)[end-1:end]], label = "", arrow = true, color = :red)
			end
			scatter!([br.branch[1].param for br in branches], [br.branch[1][1] for br in branches], marker = :cross, color=:green, label = "") |> display
		end

		# only look for new branches if the number of active branches is too small
		if mod(ii, seekEveryStep) == 0
			if nactive < maxBranches
				n_active = 0
				# we restrict to 1:nc because we don't want to update the newly found branches
				for (idb, state) in enumerate(states[1:nc])
					# TODO this could be done in parallel
					if isActive(state) && (n_active < maxBranches)
						n_active += 1
						_success = true
						verbosity >= 2 && println("----> Deflation for branch $idb")
						while _success
							sol1, hist, _success, itnewton = getNewSolution(state, current_param, idb)
							if _success && normN(sol1 - getx(state)) < optnewton.tol
								@error "Same solution found for identical parameter value!!"
								_success = false
							end
							if _success
								verbosity>=1 && printstyled(color=:green, "--> new solution for branch $idb \n")
								push!(deflationOp.roots, sol1)
								push!(states, DCState(sol1, iterate(setproperties(contIt; x0 = sol1, par = set(par,lens,current_param)))[1]))
								push!(branches, ContResult(contIt, states[end].contState))
							end
						end
					end
				end
			end
		end

		ii += 1
	end
	showplot && display(plot(branches..., label = ""))
	return branches, contIt, [getx(c.contState) for c in states if isActive(c)], current_param
end


text:
#remarks: the J matrix is skewsymmetric in two cases:
#        1) both the energy and co-energy variables are discretized
#           at Gauss-Lagrange collocation points
#        2) at least the energy variables are discretized at
#          Gauss-Lagrange... why?!
function discrete_phs_closed(N,a,b)
	Nx = N   	#    energy variables: x
	Ne = Nx   # co-energy variables: e

	ei,we = lgwt(Ne,a,b)    # discretization of co-energy variables
	ei1 = [a; ei; b]
	#ei1 = collect(linspace(a,b,Ne,))
	#ei = ei1[2:end-1]
	if Nx > 1
		xi,wx = lgwt(Nx,a,b)
		#xi = collect(linspace(a,b,Nx))
	else
		xi,wx = (a+b)/2, (b-a);
	end
	#M = massmatrix(ei,xi,we);
	xquad, wquad = lgwt(Ne+2,a,b)
	M2 = massmatrix(ei1,xi, xquad, wquad);
	M2 = M2[:,2:end-1]
	M1 = massmatrix(ei,xi, xquad, wquad);
	D2 = dermatrix(ei,xi,1);
	D1 = dermatrix(ei1,xi,1);
	D1 = D1[:,2:end-1]

    phJ = -[zeros(N,N) D1 / (M2);
		          D2 / (M1) zeros(N,N)];
    phB = zeros(2*N,1);
    phD = [];

	if N > 1
		Q = massmatrix(xi,xi,xquad,wquad)
		#Q = sparse(diagm(wx[:],0));
		Q = blkdiag(Q,Q);
		else
		#Q = diagm([wx; wx][:]);
	end

	pflow = Collocation_points(xi, wx)
	peffort = Collocation_points(ei, we)

	p = Phs(phJ, phB, phD, Q)
	#p.Q = Q
	p.disc_data = Discrete_domain(pflow, peffort)
	#func(x :: Array) = (0.5*x'*p.Q*x)[1]
	#p.Hamiltonian = func
	#p.GradHam = ForwardDiff.gradient(p.Hamiltonian)
	p
end

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
export
		add_arc,
		add_node,
		get_all_nodes,
		get_children,
		get_first_node,
		get_next_node,
		get_node,
		get_number_of_nodes,
		get_parents,
		find_node,
		is_acyclic,
		read_file,
		write_file,
		set_default_BN_algorithm,
		set_default_ID_algorithm,
		set_default_HBN_algorithm,
		get_default_BN_algorithm,
		get_default_ID_algorithm,
		get_default_HBN_algorithm,
		DSL_ALG_BN_LAURITZEN, DSL_ALG_BN_HENRION, DSL_ALG_BN_PEARL, DSL_ALG_BN_LSAMPLING,
		DSL_ALG_BN_SELFIMPORTANCE, DSL_ALG_BN_HEURISTICIMPORTANCE, DSL_ALG_BN_BACKSAMPLING,
		DSL_ALG_BN_AISSAMPLING, DSL_ALG_BN_EPISSAMPLING, DSL_ALG_BN_LBP, DSL_ALG_BN_LAURITZEN_OLD,
		DSL_ALG_BN_RELEVANCEDECOMP, DSL_ALG_BN_RELEVANCEDECOMP2, DSL_ALG_HBN_HEPI,
		DSL_ALG_HBN_HLW, DSL_ALG_HBN_HLBP, DSL_ALG_HBN_HLOGICSAMPLING,
		update_beliefs,
		invalidate_all_beliefs,
		call_ID_algorithm,
		call_BN_algorithm,
		call_Eq_algorithm,
		call_HBN_algorithm,
		clear_all_evidence,
		clear_all_decision,
		clear_all_propagated_evidence,
		is_there_any_evidence,
		is_there_any_decision,
		partial_ordering

function add_arc( net::Network, theParent::Integer, theChild::Integer )
	err_int = ccall( (:network_AddArc, LIB_SMILE), Int32, (Ptr{Void}, Int32, Int32), net.ptr, theParent, theChild)
	if err_int < 0
		error( "network_AddArc: code: ", err_int )
	end
end

function add_node( net::Network, nodeType::Integer, nodeName::AbstractString )

	if !is_nodetype_id(nodeType)
		error( "network_AddNode: nodeType not a recognized type" )
	end

	error_buffer = Array(UInt8, 128)
	nodeId = ccall( (:network_AddNode, LIB_SMILE), Int32, (Ptr{Void}, Int32, Ptr{UInt8}, Ptr{UInt8}), net.ptr, nodeType, bytestring(nodeName), error_buffer )
	if nodeId < 0
		error_str = bytestring(error_buffer)
		error( "network_AddNode: ", error_str )
	end
	return nodeId
end

function get_all_nodes( net::Network )

	n = get_number_of_nodes(net)
	if n == 0
		return Int[]
	end

	nodeHandles = Array(Int,n)
	nodeHandles[1] = get_first_node( net )
	for i = 2 : n
		nodeHandles[i] = get_next_node( net, nodeHandles[i-1] )
	end

	return nodeHandles
end

function get_children( net::Network, ofThisNode::Integer, len_start = 10 )
	len = Array(UInt32,1)
	len[1] = len_start
	arr = Array(Int32, len_start)
	ccall( (:network_GetChildren, LIB_SMILE), Void, (Ptr{Void},Int32,Ptr{Int32},Ptr{UInt32}), net.ptr, ofThisNode, arr, len )
	len_end = len[1]

	if len_end == 0
		return []
	end
	if len_end > len_start
		error("get_children:: not enough spaces allocated in array")
	end
	return arr[1:len_end]
end

function get_first_node( net::Network )
	retval = ccall( (:network_GetFirstNode, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
	return retval
end

function get_next_node( net::Network, ofThisNode::Integer )
	retval = ccall( (:network_GetNextNode, LIB_SMILE), Int32, (Ptr{Void},Int32), net.ptr, ofThisNode )
	return retval
end

function get_node( net::Network, nodeHandle::Integer )
	ptr = ccall( (:network_GetNode, LIB_SMILE), Ptr{Void}, (Ptr{Void},Int32), net.ptr, nodeHandle )
	return Node(ptr)
end

function get_number_of_nodes( net::Network )
	retval = ccall( (:network_GetNumberOfNodes, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
	return retval
end

function get_parents( net::Network, ofThisNode::Integer, len_start = 10 )
	len = Array(UInt32,1)
	len[1] = len_start
	arr = Array(Int32, len_start)
	ccall( (:network_GetParents, LIB_SMILE), Void, (Ptr{Void},Int32,Ptr{Int32},Ptr{UInt32}), net.ptr, ofThisNode, arr, len )
	len_end = len[1]

	if len_end == 0
		return []
	end
	if len_end > len_start
		error("get_parents:: not enough spaces allocated in array")
	end
	return arr[1:len_end]
end

find_node( net::Network, withThisID::AbstractString ) =
	ccall( (:network_FindNode, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8}), net.ptr, bytestring(withThisID) )

function is_acyclic( net::Network )
	retval = ccall( (:network_IsAcyclic, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
	return retval != 0 # true if non zero
end

function read_file( net::Network, fileName::AbstractString )
	# TODO: check that fileName ends with a valid file type
	retval = ccall( (:network_ReadFile, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8}), net.ptr, bytestring(fileName) )
	return retval
end

function write_file( net::Network, fileName::AbstractString )
	ext = splitext(fileName)[2]
	@assert(ext == ".dsl" || ext == ".xdsl")

	retval = ccall( (:network_WriteFile, LIB_SMILE), Int32, (Ptr{Void},Ptr{UInt8}), net.ptr, bytestring(fileName) )
	return retval
end

# BN algorithms
const DSL_ALG_BN_LAURITZEN            = Int32( 0)
const DSL_ALG_BN_HENRION              = Int32( 1)
const DSL_ALG_BN_PEARL                = Int32( 2)
const DSL_ALG_BN_LSAMPLING            = Int32( 3)
const DSL_ALG_BN_SELFIMPORTANCE       = Int32( 4)
const DSL_ALG_BN_HEURISTICIMPORTANCE  = Int32( 5)
const DSL_ALG_BN_BACKSAMPLING         = Int32( 6)
const DSL_ALG_BN_AISSAMPLING          = Int32( 7)
const DSL_ALG_BN_EPISSAMPLING         = Int32( 8)
const DSL_ALG_BN_LBP				          = Int32( 9)
const DSL_ALG_BN_LAURITZEN_OLD        = Int32(10)
const DSL_ALG_BN_RELEVANCEDECOMP      = Int32(11)
const DSL_ALG_BN_RELEVANCEDECOMP2     = Int32(12)
const DSL_ALG_HBN_HEPIS				        = Int32(13)
const DSL_ALG_HBN_HLW                 = Int32(14)
const DSL_ALG_HBN_HLBP				        = Int32(15)
const DSL_ALG_HBN_HLOGICSAMPLING	    = Int32(16)

# ID algorithms
const DSL_ALG_ID_COOPERSOLVING        = Int32(0)
const DSL_ALG_ID_SHACHTER             = Int32(1)

set_default_BN_algorithm(net::Network,  algorithm::Integer) = ccall( (:network_SetDefaultBNAlgorithm,  LIB_SMILE), Void, (Ptr{Void},Int32), net.ptr, algorithm )
set_default_ID_algorithm(net::Network,  algorithm::Integer) = ccall( (:network_SetDefaultIDAlgorithm,  LIB_SMILE), Void, (Ptr{Void},Int32), net.ptr, algorithm )
set_default_HBN_algorithm(net::Network, algorithm::Integer) = ccall( (:network_SetDefaultHBNAlgorithm, LIB_SMILE), Void, (Ptr{Void},Int32), net.ptr, algorithm )
get_default_BN_algorithm(net::Network)  = ccall( (:network_GetDefaultBNAlgorithm,  LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
get_default_ID_algorithm(net::Network)  = ccall( (:network_GetDefaultIDAlgorithm,  LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
get_default_HBN_algorithm(net::Network) = ccall( (:network_GetDefaultHBNAlgorithm, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )

# Inference
update_beliefs(net::Network)         = ccall( (:network_UpdateBeliefs,        LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
invalidate_all_beliefs(net::Network) = ccall( (:network_InvalidateAllBeliefs, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
call_ID_algorithm(net::Network)      = ccall( (:network_CallIDAlgorithm,      LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
call_BN_algorithm(net::Network)      = ccall( (:network_CallBNAlgorithm,      LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
call_Eq_algorithm(net::Network)      = ccall( (:network_CallEqAlgorithm,      LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
call_HBN_algorithm(net::Network)     = ccall( (:network_CallHBNAlgorithm,     LIB_SMILE), Int32, (Ptr{Void},), net.ptr )

clear_all_evidence(net::Network)            = ccall( (:network_ClearAllEvidence, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
clear_all_decision(net::Network)            = ccall( (:network_ClearAllDecision, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
clear_all_propagated_evidence(net::Network) = ccall( (:network_ClearAllPropagatedEvidence, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
is_there_any_evidence(net::Network)         = ccall( (:network_IsThereAnyEvidence, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )
is_there_any_decision(net::Network)         = ccall( (:network_IsThereAnyDecision, LIB_SMILE), Int32, (Ptr{Void},), net.ptr )

partial_ordering(net::Network) = IntArray(ccall( (:network_PartialOrdering, LIB_SMILE), Ptr{Void}, (Ptr{Void},), net.ptr ))


text:
using InfrastructureModels

IM = InfrastructureModels

file = ARGS[1]

mp_data = open(file) do io
	data_string = read(io, String)
	# We need to manually patch `data_string` by inserting '%column_names%' where needed.
	# This is because `PowerModels.export_matpower()` doesn't write it,
	# whereas `InfrastructureModels.parse_matlab_string()` expects it...
	data_lines = split(data_string, '\n')
	last = length(data_lines)
	for (k, line) in enumerate(data_lines)
		# This below is just a dirty condition that works for the files in this project
		# I don't guarantee that it will work for *any* .m file
		if k > 1 && startswith(line, "%") && startswith(data_lines[k - 1], "%%")
			data_lines[k] = "%column_names%" * line[2:end]
		end
		if startswith(line, "%% convert branch impedances")
			last = k - 1
			break
		end
	end
	data_lines = data_lines[1:last]
	data_string = join(data_lines, '\n')
	matlab_data, func_name, colnames = IM.parse_matlab_string(data_string, extended=true)
	# write CSVs
	for key in ("bus", "gen", "branch")
		mpckey = "mpc.$key"
		outfile = "$key.csv"
		open(outfile, "w") do f
			column_names = colnames[mpckey]
			if key in ("gen", "branch")
				pushfirst!(column_names, "$(key)_name")
			end
			println(f, join(column_names, ", "))
			for (i, row) in enumerate(matlab_data[mpckey])
				if key in ("gen", "branch")
					row = string.(row)
					pushfirst!(row, "$(key)_$i")
				end
				str_row = join(row, ", ")
				println(f, str_row)
			end
		end
	end
end

max_stars_repo_name: Energy-Reform/Mahdavi-test-systems

max_issues_repo_name: Energy-Reform/Mahdavi-test-systems

max_forks_repo_name: Energy-Reform/Mahdavi-test-systems


text:
abstract type AbstractAbundanceTable{T <: Real,
                                     F <: AbstractFeature,
                                     S <: AbstractSample} <: EcoBase.AbstractAssemblage{T, F, S}
end

"""
    CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}

An `AbstractAssemblage` from [EcoBase.jl](https://github.com/EcoJulia/EcoBase.jl)
that uses an `AxisArray` of a `SparseMatrixCSC` under the hood.

`CommunityProfile`s are tables with `AbstractFeature`-indexed rows and
`AbstractSample`-indexed columns.
Note - we can use the `name` of samples and features to index.

```jldoctest community
julia> txs = [Taxon("taxon\$i") for i in 1:10];

julia> mss = [MicrobiomeSample("sample\$i") for i in 1:5];

julia> mat = spzeros(10,5);

julia> for i in 1:5; mat[i,i] = 1.; end

julia> comm = CommunityProfile(mat, txs, mss)
CommunityProfile{Float64, Taxon, MicrobiomeSample} with 10 things in 5 places

Thing names:
taxon1, taxon2, taxon3...taxon9, taxon10

Place names:
sample1, sample2, sample3, sample4, sample5

julia> comm["taxon1", "sample1"]
1.0

julia> comm[:,["sample1", "sample5"]]
CommunityProfile{Float64, Taxon, MicrobiomeSample} with 10 things in 2 places

Thing names:
taxon1, taxon2, taxon3...taxon9, taxon10

Place names:
sample1, sample5

julia> comm[Taxon("taxon3", :kingdom), "sample1"]
0.0
```
"""
mutable struct CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}
    aa::NamedAxisArray

    function CommunityProfile(aa::NamedAxisArray)
        @assert dimnames(aa) == (:features, :samples)
        T = eltype(parent(aa))
        F = eltype(keys(axes(aa, 1)))
        S = eltype(keys(axes(aa, 2)))
        return new{T, F, S}(aa)
    end
end

function CommunityProfile(tab::SparseMatrixCSC{<:Real},
                          features::AbstractVector{<:AbstractFeature},
                          samples::AbstractVector{<:AbstractSample})
    return CommunityProfile(NamedAxisArray(tab, features=features, samples=samples))
end

function CommunityProfile{T, F, S}(tab::SparseMatrixCSC{<:T},
                                   features::AbstractVector{F},
                                   samples::AbstractVector{S}) where {T, F, S}
    return CommunityProfile(tab, features, samples)
end
## -- Convienience functions -- ##

"""
    features(at::AbstractAbundanceTable)

Returns features in `at`. To get featurenames instead, use [`featurenames`](@ref).
"""
features(at::AbstractAbundanceTable) = axes(at.aa, 1) |> keys

"""
    samples(at::AbstractAbundanceTable)

Returns samples in `at`. To get samplenames instead, use [`samplenames`](@ref).
"""
samples(at::AbstractAbundanceTable) = axes(at.aa, 2) |> keys

profiletype(at::AbstractAbundanceTable) = eltype(features(at))
clades(at::AbstractAbundanceTable) = clade.(features(at))

Base.size(at::AbstractAbundanceTable, dims...) = size(at.aa, dims...)

# -- Indexing -- #

function Base.getindex(at::CommunityProfile, inds...)
    idx = at.aa[inds...]

    # single value - return that value
    ndims(idx) == 0 && return idx
    # another table - return a new CommunityProfile with that table
    ndims(idx) == 2 && return CommunityProfile(idx)
    # a row or a column, figure out which, and make it 2D
    if ndims(idx) == 1
        dn = dimnames(idx)[1]
        # if it's a row...
        if dn == :samples
            return at[[inds[1]], inds[2]]
        # if it's a column
        elseif dn == :features
            return at[inds[1], [inds[2]]]
        end
    end
end

## -- EcoBase Translations -- ##
# see src/ecobase.jl for Microbiome function names
# thing => feature
# place => sample
# occurrences => abundances (or totals)

EcoBase.thingnames(at::AbstractAbundanceTable) = name.(features(at))
EcoBase.placenames(at::AbstractAbundanceTable) = name.(samples(at))
EcoBase.occurrences(at::AbstractAbundanceTable) = parent(parent(at.aa)) # first parent is the unnamed AxisArray
EcoBase.nthings(at::AbstractAbundanceTable) = size(at, 1)
EcoBase.nplaces(at::AbstractAbundanceTable) = size(at, 2)
# # todo
# EcoBase.thingoccurrences(at::AbstractAbundanceTable, things) = nothing
# EcoBase.placeoccurrences(at::AbstractAbundanceTable, places) = nothing

"""
    featuretotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a nfeatures x 1 `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(featuretotals(at))`.
"""
featuretotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=2)

"""
    sampletotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a 1 x nsamples `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(sampletotals(at))`.
"""
sampletotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=1)

## -- Tables Interface -- ##

Tables.istable(::AbstractAbundanceTable) = true
Tables.columnaccess(::AbstractAbundanceTable) = true
Tables.rowaccess(::AbstractAbundanceTable) = true

Tables.getcolumn(at::AbstractAbundanceTable, i::Int) = i == 1 ? featurenames(at) : abundances(at[:, i-1])
Tables.getcolumn(at::AbstractAbundanceTable, i::AbstractString) = i == "features" ? featurenames(at) : abundances(at[:, i])
Tables.getcolumn(at::AbstractAbundanceTable, i::Symbol) = Tables.getcolumn(at, string(i))

Tables.columnnames(at::AbstractAbundanceTable) = [:features, Symbol.(samplenames(at))...]

function Tables.schema(at::AbstractAbundanceTable)
    elt = eltype(abundances(at))
    coltypes = [eltype(features(at)), (elt for _ in 1:nsamples(at))...]
    return Tables.Schema(Tables.columnnames(at), coltypes)
end

Tables.columns(at::AbstractAbundanceTable) = (; (col => Tables.getcolumn(at, col) for col in Tables.columnnames(at))...)

function _makerow(row::AbstractAbundanceTable)
    size(row, 1) == 1 || error("Can't make row from table of size $(size(row))")
    NamedTuple{(:features, Symbol.(samplenames(row))...)}((name(features(row)[1]), abundances(row)...))
end

Tables.rows(at::AbstractAbundanceTable) = (_makerow(at[i, :]) for i in 1:nfeatures(at))

## -- Methods for absolute and relative abundances -- ##

# """
#     filterabund(abun::AbstractComMatrix, n::Int=minimum(10, nfeatures(abun)))

# Filter an abundance table to the top `n` features accross all samples

# This function also adds a row for "other", which sums the abundances of the
# remaining features.
# """
# function filterabund(abun::AbstractAbundanceTable, n::Int=minimum(10, nfeatures(abun)))
#     # TODO: add prevalence filter

#     totals = featuretotals(abun)

#     srt = sortperm(totals, rev=true)

#     newabun = getfeature(abun, srt[1:n])

#     remainder = [sum(occurrences(abun)[srt[n+1:end], i]) for i in 1:size(abun, 2)]'
#     newabun = vcat(newabun, remainder)
#     newrows = cat(featurenames(abun)[srt[1:n]], ["other"], dims=1)

#     return abundancetable(newabun, samplenames(abun), newrows)
# end

"""
    relativeabundance!(a::AbstractAbundanceTable; kind::Symbol=:fraction)

Normalize each sample in AbstractAbundanceTable to the sum of the sample.

By default, columns sum to 1.0.
Use `kind=:percent` for columns to sum to 100.
"""
function relativeabundance!(at::AbstractAbundanceTable; kind::Symbol=:fraction)
    in(kind, [:percent, :fraction]) || throw(ArgumentError("Invalid kind: $kind"))
    eltype(abundances(at)) <: AbstractFloat || throw(ArgumentError("relativeabundance! requires profile to have AbstractFloat eltype. Try relativeabundance instead"))
    abund = abundances(at)
    abund ./= sampletotals(at)
    kind == :percent && (abund .*= 100)
    dropzeros!(abund) # shouldn't need dropzeros - https://github.com/JuliaLang/julia/issues/39018
    return at
end

"""
    relativeabundance(at::AbstractAbundanceTable, kind::Symbol=:fraction)

Like [`relativeabundance!`](@ref), but does not mutate original.
"""
function relativeabundance(at::T, kind::Symbol=:fraction) where T <: AbstractAbundanceTable
    comm = T(float.(abundances(at)), deepcopy(features(at)), deepcopy(samples(at)))
    relativeabundance!(comm)
end

"""
    present(t::Union{Real, Missing}, minabundance::Real=0.0)
    present(at::AbstractAbundanceTable, minabundance::Real=0.0)

Check if a given (non-zero) value is greater than or equal to a minimum value.
If the minimum abundance is 0, just checks if value is non-zero.

If used on an `AbstractAbundanceTable`, returns a sparse boolean matrix of the same size.
"""
function present(t::Real, minabundance::Real=0.0)
    (minabundance >= 0 && t >= 0) || throw(DomainError("Only defined for positive values"))
    t == 0 ? false : t >= minabundance
end

present(::Missing, m) = missing

function present(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mat = spzeros(Bool, Ssize(at)...)
    for i in eachindex(mat)
        mat[i] = present(at[i], minabundance)
    end
    return mat
end

"""
    prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0)
    prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)

Return the fraction of values that are greater than or equal to a minimum.
If the minimum abundance is 0, returns the fraction of non-zero values.

If used on an `AbstractAbundanceTable`,
returns a prevalence value for each `feature` accross the `sample`s.
"""
prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0) = mean(x-> present(x, minabundance), a)

# makes it work for any iterable
prevalence(a, minabundance::Real=0.0) = mean(x-> present(x, minabundance), (y for y in a))

function prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mean(x-> present(x, minabundance), abundances(at), dims=2)
end


text:
using BinDeps
using Compat

@BinDeps.setup

sqlite3 = library_dependency("sqlite3",
    aliases=["sqlite3", "libsqlite3", "sqlite3-64", "libsqlite3-0"])

@osx_only begin
  using Homebrew
  provides(Homebrew.HB, @compat Dict("sqlite3"=>sqlite3))
end

@unix_only begin
  provides(AptGet, @compat Dict("sqlite3"=>sqlite3))
  try
    py_info = readall(`pip show sqlalchemy`)
    m = match(r"Version: (.*)", py_info)
    version = @compat VersionNumber(m.captures[1])
    if version < v"1"
      info("Attempting to install sqlalchemy python package via pip")
      try
        run(`pip install --upgrade sqlalchemy`)
      catch
        run(`sudo pip install --upgrade sqlalchemy`)
      end
    end
  catch err
    warn("Couldn't automatically install sqlalchemy. Run 'pip install sqlalchemy' manually.\nError was $err")
  end
end

@windows_only begin
    using WinRPM
    provides(WinRPM.RPM, @compat Dict("libsqlite3-0"=>sqlite3))
end

@BinDeps.install


text:
## ----------  Generic

"""
	$(SIGNATURES)

Specifies dropout rules for a set colleges.
"""
abstract type DropoutRuleSet <: ModelObject end

"""
	$(SIGNATURES)

Switches used to construct a `DropoutRuleSet`.
"""
abstract type DropoutRuleSwitches end

"""
	$(SIGNATURES)

Dropout rule for one college. Constructed from the `DropoutRuleSet`.
"""
abstract type DropoutRule end

# Stub
validate_drop_set(dr :: DropoutRuleSet) = true;

Base.show(io :: IO,  ds :: DropoutRuleSet) =
    print(io, "DropoutRuleSet");

"""
	$(SIGNATURES)

Table with settings. Columns are explanation and values.
All formatted into strings.
"""
function settings_table(ds :: DropoutRuleSet)
    timeV = Int.(college_durations(ds));
    return [
        "DropoutRuleSet"  " ";
        "Must drop at end of"  "$timeV"
    ]
end

StructLH.describe(ds :: DropoutRuleSet) = settings_table(ds);

ModelParams.has_pvector(ds :: DropoutRuleSet) = false;
college_durations(d :: DropoutRuleSwitches) = d.tMaxV;
college_durations(d :: DropoutRuleSet) = college_durations(d.switches);
college_duration(d :: DropoutRuleSwitches, ic) = d.tMaxV[ic];
college_duration(ds :: DropoutRuleSet, ic) = college_duration(ds.switches, ic);
college_duration(d :: DropoutRule) = d.tMustDrop;

# -----------------  Simple DropoutRule
# Nobody is forced to drop out until last year in college.
# Except for a small fraction for numerical reasons.

"""
	$(SIGNATURES)

Switches for constructing a simple `DropoutRuleSet`.
"""
mutable struct DropoutRuleSwitchesSimple <: DropoutRuleSwitches
    "Students can stay at most this long in each college."
    tMaxV :: Vector{TimeInt}
    # Probability that student must drop out before terminal year
    # Small. There for numerical reasons.
    probMin :: Double
end

"""
	$(SIGNATURES)

Simple `DropoutRuleSet` where students must drop out after a fixed number of years in college.
"""
mutable struct DropoutRuleSetSimple <: DropoutRuleSet
    objId :: ObjectId
    switches :: DropoutRuleSwitchesSimple
end

"""
	$(SIGNATURES)

Simple `DroputRule`.
"""
struct DropoutRuleSimple <: DropoutRule
    # Must drop out (if not graduated) after this many years in college
    tMustDrop :: TimeInt
    probMin :: Double
end

## ------------  Init

"""
	$(SIGNATURES)

Construct a `DropoutRuleSet` from its switches.
"""
make_dropout_set(objId :: ObjectId,  switches :: DropoutRuleSwitchesSimple) =
    DropoutRuleSetSimple(objId, switches);

"""
	$(SIGNATURES)

Construct the `DropoutRule` for one college.
"""
function make_dropout_rule(ds :: DropoutRuleSetSimple, iCollege :: Integer)
    return DropoutRuleSimple(college_duration(ds, iCollege), ds.switches.probMin);
end

# --------------------  For one college

"""
    $(SIGNATURES)

Dropout probability at end of t.
Inputs are h, n at END of t.
"""
function drop_prob(dr :: DropoutRuleSimple,  t :: Integer,
    h :: Double,  n :: Integer)

    if t >= dr.tMustDrop
        dProb = 1.0;
    else
        dProb = dr.probMin;
    end
    return dProb
end

# Array inputs
function drop_prob(g :: DropoutRuleSimple, t :: Integer,
    hM :: Array{Double},  nM :: Array{T1}) where T1 <: Integer

    if t >= g.tMustDrop
        probM = ones(Double, size(hM));
    else
        sizeV = size(hM);
        @assert size(nM) == sizeV
        probM = Array{Double}(undef, sizeV);
        for i1 in eachindex(hM)
            probM[i1] = drop_prob(g, t, hM[i1], nM[i1]);
        end
    end
    return probM :: Array{Double}
end

"""
	$(SIGNATURES)

Dropout probability on an [h, n] grid; at end of t.
"""
function drop_prob_grid(g :: DropoutRuleSimple,  t :: Integer,  hV :: Vector{Double},
    nTakenV :: Vector{TI})  where  TI <: Integer

    nh = length(hV);
    nn = length(nTakenV);
    if t >= g.tMustDrop
        probM = ones(Double, nh, nn);
    else
        probM = Matrix{Double}(undef, nh, nn);
        for i_h = 1 : nh
            for i_n = 1 : nn
                probM[i_h, i_n] = drop_prob(g, t, hV[i_h], nTakenV[i_n]);
            end
        end
    end
    validate_drop_prob_grid(probM)
    return probM :: Matrix{Double}
end

function validate_drop_prob_grid(prob_hnM :: Matrix{Double})
    isValid = true;
    isValid = isValid && all_at_least(prob_hnM, 0.0) && all_at_most(prob_hnM, 1.0);
    isValid = isValid && all_at_most(diff(prob_hnM, dims = 2), 0.0);
    isValid = isValid && all_at_most(diff(prob_hnM, dims = 1), 0.0);
end

# ---------------------


text:
"""
    positive(val::Real, transform=exp, ε=sqrt(eps(typeof(val))))

Return a `Positive`.
The `value` of a `Positive` is a `Real` number that is constrained to be positive.
This is represented in terms of a `transform` that maps an `unconstrained_value` to the
positive reals.
Satisfies `val ≈ transform(unconstrained_value)`.
"""
function positive(val::Real, transform=exp, ε=sqrt(eps(typeof(val))))
    val > 0 || throw(ArgumentError("Value ($val) is not positive."))
    val > ε || throw(ArgumentError("Value ($val) is too small, relative to ε ($ε)."))
    unconstrained_value = inverse(transform)(val - ε)
    return Positive(unconstrained_value, transform, convert(typeof(unconstrained_value), ε))
end

struct Positive{T<:Real,V,Tε<:Real} <: AbstractParameter
    unconstrained_value::T
    transform::V
    ε::Tε
end

value(x::Positive) = x.transform(x.unconstrained_value) + x.ε

function flatten(::Type{T}, x::Positive) where {T<:Real}
    v, unflatten_to_Real = flatten(T, x.unconstrained_value)

    function unflatten_Positive(v_new::Vector{T})
        return Positive(unflatten_to_Real(v_new), x.transform, x.ε)
    end

    return v, unflatten_Positive
end

"""
    bounded(val::Real, lower_bound::Real, upper_bound::Real)

Constructs a `Bounded`.
The `value` of a `Bounded` is a `Real` number that is constrained to be within the interval
(`lower_bound`, `upper_bound`), and is equal to `val`.
This is represented internally in terms of an `unconstrained_value` and a `transform` that
maps any `Real` to this interval.
"""
function bounded(val::Real, lower_bound::Real, upper_bound::Real)
    lb = convert(typeof(val), lower_bound)
    ub = convert(typeof(val), upper_bound)

    # construct open interval
    ε = convert(typeof(val), 1e-12)
    lb_plus_ε = lb + ε
    ub_minus_ε = ub - ε

    if val > ub_minus_ε || val < lb_plus_ε
        throw(
            ArgumentError(
                "Value, $val, outside of specified bounds ($lower_bound, $upper_bound)."
            ),
        )
    end

    length_interval = ub_minus_ε - lb_plus_ε
    unconstrained_val = logit((val - lb_plus_ε) / length_interval)
    transform(x) = lb_plus_ε + length_interval * logistic(x)

    return Bounded(unconstrained_val, lb, ub, transform, ε)
end

struct Bounded{T<:Real,V,Tε<:Real} <: AbstractParameter
    unconstrained_value::T
    lower_bound::T
    upper_bound::T
    transform::V
    ε::Tε
end

value(x::Bounded) = x.transform(x.unconstrained_value)

function flatten(::Type{T}, x::Bounded) where {T<:Real}
    v, unflatten_to_Real = flatten(T, x.unconstrained_value)

    function unflatten_Bounded(v_new::Vector{T})
        return Bounded(
            unflatten_to_Real(v_new), x.lower_bound, x.upper_bound, x.transform, x.ε
        )
    end

    return v, unflatten_Bounded
end


text:
push!(LOAD_PATH,"../src/")
using Documenter
using StreamGraphs

makedocs(
         modules  = [StreamGraphs],
         format   = :html,
         sitename = "StreamGraphs.jl",
         doctest  = false,
         pages    = Any[
                        "Getting Started" => "index.md",
                        "LinkStreams vs. StreamGraphs" => "streamtypes.md"
                        ]
         )

deploydocs(
           deps   = nothing,
           make   = nothing,
           repo   = "github.com/NicolasGensollen/StreamGraphs.jl.git",
           branch = "gh-pages",
           target = "build",
           )


text:
using CodecZlib: GzipDecompressorStream

using CSV: File

using DataFrames: DataFrame, rename

function make_vcf_dataframe(path::String)::DataFrame

    io = open(path)

    io2 = GzipDecompressorStream(io)

    df = DataFrame(File(io2; comment="##", delim='\t', header=1))

    df = rename(df, ("#CHROM" => "CHROM"))

    close(io)

    return df

end

export vcf_dataframe


text:
# # Diffusion Equation in a Vertical Column
# ## Setup
# - A vertical diffusive column split over two Cartesian components (atmos & ocean columns)
# - Boundary conditions used to set diffusive flux in other components have different vertical discretizations and timesteps
# - The coupler maps fields between the two components

# ## Equations
# The RHS evaluation follows:
# ```math
# \frac{\partial \theta}{\partial t} = - \nabla \cdot [\kappa(\phi_{init}) \nabla \phi]
# ```
# where
#  -  $\theta$ is the tracer (e.g. potential temperature)
#  -  $\kappa$ is the diffusivity tensor

# # Import packages
using ClimateMachine
using MPI
using Statistics

## To couple
using ClimaCoupler
using Unitful, Dates

## To create meshes & grids
using ClimateMachine.Ocean.Domains
using ClimateMachine.Grids
import ClimateMachine.DGMethods.NumericalFluxes: NumericalFluxSecondOrder

## To setup some callbacks
using ClimateMachine.GenericCallbacks

## To invoke timestepper
using ClimateMachine.ODESolvers
using ClimateMachine.ODESolvers: solve!
using ClimateMachine.MPIStateArrays: weightedsum

ClimateMachine.init()
const FT = Float64;

if !(:CplTestingBL in names(Main))
    include("CplTestingBL.jl") # allows re-run of script without restarting julia
end
using .CplTestingBL

# # Set simulation parameters

couple_dt = 3600.0 # timestep at which coupled components sync
nstepsA = 10 # atmos steps per coupled timestep
nstepsO = 5 # ocean steps per coupled timestep

##  Haney like relaxation time scale and a length scale (Haney, 1971).
##  Air-sea exchange vigor is governed by length/time-scale.
const τ_airsea = FT(60 * 86400)
const L_airsea = FT(500)
const λ_airsea = FT(L_airsea / τ_airsea)
function coupling_lambda()
    return (λ_airsea)
end;

##  Background atmos and ocean diffusivities
const κᵃʰ = FT(1e4) * 0.0
const κᵃᶻ = FT(1e-1)
const κᵒʰ = FT(1e3) * 0.0
const κᵒᶻ = FT(1e-4)

# # Set up coupled model
# Define component models and initialize the coupler

function main(::Type{FT}) where {FT}
    ## Domain
    Np = 4
    ΩA = RectangularDomain(
        Ne = (10, 10, 5),
        Np = Np,
        x = (0, 1e6),
        y = (0, 1e6),
        z = (0, 1e5),
        periodicity = (true, true, false),
    )
    ΩO = RectangularDomain(
        Ne = (10, 10, 4),
        Np = Np,
        x = (0, 1e6),
        y = (0, 1e6),
        z = (-4e3, 0),
        periodicity = (true, true, false),
    )

    ## Grid
    btags = ((0, 0), (0, 0), (1, 2))
    gridA = DiscontinuousSpectralElementGrid(ΩA; boundary_tags = btags)
    gridO = DiscontinuousSpectralElementGrid(ΩO; boundary_tags = btags)

    ## Numerics-specific options
    numerics = (NFsecondorder = CplTestingBL.PenaltyNumFluxDiffusive(),)

    ## Callbacks (TODO)
    callbacks = ()

    ## Collect spatial info, timestepping, balance law and DGmodel for the two components

    ## 1. Atmos component
    mA = CplModel(;
        grid = gridA,
        equations = CplTestBL(bl_propA, (CoupledPrimaryBoundary(), ExteriorBoundary())),
        nsteps = nstepsA,
        dt = couple_dt / nstepsA,
        numerics...,
    )

    ## 2. Ocean component
    mO = CplModel(;
        grid = gridO,
        equations = CplTestBL(bl_propO, (ExteriorBoundary(), CoupledSecondaryBoundary())),
        nsteps = nstepsO,
        dt = couple_dt / nstepsO,
        numerics...,
    )
    #+
    # Create the coupler object for holding import/export fields and performs mappings
    # and instantiate the coupled timestepper:
    coupler = CplState()
    coupler_register!(coupler, :Ocean_SST, deepcopy(mO.state.θ[mO.boundary]), mO.grid, DateTime(0), u"°C")
    coupler_register!(
        coupler,
        :Atmos_MeanAirSeaθFlux,
        deepcopy(mA.state.F_accum[mA.boundary]),
        mA.grid,
        DateTime(0),
        u"°C",
    )

    compA = (pre_step = preatmos, component_model = mA, post_step = postatmos)
    compO = (pre_step = preocean, component_model = mO, post_step = postocean)
    component_list = (atmosphere = compA, ocean = compO)
    cpl_solver = CplSolver(component_list = component_list, coupler = coupler, coupling_dt = couple_dt, t0 = 0.0)

    return cpl_solver, callbacks
end

function run(cpl_solver, numberofsteps, cbvector)
    solve!(nothing, cpl_solver; numberofsteps = numberofsteps, callbacks = cbvector)
end

# # Define `pre_step` and `post_step` functions
# Each component model must define `pre_step` and `post_step` functions.
# In the `pre_step`, a component imports necessary boundary state and flux data from the coupler.
# In the `post_step`, a component exports boundary data to the coupler to be later received by other components.

function get_components(csolver)
    mA = csolver.component_list.atmosphere.component_model
    mO = csolver.component_list.ocean.component_model
    return mA, mO
end

function preatmos(csolver)
    mA, mO = get_components(csolver)

    ## Set boundary SST used in atmos to SST of ocean surface at start of coupling cycle.
    mA.discretization.state_auxiliary.θ_secondary[mA.boundary] .=
        coupler_get(csolver.coupler, :Ocean_SST, mA.grid, DateTime(0), u"°C")
    ## Set atmos boundary flux accumulator to 0.
    mA.state.F_accum .= 0

    @info(
        "preatmos",
        time = csolver.t,
        total_θ_atmos = weightedsum(mA.state, 1),
        total_θ_ocean = weightedsum(mO.state, 1),
        total_θ = weightedsum(mA.state, 1) + weightedsum(mO.state, 1),
        atmos_θ_surface_max = maximum(mA.state.θ[mA.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
    )
end

function postatmos(csolver)
    mA, mO = get_components(csolver)

    ## Pass atmos exports to "coupler" namespace
    ## 1. Save mean θ flux at the Atmos boundary during the coupling period
    coupler_put!(
        csolver.coupler,
        :Atmos_MeanAirSeaθFlux,
        mA.state.F_accum[mA.boundary] ./ csolver.dt,
        mA.grid,
        DateTime(0),
        u"°C",
    )

    @info(
        "postatmos",
        time = time = csolver.t + csolver.dt,
        total_θ_atmos = weightedsum(mA.state, 1),
        total_θ_ocean = weightedsum(mO.state, 1),
        total_F_accum = mean(mA.state.F_accum[mA.boundary]) * 1e6 * 1e6,
        total_θ = weightedsum(mA.state, 1) + weightedsum(mO.state, 1) + mean(mA.state.F_accum[mA.boundary]) * 1e6 * 1e6,
        F_accum_max = maximum(mA.state.F_accum[mA.boundary]),
        F_avg_max = maximum(mA.state.F_accum[mA.boundary] ./ csolver.dt),
        atmos_θ_surface_max = maximum(mA.state.θ[mA.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
    )
end

function preocean(csolver)
    mA, mO = get_components(csolver)

    ## Set mean air-sea theta flux
    mO.discretization.state_auxiliary.F_prescribed[mO.boundary] .=
        coupler_get(csolver.coupler, :Atmos_MeanAirSeaθFlux, mO.grid, DateTime(0), u"°C")
    ## Set ocean boundary flux accumulator to 0. (this isn't used)
    mO.state.F_accum .= 0

    @info(
        "preocean",
        time = csolver.t,
        F_prescribed_max = maximum(mO.discretization.state_auxiliary.F_prescribed[mO.boundary]),
        F_prescribed_min = maximum(mO.discretization.state_auxiliary.F_prescribed[mO.boundary]),
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
        ocean_θ_surface_min = maximum(mO.state.θ[mO.boundary]),
    )
end

function postocean(csolver)
    mA, mO = get_components(csolver)
    @info(
        "postocean",
        time = csolver.t + csolver.dt,
        ocean_θ_surface_max = maximum(mO.state.θ[mO.boundary]),
        ocean_θ_surface_min = maximum(mO.state.θ[mO.boundary]),
    )

    ## Pass ocean exports to "coupler" namespace
    ##  1. Ocean SST (value of θ at z=0)
    coupler_put!(csolver.coupler, :Ocean_SST, mO.state.θ[mO.boundary], mO.grid, DateTime(0), u"°C")
end

# # Specify balance law

## Set atmosphere initial state function
function atmos_init_theta(xc, yc, zc, npt, el)
    return 30.0
end
## Set atmosphere shadow boundary flux function
function atmos_theta_shadow_boundary_flux(θᵃ, θᵒ, npt, el, xc, yc, zc)
    if zc == 0.0
        tflux = (1.0 / τ_airsea) * (θᵃ - θᵒ)
    else
        tflux = 0.0
    end
    return tflux
end
## Set atmsophere diffusion coeffs
function atmos_calc_kappa_diff(_...)
    return κᵃʰ, κᵃʰ, κᵃᶻ
end
## Set atmos source!
function atmos_source_theta(θᵃ, npt, el, xc, yc, zc, θᵒ)
    tsource = 0.0
    if zc == 0.0
        ## tsource = -(1. / τ_airsea)*( θᵃ-θᵒ )
    end
    return tsource
end
## Set penalty term tau (for debugging)
function atmos_get_penalty_tau(_...)
    return FT(3.0 * 0.0)
end
## Create atmos component
bl_propA = CplTestingBL.prop_defaults()

bl_propA = (; bl_propA..., init_theta = atmos_init_theta, theta_shadow_boundary_flux = atmos_theta_shadow_boundary_flux)
bl_propA = (bl_propA..., init_theta = atmos_init_theta)
bl_propA = (bl_propA..., theta_shadow_boundary_flux = atmos_theta_shadow_boundary_flux)
bl_propA = (bl_propA..., calc_kappa_diff = atmos_calc_kappa_diff)
bl_propA = (bl_propA..., source_theta = atmos_source_theta)
bl_propA = (bl_propA..., get_penalty_tau = atmos_get_penalty_tau)
bl_propA = (bl_propA..., coupling_lambda = coupling_lambda)

## Set initial temperature profile
function ocean_init_theta(xc, yc, zc, npt, el)
    return 20.0
end
## Set boundary source imported from atmos
function ocean_source_theta(θ, npt, el, xc, yc, zc, air_sea_flux_import)
    sval = 0.0
    if zc == 0.0
        ## sval=air_sea_flux_import
    end
    return sval
end
## Set ocean diffusion coeffs
function ocean_calc_kappa_diff(_...)
    ## return κᵒʰ,κᵒʰ,κᵒᶻ*FT(100.)
    return κᵒʰ, κᵒʰ, κᵒᶻ # m^2 s^-1
end
## Set penalty term tau (for debugging)
function ocean_get_penalty_tau(_...)
    return FT(0.15 * 0.0)
end
## Create ocean component
bl_propO = CplTestingBL.prop_defaults()
bl_propO = (bl_propO..., init_theta = ocean_init_theta)
bl_propO = (bl_propO..., source_theta = ocean_source_theta)
bl_propO = (bl_propO..., calc_kappa_diff = ocean_calc_kappa_diff)
bl_propO = (bl_propO..., get_penalty_tau = ocean_get_penalty_tau)
bl_propO = (bl_propO..., coupling_lambda = coupling_lambda)

# # Run simulation
simulation, cbvector = main(Float64);
nsteps = 10
println("Initialized. Running...")
run(simulation, nsteps, cbvector)

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
## sample space/domain

@compat abstract type VariateForm end
type Univariate    <: VariateForm end
type Multivariate  <: VariateForm end
type Matrixvariate <: VariateForm end

@compat abstract type ValueSupport end
type Discrete   <: ValueSupport end
type Continuous <: ValueSupport end

Base.eltype(::Type{Discrete}) = Int
Base.eltype(::Type{Continuous}) = Float64

## Sampleable

@compat abstract type Sampleable{F<:VariateForm,S<:ValueSupport} end

Base.length(::Sampleable{Univariate}) = 1
Base.length(s::Sampleable{Multivariate}) = throw(MethodError(length, (s,)))
Base.length(s::Sampleable) = prod(size(s))

Base.size(s::Sampleable{Univariate}) = ()
Base.size(s::Sampleable{Multivariate}) = (length(s),)

Base.eltype{F,S}(s::Sampleable{F,S}) = eltype(S)
Base.eltype{F}(s::Sampleable{F,Discrete}) = Int
Base.eltype{F}(s::Sampleable{F,Continuous}) = Float64

nsamples{D<:Sampleable{Univariate}}(::Type{D}, x::Number) = 1
nsamples{D<:Sampleable{Univariate}}(::Type{D}, x::AbstractArray) = length(x)
nsamples{D<:Sampleable{Multivariate}}(::Type{D}, x::AbstractVector) = 1
nsamples{D<:Sampleable{Multivariate}}(::Type{D}, x::AbstractMatrix) = size(x, 2)
nsamples{D<:Sampleable{Matrixvariate}}(::Type{D}, x::Number) = 1
nsamples{D<:Sampleable{Matrixvariate},T<:Number}(::Type{D}, x::Array{Matrix{T}}) = length(x)

## Distribution <: Sampleable

@compat abstract type Distribution{F<:VariateForm,S<:ValueSupport} <: Sampleable{F,S} end

@compat const UnivariateDistribution{S<:ValueSupport}   = Distribution{Univariate,S}
@compat const MultivariateDistribution{S<:ValueSupport} = Distribution{Multivariate,S}
@compat const MatrixDistribution{S<:ValueSupport}       = Distribution{Matrixvariate,S}
const NonMatrixDistribution = Union{UnivariateDistribution, MultivariateDistribution}

@compat const DiscreteDistribution{F<:VariateForm}   = Distribution{F,Discrete}
@compat const ContinuousDistribution{F<:VariateForm} = Distribution{F,Continuous}

const DiscreteUnivariateDistribution     = Distribution{Univariate,    Discrete}
const ContinuousUnivariateDistribution   = Distribution{Univariate,    Continuous}
const DiscreteMultivariateDistribution   = Distribution{Multivariate,  Discrete}
const ContinuousMultivariateDistribution = Distribution{Multivariate,  Continuous}
const DiscreteMatrixDistribution         = Distribution{Matrixvariate, Discrete}
const ContinuousMatrixDistribution       = Distribution{Matrixvariate, Continuous}

variate_form{VF<:VariateForm,VS<:ValueSupport}(::Type{Distribution{VF,VS}}) = VF
variate_form{T<:Distribution}(::Type{T}) = variate_form(supertype(T))

value_support{VF<:VariateForm,VS<:ValueSupport}(::Type{Distribution{VF,VS}}) = VS
value_support{T<:Distribution}(::Type{T}) = value_support(supertype(T))

## TODO: the following types need to be improved
@compat abstract type SufficientStats end
@compat abstract type IncompleteDistribution end

@compat const DistributionType{D<:Distribution} = Type{D}
const IncompleteFormulation = Union{DistributionType,IncompleteDistribution}


text:
mutable struct TestMeta <: CachedArrays.AbstractMetadata
    dirty::Bool
    id::UInt
    pool::CachedArrays.Pool
    sibling::Union{Nothing,TestMeta}
end

TestMeta(id) = TestMeta(false, id, CachedArrays.Local, nothing)

# Implement the API
CachedArrays.setdirty!(T::TestMeta, flag) = (T.dirty = flag)
CachedArrays.isdirty(T::TestMeta) = (T.dirty)
CachedArrays.getid(T::TestMeta) = T.id
CachedArrays.getpool(T::TestMeta) = T.pool
CachedArrays.getsibling(T::TestMeta) = T.sibling
CachedArrays.setsibling!(A::TestMeta, B::TestMeta) = (A.sibling = B)

struct ArrayWrapper
    array::Array
end

const TRACKER = IdDict{ArrayWrapper,TestMeta}

CachedArrays.metastyle(::ArrayWrapper) = TestMeta
CachedArrays.metadata(x::ArrayWrapper, ::Type{TestMeta})::TestMeta = TRACKER[x]

@testset "Testing Metadata API" begin
    x = rand(Float32, 100)
    A = ArrayWrapper(x)
    TRACKER[A] = TestMeta(0)

    @test CachedArrays.isdirty(A) == false
    CachedArrays.setdirty!(A, true)
    @test CachedArrays.isdirty(A) == true

    @test CachedArrays.getid(A) == zero(UInt)
    @test CachedArrays.getpool(A) == CachedArrays.Local
    @test isnothing(CachedArrays.getsibling(A))

    # Make another array and make a sibling for the array.
    B = ArrayWrapper(rand(Float32, 100))
    CachedArrays.setsibling!(A, B)

    @test CachedArrays.getsibling(A) === B

    # cleanup
    empty!(TRACKER)
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
"""
    struct PixBox
        x::Int32
        y::Int32
        w::Int32
        h::Int32
    end

The PixBox structure is used to identify a region in an image.

__Values:__

| Name | Description
| :--- | :----------
| x    | The left side of the region identified by this box.
| y    | The top edge of the region identified by this box.
| w    | The width of the box.
| h    | The height of the box.
"""
struct PixBox
    x::Int32
    y::Int32
    w::Int32
    h::Int32
end

# =================================================================================================
"""
    show(
        io::IO,
        box::PixBox
    )::Nothing

Write some summary information about the PixBox to the IO stream.

__Arguments:__

| T | Name | Default | Description
|---| :--- | :------ | :----------
| R | io   |         | The stream to write the information to.
| R | box  |         | The box to display information about.
"""
function Base.show(
            io::IO,
            box::PixBox
        )::Nothing
    println(io, "Box at: ($(box.x), $(box.y)) Dims: ($(box.w), $(box.h))")
    nothing
end


text:
using SEMPP
using Test, Random, Distributions, EGPD, CSV, DataFrames
using Dates

Random.seed!(12)

@testset "SEMPP.jl" begin
    include("Structures_tests.jl")
    include("ParametersEstimation_tests.jl")
    include(joinpath("Plots_tests", "PlotsData_tests.jl"))
end


text:
using FiniteDepthParametricMachines: dense_filtrations,
                          conv_filtrations,
                          recur_filtrations,
                          filtrations,
                          DenseMachine,
                          ConvMachine,
                          RecurMachine,
                          clean,
                          clean!,
                          solve,
                          solve!,
                          derivative,
                          flip

using NNlib, ChainRulesCore, ChainRulesTestUtils
using Test

@testset "utils" begin
    W = rand(3, 3)
    @test flip(W) == permutedims(W)
    W = rand(5, 3, 3)
    @test flip(W) == permutedims(W[end:-1:1, :, :], (1, 3, 2))
    W = rand(5, 5, 3, 3)
    @test flip(W) == permutedims(W[end:-1:1, end:-1:1, :, :], (1, 2, 4, 3))
end

@testset "dense" begin
    # non linear machine
    minibatch = 6
    dims = [2, 4, 7, 1]
    y₀, z₀ = rand(sum(dims), minibatch), rand(sum(dims), minibatch)
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(sum(dims), sum(dims)), tanh
    (forward, backward) = dense_filtrations(y₀, dims)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    # test that forward pass is a machine
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ.(y) + z₀ # z = σ(y) + z₀
    @test y ≈ W̃' * z + y₀ # y = y₀ + W̃' z
    # test that backward pass is a machine
    Dσ = derivative.(z .- z₀, σ, y)
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀ # z = J ⊙ y + z₀
    @test y ≈ W̃ * z + y₀ # y = y₀ + W̃ z
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # linear machine
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(sum(dims), sum(dims)), rand(sum(dims), minibatch)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    y, z = copy(y₀), copy(z₀)
    # test that forward pass is a machine
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ .* y + z₀
    @test y ≈ W̃' * z + y₀
    # test that backward pass is a machine
    Dσ = σ
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ W̃ * z + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # machine versus `solve`
    machine1 = DenseMachine(dims, tanh)
    machine2 = DenseMachine(tanh, dims)
    for m in (machine1, machine2)
        @test m == DenseMachine(m.W, m.σ, m.dims)
        @test m(y₀) ≈ solve(y₀, m.W, m.σ, dense_filtrations(y₀, m.dims))
    end
end

@testset "conv 1D" begin
    # nonlinear machine
    minibatch = 6
    dims = [2, 4, 7, 1]
    pad = (1, 3)
    kernelsize = (5,)
    datasize = (12,)
    y₀, z₀ = rand(datasize..., sum(dims), minibatch), rand(datasize..., sum(dims), minibatch)
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), tanh
    (forward, backward) = conv_filtrations(y₀, dims; pad)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ.(y) + z₀
    @test y ≈ conv(z, W̃; pad) + y₀
    # test that backward pass is a machine
    Dσ = derivative.(z .- z₀, σ, y)
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=pad)) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # linear machine
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), rand(datasize..., sum(dims), minibatch)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ .* y + z₀
    @test y ≈ conv(z, W̃; pad) + y₀
    # test that backward pass is a machine
    Dσ = σ
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=pad)) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # machine versus `solve`
    machine1 = ConvMachine(dims, tanh; pad, init=rand)
    machine2 = ConvMachine(tanh, dims; pad, init=rand)
    for m in (machine1, machine2)
        @test m == ConvMachine(m.W, m.σ, m.dims; pad)
        @test m(y₀) ≈ solve(y₀, m.W, m.σ, conv_filtrations(y₀, m.dims; m.pad))
    end
end

@testset "conv 2D" begin
    # nonlinear machine
    minibatch = 3
    dims = [2, 1, 3, 1]
    pad = (1, 3, 2, 0)
    kernelsize = (pad[1] + pad[2] + 1, pad[3] + pad[4] + 1)
    datasize = (6, 9)
    y₀, z₀ = rand(datasize..., sum(dims), minibatch), rand(datasize..., sum(dims), minibatch)
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), tanh
    (forward, backward) = conv_filtrations(y₀, dims; pad)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ.(y) + z₀
    @test y ≈ conv(z, W̃; pad) + y₀
    # test that backward pass is a machine
    Dσ = derivative.(z .- z₀, σ, y)
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=pad)) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # linear machine
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), rand(datasize..., sum(dims), minibatch)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ .* y + z₀
    @test y ≈ conv(z, W̃; pad) + y₀
    # test that backward pass is a machine
    Dσ = σ
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=pad)) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # machine versus `solve`
    machine1 = ConvMachine(dims, tanh; pad, init=rand)
    machine2 = ConvMachine(tanh, dims; pad, init=rand)
    for m in (machine1, machine2)
        @test m == ConvMachine(m.W, m.σ, m.dims; pad)
        @test m(y₀) ≈ solve(y₀, m.W, m.σ, conv_filtrations(y₀, m.dims; m.pad))
    end
end

@testset "recurrent" begin
    # nonlinear machine
    minibatch = 6
    dims = [2, 4, 7, 1]
    pad = 6
    kernelsize = (pad + 1,)
    timeblock = 5
    datalength = 32 # test `datalength % timeblock != 0`
    y₀, z₀ = rand(datalength, sum(dims), minibatch), rand(datalength, sum(dims), minibatch)
    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), tanh
    (forward, backward) = recur_filtrations(y₀, dims; pad, timeblock)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ.(y) + z₀
    @test y ≈ conv(z, W̃; pad=(pad, 0)) + y₀
    # test that backward pass is a machine
    Dσ = derivative.(z .- z₀, σ, y)
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=(pad, 0))) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    y, z = copy(y₀), copy(z₀)
    W, σ = randn(kernelsize..., sum(dims), sum(dims)), rand(datalength, sum(dims), minibatch)
    @test flip(clean(W, forward)) ≈ clean(flip(W), backward)
    W̃ = clean(W, forward)
    solve!(y, z, W̃, σ, forward)
    @test z ≈ σ .* y + z₀
    @test y ≈ conv(z, W̃; pad=(pad, 0)) + y₀
    # test that backward pass is a machine
    Dσ = σ
    y, z = copy(y₀), copy(z₀)
    solve!(y, z, flip(W̃), Dσ, backward)
    @test z ≈ Dσ .* y + z₀
    @test y ≈ ∇conv_data(z, W̃, DenseConvDims(y, W̃, padding=(pad, 0))) + y₀
    test_rrule(solve, y₀, z₀, W, σ, (forward, backward) ⊢ NoTangent()) # test derivatives

    # machine versus `solve`
    machine1 = RecurMachine(dims, tanh; pad, timeblock, init=rand)
    machine2 = RecurMachine(tanh, dims; pad, timeblock, init=rand)
    for m in (machine1, machine2)
        @test m == RecurMachine(m.W, m.σ, m.dims; pad, timeblock)
        @test m(y₀) ≈ solve(y₀, m.W, m.σ, recur_filtrations(y₀, m.dims; m.pad, m.timeblock))
    end
end

@testset "filtrations nograd" begin
    dimensions = [2, 2, 2]
    minibatch = 32
    y = rand(sum(dimensions), minibatch)
    machine = DenseMachine(dimensions, tanh)
    (forward, backward) = filtrations(machine, y)
    res, pb = rrule(filtrations, machine, y)
    @test res.forward.shapes == forward.shapes
    @test res.backward.shapes == backward.shapes
    @test pb(nothing, nothing) == (NoTangent(), NoTangent(), NoTangent())

    duration, minibatch = 12, 32
    pad = (3, 3)
    y = rand(duration, sum(dimensions), minibatch)
    machine = ConvMachine(dimensions, tanh; pad)
    (forward, backward) = filtrations(machine, y)
    res, pb = rrule(filtrations, machine, y)
    @test res.forward.shapes == forward.shapes
    @test res.backward.shapes == backward.shapes
    @test pb(nothing, nothing) == (NoTangent(), NoTangent(), NoTangent())

    duration, minibatch = 12, 32
    pad, timeblock = 3, 5
    y = rand(duration, sum(dimensions), minibatch)
    machine = RecurMachine(dimensions, tanh; pad, timeblock)
    (forward, backward) = filtrations(machine, y)
    res, pb = rrule(filtrations, machine, y)
    @test res.forward.shapes == forward.shapes
    @test res.backward.shapes == backward.shapes
    @test pb(nothing, nothing) == (NoTangent(), NoTangent(), NoTangent())
end


text:
struct RelaxInterval <: AbstractCacheAttribute end

Base.@kwdef mutable struct IntervalCache{T<:Real} <: AbstractCache
    v::VariableValues{T}                       = VariableValues{T}()
    _set::Vector{Interval{T}}                  = Interval{T}[]
    _subexpression_set::Dict{Int,Interval{T}}  = Dict{Int,Interval{T}}()
    _set_mv_buffer::Vector{Interval{T}}        = Interval{T}[]
    first_eval::Bool                           = true
end
function IntervalCache{T}(::RelaxInterval, n::Int, m::Int, p::Int) where T<:Real
    IntervalCache{T}(_set = zeros(Interval{T}, n), _set_mv_buffer = zeros(Interval{T}, p))
end
function initialize!(c::IntervalCache{T}, g::DirectedTree) where T<:Real
    c._set                 = zeros(Interval{T}, node_count(g))
    c._set_mv_buffer       = zeros(Interval{T}, length(sparsity(g, 1)))
    return
end

set(::RelaxInterval, b::IntervalCache{T}, i) where T<:Real = b._set[i]
num(::RelaxInterval, b::IntervalCache{T}, i) where T<:Real = b._set[i].lo
is_num(::RelaxInterval, b::IntervalCache{T}, i) where T<:Real = b._set[i]
subexpression_set(::RelaxInterval, b::IntervalCache{T}, i) where T<:Real = b._subexpression_set[i]
set_mv_buffer(::RelaxInterval, b::IntervalCache{T}, i) where T<:Real = b._set_mv_buffer[i]

val(b::IntervalCache{T}, i::Int) where T<:Real = val(b.v, i)
lbd(b::IntervalCache{T}, i::Int) where T<:Real = lbd(b.v, i)
ubd(b::IntervalCache{T}, i::Int) where T<:Real = ubd(b.v, i)

function Base.setindex!(b::IntervalCache{T}, v::Interval{T}, i::Int) where T<:Real
    if first_eval
        b._set[i] = v
    else
        b._set[i] = b._set[i] ∩ v
    end
    nothing
end

include(joinpath(@__DIR__, "forward.jl"))
include(joinpath(@__DIR__, "reverse.jl"))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
using Distributions, LinearAlgebra, LaTeXStrings, Random, Plots; pyplot()
Random.seed!(1)

N = 10^5

SigY = [ 6 4 ;
         4 9]
muY = [15 ;
       20]
A = cholesky(SigY).L

rngGens = [()->rand(Normal()),
           ()->rand(Uniform(-sqrt(3),sqrt(3))),
           ()->rand(Exponential())-1]

rv(rg) = A*[rg(),rg()] + muY

data = [[rv(r) for _ in 1:N] for r in rngGens]

stats(data) = begin
    data1, data2 = first.(data),last.(data)
    println(round(mean(data1),digits=2), "\t",round(mean(data2),digits=2),"\t",
            round(var(data1),digits=2), "\t", round(var(data2),digits=2), "\t",
            round(cov(data1,data2),digits=2))
end

println("Mean1\tMean2\tVar1\tVar2\tCov")
for d in data
    stats(d)
end

scatter(first.(data[1]), last.(data[1]), c=:blue, ms=1, msw=0, label="Normal")
scatter!(first.(data[2]), last.(data[2]), c=:red, ms=1, msw=0, label="Uniform")
scatter!(first.(data[3]),last.(data[3]),c=:green, ms=1,msw=0,label="Exponential",
	xlims=(0,40), ylims=(0,40), legend=:bottomright, ratio=:equal,
    xlabel=L"X_1", ylabel=L"X_2")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))" --color=yes --startup-file=no -e 'include(popfirst!(ARGS))' \
    "${BASH_SOURCE[0]}" "$@"
    =#

abstract type Machine end
abstract type TuringMachine <: Machine end
abstract type MachineComponent end
abstract type Programme end
abstract type ProgrammeCompoment end


text:
# Tilting (addition with affine function)

"""
  Tilt(f::ProximableFunction, a::AbstractArray, b::Real)

Given function `f`, returns `g(x) = f(x) + <a,x> + b`.
"""

immutable Tilt{T <: ProximableFunction, S <: AbstractArray, R <: Real} <: ProximableFunction
  f::T
  a::S
  b::R
end

Tilt{T <: ProximableFunction, S <: AbstractArray}(f::T, a::S) = Tilt{T, S, eltype(a)}(f, a, 0.0)

function (g::Tilt){T <: RealOrComplex}(x::AbstractArray{T})
  return g.f(x) + vecdot(g.a, x) + g.b
end

function prox!{T <: RealOrComplex}(g::Tilt, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  v = prox!(g.f, x - gamma*g.a, y, gamma)
  return v + vecdot(g.a, y) + g.b
end

function prox_naive{T <: RealOrComplex}(g::Tilt, x::AbstractArray{T}, gamma::Real=1.0)
  y, v = prox_naive(g.f, x - gamma*g.a, gamma)
  return y, v + vecdot(g.a, y) + g.b
end

is_prox_accurate(f::Tilt) = is_prox_accurate(f.f)

fun_name(f::Tilt) = string("Tilted ", fun_name(f.f))
fun_dom(f::Tilt) = fun_dom(f.f)
fun_expr(f::Tilt) = string(fun_expr(f.f)," + a'x + b")
fun_params(f::Tilt) = string(fun_expr(f.f), ", a = $(typeof(f.a))", ", b = $(f.b)")


text:
module DenseJacobian

"""Finite Differences are far from the best solution to find jacobians but are the easiest to
implement. Hence this solution.

    getJacobian!(jacobian::Array{Float64,2}, f::Function, x::Array{Float64,1})

Here "jacobian" must have row size equal to length of vector returned by Function "f" and column size equal to
length of vector "x"
"""
function getJacobian!(jacobian::Array{Float64,2}, f::Function, x::Array{Float64,1})
    f_col = f(x)
    rows::Int64 = length(f_col)
    cols::Int64 = length(x)
    @assert size(jacobian) == (rows, cols) "The rows must be the length of the function output and the cols the length of the input function"
    h::Float64 = 0.0
    for col ∈ 1:cols
        h = x[col] == 0.0 ? sqrt(eps(1.0)) : sqrt(eps(x[col]))*x[col]
        jacobian[:, col] .= (f(x+[zeros(col-1); h; zeros(cols-col)]) - f_col)/h
        #println((x+[zeros(col-1); h; zeros(cols-col)]))

    end
    return nothing
end

"""getJacobian is an easier to use version of denseJacobian!. This is less efficient than denseJacobian!.
Hence whenever possible denseJacobian! must be used."""
function getJacobian(f::Function, x::Array{Float64,1})
    f_col = f(x)
    rows::Int64 = length(f_col)
    cols::Int64 = length(x)
    jacobian::Array{Float64,2} = zeros(rows, cols)
    getJacobian!(jacobian, f, x)
    return jacobian
end

export getJacobian!, getJacobian

end # module


text:
module UOT_toolbox

using LinearAlgebra

function line(x1::Array{Float64, 1}, x2::Array{Float64, 1})
    t = collect(LinRange(0, 1, 100));
    return [(1-tt)*x1[1] + tt*x2[1] for tt in t], [(1-tt)*x1[2] + tt*x2[2] for tt in t];
end

end

max_stars_repo_name: qdenoyelle/OTMSL

max_issues_repo_name: qdenoyelle/OTMSL

max_forks_repo_name: qdenoyelle/OTMSL


text:
# TODO:
const BIG_ENDIAN = ENDIAN_BOM == 0x01020304

# TODO: endianness

const INIT_STATE = UInt32[0x67452301,
                          0xEFCDAB89,
                          0x98BADCFE,
                          0x10325476,
                          0xC3D2E1F0]

const K0  = 0x00000000
const K1  = 0x5A827999
const K2  = 0x6ED9EBA1
const K3  = 0x8F1BBCDC
const K4  = 0xA953FD4E
const KK0 = 0x50A28BE6
const KK1 = 0x5C4DD124
const KK2 = 0x6D703EF3
const KK3 = 0x7A6D76E9
const KK4 = 0x00000000

@inline ROTL32(num, shift::T) where T = (num << shift) | (num >> (T(32) - shift))
@inline F0(x, y, z) = x ⊻ y ⊻ z
@inline F1(x, y, z) = (x & y) | (~x & z)
@inline F2(x, y, z) = (x | ~y) ⊻ z
@inline F3(x, y, z) = (x & z) | (y & ~z)
@inline F4(x, y, z) = x ⊻ (y | ~z)

const left_p = ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
                 8, 5, 14, 2, 11, 7, 16, 4, 13, 1, 10, 6, 3, 15, 12, 9,
                 4, 11, 15, 5, 10, 16, 9, 2, 3, 8, 1, 7, 14, 12, 6, 13,
                 2, 10, 12, 11, 1, 9, 13, 5, 14, 4, 8, 16, 15, 6, 7, 3,
                 5, 1, 6, 10, 8, 13, 3, 11, 15, 2, 4, 9, 12, 7, 16, 14 )
const left_q = ( 11, 14, 15, 12, 5, 8, 7, 9, 11, 13, 14, 15, 6, 7, 9, 8,
                 7, 6, 8, 13, 11, 9, 7, 15, 7, 12, 15, 9, 11, 7, 13, 12,
                 11, 13, 6, 7, 14, 9, 13, 15, 14, 8, 13, 6, 5, 12, 7, 5,
                 11, 12, 14, 15, 14, 15, 9, 8, 9, 14, 5, 6, 8, 6, 5, 12,
                 9, 15, 5, 11, 6, 8, 13, 12, 5, 12, 13, 14, 11, 8, 5, 6 )
const right_p = ( 6, 15, 8, 1, 10, 3, 12, 5, 14, 7, 16, 9, 2, 11, 4, 13,
                  7, 12, 4, 8, 1, 14, 6, 11, 15, 16, 9, 13, 5, 10, 2, 3,
                  16, 6, 2, 4, 8, 15, 7, 10, 12, 9, 13, 3, 11, 1, 5, 14,
                  9, 7, 5, 2, 4, 12, 16, 1, 6, 13, 3, 14, 10, 8, 11, 15,
                  13, 16, 11, 5, 2, 6, 9, 8, 7, 3, 14, 15, 1, 4, 10, 12 )
const right_q = ( 8, 9, 9, 11, 13, 15, 15, 5, 7, 7, 8, 11, 14, 14, 12, 6,
                  9, 13, 15, 7, 12, 8, 9, 11, 7, 7, 12, 7, 6, 15, 13, 11,
                  9, 7, 15, 11, 8, 6, 6, 14, 12, 13, 5, 14, 13, 13, 7, 5,
                  15, 5, 8, 11, 14, 14, 6, 14, 6, 9, 12, 9, 12, 5, 15, 8,
                  8, 5, 12, 9, 12, 5, 14, 6, 8, 13, 6, 5, 15, 13, 11, 11 )

# Left lane generator macro
macro L(i)
    # @assert i <= 80
    # @assert i >= 1

    # Rotate the words
    ww = (:a, :b, :c, :d, :e)
    a = ww[((81 - i) % 5) + 1]
    b = ww[((82 - i) % 5) + 1]
    c = ww[((83 - i) % 5) + 1]
    d = ww[((84 - i) % 5) + 1]
    e = ww[((85 - i) % 5) + 1]

    f = Symbol("F", div(i - 1, 16))
    k = Symbol("K", div(i - 1, 16))

    r = left_p[i]
    s = left_q[i]

    return esc(quote
        t = $a + $f($b, $c, $d) + $k + unsafe_load(buf, $r)
        $a = ROTL32(UInt32( t), UInt8($s)) + $e
        $c = ROTL32(UInt32($c), UInt8(10))
    end)
end

# Right lane generator macro
macro R(i)
    # @assert i <= 80
    # @assert i >= 1

    # Rotate the words
    ww = (:a, :b, :c, :d, :e)
    a = ww[((81 - i) % 5) + 1]
    b = ww[((82 - i) % 5) + 1]
    c = ww[((83 - i) % 5) + 1]
    d = ww[((84 - i) % 5) + 1]
    e = ww[((85 - i) % 5) + 1]

    f = Symbol("F",  4 - div(i - 1, 16))
    k = Symbol("KK",     div(i - 1, 16))

    r = right_p[i]
    s = right_q[i]

    return esc(quote
        t = $a + $f($b, $c, $d) + $k + unsafe_load(buf, $r)
        $a = ROTL32(UInt32( t), UInt8($s)) + $e
        $c = ROTL32(UInt32($c), UInt8(10))
    end)
end