text:
mutable struct Lininterp

	d          :: Array{Int}    # number of points in each dimension
	n          :: Int 	        # number  of dims
	nfunc      :: Int 	        # number  of different functions to evaluate on the same point
	ifunc      :: Vector{Int} 	# indices of which funcs to evaluate
	hascache   :: Bool 	        # has active cache
	cache      :: Array{Int}	# current index of lower bound of bracket containing x
	infs       :: Array{Float64}	    # current value of inf: lower bound of bracket in xgrid
	sups       :: Array{Float64}	   	# current value of sup: upper bound of bracket in xgrid
	z          :: Array{Float64}    	# current position in bracket: z = (x-x[inf]) / (x[sup] - x[inf])
	hits       :: Int 			# count of cache hits
	miss       :: Int 			# count of cache misses
	grids      :: Array{Vector{Float64}} 	# grids for each dimension
	vals       :: Array{Array{Float64}} 	# arrays with function values on the grid.
	vertex     :: Array{Array{Float64}}  # arrays with current vertices

	# constructor for multiple functions to approximate
	function Lininterp(v::Array{Array{Float64}},g::Array{Vector{Float64}})
		d = [size(v[1])...]
		if !all(map(z->size(z)==size(v[1]),v))
			throw(ArgumentError("all arrays in v must have the same shape"))
		end
		n = length(d)
		nfunc = length(v)
		ifunc = collect(1:nfunc)
		if length(d) != length(g)
			throw(ArgumentError("v must have as many dims as g has grids"))
		end
		for i in 1:length(g)
			if d[i] != length(g[i])
				throw(ArgumentError("size(v,$i) must have the same length as g[$i]"))
			end
			if !issorted(g[i])
				throw(ArgumentError("g[$i] must be sorted!"))
			end
		end
		if n > 4
			throw(ArgumentError("currently only up to 4D implemented"))
		end
		cache = ones(Int,n)
		vertex = Array{Float64}[]
		for i=1:nfunc
			push!(vertex,zeros(2^n))
		end
		return new(d,n,nfunc,ifunc,false,cache,zeros(n),zeros(n),zeros(n),0,0,g,v,vertex)
	end

end

# additional outer constructors
# constructor for single function
function Lininterp(v::Array{Float64},g::Array{Vector{Float64}})
	v1 = Array{Float64}[]
	push!(v1,v)
	Lininterp(v1,g)
end
# constructor for two functions
function Lininterp(v1::Array{Float64},v2::Array{Float64},g::Array{Vector{Float64}})
	v = Array{Float64}[]
	push!(v,v1)
	push!(v,v2)
	Lininterp(v,g)
end

# general utility methods
# =======================

getDims(l::Lininterp) = l.d
getGrids(l::Lininterp) =  l.grids
function resetCache!(l::Lininterp)
	fill!(l.cache,0)
	l.hits = 0
	l.miss = 0
	l.hascache=false
	return nothing
end
function getCache(l::Lininterp)
	l.cache
end
function getCache(l::Lininterp,i::Int)
	l.cache[i]
end
function getCachedVal(l::Lininterp,i::Int)
	l.grids[i][l.cache[i]]
end
function getNextCachedVal(l::Lininterp,i::Int)
	l.grids[i][l.cache[i] + 1]
end

function setValue!(l::Lininterp,which::Int,v::Array)
	if size(v) != l.d
		throw(ArgumentError("size(v) must be equal to size d"))
	end
	vals[which] = v
	return nothing
end

function setGrid!(l::Lininterp,i::Int,g::Vector{Float64})
	if length(g) != l.d[i]
		throw(ArgumentError("new grid g incompatible with dim $i = $(l.d[i])"))
	end
	if !issorted(g)
		throw(ArgumentError("g must be sorted!"))
	end
	l.grids[i] = g
	return nothing
end

function hitmiss(l::Lininterp)
	return (l.hits,l.miss)
end

function getValue(l::Lininterp,x::Vector{Float64})
	# make sure all get evaluated
	l.ifunc = collect(1:l.nfunc)
	if l.n == 1
		eval1D(l,x)
	elseif l.n == 2
		eval2D(l,x)
	elseif l.n == 3
		eval3D(l,x)
	elseif l.n ==4
		eval4D(l,x)
	else
		warn("only up to 4D implemented so far")
	end
end
function getValue(l::Lininterp,y::Float64)
	# make sure all get evaluated
	x = [y]
	l.ifunc = collect(1:l.nfunc)
	if l.n == 1
		eval1D(l,x)
	elseif l.n == 2
		eval2D(l,x)
	elseif l.n == 3
		eval3D(l,x)
	elseif l.n ==4
		eval4D(l,x)
	else
		warn("only up to 4D implemented so far")
	end
end

# method with preallocated output vector y
# and index `which` indicating which functions to evaluate
function getValue!(y::Vector{Float64},l::Lininterp,x::Vector{Float64},which::Vector{Int})
	if maximum(which) > l.nfunc
		throw(ArgumentError("which contains a higher index than there are functions"))
	end
	if length(which) != length(y)
		throw(ArgumentError("output y must be as long as `which` index"))
	end
	# only which get evaluated
	l.ifunc = which
	if l.n == 1
		eval1D!(y,l,x)
	elseif l.n == 2
		eval2D!(y,l,x)
	elseif l.n == 3
		eval3D!(y,l,x)
	elseif l.n ==4
		eval4D!(y,l,x)
	else
		warn("only up to 4D implemented so far")
	end
end

# evaluation methods
# ==================

# finds the inf of x in it's grid for each dimension
# uses caching
# 1. remembers last evaluation, and corresponding locations in grid
# 2. finds the bracket of grid values x in is
# 3. finds position of x in the bracket
function findBracket!(l::Lininterp,x::Vector{Float64})

	# if does not have active cache, search entire interval
	if !l.hascache
		for i in 1:length(x)
			l.miss += 1
			# deal with values out of grid: set to grid bounds
			if x[i] <= l.grids[i][1]
				x[i] = l.grids[i][1]
				l.cache[i] = 1
				l.infs[i] = l.grids[i][1]
				l.sups[i] = l.grids[i][2]
				l.z[i]     = 0.0

			elseif x[i] >= l.grids[i][end]
				x[i] = l.grids[i][end]
				l.cache[i] = l.d[i]-1
				l.infs[i] = l.grids[i][end-1]
				l.sups[i] = l.grids[i][end]
				l.z[i]     = 1.0
			else
				l.cache[i] = searchsortedlast(l.grids[i],x[i],1,l.d[i]-1,Base.Forward)
				l.infs[i]  = l.grids[i][l.cache[i]]
				l.sups[i]  = l.grids[i][l.cache[i] + 1]
				l.z[i]     = (x[i] - l.infs[i]) / (l.sups[i] - l.infs[i])
			end
		end
		# now l has a cache
		l.hascache = true

	# else, search cache
	else
		for i in 1:length(x)
			# deal with values out of grid: set to grid bounds
			if x[i] <= l.grids[i][1]
				l.miss += 1
				x[i] = l.grids[i][1]
				l.cache[i] = 1
				l.infs[i] = l.grids[i][1]
				l.sups[i] = l.grids[i][2]
				l.z[i]     = 0.0

			elseif x[i] >= l.grids[i][end]
				l.miss += 1
				x[i] = l.grids[i][end]
				l.cache[i] = l.d[i]-1
				l.infs[i] = l.grids[i][end-1]
				l.sups[i] = l.grids[i][end]
				l.z[i]     = 1.0
			# if x is below current lower bound of bracket, search below
			elseif x[i] < getCachedVal(l,i)
				l.miss += 1
				l.cache[i] = searchsortedlast(l.grids[i],x[i],1,l.cache[i]-1,Base.Forward)
				l.infs[i]  = l.grids[i][l.cache[i]]
				l.sups[i]  = l.grids[i][l.cache[i] + 1]
				l.z[i]     = (x[i] - l.infs[i]) / (l.sups[i] - l.infs[i])
			# if x is above current upper bound of bracket, search above
			elseif x[i] >= getNextCachedVal(l,i)
				l.miss += 1
				l.cache[i] = searchsortedlast(l.grids[i],x[i],l.cache[i]+1,l.d[i]-1,Base.Forward)
				l.infs[i]  = l.grids[i][l.cache[i]]
				l.sups[i]  = l.grids[i][l.cache[i] + 1]
				l.z[i]     = (x[i] - l.infs[i]) / (l.sups[i] - l.infs[i])
			# x is in current bracket. cool!
			else
				l.hits += 1
				l.z[i]   = (x[i] - l.infs[i]) / (l.sups[i] - l.infs[i])
			end
		end
	end
	return nothing
end

# 1D evaluation
# =============

function find1DVertices!(l::Lininterp)
	for i in l.ifunc
		l.vertex[i][1] = l.vals[i][l.cache[1]   ]	# v0
		l.vertex[i][2] = l.vals[i][l.cache[1] + 1]	# v1
	end
end

function eval1D(l::Lininterp,x::Vector{Float64})

	if length(x) != 1
		throw(ArgumentError("x needs 1 elements: one for each D in 1D!"))
	end
	if l.n != 1
		throw(ArgumentError("must supply a Lininterp with 1D"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find1DVertices!(l)

	# build up linear combinations
	out = Float64[]
	cc = 0.0
	for i in l.ifunc
		cc = (1.0-l.z[1])* l.vertex[i][1] +   # v0
			 (    l.z[1])* l.vertex[i][2]    # v1
        push!(out,cc)
    end
    return out
end

function eval1D!(y::Vector{Float64},l::Lininterp,x::Vector{Float64})

	if length(x) != 1
		throw(ArgumentError("x needs 1 elements: one for each D in 1D!"))
	end
	if l.n != 1
		throw(ArgumentError("must supply a Lininterp with 1D"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find1DVertices!(l)

	# build up linear combinations
	yi = 0
	for i in l.ifunc
		yi += 1
		cc = (1.0-l.z[1])* l.vertex[i][1] +   # v0
			 (    l.z[1])* l.vertex[i][2]    # v1
        y[yi] = cc
    end
    return nothing
end

# 2D evaluation
# =============

function find2DVertices!(l::Lininterp)
	for i in l.ifunc
		l.vertex[i][1] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1)]	# v00
		l.vertex[i][2] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1)]	# v10
		l.vertex[i][3] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]  )]	# v01
		l.vertex[i][4] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]  )]	# v11
	end
end

function eval2D(l::Lininterp,x::Array{Float64,1})

	if length(x) != 2
		throw(ArgumentError("x needs 2 elements: one for each D in 2D!"))
	end
	if l.n != 2
		throw(ArgumentError("must supply a Lininterp with 2D"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find2DVertices!(l)

	# build up linear combinations
	out = Float64[]
	cc = 0.0
	for i in l.ifunc
		cc = (1.0-l.z[1])*(1.0-l.z[2]) * l.vertex[i][1] +   # v00
			 (    l.z[1])*(1.0-l.z[2]) * l.vertex[i][2] +   # v10
			 (1.0-l.z[1])*(    l.z[2]) * l.vertex[i][3] +   # v01
	         (    l.z[1])*(    l.z[2]) * l.vertex[i][4]     # v11
        push!(out,cc)
    end
    return out
end

function eval2D!(y::Vector{Float64},l::Lininterp,x::Array{Float64,1})

	if length(x) != 2
		throw(ArgumentError("x needs 2 elements: one for each D in 2D!"))
	end
	if l.n != 2
		throw(ArgumentError("must supply a Lininterp with 2D"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find2DVertices!(l)

	# build up linear combinations
	yi = 0
	for i in l.ifunc
		yi += 1
		cc = (1.0-l.z[1])*(1.0-l.z[2]) * l.vertex[i][1] +   # v00
			 (    l.z[1])*(1.0-l.z[2]) * l.vertex[i][2] +   # v10
			 (1.0-l.z[1])*(    l.z[2]) * l.vertex[i][3] +   # v01
	         (    l.z[1])*(    l.z[2]) * l.vertex[i][4]     # v11
        y[yi] = cc
    end
    return nothing
end

# 3D evaluation
# =============

# finds the function values at cartesian grid of indices
# xinf, xsup, yinf, ysup, zinf, zsup
function find3DVertices!(l::Lininterp)
	for i in l.ifunc
		l.vertex[i][1] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1 + l.d[2]*(l.cache[3]-1))]	# v000
		l.vertex[i][2] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*(l.cache[3]-1))]	# v100
		l.vertex[i][3] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*(l.cache[3]-1))]	# v010
		l.vertex[i][4] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1 + l.d[2]*(l.cache[3]  ))]	# v001
		l.vertex[i][5] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*(l.cache[3]-1))]	# v110
		l.vertex[i][6] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*(l.cache[3]  ))]	# v101
		l.vertex[i][7] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*(l.cache[3]  ))]	# v011
		l.vertex[i][8] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*(l.cache[3]  ))]	# v111
	end
