text:
module KernelStat
#-----------------------------------------------------------
#PURPOSE
#	Colin T. Bowers module for kernel functions and statistics
#NOTES
#	Currently this module only really contains a function for estimating bandwidths (needed for estimating block lengths in dependent bootstraps)
#LICENSE
#	MIT License (see github repository for more detail: https://github.com/colintbowers/KernelStat.jl.git)
#-----------------------------------------------------------

#Load any entire modules that are needed (use import ModuleName1, ModuleName2, etc)
using 	StatsBase

#Load any specific variables/functions that are needed (use import ModuleName1.FunctionName1, ModuleName2.FunctionName2, etc)
import 	Base.string,
		Base.show,
		Base.copy

#Specify the variables/functions to export (use export FunctionName1, FunctionName2, etc)
export	evaluate, #Evaluate kernel function at given input
		activedomain, #Function that returns a UnitRange indicating the range over which kernel output is non-zero
		paramdomain, #Function that returns a UnitRange indicating the feasible domain kernel parameters
		setbandwidth!, #Function to adjust bandwidth of input kernel function
		KernelFunction, #Abstract type that nests all kernel function types
		KernelDummy, #Dummy type (useful for update! routines)
		KernelUniform, #Kernel function type
		KernelBartlett, #Kernel function type
		KernelTriangular, #Identical to KernelBartlett
		KernelEpanechnikov, #Kernel function type
		KernelGaussian, #Kernel function type
		KernelParzen2009RK, #Kernel function type
		KernelPR1993FlatTop, #Kernel function type
		KernelP2003FlatTop, #Kernel function type
		KernelPP2002Trap, #Kernel function type
		KernelPP2002Smooth, #Kernel function type
		KernelPR1994SB, #Kernel function for variance estimator from Politis, Romano (1994) "The Stationary Bootstrap"
		hacvariance, #Function for estimating heteroskedasticity and auto-correlation consistent variance
		HACVarianceMethod, #Abstract type for nesting all HAC variance methods
		HACVarianceBasic, #Basic HAC variance estimator
		bandwidth, #Function for estimating bandwidth
		BandwidthMethod, #Abstract type that nests all bandwidth method types
		BandwidthDummy, #Dummy type (useful for update! functions)
		BandwidthMax, #Maximum possible bandwidth
		BandwidthWhiteNoise, #Bandwidth type that uses Gaussian white noise assumption to get correlation confidence bounds
		BandwidthBartlett, #Bandwidth type that uses Bartlett formula to get correlation confidence bounds
		BandwidthP2003 #Bandwidth type implements the bandwidth estimation method of Politis (2003) "Adaptive Bandwidth Choice"

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

#----------------------------------------------------------
#SET CONSTANTS FOR MODULE
#----------------------------------------------------------
const gaussianPremult = (1 / sqrt(2*pi))::Float64

#----------------------------------------------------------
#TYPE
#	abstract super-type and one type for each kernel function
#PURPOSE
#	These types are used to define different kernel functions, parameterised by the fields of the type, and all of which can be called via the evaluate function.
#----------------------------------------------------------
#Abstract super type
abstract KernelFunction
#Dummy type for evaluate function that ditches the indicator function check
type KernelNoIndCheck; end
#General dummy type
type KernelDummy <: KernelFunction ; end
#------- TYPE DEFINITIONS ---------------
#Common kernel functions
type KernelUniform <: KernelFunction
	LB::Float64
	UB::Float64
	function KernelUniform(LB::Float64, UB::Float64)
		LB >= UB && error("lower bound greater than upper bound")
		new(LB, UB)
	end
end
KernelUniform() = KernelUniform(-1.0, 1.0)
KernelUniform(LB::Number, UB::Number) = KernelUniform(convert(Float64, LB), convert(Float64, UB))
type KernelBartlett <: KernelFunction
	bandwidth::Float64
	KernelBartlett(bandwidth::Float64) = bandwidth <= 0 ? error("Bandwidth must be strictly positive") : new(bandwidth)
end
KernelBartlett() = KernelBartlett(1.0)
KernelBartlett(bw::Number) = KernelBartlett(convert(Float64, bw))
type KernelTriangular <: KernelFunction #Identical to KernelBartlett
	bandwidth::Float64
	KernelTriangular(bandwidth::Float64) = bandwidth <= 0 ? error("Bandwidth must be strictly positive") : new(bandwidth)
end
KernelTriangular() = KernelTriangular(1.0)
KernelTriangular(bw::Number) = KernelTriangular(convert(Float64, bw))
type KernelEpanechnikov <: KernelFunction
	bandwidth::Float64
	KernelEpanechnikov(bandwidth::Float64) = bandwidth <= 0 ? error("Bandwidth must be strictly positive") : new(bandwidth)
end
KernelEpanechnikov() = KernelEpanechnikov(1.0)
KernelEpanechnikov(bw::Number) = KernelEpanechnikov(convert(Float64, bw))
type KernelGaussian <: KernelFunction
	bandwidth::Float64
	KernelGaussian(bandwidth::Float64) = (bandwidth <= 0) ? error("Bandwidth must be strictly positive") : new(bandwidth)
end
KernelGaussian() = KernelGaussian(1.0) #Constructor for the most common use-case
KernelGaussian(bw::Number) = KernelGaussian(convert(Float64, bw))
#Kernel functions from specific papers
#Barndorff-Nielsen, Hansen, Lunde, Shephard (2009) "Realized Kernels in Practice"
type KernelParzen2009RK <: KernelFunction
end
#Politis, Romano (1993) "On a family of smoothing kernels of infinite order"
type KernelPR1993FlatTop <: KernelFunction
	m::Float64
	M::Float64
	function KernelPR1993FlatTop(m::Float64, M::Float64)
		(m <= 0 || M <= 0) && error("Parameters of this kernel function must be strictly positive")
		(m >= M) && error("First parameter of this kernel function must be stricty smaller than second parameter")
		new(m, M)
	end
end
KernelPR1993FlatTop() = KernelPR1993FlatTop(0.5, 1.0) #A popular specification of this kernel function
#Politis, Romano (2003) "Adaptive Bandwidth Choice" [NOTE, equivalent to KernelPR1993FlatTop with m = 1/2, and M = 1]
type KernelP2003FlatTop <: KernelFunction; end
#Paparoditis, Politis (2002) "The tapered block bootstrap for general statistics from stationary sequences"
type KernelPP2002Trap <: KernelFunction
	p1::Float64
	KernelPP2002Trap(p1::Float64) = !(0 < p1 <= 0.5) ? error("Invalid kernel function parameter") : new(p1)
end
KernelPP2002Trap() = KernelPP2002Trap(0.43) #Optimal value from Paparoditis, Politis (2002)
#Paparoditis, Politis (2002) "The tapered block bootstrap for general statistics from stationary sequences"
type KernelPP2002Smooth <: KernelFunction
	p1::Float64
	KernelPP2002Smooth(p1::Float64) = !(1 <= p1 < Inf) ? error("Invalid kernel function parameter") : new(p1)
end
KernelPP2002Smooth(p1::Int) = KernelPP2002Smooth(convert(Float64, p1))
KernelPP2002Smooth() = KernelPP2002Smooth(1.3) #Optimal value from Paparoditis, Politis (2002)
#Politis, Romano (1994) "The Stationary Bootstrap" (kernel function for variance estimator, see eqn 7)
type KernelPR1994SB <: KernelFunction
	bandwidth::Float64 #The upper bound (usually the number of observations)
	p::Float64 #parameter of geometric distribution (1 / expected block length)
	function KernelPR1994SB(bandwidth::Float64, p::Float64)
		bandwidth <= 0 && error("Bandwidth must be strictly positive")
		!(0 < p < 1) && error("Geometric distribution parameter must lie on the interval (0, 1)")
		new(bandwidth, p)
	end
end
KernelPR1994SB() = KernelPR1994SB(1.0, 0.1)
KernelPR1994SB(bandwidth::Number) = KernelPR1994SB(convert(Float64, bandwidth), 0.1)
KernelPR1994SB(bandwidth::Number, p::Number) = KernelPR1994SB(convert(Float64, bandwidth), convert(Float64, p))
#------- METHODS ------
#string method gets string representation of type
string(kT::KernelDummy) = "dummy"
string(kT::KernelUniform) = "uniform"
string(kT::KernelBartlett) = "bartlett"
string(kT::KernelTriangular) = "triangular"
string(kT::KernelEpanechnikov) = "epanechnikov"
string(kT::KernelGaussian) = "gaussian"
string(kT::KernelParzen2009RK) = "parzen2009RK"
string(kT::KernelPR1993FlatTop) = "PR1993FlatTop"
string(kT::KernelP2003FlatTop) = "P2003FlatTop"
string(kT::KernelPP2002Trap) = "PP2002Trap"
string(kT::KernelPP2002Smooth) = "PP2002Smooth"
string(kT::KernelPR1994SB) = "PR1994SB"
#show method
function show(io::IO, kF::KernelFunction)
	println(io, "Kernel function = " * string(kF))
	fieldNames = names(kF)
	for n = 1:length(fieldNames)
		println("    field " * string(fieldNames[n]) * " = " * string(getfield(kF, n)))
	end
end
show{T<:KernelFunction}(k::T) = show(STDOUT, k)
#deepcopy method
# function deepcopy(kF::KernelFunction)
# 	tempArgs = [ deepcopy(getfield(kF, i)) for i = 1:length(names(kF)) ]
# 	return(eval(parse(string(typeof(kF)) * "(tempArgs...)")))
# end
#------ EVALUATE method for evaluating kernel function at a given value ------------------------------
#Common Kernel functions
evaluate(x::Number, kT::KernelUniform) = indicator(x, UnitRange(kT.LB, kT.UB)) * (1 / (kT.UB - kT.LB))
evaluate(x::Number, kT::KernelUniform, ::KernelNoIndCheck) = 1 / (kT.UB - kT.LB)
evaluate{T<:Union(KernelBartlett, KernelTriangular)}(x::Number, kT::T) = indicator(x / kT.bandwidth, UnitRange(-1, 1)) * (1 - abs(x / kT.bandwidth))
evaluate{T<:Union(KernelBartlett, KernelTriangular)}(x::Number, kT::T, ::KernelNoIndCheck) = 1 - abs(x / kT.bandwidth)
evaluate(x::Number, kT::KernelEpanechnikov) = indicator(x / kT.bandwidth, UnitRange(-1, 1)) * (0.75 / kT.bandwidth) * (1 - (x^2 / kT.bandwidth^2))
evaluate(x::Number, kT::KernelEpanechnikov, ::KernelNoIndCheck) = (0.75 / kT.bandwidth) * (1 - (x^2 / kT.bandwidth^2))
evaluate(x::Number, kT::KernelGaussian) = (1 / kT.bandwidth) * gaussianPremult * exp(-1 * (x^2 / (2 * kT.bandwidth^2)))
evaluate(x::Number, kT::KernelGaussian, ::KernelNoIndCheck) = evaluate(x, kT)
function evaluate{T<:Number}(x::T, kT::KernelParzen2009RK)
	if x < 0; return(zero(T))
	elseif x <= 0.5; return(1 - 6*x^2 + 6*x^3)
	elseif x <= 1; return(2*(1 - x)^3)
	else return(zero(T))
	end
end
function evaluate{T<:Number}(x::T, kT::KernelParzen2009RK, ::KernelNoIndCheck)
	if x <= 0.5; return(1 - 6*x^2 + 6*x^3)
	else; return(2*(1 - x)^3)
	end
end
#Kernel functions from specific papers (note, may not satisfy standard kernel function properties)
#KernelPR1993FlatTop
function evaluate{T<:Number}(x::T, kT::KernelPR1993FlatTop)
	if abs(x) <= kT.m; return(one(T))
	elseif abs(x) <= kT.M; return(1 - (abs(x) - kT.m) / (kT.M - kT.m))
	else; return(zero(T))
	end
end
evaluate(x::Number, kT::KernelPR1993FlatTop, ::KernelNoIndCheck) = evaluate(x, kT)
#KernelP2003FlatTop (equivalent to KernelPR1993FlatTop with m=0.5 and M=1)
function evaluate{T<:Number}(x::T, kT::KernelP2003FlatTop)
	if abs(x) <= 0.5; return(one(T))
	elseif abs(x) <= 1; return(2 * (1 - abs(x)))
	else; return(zero(T))
	end
end
evaluate(x::Number, kT::KernelP2003FlatTop, ::KernelNoIndCheck) = evaluate(x, kT)
#KernelPP2002Trap
function evaluate{T<:Number}(x::T, kT::KernelPP2002Trap)
	if x < 0; return(zero(T))
	elseif x < kT.p1; return(x / kT.p1)
	elseif x < 1 - kT.p1; return(one(T))
	elseif x < 1; return((1 - x) / kT.p1)
	else; return(zero(T))
	end
end
evaluate(x::Number, kT::KernelPP2002Trap, ::KernelNoIndCheck) = evaluate(x, kT)
#KernelPP2002Smooth
function evaluate{T<:Number}(x::T, kT::KernelPP2002Smooth)
	if 0 <= x <= 1; return(1 - abs(2*x - 1)^kT.p1)
	else; return(zero(T))
	end
end
evaluate(x::Number, kT::KernelPP2002Smooth, ::KernelNoIndCheck) = 1 - abs(2*x - 1)^kT.p1
#KernelPR1994SB
evaluate(x::Number, kT::KernelPR1994SB) = indicator(x, UnitRange(0.0, kT.bandwidth)) * ((1 - x/kT.bandwidth)*(1 - kT.p)^x + (x/kT.bandwidth)*(1 - kT.p)^(kT.bandwidth - x))
evaluate(x::Number, kT::KernelPR1994SB, ::KernelNoIndCheck) = (1 - x/kT.bandwidth)*(1 - kT.p)^x + (x/kT.bandwidth)*(1 - kT.p)^(kT.bandwidth - x)
#Array input wrappers
evaluate{T<:Number}(x::AbstractVector{T}, kT::KernelFunction) = [ evaluate(x[n], kT) for n = 1:length(x) ]
evaluate{T<:Number}(x::AbstractMatrix{T}, kT::KernelFunction) = [ evaluate(x[n, m], kT) for n = 1:size(x, 1), m = 1:size(x, 2) ]
#-----kernelweight ------------------------
#This function transforms a kernelFunction into a set of weights
function kernelweight(kF::KernelFunction, numWeight::Int)
	r = activedomain(kF)
	r.start > 0 && error("Kernel function active domain start point too high to construct weights")
	r.stop < numWeight - 1 && error("Kernel function active domain end point not high enough to construct weights")
	w = [ evaluate(n, kF) for n = 0:numWeight-1 ]
	w *= (1 / sum(w)) #scale weights so they sum to unity
	return(w)
