text:
module Regularizer

export regularizer, createElasticOperatorCentered, createElasticOperatorStaggered, createDiffusiveOperatorCentered
export createCurvatureOperatorCentered

function regularizer(deformationField::Array{Float64,1}, operator)

    d2functionValue = operator'*operator
    dfunctionValue  = d2functionValue*deformationField
    functionValue   = (0.5 .* deformationField'*dfunctionValue)[1]
    return [functionValue, dfunctionValue, d2functionValue]

end

function createElasticOperatorCentered(h::Array{Float64,1},m::Array{Int64,1}; mu::Number = 1, lambda::Number = 0)

  m = m.-1

  a   = sqrt(mu)
  b   = sqrt(mu+lambda)
  dx(k) = spdiagm((-1.*ones(m[k]),ones(m[k])),0:1,m[k],m[k]+1)/h[k]
  av(k) = spdiagm((ones(m[k])/2,ones(m[k])/2),0:1,m[k],m[k]+1)
  D1   = kron(speye(m[2]+1),dx(1))
  D2   = kron(dx(2),speye(m[1]+1))
  A1   = kron( av(2) , speye(m[1]))
  A2   = kron( speye(m[2]),  av(1))
  p1,p2 = size(D1)
  p3,p4 = size(D2)

  B = [ a*D1 spzeros(p1,p2);
       a*D2 spzeros(p3,p4);
       spzeros(p1,p2) a*D1;
       spzeros(p3,p4) a*D2;
       b*A1*D1 b*A2*D2 ]

  return B

end

# function createElasticOperatorStaggered(h::Array{Float64,1},m::Array{Int64,1}; mu::Number = 1, lambda::Number = 0)
#
#     a   = sqrt(mu)
#     b   = sqrt(mu+lambda)
#
#     dx(h,m) = spdiagm((-ones(m,1),ones(m,1)),[0,1],m,m+1)/h
#
#     D11 = kron(dx(h[2],m[2]),speye(m[1]))
#     D21 = kron(speye(m[2]+1),dx(h[1],m[1]-1))
#     D12 = kron(dx(h[2],m[2]-1),speye(m[1]+1))
#     D22 = kron(speye(m[2]),dx(h[1],m[1]))
#
#     B = spzeros(2*size(D11,1)+size(D21,1)+size(D12,1)+size(D22,1),size(D11,2)+size(D12,2))
#     B[1:size(D11,1), 1:size(D11,2)]                         = a*D11
#     B[size(D11,1)+1:size(D11,1)+size(D21,1), 1:size(D11,2)] = a*D21
#     B[size(D11,1)+size(D21,1)+1:size(D11,1)+size(D21,1)+size(D12,1), size(D11,2)+1:end]                         = a*D12
#     B[size(D11,1)+size(D21,1)+size(D12,1)+1:size(D11,1)+size(D21,1)+size(D12,1)+size(D22,1), size(D11,2)+1:end] = a*D22
#     B[size(D11,1)+size(D21,1)+size(D12,1)+size(D22,1)+1:end, 1:size(D11,2)]     = b*D11
#     B[size(D11,1)+size(D21,1)+size(D12,1)+size(D22,1)+1:end, size(D11,2)+1:end] = b*D22
#
#     return B
#
# end

function createDiffusiveOperatorCentered(h::Array{Float64,1},m::Array{Int64,1})

  d(k) = spdiagm((-ones(m[k]-1,1),ones(m[k]-1,1)),[0 1],m[k]-1,m[k])/(h[k])
  dx = d(1)
  dy = d(2)
  D1 = kron(speye(m[2]),dx)
  D2 = kron(dy,speye(m[1]))
  p1,p2 = size(D1)
  p3,p4 = size(D2)
  B = [ D1 spzeros(p1,p2);
        D2 spzeros(p3,p4);
        spzeros(p1,p2) D1;
        spzeros(p3,p4) D2
      ]
end

function createCurvatureOperatorCentered(h::Array{Float64,1},m::Array{Int64,1})

  d(k) = spdiagm((-ones(m[k]-2,1),2*ones(m[k]-2,1), -ones(m[k]-2,1)),[0 1 2],m[k]-2,m[k])/(h[k])
  dx = d(1); #dx[1,1]=1; dx[end, end]=1
  dy = d(2); #dy[1,1]=1; dy[end, end]=1
  D1 = kron(speye(m[2]),dx)
  D2 = kron(dy,speye(m[1]))
  p1,p2 = size(D1)
  p3,p4 = size(D2)
  B = [ D1 spzeros(p1,p2);
        D2 spzeros(p3,p4);
        spzeros(p1,p2) D1;
        spzeros(p3,p4) D2
      ]
end

end#


text:
module SphericalHarmonicModes

export LM,ML,L₂L₁
export modeindex,l_range,m_range,l₂_range,l₁_range

abstract type ModeRange end
abstract type SHModeRange <: ModeRange end
abstract type ModeProduct <: ModeRange end

Base.eltype(::SHModeRange) = Tuple{Int,Int}

struct LM <: SHModeRange
	l_min :: Int
	l_max :: Int
	m_min :: Int
	m_max :: Int

	function LM(l_min,l_max,m_min,m_max)
		check_if_st_range_is_valid(l_min,l_max,m_min,m_max)
		l_min = max(l_min,max(m_min,-m_max))
		new(l_min,l_max,m_min,m_max)
	end
end

struct ML <: SHModeRange
	l_min :: Int
	l_max :: Int
	m_min :: Int
	m_max :: Int

	function ML(l_min,l_max,m_min,m_max)
		check_if_st_range_is_valid(l_min,l_max,m_min,m_max)
		l_min = max(l_min,max(m_min,-m_max))
		new(l_min,l_max,m_min,m_max)
	end
end

struct L₂L₁ <: ModeProduct
	l₁_min :: Int
	l₁_max :: Int
	Δl_max :: Int
	l₂_min :: Int
	l₂_max :: Int

	function L₂L₁(l_min::Integer,l_max::Integer,
		Δl_max::Integer,l₂_min::Integer,l₂_max::Integer)

		map(check_if_non_negative,(l_min,l_max,Δl_max,l₂_min,l₂_max))

		# Alter the ranges if necessary
		l_min = max(l_min,l₂_min-Δl_max,0)
		l₂_min = max(l_min-Δl_max,0,l₂_min)
		l₂_max = max(min(max(l_max+Δl_max,0),l₂_max),l₂_min)
		l_max = min(l_max,l₂_max+Δl_max)

		new(l_min,l_max,Δl_max,l₂_min,l₂_max)
	end
end

Base.eltype(::L₂L₁) = Tuple{Int,Int}

# Throw errors if values are invalid
struct OrderError{T} <: Exception
	var :: String
	low :: T
	high :: T
end

struct tRangeError <: Exception
	l_max :: Integer
	m :: Integer
end

struct ModeMissingError{M,T} <: Exception
	firstmode :: T
	secondmode :: T
	itr :: M
end

struct NegativeDegreeError{T<:Integer} <: Exception
	l :: T
end

struct InvalidModeError{ML<:Integer,Tt<:Integer} <: Exception
	l :: ML
	m :: Tt
end

struct NonContiguousError <: Exception end

Base.showerror(io::IO, e::tRangeError) = print(io," m = ", e.m,
		" does not satisfy ",-e.l_max," ⩽ m ⩽ ",e.l_max)

Base.showerror(io::IO, e::OrderError) = print(io,e.var,"min = ",e.low,
	" is not consistent with ",e.var,"max = ",e.high)

Base.showerror(io::IO, e::ModeMissingError{<:SHModeRange}) =
	print(io,"Mode with (l=",e.firstmode,",m=",e.secondmode,")",
	" is not included in the range given by ",e.itr)

Base.showerror(io::IO, e::ModeMissingError{L₂L₁}) =
	print(io,"Mode with (l₂=",e.firstmode,",l₁=",e.secondmode,")",
	" is not included in the range given by ",e.itr)

Base.showerror(io::IO, e::NegativeDegreeError) = print(io,"l = ",e.l,
	" does not correspond to a valid mode")

Base.showerror(io::IO, e::InvalidModeError) = print(io,"(l=",e.l,",m=",e.m,")",
	" is not a valid mode. |m| <= l is not satisfied")

Base.showerror(io::IO, e::NonContiguousError) = print(io,
	"ModeRanges only support contiguous indexing")

function check_if_non_negative(l::Integer)
	l >= 0 || throw(NegativeDegreeError(l))
end

function check_if_st_range_is_valid(l_min,l_max,m_min,m_max)
	check_if_non_negative.((l_min,l_max))
	if m_min > m_max
		throw(OrderError("m",m_min,m_max))
	end
	if l_min > l_max
		throw(OrderError("l",l_min,l_max))
	end
	if abs(m_max) > l_max
		throw(tRangeError(l_max,m_max))
	end
	if abs(m_min) > l_max
		throw(tRangeError(l_max,m_min))
	end
end

@inline function check_if_valid_mode(l::Integer,m::Integer)
	abs(m) <= l || throw(InvalidModeError(l,m))
end

@inline function check_if_mode_present(mr,l,m)
	(l,m) in mr || throw(ModeMissingError(l,m,mr))
end

# Constructors. Both ML and LM are constructed identically,
# so we can dispatch on the supertype
(::Type{T})(l_min::Integer,l_max::Integer) where {T<:SHModeRange} = T(l_min,l_max,-l_max,l_max)
(::Type{T})(l::Integer) where {T<:SHModeRange} = T(l,l)

(::Type{T})(l_range::AbstractUnitRange{<:Integer},
	m_range::AbstractUnitRange{<:Integer}) where {T<:SHModeRange} =
	T(extrema(l_range)...,extrema(m_range)...)

(::Type{T})(l_range::AbstractUnitRange{<:Integer},m::Integer) where {T<:SHModeRange} =
	T(extrema(l_range)...,m,m)

(::Type{T})(l::Integer,m_range::AbstractUnitRange{<:Integer}) where {T<:SHModeRange} =
	T(l,l,extrema(m_range)...)

(::Type{T})(l_range::AbstractUnitRange{<:Integer}) where {T<:SHModeRange} =
	T(extrema(l_range)...)

function L₂L₁(l_min::Integer,l_max::Integer,Δl_max::Integer,
	l₂_min::Integer=max(l_min-Δl_max,0))

	L₂L₁(l_min,l_max,Δl_max,l₂_min,l_max+Δl_max)
end

L₂L₁(l_min::Integer,l_max::Integer,mr::SHModeRange,args...) = L₂L₁(l_min,l_max,mr.l_max,args...)

L₂L₁( l_range::AbstractUnitRange{<:Integer},Δl_max::Integer,
	l₂_range::AbstractUnitRange{<:Integer}) =
	L₂L₁(extrema(l_range)...,Δl_max,extrema(l₂_range)...)

L₂L₁(l_min::Integer,l_max::Integer,Δl_max::Integer,
	l₂_range::AbstractUnitRange{<:Integer}) =
	L₂L₁(l_min,l_max,Δl_max,extrema(l₂_range)...)

L₂L₁( l_range::AbstractUnitRange{<:Integer},Δl_max::Union{Integer,<:SHModeRange},
	args...) = L₂L₁(extrema(l_range)...,Δl_max,args...)

# Get the ranges of the modes

@inline l_range(mr::SHModeRange) = mr.l_min:mr.l_max
@inline m_range(mr::SHModeRange) = mr.m_min:mr.m_max

@inline l₁_range(mr::L₂L₁) = mr.l₁_min:mr.l₁_max
@inline l₂_range(mr::L₂L₁) = mr.l₂_min:mr.l₂_max

@inline l_range(mr::SHModeRange,m::Integer) = max(abs(m),mr.l_min):mr.l_max

@inline m_range(mr::SHModeRange,l::Integer) = max(-l,mr.m_min):min(l,mr.m_max)

@inline function l₂_range(mr::L₂L₁,l::Integer)
	max(l - mr.Δl_max,mr.l₂_min):min(l + mr.Δl_max,mr.l₂_max)
end

function _length(mr::LM)

	# We evaluate Sum[l_max - Max[Abs[m], l_min] + 1] piecewise in
	# Mathematica. There are four cases to consider
	# depending on which of |m| and l_min is larger

	l_min,l_max,m_min,m_max = mr.l_min,mr.l_max,mr.m_min,mr.m_max

	N = 0

	up = min(-1, m_max)
	lo = max(m_min,-l_min)
	if up >= lo
		N += (1 + l_max - l_min)*(1 - lo + up)
	end

	up = min(m_max, l_min)
	lo = max(m_min, 0)
	if up >= lo
		N += (1 + l_max - l_min)*(1 - lo + up)
	end

	up = m_max
	lo = max(m_min, l_min + 1)
	if up >= lo
		N += div(-((1 + up - lo)*(-2 - 2l_max + up + lo)),2)
	end

	up = min(-1, -1 - l_min, m_max)
	lo = max(-l_max, m_min)
	if up >= lo
		N += div(-((-1 + lo - up)*(2 + 2l_max + lo + up)),2)
	end

	return N
end

function _length(mr::ML)

	# We evaluate Sum[min(l,m_max)-max(-l,m_min)+1] piecewise in
	# Mathematica. There are four cases to consider
	# depending on l, m_max and m_min

	l_min,l_max,m_min,m_max = mr.l_min,mr.l_max,mr.m_min,mr.m_max

	N = 0

	up = min(-m_min, m_max,l_max)
	lo = l_min
	if up >= lo
		N += -(-1 + lo - up)*(1 + lo + up)
	end

	up = min(-m_min, l_max)
	lo = max(l_min, m_max + 1, -m_max)
	if up >= lo
		N += div(-((-1 + lo - up)*(2 + 2m_max + lo + up)),2)
	end

	up = min(m_max, l_max)
	lo = max(l_min, -m_min + 1, m_min)
	if up >= lo
		N += div(-((-1 + lo - up)*(2 - 2m_min + lo + up)),2)
	end

	up = l_max
	lo = max(-m_min + 1, m_max + 1,l_min)
	if up >= lo
		N += (1 + m_max - m_min)*(1 + up - lo)
	end

	return N
end

function _length(mr::L₂L₁)

	# Number of modes is given by count(l₂_range(mr,l) for l in l_range(mr))
	l_min,l_max,l₂_min,l₂_max = mr.l₁_min,mr.l₁_max,mr.l₂_min,mr.l₂_max
	Δl_max = mr.Δl_max

    N = 0

    up = min(l_max, l₂_max - Δl_max)
    lo = max(l_min, l₂_min + Δl_max)
    if up >= lo
    	N += (1 + 2*Δl_max)*(1 - lo + up)
    end

    up = min(l_max,l₂_max + Δl_max)
    lo = max(l_min, l₂_min + Δl_max, l₂_max - Δl_max + 1)
    if up >= lo
    	N += div((-1 + lo - up)*(-2*(1 + Δl_max + l₂_max) + lo + up),2)
    end

    up = min(l_max, l₂_min + Δl_max - 1, l₂_max - Δl_max)
    lo = max(l_min,l₂_min - Δl_max)
    if up >= lo
    	N += div(-((-1 + lo - up)*(2 + 2Δl_max - 2l₂_min + lo + up)),2)
    end

    up = min(l_max, l₂_min + Δl_max - 1)
    lo = max(l_min,l₂_max - Δl_max + 1)
    if up >= lo
    	N += (1 + l₂_max - l₂_min)*(1 - lo + up)
    end

    return N
end

Base.length(mr::ModeRange) = _length(mr)

Base.firstindex(mr::ModeRange) = 1
Base.lastindex(mr::ModeRange) = length(mr)

# Size and axes behave similar to vectors
@inline Base.size(mr::ModeRange) = (length(mr),)
@inline Base.size(mr::ModeRange,d::Integer) = d == 1 ? length(mr) : 1

@inline Base.axes(mr::ModeRange) = (Base.OneTo(length(mr)),)
@inline Base.axes(mr::ModeRange,d::Integer) = d == 1 ? Base.OneTo(length(mr)) : Base.OneTo(1)

# Convenience function to generate the first step in the iteration
@inline first_m(mr::LM) = mr.m_min
@inline first_l(mr::LM) = first(l_range(mr,first_m(mr)))
@inline first_l(mr::ML) = mr.l_min
@inline first_m(mr::ML) = first(m_range(mr,first_l(mr)))

@inline first_l₁(mr::L₂L₁) = max(mr.l₁_min,mr.l₂_min-mr.Δl_max)
@inline first_l₂(mr::L₂L₁) = first(l₂_range(mr,first_l₁(mr)))

function Base.iterate(mr::LM, state=((first_l(mr),first_m(mr)), 1))

	(l,m), count = state

	if count > _length(mr)
		return nothing
	end

	if l == mr.l_max
		next_m = m + 1
		next_l = max(mr.l_min,abs(next_m))
	else
		next_m = m
		next_l = l + 1
	end

	return (l,m), ((next_l,next_m), count + 1)
end

function Base.iterate(mr::ML, state=((first_l(mr),first_m(mr)), 1))

	(l,m), count = state

	if count > length(mr)
		return nothing
	end

	if m == min(mr.m_max,l)
		next_l = l + 1
		next_m = max(mr.m_min,-next_l)
	else
		next_l = l
		next_m = m + 1
	end

	return (l,m), ((next_l,next_m), count + 1)
end

function Base.iterate(mr::L₂L₁,state=((first_l₂(mr),first_l₁(mr)), 1))

	(l₂,l₁),count = state
	if count > length(mr)
		return nothing
	end

	l₂_range_l₁ = l₂_range(mr,l₁)

	if l₂ == last(l₂_range_l₁)
		next_l₁ = l₁ + 1
		next_l₂ = first(l₂_range(mr,next_l₁))
	else
		next_l₁ = l₁
		next_l₂ = l₂ + 1
	end

	return (l₂,l₁),((next_l₂,next_l₁),count+1)
end

function _in((l,m)::Tuple{<:Integer,<:Integer},mr::ML)
	(abs(m)<=l) && (mr.l_min <= l <= mr.l_max) && (mr.m_min <= m <= mr.m_max) &&
	(m in m_range(mr,l))
end

function _in((l,m)::Tuple{<:Integer,<:Integer},mr::LM)
	(abs(m)<=l) && (mr.l_min <= l <= mr.l_max) && (mr.m_min <= m <= mr.m_max) &&
	(l in l_range(mr,m))
end

function _in((l₂,l₁)::Tuple{<:Integer,<:Integer},mr::L₂L₁)
	(mr.l₁_min <= l₁ <= mr.l₁_max) &&
	(mr.l₂_min <= l₂ <= mr.l₂_max) &&
	(l₂ in l₂_range(mr,l₁))
end

Base.in(T::Tuple,mr::ModeRange) = _in(T,mr)

Base.last(mr::LM) = (last(l_range(mr,mr.m_max)),mr.m_max)
Base.last(mr::ML) = (mr.l_max,last(m_range(mr,mr.m_max)))
Base.last(mr::L₂L₁) = (last(l₂_range(mr,mr.l₁_max)), mr.l₁_max)

function modeindex(mr::LM,l::Integer,m::Integer)
	# Check if the l and m supplied correspond to valid modes
	check_if_non_negative(l)
	check_if_valid_mode(l,m)
	check_if_mode_present(mr,l,m)

	Nskip = 0

	l_min,l_max = mr.l_min,mr.l_max
	m_min,m_max = mr.m_min,mr.m_max

	# Nskip = sum(length(l_range(mr,m_i)) fot m_i=mr.m_min:m-1)

	up = min(m - 1, -1)
	lo = max(m_min, -l_min)
	if up >= lo
		Nskip += (1 + l_max - l_min)*(1 - lo + up)
	end

	up = min(m - 1, l_min)
	lo = max(m_min, 0)
	if up >= lo
		Nskip += (1 + l_max - l_min)*(1 - lo + up)
	end

	up = min(m - 1, m_max, l_max)
	lo = max(m_min, l_min + 1)
	if up >= lo
		Nskip += div((-1 + lo - up)*(-2*(1 + l_max) + lo + up),2)
	end

	up = min(m - 1, -1, -l_min - 1)
	lo = max(m_min, -l_max)
	if up >= lo
		Nskip += div(-((-1 + lo - up)*(2 + 2l_max + lo + up)),2)
	end

	ind_l = l - first(l_range(mr,m)) + 1
	Nskip + ind_l
end

function modeindex(mr::LM,l::AbstractUnitRange{<:Integer},m::Integer)
	last(l) in l_range(mr,m) || throw(ModeMissingError(last(l),m,mr))
	start = modeindex(mr,first(l),m)
	start:start + length(l) - 1
end

modeindex(mr::LM,::Colon,m::Integer) = modeindex(mr,l_range(mr,m),m)

function modeindex(mr::ML,l::Integer,m::Integer)
	# Check if the l and m supplied correspond to valid modes
	check_if_non_negative(l)
	check_if_valid_mode(l,m)
	check_if_mode_present(mr,l,m)

	Nskip = 0

	l_min,l_max = mr.l_min,mr.l_max
	m_min,m_max = mr.m_min,mr.m_max

	# Nskip = sum(length(m_range(mr,l_i)) for l_i=mr.l_min:l-1)

	up = min(-m_min, m_max, l - 1)
	lo = l_min
	if up >= lo
		Nskip += -(-1 + lo - up)*(1 + lo + up)
	end

	up = min(-m_min, l - 1)
	lo = max(l_min, m_max + 1, -m_max)
	if up >= lo
		Nskip += div(-((-1 + lo - up)*(2 + 2m_max + lo + up)),2)
	end

	up = min(m_max, l - 1)
	lo = max(l_min, -m_min + 1, m_min)
	if up >= lo
		Nskip += div(-((-1 + lo - up)*(2 - 2m_min + lo + up)),2)
	end

	up = l - 1
	lo = max(-m_min + 1, m_max + 1, l_min)
	if up >= lo
		Nskip += (1 + m_max - m_min)*(1 - lo + up)
	end

	ind_m = m - first(m_range(mr,l)) + 1
	Nskip + ind_m
end

function modeindex(mr::ML,l::Integer,m::AbstractUnitRange{<:Integer})
	last(m) in m_range(mr,l) || throw(ModeMissingError(l,last(m),mr))
	start = modeindex(mr,l,first(m))
	start:start + length(m) - 1
end

modeindex(mr::ML,l::Integer,::Colon) = modeindex(mr,l,m_range(mr,l))

function modeindex(mr::L₂L₁,l₂::Integer,l₁::Integer)
	# Check if l₁ and l₂ supplied correspond to valid modes
	map(check_if_non_negative,(l₁,l₂))
	check_if_mode_present(mr,l₂,l₁)

	Nskip = 0

	l_min,l_max,l₂_min,l₂_max = mr.l₁_min,mr.l₁_max,mr.l₂_min,mr.l₂_max
	Δl_max = mr.Δl_max

	up = min(l₁ - 1, l₂_max - Δl_max)
	lo = max(l_min, l₂_min + Δl_max)
	if up >= lo
		Nskip += (1 + 2Δl_max)*(1 - lo + up)
	end

	up = min(l₁ - 1, l₂_max + Δl_max)
	lo = max(l_min, l₂_min + Δl_max, l₂_max - Δl_max + 1)
	if up >= lo
		Nskip += div((-1 + lo - up)*(-2*(1 + Δl_max + l₂_max) + lo + up),2)
	end

	up = min(l₁ - 1, l₂_min + Δl_max - 1, l₂_max - Δl_max)
	lo = max(l_min, l₂_min - Δl_max)
	if up >= lo
		Nskip += div(-((-1 + lo - up)*(2 + 2Δl_max - 2l₂_min + lo + up)),2)
	end

	up = min(l₁ - 1, l₂_min + Δl_max - 1)
	lo = max(l_min, l₂_max - Δl_max + 1)
	if up >= lo
		Nskip += (1 + l₂_max - l₂_min)*(1 - lo + up)
	end

	ind_l₂ = l₂ - first(l₂_range(mr,l₁)) + 1
	Nskip + ind_l₂
end

function modeindex(mr::L₂L₁,l₂::AbstractUnitRange{<:Integer},l₁::Integer)
	last(l₂) in l₂_range(mr,l₁) || throw(ModeMissingError(last(l₂),l₁,mr))
	start = modeindex(mr,first(l₂),l₁)
	start:start + length(l₂) - 1
end

modeindex(mr::L₂L₁,::Colon,l₁::Integer) = modeindex(mr,l₂_range(mr,l₁),l₁)

modeindex(mr::ModeRange,T::Tuple{<:Any,<:Any}) = modeindex(mr,T...)
modeindex(mr::ModeRange,T::Vararg{<:Any,2}) = modeindex(mr,T...)
modeindex(mr::ModeRange,::Colon,::Colon) = Base.OneTo(length(mr))

# modeindex with two ranges works only if the values are contiguously stored in the iterator
# If the values are contiguous then it returns modeindex(l_min,m_min):modeindex(l_max,m_max)
# All (l,m) values need to be present in the iterator
# Some invalid values might be caught by internal modeindex calls
function modeindex(mr::ModeRange,l_range::AbstractUnitRange,m_range::AbstractUnitRange)
	m_min,m_max = extrema(m_range)
	l_min,l_max = extrema(l_range)
	firstind = modeindex(mr,l_min,m_min)
	lastind = modeindex(mr,l_max,m_max)

	if (lastind - firstind + 1) != length(l_range)*length(m_range)
		throw(NonContiguousError())
	end
	firstind:lastind
end

# Return the range of indices of m containing the entirety of mpart
# If mpart is contiguously stored in m then collect(m)[modeindex(mr,mpart)] == collect(mpart)
function modeindex(mr::T,mpart::T) where {T<:ModeRange}
	mode_start = first(mpart)
	mode_end = last(mpart)
	modeindex(mr,mode_start):modeindex(mr,mode_end)
end

# Display the iterators

function Base.show(io::IO, mr::SHModeRange)
	print(io,"(l=",l_range(mr),",m=",m_range(mr),")")
end

function Base.show(io::IO, mr::L₂L₁)
	print(io,"(l₁=",l₁_range(mr),",Δl_max=",mr.Δl_max,
		",l₂=",l₂_range(mr),")")
end

function Base.show(io::IO, ::MIME"text/plain", mr::LM)
	println("Spherical harmonic modes with l increasing faster than m")
	print(io,"(l_min = ",mr.l_min,", l_max = ",mr.l_max,
		", m_min = ",mr.m_min,", m_max = ",mr.m_max,")")
end

function Base.show(io::IO, ::MIME"text/plain", mr::ML)
	println("Spherical harmonic modes with m increasing faster than l")
	print(io,"(l_min = ",mr.l_min,", l_max = ",mr.l_max,
		", m_min = ",mr.m_min,", m_max = ",mr.m_max,")")
end

function Base.show(io::IO, ::MIME"text/plain", mr::L₂L₁)
	println("Spherical harmonic modes (l₂,l₁) where |l₁-Δl| ⩽ l₂ ⩽ l₁+Δl "*
		"for 0 ⩽ Δl ⩽ Δl_max, l₁_min ⩽ l₁ ⩽ l₁_max, and l₂_min ⩽ l₂ ⩽ l₂_max")
	print(io,"(",mr.l₂_min," ⩽ l₂ ⩽ ",mr.l₂_max," and ",
		mr.l₁_min," ⩽ l₁ ⩽ ",mr.l₁_max,", with Δl_max = ",mr.Δl_max,")")
end

end # module

ext: jl


text:
# Example Julia to export EndNote/BibTex citations from a PMID or a list of
# PMIDs. If you need to search Entrez/PubMed and save the results as citations,
# refer to Examples/pubmed_search_and_save.jl
# Date: September 8, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using BioMedQuery.Processes

#************************ LOCALS TO CONFIGURE!!!! **************************
results_dir = "./results"

 if !isdir(results_dir)
     mkdir(results_dir)
 end

email= ENV["NCBI_EMAIL"] #This is an enviroment variable that you need to setup
citation_type="endnote"
pmid = 11748933
output_file1=results_dir*"/11748933.enw"

pmid_list = [24008025, 24170597]
output_file2 = results_dir* "/pmid_list.enw"

#***************************************************************************

export_citation(email, pmid, citation_type, output_file1)
export_citation(email, pmid_list, citation_type, output_file2)


text:
module PLEXOSUtils

import EzXML: Document, eachelement, namespace, Node, nodecontent, parsexml
import InfoZIP: Archive, open_zip

export open_plexoszip, PLEXOSSolutionDataset, PLEXOSSolutionDatasetSummary

struct PLEXOSTable

    name::String
    fieldname::Symbol
    fieldtype::Symbol
    loadorder::Int
    identifier::Union{String,Nothing}
    timestampfield::Union{Symbol,Nothing}
    indexoffset::Int

    PLEXOSTable(name::String, fieldname::Symbol, fieldtype::Symbol,
                loadorder::Int,
                identifier::Union{String,Nothing}=nothing,
                timestampfield::Union{Symbol,Nothing}=nothing,
                indexoffset::Int=0
    ) = new(name, fieldname, fieldtype, loadorder,
            identifier, timestampfield, indexoffset)

end

plexostables = [

    PLEXOSTable("t_config", :configs, :PLEXOSConfig, 1),
    PLEXOSTable("t_unit", :units, :PLEXOSUnit, 1, "unit_id", nothing, 1),
    PLEXOSTable("t_timeslice", :timeslices, :PLEXOSTimeslice, 1, "timeslice_id", nothing, 1),
    PLEXOSTable("t_model", :models, :PLEXOSModel, 1, "model_id"),

    PLEXOSTable("t_sample", :samples, :PLEXOSSample, 1, "sample_id", nothing, 4),
    PLEXOSTable("t_sample_weight", :sampleweights, :PLEXOSSampleWeight, 2),

    PLEXOSTable("t_class_group", :classgroups, :PLEXOSClassGroup, 1, "class_group_id"),
    PLEXOSTable("t_class", :classes, :PLEXOSClass, 2, "class_id"),
    PLEXOSTable("t_category", :categories, :PLEXOSCategory, 3, "category_id"),

    PLEXOSTable("t_attribute", :attributes, :PLEXOSAttribute, 3, "attribute_id"),
    PLEXOSTable("t_collection", :collections, :PLEXOSCollection, 3, "collection_id"),
    PLEXOSTable("t_property", :properties, :PLEXOSProperty, 4, "property_id"),

    PLEXOSTable("t_object", :objects, :PLEXOSObject, 4, "object_id"),
    PLEXOSTable("t_membership", :memberships, :PLEXOSMembership, 5, "membership_id"),
    PLEXOSTable("t_attribute_data", :attribute_data, :PLEXOSAttributeData, 5),

    PLEXOSTable("t_custom_column", :customcolumns, :PLEXOSCustomColumn, 3, "column_id"),
    PLEXOSTable("t_memo_object", :memos, :PLEXOSMemoObject, 5),

    PLEXOSTable("t_period_0", :intervals, :PLEXOSPeriod0, 1, "interval_id", :datetime),
    PLEXOSTable("t_period_1", :days, :PLEXOSPeriod1, 1, "day_id", :date),
    PLEXOSTable("t_period_2", :weeks, :PLEXOSPeriod2, 1, "week_id", :week_ending),
    PLEXOSTable("t_period_3", :months, :PLEXOSPeriod3, 1, "month_id", :month_beginning),
    PLEXOSTable("t_period_4", :years, :PLEXOSPeriod4, 1, "fiscal_year_id", :year_ending),
    PLEXOSTable("t_period_6", :hours, :PLEXOSPeriod6, 1, "hour_id", :datetime),
    PLEXOSTable("t_period_7", :quarters, :PLEXOSPeriod7, 1, "quarter_id", :quarter_beginning),

    PLEXOSTable("t_phase_1", :lt, :PLEXOSPhase1, 2),
    PLEXOSTable("t_phase_2", :pasa, :PLEXOSPhase2, 2),
    PLEXOSTable("t_phase_3", :mt, :PLEXOSPhase3, 2),
    PLEXOSTable("t_phase_4", :st, :PLEXOSPhase4, 2),

    PLEXOSTable("t_key", :keys, :PLEXOSKey, 6, "key_id"),
    PLEXOSTable("t_key_index", :keyindices, :PLEXOSKeyIndex, 7)

]

plexostables_lookup = Dict(x.name => x for x in plexostables)

include("types.jl")
include("PLEXOSSolutionDatasetSummary.jl")
include("PLEXOSSolutionDataset.jl")
include("utils.jl")

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

abstract type SolutionMethod end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	configure_solver(solver::String, solver_settings_path::String)

This method returns a solver-specific MathOptInterface OptimizerWithAttributes optimizer instance to be used in the GenX.generate\_model() method.

The "solver" argument is a string which specifies the solver to be used. It is not case sensitive.
Currently supported solvers include: "Gurobi", "CPLEX", "Clp", "Cbc", or "SCIP"

The "solver\_settings\_path" argument is a string which specifies the path to the directory that contains the settings YAML file for the specified solver.

"""
function configure_solver(solver::String, solver_settings_path::String)

    println("Configuring Solver")

    solver = lowercase(solver)

    if solver == "gurobi"
        gurobi_settings_path = joinpath(solver_settings_path, "gurobi_settings.yml")
        OPTIMIZER = configure_gurobi(gurobi_settings_path)
    elseif solver == "cplex"
        cplex_settings_path = joinpath(solver_settings_path, "cplex_settings.yml")
        OPTIMIZER = configure_cplex(cplex_settings_path)
    elseif solver == "clp"
        clp_settings_path = joinpath(solver_settings_path, "clp_settings.yml")
        OPTIMIZER = configure_clp(clp_settings_path)
    elseif solver == "cbc"
        cbc_settings_path = joinpath(solver_settings_path, "cbc_settings.yml")
        OPTIMIZER = configure_cbc(cbc_settings_path)
    elseif solver == "scip"
        scip_settings_path = joinpath(solver_settings_path, "scip_settings.yml")
        OPTIMIZER = configure_scip(scip_settings_path)
    elseif solver == "glpk"
        glpk_settings_path = joinpath(solver_settings_path, "glpk_settings.yml")
        OPTIMIZER = configure_glpk(glpk_settings_path)
    end

    return OPTIMIZER
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text: for i = 1 : (39 + "three") end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test
using Nord, ColorTypes, FixedPointNumbers

@test Nord.nord0 == RGB{N0f8}(46/255, 52/255, 64/255)


text:
using NiLang
using Test

@testset "⊕, ⊖" begin
    x, y = 0.2, 0.5
    @instr x ⊕ y
    @test x == 0.7 && y==0.5
end

@testset "*, /" begin
    x, y, out = 2.0, 2.0, 1.0
    @instr out += x * y
    @test x == 2.0 && y == 2.0 && out == 5.0
    x, y, out = 2.0, 2.0, 1.0
    @instr out += x / y
    @test x == 2.0 && y == 2.0 && out == 2.0
end

@testset "XOR, SWAP" begin
    x, y = 1, 2
    @instr XOR(x, y)
    @test x == 3 && y == 2
    x, y = 1, 2
    @instr SWAP(x, y)
    @test x == 2 && y == 1
end

@testset "CONJ, NEG" begin
    x = 0.3 + 2im
    @instr CONJ(x)
    @test x == 0.3-2im
    @instr NEG(x)
    @test x == -0.3+2im
    @test check_inv(CONJ, (x,))
    @test check_inv(NEG, (x,))
end

@testset "ROT" begin
    x, y, θ = 0.0, 1.0, π
    #@instr ROT(x, y, θ)
    @test check_inv(ROT, (x, y, θ); verbose=true)
    @test check_inv(IROT, (x, y, θ); verbose=true)
end


text:
# # [Negative acknowledgement (NACK)](@id nack)

using Reagents
using Reagents: WithNack, Return
using Test

# Reagents.jl provides so-called *negative acknowledgement* (NACK) reagent
# [`Reagents.WithNack`](@ref) which is taken from [Concurrent
# ML](https://en.wikipedia.org/wiki/Concurrent_ML).  This is useful for writing
# "client-server" style of code where the client can abort the request.
#
# !!! note
#
#     Concurrent ML provides composable synchronizable operations called
#     *events* which are similar to [`Reagent`](@ref Reagents.Reagent).  Turon
#     (2012) discussed influences of Concurrent ML on reagents.
#
#     See also:
#     * [Concurrent ML's manual on `withNack`](http://cml.cs.uchicago.edu/pages/cml.html#SIG:CML.withNack:VAL)
#     * [Racket's manual on `nack-guard-evt`](https://docs.racket-lang.org/reference/sync.html#%28def._%28%28quote._~23~25kernel%29._nack-guard-evt%29%29)
#
# ## How it works
#
# Let us set up a demo.

function nack_demo()
    #=
    [`Reagents.WithNack`](@ref) has non-trivial effect only when used inside the
    choice combinator which possibly blocking branches.  Thus, to selectively
    trigger two branches in the choice combinator, we create two channels:
    =#
    s1, r1 = Reagents.channel()
    s2, r2 = Reagents.channel()
    #=
    To receive the negative acknowledgement, we craete one more channel:
    =#
    send_gotnack, receive_gotnack = Reagents.channel()
    #=
    The first branch `br1` (below) uses [`Reagents.WithNack`](@ref).  It passes
    the negative acknowledgement reagent `nack` to the user-defined function
    (the `do` block; it returns a reagent).  The reagent `nack` blocks until
    this branch `br1` is cancelled (i.e., another branch of `|` is chosen).
    =#
    br1 = WithNack() do nack
        @async (nack ∘ send_gotnack)(:gotnack)
        return r1
    end
    #=
    We just use a channel endpoint for another branch:
    =#
    br2 = r2
    #=
    These two reagents are composed with the choice combinator `|`:
    =#
    choice = br1 | br2
    #=
    Returning the reagents so that they can be invoked differently for trying
    differnt scenarios:
    =#
    return (; choice, s1, s2, receive_gotnack)
end

# ### Scenario 1: `nack` is triggered

function test_nack_demo_1()
    (; choice, s2, receive_gotnack) = nack_demo()
    @sync begin
        #=
        Let us choose the second branch `br2` which does *not* include
        `WithNack`:
        =#
        @async s2(222)
        @test choice() == 222
    end
    #=
    Since the branch `br1` with `WithNack` is not chosen, we get the negative
    acknowledgement:
    =#
    @test receive_gotnack() == :gotnack
end

# ### Scenario 2: `nack` is not triggered

function test_nack_demo_2()
    (; choice, s1, receive_gotnack) = nack_demo()
    @sync begin
        #=
        This time, we choose the first branch `br1` which includes `WithNack`:
        =#
        @async s1(111)
        @test choice() == 111
    end
    #=
    Since we chose the `WithNack`'s branch, `nack` is not triggered this time:
    =#
    @test Reagents.trysync!(receive_gotnack) === nothing
end

# ## Client-server pattern
#
# `WithNack` is useful for writing "client-server" pattern. As an example, we'll
# create an in-process "server" that issues unique IDs. That is to say, we'd
# like to have the following API:

function test_unique_id_provider_api()
    with_unique_id_provider() do unique_id
        @test unique_id() == 0
        @test unique_id() == 1
    end
end

# Here, `unique_id` is a reagent for communicating with a server created in
# `with_unique_id_provider`.
#
# ### `unique_id_provider!`
#
# Let us start from the event loop of the server.  The server listens to ID
# requests from `request_receive` and a shutdown request from
# `shutdown_receive`.

function unique_id_provider!(request_receive, shutdown_receive)
    #=
    It keeps the current available ID as its local variable:
    =#
    id = 0
    while true
        #=
        First, the server listens to both `request_receive` and
        `shutdown_receive`. The latter returns `nothing` upon shutdown request.
        =#
        receive_request_or_shutdown = request_receive | shutdown_receive
        #=
        When the `shutdown_receive` reagent is chosen (i.e., the reaction result
        is `nothing`), the short-circuting `@something` evaluates the `break`
        statement so that the server exits the loop:
        =#
        (; reply, abort) = @something(receive_request_or_shutdown(), break)
        #=
        The client (see below) sends `reply` and `abort` channel endpoints.  The
        server tries to send the ID with `Return(id) ⨟ reply` while also
        listening to the abort (NACK) and shutdown requests:
        =#
        try_reply = (
            (Return(id) ⨟ reply ⨟ Return(true)) |  # try sending the id
            (abort ⨟ Return(false)) |              # or wait for the abort (NACK)
            shutdown_receive                       # or wait for shutdown
        )
        #=
        The server only increments the ID when the client received the ID.
        =#
        if @something(try_reply(), break)
            id += 1
        end
    end  # while true
end  # function unique_id_provider!

# (For an ID server, this property is probably not required. But consider, e.g.,
# a lock server, where it is important to know that the client received the
# reply.)
#
# ### `with_unique_id_provider`
#
# The channels connecting the server and client are set up in the function
# below. The client API can be invoked inside the function `f` passed as the
# argument:

function with_unique_id_provider(f)
    request_send, request_receive = Reagents.channel()
    shutdown_send, shutdown_receive = Reagents.channel(Nothing)

    #=
    For each request, the client creates the channel (`reply`) for receiving the
    ID and also the negative acknowledgement reaagent `abort` for communicating
    that the request is aborted:
    =#
    unique_id = WithNack() do abort
        reply, receive = Reagents.channel(Int, Nothing)
        request_send((; reply, abort))
        return receive
    end

    #=
    Finally, we start the server in a task and execute the client's code `f`:
    =#
    @sync begin
        @async unique_id_provider!(request_receive, shutdown_receive)
        try
            f(unique_id)
        finally
            shutdown_send()
        end
    end
end

# ### Testing the ID server

function test_unique_id_provider()
    with_unique_id_provider() do unique_id
        #=
        When used alone, `unique_id` simply sends a request and wait for a
        reply from the ID server:
        =#
        @test unique_id() == 0
        @test unique_id() == 1

        #=
        Demonstrating the behavior of aborting the request is a bit more
        involved. First, we create a task that tries to send the "cancellation"
        request via a channel:
        =#
        send, receive = Reagents.channel(Nothing)
        canceller = @task send()
        yield(canceller)
        #=
        Since we don't know when `send()` will be invoked, we'll try it in a
        loop. The variable `prev` keeps track of the last `id` issued by the
        server:
        =#
        prev = unique_id()
        while true
            #=
            Then invoke `unique_id` and `receive` together. If this reaction
            takes choose the branch of `receive` ("cancellation"), it returns a
            `nothing`:
            =#
            ans = (unique_id | receive)()
            if ans === nothing
                #=
                Here, we have attempted to invoke the `unique_id` reagent but it
                was aborted by another reagent `receive`.  Since this triggers
                the nack reagent `abort`, this reaction did not update the
                server's state (the variable `id`). So, the next call to
                `unique_id` should increment the ID only by one:
                =#
                @test unique_id() == prev + 1
                break
                #=
                If `receive` was not tirggered, we keep the id `ans` so that it
                can be used in the next iteration:
                =#
            else
                prev = ans::Int
            end
        end
    end
end


text:
function evaluate_individual_view!(p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig,results::Union{Array{Float64,3}, Array{Float64,4}})

    grn = initialise_grn(p_float,p_int,grn_config)

    environment = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,2 + size(p_float,1)))

    zygote = vcat([1.0,4.0], ones(Float64,size(p_float,1)))

    for i in 1:size(zygote,1)
        if i in 2 .+ grn.rec_id[1:3]
            zygote[i] = 0.
        end
    end

    environment[Int(round(cell_config.grid_size/2)),Int(round(cell_config.grid_size/2)),:] = zygote

    J = create_differential_adhesion_matrix(grn,grn_config)

    dc = zeros(Float64,size(p_float,1))

    Ax = zeros(Float64,cell_config.grid_size,cell_config.grid_size)
    Ay = zeros(Float64,cell_config.grid_size,cell_config.grid_size)

    diffusion_constants = create_diffusion_constants(grn,grn_config,exp_config)

    run_life_view!(results,environment,grn,J,diffusion_constants,
            Ax,Ay,dc,cell_config,exp_config,grn_config)

end

function evaluate_individual_view!(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig,results::Array{Float64,3})

    grn = initialise_grn(p_float,p_int,grn_config)

    environment = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,2 + size(p_float,1)))

    zygote = vcat([1.0,4.0], ones(Float64,size(p_float,1)))

    for i in 1:size(zygote,1)
        if i in 2 .+ grn.rec_id[1:3]
            zygote[i] = 0.
        end
    end

    environment[Int(round(cell_config.grid_size/2)),Int(round(cell_config.grid_size/2)),:] = zygote

    J = create_differential_adhesion_matrix(grn,grn_config)

    dc = zeros(Float64,size(p_float,1))

    Ax = zeros(Float64,cell_config.grid_size,cell_config.grid_size)
    Ay = zeros(Float64,cell_config.grid_size,cell_config.grid_size)

    diffusion_constants = create_diffusion_constants(grn,grn_config,exp_config)

    run_life_view!(damage_times,damage_percentage,results,environment,grn,J,diffusion_constants,
            Ax,Ay,dc,cell_config,exp_config,grn_config)

end

function evaluate_individual_view!(damage_times::Vector{Int64},damage_percentage::Float64,
    p_float::Array{Float64,2},p_int::Vector{Int64},cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig,energies::Array{Float64,1},results::Array{Float64,3})

    grn = initialise_grn(p_float,p_int,grn_config)

    environment = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,2 + size(p_float,1)))

    zygote = vcat([1.0,4.0], ones(Float64,size(p_float,1)))

    for i in 1:size(zygote,1)
        if i in 2 .+ grn.rec_id[1:3]
            zygote[i] = 0.
        end
    end

    environment[Int(round(cell_config.grid_size/2)),Int(round(cell_config.grid_size/2)),:] = zygote

    J = create_differential_adhesion_matrix(grn,grn_config)

    dc = zeros(Float64,size(p_float,1))

    Ax = zeros(Float64,cell_config.grid_size,cell_config.grid_size)
    Ay = zeros(Float64,cell_config.grid_size,cell_config.grid_size)

    diffusion_constants = create_diffusion_constants(grn,grn_config,exp_config)

    run_life_view!(damage_times,damage_percentage,energies,results,environment,grn,J,diffusion_constants,
            Ax,Ay,dc,cell_config,exp_config,grn_config)

end

function evaluate_individual_fitness(p_float::Array{Float64,2},p_int::Vector{Int64},evo_config::EMAPConfig,cell_config::CellularConfig,
    exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    grn = initialise_grn(p_float,p_int,grn_config)

    environment = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,2 + size(p_float,1)))

    zygote = vcat([1.0,4.0], ones(Float64,size(p_float,1)))

    for i in 1:size(zygote,1)
        if i in 2 .+ grn.rec_id[1:3]
            zygote[i] = 0.
        end
    end

    environment[Int(round(cell_config.grid_size/2)),Int(round(cell_config.grid_size/2)),:] = zygote

    J = create_differential_adhesion_matrix(grn,grn_config)

    diffusion_constants = create_diffusion_constants(grn,grn_config,exp_config)

    Ax = zeros(Float64,cell_config.grid_size,cell_config.grid_size)
    Ay = zeros(Float64,cell_config.grid_size,cell_config.grid_size)

    dc = zeros(Float64,size(p_float,1))

    fingerprints = zeros(Int,16)
    fingerprint_fitnesses = zeros(Float64,16)

    n_zero = 0
    center = Int(round(cell_config.grid_size/2))

    for trial in 1:evo_config.n_evaluation

        environment .= zero(Float64)

        environment[center,center,:] = zygote

        target_cell_areas = Dict{Float64,Float64}()
        n_div = Dict{Float64,Integer}()

        all_cells::Array{Float64,1} = [1.]
        n_div[1.] = 0
        target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

        run_life_fitness!(evo_config.time_steps,environment,all_cells,target_cell_areas,n_div,
                grn,J,diffusion_constants,
                Ax,Ay,dc,cell_config,exp_config,grn_config)

        f,fp = evaluate_fitness(environment[:,:,2],evo_config)

        for i in 1:16
            if fp == i
                fingerprints[i] += 1
                fingerprint_fitnesses[i] += f
            end
        end

    end

    winner = argmax(fingerprints)
    fitness = fingerprint_fitnesses[winner]/fingerprints[winner]

    return fitness,winner

end

function evaluate_individual_fitness(p_float::Array{Float64,2},p_int::Vector{Int64},evo_config::POOLConfig,
    cell_config::CellularConfig,exp_config::ExperimentalConfig,grn_config::FGRNConfig)

    grn = initialise_grn(p_float,p_int,grn_config)

    environment = zeros(Float64,(cell_config.grid_size,cell_config.grid_size,2 + size(p_float,1)))

    zygote = vcat([1.0,4.0], ones(Float64,size(p_float,1)))

    for i in 1:size(zygote,1)
        if i in 2 .+ grn.rec_id[1:3]
            zygote[i] = 0.
        end
    end

    environment[Int(round(cell_config.grid_size/2)),Int(round(cell_config.grid_size/2)),:] = zygote

    J = create_differential_adhesion_matrix(grn,grn_config)

    diffusion_constants = create_diffusion_constants(grn,grn_config,exp_config)

    Ax = zeros(Float64,cell_config.grid_size,cell_config.grid_size)
    Ay = zeros(Float64,cell_config.grid_size,cell_config.grid_size)

    dc = zeros(Float64,size(p_float,1))

    fitnesses = zeros(evo_config.n_evaluation)

    n_zero = 0
    center = Int(round(cell_config.grid_size/2))

    for trial in 1:evo_config.n_evaluation

        environment .= zero(Float64)

        environment[center,center,:] = zygote

        target_cell_areas = Dict{Float64,Float64}()
        n_div = Dict{Float64,Integer}()

        all_cells::Array{Float64,1} = [1.]
        n_div[1.] = 0
        target_cell_areas[1.] = 2*0.25*cell_config.cell_size_upper_limit

        run_life_fitness!(evo_config.time_steps,environment,
                all_cells,target_cell_areas,n_div,
                grn,J,diffusion_constants,Ax,Ay,dc,cell_config,exp_config,grn_config)

        f = evaluate_fitness(environment[:,:,2],evo_config)

        fitnesses[trial] = f

    end

    return mean(fitnesses)

end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
export SampleList

mutable struct SampleList <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function SampleList(out, s, w; id=generateId(SampleList))
        @ensureVariables(out, s, w)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out] = self.interfaces[1] = associate!(Interface(self), out)
        self.i[:s] = self.interfaces[2] = associate!(Interface(self), s)
        self.i[:w] = self.interfaces[3] = associate!(Interface(self), w)

        return self
    end
end

slug(::Type{SampleList}) = "SampleList"

format(dist::ProbabilityDistribution{V, SampleList}) where V<:VariateType = "$(slug(SampleList))(s=$(format(dist.params[:s])),w=$(format(dist.params[:w])))"

ProbabilityDistribution(::Type{Univariate}, ::Type{SampleList}; s=[0.0], w=[1.0]) = ProbabilityDistribution{Univariate, SampleList}(Dict{Symbol, Any}(:s=>s, :w=>w))
ProbabilityDistribution(::Type{Multivariate}, ::Type{SampleList}; s=[[0.0]], w=[1.0]) = ProbabilityDistribution{Multivariate, SampleList}(Dict{Symbol, Any}(:s=>s, :w=>w))
ProbabilityDistribution(::Type{MatrixVariate}, ::Type{SampleList};s=[mat(0.0)], w=[1.0]) = ProbabilityDistribution{MatrixVariate, SampleList}(Dict{Symbol,Any}(:s=>s,:w=>w))

dims(dist::ProbabilityDistribution{<:VariateType, SampleList}) = size(dist.params[:s][1])

function vague(::Type{SampleList})
    n_samples = default_n_samples # Fixed number of samples

    return ProbabilityDistribution(Univariate, SampleList, s=rand(n_samples), w=ones(n_samples)/n_samples)
end

function vague(::Type{SampleList}, dims::Tuple)
    n_samples = default_n_samples # Fixed number of samples

    s_list = Vector(undef, n_samples)
    for n=1:n_samples
        s_list[n] = randn(dims)
    end

    return ProbabilityDistribution(Multivariate, SampleList, s=s_list, w=ones(n_samples)/n_samples)
end

unsafeMean(dist::ProbabilityDistribution{V, SampleList}) where V<:VariateType = sum(dist.params[:s].*dist.params[:w])

unsafeLogMean(dist::ProbabilityDistribution{Univariate, SampleList}) = sum(log.(dist.params[:s]).*dist.params[:w])

unsafeMeanLogMean(dist::ProbabilityDistribution{Univariate, SampleList}) = sum(dist.params[:s].*log.(dist.params[:s]).*dist.params[:w])

function unsafeLogMean(dist::ProbabilityDistribution{Multivariate, SampleList})
    sum = zeros(length(dist.params[:s][1]))
    for i=1:length(dist.params[:s])
        sum = sum .+ log.(dist.params[:s][i]).*dist.params[:w][i]
    end
    return sum
end

# Unbiased (co)variance estimates
function unsafeVar(dist::ProbabilityDistribution{Univariate, SampleList})
    samples = dist.params[:s]
    n_samples = length(samples)

    return (n_samples/(n_samples - 1))*sum(dist.params[:w].*(samples .- unsafeMean(dist)).^2)
end

unsafeCov(dist::ProbabilityDistribution{Univariate, SampleList}) = unsafeVar(dist)

unsafeVar(dist::ProbabilityDistribution{Multivariate, SampleList}) = diag(unsafeCov(dist))
function unsafeCov(dist::ProbabilityDistribution{Multivariate, SampleList})
    samples = dist.params[:s]
    weights = dist.params[:w]

    n_samples = length(samples)
    m = unsafeMean(dist)
    d = dims(dist)[1]
    tot = zeros(d, d)
    for i = 1:n_samples
        tot += (samples[i] .- m)*transpose(samples[i] .- m).*weights[i]
    end

    return (n_samples/(n_samples - 1)).*tot
end

function unsafeCov(dist::ProbabilityDistribution{MatrixVariate, SampleList})
    samples = dist.params[:s]
    weights = dist.params[:w]

    n_samples = length(samples)
    m = unsafeMean(dist)
    (d1, d2) = dims(dist)
    cov1 = zeros(d1, d1)
    cov2 = zeros(d2, d2)

    for i = 1:n_samples
        cov1 += ((samples[i] .- m))*transpose((samples[i] .- m)).*weights[i]
        cov2 += transpose((samples[i] .- m))*((samples[i] .- m)).*weights[i]
    end
    cov1 = (n_samples/(n_samples - 1)).*cov1
    cov2 = (n_samples/(n_samples - 1)).*cov2

    return kron(cov1, cov2)
end

unsafeMeanCov(dist::ProbabilityDistribution{<:VariateType, SampleList}) = (unsafeMean(dist), unsafeCov(dist))

function unsafeMirroredLogMean(dist::ProbabilityDistribution{Univariate, SampleList})
    all(0 .<= dist.params[:s] .< 1) || error("unsafeMirroredLogMean does not apply to variables outside of the range [0, 1]")

    return sum(log.(1 .- dist.params[:s]) .* dist.params[:w])
end

unsafeMeanVector(dist::ProbabilityDistribution{<:VariateType, SampleList}) = sum(dist.params[:s].*dist.params[:w])

isProper(dist::ProbabilityDistribution{<:VariateType, SampleList}) = abs(sum(dist.params[:w]) - 1) < 0.001

# prod of a pdf (or distribution) message and a SampleList message
# this function is capable to calculate entropy with SampleList messages in VMP setting
@symmetrical function prod!(
    x::ProbabilityDistribution{V}, # Includes function distributions
    y::ProbabilityDistribution{V, SampleList},
    z::ProbabilityDistribution{V, SampleList}=ProbabilityDistribution(V, SampleList, s=[0.0], w=[1.0])) where V<:VariateType

    # Suppose in the previous time step m1(pdf) and m2(pdf) messages collided.
    # The resulting collision m3 (sampleList) = m1*m2 is supposed to carry
    # the proposal (m1) and integrand (m2) distributions. m1 is the message from which
    # the samples are drawn. m2 is the message on which the samples are evaluated and
    # weights are calculated. In case Particle Filtering (BP), entropy will not be calculated
    # and in the first step there won't be any integrand information.
    if haskey(y.params, :logintegrand)
        # recall that we are calculating m3*m4. If m3 consists of integrand information
        # update it: new_integrand = m2*m3. This allows us to collide arbitrary number of beliefs
        # to approximate posterior and yet estimate the entropy.
        logIntegrand = (samples) -> y.params[:logintegrand](samples) .+ logPdf.([x], samples)
    else
        # If there is no integrand information before, set it to m4
        logIntegrand = (samples) -> logPdf.([x], samples)
    end

    samples = y.params[:s] # samples come from proposal (m1)
    n_samples = length(samples) # number of samples
    log_samples_x = logPdf.([x], samples) # evaluate samples in logm4, i.e. logm4(s)

    # Compute sample weights
    w_raw_x = exp.(log_samples_x) # m4(s)
    w_prod = w_raw_x.*y.params[:w] # update the weights of posterior w_unnormalized = m4(s)*w_prev
    weights = w_prod./sum(w_prod) # Normalize weights

    # Resample if required
    n_eff = 1/sum(weights.^2) # Effective number of particles
    if n_eff < n_samples/10
        samples = sample(samples, Weights(weights), n_samples) # Use StatsBase for resampling
        weights = ones(n_samples)./n_samples
    end

    # resulting posterior or message
    z.params[:w] = weights # set adjusted weights
    z.params[:s] = samples # samples are still coming from the same proposal
    z.params[:logintegrand] = logIntegrand # set integrand
    if haskey(y.params, :logproposal) && haskey(y.params, :unnormalizedweights)
        z.params[:unnormalizedweights] = w_raw_x.*y.params[:unnormalizedweights] # m4(s)*m2(s)
        logProposal = y.params[:logproposal] # m1
        z.params[:logproposal] = logProposal # m1
        # calculate entropy
        H_y = log(sum(w_raw_x.*y.params[:unnormalizedweights])) - log(n_samples) # log(sum_i(m4(s_i)*m2(s_i))/N)
        # -sum_i(w_i*log(m1(s_i)*m2(s_i)*m4(s_i)))
        H_x = -sum( weights.*(logProposal(samples) + log.(y.params[:unnormalizedweights]) + log_samples_x) )
        entropy = H_x + H_y
        z.params[:entropy] = entropy
    end

    return z
end

# Disambiguate beteen SampleList product and nonlinear product of any distribution with a Gaussian
# Following two definitions must be parameterized on separate Univariate and Multivariate types
@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate, <:Gaussian},
    y::ProbabilityDistribution{Univariate, SampleList})

    z = ProbabilityDistribution(Univariate, SampleList, s=[0.0], w=[1.0])

    return prod!(x, y, z) # Return a SampleList
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Multivariate, <:Gaussian},
    y::ProbabilityDistribution{Multivariate, SampleList})

    z = ProbabilityDistribution(Multivariate, SampleList, s=[[0.0]], w=[1.0])

    return prod!(x, y, z) # Return a SampleList
end

function sampleWeightsAndEntropy(x::ProbabilityDistribution{<:VariateType, <:Function}, y::ProbabilityDistribution)
    sampleWeightsAndEntropy(y, x)
end

function sampleWeightsAndEntropy(x::ProbabilityDistribution, y::ProbabilityDistribution)
    n_samples = default_n_samples # Number of samples is fixed
    samples = sample(x, n_samples)

    # Apply log-pdf functions to the samples
    log_samples_x = logPdf.([x], samples)
    log_samples_y = logPdf.([y], samples)

    # Extract the sample weights
    w_raw = exp.(log_samples_y) # Unnormalized weights
    w_sum = sum(w_raw)
    weights = w_raw./w_sum # Normalize the raw weights

    # Compute the separate contributions to the entropy
    H_y = log(w_sum) - log(n_samples)
    H_x = -sum( weights.*(log_samples_x + log_samples_y) )
    entropy = H_x + H_y

    # Inform next step about the proposal and integrand to be used in entropy calculation in smoothing
    logproposal = (samples) -> logPdf.([x], samples)
    logintegrand = (samples) -> logPdf.([y], samples)

    return (samples, weights, w_raw, logproposal, logintegrand, entropy)
end

# General product definition that returns a SampleList
function prod!(
    x::ProbabilityDistribution{V},
    y::ProbabilityDistribution{V},
    z::ProbabilityDistribution{V, SampleList} = ProbabilityDistribution(V, SampleList, s=[0.0], w=[1.0])) where {V<:VariateType}

    (samples, weights, unnormalizedweights, logproposal, logintegrand, entropy) = sampleWeightsAndEntropy(x, y)

    z.params[:s] = samples
    z.params[:w] = weights
    z.params[:unnormalizedweights] = unnormalizedweights
    z.params[:logproposal] = logproposal
    z.params[:logintegrand] = logintegrand
    z.params[:entropy] = entropy

    return z
end

# Specialized product definition that accounts for the Categorical VariateType
@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate, Categorical},
    y::ProbabilityDistribution{Multivariate, Function},
    z::ProbabilityDistribution{Univariate, SampleList} = ProbabilityDistribution(Univariate, SampleList, s=[0.0], w=[1.0]))

    (samples, weights, unnormalizedweights, logproposal, logintegrand, entropy) = sampleWeightsAndEntropy(x, y)

    z.params[:s] = samples
    z.params[:w] = weights
    z.params[:unnormalizedweights] = unnormalizedweights
    z.params[:logproposal] = logproposal
    z.params[:logintegrand] = logintegrand
    z.params[:entropy] = entropy

    return z
end

# Bootstrap samples
function bootstrap(dist_mean::ProbabilityDistribution{Univariate, SampleList}, dist_var::ProbabilityDistribution{Univariate, PointMass})
    s_m = dist_mean.params[:s] # Samples representing the mean
    N = length(s_m)
    v = dist_var.params[:m] # Fixed variance

    return sqrt(v)*randn(N) .+ s_m # New samples
end

function bootstrap(dist_mean::ProbabilityDistribution{Multivariate, SampleList}, dist_var::ProbabilityDistribution{MatrixVariate, PointMass})
    s_m = dist_mean.params[:s] # Samples representing the mean
    N = length(s_m)
    V = dist_var.params[:m] # Fixed variance
    U = (cholesky(default_cholesky_mode, V)).U # Precompute Cholesky

    return [U' *randn(dims(dist_mean)) + s_m[i] for i in 1:N] # New samples
end

function bootstrap(dist_mean::ProbabilityDistribution{Univariate, <:Gaussian}, dist_var::ProbabilityDistribution{Univariate, SampleList})
    s_v = dist_var.params[:s] # Samples representing the variance
    N = length(s_v)
    (m, v) = unsafeMeanCov(dist_mean)
    s_u = sqrt.(s_v .+ v) # Standard deviation for each variance sample

    return s_u.*randn(N) .+ m # New samples
end

function bootstrap(dist_mean::ProbabilityDistribution{Multivariate, <:Gaussian}, dist_var::ProbabilityDistribution{MatrixVariate, SampleList})
    s_V = dist_var.params[:s] # Samples representing the covariance
    N = length(s_V)
    (m, V) = unsafeMeanCov(dist_mean)
    s_U = [(cholesky(default_cholesky_mode, s_V[i] + V)).U for i in 1:N] # Precompute Cholesky for each covariance sample; this can be expensive

    return [s_U[i]' *randn(dims(dist_mean)) + m for i in 1:N] # New samples
end

#Sampling from a distribution in ForneyLab returns equally weighted samples from the distribution
#To retain the unified standard procedure, we allow sampling from SampleList not through directly returning
#sample and weight parameters but drawing samples from the existing list of samples according to weights.
# Inverse-transform sampling
sample(dist::ProbabilityDistribution{<:VariateType, SampleList}) = sample(dist.params[:s], Weights(dist.params[:w]))

# Differential entropy for SampleList
function differentialEntropy(dist::ProbabilityDistribution{<:VariateType, SampleList})
    haskey(dist.params, :entropy) || error("Missing entropy for SampleList; quantity is requested but not computed")

    return dist.params[:entropy] # Entropy is pre-computed during computation of the marginal
end


text:
# Takes full dataset and thins by a factor of 100:

using JLD
@load "tyc3559_mcmc_norv_fixedkernel_v03.jld"

nthin = 100
par_mcmc_thin = zeros(nwalkers,nthin,nparam)
log_like_mcmc_thin =  zeros(nwalkers,nthin)

for i=1:nwalkers
  for j=1:nthin
    par_mcmc_thin[i,j,:]=par_mcmc[i,(j-1)*100+1,:]
    log_like_mcmc_thin[i,j,:]=log_like_mcmc[i,(j-1)*100+1]
  end
end

save("tyc3559_mcmc_norv_fixedkernel_v03_thin.jld","par_mcmc",par_mcmc_thin,"log_like_mcmc",log_like_mcmc_thin,"nwalkers",nwalkers,"nsteps",nthin,"nparam",nparam,"nvary",nvary,"ivary",ivary,"time_phot",time_phot,"flux",flux,"nlast",nlast)


text:
using Oceananigans.Utils: with_tracers
using Oceananigans.Operators: assumed_field_location

@inline zeroforcing(args...) = 0

"""
    regularize_forcing(forcing, field_name, model_field_names)

"Regularizes" or "adds information" to
user-defined forcing objects that are passed to the `IncompressibleModel`
constructor. `regularize_forcing` is called inside `model_forcing`.
We need `regularize_forcing` because it is only until `model_forcing` is called that
the *names* of various forcing fields are available. The `field_name` can be used to infer
the location at which the forcing is applied, or to add a field dependency
to a special forcing object, as for `Relxation`.
"""
regularize_forcing(forcing, field_name, model_field_names) = forcing # fallback

""" Wrap `forcing` in a `ContinuousForcing` at the location of `field_name`. """
function regularize_forcing(forcing::Function, field_name, model_field_names)
    X, Y, Z = assumed_field_location(field_name)
    return ContinuousForcing{X, Y, Z}(forcing)
end

regularize_forcing(::Nothing, field_name, model_field_names) = zeroforcing

"""
    model_forcing(; u=zeroforcing, v=zeroforcing, w=zeroforcing, tracer_forcings...)

Return a named tuple of forcing functions for each solution field, wrapping
forcing function in `ContinuousForcing`s and ensuring that
`ContinuousForcing`s are located correctly for velocity fields.
"""
function model_forcing(tracer_names; u=nothing, v=nothing, w=nothing, tracer_forcings...)

    model_field_names = tuple(:u, :v, :w, tracer_names...)

    u = regularize_forcing(u, :u, model_field_names)
    v = regularize_forcing(v, :v, model_field_names)
    w = regularize_forcing(w, :w, model_field_names)

    # Build tuple of user-specified tracer forcings
    specified_tracer_forcings_tuple = Tuple(regularize_forcing(f.second, f.first, model_field_names) for f in tracer_forcings)
    specified_tracer_names = Tuple(f.first for f in tracer_forcings)

    specified_forcings = NamedTuple{specified_tracer_names}(specified_tracer_forcings_tuple)

    # Re-build with defaults for unspecified tracer forcing
    tracer_forcings = with_tracers(tracer_names, specified_forcings, (name, initial_tuple) -> zeroforcing)

    return merge((u=u, v=v, w=w), tracer_forcings)
end


text:
__precompile__()

module DroneHunter

# package code goes here
include("mcts.jl")

export MCTSPolicy

include("batch.jl")

export my_parsim

end # module


text:
const Maybe{T} = Union{T, Nothing} where T

const Symbolic = Union{Expr, Symbol}

struct Bounds
    a :: Number
    b :: Number
end

struct Integral
    expr   :: Union{Symbolic, Number}
    diff   :: Symbolic
    bounds :: Maybe{Bounds}
end

createIntegral(expr, diff::Symbol, bounds = nothing) = Integral(expr, diff, bounds)
function createIntegral(expr, diff::Expr, bounds = nothing)
    # We need to validate that the differential contains ONE AND ONLY ONE variable
    vars = 0
    f(x) = nothing
    f(x::Symbol) = (vars += 1)
    walk(f, diff)

    if vars ≠ 1
        error("the differential must contain strictly one variable.")
    end

    Integral(expr, diff, bounds)
end

struct Thunk
    expr :: Expr
end

createThunk(x) = x

# Eagerly evaluate eval() calls
function createThunk(ex::Expr)
    new_ex = Context.eagerCalc(ex)
    root_eager = ex.head == :call && ex.args[1] ∈ [:eval, :eval!]
    root_eager ? new_ex : Thunk(new_ex)
end

const Value = Union{Symbolic, Number, Integral, Thunk}

struct Assignment
    lvalue :: Union{Symbolic, Integral}
    rvalue :: Value
end

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

# string conversion overload
using Base: string
Base.string(t::Thunk) = Base.string(t.expr)

Base.string(bounds::Bounds) = "_$(wrap(bounds.a))^$(wrap(bounds.b))"

function Base.string(int::Integral)
    bounds = int.bounds ≠ nothing ? Base.string(int.bounds) : ""
    "int$(bounds) $(wrap(int.expr))*d$(wrap(int.diff))"
end


text:
include("../common.jl")

name = "SuiteSparse_GPU"

sources = [
    sources;
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SuiteSparse

# Apply Jameson's shlib patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/SuiteSparse-shlib.patch

# Disable OpenMP as it will probably interfere with blas threads and Julia threads
FLAGS+=(INSTALL="${prefix}" INSTALL_LIB="${libdir}" INSTALL_INCLUDE="${prefix}/include" CFOPENMP=)

if [[ ${target} == *mingw32* ]]; then
    FLAGS+=(UNAME=Windows)
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir} -shared")
else
    FLAGS+=(UNAME="$(uname)")
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir}")
fi

BLAS_NAME=blastrampoline
if [[ ${nbits} == 64 ]]; then
    SUN="-DSUN64 -DLONGBLAS='long long'"
fi

FLAGS+=(BLAS="-l${BLAS_NAME}" LAPACK="-l${BLAS_NAME}")

# CUDA
FLAGS+=(CUDA_PATH="$prefix/cuda")

# METIS
FLAGS+=(MY_METIS_LIB="-lmetis" MY_METIS_INC="${prefix}/include")

make -j${nproc} -C SuiteSparse_config "${FLAGS[@]}" library config

for proj in SuiteSparse_config SuiteSparse_GPURuntime GPUQREngine AMD BTF CAMD CCOLAMD COLAMD CHOLMOD LDL KLU UMFPACK RBio SPQR SLIP_LU Mongoose; do
    make -j${nproc} -C $proj "${FLAGS[@]}" library CFOPENMP="$CFOPENMP"
    make -j${nproc} -C $proj "${FLAGS[@]}" install CFOPENMP="$CFOPENMP"
done

# For now, we'll have to adjust the name of the OpenBLAS library on macOS and FreeBSD.
# Eventually, this should be fixed upstream
if [[ ${target} == *-apple-* ]] || [[ ${target} == *freebsd* ]]; then
    echo "-- Modifying library name for OpenBLAS"

    for nm in libcholmod libspqr libumfpack; do
        # Figure out what version it probably latched on to:
        if [[ ${target} == *-apple-* ]]; then
            OPENBLAS_LINK=$(otool -L ${libdir}/${nm}.dylib | grep lib${BLAS_NAME} | awk '{ print $1 }')
            install_name_tool -change ${OPENBLAS_LINK} @rpath/lib${BLAS_NAME}.dylib ${libdir}/${nm}.dylib
        elif [[ ${target} == *freebsd* ]]; then
            OPENBLAS_LINK=$(readelf -d ${libdir}/${nm}.so | grep lib${BLAS_NAME} | sed -e 's/.*\[\(.*\)\].*/\1/')
            patchelf --replace-needed ${OPENBLAS_LINK} lib${BLAS_NAME}.so ${libdir}/${nm}.so
        fi
    done
fi

# Delete the extra soversion libraries built. https://github.com/JuliaPackaging/Yggdrasil/issues/7
if [[ "${target}" == *-mingw* ]]; then
    rm -f ${libdir}/lib*.*.${dlext}
    rm -f ${libdir}/lib*.*.*.${dlext}
fi

install_license LICENSE.txt

# Compile SuiteSparse_wrapper shim
cd $WORKSPACE/srcdir/SuiteSparse_wrapper
"${CC}" -O2 -shared -fPIC -I${prefix}/include SuiteSparse_wrapper.c -o ${libdir}/libsuitesparse_wrapper.${dlext} -L${libdir} -lcholmod
"""

# Override the default platforms
platforms = [
    Platform("x86_64", "linux"),
]

append!(products, [
    LibraryProduct("libsliplu",                 :libsliplu),
#    LibraryProduct("libmongoose",               :libmongoose),
    LibraryProduct("libGPUQREngine",            :libGPUQREngine),
    LibraryProduct("libSuiteSparse_GPURuntime", :libSuiteSparse_GPURuntime),
])

dependencies = [
    dependencies;
    Dependency("METIS_jll")
    Dependency("MPFR_jll")
    Dependency("GMP_jll")
    BuildDependency(PackageSpec(name="CUDA_full_jll", version=v"10.0.130"))
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"6")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
"""
Types marked with this trait need to implement the following queries

- [`get_free_values(object)`]
- [`get_fe_space(object)`]

"""
FEFunctionStyle(::Type{T}) where T = Val{false}()

FEFunctionStyle(fe_function) = FEFunctionStyle(typeof(fe_function))

"""
"""
function is_a_fe_function(fe_function)
  v = FEFunctionStyle(fe_function)
  get_val_parameter(v)
end

function is_a_fe_function(::Type)
  @unreachable "is_a_fe_function cannot be called on types"
end

"""
"""
function get_free_values(object)
  @abstractmethod
end

"""
"""
function get_cell_values(object)
  @abstractmethod
end

function get_cell_values(object,cellids::AbstractArray)
  reindex(get_cell_values(object),cellids)
end

function get_cell_values(object,cellids::SkeletonPair)
  vals = reindex(get_cell_values(object),cellids)
  f = get_fe_space(object)
  axs = reindex(get_cell_axes(f),cellids)
  merge_cell_dofs_at_skeleton(vals.left,vals.right,axs.left,axs.right)
end

"""
"""
function get_fe_space(object)
  @abstractmethod
end

"""
"""
function test_fe_function(f)
  @test is_a_fe_function(f)
  free_values = get_free_values(f)
  fe_space = get_fe_space(f)
  @test length(free_values) == num_free_dofs(fe_space)
  cell_vals = get_cell_values(f)
end


text:
  ci,fmod =   2
# Can we approximate power-2 limb-darkening with polynomial?

using PyPlot
#include("/Users/ericagol/Computer/Julia/regress.jl")
include("get_coeff.jl")

alpha = big.([0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.99,1.01,1.05,1.1,1.15,1.2,1.25,1.3,1.35,1.4,1.45,1.5,1.55,1.6,1.65,1.7,1.75,1.8])
nmu = 100; mu = logspace(-big(3.0),big(0.0),nmu)
nalpha = length(alpha)
#nmu = 1000; mu = linspace(1e-3,1.0,nmu)
nord = 50
coeff = zeros(BigFloat,nalpha,nord,nord-1)
scat = zeros(BigFloat,nalpha,nord)
for ialpha=1:nalpha
  fnon = mu.^alpha[ialpha]-mu
  for i=2:nord
    fmod = zeros(BigFloat,nmu)
    ci,fmod = get_coeff!(mu,fnon,i,nmu)
    clf()
    plot(mu,convert(Array{Float64,1},fnon))
    plot(mu,convert(Array{Float64,1},fmod))
    plot(mu,convert(Array{Float64,1},fmod-fnon))
    scat[ialpha,i]=convert(Float64,std(fmod-fnon))
    coeff[ialpha,i,1:i-1]=convert(Array{Float64,1},ci)
    println("alpha: ",convert(Float64,alpha[ialpha])," i: ",i," scatter: ",scat[ialpha,i]," max: ",maximum(abs,convert(Array{Float64,1},fmod-fnon)))
  end
end


text:
"""

julia dcgan.jl --outdir ~/dcgan-out
julia dcgan.jl -h # to see all other script options

This example implements a DCGAN (Deep Convolutional Generative Adversarial Network) on MNIST dataset. This implemented model is not identical with the original model. LeNet is used as a base to adapt DCGAN to original MNIST data.1

* Paper url: https://arxiv.org/abs/1511.06434

"""
module DCGAN
using Knet
using Knet: params
using Sloth
using Images
using ArgParse
using JLD2, FileIO
using Printf
using Dates
using Statistics
using Random

include(Knet.dir("data","mnist.jl"))
include(Knet.dir("data","imagenet.jl"))

function main(args)
    o = parse_options(args)
    o[:seed] > 0 && Knet.setseed(o[:seed])

    # load models, data, optimizers
    D = Discriminator(o[:atype])
    G = Generator(o[:zdim], o[:atype])
    xtrn,ytrn,xtst,ytst = mnist()
    dtrn = minibatch(xtrn, ytrn, o[:batchsize]; shuffle=true, xtype=o[:atype])
    init_optimizers!(D, o[:optim])
    init_optimizers!(G, o[:optim])
    z = sample_noise(o[:atype],o[:zdim],prod(o[:gridsize]))

    if o[:outdir] != nothing && !isdir(o[:outdir])
        mkpath(o[:outdir])
        mkpath(joinpath(o[:outdir],"models"))
        mkpath(joinpath(o[:outdir],"generations"))
    end

    # training
    println("training started..."); flush(stdout)
    for epoch = 1:o[:epochs]
        dlossval = glossval = 0
        @time for (x,y) in dtrn
            input = 2 .* x .- 1
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            dlossval += train_discriminator!(D,G,input,noise)
            noise = sample_noise(o[:atype],o[:zdim],length(y))
            glossval += train_generator!(G,D,noise)
        end
        dlossval /= length(dtrn); glossval /= length(dtrn)
        println((:epoch,epoch,:dloss,dlossval,:gloss,glossval))
        flush(stdout)

        # save models and generations
        if o[:outdir] != nothing
            filename = @sprintf("%04d.png",epoch)
            filepath = joinpath(o[:outdir],"generations",filename)
            plot_generations(G; z=z, savefile=filepath,
                             scale=o[:gridscale], gridsize=o[:gridsize])

            filename = @sprintf("%04d.jld2",epoch)
            filepath = joinpath(o[:outdir],"models",filename)
        end
    end

    return D, G
end

function parse_options(args)
    s = ArgParseSettings()
    s.description =
        "Deep Convolutional Generative Adversarial Networks on MNIST."

    @add_arg_table s begin
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}");
         help="array and float type to use")
        ("--batchsize"; arg_type=Int; default=100; help="batch size")
        ("--zdim"; arg_type=Int; default=100; help="noise dimension")
        ("--epochs"; arg_type=Int; default=20; help="# of training epochs")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
        ("--gridsize"; arg_type=Int; nargs=2; default=[8,8])
        ("--gridscale"; arg_type=Float64; default=2.0)
        ("--optim"; default="Adam(;lr=0.0002, beta1=0.5)")
        ("--loadfile"; default=nothing; help="file to load trained models")
        ("--outdir"; default=nothing; help="output dir for models/generations")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:atype] = eval(Meta.parse(o[:atype]))
    if o[:outdir] != nothing
        o[:outdir] = abspath(o[:outdir])
    end
    return o