end

function eval3D(l::Lininterp,x::Array{Float64,1})

	if length(x) != 3
		throw(ArgumentError("x needs 3 elements: one for each D in 3D!"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)
	find3DVertices!(l)

	# build up linear combinations
	out = Float64[]
	cc = 0.0
	for i in l.ifunc
		cc = (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3]) * l.vertex[i][1] +   # v000
			 (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3]) * l.vertex[i][2] +   # v100
			 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3]) * l.vertex[i][3] +   # v010
	         (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3]) * l.vertex[i][4] +   # v001
	         (    l.z[1])*(    l.z[2])*(1.0-l.z[3]) * l.vertex[i][5] +   # v110
	         (    l.z[1])*(1.0-l.z[2])*(    l.z[3]) * l.vertex[i][6] +   # v101
	         (1.0-l.z[1])*(    l.z[2])*(    l.z[3]) * l.vertex[i][7] +   # v011
	         (    l.z[1])*(    l.z[2])*(    l.z[3]) * l.vertex[i][8]     # v111
        push!(out,cc)
	end
	return out
end

function eval3D!(y::Vector{Float64},l::Lininterp,x::Array{Float64,1})

	if length(x) != 3
		throw(ArgumentError("x needs 3 elements: one for each D in 3D!"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)
	find3DVertices!(l)

	# build up linear combinations
	yi = 0
	for i in l.ifunc
		yi += 1
		cc = (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3]) * l.vertex[i][1] +   # v000
			 (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3]) * l.vertex[i][2] +   # v100
			 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3]) * l.vertex[i][3] +   # v010
	         (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3]) * l.vertex[i][4] +   # v001
	         (    l.z[1])*(    l.z[2])*(1.0-l.z[3]) * l.vertex[i][5] +   # v110
	         (    l.z[1])*(1.0-l.z[2])*(    l.z[3]) * l.vertex[i][6] +   # v101
	         (1.0-l.z[1])*(    l.z[2])*(    l.z[3]) * l.vertex[i][7] +   # v011
	         (    l.z[1])*(    l.z[2])*(    l.z[3]) * l.vertex[i][8]     # v111
        y[yi] = cc
	end
	return nothing
end

# 4D evaluation
# =============

function find4DVertices!(l::Lininterp)
	for i in l.ifunc
		l.vertex[i][1]  = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]-1)))]	# v0000
		l.vertex[i][2]  = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]-1)))]	# v1000
		l.vertex[i][3]  = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]-1)))]	# v0100
		l.vertex[i][4]  = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]-1)))]	# v0010
		l.vertex[i][5]  = l.vals[i][l.cache[1] +   + l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]  )))]	# v0001
		l.vertex[i][6]  = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]-1)))]	# v1100
		l.vertex[i][7]  = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]-1)))]	# v1010
		l.vertex[i][8]  = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]  )))]	# v1001
		l.vertex[i][9]  = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]-1)))]	# v0110
		l.vertex[i][10] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]  )))]	# v0101
		l.vertex[i][11] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]  )))]	# v0011
		l.vertex[i][12] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]-1)))]	# v1110
		l.vertex[i][13] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]-1) + l.d[3]*(l.cache[4]  )))]	# v1101
		l.vertex[i][14] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]-1 + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]  )))]	# v1011
		l.vertex[i][15] = l.vals[i][l.cache[1] +     l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]  )))]	# v0111
		l.vertex[i][16] = l.vals[i][l.cache[1] + 1 + l.d[1]*(l.cache[2]   + l.d[2]*((l.cache[3]  ) + l.d[3]*(l.cache[4]  )))]	# v1111
	end
end

function eval4D(l::Lininterp,x::Array{Float64,1})

	if length(x) != 4
		throw(ArgumentError("x needs 4 elements: one for each D in 4D!"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find4DVertices!(l)

	wgts = hcat((1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]))

	# build up linear combinations
	m = length(l.ifunc)
	out = zeros(m)
	for i in 1:m
		# cc = (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][1]  +   # v0000
		# 	 (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][2]  +   # v1000
		# 	 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][3]  +   # v0100
	 #         (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][4]  +   # v0010
	 #         (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][5]  +   # v0001
	 #         (    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][6]  +   # v1100
	 #         (    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][7]  +   # v1010
	 #         (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][8]  +   # v1001
		#      (1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][9]  +   # v0110
		# 	 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][10] +   # v0101
		# 	 (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][11] +   # v0011
	 #         (    l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][12] +   # v1110
	 #         (    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][13] +   # v1101
	 #         (    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][14] +   # v1011
	 #         (1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][15] +   # v0111
	 #         (    l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][16]     # v1111
        # push!(out,cc)
      	out[i] = (wgts * l.vertex[l.ifunc[i]])[1]
	end
	return out
end

function eval4D!(y::Vector{Float64},l::Lininterp,x::Array{Float64,1})

	if length(x) != 4
		throw(ArgumentError("x needs 4 elements: one for each D in 4D!"))
	end

	# find in which bracket of grid values x is in.
	# using cached values
	# finds the relative position of x inside the bracket
	findBracket!(l,x)

	# get the function values on the box
	find4DVertices!(l)
	wgts = hcat((1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]),
	(    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]),
	(    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]),
	(1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]),
	(    l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]))

	# build up linear combinations
	yi = 0
	for i in l.ifunc
		yi += 1
		# cc = (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][1]  +   # v0000
		# 	 (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][2]  +   # v1000
		# 	 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][3]  +   # v0100
	 #         (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][4]  +   # v0010
	 #         (1.0-l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][5]  +   # v0001
	 #         (    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(1.0-l.z[4]) * l.vertex[i][6]  +   # v1100
	 #         (    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][7]  +   # v1010
	 #         (    l.z[1])*(1.0-l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][8]  +   # v1001
		#      (1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][9]  +   # v0110
		# 	 (1.0-l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][10] +   # v0101
		# 	 (1.0-l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][11] +   # v0011
	 #         (    l.z[1])*(    l.z[2])*(    l.z[3])*(1.0-l.z[4]) * l.vertex[i][12] +   # v1110
	 #         (    l.z[1])*(    l.z[2])*(1.0-l.z[3])*(    l.z[4]) * l.vertex[i][13] +   # v1101
	 #         (    l.z[1])*(1.0-l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][14] +   # v1011
	 #         (1.0-l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][15] +   # v0111
	 #         (    l.z[1])*(    l.z[2])*(    l.z[3])*(    l.z[4]) * l.vertex[i][16]     # v1111
        y[yi] = (wgts * l.vertex[i])[1]
	end
	return nothing
end

function show(io::IO, l::Lininterp)
	print(io,"linear interpolation object\n")
	print(io,"dimensions: $(l.d)\n")
	print(io,"approximates $(l.nfunc) functions at given point\n")
	print(io,"has active cache: $(l.hascache)\n" )
	print(io,"currently evaluates functions $(l.ifunc)\n" )
	if l.hits+l.miss == 0
		print(io,"accelerator hits 0% of attemps\n")
	else
		print(io,"accelerator hits $(round(Int,100*l.hits/(l.hits+l.miss))) % of attemps\n")
	end
end


text:
function CroutLU(A::Array)

  #A = the matrix of coefficients 'A' must be squared

    m,n =size(A)

    if m != n
        error("Matrix must be squared")
    end

    # System of equations
    L=zeros(m,m)
    U=Matrix{Float64}(I,m,m)

    for i=1:m
        L[i,1]=A[i,1]
    end

    for j=2:m
       U[1,j]=A[1,j]/L[1,1]
    end

    for j=2:m-1
       for i=j:m
           L[i,j]=A[i,j]- (L[i,1:j-1]'*U[1:j-1,j])
        end

        for k=j:m
           U[j,k]=(A[j,k]-(L[j,1:j-1]'*U[1:j-1,k]))/L[j,j]
        end
    end

       L[m,m]=A[m,m]-(L[m,1:m-1]'*U[1:m-1,m])

    return L, U
end

ext: jl


text:
function StatsBase.addcounts!(cm::Dict{U,Int}, x::AbstractDataArray{T}) where {T,U}
    for v in x
        cm[v] = get(cm, v, 0) + 1
    end
    return cm
end

function StatsBase.addcounts!(cm::Dict{U,W}, x::AbstractDataArray{T}, wv::Weights{W}) where {T,U,W}
    n = length(x)
    length(wv) == n || raise_dimerror()
    w = values(wv)
    z = zero(W)

    @inbounds for i in 1:n
        xi = x[i]
        wi = w[i]
        cm[xi] = get(cm, xi, z) + wi
    end
    return cm
end

function StatsBase.countmap(x::AbstractDataArray{T}) where {T}
    addcounts!(Dict{Union{T,Missing}, Int}(), x)
end

function StatsBase.countmap(x::AbstractDataArray{T}, wv::Weights{W}) where {T,W}
    addcounts!(Dict{Union{T,Missing}, W}(), x, wv)
end

"""
    cut(x::AbstractVector, breaks::Vector) -> PooledDataArray
    cut(x::AbstractVector, ngroups::Integer) -> PooledDataArray

Divide the range of `x` into intervals based on the cut points specified in `breaks`,
or into `ngroups` intervals of approximately equal length.

# Examples

```jldoctest
julia> cut([1, 2, 3, 4], [1, 3])
4-element DataArrays.PooledDataArray{String,UInt32,1}:
 "[1,3]"
 "[1,3]"
 "[1,3]"
 "(3,4]"
```
"""
function cut(x::AbstractVector{S}, breaks::Vector{T}) where {S, T}
    if !issorted(breaks)
        sort!(breaks)
    end
    min_x, max_x = minimum(x), maximum(x)
    if breaks[1] > min_x
        pushfirst!(breaks, min_x)
    end
    if breaks[end] < max_x
        push!(breaks, max_x)
    end
    refs = fill(zero(DEFAULT_POOLED_REF_TYPE), length(x))
    for i in 1:length(x)
        if x[i] == min_x
            refs[i] = 1
        else
            refs[i] = searchsortedfirst(breaks, x[i]) - 1
        end
    end
    n = length(breaks)
    from = map(x -> sprint(showcompact, x), breaks[1:(n - 1)])
    to = map(x -> sprint(showcompact, x), breaks[2:n])
    pool = Vector{String}(n - 1)
    if breaks[1] == min_x
        pool[1] = string("[", from[1], ",", to[1], "]")
    else
        pool[1] = string("(", from[1], ",", to[1], "]")
    end
    for i in 2:(n - 1)
        pool[i] = string("(", from[i], ",", to[i], "]")
    end
    PooledDataArray(RefArray(refs), pool)
end

cut(x::AbstractVector, ngroups::Integer) = cut(x, quantile(x, collect(1 : ngroups - 1) / ngroups))

function Base.repeat(A::DataArray{T,N};
                     inner = ntuple(x->1, ndims(A)),
                     outer = ntuple(x->1, ndims(A))) where {T,N}
    DataArray{T,N}(repeat(A.data; inner=inner, outer=outer),
                   BitArray(repeat(A.na; inner=inner, outer=outer)))
end

function Base.repeat(A::PooledDataArray{T,R,N};
                     inner = ntuple(x->1, ndims(A)),
                     outer = ntuple(x->1, ndims(A))) where {T,R,N}
    PooledDataArray(RefArray{R,N}(repeat(A.refs; inner=inner, outer=outer)),
                    A.pool)
end


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project

using DataFrames, CSV, StatsBase

df = CSV.read("../externals/Core.Math.Data/data/Pejcic_318.csv", copycols = true)

println(corspearman(df.ATV, df.ATT))


text:
# Linear Programming
# ==================
const ElLPApproach = Cuint

const LP_ADMM = ElLPApproach(0)
const LP_MEHROTRA = ElLPApproach(1)

const ElKKTSystem = Cuint

const FULL_KKT = ElKKTSystem(0)
const AUGMENTED_KKT = ElKKTSystem(1)
const NORMAL_KKT = ElKKTSystem(2)

struct MehrotraCtrl{T<:ElFloatType}
    primalInit::ElBool
    dualInit::ElBool
    minTol::T
    targetTol::T
    maxIts::ElInt
    maxStepRatio::T
    system::ElKKTSystem
    mehrotra::ElBool
    forceSameStep::ElBool
    solveCtrl::RegSolveCtrl{T}
    resolveReg::ElBool
    outerEquil::ElBool
    basisSize::ElInt
    print::ElBool
    time::ElBool
    wSafeMaxNorm::T
    wMaxLimit::T
    ruizEquilTol::T
    ruizMaxIter::ElInt
    diagEquilTol::T
    checkResiduals::ElBool

    reg0Tmp::T
    reg1Tmp::T
    reg2Tmp::T
    reg0Perm::T
    reg1Perm::T
    reg2Perm::T
end
function MehrotraCtrl(::Type{T};
    primalInit::Bool = false,
    dualInit::Bool = false,
    minTol = eps(T)^0.3,
    targetTol = eps(T)^0.5,
    maxIts = 1000,
    maxStepRatio = 0.99,
    system = FULL_KKT,
    mehrotra = true,
    forceSameStep = true,
    solveCtrl = RegSolveCtrl(T),
    resolveReg = true,
    outerEquil::Bool = true,
    basisSize = 6,
    print = false,
    time = false,
    wSafeMaxNorm = eps(T)^(-0.15),
    wMaxLimit = eps(T)^(-0.4),
    ruizEquilTol = eps(T)^(-0.25),
    ruizMaxIter = 3,
    diagEquilTol = eps(T)^(-0.15),
    checkResiduals = false,

    reg0Tmp = eps(T)^0.25,
    reg1Tmp = eps(T)^0.25,
    reg2Tmp = eps(T)^0.25,

    reg0Perm = eps(T)^0.35,
    reg1Perm = eps(T)^0.35,
    reg2Perm = eps(T)^0.35) where {T<:ElFloatType}

    MehrotraCtrl{T}(ElBool(primalInit),
                    ElBool(dualInit),
                    T(minTol),
                    T(targetTol),
                    ElInt(maxIts),
                    T(maxStepRatio),
                    ElKKTSystem(system),
                    ElBool(mehrotra),
                    ElBool(forceSameStep),
                    solveCtrl,
                    ElBool(resolveReg),
                    ElBool(outerEquil),
                    ElInt(basisSize),
                    ElBool(print),
                    ElBool(time),
                    T(wSafeMaxNorm),
                    T(wMaxLimit),
                    T(ruizEquilTol),
                    ElInt(ruizMaxIter),
                    T(diagEquilTol),
                    ElBool(checkResiduals),
                    T(reg0Tmp),
                    T(reg1Tmp),
                    T(reg2Tmp),
                    T(reg0Perm),
                    T(reg1Perm),
                    T(reg2Perm))
end

struct LPAffineCtrl{T<:ElFloatType}
    approach::Cuint
    mehrotraCtrl::MehrotraCtrl{T}
end

function LPAffineCtrl(::Type{T};
    approach::Cuint = LP_MEHROTRA,
    mehrotraCtrl::MehrotraCtrl = MehrotraCtrl(T)) where {T<:ElFloatType}

    LPAffineCtrl{T}(approach, mehrotraCtrl)
end

for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d))
    @eval begin
        function LPAffine(
            A::DistSparseMatrix{$elty},
            G::DistSparseMatrix{$elty},
            b::DistMultiVec{$elty},
            c::DistMultiVec{$elty},
            h::DistMultiVec{$elty},
            x::DistMultiVec{$elty},
            y::DistMultiVec{$elty},
            z::DistMultiVec{$elty},
            s::DistMultiVec{$elty},
            ctrl::LPAffineCtrl=SOCPAffineCtrl($elty))
            ElError(ccall(($(string("ElLPAffine_", ext)), libEl), Cuint,
                (Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},
                 Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},
                 LPAffineCtrl{$elty}),
                A.obj, G.obj, b.obj, c.obj, h.obj,
                x.obj, y.obj, z.obj, s.obj, ctrl))
            return nothing
        end
    end