end
#-----activedomain and paramdomain ------------------------
#The purpose of activedomain is to return the domain (of kernel function input) over which the output of the kernel function is active
#The purpose of paramdomain is to return the domain over which the indicated parameter can be defined.
activedomain(kT::KernelUniform) = UnitRange(kT.LB, kT.UB)
activedomain(kT::KernelEpanechnikov) = UnitRange(-kT.bandwidth, kT.bandwidth)
activedomain{T<:Union(KernelBartlett, KernelTriangular)}(kT::T) = UnitRange(-kT.bandwidth, kT.bandwidth)
activedomain(kT::KernelGaussian) = UnitRange(-Inf, Inf)
activedomain(kT::KernelParzen2009RK) = UnitRange(0, 1)
activedomain(kT::KernelPR1993FlatTop) = UnitRange(-kT.M, kT.M)
activedomain(kT::KernelP2003FlatTop) = UnitRange(-1, 1)
activedomain(kT::KernelPP2002Trap) = UnitRange(0, 1)
activedomain(kT::KernelPP2002Smooth) = UnitRange(0, 1)
activedomain(kT::KernelPR1994SB) = UnitRange(0, kT.bandwidth)
function paramdomain(kT::KernelUniform, pNum::Int)
	pNum == 1 && return(UnitRange(-Inf, kT.UB))
	pNum == 2 && return(UnitRange(kT.LB, Inf))
	error("Invalid parameter number")
end
paramdomain(kT::KernelEpanechnikov, pNum::Int) = pNum == 1 ? UnitRange(nextfloat(0.0), prevfloat(Inf)) : error("Invalid parameter number")
paramdomain{T<:Union(KernelBartlett, KernelTriangular)}(kT::T, pNum::Int) = pNum == 1 ? UnitRange(nextfloat(0.0), prevfloat(Inf)) : error("Invalid parameter number")
paramdomain(kT::KernelGaussian, paramNum::Int) = pNum == 1 ? UnitRange(nextfloat(0.0), prevfloat(Inf)) : error("Invalid parameter number")
paramdomain(kT::KernelParzen2009RK, pNum::Int) = error("Kernel function has no parameters")
function paramdomain(kT::KernelPR1993FlatTop, paramNum::Int)
	paramNum == 1 && return(UnitRange(nextfloat(0.0), Inf))
	paramNum == 2 && return(UnitRange(nextfloat(0.0), Inf))
	error("Invalid parameter number")
end
paramdomain(kT::KernelP2003FlatTop, pNum::Int) = error("Kernel function has no parameters")
paramdomain(kT::KernelPP2002Trap, pNum::Int) = (pNum == 1) ? UnitRange(nextfloat(0.0), 0.5) : error("Invalid parameter number")
paramdomain(kT::KernelPP2002Smooth) = (pNum == 1) ? UnitRange(1, Inf) : error("Invalid parameter number")
function paramdomain(kT::KernelPR1994SB, pNum::Int)
	pNum == 1 && return(UnitRange(0, Inf))
	pNum == 2 && return(UnitRange(nextfloat(0.0), prevfloat(1.0)))
	error("Invalid parameter number")
end
function setbandwidth!(kT::KernelUniform, LB::Number, UB::Number)
	kT.LB = convert(Float64, LB)
	kT.UB = convert(Float64, UB)
end
function setbandwidth!(kT::KernelUniform, bandwidth::Number)
	kT.LB = 0.0
	kT.UB = convert(Float64, bandwidth)
end
function setbandwidth!{T<:Union(KernelEpanechnikov, KernelGaussian, KernelPR1994SB, KernelTriangular, KernelBartlett)}(kT::T, bandwidth::Number)
	kT.bandwidth = convert(Float64, bandwidth)
end
setbandwidth!(kT::KernelFunction, bandwidth::Number) = error("The input kernel function does not have a natural bandwidth measure")

#---------- TYPES FOR DIFFERENT METHODS OF BANDWIDTH ESTIMATION
#Abstract super-type
abstract BandwidthMethod
#General dummy type (useful for update! functions)
type BandwidthDummy <: BandwidthMethod; end
#Dummy type for using maximum possible bandwidth
type BandwidthMax <: BandwidthMethod; end
#type definitions
type BandwidthStockWatsonDefault
	adjustmentTerm::Float64
	BandwidthStockWatsonDefault(adjustmentTerm::Float64) = (adjustmentTerm <= 0) ? error("Adjustment scalar must be strictly greater than zero") : new(adjustmentTerm)
end
BandwidthStockWatsonDefault() = BandwidthStockWatsonDefault(1.0)
type BandwidthWhiteNoise <: BandwidthMethod
	adjustmentTerm::Float64
	BandwidthWhiteNoise(adjustmentTerm::Float64) = (adjustmentTerm <= 0) ? error("Adjustment scalar must be strictly greater than zero") : new(adjustmentTerm)
end
BandwidthWhiteNoise() = BandwidthWhiteNoise(1.0)
type BandwidthBartlett <:BandwidthMethod
	adjustmentTerm::Float64
	BandwidthBartlett(adjustmentTerm::Float64) = (adjustmentTerm <= 0) ? error("Adjustment scalar must be strictly greater than zero") : new(adjustmentTerm)
end
BandwidthBartlett() = BandwidthBartlett(1.0)
type BandwidthP2003 <:BandwidthMethod
	adjustmentTerm::Float64
	c::Float64
	K::Int
	function BandwidthP2003(adjustmentTerm::Float64, c::Float64, K::Int)
		adjustmentTerm <= 0 && error("Adjustment scalar must be strictly greater than zero")
		c <= 0 && error("c must be strictly positive")
		K < 1 && error("K must be strictly positive")
		new(adjustmentTerm, c, K)
	end
end
BandwidthP2003() = BandwidthP2003(1.0, 2.0, 5)
BandwidthP2003(numObs::Int) = BandwidthP2003(1.0, 2.0, max(5, convert(Int, ceil(sqrt(log(10, numObs))))))
#-------------- METHODS---------
#string
string(x::BandwidthMax) = "maximumBandwidth"
string(x::BandwidthStockWatsonDefault) = "stockWatsonDefault"
string(x::BandwidthWhiteNoise) = "whiteNoise"
string(x::BandwidthBartlett) = "bartlett"
string(x::BandwidthP2003) = "P2003"
#show method for BandwidthMethod
function show(io::IO, b::BandwidthMethod)
	println(io, "bandwidth method = " * string(b))
	fieldNames = names(b)
	for n = 1:length(fieldNames)
		println("    field " * string(fieldNames[n]) * " = " * string(getfield(b, n)))
	end
end
show{T<:BandwidthMethod}(b::T) = show(STDOUT, b)
#deepcopy
function deepcopy(x::BandwidthMethod)
	tempArgs = [ deepcopy(getfield(x, i)) for i = 1:length(names(x)) ]
	return(eval(parse(string(typeof(x)) * "(tempArgs...)")))
end
#---------- TYPES FOR DIFFERENT METHODS OF HAC-VARIANCE ESTIMATION
#Abstract super-type
abstract HACVarianceMethod
#type definitions
type HACVarianceBasic <: HACVarianceMethod
	kernelFunction::KernelFunction
	bandwidthMethod::BandwidthMethod
	function HACVarianceBasic(kernelFunction::KernelFunction, bandwidthMethod::BandwidthMethod)
		!(typeof(kernelFunction) <: Union(KernelUniform, KernelGaussian, KernelEpanechnikov, KernelPR1994SB, KernelBartlett, KernelTriangular)) && error("This kernel function does not have appropriate properties for HAC variance estimation")
		new(kernelFunction, bandwidthMethod)
	end
end
HACVarianceBasic() = HACVarianceBasic(KernelEpanechnikov(), BandwidthP2003())
HACVarianceBasic(numObs::Int) = HACVarianceBasic(KernelEpanechnikov(), BandwidthP2003(numObs))
HACVarianceBasic{T<:Number}(x::Vector{T}) = HACVarianceBasic(KernelEpanechnikov(), BandwidthP2003(length(x)))
#---------- METHODS ----------------
string(::HACVarianceBasic) = "hacVarianceBasic"
function show(io::IO, x::HACVarianceBasic)
	println(io, "HAC variance estimator type. Parameters are:")
	show(io, x.kernelFunction)
	show(io, x.bandwidthMethod)
end
#deepcopy
function deepcopy(x::HACVarianceMethod)
	tempArgs = [ deepcopy(getfield(x, i)) for i = 1:length(names(x)) ]
	return(eval(parse(string(typeof(x)) * "(tempArgs...)")))
end

#----------------------------------------------------------
#FUNCTION
#	bandwidth
#INPUT
#	(x::Vector{T<:Number}, method::BandwidthMethod): Estimate bandwidth using specified method on input data x
#OUTPUT
#	Output is the tuple (M::Int, xVar::Float64, covVec::Vector{Float64})
#		M is the bandwidth estimator. Minimum possible estimate is 1, maximum possible estimate is length(x) - 1.
#		xVar is the variance of the input data x but scaled using (1/length(x)) instead of (1/(length(x)-1))
#		covVec is the vector of sample autocovariances (NOT autocorrelations) from lag 1 to lag K of input x. K is determined within the function.
#	We pass all three of these variables as output because xVar and covVec are frequently used in the calling function, and so this saves duplicating computations.
#NOTES
#	M can be close to length(x). Typically, this is not a good situation, so the user will probably want to control for this in the calling function and bind maximum M to a more reasonable value.
#	M has an enforced minimum value of 2. Allowing results of M < 2 typically result in kernel function evaluations that always equal 0.
#	Autocorrelations are calculated using autocor (e.g. fft methods) in blocks of 20 (i.e. if bandwidth estimate is not found using first 20 autocorrelations, the next 20 will be calculated, and so on)
#----------------------------------------------------------
function bandwidth{T<:Number}(x::AbstractVector{T}, method::BandwidthMax)
	length(x) < 2 && error("Input data must have at least two observations")
	return(length(x)-1, ((length(x)-1) / length(x)) * var(x), Array(Float64, 0))
end
function bandwidth{T<:Number}(x::AbstractVector{T}, method::BandwidthStockWatsonDefault)
	length(x) < 2 && error("Input data must have at least two observations")
	return(0.75 * length(x)^(1/3), ((length(x)-1) / length(x)) * var(x), Array(Float64, 0))
end
function bandwidth{T<:Number}(x::AbstractVector{T}, method::BandwidthWhiteNoise)
	length(x) < 2 && error("Input data must have at least two observations")
	mHat = 1
	corVec = Array(Float64, 0)
	append!(corVec, autocor(x, 1:min(20, length(x)-1)))
	whiteNoiseBound = sqrt(1 / length(x))
	nullRejected = false
	if abs(corVec[1]) > 1.96 * whiteNoiseBound #"1.96*" implies we are hard-coding hypothesis test at ~95% confidence level. Note, if this condition is not triggered, then we detect no autocorrelation and so leave mHat at 1
		mSt = 1
		while mSt < length(x)
			for m = mSt:length(corVec)-1
				if abs(corVec[m+1]) < 1.96 * whiteNoiseBound
					nullRejected = true
					mHat = m + 1
					break
				end
			end
			if nullRejected == true
				break
			else
				mSt = length(corVec)
				append!(corVec, autocor(x, mSt+1:min(mSt+20, length(x)))) #Grow auto-correlation vector by 20 and try again
			end
		end
		if nullRejected == false
			mHat = length(x)
		end
	end
	M = convert(Int, ceil(method.adjustmentTerm * 2 * mHat)) #"2*" is a standard rule, see e.g. Politis (2003) "Adaptive Bandwidth Choice". method.adjustmentTerm allows user to tinker with output.
	if M > length(x) - 1
		M = length(x) - 1 #Set maximum value for M
	elseif M < 2
		M = 2 #Set minimum value for M
	end
	xVar = ((length(x)-1) /  length(x)) * var(x) #Used to scale autocorrelations to autocovariances for return argument
	return(M, xVar, xVar * corVec)
end
function bandwidth{T<:Number}(x::AbstractVector{T}, method::BandwidthBartlett)
	length(x) < 2 && error("Input data must have at least two observations")
	mHat = 1
	corVec = Array(Float64, 0)
	append!(corVec, autocor(x, 1:min(20, length(x)-1)))
	bartlettBound = sqrt(1 / length(x))
	nullRejected = false
	if abs(corVec[1]) > 1.96 * bartlettBound #"1.96*" implies we are hard-coding hypothesis test at ~95% confidence level. Note, if this condition is not triggered, then we detect no autocorrelation and so leave mHat at 1
		corSum = 0.0 #Sum of squared autocorrelations. We add lags to the overall sum iteratively within the loop
		mSt = 1
		while mSt < length(x)
			for m = mSt:length(corVec)-1
				corSum += 2*corVec[m]^2
				bartlettBound = sqrt((1 + corSum) / length(x))
				if abs(corVec[m+1]) < 1.96 * bartlettBound
					nullRejected = true
					mHat = m + 1
					break
				end
			end
			if nullRejected == true
				break
			else
				mSt = length(corVec)
				append!(corVec, autocor(x, mSt+1:min(mSt+20, length(x)))) #Grow auto-correlation vector by 20 and try again
			end
		end
		if nullRejected == false
			mHat = length(x)
		end
	end
	M = convert(Int, ceil(method.adjustmentTerm * 2 * mHat)) #"2*" is a standard rule, see e.g. Politis (2003) "Adaptive Bandwidth Choice". method.adjustmentTerm allows user to tinker with output.
	if M > length(x) - 1
		M = length(x) - 1 #Set maximum value for M
	elseif M < 2
		M = 2 #Set minimum value for M
	end
	xVar = ((length(x)-1) /  length(x)) * var(x) #Used to scale autocorrelations to autocovariances for return argument
	return(M, xVar, xVar * corVec)