end

function leaky_relu(x, alpha=0.2)
    pos = max.(0,x)
    neg = min.(0,x) .* alpha
    return pos + neg
end

function sample_noise(atype,zdim,nsamples,mu=0.5,sigma=0.5)
    noise = convert(atype, randn(zdim,nsamples))
    normalized = (noise .- mu) ./ sigma
end

mutable struct DLayer1
    conv
    bnorm
end

function (l::DLayer1)(x; stride=1, padding=0, alpha=0.2, training=true)
    y = l.conv(x; stride=stride, padding=padding)
    y = l.bnorm(y; training=training)
    y = pool(y; mode=2)
end

function DLayer1(ci::Int, co::Int, k::Int; atype=Array{Float64})
    conv = Conv(ci, co, k; atype=atype)
    bnorm = BatchNorm(co; atype=atype)
    return DLayer1(conv, bnorm)
end

mutable struct DLayer2
    linear
    bnorm
end

function (l::DLayer2)(x; alpha=0.2, training=true)
    y = l.linear(x)
    y = l.bnorm(y; training=training)
    y = leaky_relu(y, alpha)
end

function DLayer2(input_dim::Int, output_dim::Int; atype=Array{Float64})
    linear = Linear(input_dim, output_dim; atype=atype, bias=false)
    bnorm = BatchNorm(output_dim; atype=atype)
    return DLayer2(linear, bnorm)