end

# Second-Order Cone Programming
# =============================
const SOCP_ADMM = Cuint(0)
const SOCP_IPF = Cuint(1)
const SOCP_IPF_SELFDUAL = Cuint(2)
const SOCP_MEHROTRA = Cuint(3)
const SOCP_MEHROTRA_SELFDUAL = Cuint(4)

struct SOCPAffineCtrl{T<:ElFloatType}
    approach::Cuint
    mehrotraCtrl::MehrotraCtrl{T}
end
function SOCPAffineCtrl(::Type{T};
                      approach::Cuint = SOCP_MEHROTRA,
                      mehrotraCtrl::MehrotraCtrl = MehrotraCtrl(T)) where {T<:ElFloatType}

    SOCPAffineCtrl{T}(approach, mehrotraCtrl)
end

for (elty, ext) in ((:Float32, :s),
                    (:Float64, :d))
    @eval begin
        function SOCPAffine(
            A::DistSparseMatrix{$elty},
            G::DistSparseMatrix{$elty},
            b::DistMultiVec{$elty},
            c::DistMultiVec{$elty},
            h::DistMultiVec{$elty},
            orders::DistMultiVec{ElInt},
            firstInds::DistMultiVec{ElInt},
            labels::DistMultiVec{ElInt},
            x::DistMultiVec{$elty},
            y::DistMultiVec{$elty},
            z::DistMultiVec{$elty},
            s::DistMultiVec{$elty},
            ctrl::SOCPAffineCtrl=SOCPAffineCtrl($elty))

            ElError(ccall(($(string("ElSOCPAffine_", ext)), libEl), Cuint,
                (Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},
                 Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},
                 Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},
                 SOCPAffineCtrl{$elty}),
                A.obj, G.obj, b.obj, c.obj, h.obj,
                orders.obj, firstInds.obj, labels.obj,
                x.obj, y.obj, z.obj, s.obj, ctrl))
            return nothing
        end
    end
end


text:
using Test

using AutomotiveDrivingModels
using EzXML
using AutoUrban

include("test_convert2xodr.jl")
# broken
# include("test_drivermodel.jl")


text:
#=

  Nontransitive dice in ConstraintSolver.jl

  From
  http://en.wikipedia.org/wiki/Nontransitive_dice
  """
  A set of nontransitive dice is a set of dice for which the relation
  "is more likely to roll a higher number" is not transitive. See also
  intransitivity.

  This situation is similar to that in the game Rock, Paper, Scissors,
  in which each element has an advantage over one choice and a
  disadvantage to the other.
  """

  From the Wikipedia page:
    A > B = 5/9   B > C = 5/9  C > A = 5/9
     [2,2,4,4,9,9] die A
     [1,1,6,6,8,8] die B
     [3,3,5,5,7,7] die C

  Here are some examples:

  * The minimum die value for a 3 x 6 dice problem is 4, e.g.
    dice:
    [1, 1, 3, 4, 4, 4]
    [3, 3, 3, 3, 3, 4]
    [2, 2, 2, 4, 4, 4]
    competitions:
    1 vs 2:[15, 13]
    2 vs 3:[18, 15]
    3 vs 1:[15, 12]
    max_dice_val:4
    max_win:18

  * Here's an example of a 6 x 6 dice configuration:
  dice:
    [1, 5, 5, 6, 6, 6]
    [5, 5, 5, 5, 6, 6]
    [4, 4, 5, 6, 6, 6]
    [3, 3, 6, 6, 6, 6]
    [2, 2, 6, 6, 6, 6]
    [1, 4, 6, 6, 6, 6]
    competitions:
    1 vs 2:[12, 10]
    2 vs 3:[14, 12]
    3 vs 4:[12, 11]
    4 vs 5:[12, 8]
    5 vs 6:[10, 9]
    6 vs 1:[13, 10]
    max_dice_val:6
    max_win:14

   * And 3 dice with 12 sides (75s):
    dice:
    [1, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 6]
    [3, 3, 3, 4, 5, 6, 6, 6, 6, 6, 6, 6]
    [2, 2, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6]
    competitions:
    1 vs 2:[47, 44]
    2 vs 3:[52, 46]
    3 vs 1:[52, 50]
    max_dice_val:6
    max_win:52

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

# m: number of dice
# n: number of sides of each die
function nontransitive_dice(m=3,n=6,max_dice_val=6,print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    # The dice
    # m: number dice
    # n: sides per die
    println("number of dice: $m")
    println("number of sides per die: $n")
    @variable(model, 1 <= dice[1:m,1:n] <= max_dice_val, Int)

    # the competitions:
    # (die 1 vs die 2, die 2 vs die 1), ... (die m vs die 1, die 1 vs die m)
    # And the first die must beat the second in each round.
    # Note the last wrap around which breaks the transitivity.
    @variable(model, 0 <= comp[1:m,1:2] <= n*n, Int)

    # Symmetry breaking
    @constraint(model, dice[1,1] == 1)

    # Symmetry breaking: order the number of each die (increasing)
    for dd in 1:m
        increasing(model, dice[dd,:])
    end

    # and now we roll...
    # Number of wins for [A vs B, B vs C, C vs A]
    b1 = @variable(model, [1:m,1:n,1:n], binary=true)
    b2 = @variable(model, [1:m,1:n,1:n], binary=true)
    for dd in 1:m
        # which dice to compare?
        d1 = dd
        d2 = 1+mod(dd,m)

        # This don't work (yet)!
        # @constraint(model, comp[dd,1] == sum([dice[d1, r1] > dice[d2, r2]
        #                for r1 in 1:n, r2 in 1:n]))

        # @constraint(model, comp[dd,2] == sum([dice[d2, r1] > dice[d1, r2]
        #                                for r1 in 1:n, r2 in 1:n]))

        # Rewrite:
        println("match $dd: d1:$d1 vs d2:$d2")
        for r1 in 1:n, r2 in 1:n
            # d1 > d2
            @constraint(model, b1[dd,r1,r2] := {dice[d1, r1] > dice[d2, r2]})
            # d2 > d1
            @constraint(model, b2[dd,r1,r2] := {dice[d2, r1] > dice[d1, r2]})
        end
        @constraint(model, comp[dd,1] == sum(b1[dd,:,:]))
        @constraint(model, comp[dd,2] == sum(b2[dd,:,:]))
    end

    # non-transitivity
    # All dice 1..m-1 must beat the follower, and die m must beat die 1
    for dd in 1:m
        @constraint(model, comp[dd,1] > comp[dd,2])
    end
    @constraint(model, comp[m,1] > comp[1,2])

    # symmetry breaking
    # TODO! (lex_lesseq is defined in constraints_util.jl but either wrong or too slow!)
    # for dd in 2:m
    #    lex_less_eq(model,dice[dd-1,:],dice[dd,:])
    # end

    # Too slow!
    # @variable(model, 1 <= max_val <= max_dice_val, Int)
    # @constraint(model, max_val .>= dice)
    # @objective(model, Min, max_val)

    # Solve the problem
    println("solve")
    optimize!(model)

    status = JuMP.termination_status(model)
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # max_val_val = convert.(Integer,JuMP.value.(max_val; result=sol))
                # println("max_val:$max_val_val", )
                dice_val = convert.(Integer,JuMP.value.(dice; result=sol))
                println("dice:")
                for r in eachrow(dice_val)
                    println(r)
                end
                println("competitions:")
                comp_val = convert.(Integer,JuMP.value.(comp; result=sol))
                for i in 1:m
                    println("$i vs $(1+mod(i,m)):", comp_val[i,:])
                end
                println("max_dice_val:", maximum(dice_val))
                println("max_win:", maximum(comp_val))
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time nontransitive_dice(3,6,6,true,false,6)
# @time nontransitive_dice(3,6,4,true,false,6)
# @time nontransitive_dice(4,6,6,true,false,10)
# @time nontransitive_dice(5,6,6,true,false,20)
@time nontransitive_dice(6,6,6,true,false,20)
# @time nontransitive_dice(3,12,6,true,false,120) # 3 12 sided dice

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using CellwiseEditingDifferentiation
using Distributions
using Lora

ncells = 20

λ = Hyperparameter(:λ)

μ = Data(:μ)
coverage = Data(:coverage)
edited = Data(:edited)
rate = Data(:rate)

function wlogtarget(w::Float64, s::Vector)
  v = s[2]*(1-s[2])*inv_logit(w)
  a, b = beta_pars_from_mv(s[2], v)
  mapreduce(x -> logpdf(Beta(a, b), x), +, s[7:(6+ncells)])+vpcprior(s[2], v, a, b, s[1])+logdvdw(w, s[2])
end

w = BasicContUnvParameter(:w, logtarget=wlogtarget, nkeys=6+ncells)

psetpdf = Array(Function, ncells)
for i in 1:ncells
  psetpdf[i] = function (p::Float64, s::Vector)
    a, b = beta_pars_from_mv(s[2], s[2]*(1-s[2])*inv_logit(s[6]))
    Beta(a+s[4][i], b+s[3][i]-s[4][i])
  end
end

p = Array(BasicContUnvParameter, ncells)
for i in 1:ncells
  p[i] = BasicContUnvParameter(Symbol('p', i), setpdf=psetpdf[i], nkeys=6+ncells)
end

model = GenericModel([λ, μ, coverage, edited, rate, w, p...], isindexed=false)

init = Dict(
  :λ => 10.,
  :μ => 0.5,
  :coverage => fill(10.0, 20),
  :edited => rand([1., 2., 3., 8., 9.], 20),
  :w => 0.
)
init[:rate] = float(init[:edited]./init[:coverage])
for i in 1:ncells
  init[Symbol('p', i)] =
    rand(psetpdf[i](NaN, [[init[k] for k in [:λ, :μ, :coverage, :edited, :rate, :w]]; fill(NaN, ncells)]))
end

wjob = BasicMCJob(
  model,
  RAM(),
  BasicMCRange(nsteps=100, burnin=10),
  init;
  pindex=6,
  resetpstate=false,
  verbose=false
)

mcrange = BasicMCRange(nsteps=1100, burnin=100)

outopts = fill(Dict(:monitor => [:value]), ncells+1)

gibbsjob = GibbsJob(model, Dict(:w => wjob), mcrange, init)

@time run(gibbsjob)

output(gibbsjob)

chains = Dict(gibbsjob)

v = map(x -> init[:μ]*(1-init[:μ])*inv_logit(x), chains[:w].value)

for i in 1:ncells
  writecsv(string("pchain", i, ".csv"), chains[Symbol('p', i)].value)
end
writecsv("vchain.csv", v)


text:
@testset "Krylov K3" begin
  for kmethod in [:bicgstab, :usymqr, :qmr]
    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(
        sp = K3KrylovParams(kmethod = kmethod),
        solve_method = :IPF,
        history = true,
      ),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K3KrylovParams(uplo = :U, kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-4, ϵ_rb = 1.0e-4, ϵ_pdd = 1.0e-4),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "KrylovK3_5" begin
  for kmethod in [:minres, :minres_qlp]
    stats1 = ripqp(
      QuadraticModel(qps1),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(kmethod = kmethod), history = true),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats1.objective, -1.59078179, atol = 1e-1)
    @test stats1.status == :acceptable

    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(uplo = :U, kmethod = kmethod), solve_method = :IPF),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K3_5KrylovParams(kmethod = kmethod)),
      itol = InputTol(max_iter = 50, max_time = 20.0, ϵ_rc = 1.0e-2, ϵ_rb = 1.0e-2, ϵ_pdd = 1.0e-2),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end