end
function bandwidth{T<:Number}(x::AbstractVector{T}, method::BandwidthP2003)
	length(x) < 2 && error("Input data must have at least two observations")
	mHat = 1
	corVec = Array(Float64, 0)
	append!(corVec, autocor(x, 1:min(20, length(x)-1))) #Add first block of autocorrelations
	corBound = method.c * sqrt(log(10, length(x)) / length(x)) #Note, use of log base 10 is deliberate and recommended in Politis (2003)
	KCounter = 0
	mHatFound = false
	for k = 1:length(x)-2
		if abs(corVec[k]) < corBound
			KCounter += 1 #Bound is satisfied so add one to counter
		else
			KCounter = 0 #Bound is not satisfied, so reset the counter
		end
		if KCounter >= method.K #We found K autocorrelations in a row that satisfy the bound, so break.
			mHat = k - method.K + 1
			mHatFound = true
			break
		end
		if k == length(corVec) #We've run out of autocorrelations to check, so we need to add another block of them to corVec
			append!(corVec, autocor(x, length(corVec)+1:min(length(corVec)+20, length(x)-1)))
		end
	end
	if mHatFound == false
		mHat = length(x) - 1
	end
	M = convert(Int, ceil(method.adjustmentTerm * 2 * mHat)) #"2*" is a standard rule, see e.g. Politis (2003) "Adaptive Bandwidth Choice". method.adjustmentTerm allows user to tinker with output.
	if M > length(x) - 1
		M = length(x) - 1 #Set maximum value for M
	elseif M < 2
		M = 2 #Set minimum value for M
	end
	xVar = ((length(x)-1) /  length(x)) * var(x) #Used to scale autocorrelations to autocovariances for return argument
	return(M, xVar, xVar * corVec)
end

#----------------------------------------------------------
#FUNCTION
#	hacvariance
#INPUT
#	(x::Vector{T<:Number}, method::HACVarianceMethod): Estimate variance using specified method
#OUTPUT
#	Output is a Float64 variance estimator
#NOTES
#----------------------------------------------------------
function hacvariance{T<:Number}(x::AbstractVector{T}, method::HACVarianceBasic)
	(M, xVar, xCov) = bandwidth(x, method.bandwidthMethod)
	length(xCov) < M && append!(xCov, autocov(x, length(xCov)+1:M)) #Get any additional autocovariances that we might need
	setbandwidth!(method.kernelFunction, M)
	kernelAdjTerm = 1 / evaluate(0.0, method.kernelFunction) #Used to scale kernel functions that don't satisfy k(0) = 1
	v = xVar
	noIndCheck = KernelNoIndCheck() #Since we've set the bandwidth to M, we can safely ignore indicator checks
	for m = 1:M
		v += 2 * kernelAdjTerm * evaluate(m, method.kernelFunction, noIndCheck) * xCov[m]
	end
	v = convert(Float64, max(v, 0.0)) #variance estimators should be non-negative
	return(v, M)
end

#Non-exported indicator function for if a value lies in a specified range
indicator{T<:Number}(x::T, LB::T, UB::T) = (LB <= x <= UB) ? 1 : 0
indicator{T<:Number}(x::T, r::UnitRange) = (r.start <= x <= r.stop) ? 1 : 0

end # module


text:
using Notes
using Test

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

ext: jl


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    SOCtoRSOCBridge{T} <: Bridges.Variable.AbstractBridge

`SOCtoRSOCBridge` implements the following reformulation:

 * ``||x||_2 \\le t`` into ``2uv \\ge ||w||_2^2``, with the substitution rules
   ``t = \\frac{u}{\\sqrt 2} + \\frac{v}{\\sqrt 2}``,
   ``x = (\\frac{u}{\\sqrt 2} - \\frac{v}{\\sqrt 2}, w)``.

## Assumptions

 * `SOCtoRSOCBridge` assumes that ``|x| \\ge 1``.

## Source node