end

mutable struct Discriminator
    dlayer1a
    dlayer1b
    dlayer2a
    softlayer
end

function (model::Discriminator)(x; alpha=0.2, training=true)
    y = x
    y = model.dlayer1a(y; training=training)
    y = model.dlayer1b(y; training=training)
    y = reshape(y, 800, size(y, 4))
    y = model.dlayer2a(y; training=training)
    y = sigm.(model.softlayer(y))
end

function Discriminator(atype)
    dlayer1a = DLayer1(1, 20, 5; atype=atype)
    dlayer1b = DLayer1(20, 50, 5; atype=atype)
    dlayer2a = DLayer2(800, 500; atype=atype)
    softlayer = Linear(500, 1; atype=atype)

    return Discriminator(dlayer1a, dlayer1b, dlayer2a, softlayer)
end

function dloss(D::Discriminator, real_images, fake_images)
    yreal = D(real_images)
    real_loss = -mean(log.(yreal))
    yfake = D(fake_images)
    fake_loss = -mean(log.(1 .- yfake))
    return real_loss+fake_loss
end

mutable struct GLayerA
    linear
    bnorm
end

function (l::GLayerA)(x; training=true)
    y = l.linear(x)
    y = l.bnorm(y; training=true)
    y = relu.(y)
end