@testset "K3_5 structured" begin
  for kmethod in [:tricg, :trimr]
    stats1 = ripqp(
      QuadraticModel(qps1),
      display = false,
      iconf = InputConfig(
        sp = K3_5StructuredParams(kmethod = kmethod),
        solve_method = :IPF,
        history = true,
      ),
      itol = InputTol(
        max_iter = 100,
        max_time = 20.0,
        ϵ_rc = 1.0e-2,
        ϵ_rb = 1.0e-2,
        ϵ_pdd = 1.0e-2,
      ),
    )
    @test isapprox(stats1.objective, -1.59078179, atol = 1e-1)
    @test stats1.status == :acceptable

    stats2 = ripqp(
      QuadraticModel(qps2),
      display = false,
      iconf = InputConfig(sp = K3_5StructuredParams(kmethod = kmethod), solve_method = :IPF),
      itol = InputTol(
        max_iter = 100,
        max_time = 20.0,
        ϵ_rc = 1.0e-2,
        ϵ_rb = 1.0e-2,
        ϵ_pdd = 1.0e-2,
      ),
    )
    @test isapprox(stats2.objective, -9.99599999e1, atol = 1e-1)
    @test stats2.status == :acceptable

    stats3 = ripqp(
      QuadraticModel(qps3),
      display = false,
      iconf = InputConfig(sp = K3_5StructuredParams(kmethod = kmethod), solve_method = :IPF),
      itol = InputTol(
        max_iter = 100,
        max_time = 20.0,
        ϵ_rc = 1.0e-2,
        ϵ_rb = 1.0e-2,
        ϵ_pdd = 1.0e-2,
      ),
    )
    @test isapprox(stats3.objective, 5.32664756, atol = 1e-1)
    @test stats3.status == :acceptable
  end
end


text:
using Test, Documenter
using IndividualDisplacements, MeshArrays, OrdinaryDiffEq, OceanStateEstimation
include(joinpath(dirname(pathof(IndividualDisplacements)),"../examples/helper_functions.jl"))

IndividualDisplacements.flt_example_download()
OceanStateEstimation.get_ecco_velocity_if_needed()
OceanStateEstimation.get_occa_velocity_if_needed()

@testset "test1" begin
    uvetc,sol=test1_setup()
    @test isapprox(sol[1,end],23.4474; atol=0.01)
    @test isapprox(sol[2,end],8.0896; atol=0.01)
end

@testset "test2" begin
    df,𝑃=test2_periodic_domain()
    @test prod(isapprox.(df[end-35:6:end,:y],12*(0.4:0.04:0.6)))
    @test prod(isapprox.(df[end-5:end,:x],12*(0.4:0.04:0.6).+4.0))
end

@testset "test3" begin
    p=dirname(pathof(IndividualDisplacements))
    include(joinpath(p,"../examples/basics/random_flow_field.jl"))

    show(𝐼)
    diff(𝐼)
    size(𝐼)
    𝐽=similar(𝐼)
    @test isa(𝐽,Individuals)

    tmp1=randn_lonlat(10)
    𝐺=convert_to_FlowFields(u,v,10.0)
    tmp2=nearest_to_xy(𝐺.u0,3.,3.,1.)
    @test isa(tmp2,Array)
    tmp3=nearest_to_xy(𝐹.u0,3.,3.)
    @test isa(tmp3,Array)
end

@testset "doctests" begin
    doctest(IndividualDisplacements; manual = false)
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    numbersCalled = []
    boards = []

    open(fname, "r") do f
        numbersCalled = map(x -> parse(Int, x), split(readline(f), ","))

        # skip the empty line
        readline(f)

        while !eof(f)
            currentBoard = zeros(Int, 5,5)

            for i=1:5
                currentBoard[i,:] = map(x -> parse(Int, x), split(readline(f), " ", keepempty=false))
            end

            push!(boards, currentBoard)

            readline(f)
        end
    end

    return numbersCalled, boards

end

function WinningBoard(board::Matrix)
    # a board has won when any row/column sums to -5 (after matching entries are set to -1)
    for i = 1:5
        if sum(board[i,:]) == -5 || sum(board[:,i]) == -5
            return true
        end
    end
    return false
end

function UpdateBoard!(board, numberCalled)
    # when a number is called, set any matching board elements to -1
    ndx = findfirst(x -> x==numberCalled, board)
    if ndx != nothing
        setindex!(board, -1, ndx)
    end
end

function Part1(data)
    numbersCalled, boards = data

    winner = nothing
    ndx = 0
    while winner == nothing
        ndx += 1
        map(b -> UpdateBoard!(b, numbersCalled[ndx]), boards)
        winner = findfirst(WinningBoard, boards)
    end

    winningBoard = boards[winner]
    return (sum(winningBoard) + count(x -> x==-1, winningBoard)) * numbersCalled[ndx]
end

function Part2(data)
    numbersCalled, boards = data

    lastWinner = zeros(Int, 5,5)
    lastNumber = -1
    ndx = 0
    while ndx < length(numbersCalled)
        winner = nothing
        while winner == nothing && ndx < length(numbersCalled)
            ndx += 1
            map(b -> UpdateBoard!(b, numbersCalled[ndx]), boards)
            winner = findfirst(WinningBoard, boards)
        end

        if winner != nothing
            lastWinner = boards[winner]
            lastNumber = numbersCalled[ndx]
            filter!(x -> !WinningBoard(x), boards)
        end
    end

    return (sum(lastWinner) + count(x -> x==-1, lastWinner)) * lastNumber
end

data = ReadData("../data/test4.txt")
data = ReadData("../data/day4.txt")

println(Part1(data))

println(Part2(data))


text:
export LoadCamvid, GetTVvalues

"""
Obtain the anisotropic total-variation per image. If the noise factor is larger than zero,
the function return a number from the uniform distribution [trueTV * (1-noise_factor) - trueTV * (1+noise_factor)]
"""
function GetTVvalues(grid,labels,noise_factor,compgrid)

  (TV,dummy1,dummy2,dummy3) = get_TD_operator(compgrid,"TV",Float32)

  true_TV_list = zeros(Float32,length(labels))

  for i=1:length(labels)
    true_TV_list[i] = norm(TV*vec(labels[i][:,:,1]),1)
  end

    #add 'noise' to true tv list
  if noise_factor>0.0
    for i=1:length(labels)
      d = Uniform((1-noise_factor)*true_TV_list[i],(1+noise_factor)*true_TV_list[i])
      true_TV_list[i] = rand(d,1)[1]
    end
  end

  return true_TV_list
end

function LoadCamvid(basepath::String,nr_output_chan::Int)

  current_dir = pwd()

  #set data directories
  data_dir_train  = joinpath(basepath,"Camvid/train/")
  data_dir_val    = joinpath(basepath,"Camvid/val/")
  label_dir_train = joinpath(basepath,"Camvid/trainannot/")
  label_dir_val   = joinpath(basepath,"Camvid/valannot/")

  #read file names
  all_train_data = readdir(data_dir_train);
  all_val_data   = readdir(data_dir_val);

  all_train_labels = readdir(label_dir_train);
  all_val_labels   = readdir(label_dir_val);

  #allocate tensors for data/label storage
  train_data     = zeros(Float32,3,360,480,62);
  val_data       = zeros(Float32,3,360,480,10);

  train_labels   = zeros(Float32,360,480,1,62);
  val_labels     = zeros(Float32,360,480,1,10);

  #read training data (hardcoded for getting the first video only)
  cd(data_dir_train)
  for i=1:62#length(all_train_data)
      temp_img = Images.load(all_train_data[i])
      train_data[1:3,:,:,i] = Images.channelview(temp_img)
      println(i)
   end

  #permute to x-y-n_chan-n_examples
  train_data    = permutedims(train_data,[2,3,1,4]);

  #repeat data to increase channel count explicitly for input-output
  repeat_factor = Int(nr_output_chan/3)
  train_data    = repeat(train_data,outer=[1,1,repeat_factor,1])

  #read val data
  cd(data_dir_val)
  for i=1:10#length(all_val_data)
      temp_img = Images.load(all_val_data[i])
      val_data[1:3,:,:,i] = Images.channelview(temp_img)
      println(i)
   end
  val_data = permutedims(val_data,[2,3,1,4]);;
  val_data=repeat(val_data,outer=[1,1,16,1])

  #read training labels
  cd(label_dir_train)
  counter = 1
  for i=1:62#length(all_train_labels)
      temp_img = Images.load(all_train_labels[i])
      train_labels[:,:,1,i] = Images.channelview(temp_img)
      println(i)
   end

   #read validation labels
   cd(label_dir_val)
   counter = 1
   for i=1:10#length(all_val_labels)
       temp_img = Images.load(all_val_labels[i])
       val_labels[:,:,1,i] = Images.channelview(temp_img)
       println(i)
    end

    #reduce resolution and make sure each dimension can be divided by 2 at least 2 times (to be able to take wavelet transforms)
    train_labels = train_labels[1:2:end,1:2:end,:,:]; train_labels = train_labels[1:176,:,:,:];
    val_labels   = val_labels[1:2:end,1:2:end,:,:];   val_labels   = val_labels[1:176,:,:,:];
    val_data     = val_data[1:2:end,1:2:end,:,:];     val_data     = val_data[1:176,:,:,:];
    train_data   = train_data[1:2:end,1:2:end,:,:];   train_data   = train_data[1:176,:,:,:];

    #change labels to 1-hot encodings
    n_class             = length(unique(train_labels))
    class_values        = sort(unique(train_labels))
    train_labels_OneHot = zeros(Int,size(train_labels,1),size(train_labels,2),length(class_values),size(train_labels,4));
    for i=1:length(class_values)
        println(i)
        class_indices = findall(train_labels[:,:,1,:] .== class_values[i])
        class_image = zeros(Int,size(train_labels)[1],size(train_labels)[2],size(train_labels)[4])
        class_image[class_indices] .= 1
        train_labels_OneHot[:,:,i,:] .= class_image
    end

    val_labels_OneHot = zeros(Int,size(val_labels,1),size(val_labels,2),length(class_values),size(val_labels,4));
    for i=1:length(class_values)
        class_indices = findall(val_labels[:,:,1,:] .== class_values[i])
        class_image = zeros(Int,size(val_labels)[1],size(val_labels)[2],size(val_labels)[4])
        class_image[class_indices] .= 1
        val_labels_OneHot[:,:,i,:] .= class_image
    end

    #print % of pixels in each class for the first video
    println("% of pixels in each class, for first video in CamVid dataset")
    for i=1:12;
      println("class: ",i)
      println(sum(train_labels_OneHot[:,:,i,:])/prod(size(train_labels_OneHot[:,:,1,:])).*100);
    end

    # selected_classes = [1, 2, 4, 5, 6, 9, 12]
    # combine_classes = [3,7,8,10,11]
    #
    # val_labels_OneHot[:,:,[combine_classes[1]],:]  .= sum(val_labels_OneHot[:,:,combine_classes,:],dims=3)
    # train_labels_OneHot[:,:,[combine_classes[1]],:].= sum(train_labels_OneHot[:,:,combine_classes,:],dims=3)

    #class 9 corresponds to vehicles
    selected_classes=[9]
    #val_labels_OneHot = val_labels_OneHot[:,:,selected_classes,:]
    train_labels_OneHot = train_labels_OneHot[:,:,selected_classes,:]

    #add the second channel (1 - first channel) (this is unnecessary if n_classes > 3)
    train_labels_OneHot = repeat(train_labels_OneHot,outer=[1,1,2,1])
    train_labels_OneHot[:,:,2,:] .= 1 .- train_labels_OneHot[:,:,1,:]

    #create a vector with one label per entry
    train_labels   = Vector{Any}(undef,size(train_data)[end])
    val_labels     = Vector{Any}(undef,size(val_data)[end])

    for i=1:length(train_labels)
      train_labels[i] = train_labels_OneHot[:,:,:,i]
    #   train_labels[i] = zeros(Float32,size(train_labels_OneHot,1),size(train_labels_OneHot,2),2,1)
    #   train_labels[i][:,:,1,1] .= 1f0 .*train_labels_OneHot[:,:,9,i]
    #   train_labels[i][:,:,2,1] .= 1f0 .- train_labels_OneHot[:,:,9,i]
    end

    for i=1:length(val_labels)
      val_labels[i] = val_labels_OneHot[:,:,:,i]
    #   val_labels[i] = zeros(Float32,size(val_labels_OneHot,1),size(val_labels_OneHot,2),2,1)
    #   val_labels[i][:,:,1,1] .= val_labels_OneHot[:,:,9,i]
    #   val_labels[i][:,:,2,1] .= 1f0 .- val_labels_OneHot[:,:,9,i]
    end

    #change data to a vector of examples
    dataL    = Vector{Array{Float32,4}}(undef,size(train_data,4))
    datavalL = Vector{Array{Float32,4}}(undef,size(val_data,4))
    for i=1:length(dataL)
        dataL[i] = train_data[:,:,:,[i]]
    end
    for i=1:length(datavalL)
        datavalL[i] = val_data[:,:,:,[i]]
    end
    train_data = deepcopy(dataL)
    dataL      = []
    val_data   = deepcopy(datavalL)
    datavalL   = []

    cd(current_dir) #switch back to current working dir

    #repartition training and validation:
    val_data   = deepcopy(train_data[47:end])
    train_data = train_data[1:46]

    val_labels = deepcopy(train_labels[47:end])
    train_labels = train_labels[1:46]

    return train_data, train_labels, val_data, val_labels