`SOCtoRSOCBridge` supports:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.SecondOrderCone`](@ref)

## Target node

`SOCtoRSOCBridge` creates:

 * [`MOI.VectorOfVariables`](@ref) in [`MOI.RotatedSecondOrderCone`](@ref)
"""
struct SOCtoRSOCBridge{T} <:
       SetMapBridge{T,MOI.RotatedSecondOrderCone,MOI.SecondOrderCone}
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{
        MOI.VectorOfVariables,
        MOI.RotatedSecondOrderCone,
    }
end

const SOCtoRSOC{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoRSOCBridge{T},OT}

function MOI.Bridges.map_set(
    ::Type{<:SOCtoRSOCBridge},
    set::MOI.RotatedSecondOrderCone,
)
    return MOI.SecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.inverse_map_set(
    ::Type{<:SOCtoRSOCBridge},
    set::MOI.SecondOrderCone,
)
    if MOI.dimension(set) == 1
        error(
            "Unable to reformulate a `SecondOrderCone` into a " *
            "`RotatedSecondOrderCone` because the dimension of `1` is too " *
            "small",
        )
    end
    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
end

function MOI.Bridges.map_function(
    ::Type{<:SOCtoRSOCBridge{T}},
    func,
    i::MOI.Bridges.IndexInVector,
) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    if i.value > 2
        return scalars[i.value]
    end
    t, u = scalars[1], scalars[2]
    ts = MOI.Utilities.operate!(/, T, t, sqrt(T(2)))
    us = MOI.Utilities.operate!(/, T, u, sqrt(T(2)))
    if i.value == 1
        return MOI.Utilities.operate!(+, T, ts, us)
    else
        return MOI.Utilities.operate!(-, T, ts, us)
    end
end

function MOI.Bridges.map_function(::Type{<:SOCtoRSOCBridge{T}}, func) where {T}
    scalars = MOI.Utilities.eachscalar(func)
    t, u, x = scalars[1], scalars[2], scalars[3:end]
    ts = MOI.Utilities.operate!(/, T, t, sqrt(T(2)))
    us = MOI.Utilities.operate!(/, T, u, sqrt(T(2)))
    return MOI.Utilities.operate(vcat, T, ts + us, ts - us, x)
end

# The map is an involution
function MOI.Bridges.inverse_map_function(BT::Type{<:SOCtoRSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is symmetric
function MOI.Bridges.adjoint_map_function(BT::Type{<:SOCtoRSOCBridge}, func)
    return MOI.Bridges.map_function(BT, func)
end

# The map is a symmetric involution
function MOI.Bridges.inverse_adjoint_map_function(
    BT::Type{<:SOCtoRSOCBridge},
    func,
)
    return MOI.Bridges.map_function(BT, func)
end


text:
# ---
# title: 105. Construct Binary Tree from Preorder and Inorder Traversal
# id: problem105
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Tree, Depth-first Search
# link: <https://leetcode.com/problems/construct-binary-tree-from-preorder-and-inorder-traversal/description/>
# hidden: true
# ---
#
# Given preorder and inorder traversal of a tree, construct the binary tree.
#
# **Note:**
# You may assume that duplicates do not exist in the tree.
#
# For example, given
#
#
#
#     preorder = [3,9,20,15,7]
#     inorder = [9,3,15,20,7]
#
# Return the following binary tree:
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#
## @lc code=start
using LeetCode

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


text:
using Test, AdvancedHMC, LinearAlgebra, UnicodePlots, Random
using Statistics: mean, var, cov
include("common.jl")

@testset "Matrix mode" begin
    n_chains_max = 20
    n_chains_list = collect(1:n_chains_max)
    θ_init_list = [rand(D, n_chains) for n_chains in n_chains_list]
    ϵ = 0.2
    lf = Leapfrog(ϵ)
    i_test = 5
    lfi = Leapfrog(fill(ϵ, i_test))
    lfi_jittered = JitteredLeapfrog(fill(ϵ, i_test), 1.0)
    n_steps = 10
    n_samples = 12_000
    n_adapts = 2_000

    for metricT in [
        UnitEuclideanMetric,
        DiagEuclideanMetric,
        # DenseEuclideanMetric  # not supported at the moment
    ], τ in [
        StaticTrajectory(lfi, n_steps),
        StaticTrajectory(lfi_jittered, n_steps),
        HMCDA(lf, ϵ * n_steps)
    ]
        n_chains = n_chains_list[i_test]
        metric = metricT((D, n_chains))
        h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
        # NoAdaptation
        samples, stats = sample(h, τ, θ_init_list[i_test], n_samples; verbose=false)
        @test mean(samples) ≈ zeros(D, n_chains) atol=RNDATOL * n_chains
        # Adaptation
        for adaptor in [
            Preconditioner(metric),
            NesterovDualAveraging(0.8, lfi.ϵ),
            NaiveHMCAdaptor(
                Preconditioner(metric),
                NesterovDualAveraging(0.8, lfi.ϵ),
            ),
            StanHMCAdaptor(
                Preconditioner(metric),
                NesterovDualAveraging(0.8, lfi.ϵ),
            ),
        ]
            τ isa HMCDA && continue
            samples, stats = sample(h, τ, θ_init_list[i_test], n_samples, adaptor, n_adapts; verbose=false, progress=false)
            @test mean(samples) ≈ zeros(D, n_chains) atol=RNDATOL * n_chains
        end
        # Passing a vector of same RNGs
        rng = [MersenneTwister(1) for _ in 1:n_chains]
        h = Hamiltonian(metricT((D, n_chains)), ℓπ, ∂ℓπ∂θ)
        θ_init = repeat(rand(D), 1, n_chains)
        samples, stats = sample(rng, h, τ, θ_init, n_samples; verbose=false)
        all_same = true
        for i_sample in 2:10
            for j in 2:n_chains
                all_same = all_same && samples[i_sample][:,j] == samples[i_sample][:,1]
            end
        end
        @test all_same
    end
    @info "Adaptation tests for HMCDA with NesterovDualAveraging are skipped"

    # Simple time benchmark
    let metricT=UnitEuclideanMetric
        τ = StaticTrajectory(lf, n_steps)

        time_mat = Vector{Float64}(undef, n_chains_max)
        for (i, n_chains) in enumerate(n_chains_list)
            h = Hamiltonian(metricT((D, n_chains)), ℓπ, ∂ℓπ∂θ)
            t = @elapsed samples, stats = sample(h, τ, θ_init_list[i], n_samples; verbose=false)
            time_mat[i] = t
        end

        # Time for multiple runs of single chain
        time_seperate = Vector{Float64}(undef, n_chains_max)

        for (i, n_chains) in enumerate(n_chains_list)
            t = @elapsed for j in 1:n_chains
                h = Hamiltonian(metricT(D), ℓπ, ∂ℓπ∂θ)
                samples, stats = sample(h, τ, θ_init_list[i][:,j], n_samples; verbose=false)
            end
            time_seperate[i] = t
        end

        # Make plot
        fig = lineplot(
            collect(1:n_chains_max),
            time_mat;
            title="Scalabiliry of multiple chains",
            name="matrix parallism",
            xlabel="Num of chains",
            ylabel="Time (s)"
        )
        lineplot!(fig, collect(n_chains_list), time_seperate; color=:blue, name="seperate")
        println(); show(fig); println(); println()
    end
end


text:
using Random
using Images, CoordinateTransformations

#taken from
#https://github.com/FluxML/Metalhead.jl/blob/b3c09e58ffa907ab1ef8946a6f9e45858bea0a93/src/utils.jl#L28
function preprocess(img)
    img = channelview(imresize(img, (224, 224))) .* float64(1)
    img = permutedims(img, (3,2,1))
    img
end

function copyimg(img)
    # to reset index of imageview
    imresize(img,size(img))
end

struct Dataset
    len::Int
    data::Array{Tuple{String,Int64},1}
    augment::Bool
    image_cache::Dict{Int,Array{RGB{Normed{UInt8,8}},2}}
    use_cache::Bool
    function Dataset(data; train=true)
        augment=train
        use_cache=train
        image_cache = Dict{Int,Array{RGB{Normed{UInt8,8}},2}}()
        new(length(data), data, augment, image_cache, use_cache)
    end
end

function get_example(dataset::Dataset, i::Int)
    path, label = dataset.data[i]
    if dataset.use_cache && haskey(dataset.image_cache, i)
        img = dataset.image_cache[i]
    else
        img = load(path)
        dataset.image_cache[i] = img
    end
    img = copyimg(img)
    if dataset.augment
        #scale transform
        img = imresize(img, (rand(368:512), rand(368:512)))
        #crop
        h, w = size(img)
        b_h = rand(1:h - 224)
        b_w = rand(1:w - 224)
        img = img[b_h:b_h+224, b_w:b_w+224]
        tfm = LinearMap(RotMatrix(rand(range(-pi/12,length=50,stop=pi/12))))
        img = warp(img, tfm)
        img = imresize(img,(224, 224))

        domirror = rand([true, false])
        if domirror
            img = channelview(img)
            newimg = zeros(eltype(img), size(img))
            x_range = size(img)[3]
            for i in 1:x_range
                newimg[:,:,x_range - i + 1] = img[:, :, i]
            end
            img = colorview(RGB, newimg)
        end
    end
    img = preprocess(img)
    return img, label
end

Base.length(dataset::Dataset) = dataset.len


text:
const Δt_sec = 86400.0
const ldd_mv = 255

# read the staticmaps into memory
nc = NCDataset(staticmaps_rhine_path)
ldd_2d = Wflow.ncread(nc, "ldd")

inds, _ = Wflow.active_indices(ldd_2d, ldd_mv)
n = length(inds)

# take out only the active cells
ldd = ldd_2d[inds]
slope = Wflow.ncread(nc, "slope"; sel = inds)
N = Wflow.ncread(nc, "N"; sel = inds)
Qold = Wflow.ncread(nc, "Qold"; sel = inds)
Bw = Wflow.ncread(nc, "Bw"; sel = inds)
waterlevel = Wflow.ncread(nc, "waterlevel"; sel = inds)
DCL = Wflow.ncread(nc, "DCL"; sel = inds)
close(nc)

# create the directed acyclic graph from the drainage direction array
graph = Wflow.flowgraph(ldd, inds, Wflow.pcrdir)
# a topological sort is used for visiting nodes in order from upstream to downstream
toposort = topological_sort_by_dfs(graph)
sink = toposort[end]
@test ldd[sink] == 5  # the most downstream node must be a sink

# calculate parameters of kinematic wave
const q = 0.000001
const β = 0.6
const AlpPow = (2.0 / 3.0) * β
AlpTermR = (N ./ sqrt.(slope)) .^ β
P = Bw + (2.0 * waterlevel)
α = AlpTermR .* P .^ AlpPow

Q = zeros(n)
Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt_sec)

@testset "flow rate" begin
    @test sum(Q) ≈ 2.957806043289641e6
    @test Q[toposort[1]] ≈ 0.00021453683608501235
    @test Q[toposort[n-100]] ≈ 4054.9507466731234
    @test Q[sink] ≈ 4131.101474418251
end


text:
using ControlSystems
using MonteCarloMeasurements
using RobustAndOptimalControl.Weights
w = neglected_delay(1)
@test dcgain(w)[] < 1e-6
@test evalfr(w, 10000im)[] ≈ 2 atol=1e-3

w = gain_and_delay_uncertainty(1, 2, 1)

w = Weights.makeweight(0.1, 2)
@test dcgain(w)[] ≈ 0.1
@test evalfr(w, 10000im)[] ≈ 2 atol=1e-3

w = neglected_lag(1)
@test dcgain(w)[] < 1e-6
@test evalfr(w, 10000im)[] ≈ 1 atol=1e-3

P = tf(1 ± 0.1, [1, 1 ± 0.2, 1])
w = 2π .* exp10.(LinRange(-1, 1, 200))
centers, radii = fit_complex_perturbations(P, w; relative=true, nominal=:mean)

@test 20 <= argmax(radii) <= 22

@test maximum(radii) > 1.4
@test minimum(radii) < 0.5

##
w = 2π .* exp10.(LinRange(-2, 1, 200))
centers, radii = fit_complex_perturbations(P, w; relative=false, nominal=:mean)
nyquistplot(P, w, ylims=(-2,2), xlims=(-2,2))
nyquistcircles!(w, centers, radii)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Add <: Edit
        col::Symbol
        val::Any
    end

Add new column `col` filled with `val`

# Arguments
- `col::Symbol`: name of new column
- `val::Any`: value to add to new column
"""
mutable struct Add <: Edit
    "name of new column"
    col::Symbol
    "value to add to new column"
    val::Any
end

function Add(; col, val, )
    Add(col, val, )
end

"""Get [`Add`](@ref) `col`."""
get_col(value::Add) = value.col
"""Get [`Add`](@ref) `val`."""
get_val(value::Add) = value.val

"""Set [`Add`](@ref) `col`."""
set_col!(value::Add, val) = value.col = val
"""Set [`Add`](@ref) `val`."""
set_val!(value::Add, val) = value.val = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
function process_sweeps(s::Sweeps)
  return (;
    nsweeps=s.nsweep, maxdim=s.maxdim, mindim=s.mindim, cutoff=s.cutoff, noise=s.noise
  )
end

function tdvp(H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...)
  return tdvp(H, t, psi0; process_sweeps(sweeps)..., kwargs...)
end

function tdvp(solver, H, t::Number, psi0::MPS, sweeps::Sweeps; kwargs...)
  return tdvp(solver, H, t, psi0; process_sweeps(sweeps)..., kwargs...)
end

function dmrg(H, psi0::MPS, sweeps::Sweeps; kwargs...)
  return dmrg(H, psi0; process_sweeps(sweeps)..., kwargs...)
end


text:
immutable FIRFiltType
    kind::Cint
end

const unknown   = FIRFiltType(0 )   # LIQUID_FIRFILT_UNKNOWN    unknown filter type
const kaiser    = FIRFiltType(1 )   # LIQUID_FIRFILT_KAISER     Nyquist Kaiser filter
const pm        = FIRFiltType(2 )   # LIQUID_FIRFILT_PM         Parks-McClellan filter
const rcos      = FIRFiltType(3 )   # LIQUID_FIRFILT_RCOS       raised-cosine filter
const fexp      = FIRFiltType(4 )   # LIQUID_FIRFILT_FEXP       flipped exponential
const fsech     = FIRFiltType(5 )   # LIQUID_FIRFILT_FSECH      flipped hyperbolic secant
const farcsech  = FIRFiltType(6 )   # LIQUID_FIRFILT_FARCSECH   flipped arc-hyperbolic secant
const arkaiser  = FIRFiltType(7 )   # LIQUID_FIRFILT_ARKAISER   root-Nyquist Kaiser (approximate optimum)
const rkaiser   = FIRFiltType(8 )   # LIQUID_FIRFILT_RKAISER    root-Nyquist Kaiser (true optimum)
const rrc       = FIRFiltType(9 )   # LIQUID_FIRFILT_RRC        root raised-cosine
const hm3       = FIRFiltType(10)   # LIQUID_FIRFILT_hM3        harris-Moerder-3 filter
const gmsktx    = FIRFiltType(11)   # LIQUID_FIRFILT_GMSKTX     GMSK transmit filter
const gmskrx    = FIRFiltType(12)   # LIQUID_FIRFILT_GMSKRX     GMSK receive filter
const rfexp     = FIRFiltType(13)   # LIQUID_FIRFILT_RFEXP      flipped exponential
const rfsech    = FIRFiltType(14)   # LIQUID_FIRFILT_RFSECH     flipped hyperbolic secant
const rfarcsech = FIRFiltType(15)   # LIQUID_FIRFILT_RFARCSECH  flipped arc-hyperbolic secant

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
import SymPy.fn_map
SymPy.fn_map["Mul"] = :__prod__

function Quote(ex)
    if ex isa Real
        return ex
    else
        return convert(Expr, ex)
    end
end

function QuoteFn(name,ex,vars=free_symbols(ex))
    quot= Expr(:function, Expr(:call, Symbol(name), map(Symbol,vars)...), Quote(ex))
    if name == "_Eh"
        println(quot)
    end
    return quot
end

function QuoteFnCSE(name,cse_ex,vars=free_symbols(cse_ex[2]))
    quot = Expr(:block)
    quot.head = :function
    push!(quot.args, :($(Symbol(name))($(map(Symbol,vars)...) )))
    fnbody = QuoteCSE(cse_ex)
    push!(quot.args,fnbody)
    return quot
end

function QuoteCSE(cse_ex)
    cse,ex = cse_ex
    ex=ex[1]
    fnbody=Expr(:block)
    for (k,v) in cse
        push!(fnbody.args, :($(Quote(k)) = $(Quote(v))) )
    end
    push!(fnbody.args,Quote(ex))
    return fnbody
end

function QuoteArrCSE(cse_ex)
    # assume we did simultaneous CSE on the array
    # cse is common list of substitutions
    # ex is array of cse'd expressions
    cse,ex = cse_ex
    fnbody=Expr(:block)
    for (k,v) in cse
        push!(fnbody.args, :($(Quote(k)) = $(Quote(v))) )
    end
    if length(ex) == 1
        ex = ex[1]
    end
    push!(fnbody.args,:([$(Quote.(ex)...)]))
    return fnbody
end
function QuoteArr(ex)
    fnbody=Expr(:block)
    push!(fnbody.args,:([$(Quote.(ex)...)]))
    return fnbody
end

QuoteFnArr(name,ex, vars=free_symbols(ex)) = Expr(:function,
                          Expr(:call, Symbol(name), map(Symbol,vars)...),
                         :([$(Quote.(ex)...)]))

function QuoteFnArrCSE(name, cse_ex, vars=free_symbols(ex))
    quot = Expr(:block)
    quot.head = :function
    push!(quot.args, :($(Symbol(name))($(map(Symbol,vars)...) )))
    fnbody = QuoteArrCSE(cse_ex)
    push!(quot.args,fnbody)
    return quot
end


text:
using Lint

named_pipe = ARGS[1]
notinstalled = Pkg.installed("Lint") == nothing

if notinstalled
    print(STDERR, "linter-julia-installing-lint")
    try
        Pkg.add("Lint")
    catch
        print(STDERR, "linter-julia-msg-install")
        rethrow()
    end
else
    if Pkg.installed("Lint") < v"0.3.0"
        print(STDERR, "linter-julia-updating-lint")
        try
            Pkg.update("Lint")
        catch
            print(STDERR, "linter-julia-msg-update")
            rethrow()
        end
    else # start the server
        try
            lintserver(named_pipe,"standard-linter-v2")
        catch
            print(STDERR, "linter-julia-msg-load")
            rethrow()
        end
    end
end

max_stars_repo_name: keplersj/linter-julia

max_issues_repo_name: keplersj/linter-julia

max_forks_repo_name: keplersj/linter-julia


text:
module E

export f_1, f_2, f_3

"f_1"
f_1(x) = x

"f_2"
f_2(x) = x

"f_3"
f_3(x) = x

"g_1"
g_1(x) = x

"g_2"
g_2(x) = x

"g_3"
g_3(x) = x

export T_1

"T_1"
type T_1 end

"T_2"
type T_2 end

end


text:
# Stokes
# ======
a_ST((v,p),(φ,q),μ,dΩ) = ∫( ε(φ) ⊙ σᵥ_Ωf(μ,v) - (∇⋅φ) * p + q * (∇⋅v) )dΩ
l_ST((φ,q),f,dΩ) = ∫( φ⋅f )dΩ
function a_ST_Γd((v,p),(φ,q),n,μ,γ,h,dΓ)
  ∫( γ*μ/h*φ⋅v - φ ⋅ ( n⋅σᵥ_Ωf(μ,v) - n*p ) + ( n⋅σᵥ_Ωf(μ,φ) - n*q) ⋅ v )dΓ
end
function l_ST_Γd((φ,q),n,μ,γ,h,vD,dΓ)
  ∫( γ*μ/h*φ⋅vD + ( n⋅σᵥ_Ωf(μ,φ) - n*q) ⋅ vD )dΓ
end

# Navier-Stokes
# =============
function a_NS((v,p),(vt,pt),(φ,q),μ,ρ,dΩ)
  ∫( ρ*(ϕ ⋅ vt) + ρ*( conv∘(v,∇(v))) + ε(φ) ⊙ σᵥ_Ωf(μ,v) - (∇⋅φ) * p + q * (∇⋅v) )dΩ
end
function da_NS_dx((v,p),(dv,dp),(φ,q),μ,ρ,dΩ)
  ∫( 0.0*ρ*(ϕ ⋅ dv) + ρ*( dconv∘(dv,∇(dv),v,∇(v))) + ε(φ) ⊙ (σᵥ_Ωf(μ)∘(ε(dv))) - (∇⋅φ) * dp + q * (∇⋅dv) )dΩ
end
function da_NS_dxt((dvt,dpt),(φ,q),ρ,dΩ)
  ∫( ρ*(ϕ ⋅ dvt) )dΩ
end
l_NS((φ,q),f,dΩ) = ∫( φ⋅f )dΩ

# Navier-Stokes ALE
# =================
# LHS terms
function a_NS_ALE((u,v,p),(ut,vt,pt),(ϕ,φ,q),μ,ρ,dΩ)
  ∫( φ ⋅ ( (J∘∇(u)) * ρ * vt ) +
     φ ⋅ ( (J∘∇(u)) * ρ * (conv∘( (Finv∘∇(u))⋅(v-ut), ∇(v) ) )) +
     ∇(φ) ⊙ Pᵥ_Ωf(μ,u,v) +
     (∇⋅φ) * Pₚ_Ωf(u,p) +
     q * ((J∘∇(u)) * (∇(v) ⊙ (FinvT∘∇(u)))) )dΩ
end
function a_NS_ALE_ΓD((u,v,p),(ut,vt,pt),(ϕ,φ,q),t,n,μ,γ,h,dΓ)
  ∫( 0.0*(ϕ⋅u) + γ*μ/h*(ϕ⋅ut) + γ*μ/h*(φ⋅v) +
    -(φ ⋅ (n⋅Pᵥ_Ωf(μ,u,v))) - (φ ⋅ n)*Pₚ_Ωf(u,p) +
    (n⋅Pᵥ_Ωf(μ,u,φ)) ⋅ v + (n*Pₚ_Ωf(u,q)) ⋅ v )dΓ
end

# LHS linearized terms
function da_NS_ALE_dx((u,v,p),(ut,vt,pt),(du,dv,dp),(ϕ,φ,q),μ,ρ,dΩ)
  ∫( φ ⋅ ( (dJ∘(∇(u),∇(du))) * ρ * vt ) +
     φ ⋅ ( ( (dJ∘(∇(u),∇(du))) * ρ * (conv∘( (Finv∘∇(u))⋅(v-ut), ∇(v))) ) +
           ( (J∘∇(u)) * ρ * (conv∘( (dFinv∘(∇(u),∇(du)))⋅(v-ut), ∇(v))) ) ) +
     ∇(φ) ⊙ dPᵥ_Ωf_du(μ,u,du,v) +
     (∇⋅φ) * dPₚ_Ωf_du(u,du,p) +
     φ ⋅ ( (J∘∇(u)) * ρ * (dconv∘( (Finv∘∇(u))⋅dv, ∇(dv), (Finv∘∇(u))⋅(v-ut) , ∇(v))) ) +
     ∇(φ) ⊙ Pᵥ_Ωf(μ,u,dv) +
     (∇⋅φ) * dPₚ_Ωf_dp(u,dp) +
     q * ( (dJ∘(∇(u),∇(du))) * (∇(v) ⊙ (FinvT∘∇(u))) +
           (J∘∇(u)) * (∇(v) ⊙ (dFinvT∘(∇(u),∇(du)))) ) +
     q * ( (J∘∇(u)) * (∇(dv) ⊙ (FinvT∘∇(u))) ) )dΩ
end
function da_NS_ALE_dxt((u,v,p),(dut,dvt,dpt),(ϕ,φ,q),ρ,dΩ)
  ∫( - φ ⋅ ( (J∘∇(u)) * ρ * (conv∘((Finv∘∇(u))⋅dut, ∇(v))) ) +
     φ ⋅ ( (J∘∇(u)) * ρ * dvt ) )dΩ
end
function da_NS_ALE_ΓD_dx((u,v,p),(ut,vt,pt),(du,dv,dp),(ϕ,φ,q),n,μ,γ,h,dΓ)
  dP_tensor(u,du,v,dv,p) = dPᵥ_Ωf_du(μ,u,du,v) + Pᵥ_Ωf_dv(μ,u,dv) + dPₚ_Ωf_du(u,du,p)
  dP_scalar(u,dp) = Pₚ_Ωf_dp(u,dp)
  ∫( 0.0*(ϕ⋅du) + γ*μ/h*(φ⋅dv) +
     - φ ⋅ ( n⋅dP_tensor(u,du,v,dv,p) + n*dP_scalar(u,dp) ) +
     (n⋅Pᵥ_Ωf(μ,u,φ)) ⋅ dv + (n⋅dPᵥ_Ωf_du(μ,u,du,φ)) ⋅ v +
     (n*Pₚ_Ωf(u,q)) ⋅ dv + (n⋅dPₚ_Ωf_du(u,du,q)) ⋅ v )dΓ
end
function da_NS_ALE_ΓD_dxt((dut,dvt,dpt),(ϕ,φ,q),μ,γ,h,dΓ)
  ∫( γ*μ/h*(ϕ⋅dut) )dΓ
end

# RHS terms
function l_NS_ALE((ϕ,φ,q),f,t,dΩ)
  ∫( φ⋅f(t) )dΩ
end
function l_NS_ALE_ΓD((ϕ,φ,q),t,vD,n,μ,γ,h,dΓ)
  ∫( γ*μ/h*(ϕ⋅vD(t)) + γ*μ/h*(φ⋅vD(t)) +
    (n⋅Pᵥ_Ωf(μ,u,φ)) ⋅ vD(t) + (n*Pₚ_Ωf(u,q)) ⋅ vD(t) )dΓ
end

# Monolithic mesh motion
# ======================
function a_mesh(strategy::MeshStrategy{:laplacian},(u,v),(ϕ,φ),α,dΩ)
  ∫( α * (∇(ϕ) ⊙ ∇(u)) )dΩ
end
function a_mesh(strategy::MeshStrategy{:linearElasticity},(u,v),(ϕ,φ),λ,μ,dΩ)
  ∫( ε(ϕ) ⊙ (σₘ(λ,μ)∘ε(u)) )dΩ
end
function a_mesh(strategy::MeshStrategy{:biharmonic},(w,u,v),(ψ,ϕ,φ),α₁,α₂,dΩ)
  ∫( α₁ * ( - (ψ ⋅ w) + (∇(ψ) ⊙ ∇(u)) ) + α₂ * ( (∇(ϕ) ⊙ ∇(w)) ) )dΩ
end
function a_mesh_Γi(strategy::MeshStrategy{:laplacian},(u,v),(ϕ,φ),n,α,dΓ)
  ∫( - α * (ϕ ⋅ (n.⁺⋅∇(u))) )dΓ
end
function a_mesh_Γi(strategy::MeshStrategy{:linearElasticity},(u,v),(ϕ,φ),n,λ,μ,dΓ)
  ∫( - (ϕ ⋅  (n.⁺⋅(σₘ(λ,μ)∘ε(u))) ) )dΓ
end
function a_mesh_Γi(strategy::MeshStrategy{:biharmonic},(w,u,v),(ψ,ϕ,φ),n,α₁,α₂,dΓ)
  ∫( - α₁ * (ψ ⋅  (n.⁺⋅∇(u))) - α₂ * (ϕ ⋅  (n.⁺⋅∇(w))) )dΓ
end
function l_mesh(strategy::MeshStrategy,(ϕ,φ),f,t,dΩ)
  ∫( ϕ⋅f(t) )dΩ
end
function l_mesh(strategy::MeshStrategy{:biharmonic},(ψ,ϕ,φ),f,t,dΩ)
  ∫( ψ⋅f(t) )dΩ
end


text:
module TestGaussianDistributions
using Test
using LinearAlgebra
using GaussianDistributions

@testset "basic properties" begin
    @testset "univariate" begin
        μ, σ = randn(), exp(randn())
        G = Gaussian(μ, σ^2)
        @test G isa Univariate
        @test mean(G) == μ
        @test cov(G) == σ^2
        @test var(G) == σ^2
        @test std(G) == σ
        # indexing
        @test length(G) == 1
        @test G[1] == G # equality
        @test G[1] ≈ G
    end

    @testset "multivariate" begin
        d = 3
        μ = randn(d)
        Σ = randn(d, d)
        Σ = Σ'Σ
        G = Gaussian(μ, Σ)

        @test G isa Multivariate
        @test mean(G) == μ
        @test cov(G) == Σ
        @test var(G) == diag(Σ)
        @test std(G) == sqrt.(diag(Σ))

        # indexing
        @test length(G) == d
        for i in 1:d
            @test G[i] isa Univariate
            @test G[i] == Gaussian(μ[i], Σ[i, i])
        end
        @test G[1:2] isa Multivariate
        @test G[1:2] == Gaussian(μ[1:2], Σ[1:2, 1:2])

        # convenience constructor for diagonal covariance
        σ² = exp(randn())
        G = Gaussian(μ, σ²)
        @test G isa Multivariate
        @test mean(G) ≈ μ
        @test cov(G) ≈ σ² * I(d)
    end

    @testset "process" begin
        n = 3
        μ(x) = zero(x)
        Σ(x, y) = exp(-(x-y)^2/2)
        G = Gaussian(μ, Σ)
        @test !(G isa Univariate)
        @test !(G isa Multivariate)
        @test mean(G) == μ
        @test cov(G) == Σ

        x = randn()
        @test var(G)(x) == Σ(x, x)
        @test std(G)(x) == sqrt(Σ(x, x))

        Gx = G(x)
        @test Gx isa Univariate
        @test mean(Gx) == μ(x)
        @test var(Gx) == Σ(x, x)

        x = randn(3)
        Gx = marginal(G, x)
        @test Gx isa Multivariate
        @test mean(Gx) == μ.(x)
        @test cov(Gx) == Σ.(x, x')
        @test std(Gx) == sqrt.(Σ.(x, x))
    end
end

end


text:
export
ruleSPDotProductOutNPG,
ruleSPDotProductOutNGP,
ruleSPDotProductIn2GPN,
ruleSPDotProductIn1GNP

function ruleSPDotProductOutNPG(msg_out::Nothing,
                                msg_in1::Message{PointMass, Multivariate},
                                msg_in2::Message{F, Multivariate}) where F<:Gaussian

    x = msg_in1.dist.params[:m]
    β = convert(ProbabilityDistribution{Multivariate, GaussianMeanVariance}, msg_in2.dist)
    out_m = x' * β.params[:m]
    out_v = x' * β.params[:v] * x

    return Message(Univariate, GaussianMeanVariance, m=out_m, v=out_v)

end

function ruleSPDotProductOutNGP(msg_out::Nothing,
                                msg_in1::Message{F, Multivariate},
                                msg_in2::Message{PointMass, Multivariate}) where F<:Gaussian

    ruleSPDotProductOutNPG(msg_out, msg_in2, msg_in1)
end

function ruleSPDotProductIn2GPN(msg_out::Message{F, Univariate},
                                msg_in1::Message{PointMass, Multivariate},
                                msg_in2::Nothing) where F<:Gaussian

    x = msg_in1.dist.params[:m] # We'll call in1 x
    d = length(x)

    y = convert(ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}, msg_out.dist)
    xi = x * y.params[:xi]
    w = x * y.params[:w] * x'
    w += tiny*diageye(size(w)[1]) # Ensure w is invertible

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

function ruleSPDotProductIn1GNP(msg_out::Message{F, Univariate},
                                msg_in1::Nothing,
                                msg_in2::Message{PointMass, Multivariate}) where F<:Gaussian

    ruleSPDotProductIn2GPN(msg_out, msg_in2, msg_in1)
end


text:
monitors_default() = Vector{Vector{Symbol}}([])

"""
Default callback function for `gibbs` does nothing. You can provide your own
callback which has the same signature. One use for this is computing summary
statistics like log-likelihood after each MCMC iteration. If summary statistics
are returned, they will be appended and returned at the end of the MCMC.

For example,
```julia
function callback_fn(state, iter, pbar)
  out = nothing
  if iter % 10 == 0
    ll = round(compute_loglike(state), digits=3)
    set_description(pbar, "loglike: \$(ll)")
    out = Dict(:ll => ll)
  end

  return out