function GLayerA(input_dim::Int, output_dim::Int; atype=Array{Float64})
    linear = Linear(input_dim, output_dim; atype=atype, bias=false)
    bnorm = BatchNorm(output_dim; atype=atype)
    return GLayerA(linear, bnorm)
end

mutable struct GLayerB
    deconv
    bnorm
end

function (l::GLayerB)(x; stride=2, training=true)
    y = l.deconv(x; stride=stride)
    y = l.bnorm(y; training=training)
end

function GLayerB(ci::Int, co::Int, k::Int; atype=Array{Float64})
    deconv = Deconv(ci, co, k; atype=atype, bias=false)
    bnorm = BatchNorm(co; atype=atype)
    return GLayerB(deconv, bnorm)
end

mutable struct Generator
    glayer1a
    glayer2a
    glayer3b
    glayer4c
    glayer5b
    deconv
end

function (model::Generator)(z; training=true)
    y = model.glayer1a(z; training=training)
    y = model.glayer2a(y; training=true)
    y = reshape(y, 4, 4, 50, size(y, 2))
    y = model.glayer3b(y; training=training)
    y = relu.(model.glayer4c(y; stride=1, training=training))
    y = model.glayer5b(y; training=training)
    y = tanh.(model.deconv(y))
    return y
end