end #end function


text:
module Runtests

using Test, DescriptorSystems

@testset "Test DescriptorSystems" begin
# test constructors
include("test_dss.jl")
include("test_rtf.jl")
include("test_polrat.jl")
include("test_connections.jl")
# test basic functions
include("test_conversions.jl")
include("test_ordred.jl")
include("test_analysis.jl")
# test factorizations
include("test_pscf.jl")
include("test_cfid.jl")
include("test_iofac.jl")
# test advanced operations
include("test_gsdec.jl")
include("test_nullrange.jl")
include("test_covers.jl")
include("test_linsol.jl")
include("test_ginv.jl")
include("test_model_matching.jl")
end

end


text:
# This file was generated, do not modify it. # hide
boston[3, 5]


text:
using ITensors
using LinearAlgebra
using Printf
using Random
using TBLIS

Random.seed!(1234)

let
  nthreads = 4

  N = 100
  sites = siteinds("S=1",N)
  ampo = AutoMPO()
  for j=1:N-1
    ampo += 0.5, "S+", j, "S-", j+1
    ampo += 0.5, "S-", j, "S+", j+1
    ampo +=      "Sz", j, "Sz", j+1
  end
  H = MPO(ampo,sites)
  psi0 = randomMPS(sites,10)

  sweeps_compile = Sweeps(2)
  maxdim!(sweeps_compile, 10)
  cutoff!(sweeps_compile, 1E-15)

  sweeps = Sweeps(6)
  maxdim!(sweeps, 20,100,200,300)
  cutoff!(sweeps, 0.0)
  @show sweeps

  #
  # Using BLAS backend
  #

  disable_tblis!()
  BLAS.set_num_threads(nthreads)

  # Compile
  dmrg(H, psi0, sweeps_compile; outputlevel = 0)

  println("Using BLAS with $nthreads threads\n")
  energy, psi = dmrg(H, psi0, sweeps)
  @printf("Final energy = %.12f\n",energy)
  println()

  #
  # Using TBLIS backend
  #

  enable_tblis!()
  BLAS.set_num_threads(1)
  TBLIS.set_num_threads(nthreads)

  # Compile
  dmrg(H, psi0, sweeps_compile; outputlevel = 0)

  println("Using TBLIS with $(TBLIS.get_num_threads()) threads (and 1 BLAS thread)\n")
  energy, psi = dmrg(H, psi0, sweeps)
  @printf("Final energy = %.12f\n",energy)
end


text:
module NMRTools

using Colors
using DimensionalData
using LsqFit
using Measurements
using OffsetArrays
using Optim
using Plots
using SimpleTraits
using SpecialFunctions
using Statistics

import Base.getindex, Base.setindex!
import DimensionalData.dims, DimensionalData.refdims, DimensionalData.data, DimensionalData.name, DimensionalData.metadata, DimensionalData.label, DimensionalData.rebuild
import DimensionalData.X, DimensionalData.Y, DimensionalData.Z, DimensionalData.Ti
import DimensionalData.Between, DimensionalData.At, DimensionalData.Near

export NMRData
export dims, X, Y, Z, Ti
export xval, yval, zval, tval, settval
export At, Near, Between
export haspseudodimension, HasPseudoDimension
export NMRToolsException
export loadnmr
export data, acqus, label, label!, metadata, metadatahelp, scale
export WindowFunction, NullWindow, UnknownWindow, ExponentialWindow, SineWindow, GaussWindow
export GeneralSineWindow, CosWindow, Cos²Window, GeneralGaussWindow, LorentzToGaussWindow
export estimatenoise!
export fitexp, fitdiffusion

# type definitions
include("util.jl")
include("exceptions.jl")
include("nmrdata.jl")
include("windows.jl")

# routines
include("noise.jl")
include("loadnmr.jl")
include("plotrecipes.jl")
include("fittingpseudo2d.jl")

end # module


text: println((sum(1:100) ^ 2) - sum(map(x -> x ^ 2, 1:100)))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
function initialprobs(
    lower::Array{T,1},
    upper::Array{T,1};
    maxsamples = 10000,
) where {T<:Number}
    p = length(lower)
    @assert p == length(upper)

    bitlen = p * 32
    mat = Array{Int8,2}(undef, maxsamples, bitlen)
    for tries = 1:maxsamples
        randvalues = map((x, y) -> x + rand() * (y - x), lower, upper)
        mat[tries, :] .= bits(randvalues)
    end
    return colmeans(mat)
end

function sample(probs::Array{T,1}) where {T<:Number}
    return map(x -> if rand() < x
        1
    else
        0
    end, probs)
end

function mccga(;
    lower::Array{T,1},
    upper::Array{T,1},
    costfunction::Function,
    popsize::Int,
    maxsamples = 10000,
) where {T<:Number}
    probvector = initialprobs(lower, upper, maxsamples = maxsamples)
    initial_probvector = copy(probvector)
    chsize = length(probvector)
    mutation = 1.0 / convert(Float64, popsize)
    numiters = 0
    while !(all(x -> (x <= mutation) || (x >= 1.0 - mutation), probvector))
        numiters += 1
        ch1 = sample(probvector)
        ch2 = sample(probvector)
        cost1 = costfunction(floats(ch1))
        cost2 = costfunction(floats(ch2))
        winner = ch1
        loser = ch2
        if (cost2 < cost1)
            winner = ch2
            loser = ch1
        end
        for i = 1:chsize
            if winner[i] != loser[i]
                if winner[i] == 1
                    probvector[i] += mutation
                else
                    probvector[i] -= mutation
                end
                if probvector[i] > 1
                    probvector[i] = 1
                elseif probvector[i] < 0
                    probvector[i] = 0
                end
            end
        end
    end

    sampledvector = sample(probvector)
    initial_solution = floats(sampledvector)

    local_result = Optim.optimize(costfunction, initial_solution, Optim.NelderMead())

    resultdict = Dict(
        "initial_prob_vector" => initial_probvector,
        "final_prob_vector" => sampledvector,
        "iterations" => numiters,
        "initial_solution" => floats(sampledvector),
        "final_solution" => local_result.minimizer,
        "initial_minimum" => costfunction(floats(sampledvector)),
        "final_minimum" => costfunction(local_result.minimizer),
    )

    return resultdict
end

max_stars_repo_name: jbytecode/MCCGA

max_issues_repo_name: jbytecode/MCCGA

max_forks_repo_name: jbytecode/MCCGA


text:
using Images, FileIO
using ImageTransformations
using CairoMakie
using JLD2
using WindMixing: plot_T_profiles_final_diurnal

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot\\z.png")))
T_3D_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T_3D.png"))))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
    T_3D = T_3D_img,
)

RESULTS_DIR = "C:\\Users\\xinle\\Documents\\OceanParameterizations.jl\\final_results"
NDE_DIR = "3sim_diurnal"

function open_NDE_profile(filepath)
    file = jldopen(filepath)
    NDE_profile = file["NDE_profile"]
    close(file)
    return NDE_profile
end

files_training = [
    "train_wind_-5e-4_diurnal_5e-8"
    "train_wind_-3.5e-4_diurnal_5e-8"
    "train_wind_-2e-4_diurnal_5e-8"
]

NDE_profiles_training = [
    open_NDE_profile(joinpath(RESULTS_DIR, NDE_DIR, file, "profiles_fluxes_oceananigans.jld2")) for file in files_training
]

subplot_titles_training = [
    "Strong Wind, Strong Cycle"
    "Moderate Wind, Strong Cycle"
    "Weak Wind, Strong Cycle"
]

frame = 1009

plot_T_profiles_final_diurnal(NDE_profiles_training, frame, subplot_titles_training, axis_images, "final_results\\diurnal_training_T.png")

function calculate_losses(NDE_profiles)
    @inline T_loss(data) = sum(data) / 1153

    return [
        (
            NN = T_loss(profile["T_losses"]),
            mpp = T_loss(profile["T_losses_modified_pacanowski_philander"]),
            kpp = T_loss(profile["T_losses_kpp"]),
        ) for profile in NDE_profiles
    ]
end

calculate_losses(NDE_profiles_training)

files_interpolating = [
    "test_wind_-4.5e-4_diurnal_4e-8"
    "test_wind_-3e-4_diurnal_4e-8"
    "test_wind_-3e-4_diurnal_2e-8"
]

NDE_profiles_interpolating = [
    open_NDE_profile(joinpath(RESULTS_DIR, NDE_DIR, file, "profiles_fluxes_oceananigans.jld2")) for file in files_interpolating
]

subplot_titles_interpolating = [
    "Strong Wind, Strong Cycle"
    "Weak Wind, Strong Cycle"
    "Weak Wind, Weak Cycle"
]

frame = 1009

plot_T_profiles_final_diurnal(NDE_profiles_interpolating, frame, subplot_titles_interpolating, axis_images, "final_results\\diurnal_interpolating_T.png")

calculate_losses(NDE_profiles_interpolating)

files_extrapolating = [
    "test_wind_-5.5e-4_diurnal_5.5e-8"
    "test_wind_-1.5e-4_cooling_3.5e-8"
    "test_wind_-5.5e-4_new"

    "test_wind_-1.5e-4_diurnal_5.5e-8"
    "test_wind_-1.5e-4_heating_-3.5e-8"
    "test_wind_-5.5e-4_cooling_3.5e-8"
]