end
```
"""
function default_callback_fn(state::T, iter::Int, pbar::ProgressBar) where T
  return
end

"""
Checks if a field is a subtype.
If a field has the form "a__b", the name
of the field is `a`, and it has a field `b`.
"""
issubtype(x::Symbol)::Bool = occursin("__", String(x))

function deepcopyFields(state::T, fields::Vector{Symbol}) where T
  substate = Dict{Symbol, Any}()

  # Partition fields into subtypes and regular fields
  subtypefields, regfields = partition(issubtype, fields)

  # Get level1 fields
  for field in regfields
    substate[field] = deepcopy(getfield(state, field))
  end

  # Get level2 fields
  for field in subtypefields
    statename, _field = split(String(field), "__")
    _state = getfield(state, Symbol(statename))
    substate[field] = deepcopy(getfield(_state, Symbol(_field)))
  end

  return substate
end

showtime() = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")

function gibbs(init::T, update!::Function;
               monitors::Vector{Vector{Symbol}}=monitors_default(),
               nsamps::Vector{Int}=[1000], nburn::Int=0,
               thin::Int=1, callback_fn::Function=default_callback_fn,
               verbose::Int=1) where T
  # Assert that the first monitor is the primary one.
  # The first sample should collect the most samples.
  @assert issorted(nsamps, rev=true)

  # Make a copy of the initial state.
  state = deepcopy(init)

  # Checking number of monitors.
  num_monitors = length(monitors)
  println("Number of monitors: $(num_monitors)"); flush(stdout)
  @assert num_monitors == length(nsamps)

  # Total number of MCMC iterations
  nmcmc = nburn + thin * nsamps[1]
  verbose > 0 && println("Total number of MCMC iterations: $(nmcmc)")

  # Thin factor for each monitor
  thins = div.(nmcmc - nburn, nsamps)
  println("Thinning factors: ", thins)

  # Check monitor
  if num_monitors == 0
    verbose > 0 && println("Using default monitor.")
    fnames = fieldnames(typeof(init))
    append!(monitors, [fnames])
    num_monitors = 1
  end

  # Ensure parameters unique in each monitor.
  monitors = unique.(monitors)

  verbose > 0 && println("Preallocating memory...")

  # Create object to return
  @time chain = [fill(deepcopyFields(state, monitors[i]), nsamps[i])
                 for i in 1:num_monitors]

  # Initialize empty summary statistics.
  summary_stats = []

  # Initialize sample index tracker for each monitor.
  counters = zeros(Int, num_monitors)

  # Initialize progress bar.
  pbar = ProgressBar(1:nmcmc)

  # Print time at beginning of MCMC.
  println("Start time: ", showtime()); flush(stdout)

  # Gibbs loop.
  for i in pbar
    update!(state)
    if i > nburn
      for j in 1:num_monitors
        if i % thins[j] == 0
          substate = deepcopyFields(state, monitors[j])
          counters[j] += 1
          chain[j][counters[j]] = substate
        end
      end
    end

    # Execute callback function.
    summary_stat = callback_fn(state, i, pbar)

    # Append summary statistics if callback returns something.
    if summary_stat != nothing
      append!(summary_stats, [summary_stat])
    end

    # This is needed so that `/usr/bin/tail` in works properly
    # in a BASH terminal.
    flush(stdout)
  end

  # Print end time.
  println("End time: ", showtime()); flush(stdout)

  # Return chain, last state, and summary statistics.
  return (chain, state, summary_stats)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module Graph

using DocStringExtensions

import ..Delayed

"""
$(TYPEDSIGNATURES)

Given an delayed object, return an array of argument values.
"""
function delayedproperties(item::Delayed.AbstractDelayed)::Vector

    return [
        getproperty(item, propertyname)
            for propertyname in propertynames(item)
    ]

end
delayedproperties(a::Delayed.AbstractSwitch) = [a.choice, a.options...]
delayedproperties(a::Delayed.CategoricalIndex) = a.probabilities

"""
$(TYPEDSIGNATURES)

Depth-first search
Unrolls a graph into an array of all the nodes.
"""
dfs(space::Dict{Symbol, T} where T)::Vector = dfs!(Delayed.AbstractDelayed[], space)

"""
$(TYPEDSIGNATURES)

Unrolls a nested space into a vector of all the Delayed.AbstractDelayed nodes.
"""
function dfs!(seq::Vector{Delayed.AbstractDelayed}, item::Delayed.AbstractDelayed)::Vector
    # For every input (arg) of the object, add those to the list too.
    for prop in delayedproperties(item)
        dfs!(seq, prop)
    end
    push!(seq, item)

    return seq
end
function dfs!(seq::Vector{Delayed.AbstractDelayed}, item::Dict)::Vector
    for v in values(item)
        dfs!(seq, v)
    end

    return seq
end
function dfs!(seq::Vector{Delayed.AbstractDelayed}, item::Union{Tuple, Vector})::Vector
    for v in item
        dfs!(seq, v)
    end

    return seq
end
dfs!(seq::Vector{Delayed.AbstractDelayed}, item::Any)::Vector = seq

function checklabel!(labels::Vector{Symbol}, node::Delayed.AbstractParam)::Nothing
    if node.label in labels
        throw(DomainError(node.label, "Your space has a duplicate label :$(node.label)"))
    end
    push!(labels, node.label)

    return nothing
end
checklabel!(labels::Vector{Symbol}, node::T where T)::Nothing = nothing

"""
$(TYPEDSIGNATURES)

Ensure the user hasn't submitted any duplicate labels in their space.
"""
function checkspace(space::Dict{Symbol, T})::Dict{Symbol, T} where T
    labels = Symbol[]
    for node in Graph.dfs(space)
        checklabel!(labels, node)
    end

    return space
end

end # module Graph


text:
using LinearAlgebra
using MixedModels
using Random
using SparseArrays
using StatsModels
using Test

const LMM = LinearMixedModel

@testset "UBlk" begin
    ex22 = UniformBlockDiagonal(reshape(Vector(1.0:12.0), (2, 2, 3)))
    Lblk = UniformBlockDiagonal(fill(0., (2,2,3)))
    ds = (Y = rand(12), A = repeat(['N','Y'], outer=6), G = repeat('a':'c', inner=4),
        H = repeat('A':'B', outer=6), U = repeat([-1,0,1], inner=2, outer=2))
    sch = schema(ds, Dict(:A=>EffectsCoding()))
    vf1 = modelcols(apply_schema(@formula(Y ~ 1 + A + (1+A|G)), sch, LMM), ds)[2][2]
    vf2 = modelcols(apply_schema(@formula(Y ~ 1 + U + (1+U|H)), sch, LMM), ds)[2][2]
    prd = vf2'vf1

    @testset "size" begin
        @test size(ex22) == (6, 6)
        @test size(ex22, 1) == 6
        @test size(ex22, 2) == 6
        @test size(ex22.data) == (2, 2, 3)
    #    @test length(ex22.facevec) == 3
        @test size(vf1) == (12, 6)
        @test size(vf2) == (12, 4)
        @test size(prd) == (4, 6)
    end

    @testset "elements" begin
        @test ex22[1, 1] == 1
        @test ex22[2, 1] == 2
        @test ex22[3, 1] == 0
        @test ex22[2, 2] == 4
        @test ex22[3, 3] == 5
        @test ex22[:, 3] == [0,0,5,6,0,0]
        @test ex22[5, 6] == 11
    end

    @testset "facevec" begin
        @test view(ex22.data, :, :, 3) == reshape(9:12, (2,2))
    end

    @testset "scaleinflate" begin
        MixedModels.scaleinflate!(copyto!(Lblk, ex22), vf1)
        @test view(Lblk.data, :, :, 1) == [2. 3.; 2. 5.]
        setθ!(vf1, [1.,1.,1.])
        Λ = vf1.λ
        MixedModels.scaleinflate!(copyto!(Lblk, ex22), vf1)
        target = Λ'view(ex22.data, :, :, 1)*Λ + I
        @test view(Lblk.data, :, :, 1) == target
    end

    @testset "updateL" begin
        @test ones(2, 2) == MixedModels.rankUpdate!(Hermitian(zeros(2, 2)), ones(2), 1., 1.)
        d3 = MixedModels.dataset(:d3)
        sch = schema(d3)
        vf1 = modelcols(apply_schema(@formula(y ~ 1 + u + (1+u|g)), sch, LMM), d3)[2][2]
        vf2 = modelcols(apply_schema(@formula(y ~ 1 + u + (1+u|h)), sch, LMM), d3)[2][2]
        @test vf1.λ == LowerTriangular(Matrix(I, 2, 2))
        setθ!(vf2, [1.75, 0.0, 1.0])
        A11 = vf1'vf1
        L11 = MixedModels.cholUnblocked!(MixedModels.scaleinflate!(copyto!(UniformBlockDiagonal(fill(0., size(A11.data))), A11), vf1), Val{:L})
        L21 = vf2'vf1
        @test isa(L21, BlockedSparse)
        @test L21[1,1] == 30.0
        @test size(L21) == (344, 9452)
        @test size(L21, 1) == 344
        MixedModels.lmulΛ!(vf2', MixedModels.rmulΛ!(L21, vf1))
        @test size(Matrix(L21)) == size(sparse(L21))
#        L21cb1 = copy(L21.colblocks[1])
#        @test L21cb1 == Vf2.Λ * A21cb1 * Vf1.Λ
#        rdiv!(L21, adjoint(LowerTriangular(L11)))
#        @test_broken L21.colblocks[1] == rdiv!(L21cb1, adjoint(LowerTriangular(L11.facevec[1])))
         A22 = vf2'vf2
         L22 = MixedModels.scaleinflate!(copyto!(UniformBlockDiagonal(fill(0., size(A22.data))), A22), vf2)
    end
end


text:
using Weave

weave("$(@__DIR__)/three_compartment_model/three_compartment_model.jmd",
    doctype = "github",
    out_path = "$(@__DIR__)/three_compartment_model/README.md")

weave("$(@__DIR__)/simple_lung/simple_lung.jmd",
    doctype = "github",
    out_path = "$(@__DIR__)/simple_lung/README.md")

weave("$(@__DIR__)/parameter_estimation/parameter_estimation.jmd",
    doctype = "github",
    out_path = "$(@__DIR__)/parameter_estimation/README.md")


text:
export vec_bs, spike_count_win

#Spike Counts

function spike_count_win(rate::SpikeTrain,t1::Float64,t2::Float64,out=zeros(Int64,length(rate.trials)))

    @inbounds for i=1:length(out)
        mycent=rate.center[i,1]
        for j in rate.trials[i].inds
            if rate.ts[j]-mycent > t1
                if rate.ts[j]-mycent < t2
                    out[i]+=1
                end
            elseif rate.ts[j]-mycent > t2
                break
            end
        end
    end

    out
end

function spike_count_win(rate::SpikeTrain,t1::Array{Float64,1},t2::Array{Float64,1},out=zeros(Float64,length(rate.trials)))

    for i=1:length(out)
        for j in rate.trials[i].inds
            if rate.ts[j] > t1[i]
                if rate.ts[j] < t2[i]
                    out[i]+=1
                end
            elseif rate.ts[j] > t2[i]
                break
            end
        end
        out[i]=out[i]/(t2[i]-t1[i])
    end

    out
end

function spike_count_win_ind(rate::SpikeTrain,t1::Float64,t2::Float64,inds::Array{Int64,1},out=zeros(Int64,length(inds)))

    @inbounds for i=1:length(out)
        mycent=rate.center[inds[i],1]
        for j in rate.trials[inds[i]].inds
            if rate.ts[j]-mycent > t1
                if rate.ts[j]-mycent < t2
                    out[i]+=1
                end
            elseif rate.ts[j]-mycent > t2
                break
            end
        end
    end

    out
end

function spike_count_win(rate::SpikeTrain,ts::Array{Float64,1},out=zeros(Int64,length(rate.trials),length(ts)-1))

    @inbounds for i=1:length(out)
        out[i]=0.0
    end
   @inbounds for i in 1:length(rate.trials)
        count=1
        mycent=rate.center[i,1]
        for j in rate.trials[i].inds
            myt=rate.ts[j]-mycent
            if myt>ts[count]
                k=count+1
                while k < length(ts)+1
                    if myt<ts[k]
                        out[i,k-1]+=1.0
                        count=k-1
                        break
                    end
                    k+=1
                end
            end
        end
    end
    out
end

#=
Modulation Index
Vector Method
=#

function vec_bs(rate::SpikeTrain,t1::Array{Float64,1},t2::Array{Float64,1},cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})
    myvec1=spike_count_win(rate,t1,t2)

    vec_bs(myvec1,cons,t_trials,t_dir,inds)
end

function vec_bs(rate::SpikeTrain,t1::Float64,t2::Float64,cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})
    # cons - trial type identifier for each trial of interest
    # t_trials - Maps trial type identifier to index in dictionary
    # t_dir - angle for each trial type
    # inds - maps index of t_types to index of s (first trial of interest (index 1 in t_types) may be the 4th absolute trial, so inds[1]=4)

    myvec1=spike_count_win(rate,t1,t2)

    vec_bs(myvec1,cons,t_trials,t_dir,inds)
end

function vec_bs{T}(myvec1::Array{T,1},cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})

    mv1=PD_vec(myvec1,cons,t_dir,t_trials,inds)

    miss=0.0
    for i=1:1000
        cons2=deepcopy(cons)
        shuffle!(cons2)
        mv2=PD_vec(myvec1,cons2,t_dir,t_trials,inds)
        if mv2>mv1
            miss+=1.0
        end
    end

    miss/1000
end

function PD_vec{T}(s::Array{T,1},cons::Array{Int64,1},t_dir::Array{Float64,1},t_trials::Dict{Int64,Int64},inds::Array{Int64,1})
    # s - Spike counts per trial
    # cons - condition identifier for each trial of interest
    # t_dir - angle for each trial type
    # t_trials - Maps trial type identifier to index in dictionary
    # inds - maps index of t_types to index of s (first trial of interest (index 1 in t_types) may be the 4th absolute trial, so inds[1]=4)

    dirs=zeros(Float64,length(t_dir))
    dirn=zeros(Int64,length(t_dir))

    for i=1:length(cons)
        dirs[t_trials[cons[i]]]+=s[inds[i]]
        dirn[t_trials[cons[i]]]+=1
    end

    #normalize for unequal trial types
    for i=1:length(dirs)
        dirs[i]=dirs[i]/dirn[i]
    end

    largest_dir=maximum(dirs)
    x=0.0
    y=0.0

    for i=1:length(dirs)
        #dirs[i]=dirs[i]/largest_dir
        x+=cosd(t_dir[i])*dirs[i]
        y+=sind(t_dir[i])*dirs[i]
    end

    sqrt(x*x+y*y)
end

function mod_bs(rate::SpikeTrain,t1::Float64,t2::Float64,cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})

    myvec1=spike_count_win(rate,t1,t2)
    mod_bs(myvec1,cons,t_trials,t_dir,inds)
end

function mod_bs{T}(myvec1::Array{T,1},cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})

    mv1=PD_mod(s,cons,t_trials,t_dir,inds)

    miss=0.0
    for i=1:1000
        cons2=deepcopy(cons)
        shuffle!(cons2)
        mv2=PD_mod(myvec1,cons2,t_trials,t_dir,inds)
        if mv2>mv1
            miss+=1.0
        end
    end

    miss/1000
end

function PD_mod{T}(s::Array{T,2},cons::Array{Int64,1},t_trials::Dict{Int64,Int64},t_dir::Array{Float64,1},inds::Array{Int64,1})

end


text:
valueof(::Val{x}) where {x} = x

# @inline fieldindex(::T, field::Val) where {T} = fieldindex(T, field)
@generated function fieldindex(::Type{T}, ::Val{field}) where {T,field}
    field isa Symbol || return :(error("`field` must be a symbol; given: $field"))
    return findfirst(n -> n === field, fieldnames(T))
end

if !isdefined(Base, :ismutabletype)
    function ismutabletype(@nospecialize(t::Type))
        t = Base.unwrap_unionall(t)
        # TODO: what to do for `Union`?
        return isa(t, DataType) && t.mutable
    end
end

function noinline(f)
    @noinline wrapper() = f()
    wrapper()
end

alignto(ptr, ::Nothing) = ptr
function alignto(ptr, align)
    ispow2(align) || noinline() do
        error("not a power of 2: `align = $align`")
    end
    # return ptr + (align - mod(UInt(ptr), align))
    return ptr + (align - (UInt(ptr) & (align - 1)))
end

struct UnionValue{T}
    value::T
end

unwrap_union_value(x) = x
unwrap_union_value(x::UnionValue) = x.value

function wrap_union_type(::Type{T}) where {T}
    if T isa Union
        return UnionValue{T}
    else
        return T
    end
end

@inline function field_convert(::Type{T}, x) where {T}
    if T isa Union
        return UnionValue{T}(convert(T, x))
    else
        return convert(T, x)
    end
end

function define_docstrings()
    docstrings = [:RecordArrays => joinpath(dirname(@__DIR__), "README.md")]
    #=
    docsdir = joinpath(@__DIR__, "docs")
    for filename in readdir(docsdir)
        stem, ext = splitext(filename)
        ext == ".md" || continue
        name = Symbol(stem)
        name in names(RecordArrays, all=true) || continue
        push!(docstrings, name => joinpath(docsdir, filename))
    end
    =#
    for (name, path) in docstrings
        include_dependency(path)
        doc = read(path, String)
        doc = replace(doc, r"^```julia"m => "```jldoctest $name")
        doc = replace(doc, "<kbd>TAB</kbd>" => "_TAB_")
        @eval RecordArrays $Base.@doc $doc $name
    end
end


text:
module Integration

using Compat
using FastGaussQuadrature

export qdata

function qdata(beta::Vector{Float64},nqpts::Int64)

  n = length(beta)
  qnode = zeros(nqpts,n+1)
  qwght = zeros(nqpts,n+1)
  for j = findall(beta.>-1)
    qnode[:,j],qwght[:,j] = gaussjacobi(nqpts,0,beta[j])
  end
  qnode[:,n+1],qwght[:,n+1] = gaussjacobi(nqpts,0,0)
  return qnode, qwght

end

end


text:
function load_mnist(batchsize::Int, transform = cpu)
    # Use MLDataUtils LabelEnc for natural onehot conversion
    onehot(labels_raw) =
        convertlabel(LabelEnc.OneOfK, labels_raw, LabelEnc.NativeLabels(collect(0:9)))
    # Load MNIST
    imgs, labels_raw = MNIST.traindata()
    imgs = permutedims(imgs, (2, 1, 3))
    # Process images into (H,W,C,BS) batches
    x_train_data = Float32.(reshape(imgs, 28, 28, 1, size(imgs, 3)))
    y_train_data = onehot(labels_raw)
    # Load MNIST Test
    imgs, labels_raw = MNIST.testdata()
    imgs = permutedims(imgs, (2, 1, 3))
    # Process images into (H,W,C,BS) batches
    x_test_data = Float32.(reshape(imgs, 28, 28, 1, size(imgs, 3)))
    y_test_data = onehot(labels_raw)
    return (
        # Use Flux's Flux.Data.DataLoader to automatically minibatch and shuffle the data
        Flux.Data.DataLoader(
            transform.(copy.((x_train_data, y_train_data)));
            batchsize = batchsize,
            shuffle = true,
        ),
        # Don't shuffle the test data
        Flux.Data.DataLoader(
            transform.(copy.((x_test_data, y_test_data)));
            batchsize = batchsize,
            shuffle = false,
        ),
    )
end

function load_miniboone(
    batchsize::Int,
    path::String,
    train_test_split::Float64 = 0.8,
    transform = cpu,
)
    data = Float32.(npzread(path)')
    data = (data .- mean(data, dims = 2)) ./ std(data, dims = 2)
    total_obs = size(data, 2)
    train_data, test_data = splitobs(shuffleobs(data), train_test_split)

    return (
        Flux.Data.DataLoader(
            transform(copy(train_data)),
            batchsize = batchsize,
            shuffle = true,
        ),
        Flux.Data.DataLoader(
            transform(copy(test_data)),
            batchsize = batchsize,
            shuffle = false,
        ),
    )
end

function load_physionet(
    batchsize::Int,
    path::String,
    train_test_split::Float64 = 0.8,
    transform = cpu,
)
    data = BSON.load(path)[:data]
    total_obs = size(data[:observed_data], 3)
    train_idx, test_idx = splitobs(shuffleobs(collect(1:total_obs)), train_test_split)
    train_data = []
    test_data = []
    for key in [:observed_data, :observed_mask, :data_to_predict, :mask_predicted_data]
        push!(train_data, data[key][:, :, train_idx])
        push!(test_data, data[key][:, :, test_idx])
    end
    for key in [:observed_tp, :tp_to_predict]
        push!(train_data, reshape(data[key][:, train_idx], 1, 49, :)),
        push!(test_data, reshape(data[key][:, test_idx], 1, 49, :))
    end
    return (
        Flux.Data.DataLoader(
            transform.(copy.(train_data))...,
            batchsize = batchsize,
            shuffle = true,
            partial = false,
        ),
        Flux.Data.DataLoader(
            transform.(copy.(test_data))...,
            batchsize = batchsize,
            shuffle = true,
            partial = false,
        ),
    )
end

function load_spiral2d(
    batchsize::Int,
    transform = cpu;
    nspiral = 1000,
    ntotal = 500,
    nsample = 100,
    start = 0.0,
    stop = 1.0,
    noise_std = 0.1,
    a = 0.0,
    b = 1.0,
)
    # Toy Spiral 2D dataset for testing regularization on time series
    # extrapolation problem
    # A 2D spiral is parameterized by `r = a + b * theta`
    orig_ts = range(start, stop, length = ntotal)
    samp_ts = orig_ts[1:nsample]

    # clockwise and counter clockwissse spirals in observation space
    zs_cw = stop .+ 1.0f0 .- orig_ts
    rs_cw = a .+ b .* 50.0f0 ./ zs_cw
    xs = reshape(rs_cw .* cos.(zs_cw) .- 5.0f0, 1, :)
    ys = reshape(rs_cw .* sin.(zs_cw), 1, :)
    orig_traj_cw = reshape(cat(xs, ys, dims = 1), 2, :, 1)

    zs_cc = orig_ts
    rw_cc = a .+ b .* zs_cc
    xs = reshape(rw_cc .* cos.(zs_cc) .+ 5.0f0, 1, :)
    ys = reshape(rw_cc .* sin.(zs_cc), 1, :)
    orig_traj_cc = reshape(cat(xs, ys, dims = 1), 2, :, 1)

    sample_trajectories = []
    original_trajectories = []
    for _ = 1:nspiral
        t₀ = rand(1:(ntotal-2*nsample)) + nsample - 1

        orig_traj = rand() > 0.5 ? orig_traj_cc : orig_traj_cw
        push!(original_trajectories, copy(orig_traj))

        samp_traj = copy(orig_traj[:, t₀:t₀+nsample-1, :])
        samp_traj .+= randn(size(samp_traj)) .* noise_std
        push!(sample_trajectories, samp_traj)
    end

    original_trajectories = Float32.(cat(original_trajectories..., dims = 3))
    original_tp = Float32.(reshape(repeat(collect(orig_ts), nspiral), :, nspiral))
    sampled_trajectories = Float32.(cat(sample_trajectories..., dims = 3))
    sampled_tp = Float32.(reshape(repeat(collect(samp_ts), nspiral), :, nspiral))

    return (
        Flux.Data.DataLoader(
            transform.((sampled_trajectories, sampled_tp)),
            batchsize = batchsize,
            shuffle = true,
        ),
        Flux.Data.DataLoader(
            transform.((original_trajectories, original_tp)),
            batchsize = batchsize,
            shuffle = true,
        ),
    )
end

function load_gaussian_mixture(
    batchsize,
    transform = cpu,
    train_test_split = 0.75;
    nsamples = 1000,
    ngaussians = 6,
    dim = 2,
    radius = 5.0f0,
    σ = 0.1f0,
    noise = 0.3f0,
)
    samples_per_gaussian = nsamples ÷ ngaussians
    μ = zeros(Float32, dim)
    θ = 0
    X = Array{Float32}(undef, 2, samples_per_gaussian * ngaussians)
    for i = 1:ngaussians
        θ += Float32(2π / ngaussians)
        μ[1] = cos(θ) * radius
        μ[2] = sin(θ) * radius

        dist = MvNormal(μ, σ)
        samples = rand(dist, samples_per_gaussian)
        noise_vec = Float32.(randn(2, samples_per_gaussian)) .* noise
        X[:, (i-1)*samples_per_gaussian+1:i*samples_per_gaussian] = samples + noise_vec
    end

    X_train, X_test = splitobs(shuffleobs(X), train_test_split)

    return (
        Flux.Data.DataLoader(
            transform(copy(X_train)),
            batchsize = batchsize,
            shuffle = true,
        ),
        Flux.Data.DataLoader(
            transform(copy(X_test)),
            batchsize = batchsize,
            shuffle = false,
        ),
    )
end


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

"""
    Line(a, b)

A line passing through points `a` and `b`.

See also [`Segment`](@ref).
"""
struct Line{Dim,T} <: Primitive{Dim,T}
  a::Point{Dim,T}
  b::Point{Dim,T}
end

Line(a::Tuple, b::Tuple) = Line(Point(a), Point(b))

paramdim(::Type{<:Line}) = 1

(l::Line)(t) = l.a + t * (l.b - l.a)


text:
using Documenter
using Catalyst, ModelingToolkit

include("pages.jl")

# supposed hack to get mhchem but doesn't seem to work yet...
# const katex_version = "0.11.1"
# function Documenter.Writers.HTMLWriter.RD.mathengine!(r::Documenter.Utilities.JSDependencies.RequireJS, engine::Documenter.KaTeX)
#     push!(r, Documenter.Utilities.JSDependencies.RemoteLibrary(
#         "katex",
#         "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/$(katex_version)/katex.min.js"
#     ))
#     push!(r, Documenter.Utilities.JSDependencies.RemoteLibrary(
#         "mhchem",
#         "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/$(katex_version)/contrib/mhchem.min.js"
#     ))
#     push!(r, Documenter.Utilities.JSDependencies.RemoteLibrary(
#         "katex-auto-render",
#         "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/$(katex_version)/contrib/auto-render.min.js",
#         deps = ["katex"],
#     ))
#     push!(r, Documenter.Utilities.JSDependencies.Snippet(
#         ["jquery", "katex", "mhchem", "katex-auto-render"],
#         ["\$", "katex", "mhchem", "renderMathInElement"],
#         """
#         \$(document).ready(function() {
#           renderMathInElement(
#             document.body,
#             $(Documenter.Utilities.JSDependencies.json_jsescape(engine.config, 2))
#           );
#         })
#         """
#     ))
# end

makedocs(
    sitename = "Catalyst.jl",
    authors = "Samuel Isaacson",
    format = Documenter.HTML(analytics = "UA-90474609-3", mathengine=Documenter.Writers.HTMLWriter.MathJax(), prettyurls = (get(ENV, "CI", nothing) == "true")),
    modules = [Catalyst,ModelingToolkit],
    doctest = false,
    clean = true,
    pages = pages
)

deploydocs(
   repo = "github.com/SciML/Catalyst.jl.git";
   push_preview = true
)


text:
using OptionalArgChecks
using Documenter

DocMeta.setdocmeta!(OptionalArgChecks, :DocTestSetup, :(using OptionalArgChecks); recursive=true)

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

deploydocs(;
    repo="github.com/simeonschaub/OptionalArgChecks.jl",
)


text:
data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
    datadep"CiteSeer"
end

@testset "CiteSeer" begin
    data  = CiteSeer.dataset()

    @test data.num_nodes == 3327
    @test data.num_edges == 9104
    @test data.directed == false
    @test size(data.node_features) == (3703, data.num_nodes)
    @test size(data.node_labels) == (data.num_nodes,)
    @test size(data.train_indices) == (120,)
    @test size(data.val_indices) == (500,)
    @test size(data.test_indices) == (1015,)
    @test size(data.adjacency_list) == (data.num_nodes, )
    @test sum(length.(data.adjacency_list)) == (data.num_edges)
    if VERSION >= v"1.6.0"
        @test minimum(minimum.(data.adjacency_list; init=1000)) == 1
        @test maximum(maximum.(data.adjacency_list; init=1000)) == data.num_nodes
    end
end


text:
function normalize_node(node::R,domain::Array{T,1}) where {R<:Number,T<:AbstractFloat}

  if domain[1] == domain[2]
    norm_node = zero(T)
    return norm_node
  else
    norm_node = 2.0*(node-domain[2])/(domain[1]-domain[2])-1.0
    return norm_node
  end

end

function normalize_node(node::AbstractArray{R,1},domain::Array{T,1}) where {R<:Number,T<:AbstractFloat}

  norm_nodes = similar(node)
  for i in eachindex(node)
    norm_nodes[i] = normalize_node(node[i],domain)
  end

  return norm_nodes

end

function normalize_node!(node::AbstractArray{R,1},domain::Array{T,1}) where {R<:Number,T<:AbstractFloat}

  for i in eachindex(node)
    node[i] = normalize_node(node[i],domain)
  end

end

max_stars_repo_name: RJDennis/ChebyshevApprox

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
function get_vector_info(vecname, maxlen=Int(maxintfloat()))
    factor = 1.0
    if occursin("*", vecname)
        factor, vecname = split(vecname, "*")
        factor = parse(Float64, factor)
    end
    vec = ngGet_Vec_Info(vecname)
    vec != C_NULL || throw("Vector $(vecname) not found")
    vecinfo = unsafe_load(vec)
    vname  = unsafe_string(vecinfo.name)
    len = min(maxlen, vecinfo.length)
    typelist = Dict(0 => "notype",
                    1 => "time",
                    2 => "frequency",
                    3 => "voltage",
                    4 => "current",
                    5 => "voltage density",
                    6 => "current density",
                    7 => "sqr voltage density",
                    8 => "sqr current density",
                    9 => "sqr voltage",
                    10 => "sqr current",
                    11 => "pole",
                    12 => "zero",
                    13 => "sparam",
                    14 => "temp",
                    15 => "res",
                    16 => "impedence",
                    17 => "admittance",
                    18 => "power",
                    19 => "phase",
                    20 => "db",
                    21 => "capacitance",
                    22 => "charge",
                    )
    vtype = get(typelist, vecinfo.type, "-")
    if (vecinfo.flags & VF_REAL) != 0
        vreal = copy(unsafe_wrap(Array, vecinfo.realdata, (len,)))
        return vname, vtype, vreal*factor
    elseif (vecinfo.flags & VF_COMPLEX) != 0
        vcomplex = copy(unsafe_wrap(Array, vecinfo.compdata, (len,)))
        return vname, vtype, vcomplex*factor
    else
        error("Unknown vector type")
    end
end

function curplot()
    cur = ngSpice_CurPlot()
    cur != C_NULL || throw("No current plots")
    unsafe_string(cur)
end

function listallplots(nplots=Int(maxintfloat()))
    ppallplots = ngSpice_AllPlots()
    pallplots = unsafe_wrap(Array, ppallplots, nplots)
    allplots = String[]
    for pplot in pallplots
        pplot != C_NULL || return allplots
        plt = unsafe_string(pplot)
        push!(allplots, plt)
    end
end

function listcurvecs(plot=curplot(), nvecs=Int(maxintfloat()))
    ppcurvecs = ngSpice_AllVecs(plot)
    pcurvecs = unsafe_wrap(Array, ppcurvecs, nvecs)
    curvecs = String[]
    for pvec in pcurvecs
        pvec != C_NULL || return curvecs
        v = unsafe_string(pvec)
        push!(curvecs, v)
    end
end

function listallvecs(lplot=listallplots())
    allvecs = Dict(lplot .=> listcurvecs.(lplot))
    return allvecs
end

function getvec(name, maxlen=Int(maxintfloat()))
    get_vector_info(name, maxlen)
end

function getrealvec(name, maxlen=Int(maxintfloat()))
    _, __, data = get_vector_info(String(name), maxlen)
    real(data)
end

function getimagvec(name, maxlen=Int(maxintfloat()))
    _, __, data = get_vector_info(name, maxlen)
    imag(data)
end

function getmagnitudevec(name, maxlen=Int(maxintfloat()))
    _, __, data = get_vector_info(name, maxlen)
    abs.(data)
end

function getphasevec(name, maxlen=Int(maxintfloat))
    _, __, data = get_vector_info(name, maxlen)
    angle.(data)
end

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
using Documenter, PhyloNetworks

# Weave the .jmd
include(Pkg.dir("PhyloNetworks","docs","src", "man", "src", "make_weave.jl"))

makedocs()

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo = "github.com/crsl4/PhyloNetworks.jl.git",
    julia  = "0.6",
    osname = "linux"
)


text:
abstract type IntervalEstimator end

Base.iterate(estimator::IntervalEstimator, state = 0) = state > 0 ? nothing : (estimator, state + 1)
Base.length(estimator::IntervalEstimator) = 1

Base.show(io::IO, estimator::IntervalEstimator) = print(io, string(estimator))
Base.show(io::IO, ::MIME"application/prs.juno.inline", estimator::IntervalEstimator) = print(io, string(estimator))

Base.string(estimator::IntervalEstimator) = string(typeof(estimator))

function get_bounds(estimator::IntervalEstimator, x1::TI, x2::TI) where {TI<:Integer}

    for i in 1:size(estimator.XX, 1)
        if [x1, x2] == estimator.XX[i,:]
            return estimator.bounds[i, :]
        end
    end
    error("(x1,x2) not found in sample space, valid observation?")
end
function get_bounds(estimator::IntervalEstimator, x1::Vector{TI}, x2::Vector{TI}) where {TI<:Integer}

    return convert(Matrix{eltype(estimator.bounds)}, hcat(get_bounds.(estimator, x1, x2)...)')
end

(estimator::IntervalEstimator)(x1::TI, x2::TI) where {TI<:Integer} = get_bounds(estimator, x1, x2)
(estimator::IntervalEstimator)(x1::Vector{TI}, x2::Vector{TI}) where {TI<:Integer} = get_bounds(estimator, x1, x2)

function compatible(estimator::IntervalEstimator, design::AbstractDesign, pnull::Real)

    @assert (0 <= pnull <= 1) "pnull must be in [0,1]"

    XX                      = sample_space(design)
    design_rejects          = reject.(XX[:,1], XX[:,2], design)
    lower_boundaries        = get_bounds(estimator, XX[:,1], XX[:,2])[:,1]
    iv_rejects              = lower_boundaries .> pnull
    incompatibility_degree  = sum(design_rejects .& .!iv_rejects)

    df = DataFrames.DataFrame(
        x1 = XX[:,1],
        x2 = XX[:,2],
        design_rejects = design_rejects,
        ci_reject = iv_rejects,
        lower_boundaries = lower_boundaries
    )

    return Dict{String,Any}(
        "compatible" => incompatibility_degree == 0,
        "incompatibility degree" => incompatibility_degree,
        "details" => df
    )
end

function coverage_probability(estimator::IntervalEstimator, p::Real)

    x1, x2 = estimator.XX[:,1], estimator.XX[:,2]
    pmff   = pmf.(x2, n2.(estimator.design, x1), p) .* pmf.(x1, n1(estimator.design), p)
    lower  = sum(pmff[estimator.bounds[:,1] .<= p])
    upper  = sum(pmff[estimator.bounds[:,2] .>= p])
    joint  = sum(pmff[(estimator.bounds[:,1] .<= p) .& (estimator.bounds[:,2] .>= p)])
    return [lower, joint, upper]
end

function coverage_probability(estimator::IntervalEstimator, p::Vector{TR}) where {TR<:Real}

    return convert(Matrix{TR}, hcat(coverage_probability.(estimator, p)...)')
end

function mean_width(estimator::IntervalEstimator, p::Real)

    x1, x2 = estimator.XX[:,1], estimator.XX[:,2]
    pmff   = pmf.(x2, n2.(estimator.design, x1), p) .* pmf.(x1, n1(estimator.design), p)
    return sum( (estimator.bounds[:,2] .- estimator.bounds[:,1]) .* pmff )
end


text:
using StructArrays: StructArray
using ComponentArrays: ComponentArrays
using LinearAlgebra: diagm, mul!
using NQCDynamics: RingPolymers
using NQCModels: nstates
using RingPolymerArrays: get_centroid
using NQCDistributions: ElectronicDistribution, Adiabatic, Diabatic, density_matrix
using ..Calculators: AbstractDiabaticCalculator, DiabaticCalculator, RingPolymerDiabaticCalculator

function set_quantum_derivative! end

function calculate_density_matrix_propagator!(sim::Simulation, v)
    V = sim.method.density_propagator

    V .= diagm(sim.calculator.eigen.values)
    for I in eachindex(v)
        @. V -= im * v[I] * sim.calculator.nonadiabatic_coupling[I]
    end
    return V
end

function calculate_density_matrix_propagator!(sim::RingPolymerSimulation, v)
    V = sim.method.density_propagator
    centroid_v = get_centroid(v)

    V .= diagm(sim.calculator.centroid_eigen.values)
    for I in eachindex(centroid_v)
        @. V -= im * centroid_v[I] * sim.calculator.centroid_nonadiabatic_coupling[I]
    end
    return V
end

function commutator!(out, A, B, tmp)
    mul!(out, A, B)
    mul!(tmp, B, A)
    out .-= tmp
    return nothing
end

get_quantum_subsystem(u::ComponentArrays.ComponentVector{T}) where {T} =
    StructArray{Complex{T}}((u.σreal, u.σimag))

function initialise_adiabatic_density_matrix(
    electronics::ElectronicDistribution{Diabatic},
    calculator::AbstractDiabaticCalculator,
    r
)

    diabatic_density = density_matrix(electronics, nstates(calculator))
    return transform_density!(diabatic_density, calculator, r, :to_adiabatic)
end

function initialise_adiabatic_density_matrix(
    electronics::ElectronicDistribution{Adiabatic},
    calculator::AbstractDiabaticCalculator,
    r
)

    adiabatic_density = density_matrix(electronics, nstates(calculator))
    return adiabatic_density
end

function transform_density!(
    density::AbstractMatrix, calculator::AbstractDiabaticCalculator, r, direction
)
    U = evaluate_transformation(calculator, r)
    if direction === :to_diabatic
        U = U'
    elseif !(direction === :to_adiabatic)
        throw(ArgumentError("`direction` $direction not recognised."))
    end
    density .= U' * density * U
    return density
end

function evaluate_transformation(calculator::DiabaticCalculator, r)
    eigs =  Calculators.get_eigen(calculator, r)
    return eigs.vectors
end

function evaluate_transformation(calculator::RingPolymerDiabaticCalculator, r)
    centroid_eigs = Calculators.get_centroid_eigen(calculator, r)
    return centroid_eigs.vectors
end


text:
using DelimitedFiles
using DataFrames
using FluPredictibility.BioTools

"""
	`true` if position only has one amino acid
"""
function ispoly(p)
	fref = first(p.data)[2]
	fmax, aref = findmax(fref.freq)
	if fmax > 0.95
		for f in values(p.data)
			if get(f.freq, aref, 0) <= 0.95
				# println(p.i)
				# println(aref)
				# println(f.freq[aref])
				return true
			end
		end
	else
		return true
	end
	return false
end

"""
"""
function export_strains(aln::String, outfile, seqtype=:aa)
	headers = ["strain", "?", "EPI", "date", "?", "?", "?", "?","?","lab","database"]
	fp = Flu.FluPop(aln, seqtype, headers, flulineage="h3n2", segment="ha", ignore_read_errors=true)

	out = Array{String,2}(undef, length(fp.strains), 4)
	for (i,s) in enumerate(values(fp.strains))
		out[i,:] .= [s[:strain], s[:EPI], s[:lab], s[:database]]
	end

	header=["Strain" "EPI" "Lab."  "Database"];
	writedlm(outfile, vcat(header, out), "\t");
	return out
end

"""
	compute_activity_table(intrajectories, dt, max_date)

Return a `DataFrame` containing activity statistics of input trajectories.
- `dt`: width of timebins trajectories are based on.
- `max_date`: Maximum date. Length of the table will be `div(max_date, dt)+1`.
"""
function compute_activity_table(intrajectories, dt, max_date)
    imax = div(max_date, dt) + 1
    activity_table = DataFrame(:days=>collect(0:dt:max_date), :ntot=>zeros(Int64, imax),
                    :nact=>zeros(Int64, imax), :nlost=>zeros(Int64, imax), :nfixed=>zeros(Int64, imax),
                    :fact=>zeros(Float64, imax), :flost=>zeros(Float64, imax), :ffixed=>zeros(Float64, imax))

    for traj in intrajectories
        for (i,t) in enumerate([x.value for x in traj.t[1:end-1]])
            activity_table[i, :nact] += 1
            activity_table[i, :ntot] += 1
        end
        for i in (length(traj.t)):imax
            if traj.fixation == :fixed
                activity_table[i,:nfixed] += 1
                activity_table[i, :ntot] += 1
            elseif traj.fixation == :lost
                activity_table[i,:nlost] += 1
                activity_table[i, :ntot] += 1
            end
        end
    end
    activity_table[!, :fact] .= activity_table[!, :nact] ./ activity_table[!, :ntot]
    activity_table[!, :flost] .= activity_table[!, :nlost] ./ activity_table[!, :ntot]
    activity_table[!, :ffixed] .= activity_table[!, :nfixed] ./ activity_table[!, :ntot]
    return activity_table
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
# Scalar Spaces
struct ParametricScalarFunctionSpace{dim,fdim,T,comps} <: AbstractScalarFunctionSpace{dim,fdim,T,comps}
    L::Array{T,3}
    E :: Matrix{Vector{T}}
    detJf::Array{T,3}
    normals::Array{Vec{dim,T},3}
    qr_face_weigths::Vector{T}
    qr_face_points::Vector{Vec{fdim,T}}
end
struct ParametricScalarFunctionSpace{dim,T<:Real,FE<:FiniteElement,M,N1,N2,N3} <: AbstractScalarFunctionSpace{dim,T,FE,M,N1,N2,N3}
    N::Matrix{T}
    dNdξ::Matrix{Vec{dim,T}}
    detJ::Matrix{T}
    Jinv::Matrix{Tensor{2,dim,T,M}}
    M::Matrix{T}
    qr_weights::Vector{T}
    fe::FE
    mesh::PolygonalMesh{dim,N1,N2,N3,T}
end
function _scalar_fs(::Type{T}, mesh::PolygonalMesh{dim,N1,N2,N3,T}, quad_rule::QuadratureRule{dim,shape},
    felem::FiniteElement{dim,shape,order,gorder}) where {dim, T,shape<:AbstractRefShape,N1,N2,N3,order,gorder}
    n_qpoints = length(getpoints(quad_rule))
    n_cells = getncells(mesh)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    N    = fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints)
    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints)

    # Geometry interpolation
    n_geom_basefuncs = getngeombasefunctions(felem)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints)

    for (qp, ξ) in enumerate(quad_rule.points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp], N[i, qp]  = gradient(ξ -> value(felem, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp], M[i, qp] = gradient(ξ -> geom_value(felem, i, ξ), ξ, :all)
        end
    end

    detJ = fill(T(NaN), n_cells, n_qpoints)
    Jinv = fill(zero(Tensor{2,dim,T}) * T(NaN), n_cells, n_qpoints)
    coords = fill(zero(Vec{dim,T}) * T(NaN), n_nodes_per_cell(mesh))
    #Precompute detJ and invJ
    for k in 1:n_cells
        x = get_cell_coordinates!(coords,k, mesh)
        @inbounds for i in 1:n_qpoints
            fe_J = zero(Tensor{2,dim})
            for j in 1:n_geom_basefuncs
                fe_J += x[j] ⊗ dMdξ[j, i]
            end
            detJ_c = det(fe_J)
            detJ_c > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ_c)"))
            detJ[k, i] = detJ_c
            Jinv[k,i] = inv(fe_J)
        end
    end
    MM = Tensors.n_components(Tensors.get_base(eltype(Jinv)))
    ParametricScalarFunctionSpace{dim,T,eltype(felem),MM,N1,N2,N3}(N, dNdξ, detJ, Jinv,
    M, getweights(quad_rule), felem, mesh)
end

function _sface_data(::Type{T}, mesh::PolygonalMesh{dim,N1,N2,N3,T}, quad_degree::Int,
    felem::FiniteElement{dim,shape,order,gorder}, comps = 1) where {dim, T,shape<:AbstractRefShape,N1,N2,N3,order,gorder}
    fdim = dim - 1
    f_quad_rule = QuadratureRule{fdim,shape}(DefaultQuad(), quad_degree)
    n_faceqpoints = length(getpoints(f_quad_rule))
    q_ref_facepoints = getpoints(f_quad_rule)
    q_ref_faceweights = getweights(f_quad_rule)
    face_quad_rule = create_face_quad_rule(f_quad_rule)

    n_cells = getncells(mesh)
    n_faces = n_faces_per_cell(mesh)

    # Face interpolation
    n_func_basefuncs = getnbasefunctions(felem)
    E    = fill(zeros(T,n_faces) * T(NaN), n_func_basefuncs, n_faceqpoints)

    # Geometry face Interpolation
    n_geom_basefuncs = getngeombasefunctions(felem)
    L =    fill(zero(T)          * T(NaN), n_geom_basefuncs, n_faceqpoints, n_faces)
    dLdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_faceqpoints, n_faces)

    for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
        for i in 1:n_geom_basefuncs
            dLdξ[i, qp, face], L[i, qp, face] = gradient(ξ -> geom_value(felem, i, ξ), ξ, :all)
        end
    end

    for qp in 1:n_faceqpoints
        for i in 1:n_func_basefuncs
            E_f = zeros(T, n_faces)
            for j in 1:n_faces
                #Evaluate shape function on q_point map to edge j
                E_f[j] = value(felem, i, face_quad_rule[j].points[qp])
            end
            E[i, qp] = E_f
        end
    end

    detJf = fill(T(NaN), n_cells, n_faces, n_faceqpoints)
    coords = fill(zero(Vec{dim,T}) * T(NaN), n_nodes_per_cell(mesh))
    normals = fill(zero(Vec{dim,T}) * T(NaN), n_cells, n_faces, n_faceqpoints)
    #Precompute detJ and invJ
    for k in 1:n_cells
        x = get_cell_coordinates!(coords,k, mesh)
        @inbounds for i in 1:n_faceqpoints
            for l in 1:n_faces
                fef_J = zero(Tensor{2,dim})
                for j in 1:n_geom_basefuncs
                    fef_J += x[j] ⊗ dLdξ[j, i, l]
                end
                weight_norm = weighted_normal(fef_J, l, shape, Val{dim})
                detJ_f = norm(weight_norm)
                normals[k,l,i] = weight_norm / detJ_f
                detJ_f > 0.0 || throw(ArgumentError("det(Jf) is not positive: det(Jf) = $(detJ_f)"))
                detJf[k,l,i] = detJ_f
            end
        end
    end
    ParametricScalarFunctionSpace{dim,fdim,T,comps}(L, E, detJf, normals, q_ref_faceweights, q_ref_facepoints)
end

#Data
@inline getdetJdV(fs::ParametricScalarFunctionSpace, cell::Int, q_point::Int) = fs.detJ[cell,q_point]*fs.qr_weights[q_point]
@inline shape_gradient(fs::ParametricScalarFunctionSpace, q_point::Int, base_func::Int, cell::Int) = fs.dNdξ[base_func, q_point] ⋅ fs.Jinv[cell,q_point]
@inline shape_divergence(fs::ParametricScalarFunctionSpace, q_point::Int, base_func::Int, cell::Int) = sum(fs.dNdξ[base_func, q_point] ⋅ fs.Jinv[cell,q_point])
@inline getdetJdS(fs::ParametricScalarFunctionSpace, cell::Int, face::Int, q_point::Int) = fs.detJf[cell,face, q_point]*fs.qr_face_weigths[q_point]
"""
    getnormal(fs::ScalarFunctionSpace, cell::Int, face::Int, qp::Int)
Return the normal at the quadrature point `qp` for the face `face` at
cell `cell` of the `ScalarFunctionSpace` object.
"""
@inline get_normal(fs::ParametricScalarFunctionSpace, cell::Int, face::Int, qp::Int) = fs.normals[cell, face, qp]

#Vector Space
# VectorFunctionSpace
struct ParametricVectorFunctionSpace{dim,T<:Real,FE<:FiniteElement,M,N1,N2,N3} <: AbstractVectorFunctionSpace{dim,T,FE,M,N1,N2,N3}
    n_dof::Int
    ssp::ParametricScalarFunctionSpace{dim,T,FE,M,N1,N2,N3}
end

#Constructor
function VectorFunctionSpace(mesh::PolygonalMesh, felem::FiniteElement{dim,shape,order,gorder};
    face_data = true, quad_degree = order+1) where {dim, shape, order, gorder}
    quad_rule = QuadratureRule{dim,shape}(DefaultQuad(), quad_degree)
    fs = _scalar_fs(Float64, mesh, quad_rule, felem)
    fd = face_data ? _sface_data(Float64, mesh, quad_degree, felem, dim) : nothing
    n_func_basefuncs = getnbasefunctions(felem)
    dof = n_func_basefuncs*dim
    if face_data
        return ParametricVectorFunctionSpace(dof,fs), fd
    else
        return ParametricVectorFunctionSpace(dof,fs)
    end
end

#Data
function shape_gradient(fs::ParametricVectorFunctionSpace{dim,T}, q_point::Int, base_func::Int, cell::Int)  where {dim,T}
    @assert 1 <= base_func <= fs.n_dof "invalid base function index: $base_func"
    dN_comp = zeros(T, dim, dim)
    n = size(fs.ssp.N,1)
    dN_comp[div(base_func,n+1)+1, :] = fs.ssp.dNdξ[mod1(base_func,n), q_point]
    return Tensor{2,dim,T}((dN_comp...,)) ⋅ fs.ssp.Jinv[cell,q_point]
end