function Generator(zdim, atype)
    glayer1a = GLayerA(zdim, 500; atype=atype)
    glayer2a = GLayerA(500, 800; atype=atype)
    glayer3b = GLayerB(50, 50, 2; atype=atype)
    glayer4c = GLayerB(50, 20, 5; atype=atype)
    glayer5b = GLayerB(20, 20, 2, atype=atype)
    deconv = Deconv(20, 1, 5; atype=atype)

    return Generator(glayer1a, glayer2a, glayer3b, glayer4c, glayer5b, deconv)
end

function gloss(G::Generator, D::Discriminator, noise)
    fake_images = G(noise)
    ypred = D(fake_images)
    return -mean(log.(ypred))
end

function train_discriminator!(
    D::Discriminator, G::Generator, real_images, noise)
    fake_images = G(noise; training=true)
    nsamples = div(length(real_images),784)
    J = @diff dloss(D, real_images, fake_images)
    for p in params(D)
        g = grad(J, p)
        update!(value(p), g, p.opt)
    end
    return value(J)
end

function train_generator!(G::Generator, D::Discriminator, noise)
    J = @diff gloss(G, D, noise)
    for p in params(G)
        g = grad(J, p)
        update!(value(p), g, p.opt)
    end
    return value(J)
end

function plot_generations(
    G::Generator; z=nothing, gridsize=(8,8), scale=1.0, savefile=nothing)
    if z == nothing
        nimg = prod(gridsize)
        zdim = size(G.deconv.glayer1a.w, 2)
        atype = G.deconv.w isa KnetArray ? KnetArray{Float32} : Array{Float32}
        z = sample_noise(atype, zdim, nimg)
    end
    output = Array(0.5 .* (1 .+ G(z; training=false)))
    images = map(i->output[:,:,:,i], 1:size(output,4))
    grid = make_image_grid(images; gridsize=gridsize, scale=scale)
    if savefile == nothing
        display(colorview(Gray, grid))
    else
        save(savefile, grid)
    end