NDE_profiles_extrapolating = [
    open_NDE_profile(joinpath(RESULTS_DIR, NDE_DIR, file, "profiles_fluxes_oceananigans.jld2")) for file in files_extrapolating
]

subplot_titles_extrapolating = [
    "Strong Wind, Strong Cycle"
    "Weak Wind, Strong Cooling"
    "Strong Wind Mixing"

    "Weak Wind, Strong Cycle"
    "Weak Wind, Strong Heating"
    "Strong Wind, Strong Cooling"
]

frame = 1009

plot_T_profiles_final(NDE_profiles_extrapolating, frame, subplot_titles_extrapolating, axis_images, "final_results\\diurnal_extrapolating_T.png")

calculate_losses(NDE_profiles_extrapolating)


text:
module GUI

using ..CImGui
using ..CImGui.CSyntax
using ..CImGui.CSyntax.CStatic
using ..CImGui.GLFWBackend
using ..CImGui.OpenGLBackend
using ..CImGui.GLFWBackend.GLFW
using ..CImGui.OpenGLBackend.ModernGL
using ..Gamepad, ..Zeiss, ..Pressure, ..Manipulator, ..DAQ #..Camera
using ImPlot

import ThreadPools.@tspawnat
import CImGui: ImVec2, ImVec4, ImTextureID
import ..ZEISS_TID

# Globals
global IMGUI_WINDOW #TODO: as const ref
global IMGUI_CONTEXT #TODO: as const ref
const GLSL_VERSION = 130
const BACKGROUND_COLOR = Cfloat[0.45, 0.55, 0.60, 1.00]
const PRESSURE_TID = 2
const FRAME_STEP = Threads.Condition()
const SPEED_STEP_MAX = 6

mutable struct GlobalState
    SPEED_STEP::Int64
    USE_GAMEPAD::Ref{Bool}
    GAMEPAD::GamepadState
    GAMEPAD_CACHE::GamepadState
    function GlobalState()
        speedstep = 1
        usegamepad = Ref{Bool}(false)
        gamepad = Gamepad.poll()
        return new(speedstep, usegamepad, gamepad, gamepad)
    end
end

const g = Ref{GlobalState}()
# GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"

function step_notify(c::Threads.Condition)
    lock(c)
    try
        notify(c)
    finally
        unlock(c)
    end
    return
end

function inc_speed()
    g[].SPEED_STEP == SPEED_STEP_MAX && return
    g[].SPEED_STEP += 1
end

function dec_speed()
    g[].SPEED_STEP == 1 && return
    g[].SPEED_STEP -= 1
end

function init_glfw()
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 4)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 6)
    GLFW.SetErrorCallback(error_callback)
end

function setup()
    global IMGUI_CONTEXT = CImGui.CreateContext()
    CImGui.StyleColorsDark()
    fonts_dir = joinpath(@__DIR__, "..", "fonts")
    fonts = CImGui.GetIO().Fonts
    CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "MesloLGS.ttf"), 20)
    global IMGUI_WINDOW = GLFW.CreateWindow(1280, 720, "LabRig")
    @assert IMGUI_WINDOW != C_NULL
    GLFW.MakeContextCurrent(IMGUI_WINDOW)
    GLFW.SwapInterval(1)  # enable vsync; 0 for disabled, 2 for 30fps, 3 for 20fps, etc
    ImGui_ImplGlfw_InitForOpenGL(IMGUI_WINDOW, true)
    ImGui_ImplOpenGL3_Init(GLSL_VERSION)
end

function new_frame!()
    ImGui_ImplOpenGL3_NewFrame()
    ImGui_ImplGlfw_NewFrame()
    CImGui.NewFrame()
end

function render!(window = IMGUI_WINDOW, clear_color = BACKGROUND_COLOR)
    CImGui.Render()
    GLFW.MakeContextCurrent(window)
    display_w, display_h = GLFW.GetFramebufferSize(window)
    glViewport(0, 0, display_w, display_h)
    glClearColor(clear_color...)
    glClear(GL_COLOR_BUFFER_BIT)
    ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())
    GLFW.MakeContextCurrent(window)
    GLFW.SwapBuffers(window)
end

function shutdown!(; window = IMGUI_WINDOW, ctx = IMGUI_CONTEXT)
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.DestroyWindow(window)
end
#=
# do this later...
function menubar()
    if CImGui.BeginMainMenuBar()
        if CImGui.BeginMenu("File")
=#

abstract type UIElement end
#=
mutable struct CameraUI <: UIElement
    frame::Vector{UInt16}
    pipeline::ImagePipeline
    pipeline_job::Task
    texID::UInt32
    livemode::Bool
    function CameraUI()
        livemode = false
        frame = zeros(UInt16, 1200*1200)
        pipeline = ImagePipeline(frame)
        pipeline_job = @async 0
        texID = Camera.gen_textures(1)[1]
        Camera.load_texture(texID, frame, 1200, 1200)
        camui = new(frame, pipeline, pipeline_job, texID, livemode)
        finalizer(camui) do camui
            camui.livemode && Camera.stopcont()
        end
        return camui
    end
end

function (c::CameraUI)()
    CImGui.Begin("Live View")
    if CImGui.Button("Start Live")
        try
            Camera.polled_cont!()
            c.pipeline_job = @async sleep(1.0)
            c.livemode = true
        catch e
            @warn "Failed to initialize live mode"
            c.livemode = false
        end
    end
    CImGui.SameLine()
    if CImGui.Button("Stop Live")
        Camera.stop_cont()
        c.livemode = false
    end
    if c.livemode
        Camera.reload_texture(c.texID, c.pipeline.out, 1200, 1200)
        CImGui.Image(ImTextureID(Int(c.texID)), ImVec2(1200,1200), ImVec2(0,0), ImVec2(1,1),
                     ImVec4(1.0,1.0,1.0,1.0), ImVec4(1,1,1,1))
        if istaskdone(c.pipeline_job)
            try
                c.frame .= Camera.latest_frame()
                c.pipeline_job = @tspawnat 6 c.pipeline(c.frame)
            catch e
                @warn "Frame dropped"
            end
        end
    end
    CImGui.End()
end
=#
mutable struct PressureUI <: UIElement
    cmd_pressure::Ref{Cint}
    show_history::Bool
    fs::Int64
    coax_amp::Ref{Int32}
    coax_duration::Ref{Float32}
    coax_repeat::Ref{Int32}
    rupture_amp::Ref{Int32}
    rupture_duration::Ref{Float32}
    rupture_repeat::Ref{Int32}
    function PressureUI()
        cmd_pressure = Ref{Cint}(0)
        @tspawnat PRESSURE_TID Pressure.stream_out(FRAME_STEP)
        new(cmd_pressure, false, 60, Ref{Int32}(5), Ref{Float32}(3.0), Ref{Int32}(3), Ref{Int32}(-30), Ref{Float32}(1.0), Ref{Int32}(3))
    end
end

function (p::PressureUI)()
    CImGui.Begin("Pressure control.")
    CImGui.Button("Perfusion Channel: " * (PRIMARY_SOLUTION[] ? "Main" : "Drug")) && toggle_solution()
    CImGui.Text("Current Presure: $(Pressure.CURRENT_PRESSURE[])")
    CImGui.DragInt("kPa", p.cmd_pressure, 0.5, -20, 100)
    CImGui.Button("On Cell!") && (p.cmd_pressure[] = Cint(-1))
    CImGui.SameLine()
    CImGui.Button("NEUTRAL") && (p.cmd_pressure[] = Cint(0))
    CImGui.SameLine()
    CImGui.Button("Near cell") && (p.cmd_pressure[] = Cint(1))
    CImGui.SameLine()
    CImGui.Button("Bath mode") && (p.cmd_pressure[] = Cint(20))

    #build an interface for sending ramps/transitions
    CImGui.InputInt("Coax Amp", p.coax_amp)
    CImGui.InputFloat("Coax Duration", p.coax_duration)
    CImGui.InputInt("Coax repeat", p.coax_repeat)
    CImGui.Button("Coax Cell") && Pressure.coaxcell(p.coax_amp[], p.fs, p.coax_duration[], p.coax_repeat[])

    CImGui.InputInt("Rupture Amp", p.rupture_amp)
    CImGui.InputFloat("Rupture Duration", p.rupture_duration)
    CImGui.InputInt("Rupture repeat", p.rupture_repeat)
    CImGui.Button("Rupture cell") && Pressure.attempt_break(p.rupture_amp[], p.fs, p.rupture_duration[], p.rupture_repeat[])

    if isempty(Pressure.STREAM_CHANNEL)
        put!(Pressure.STREAM_CHANNEL, [p.cmd_pressure[]])
    end
    CImGui.Button("Toggle History") && (p.show_history = !p.show_history)
    if p.show_history
        if ImPlot.BeginPlot("Pressure Sensor", "", "", ImVec2(-1, 300))
            ImPlot.PlotLine(collect(Pressure.PROBE_HISTORY), label = "kPa")
            ImPlot.EndPlot()
        end
    end
    CImGui.End()
end

mutable struct FocusUI <: UIElement
    home_pos::Int64
    work_pos::Int64
    function FocusUI()
        new(0,0)
    end
end

function follow_work(focus_zpos, fui::FocusUI)
    ump_position = Manipulator.get_position()
    diff = Cint(fui.work_pos - focus_zpos)
    target_position = Manipulator.Position(ump_position.x, ump_position.y, ump_position.z + diff, nothing)
    Manipulator.moveto(target_position)
end

function follow_home(focus_zpos, fui::FocusUI)
    ump_position = Manipulator.get_position()
    diff = Cint(fui.home_pos - focus_zpos)
    target_position = Manipulator.Position(ump_position.x, ump_position.y, ump_position.z + diff, nothing)
    Manipulator.moveto(target_position)
end

function (f::FocusUI)()
    position, velocity = fetch(@tspawnat ZEISS_TID Zeiss.get_zeiss_state())
    gamepad = g[].GAMEPAD
    gamepad_cache = g[].GAMEPAD_CACHE
    CImGui.Begin("Focus control")
    CImGui.Checkbox("Use gamepad", g[].USE_GAMEPAD)
    CImGui.Text("Position: $position")
    CImGui.Button("Go Home") && Zeiss.moveto(f.home_pos)
    CImGui.SameLine()
    CImGui.Button("Set Home: $(f.home_pos)") && (f.home_pos = position)
    CImGui.Button("Go Work") && Zeiss.moveto(f.work_pos)
    CImGui.SameLine()
    CImGui.Button("Set Work: $(f.work_pos)") && (f.work_pos = position)
    CImGui.Button("Follow Work") && follow_work(position, f)
    # Follow home is disabled atm. You would need to do it by caching Z-axis positions. The current logic fails
    # because you need to move the focus first (which alters its position) and changes the calculated offset.
    # and providing a done event that won't cause frame lag (e.g. wait for move to finish for X frames; track progress by incrementing a counter on each frame after a move)
    # In general, a return to a position where the manipulator can be retracted for pipette change would be more useful than return to home focus position!!
    # CImGui.SameLine()
    # CImGui.Button("Follow Home") && follow_home(position, f)
    CImGui.Button("STOPPP!!!!!") && Zeiss.stop!()
    if g[].USE_GAMEPAD[] # need global state for gamepad
        CImGui.Text("Z Axis Velocity: $velocity")
        CImGui.Text("Speed step: $(g[].SPEED_STEP)")
        gamepad.DPAD.left && !gamepad_cache.DPAD.left && dec_speed()
        gamepad.DPAD.right && !gamepad_cache.DPAD.right && inc_speed()
        gamepad.button.B ? Zeiss.stop!() : new_velocity = Zeiss.update_velocity(gamepad, g[].SPEED_STEP)
        Zeiss.set_velocity(new_velocity)
        Manipulator.step!(gamepad)
    end
    CImGui.End()
end

mutable struct DAQUI <: UIElement
    recording::DAQ.Recording
    data::Vector{Float64}
    xii_indexes::StepRange
    vm_indexes::StepRange
    record::Bool
    function DAQUI()
        recording = DAQ.Recording()
        xii_indexes = 1:4:length(recording.signal)
        vm_indexes = 2:4:length(recording.signal)
        data = zeros(length(recording.signal))
        return new(recording, data, xii_indexes, vm_indexes, false)
    end
end