end

splitdir(PROGRAM_FILE)[end] == "dcgan.jl" && main(ARGS)

end # module

max_stars_repo_name: ilkerkesen/adversarial

max_issues_repo_name: muhammetfaik/GAN

max_forks_repo_name: muhammetfaik/GAN


text:
using Test
using CGP
using Random

CGP.Config.init("cfg/test.yaml")

@testset "Simple fit" begin
    function simple_fit(c::Chromosome)
        process(c, [1.0])[1]
    end
    for ea in CGP.EAs
        dists = CGP.distances[:]
        crosses = CGP.crossovers[:]
        if ea!=cgpneat
            dists = CGP.distances[1:1]
        end
        if ea==oneplus
            crosses = CGP.crossovers[1:1]
        end
        for ct in CGP.chromosomes
            for mut in CGP.mutations
                for cross in crosses
                    for dist in dists
                        CGP.Config.init(Dict("mutate_method" => string("\"", mut, "\""),
                                             "crossover_method" => string("\"", cross, "\""),
                                             "distance_method" => string("\"", dist, "\"")))
                        println(ea, ", ", ct, ", ", mut, ", ", cross, ", ", dist)
                        @testset "$ea $ct $mut $cross $dist" begin
                            max_fit, genes = ea(1, 1, simple_fit; ctype=ct)
                            @test max_fit == 1.0
                            c = ct(genes, 1, 1)
                            @test process(c, [1.0])[1] == 1.0
                        end
                    end
                end
            end
        end
    end