function (d::DAQUI)()
    CImGui.Begin("Ephys Data")
    if !d.record
        if CImGui.Button("Start acquisition.")
            wait(d.recording())
            d.record = true
        end
    else
        if CImGui.Button("Stop acquisition.")
            DAQ.stop(d.recording.task)
            d.record = false
        else
            d.data .= collect(d.recording.signal)

            if ImPlot.BeginPlot("Vm Data", "","", ImVec2(-1,700))
                # downsampling to 2kHz display
                ImPlot.PlotLine(d.data[d.vm_indexes], label = "Vm")
                ImPlot.EndPlot()
            end

            if ImPlot.BeginPlot("XII Data", "", "", ImVec2(-1,500))
                #downsampling to 2kHz display
                ImPlot.PlotLine(d.data[d.xii_indexes], label = "XII")
                ImPlot.EndPlot()
            end
        end
    end
    CImGui.End()
end

mutable struct StimUI <: UIElement
    stim::DAQ.Stimulus
    user_interval::Ref{Float32}
    user_duration::Ref{Float32}
    function StimUI()
        stim = DAQ.Stimulus()
        return new(stim, Ref{Float32}(0.05), Ref{Float32}(0.03))
    end
end

function (s::StimUI)()
    CImGui.Begin("Stimulus Panel")

    CImGui.Button("Patch pulses.") && DAQ.pulsetiming!(s.stim, 0.05, 0.03)
    CImGui.SameLine()
    CImGui.Button("Ia/Ih pulses.") && DAQ.pulsetiming!(s.stim, 1.0, 0.4)
    CImGui.SameLine()
    CImGui.Button("Custom pulses.") && DAQ.pulsetiming!(s.stim, Float64(s.user_interval[]), Float64(s.user_duration[]))
    CImGui.SameLine()
    CImGui.Button("Stop pulses.") && DAQ.stop(s.stim)

    CImGui.InputFloat("Custom interval", s.user_interval)
    CImGui.InputFloat("Custom duration", s.user_duration)

    CImGui.End()
end

mutable struct ElementState
    show::Bool
    hotload::Bool
    function ElementState()
        return new(false, false)
    end
end

function build_ui()
    ui = IdDict{UIElement, ElementState}()
    for x in subtypes(UIElement)
        ui[x()] = ElementState()
    end
    return ui
end

function run_loop(; window = IMGUI_WINDOW)

    @async try
        g[] = GlobalState()
        ui = UIElement[]
        for x in [FocusUI, PressureUI, DAQUI, StimUI] #FIXME: add CameraUI back after adjusting hist
            push!(ui, x())
        end
        while !GLFW.WindowShouldClose(window)
            GLFW.PollEvents()
            new_frame!()
            g[].GAMEPAD = Gamepad.poll()

            # run UI functions...
            for el in ui
                Base.invokelatest(el)
            end
            # cache gamepad state
            g[].GAMEPAD_CACHE = g[].GAMEPAD
            # trigger pipelines on other threads
            step_notify(FRAME_STEP)
            render!()
            yield()
        end # main loop
    catch e
        @error "Error in renderloop!" exception=e
        Base.show_backtrace(stderr, catch_backtrace())
    finally
        shutdown!()
    end # try-catch-finally; @async; GC.@preserve
end # function

function launch!()
    setup()
    run_loop()
end

end # module


text:
import Base.convert

"""
    convert(::Type{UnipartiteNetwork}, N::T) where {T <: BipartiteNetwork}

Projects a deterministic bipartite network in its unipartite representation.
"""
function convert(::Type{UnipartiteNetwork}, N::T) where {T <: BipartiteNetwork}
    itype = _interaction_type(N)
    S = copy(species(N))
    B = zeros(itype, (richness(N), richness(N)))
    B[1:size(N)[1],size(N)[1]+1:richness(N)] = N.edges
    return UnipartiteNetwork(sparse(B), S)
end

"""
    convert(::Type{UnipartiteProbabilisticNetwork}, N::T) where {T <: BipartiteProbabilisticNetwork}

Projects a probabilistic bipartite network in its unipartite representation.
"""
function convert(::Type{UnipartiteProbabilisticNetwork}, N::T) where {T <: BipartiteProbabilisticNetwork}
    itype = _interaction_type(N)
    S = copy(species(N))
    B = zeros(itype, (richness(N), richness(N)))
    B[1:size(N)[1],size(N)[1]+1:richness(N)] = N.edges
    return UnipartiteProbabilisticNetwork(sparse(B), S)
end

"""
    convert(::Type{UnipartiteQuantitativeNetwork}, N::T) where {T <: BipartiteQuantitativeNetwork}

Projects a quantitative bipartite network in its unipartite representation.
"""
function convert(::Type{UnipartiteQuantitativeNetwork}, N::T) where {T <: BipartiteQuantitativeNetwork}
    itype = _interaction_type(N)
    S = copy(species(N))
    B = zeros(itype, (richness(N), richness(N)))
    B[1:size(N)[1],size(N)[1]+1:richness(N)] = N.edges
    return UnipartiteQuantitativeNetwork(sparse(B), S)
end

"""
    convert(::Type{UnipartiteNetwork}, N::T) where {T <: UnipartiteQuantitativeNetwork}

Convert a unipartite quantitative network to a unipartite binary network. This
amounts to *removing* the quantitative information.
"""
function convert(::Type{UnipartiteNetwork}, N::T) where {T <: UnipartiteQuantitativeNetwork}
    S = copy(species(N))
    B = dropzeros(N.edges)
    return UnipartiteNetwork(B .!= zero(eltype(B)), S)
end

"""
    convert(::Type{BipartiteNetwork}, N::T) where {T <: BipartiteQuantitativeNetwork}

Convert a bipartite quantitative network to a bipartite binary network. This
amounts to *removing* the quantitative information.
"""
function convert(::Type{BipartiteNetwork}, N::NT) where {NT <: BipartiteQuantitativeNetwork}
    T = copy(species(N; dims=1))
    B = copy(species(N; dims=2))
    M = dropzeros(N.edges)
    return BipartiteNetwork(M .!= zero(eltype(M)), T, B)
end

"""
    convert(::Type{AbstractUnipartiteNetwork}, N::AbstractBipartiteNetwork)

Projects any bipartite network in its unipartite representation. This function
will call the correct type-to-type `convert` function depending on the type of
the input network.

The type to be converted to *must* be `AbstractUnipartiteNetwork` -- for
example, converting a bipartite probabilistic network to a probabilistic
unipartite network is not a meaningful operation.
"""
function convert(::Type{AbstractUnipartiteNetwork}, N::AbstractBipartiteNetwork)
    if typeof(N) <: QuantitativeNetwork
        return convert(UnipartiteQuantitativeNetwork, N)
    end
    if typeof(N) <: BinaryNetwork
        return convert(UnipartiteNetwork, N)
    end
    if typeof(N) <: ProbabilisticNetwork
        return convert(UnipartiteProbabilisticNetwork, N)
    end
end

"""
    convert(::Type{BinaryNetwork}, N::QuantitativeNetwork)

Projects any bipartite network in its unipartite representation. This function
will call the correct type-to-type `convert` function depending on the type of
the input network.

This function does *not* work for probabilistic networks. The operation of
generating a deterministic network from a probabilistic one is different from a
simple conversion: it can be done either through random draws, or by selecting
only interactions with a probability greater than 0 (`N>0.0` will do this).
"""
function convert(::Type{BinaryNetwork}, N::QuantitativeNetwork)
    if typeof(N) <: BipartiteQuantitativeNetwork
        return convert(BipartiteNetwork, N)
    end
    if typeof(N) <: UnipartiteQuantitativeNetwork
        return convert(UnipartiteNetwork, N)
    end
end

# Conversion to bipartite

type_pairs = [
    (:BipartiteNetwork, :UnipartiteNetwork),
    (:BipartiteProbabilisticNetwork, :UnipartiteProbabilisticNetwork),
    (:BipartiteQuantitativeNetwork, :UnipartiteQuantitativeNetwork)
    ]

for comb in type_pairs
    t1, t2 = comb
    eval(
    quote
        """
            convert(::Type{$($t1)}, N::T) where {T <: $($t2)}

        Projects a unipartite network (specifically, a `$($t1)`) to its bipartite
        representation. The following checks are performed.

        First, the network *cannot* be degenerate, since species with no interactions
        cannot be assigned to a specific level. Second, the species cannot have both in
        and out degree. If these conditions are met, the bipartite network will be
        returned.
        """
        function convert(::Type{$t1}, N::T) where {T <: $t2}
            isdegenerate(N) && throw(ArgumentError("Impossible to convert a degenerate unipartite network into a bipartite one"))
            d1 = degree(N; dims=1)
            d2 = degree(N; dims=2)
            for s in species(N)
                (d1[s]*d2[s] == 0) || throw(ArgumentError("Species $s has both in and out degree"))
            end
            top_species = collect(keys(filter(p -> p.second == zero(eltype(N.edges)), d2)))
            bot_species = collect(keys(filter(p -> p.second == zero(eltype(N.edges)), d1)))
            A = zeros(_interaction_type(N), (length(top_species), length(bot_species)))
            B = $t1(sparse(A), top_species, bot_species)
            for s1 in species(B; dims=1), s2 in species(B; dims=2)
                B[s1,s2] = N[s1, s2]
            end
            return B
        end
    end
    )
end


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

using BAT
using Test

using LinearAlgebra, Random
using ArraysOfArrays, Distributions, PDMats, StatsBase

@testset "generic_density" begin
    mvec = [-0.3, 0.3]
    cmat = [1.0 1.5; 1.5 4.0]
    Σ = @inferred PDMat(cmat)
    mvnorm = @inferred  MvNormal(mvec, Σ)

    density = let mvnorm = mvnorm
        @inferred BAT.GenericDensity(params -> logpdf(mvnorm, params))
    end

    params = [0.0, 0.0]

    @testset "BAT.density_logval" begin
        @test @inferred(BAT.density_logval(density, params)) ≈ logpdf(mvnorm, params)
        @test @inferred(BAT.density_logval(density, params)) ≈ logpdf(mvnorm, params)
    end

    @testset "parent" begin
        @test @inferred parent(density) == density.log_f
    end
end


text:
import Base.LinAlg.BlasFloat

export eigvals_lanczos

function lanczos!{T}(K::KrylovSubspace{T})
    m = K.n
    αs = Array(T, m)
    βs = Array(T, m-1)
    for j=1:m-1
        w = nextvec(K)
        if j>1 w -= βs[j-1]*K.v[1] end
        w, y = orthogonalize(w, K, 1)
        αs[j] = y[1]
        βs[j] = convert(T, norm(w))
        append!(K, w/βs[j])
    end
    αs[m]= dot(nextvec(K), lastvec(K))
    SymTridiagonal(αs, βs)
end

function eigvals_lanczos(A, neigs::Int=size(A,1); tol::Real = size(A,1)^3*eps(real(eltype(A))), maxiter::Integer = size(A,1))
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    e1 = eiglancz_method!(history, A, neigs; maxiter=maxiter, tol=tol)
    e1, history
end

function eiglancz_method!(log::ConvergenceHistory, A, neigs::Int=size(A,1); tol::Real = size(A,1)^3*eps(real(eltype(A))), maxiter::Integer = size(A,1))
    K = KrylovSubspace(A, size(A, 1), 2) #In Lanczos, only remember the last two vectors
    initrand!(K)
    e1 = eigvals(lanczos!(K), 1:neigs)
    for iter=1:maxiter
        nextiter!(log,mvps=1)
        e0, e1 = e1, eigvals(lanczos!(K), 1:neigs)
        resnorm = norm(e1-e0)
        push!(log, :resnorm, resnorm)
        resnorm < tol && (setconv(log, resnorm>=0); break)
    end
    shrink!(log)
    e1
end


text:
struct StdSquare{T}
    half_side::T
end

struct StdCircle{T}
    radius::T
end

function get_projection(square::StdSquare, position::AbstractVector)
    half_side = square.half_side
    return clamp.(position, -half_side, half_side)
end

function is_colliding(square::StdSquare, circle::StdCircle, position::AbstractVector)
    projection = get_projection(square, position)
    vec = position .- projection
    radius = circle.radius
    return sum(vec .^ 2) < radius * radius
end

function is_player_colliding(obstacle_map, player_position_wu, player_radius_wu::T) where {T}
    height_tile_map_tu, width_tile_map_tu = size(obstacle_map)

    square = StdSquare(convert(T, 0.5))
    circle = StdCircle(player_radius_wu)

    i_player_position_tu = wu_to_tu(player_position_wu[1])
    j_player_position_tu = wu_to_tu(player_position_wu[2])

    for j in j_player_position_tu - 1 : j_player_position_tu + 1
        for i in i_player_position_tu - 1 : i_player_position_tu + 1
            tile_position_wu = similar(player_position_wu)
            tile_position_wu[1] = i - convert(T, 0.5)
            tile_position_wu[2] = j - convert(T, 0.5)
            if obstacle_map[i, j] && is_colliding(square, circle, player_position_wu .- tile_position_wu)
                return true
            end
        end
    end

    return false
end


text:
# ---
# title: 1395. Count Number of Teams
# id: problem1395
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/count-number-of-teams/description/>
# hidden: true
# ---
#
# There are `n` soldiers standing in a line. Each soldier is assigned a
# **unique** `rating` value.
#
# You have to form a team of 3 soldiers amongst them under the following rules:
#
#   * Choose 3 soldiers with index (`i`, `j`, `k`) with rating (`rating[i]`, `rating[j]`, `rating[k]`).
#   * A team is valid if:  (`rating[i] < rating[j] < rating[k]`) or (`rating[i] > rating[j] > rating[k]`) where (`0 <= i < j < k < n`).
#
# Return the number of teams you can form given the conditions. (soldiers can be
# part of multiple teams).
#
#
#
# **Example 1:**
#
#
#
#     Input: rating = [2,5,3,4,1]
#     Output: 3
#     Explanation: We can form three teams given the conditions. (2,3,4), (5,4,1), (5,3,1).
#
#
# **Example 2:**
#
#
#
#     Input: rating = [2,1,3]
#     Output: 0
#     Explanation: We can't form any team given the conditions.
#
#
# **Example 3:**
#
#
#
#     Input: rating = [1,2,3,4]
#     Output: 4
#
#
#
#
# **Constraints:**
#
#   * `n == rating.length`
#   * `1 <= n <= 200`
#   * `1 <= rating[i] <= 10^5`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ASTInterpreter

function triple_for_sess(sess)
    if isa(Gallium.getarch(sess), Gallium.X86_64.X86_64Arch)
        return "x86_64-linux-gnu"
    else
        return "i386-linux-gnu"
    end
end

function ASTInterpreter.execute_command(state, stack, ::Val{:disas}, command)
    parts = split(command, ' ')
    stacktop = 0
    if length(parts) > 1
        ip = parse(Int, parts[2], 16)
    else
        x = isa(stack, Gallium.NativeStack) ? stack.stack[end] : stack
        stacktop = (x.stacktop?0:1); ip = x.ip
    end
    base, loc, insts = get_insts(state.top_interp.session, state.top_interp.modules, ip)
    disasm_around_ip(STDOUT, insts, UInt64(ip-loc-base-stacktop);
        ipbase=base+loc, circular = false, triple = triple_for_sess(state.top_interp.session))
    return false
end

task_single_step!(timeline) = Gallium.single_step!(timeline)

function ASTInterpreter.execute_command(state, stack::Union{Gallium.NativeStack,Gallium.CStackFrame}, ::Val{:si}, command)
    task_single_step!(state.top_interp.session)
    update_stack_same_frame!(state)
    return true
end

function compute_current_line_range(state, stack)
    mod, base, ip = Gallium.modbaseip_for_stack(state, stack)
    linetab, lip = Gallium.obtain_linetable(state, stack)
    sm = start(linetab)
    local current_entry
    local newentry
    # Start by finding the entry that we're in
    while true
        newentry, sm = next(linetab, sm)
        newentry.address > lip && break
        current_entry = newentry
    end
    range = origrange = current_entry.address:(newentry.address-1)
    # Merge any subsequent entries at the same line
    while newentry.line == current_entry.line
        newentry, sm = next(linetab, sm)
        range = first(origrange):(newentry.address-1)
    end
    range += UInt64(ip-lip)
    range
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.NativeStack,Gallium.CStackFrame}, ::Val{:n}, command)
    session = state.top_interp.session
    range = compute_current_line_range(state, stack)
    step_over(session, range)
    update_stack_same_frame!(state)
    return true
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.NativeStack,Gallium.CStackFrame}, ::Val{:ip}, command)
    x = isa(stack, Gallium.NativeStack) ? stack.stack[end] : stack
    println(x.ip)
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.NativeStack,Gallium.CStackFrame}, ::Val{:nb}, command)
    session = state.top_interp.session
    # First determine the ip of the next branch
    x = isa(stack, Gallium.NativeStack) ? stack.stack[end] : stack
    base, loc, insts = get_insts(session, state.top_interp.modules, x.ip)
    ctx = DisAsmContext()
    Offset = UInt(x.ip - loc - base)
    branchip = 0
    while Offset <= sizeof(insts)
        (Inst, InstSize) = getInstruction(insts, Offset; ctx = ctx)
        if mayAffectControlFlow(Inst,ctx)
            branchip = base + loc + Offset
            break
        end
        Offset += InstSize
        free(Inst)
    end
    @assert branchip != 0
    bp = Gallium.breakpoint(session, branchip)
    task_step_until_bkpt!(session)
    Gallium.disable(bp)
    update_stack_same_frame!(state)
    return true
end

function ASTInterpreter.print_status(state, x::Gallium.CStackFrame; kwargs...)
    session = state.top_interp.session
    modules = state.top_interp.modules
    print("Stopped in function ")
    found, symb = symbolicate_frame(session, modules, x)
    println(symb)

    if x.line != 0
        code = ASTInterpreter.readfileorhist(x.file)
        if code !== nothing
            ASTInterpreter.print_sourcecode(
                code, x.line, x.declfile, x.declline)
            return
        end
    end
    ipoffset = 0
    ipbase = x.ip
    if found
        base, loc, insts = get_insts(session, modules, ipbase)
        ipbase = base+loc
        ipoffset = UInt64(x.ip-loc-base-(x.stacktop?0:1))
    else
        try
          insts = Gallium.load(session, Gallium.RemotePtr{UInt8}(x.ip), 40)
        catch
          warn("Failed to load instructions at address ($(x.ip)). Incorrect unwind info?")
          return
        end
    end
    try
        # Try disassembling at the start of the function, highlighting the
        # current ip.
        disasm_around_ip(STDOUT, insts, ipoffset; ipbase=ipbase, triple = triple_for_sess(state.top_interp.session))
    catch e
        warn("Failed to obtain instructions from object files ($e). Incorrect unwind info?")
        # This could have failed for a variety of reasons (missing unwind info,
        # self modifying code, etc). If so, get the instructions directly
        # from the target
        insts = Gallium.load(session, Gallium.RemotePtr{UInt8}(x.ip), 40)
        ipoffset = 0
        disasm_around_ip(STDOUT, insts, ipoffset; ipbase=ipbase, triple = triple_for_sess(state.top_interp.session))
    end
end

cxx"""#include <cxxabi.h>"""
function demangle(name)
    startswith(name,"_Z") || return name
    status = Ref{Cint}()
    bufsize = Ref{Csize_t}(0)
    str = icxx"""
        abi::__cxa_demangle($(pointer(name)),nullptr,
        &$bufsize, &$status);
    """
    @assert status[] == 0
    ret = unsafe_string(str)
    Libc.free(str)
    ret
end

function symbolicate_frame(session, modules, x)
    found = false
    symb = "Unknown Function"
    try
        found, symb = Gallium.Unwinder.symbolicate(session, modules, UInt64(x.ip))
        symb = demangle(symb)
        !found && (symb = "Most likely $symb")
    catch err
        (!isa(err, ErrorException) || !contains(err.msg, "found")) && rethrow(err)
    end
    found, symb
end

function ASTInterpreter.print_frame(state, io, num, x::Gallium.CStackFrame)
    session = state.top_interp.session
    modules = state.top_interp.modules
    print(io, "[$num] ")
    found, symb = symbolicate_frame(session, modules, x)
    print(io, symb, " ")
    if x.line != 0
      print(io, " at ",x.file,":",x.line)
    end
    println(io)
end

function ASTInterpreter.execute_command(state, stack, ::Val{:c}, command)
    try
        Gallium.continue!(state.top_interp.session; only_current_tgid = true)
    catch err
        !isa(err, InterruptException) && rethrow(err)
    end
    update_stack!(state)
    return true
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:reg}, command)
    ns = state.top_interp
    @assert isa(ns, Gallium.NativeStack)
    RC = ns.RCs[end-(state.level-1)]
    regname = Symbol(split(command,' ')[2])
    inverse_map = isa(Gallium.getarch(state.top_interp.session),Gallium.X86_64.X86_64Arch) ?
        Gallium.X86_64.inverse_dwarf : Gallium.X86_32.inverse_dwarf
    if !haskey(inverse_map, regname)
        print_with_color(:red, STDOUT, "No such register\n")
        return false
    end
    show(UInt(Gallium.get_dwarf(RC, inverse_map[regname])))
    println(); println()
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:regs}, command)
    ns = state.top_interp
    @assert isa(ns, Gallium.NativeStack)
    RC = ns.RCs[end-(state.level-1)]
    show(STDOUT, RC)
    println(); println()
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:unwind}, command)
    ns = state.top_interp
    newRC = Gallium.Unwinder.unwind_step(ns.session, ns.modules, ns.RCs[end-(state.level-1)];
      allow_frame_based = false, stacktop=ns.stack[end-(state.level-1)].stacktop)[2]
    @show newRC
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:symbolicate}, command)
    session = state.top_interp.session
    modules = state.top_interp.modules
    x = isa(stack, Gallium.NativeStack) ? stack.stack[end] : stack
    approximate, name = Gallium.Unwinder.symbolicate(session, modules, UInt64(x.ip))
    name = demangle(name)
    println(approximate ? "Most likely " : "", name)
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:modules}, command)
    modules = state.top_interp.modules
    modules = sort(collect(modules), by = x->x[1])
    for (base, mod) in modules
        show(UInt64(base)); print('-'); show(UInt64(base)+mod.sz); println()
    end
    return false
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:finish}, command)
    ns = state.top_interp
    session = state.top_interp.session
    @assert isa(ns, Gallium.NativeStack)
    parentRC = ns.RCs[end-(state.level)]
    theip = Gallium.ip(parentRC)
    step_to_address(session, theip)
    update_stack!(state)
    return true
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:b}, command)
    session = state.top_interp.session
    modules = state.top_interp.modules
    symbol = Symbol(split(command,' ')[2])
    bp = Gallium.breakpoint(session, modules, symbol)
    show(bp)
    return false
end

function dwarf2Cxx(dbgs, dwarfT)
    if DWARF.tag(dwarfT) == DWARF.DW_TAG_pointer_type ||
            DWARF.tag(dwarfT) == DWARF.DW_TAG_array_type
        dwarfT = get(DWARF.extract_attribute(dwarfT,DWARF.DW_AT_type))
        return Cxx.pointerTo(Cxx.instance(RemoteClang), dwarf2Cxx(dbgs, dwarfT.value))
    else
        name = DWARF.extract_attribute(dwarfT,DWARF.DW_AT_name)
        name = bytestring(get(name).value,StrTab(dbgs.debug_str))
        return cxxparse(Cxx.instance(RemoteClang),name,true)
    end
end

function iterate_frame_variables(state, stack, found_cb, not_found_cb)
    mod, base, theip = Gallium.modbaseip_for_stack(state, stack)
    lip = Gallium.compute_ip(Gallium.dhandle(mod),base,theip)
    dbgs = debugsections(Gallium.dhandle(mod))
    ns = state.top_interp
    @assert isa(ns, Gallium.NativeStack)
    RC = ns.RCs[end-(state.level-1)]

    Gallium.iterate_variables(RC, found_cb, not_found_cb, dbgs, lip)
end

function realize_remote_value(T, val, getreg)
    if isa(val, DWARF.Expressions.MemoryLocation)
        val = Gallium.load(timeline, RemotePtr{T}(val.i))
    elseif isa(val, DWARF.Expressions.RegisterLocation)
        val = reinterpret(T, [getreg(val.i)])[]
    end
    val
end

function ASTInterpreter.execute_command(state, stack::Union{Gallium.CStackFrame,Gallium.NativeStack}, ::Val{:vars}, command)
    function found_cb(dbgs, vardie, getreg, name, val)
        dwarfT = get(DWARF.extract_attribute(vardie,DWARF.DW_AT_type))
        try
            T = Cxx.juliatype(dwarf2Cxx(dbgs, dwarfT.value))
            val = realize_remote_value(T, val, getreg)
        end
        @show (name, val)
    end
    iterate_frame_variables(state, stack, found_cb, (dbgs, vardie, name)->nothing)

    return false
end


text:
"""
    struct Minimise <: CombinatorialObjective end

Minimises the natural objective function for this combinatorial problem.
(Typically, it corresponds to the cost of the solution.)
"""
struct Minimise <: CombinatorialObjective end

"""
    struct Maximise <: CombinatorialObjective end

Maximises the natural objective function for this combinatorial problem.
(Typically, it corresponds to the cost of the solution.)
"""
struct Maximise <: CombinatorialObjective end