end

@testset "Rosenbrock" begin
    # emulate the rosenbrock function
    function rosenbrock(a::Float64, b::Float64, x::Float64, y::Float64)
        (a-x)^2 + b*(y-x^2)^2
    end
    function rosenbrock_fit(c::Chromosome)
        diff = 0
        for i=1:10
            inps = collect(range(0.0, stop=i, length=4))/10
            diff += (process(c, inps)[1] - rosenbrock(inps...))^2
        end
        -diff
    end
    for ea in CGP.EAs
        for ct in CGP.chromosomes
            println(ea, ", ", ct)
            @testset "$ea $ct" begin
                randc = ct(4, 1)
                randfit = rosenbrock_fit(randc)
                max_fit, genes = ea(4, 1, rosenbrock_fit; ctype=ct)
                @test max_fit > randfit
                c = ct(genes, 4, 1)
                @test rosenbrock_fit(c) == max_fit
            end
        end
    end
end

@testset "Array function" begin
    using TestImages
    using Colors
    function find_the_mandrill(c::Chromosome, inps::Array{Array{Float64,2},1},
                               outs::Array{Float64})
        resout = 0
        for i=1:3
            perm = randperm(length(inps))
            inps = inps[perm]; outs = outs[perm]
            output = process(c, inps)
            resout += sum([(output[i] - outs[i])^2 for i in 1:3])
        end
        -resout
    end
    for ea in CGP.EAs
        for ct in CGP.chromosomes
            println(ea, ", ", ct)
            @testset "$ea $ct" begin
                mandrill = Float64.(Gray.(testimage("mandrill")))
                mandrill /= maximum(mandrill)
                inps = [rand(Float64, size(mandrill)) for i in 1:3]
                outs = zeros(3)
                correct = rand(1:3)
                inps[correct] = mandrill
                outs[correct] = 1.0
                randc = ct(3, 3)
                randfit = find_the_mandrill(randc, inps, outs)
                max_fit, genes = ea(3, 3, c->find_the_mandrill(c, inps, outs);
                                    ctype=ct)
                @test max_fit > randfit
            end
        end
    end
end
nothing


text:
mutable struct AppData

    model::ModelData

    basic_protocol::Dict{String,Any}

    # Samples holds all the sampled data in the sim or app.
    samples::Samples

    configured::Bool

    function AppData()
        o = new()

        o.model = ModelData()

        o.basic_protocol = JSON.parsefile("../data/com_protocol_basic.json")
        o.samples = Samples()
        o.configured = false
        o
    end
end

function config!(data::AppData)
    synapses = Model.synapses(data.model)
    duration = Model.duration(data.model)

    config_samples!(data.samples, synapses, duration)
    data.configured = true
end

function samples(data::AppData)
    data.samples
end

function load_data!(data::AppData)
    load!(data.model, "../data/app.json")
end

function save_data(data::AppData)
    Model.unload(data.model)
end


text:
using MeshSignedDistance
using Test

@testset "MeshSignedDistance.jl" begin
    # Write your tests here.
end

max_stars_repo_name: touste/MeshSignedDistance

max_issues_repo_name: touste/MeshSignedDistance

max_forks_repo_name: touste/MeshSignedDistance