text:
type NSBEntropy
	N::Int64
	K::Int64
	nx::Array{Int64,1} #unique counts
	kx::Array{Int64,1} #number of each unique count
	nc::Array{Int64,1} #conidience count
	kc::Array{Int64,1} #number of each unique coincidence count
	K1::Int64 # number of bins non-zero counts
	K2::Int64 #number of bins with coincidence
	B_interp::Array{Float64} #interpolation grid
	S_nsb::Float64 #nsb entropy
	dS_nsb::Float64 #standard devation
	S_ml::Float64 #maximum likelihood estimator
end

function NSBEntropy()
	K = 0
	nx = Array(Int64,0)
	kx = Array(Int64,0)
	nc = Array(Int64,0)
	kc = Array(Int64,0)
	B_interp = Array(Float64,0)
	S_nsb = 0.0
	dS_nsb = 0.0
	S_ml = 0.0
end

function NSBEntropy(n::Array{Int64,1},K::Int64)
	k = StatsBase.counts(n, 1:maximum(n))
	nx = [1:maximum(n)][k.!=0]
	kx = k[k.>0]
	idx = nx .> 1
	nc = nx[idx]
	kc = kx[idx]
	N = sum(nx.*kx)
	K1 = sum(kx)
	K2 = sum(kc)
	return NSBEntropy(N,K,nx,kx,nc,kc,K1, K2, Float64[], 0.0, 0.0, 0.0)
end

function NSBEntropy(nx::Array{Int64,1},kx::Array{Int64,1}, K::Int64)
	idx = nx .> 1
	nc = nx[idx]
	kc = kx[idx]
	N = sum(nx.*kx)
	K1 = sum(kx)
	K2 = sum(kc)
	return NSBEntropy(N,K,nx,kx,nc,kc,K1, K2, Float64[], 0.0, 0.0, 0.0)
end

function nsb_entropy(n::Array{Int64,1},K::Int64;precision::Real=1e-5)
	S_nsb = NSBEntropy(n,K)
	find_nsb_entropy(S_nsb,precision)
	return S_nsb.S_nsb, S_nsb.dS_nsb
end

function nsb_entropy(X::Dict{Int64,Int64})
	n = collect(values(X))
	K = maximum(keys(X))
	return nsb_entropy(n,K)
end

function nsb_entropy(X::Dict{Int64,Dict{Int64,Int64}})
    E = 0
    N = 0
    for (k,v) in X
        ee,n = nsb_entropy(v)
        E += ee*n
        N += n
    end
    return E/N,N
end

function nsb_entropy(X::Dict{Int64, Dict{Int64, Dict{Int64,Int64}}})
    E = 0
    N = 0
    for (k,v) in X
        ee,n = nsb_entropy(v)
        E += ee*n
        N += n
    end
    return E/N,N
end

function get_coincidence_count(S::NSBEntropy)
	if isempty(S.nc)
		idx = S.nx .> 1
		S.nc = S.nx[idx]
		S.kc = S.kx[dx]
	end
	return S.kc,S.nc
end

psi(x) = polygamma(0,x)
psi_asymp(x) = psi(x) - log(x) #FIXME is perhaps not as accurate as we want

function find_nsb_entropy(n::Array{Int64,1}, K::Integer, precision::Real;verbose::Int64=0)
	S_nsb = NSBEntropy(n,K)
	k = StatsBase.counts(n, 1:maximum(n))
	nx = [1:maximum(n)][k.!=0]
	kx = k[k.>0]
	return find_nsb_entropy(S_nsb,precision;verbose=verbose)
end

function find_nsb_entropy(S_nsb::NSBEntropy, precision::Real;verbose::Int64=0)
	integs=  [integrand_1,integrand_S,integrand_S2]
	msgs = ["normalization", "S","S^2"]
	K = S_nsb.K
	nsb_mlog_quad = zeros(1)
	edges =[precision, log(K) - precision]
	xi_lim = zeros(2)
	B_cl, xi_cl, dS_cl,errcode = max_evidence(S_nsb, precision;verbose=verbose)
	#println("size(xi_cl) = $(size(xi_cl))")
	#println("errocode = $errcode")
	if errcode > 0
		verbose > 0 && warn("FIND_ENTROPY: Switching integration over the full range of xi due to prior errors")
		#println("B_cl = $B_cl")
		if B_cl > 0
			nsb_mlog_quad = mlog_evidence(B_cl, S_nsb)
			S_cl = meanS(B_cl, S_nsb)
		else
			nsb_mlog = 0.0
			S_cl = NaN
		end
		xi_lim[:] = edges
		delta = 0.0
		dS_cl = NaN
		xi_cl = Inf
	else
		S_cl = meanS(B_cl, S_nsb)
		verbose > 0 && warn("FIND_NSB_ENTROPY: Expect S to be near $S_cl and σ near $dS_cl")
		#println("B_cl = $B_cl")
		nsb_mlog_quad = mlog_evidence(B_cl, S_nsb) #value at the saddle
		delta = erfinv(1-precision)*sqrt(2)
		verbose > 0 && println("FIND_NSB_ENTROPY: Integrating around peak")
	end

	nsb_val_quad = zeros(length(integs),1)
	ec = zeros(nsb_val_quad)
	nfun = zeros(nsb_val_quad)
	err = zeros(nsb_val_quad)
	for i in 1:length(integs)
		if delta > 0
			#println("size(xi_cl) = $(size(xi_cl))")
			#println("size(nsb_mlog_quad) = $(size(nsb_mlog_quad))")
			cent = integs[i](xi_cl, S_nsb,nsb_mlog_quad[1])
			good_enough = 0.1*cent*exp(-(delta^2)/2.0)
			wnd = [-1.0,-1.0]
			worst = [1.0,2.0]
			best = [1.0,2.0]
			limval = good_enough.*ones(2)
			finished = false
			while !finished
				if any(xi_lim[worst].==edges)
					if any(wnd[best].>10)
						wnd[best] *= 1.2
					else
						wnd[best] += 0.5
					end
				else
					if any(wnd[worst].>10)
						wnd[worst] *= 1.2
					else
						wnd[worst] += 0.5
					end
				end
				xi_lim[1] = max(edges[1],xi_cl - (delta+wnd[1])*dS_cl)[1]
				xi_lim[2] = min(edges[2],xi_cl + (delta+wnd[2])*dS_cl)[1]
				#println("xi_lim = $xi_lim")
				#println("limval = $limval")
				#println("dS_cl = $dS_cl")
				#println("delta = $delta")
				#println("nsb_mlog_quad = $nsb_mlog_quad")

				for j in 1:length(worst)
					limval[worst[j]] = integs[i](xi_lim[worst[j]],S_nsb, nsb_mlog_quad[1])[1]
				end
				for j in 1:length(best)
					limval[best[j]] = integs[i](xi_lim[best[j]],S_nsb,nsb_mlog_quad[1])[1]
				end
				tmp,worst = findmax(limval)
				tmp,best = findmin(limval)
				#println("best = $best")
				if any(limval[worst].<good_enough)
					finished = true
				end
				if any(limval[best] .< good_enough) && any(xi_lim[worst].==edges)
					finished = true
				end
				if xi_lim == edges
					finished = true
				end
			end
		end

		verbose >0 && println("FIND_NSB_ENTROPY: Doing $(msgs[i]) integral. Limits: $(xi_lim[1]) < xi < $(xi_lim[2])")

		_func(x) = integs[i](x, S_nsb,nsb_mlog_quad[1])
		nsb_val_quad[i,:],err[i] = quadgk(_func, xi_lim[1], xi_lim[2];reltol=precision)
		ec[i] = 0
		if ec[i] > 0
			verbose > 0 && warn("FIND_NSB_ENTROPY: Problem in $(msgs[i]) integral")
		end
		if err[i] > precision
			verbose > 0 && warn("FIND_NSB_ENTROPY: Precision of $precision required by only $(err[i]) achieved")
		end

	end

	S_nb  = nsb_val_quad[2]/nsb_val_quad[1]
	#println("S_nb = $S_nb")
	S_nsb.S_nsb = S_nb
	S_nsb.dS_nsb = sqrt(abs(nsb_val_quad[3]/nsb_val_quad[1] - S_nb^2))
	return S_nsb,err
end

function meanS{T<:Real}(B::T, kx::Array{Int64,1}, nx::Array{Int64,1}, K::Int64)
	return meanS([B],kx,nx,K)
end

function meanS{T<:Real}(B::Array{T}, kx::Array{Int64,1}, nx::Array{Int64,1}, K::Int64)
	@assert size(kx) == size(nx)
	if any(nx.<=0)
		ArgumentError("Bins with zero count encountered")
	end
	N = sum(kx.*nx)
	Nf = N + 0.0
	if N <= 1
		verbose > 0 && warn("Too few data samples: N = $N")
	end
	K1 = sum(kx)
	_B = B[:]

	ovrNB = 1./(Nf+_B)
	osB = ones(size(_B))
	osn = ones(size(nx))
	f = zeros(size(_B))

	f = psi(Nf+_B+1.0) - (ovrNB*osn').*(osB*nx' + (_B./K)*osn').*psi(osB*nx' + (_B./K)*osn'+1.0)*kx[:] -
			_B.*ovrNB*(1.0-K1/K).*psi(1+_B/K)
	f[!isfinite(_B)] = log(K)
	f = reshape(f, size(B))
	f
end

function meanS{T<:Real}(B::Array{T,1}, S_nsb::NSBEntropy)
	f = zeros(B)
	for i in 1:length(f)
		f[i] = meanS(B[i],S_nsb)
	end
	f
end

function meanS(B::Real, S_nsb::NSBEntropy)
	N = S_nsb.N
	Nf = N + 0.0
	K = S_nsb.K
	if !isfinite(B)
		return log(K)
	end
	K1 = S_nsb.K1
	nx = S_nsb.nx
	kx = S_nsb.kx

	ovrNB = 1./(Nf+B)

	#f = psi(Nf+_B+1.0) - (ovrNB*osn').*(osB*nx' + (_B./K)*osn').*psi(osB*nx' + (_B./K)*osn'+1.0)*kx[:] -
	#		_B.*ovrNB*(1.0-K1/K).*psi(1+_B/K)
	f = psi(N + B+1.0) - (ovrNB.*(nx + (B/K)).*psi(nx + B/K + 1.0))'*kx - B*ovrNB*(1.0 - K1/K).*psi(1+B/K)
	f[1]
end

function xi_KB{T<:Real}(K::Int64, B::T)
	return psi(B+1) - psi(1+B/K)
end

function xi_KB{T<:Real}(K::Int64,B::Array{T})
	xi = zeros(B)
	for i in 1:length(xi)
		xi[i] = psi(B[i]+1) - psi(1+B[i]/K)
	end
	#xi = psi(B+1) - psi(1+B./K)
	xi
end

function dxi_KB(K::Int64, B::Real)
	dxi = polygamma(1,B+1) - 1/K.*polygamma(1,1+B./K)
	return dxi
end

function dxi_KB{T<:Real}(K::Int64, B::Array{T})
	dxi = polygamma(1,B+1) - 1/K.*polygamma(1,1+B./K)
	return dxi
end

function max_evidence(n::Array{Int64,1}, K::Int64, precision::Real)
	kx = StatsBase.counts(n, 1:maximum(n))
	nx = [1:maximum(n)][kx.!=0]
	kx = kx[kx.!=0]
	return max_evidence(kx,nx,K,precision)
end

function get_coefficients(N::Int64)
	N2 = N*N
	N3 = N2*N
	N4 = N3*N
	N5 = N4*N
	N6 = N5*N
	N7 = N6*N
	N8 = N7*N
	N9 = N8*N
	N10 = N9*N
	N11 = N10*N
	N12 = N11*N

	Nf = float(N)
	Nm = Nf-1.0
	overN = 1/N
	Nm2 = Nm*Nm
	Nm3 = Nm2*Nm
	Nm4 = Nm3*Nm
	Nm5 = Nm4*Nm
	Nm6 = Nm5*Nm
	Nm7 = Nm6*Nm
	Nm8 = Nm7*Nm
	Nm9 = Nm8*Nm
	Nm10 = Nm9*Nm

	b = [(Nm)/(2*N), #b1
        (-2+1/N)/3, #b2
        (Nf^2 - N - 2)/(9*(Nf^2 - N)), #b3
	    (2*(2 - 3*N - 3*Nf^2 + 2*Nf^3))/(135*Nm^2*N), #b4
	    (4.0*(22.0 + 13.0*N - 12*N2 - 2.0*N3 + N4))/(405.0*Nm3*N), #b5
	    (4.0*(-40.0 + 58.0*N + 65.0*N2 - 40.0*N3 - 5.0*N4 +
	      2.0*N5))/(1701.0*Nm4*N), #b6
	    (4*(-9496 - 6912*N + 5772*Nf^2 + 2251*Nf^3 -
	      1053*Nf^4 - 87*Nf^5 + 29*Nf^6))/(42525*Nm^5*N), #b7
	    (16.0*(764.0 - 1030.0*N - 1434.0*Nf^2 + 757.0*Nf^3 + 295.0*Nf^4
	       - 111.0*Nf^5 - 7.0*Nf^6 + 2.0*Nf^7))/(18225.0*Nm^6*N), #b8
	    (16*(167000 + 142516*N - 108124*Nf^2 - 66284*Nf^3
	       + 26921*Nf^4 + 7384*Nf^5 - 2326*Nf^6 - 116*Nf^7
	       + 29*Nf^8))/(382725*Nm^7*N), #b9
	    (16.0*(-17886224.0 + 22513608.0*N + 37376676.0*Nf^2
	       - 17041380.0*Nf^3 - 11384883.0*Nf^4 + 3698262.0*Nf^5
	       + 846930.0*Nf^6 - 229464.0*Nf^7 - 9387.0*Nf^8 + 2086.0*Nf^9))/
			(37889775.0*Nm^8*N), #b10
	    (16.0*(-4166651072.0 - 3997913072.0*Nf + 2783482560.0*Nf^2 +
	       2290151964.0*Nf^3 - 803439834.0*Nf^4 - 395614251.0*Nf^5 +
	       108055443.0*Nf^6 + 20215218.0*Nf^7 - 4805712.0*Nf^8 - 165395.0*Nf^9
	       + 33079.0*Nf^10))/(795685275.0*Nm^9*N), #b11
	    (32.0*(52543486208.0 - 62328059360.0*Nf - 118489458160.0*Nf^2 +
	       47185442088.0*Nf^3 + 44875379190.0*Nf^4 - 12359832987.0*Nf^5 -
	       5400540075.0*Nf^6 + 1272974916.0*Nf^7 + 200644800.0*Nf^8 -
	       42495955.0*Nf^9 - 1255067.0*Nf^10 +
	       228194.0*Nf^11))/(14105329875.0*Nm^10*N)] #b12
   return b
end

function max_evidence(kx::Array{Int64,1}, nx::Array{Int64,1}, K::Int64, precision::Real;verbose::Int64=0)
	max_counter = 200
	errcode = 0

	if any(nx .<= 0)
		ArgumentError("Bins with zero count encountered")
	end
	N = sum(kx.*nx)
	if N <= 1
		ArgumentError("Too few data samples: N = $(N)")
	end

	ng1 = find(nx.>1)
	K1 = sum(kx)
	K2 = sum(kx[ng1])

	ep = (N-K1)/N

	if K1 >= N
		Bcl = Inf
		xicl = Inf
		errcode = 1
		dxi = NaN
		verbose > 0 && warn("MAX_EVIDENCE: No coincidence")
	elseif K1 <= 1
		Bcl = 0.0
		xicl = 0
		errcode = 2
		dxi = NaN
		verbose >0 && warn("MAX EVIDENCE: All data coincidence")
	else
		order = 10
		b = get_coefficients(N)

		B0ep = N*sum(ep.^[-1:order].*b)
		if B0ep < 0
			B0ep = precision
			verbose > 0 && warn("MAX_EVIDENCE: Series expansion for B_0 did not converge")
		end

		B0 = B0ep
		dB0 = 99999999
		counter = 0
		while ~((abs(dB0) < abs(B0*precision)) | (counter .> max_counter))
			counter += 1
			F = K1/B0 + psi(B0) - psi(B0+N)
			dF = -K1/B0^2 + polygamma(1,B0) - polygamma(1,B0+N)
			dB0 = -F/dF
			B0 += dB0
		end

		if counter > max_counter
			errcode = 3
			verbose > 0 && warn("MAX_EVIDENCE: Newton-Rhapson search for B_0 did not converge after $counter iterations")
		end

		Bcl = B0
		order_K = 4
		B = zeros(1,order_K)

		EG = -psi(1)
		pg1B0 = polygamma(1,B0)
		pg1NB0 = polygamma(1,N+B0)
		denum = K1/B0^2 - pg1B0 + pg1NB0
		pg2B0 = polygamma(2,B0)
		pg2NB0 = polygamma(2,N+B0)
		pg21 = polygamma(2,1)
		pg3B0 = polygamma(3,B0)
		pg3NB0 = polygamma(4,N+B0)
		pg4B0 = polygamma(4,B0)
		pg4NB0 = polygamma(3,N+B0)

		nxng1 = nx[ng1]
		kxng1 = kx[ng1]
		skxng1 = sum(kx[ng1])
		f0 = sum(kxng1.*polygamma(0,nxng1))
		d1f0 = sum(kxng1.*polygamma(1,nxng1))
		d2f0 = sum(kxng1.*polygamma(2,nxng1))
		d3f0 = sum(kxng1.*polygamma(3,nxng1))

		#println("f0 = $f0, d1f0 = $d1f0, d2f0 = $d2f0, d3f0 = $d3f0")

		B[1] = (B0^2*(EG*K2 + f0)) / (B0^2*denum)
		B[2] = (K2*pi^2*B0 - (6*K1*B[1]^2)/B0^3 - 3*B[1]^2*pg2B0 +
				3*B[1]^2*pg2NB0 - 6*B0*d1f0)/(-6*denum);
		B[3] = (K2*pi^2*B[1] + (6*K1*B[1]^3)/B0^4 -(12*K1*B[1]*B[2])/B0^3 +
				3*K2*B0^2*pg21 - 6*B[1]*B[2]*pg2B0 + 6*B[1]*B[2]*pg2NB0 -
				B[1]^3*pg3B0 + B[1]^3*pg3NB0 - 6*B[1]*d1f0 - 3*B0^2*d2f0)/ (-6*denum)
		B[4] = -(-(K2*pi^4*B0^3)/90 + (K1*B[1]^4)/B0^5 - (K2*pi^2*B[2])/6 -
				(3*K1*B[1]^2*B[2])/B0^4 + (K1*B[2]^2)/B0^3 +
				(2*K1*B[1]*B[3])/B0^3 - K2*B0*B[1]*pg21 + ((B[2]^2 +
				2*B[1]*B[3])*pg2B0)/2
				- ((B[2]^2 + 2*B[1]*B[3])*pg2NB0)/2 +
				(B[1]^2*B[2]*pg3B0)/2 - (B[1]^2*B[2]*pg3NB0)/2 +
				(B[1]^4*pg4B0)/ 24 - (B[1]^4*pg4NB0)/24 +  B[2]*d1f0 +
				B0*B[1]*d2f0 + (B0^3*d3f0)/6)/(-denum)

		Φ = sum(B'.*((K+0.0).^(-[1:order_K])))
		Bcl += Φ
		#println("B = $B")
		#println("Bcl = $Bcl")
		#println("Φ = $Φ")
		#ccorrect until here

		counter = 0
		dBcl = 999999999999
		F = 0.0
		dF = 0.0
		while !((abs(dBcl) < abs(Bcl*precision)) || (counter > max_counter))
			counter += 1
			F = 1/K*sum(kx[ng1].*psi(nxng1 + Bcl/K)) - K2/K*psi(1+Bcl/K) +
				K1/Bcl + psi(Bcl) - psi(Bcl+N)

			dF =  1/(K^2)*sum(kx[ng1].*polygamma(1,nxng1 + Bcl/K)) - K2/(K^2)*polygamma(1,1+Bcl/K) - K1/Bcl^2 + polygamma(1,Bcl) - polygamma(1, Bcl+N)
			#println("counter = $counter, F = $F, dF = $dF")

			dBcl = -F/dF
			if dBcl < 0.0
				verbose > 0 && warn("MAX_EVIDENCE: negative (unstable) change dBcl = $dBcl")
				dBcl = 0.0
				errcode = 4
			end
			Bcl = Bcl + dBcl
		end
		#println("Bcl = $Bcl")

		if counter > max_counter
			errcode = 3
			verbose > 0 && warn("Newton-Raphson search for Bcl did not converge after $counter iterations")
		end

		if errcode == 3 && counter < max_counter
			verbose > 0 && println("MAX_EVIDENCE. Recovered from previous error in Bcl determination")
		end

		dBcl = 1/K^2*sum(kxng1.*polygamma(1,nxng1 + Bcl/K)) -
				K2/K^2*polygamma(1,1+Bcl/K) - K1/Bcl^2 + polygamma(1,Bcl) -
				polygamma(1,Bcl + N)

		xicl = xi_KB(K,Bcl)
		verbose > 0 && println("xicl = $xicl")
		k = 0
		verbose > 0 && println("dBcl = $dBcl")
		#println(K)
		#println(dxi_KB(K,Bcl))
		dxi = 1./sqrt(-dBcl./dxi_KB(K,Bcl).^2)
		#FIXME: These numbers don't quite match with octave yet
	end
	return Bcl, xicl, dxi, errcode
end

function max_evidence(S_est::NSBEntropy, precision::Real;verbose::Int64=0)
	max_counter = 200
	errcode = 0

	N = S_est.N
	K = S_est.K
	K1 = S_est.K1
	K2 = S_est.K2
	kc = S_est.kc
	nc = S_est.nc

	ep = (N-K1)/N

	if K1 >= N
		Bcl = Inf
		xicl = Inf
		errcode = 1
		dxi = NaN
		verbose > 0 && warn("MAX_EVIDENCE: No coincidence")
	elseif K1 <= 1
		Bcl = 0.0
		xicl = 0
		errcode = 2
		dxi = NaN
		verbose > 0 && warn("MAX EVIDENCE: All data coincidence")
	else
		order = 10
		b = get_coefficients(N)

		B0ep = N*sum(ep.^[-1:order].*b)
		if B0ep < 0
			B0ep = precision
			verbose > 0 && warn("MAX_EVIDENCE: Series expansion for B_0 did not converge")
		end

		B0 = B0ep
		dB0 = 99999999
		counter = 0
		while ~((abs(dB0) < abs(B0*precision)) | (counter .> max_counter))
			counter += 1
			F = K1/B0 + psi(B0) - psi(B0+N)
			dF = -K1/B0^2 + polygamma(1,B0) - polygamma(1,B0+N)
			dB0 = -F/dF
			B0 += dB0
		end

		if counter > max_counter
			errcode = 3
			verbose >0 && warn("MAX_EVIDENCE: Newton-Rhapson search for B_0 did not converge after $counter iterations")
		end

		Bcl = B0
		order_K = 4
		B = zeros(1,order_K)

		EG = -psi(1)
		pg1B0 = polygamma(1,B0)
		pg1NB0 = polygamma(1,N+B0)
		denum = K1/B0^2 - pg1B0 + pg1NB0
		pg2B0 = polygamma(2,B0)
		pg2NB0 = polygamma(2,N+B0)
		pg21 = polygamma(2,1)
		pg3B0 = polygamma(3,B0)
		pg3NB0 = polygamma(3,N+B0)
		pg4B0 = polygamma(4,B0)
		pg4NB0 = polygamma(4,N+B0)

		nxng1 = S_est.nc
		kxng1 = S_est.kc
		skxng1 = S_est.K2
		f0 = sum(kxng1.*polygamma(0,nxng1))
		d1f0 = sum(kxng1.*polygamma(1,nxng1))
		d2f0 = sum(kxng1.*polygamma(2,nxng1))
		d3f0 = sum(kxng1.*polygamma(3,nxng1))

		#println("f0 = $f0, d1f0 = $d1f0, d2f0 = $d2f0, d3f0 = $d3f0")

		B[1] = (B0^2*(EG*K2 + f0)) / (B0^2*denum)
		B[2] = (K2*pi^2*B0 - (6*K1*B[1]^2)/B0^3 - 3*B[1]^2*pg2B0 +
				3*B[1]^2*pg2NB0 - 6*B0*d1f0)/(-6*denum);
		B[3] = (K2*pi^2*B[1] + (6*K1*B[1]^3)/B0^4 -(12*K1*B[1]*B[2])/B0^3 +
				3*K2*B0^2*pg21 - 6*B[1]*B[2]*pg2B0 + 6*B[1]*B[2]*pg2NB0 -
				B[1]^3*pg3B0 + B[1]^3*pg3NB0 - 6*B[1]*d1f0 - 3*B0^2*d2f0)/ (-6*denum)
		B[4] = -(-(K2*pi^4*B0^3)/90 + (K1*B[1]^4)/B0^5 - (K2*pi^2*B[2])/6 -
				(3*K1*B[1]^2*B[2])/B0^4 + (K1*B[2]^2)/B0^3 +
				(2*K1*B[1]*B[3])/B0^3 - K2*B0*B[1]*pg21 + ((B[2]^2 +
				2*B[1]*B[3])*pg2B0)/2
				- ((B[2]^2 + 2*B[1]*B[3])*pg2NB0)/2 +
				(B[1]^2*B[2]*pg3B0)/2 - (B[1]^2*B[2]*pg3NB0)/2 +
				(B[1]^4*pg4B0)/ 24 - (B[1]^4*pg4NB0)/24 +  B[2]*d1f0 +
				B0*B[1]*d2f0 + (B0^3*d3f0)/6)/(-denum)

		Φ = sum(B'.*((K+0.0).^(-[1:order_K])))
		Bcl += Φ
		#println("B = $B")
		#println("B0 = $B0") #correct
		#println("Bcl = $Bcl")
		#println("Φ = $Φ")
		#ccorrect until here

		counter = 0
		dBcl = 999999999999
		F = 0.0
		dF = 0.0
		while !((abs(dBcl) < abs(Bcl*precision)) || (counter > max_counter))
			counter += 1
			F = 1/K*sum(kc.*psi(nxng1 + Bcl/K)) - K2/K*psi(1+Bcl/K) +
				K1/Bcl + psi(Bcl) - psi(Bcl+N)

			dF =  1/(K^2)*sum(kc.*polygamma(1,nxng1 + Bcl/K)) - K2/(K^2)*polygamma(1,1+Bcl/K) - K1/Bcl^2 + polygamma(1,Bcl) - polygamma(1, Bcl+N)
			#println("counter = $counter, F = $F, dF = $dF")

			dBcl = -F/dF
			if dBcl < 0.0
				verbose > 0 && warn("MAX_EVIDENCE: negative (unstable) change dBcl = $dBcl")
				dBcl = 0.0
				errcode = 4
			end
			Bcl = Bcl + dBcl
		end
		#println("Bcl = $Bcl")

		if counter > max_counter
			errcode = 3
			verbose > 0 && warn("Newton-Raphson search for Bcl did not converge after $counter iterations")
		end

		if errcode == 3 && counter < max_counter
			verbose > 0 && println("MAX_EVIDENCE. Recovered from previous error in Bcl determination")
			errcode = 0
		end

		dBcl = 1/K^2*sum(kxng1.*polygamma(1,nxng1 + Bcl/K)) -
				K2/K^2*polygamma(1,1+Bcl/K) - K1/Bcl^2 + polygamma(1,Bcl) -
				polygamma(1,Bcl + N)

		xicl = xi_KB(K,Bcl)
		verbose > 0 && println("xicl = $xicl")
		verbose > 0 && println("dBcl = $dBcl")
		#println(K)
		#println(dxi_KB(K,Bcl))
		#println("Bcl = $Bcl")
		if errcode > 0
			dxi = 0.0
		else
			dxi = 1./sqrt(-dBcl./dxi_KB(K,Bcl).^2)
		end
		#FIXME: These numbers don't quite match with octave yet
		#println("dxi = $dxi")

	end
	return Bcl, xicl, dxi, errcode
end

gammaln(x) = log(gamma(x))

mlog_evidence(B::Real, kx::Array{Int64,1}, nx::Array{Int64,1}, K::Int64) = mlog_evidence([B],kx,nx,K)

function mlog_evidence{T<:Real}(B::Array{T,1}, kx::Array{Int64,1}, nx::Array{Int64,1}, K::Int64)
	if any(nx.<=0)
		ArgumentError("Bins with zero count encountered")
	end
	N = sum(kx.*nx)
	if N <= 1
		ArgumentError("Too few data samples: N = $N")
	end
	Nf = N + 0.0
	_B = B[:]
	f = zeros(_B)
	B0 = find(_B.==0)
	if length(B0) > 0
		Bn0 = find(_B.!=0)
		f[Bn0] = mlog_evidence(B[Bn0], kx, nx, K)
		f[B0] = Inf
	else
		ng1 = find(nx.>1)
		nxng1 = nx[ng1]
		K1 = sum(kx)
		K2 = sum(kx[ng1])
		if length(ng1) > 0 #coincidences
			f1 = K2*lgamma(1+_B/K)
			f2 = -lgamma(_B*ones(1,length(ng1)))/K + ones(length(B),1)*nx[ng1]'
			f3 = f2*kx[ng1]
			f[:] = f1 .+ f3[:]
			#for i in 1:length(f)
			#	f[i] = K2*lgamma(1+_B[i]/K)
			#	f2 = -lgamma(_B[i])/K
			#	f1 = 0.0
		#		f3 = 0.0
		#		for j1 in 1:length(ng1)
		#			f1 += nxng1[j1]
		#			f3 += nxng1[j1]^2
		#		end
		#		f[i] += f2*f1 + f3
		#	end
		end
	end
	#println(f)
	large = (_B .> max(100,100*N))
	if any(large) #we have large values
		nterms = max(ceil(abs(-15 - log10(N))./log10(N./_B[large]))) + 1
		#println("nterms = $nterms")
		ifac = 1./cumprod([2:(nterms+1)])
		for i in nterms:-1:1
			f[large]  += polygamma(i,_B[large]).*ifrac[i].*Nf^(i+1)
		end
	end
	#println(f)

	other = !large
	df = -K1*log(_B[other]) + lgamma(_B[other]+N) - lgamma(_B[other])
	#println(df)
	f[other] += -K1*log(_B[other]) + lgamma(_B[other]+N) - lgamma(_B[other])
end

function mlog_evidence{T<:Real}(B::Array{T,1}, S_nsb::NSBEntropy)
	f = zeros(size(B))
	for i in 1:length(f)
		f[i] = mlog_evidence(B[i],S_nsb)
	end
	return f
end

function mlog_evidence(B::Real, S_nsb::NSBEntropy)
	N = S_nsb.N
	Nf = N + 0.0
	f = 0.0
	if B == 0.0
		f = Inf
	else
		nxng1 = S_nsb.nc
		kxng1 = S_nsb.kc
		K = S_nsb.K
		K1 = S_nsb.K1
		K2 = S_nsb.K2
		if length(nxng1) > 0 #coincidences
			f1 = K2*lgamma(1.0+B/K)
			f2 = -lgamma(B/K + nxng1)
			f3 = f2'*kxng1
			#println("f1 = $f1")
			#println("f2 = $f2")
			#println("f3 = $f3")
			f = (f1 .+ f3[:])[1]
		end
	end
	#println("f = $f")
	#println("MLOG_EVIDENCE: f = $f")
	if B > max(100,100*N)
		nterms = int(ceil(abs((-15.0 - log10(N))./log10(N./B))))
		#println("nterms = $nterms")
		ifac = 1./cumprod([2:(nterms+1)])
		#println("ifac = $ifac")
		for i in nterms:-1:1
			f += polygamma(i,B).*ifac[i].*Nf^(i+1)
			#println("f = $f")
		end
		f += (N-K1)*log(B) + psi_asymp(B)*N
		#println("f = $f")
	else
	#println(f)
		#df = -K1*log(B) + lgamma(B+N) - lgamma(B)
		f += -K1*log(B) + lgamma(B+N) - lgamma(B)
	end
	#println("f = $f")
	return f
end

integrand_1{T<:Real}(xi::T,nsb_kx_quad::Array{Int64,1}, nsb_nx_quad::Array{Int64,1}, nsb_K_quad::Int64,nsb_mlog_quad::Real) = integrand_1([xi],nsb_kx_quad, nsb_nx_quad, nsb_K_quad,nsb_mlog_quad)

function integrand_1{T<:Real}(xi::Array{T},nsb_kx_quad::Array{Int64,1}, nsb_nx_quad::Array{Int64,1}, nsb_K_quad::Int64,nsb_mlog_quad::Real)
	B = B_xiK(xi, nsb_K_quad)
	_mle = mlog_evidence(B,nsb_kx_quad, nsb_nx_quad, nsb_K_quad)
	_pxi = prior_xi(xi, nsb_K_quad)
	f = exp( -_mle + nsb_mlog_quad).*_pxi
	#println("f = $f")
	return f
end

function integrand_1(xi::Real,S_nsb::NSBEntropy,nsb_mlog_quad::Real)
	B = B_xiK(xi,S_nsb)
	_mle = mlog_evidence(B,S_nsb)
	_pxi = prior_xi(xi, S_nsb.K)
	f = exp( -_mle + nsb_mlog_quad).*_pxi
	#println("f = $f")
	#println("_mle = $_mle")
	#println("nsb_mlog_quad = $nsb_mlog_quad")
	return f
end

integrand_S{T<:Real}(xi::T,nsb_kx_quad::Array{Int64,1}, nsb_nx_quad::Array{Int64,1}, nsb_K_quad::Int64,nsb_mlog_quad::Real) = integrand_S([xi],nsb_kx_quad, nsb_nx_quad, nsb_K_quad,nsb_mlog_quad)

function integrand_S{T<:Real}(xi::Array{T},nsb_kx_quad::Array{Int64,1}, nsb_nx_quad::Array{Int64,1}, nsb_K_quad::Int64,nsb_mlog_quad::Real)
	B = B_xiK(xi, nsb_K_quad)
	_mle = mlog_evidence(B,nsb_kx_quad, nsb_nx_quad, nsb_K_quad)
	_ms = meanS(B,nsb_kx_quad, nsb_nx_quad, nsb_K_quad)
	_pxi = prior_xi(xi, nsb_K_quad)
	f = exp( -_mle + nsb_mlog_quad).*_pxi.*_ms
	return f
end

function integrand_S(xi::Real,S_nsb::NSBEntropy,nsb_mlog_quad::Real)
	B = B_xiK(xi, S_nsb)
	_mle = mlog_evidence(B,S_nsb)
	_ms = meanS(B,S_nsb)
	_pxi = prior_xi(xi, S_nsb.K)
	f = exp( -_mle + nsb_mlog_quad).*_pxi.*_ms
	return f
end

function integrand_S2(xi::Real,S_nsb::NSBEntropy,nsb_mlog_quad::Real)
	B = B_xiK(xi, S_nsb)
	_mle = mlog_evidence(B,S_nsb)
	_ms2 = meanS2(B,S_nsb)
	_pxi = prior_xi(xi, S_nsb.K)
	f = exp( -_mle + nsb_mlog_quad).*_pxi.*_ms2
	return f
end

function meanS2(B::Real, S_nsb::NSBEntropy)

	N = S_nsb.N
	K = S_nsb.K
	K1 = S_nsb.K1
	nx = S_nsb.nx
	kx = S_nsb.kx
	b = B/K
	nxb = nx + b
	pnxb1 = psi(nxb+1.0)
	p0NB2 = psi(N + B + 2.0)
	p1NB2 = polygamma(1,N+B+2.0)
	f = sum(nxb.*((pnxb1-p0NB2).*kx)*(nxb.*(pnxb1-p0NB2).*kx)' - (nxb.*kx)*(nxb.*kx)'*p1NB2)
	f += sum(2*B*(1-K1/K)*nxb.*((pnxb1-p0NB2).*(psi(b+1)-p0NB2) - p1NB2).*kx)
	f += (1-K1/K)*(1-(K1+1)/K)*B*B*((psi(b+1) - p0NB2)^2-p1NB2)
	f += -sum((nxb.*(pnxb1 - p0NB2)).^2.*kx) + sum(nxb.*nxb.*kx)*p1NB2

	f += sum((nxb.*(nxb+1).*((psi(nxb+2) - p0NB2).^2 + polygamma(1,nxb+2) - p1NB2)).*kx)

	f += B*(1-K1/K)*(1+b)*((psi(2+b) - p0NB2).^2 + polygamma(1,b+2) - p1NB2)

	#normalizing
	f = f/((N+B)*(N+B+1))
	if !isfinite(f)
		f = log(K)^2
	end
	return f
end

function prior_xi(xi::Real, K::Int64)
	return 1.0
end

function prior_xi{T<:Real}(xi::Array{T}, K::Int64)
	f = ones(size(xi))
	return f
end

function B_xiK(K::Int64)
	#create the interpolation
	K_interp_in = K
	step = 1.0e-2
	b1 = 10.0
	b2 = 100.0*K
	b = zeros(int(b1/step+1 + ceil(log(b2)/step)))
	b[1:b1/step+1] = 0:step:b1
	b[b1/step+3:end] = b1*exp([step:step:log(b2)])
	#b = [0:step:b1,b1*exp([step:step:log(b2)])]
	Bxi_interp_in = zeros(2,length(b))
	Bxi_interp_in[1,:] = xi_KB(K,b)
	Bxi_interp_in[2,:] = b
	dxi = float([1.0e-3:log(K)*1.0e-3:log(K)-1.0e-3])
	Bxi_interp_in_new = zeros(2,length(dxi))
	Bxi_interp_in_new[1,:] = dxi
	Bxi_interp_in_new[2,:]  = B_xiK(Bxi_interp_in, dxi,K)
	return Bxi_interp_in_new
end

function B_xiK{T<:Real}(xi::Array{T}, S_nsb::NSBEntropy)
	if isempty(S_nsb.B_interp)
		S_nsb.B_interp = B_xiK(S_nsb.K)
	end
	return B_xiK(S_nsb.B_interp, xi, S_nsb.K)
end

function B_xiK(xi::Real, S_nsb::NSBEntropy)
	if isempty(S_nsb.B_interp)
		S_nsb.B_interp = B_xiK(S_nsb.K)
	end
	return B_xiK(S_nsb.B_interp, xi, S_nsb.K)
end

function B_xiK{T<:Real}(xi::Array{T}, K::Int64)
	if any(xi .> log(K))
		ArgumentError("Too large xi -- bigger than log(K)")
	end
	if any(xi .< 0 )
		ArgumentError("Too small xi -- smaller than 0")
	end
	#create the interpolation
	K_interp_in = K
	step = 1.0e-2
	b1 = 10.0
	b2 = 100.0*K
	b = zeros(int(b1/step+1 + ceil(log(b2)/step)))
	b[1:b1/step+1] = 0:step:b1
	b[b1/step+3:end] = b1*exp([step:step:log(b2)])
	#b = [0:step:b1,b1*exp([step:step:log(b2)])]
	Bxi_interp_in = zeros(2,length(b))
	Bxi_interp_in[1,:] = xi_KB(K,b)
	Bxi_interp_in[2,:] = b
	dxi = float([1.0e-3:log(K)*1.0e-3:log(K)-1.0e-3])
	Bxi_interp_in_new = zeros(2,length(dxi))
	Bxi_interp_in_new[1,:] = dxi
	Bxi_interp_in_new[2,:]  = B_xiK(Bxi_interp_in, dxi,K)

	return B_xiK(Bxi_interp_in_new,xi,K)
end

function B_xiK{T<:Real}(Bxi_interp_in::Array{T,2},xi::Array{T}, K::Int64)
	max_counter = 200
	if maximum(xi) > log(K)
		ArgumentError("Too large xi -- bigger than log(K)")
	end
	if minimum(xi) < 0
		ArgumentError("Too small xi -- smaller than 0")
	end
	_xi = xi[:]
	B = Bxi_interp_in[2,lookup(Bxi_interp_in[1,2:size(Bxi_interp_in,2)-1],_xi)+1]
	B = B[:]
	dB = zeros(B)
	counter = 0
    F   =  9999999.0*ones(_xi)
    dF   =  9999999.0*ones(_xi)
	#println(minimum(_xi))
	qq = abs(_xi*1e-13 + 1e-13)
	qv = _xi*1.0e-13
	while !(all(abs(F) .< qq) || counter > max_counter)
		counter += 1
		#F[:] = xi_KB(K,B) -_xi
		#println(maximum(F))
		#dF[:] = dxi_KB(K,B)
		#dB[:] = -F./dF
		for i in 1:length(F)
			F[i] = xi_KB(K,B[i]) -_xi[i]
			dF[i] = dxi_KB(K,B[i])
			dB[i] = abs(F[i]) >= qv[i] ? -F[i]/dF[i] : 0.0
		end
		#dB[abs(F).<qv] = 0.0
		B[:] += dB
	end
	#println(size(_xi))
	if counter >= max_counter
		error("Newton-Raphson root polishing did not converge after $counter iterations. Problems are likely.")
	end

	repl_large = ((B.>100*K)&(_xi.!=log(K)))
	B[repl_large] = (K-1.0)./(2*(log(K) - _xi[repl_large]))
	B[_xi.==log(K)] = Inf

	repl_small = _xi .<  (pi*pi*1e-6/6.0)
	B[repl_small] = (1.0 + 1/(K-1.0))*_xi[repl_small]*(6/pi^2)

	B = reshape(B,size(xi))
	return B
end

function B_xiK{T<:Real}(Bxi_interp_in::Array{T,2},xi::Real, K::Int64)
	max_counter = 200
	if xi > log(K)
		ArgumentError("Too large xi -- bigger than log(K)")
	end
	if xi < 0
		ArgumentError("Too small xi -- smaller than 0")
	end
	B = Bxi_interp_in[2,lookup(Bxi_interp_in[1,2:size(Bxi_interp_in,2)-1],xi)+1]
	dB = 0.0
	counter = 0
    F   =  9999999.0
    dF   =  9999999.0
	#println(minimum(_xi))
	qq = abs(xi*1e-13 + 1e-13)
	qv = xi*1.0e-13
	while !((abs(F) < qq) || counter > max_counter)
		counter += 1
		#F[:] = xi_KB(K,B) -_xi
		#println(maximum(F))
		#dF[:] = dxi_KB(K,B)
		#dB[:] = -F./dF
		F = xi_KB(K,B) -xi
		dF = dxi_KB(K,B)
		dB = abs(F) >= qv ? -F/dF : 0.0
		#dB[abs(F).<qv] = 0.0
		B += dB
	end
	#println(size(_xi))
	if counter >= max_counter
		error("Newton-Raphson root polishing did not converge after $counter iterations. Problems are likely.")
	end

	if (B > 100*K)&&(xi != log(K))
		B = (K-1.0)./(2*(log(K) - xi))
	elseif xi == log(K)
		B = Inf
	elseif xi <  (pi*pi*1e-6/6.0)
		B = (1.0 + 1/(K-1.0))*xi*(6/pi^2)
	end
	return B
end

function lookup{T<:Real}(table::Array{T,2},xi::Array{T})
	return lookup(table[:],xi)
end

function lookup{T<:Real}(table::Array{T,2},xi::Real)
	return lookup(table[:],xi)
end

function lookup{T<:Real}(table::Array{T,1},xi::Array{T})
	idx = zeros(Int64,size(xi))
	for i in 1:length(xi)
		idx[i] = last(searchsorted(table,xi[i]))
	end
	return idx
end

function lookup{T<:Real}(table::Array{T,1},xi::Real)
	idx = last(searchsorted(table,xi))
	return idx
end


text:
## tape

const MaybeFunction = Union{Function, Nothing}

mutable struct Tape
    ops::Vector{<:AbstractOp}   # linearized execution graph
    resultid::Int               # id of result variable
    derivs::Dict{Int,Int}       # derivs[var.id] == grad_var.id
    sfields::Dict{Int, Dict}    # mapping of argid -> Dict(struct field paths -> var id)
    compiled::MaybeFunction     # compiled tape or nothing
    device::AbstractDevice      # function to use for moving intermediate results to device
    Tape(device::AbstractDevice) = new(AbstractOp[], -1, Dict(), Dict(), nothing, device)
    Tape() = Tape(CPU())
end

function Base.show(io::IO, tape::Tape)
    rev_derivs = Dict((j, i) for (i, j) in tape.derivs)
    # rev_sfields = Dict((var.id, argid) for )
    println(io, "Tape")
    for (i, op) in enumerate(tape.ops)
        hint_deriv = haskey(rev_derivs, i) ? "deriv for %$(rev_derivs[i])" : ""
        # hint_struct = haskey(tape.sfields, ) ?
        hint = isempty(hint_deriv) ? "" : "\t # " * hint_deriv
        println(io, "  $op$hint")
    end
end

Base.getindex(tape::Tape, i::Integer...) = getindex(tape.ops, i...)
Base.getindex(tape::Tape, i::String...) =
    getindex(tape.ops, [parse(Int, s[2:end]) for s in i]...)
Base.getindex(tape::Tape, i::Symbol...) = getindex(tape, map(string, i)...)
Base.lastindex(tape::Tape) = lastindex(tape.ops)
Base.length(tape::Tape) = length(tape.ops)
Base.iterate(tape::Tape) = iterate(tape.ops)
Base.iterate(tape::Tape, s) = iterate(tape.ops, s)

"""
Record an operation onto a tape, assign new ID to op's var.
"""
function _record!(tape::Tape, op::AbstractOp)
    push!(tape.ops, op)
    op.var.id = length(tape)
    nothing
end

ext: jl


text:
using TropicalTensors
using Test
using Yao

@testset "test red reg" begin
    sg = Spinglass(ChimeraLattice(3,3), ones(Float32, 100), ones(Float32, 100))
    reg = TropicalTensors.red_reg(Tropical{Float32}, sg, 3, 1, 1; usecuda=false)
    @test reg isa ArrayReg{1,Tropical{Float32},Matrix{Tropical{Float32}}}
end

@testset "sg bonds" begin
    lt = ChimeraLattice(3, 2)
    bonds = sgbonds(lt)
    bonds[1] == (1, 25)
    bonds[2] == (1, 26)
    bonds[5] == (5, 1)
    bonds[37] == (5, 13)
    @test length(bonds) == 16*6 + 3*4 + 16
end

@testset "test Chimera" begin
    lt = ChimeraLattice(3, 3)
    sg = Spinglass(lt, ones(Float32, 12*4 + 9*16), zeros(Float32, 9*8))
    res = solve(sg; usecuda=false)
    @test res.n == 12*4 + 9*16
end

@testset "counting tropical" begin
    lt = ChimeraLattice(3, 3)
    sg = Spinglass(lt, ones(Float32, 12*4 + 9*16), zeros(Float32, 72))
    res = solve(CountingTropical{Float32}, sg; usecuda=false)
    @test res.n == 12*4 + 9*16
    @test res.c == 2
end


text:
#=
Example of using a module. This example uses the module contained in the
folder named MyModule.
=#

# Load the module
include("MyModule/MyModule.jl")

# Use the module and its exported function
using .MyModule

# Use the exported function
myfuncA()

# Use a function from the module
MyModule.myfuncB()

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
using LaTeXStrings
using Plots
gr()

function plot_res(result, powergrid,disturbed_node)
    ω_indices = findall(n -> isa(n, SwingEqLVS), powergrid.nodes)
    append!(ω_indices,findall(n -> isa(n, VSIMinimal), powergrid.nodes))
    #append!(ω_indices,findall(n -> isa(n, PVInverterWithFrequencyControl), powergrid.nodes))
    #append!(ω_indices,findall(n -> isa(n, WindTurbineGenType4_RotorControl), powergrid.nodes))
    append!(ω_indices,findall(n -> isa(n, FourthOrderEq), powergrid.nodes))
    append!(ω_indices,findall(n -> isa(n, ThirdOrderEq),powergrid.nodes))

    #ω_colors = reshape(Plots.get_color_palette(cgrad(:inferno), 25)[ω_indices], (1,length(ω_indices)))
    ω_labels = reshape([latexstring(string(raw"\omega", "_{$i}")) for i=ω_indices], (1, length(ω_indices)))
    p_labels = reshape([latexstring(string(raw"p", "_{$i}")) for i=1:length(powergrid.nodes)], (1, length(powergrid.nodes)))
    v_labels = reshape([latexstring(string(raw"v", "_{$i}")) for i=1:length(powergrid.nodes)], (1, length(powergrid.nodes)))

    pl_p = plot(result, :, :p, legend = (0.8, 0.95), ylabel=L"p [p.u.]", label=p_labels)
    pl_v = plot(result, :, :v, legend = (0.8, 0.7), ylabel=L"v [p.u.]", label=v_labels)
    pl_ω = plot(result, ω_indices, :ω, legend = (0.8, 0.7), ylabel=L"\omega \left[rad/s\right]", label=ω_labels)
    pl = plot(
        pl_p, pl_ω,pl_v;
        layout=(3,1),
        size = (500, 500),
        lw=3,
        xguide=L"t[s]"
    )
    display(pl)
end

max_stars_repo_name: luap-pik/PowerDynIC

max_issues_repo_name: luap-pik/PowerDynIC

max_forks_repo_name: luap-pik/PowerDynIC


text:
using PyPlot
using Jagot.plotting
plot_style("ggplot")
using PyPlotRecipes
using PyCall
PathEffects = pyimport("matplotlib.patheffects")
using Statistics
using Random
using Colors

using LinearAlgebra
using ArnoldiMethod

using CompactBases
import CompactBases: applied
using IntervalSets

function mean_color(color::String)
    c = parse(Colorant, color)
    mean([c.r,c.g,c.b])
end

lerp(a,b,t) = (1-t)*a + t*b

mean_position(x, ϕ) = ϕ'*Diagonal(x)*ϕ/(ϕ'ϕ)

function savedocfig(name,dir="figures")
    fig = gcf()
    filename = joinpath(@__DIR__, "src", dir, "$(name).svg")
    savefig(filename,
            transparent=true,
            facecolor=fig.get_facecolor())
    close(fig)
    if isfile(filename)
        println("Saved $(name) to $(filename)")
    else
        @warn "Saving $(name) to $(filename) failed"
    end
end

function logo()
    t = ArbitraryKnotSet(3, [0.0, 1, 1, 3, 4, 6], 1, 3)
    r = range(first(t), stop=last(t), length=301)
    R = BSpline(t)
    χ = R[r, :]

    cfigure("logo",figsize=(7,3)) do
        plot(r, χ, length(r) < 30 ? ".-" : "-", linewidth=4)
        margins(0.1, 0.1)
        axis("off")
    end
    savefig("docs/src/assets/logo.svg", transparent=true)
end

function simple_example()
    a,b = 0.0,1.0 # Extents
    N = 13 # Number of nodes
    k = 5 # Order of FE-DVR/B-splines
    x = range(a, stop=b, length=1000)

    Δx = (b-a)/(N+1) # Grid spacing
    # Standard, uniform finite-differences
    fd = FiniteDifferences(N, Δx)
    # Staggered, uniform finite-differences
    sfd_uni = StaggeredFiniteDifferences(N, Δx)
    # Staggered, non-uniform finite-differences
    sfd_nonuni = StaggeredFiniteDifferences(0.01, 0.5, 0.1, b)

    # Finite-element boundaries
    tf = range(a, stop=b, length=N+2)
    # We can vary the polynomial order in each element
    forder = vcat(7, fill(4,length(tf)-2))
    # By indexing the second dimension, we can implement Dirichlet0
    # boundary conditions.
    fem = FEDVR(tf, forder)[:,2:end-1]

    tb = ExpKnotSet(k, -2.0, log10(b), N+1)
    splines = BSpline(tb)[:,2:end-1]

    bases = [fd, sfd_uni, sfd_nonuni, fem, splines]
    m = length(bases)

    cfigure("simple example", figsize=(6,9)) do
        for (i,basis) in enumerate(bases)
            csubplot(m,1,i,nox=i<m) do
                plot(x, basis[x,:])
            end
        end
    end
    savedocfig("simple_example")
end

function restricted_bases()
    B = FEDVR(range(0, stop=10.0, length=5), 7)
    B̃ = B[:, 3:17]

    x = range(0, stop=10.0, length=1000)
    χ = B[x, :]
    χ̃ = B̃[x, :]

    cfigure("restricted bases") do
        csubplot(211,nox=true) do
            plot(x, χ)
            ylabel("Original basis")
        end
        csubplot(212) do
            plot(x, χ̃)
            ylabel("Restricted basis")
        end
    end
    savedocfig("restricted_bases")
end

function densities()
    f = x -> sin(2π*x)
    g = x -> x*exp(-x)
    h = x -> f(x)*g(x)

    rmin = 0.01
    rmax = 10.0
    ρ = 0.05
    N = ceil(Int, rmax/ρ)
    k = 7
    Nn = 71

    ρmin=rmin
    ρmax=0.5
    α=0.01

    coeff_centers = Vector{Vector{Float64}}()
    density_coeffs = Vector{Vector{Float64}}()
    errors = Vector{Vector{Float64}}()
    ϕ = Vector{Vector{Float64}}()

    rr = 10.0 .^ range(-2, stop=log10(rmax), length=3000)

    for R in [FiniteDifferences(N, ρ),
              StaggeredFiniteDifferences(ρmin, ρmax, α, rmax),
              FEDVR(range(0, stop=rmax, length=Nn), k),
              BSpline(LinearKnotSet(k, 0, rmax, Nn))]
        r = axes(R,1)

        cf = R \ f.(r)
        cg = R \ g.(r)
        ch = R \ h.(r)

        ρ = Density(applied(*,R,cf), applied(*,R,cg))

        push!(coeff_centers, centers(R))
        push!(density_coeffs, ρ.ρ)
        push!(errors, ρ.ρ-ch)
        push!(ϕ, R[rr,:]*ρ.ρ)
    end

    cfigure("mutual density", figsize=(7,10)) do
        csubplot(411,nox=true) do
            semilogx(rr, f.(rr), linewidth=1.0)
            semilogx(rr, g.(rr), linewidth=1.0)
            semilogx(rr, h.(rr), linewidth=2.0)
            legend([L"f(x)", L"g(x)", L"h(x)"])
        end
        csubplot(412, nox=true) do
            for ϕ in ϕ
                loglog(rr, abs.(ϕ-h.(rr)))
            end
            yl = ylim()
            ylim(max(yl[1],1e-15), yl[2])
            legend(["Finite-differences",
                    "Log–linear finite-differences",
                    "FE-DVR",
                    "B-splines"])
            ylabel("Reconstruction error")
        end
        csubplot(413, nox=true) do
            for (r,c) in zip(coeff_centers, density_coeffs)
                semilogx(r, c, ".-", linewidth=1.0)
            end
            legend(["Finite-differences",
                    "Log–linear finite-differences",
                    "FE-DVR",
                    "B-splines"])
            ylabel("Expansion coefficients")
        end
        csubplot(414) do
            for (r,e) in zip(coeff_centers, errors)
                loglog(r, abs.(e), ".-", linewidth=1.0)
            end
            legend(["Finite-differences",
                    "Log–linear finite-differences",
                    "FE-DVR",
                    "B-splines"])
            ylabel("Coefficient error")
            xlabel(L"r")
        end
    end

    savedocfig("mutual_densities")
end

function diagonal_operators()
    Z = 1.0
    rmax = 20*(2^1.05)/Z

    ρ = 0.1/Z
    N = ceil(Int, rmax/ρ)
    ufd = StaggeredFiniteDifferences(N, ρ)

    ρmin=0.1/Z
    ρmax=0.6
    α=0.002
    nufd = StaggeredFiniteDifferences(ρmin, ρmax, α, rmax)

    fedvr = FEDVR(range(0, stop=rmax, length=20),
                  vcat(10,fill(7,18)))[:,2:end-1]

    bsplines = BSpline(ExpKnotSet(5, -2.0, log10(rmax), 100))[:,2:end-1]

    R = bsplines

    r = axes(R,1)
    rr = clamp.(10.0 .^ range(min(-2,log10(r[1])), stop=log10(rmax), length=3000),
                0, rightendpoint(R))
    χ = R[rr,:]

    r = axes(R,1)
    V = r -> (0.314269680527354*r^2 + 0.20951312035157*r)*exp(-3*r/2);
    Vc = R \ V.(r);
    rc = R \ identity.(r);

    L = DiagonalOperator(applied(*, R, Vc));

    Vr = L*rc

    Vr_exact = rr .* V.(rr)

    cfigure("diagonal operators") do
        csubplot(211,nox=true) do
            semilogx(rr, χ*Vr, label=L"V(r)r")
            semilogx(rr, Vr_exact, "--", label=L"$V(r)r$ exact")
            semilogx(rr, V.(rr), label=L"V(r)")
            yl = ylim()
            semilogx(rr, rr, label=L"r")
            ylim(yl)
            legend()
        end
        csubplot(212) do
            for (R,label) in [
                (ufd, "Uniform finite-differences"),
                (nufd, "Non-uniform finite-differences"),
                (fedvr, "FE-DVR"),
                (bsplines, "B-splines")
            ]
                r = axes(R,1)
                Vc = R \ V.(r)
                rc = R \ identity.(r)

                L = DiagonalOperator(applied(*, R, Vc))

                Vr = L*rc
                χ = R[rr,:]
                loglog(rr, abs.(χ*Vr - Vr_exact), label=label)
            end
            xlabel(L"r")
            ylabel("Error")
            legend()
        end
    end

    savedocfig("diagonal_operators")
end

macro echo(expr)
    println(expr)
    :(@time $expr)
end

@info "Documentation plots"
fig_dir = joinpath(@__DIR__, "src", "figures")
mkpath(fig_dir)
@assert isdir(fig_dir)
@echo logo()
@echo simple_example()
@echo restricted_bases()
include("bspline_plots.jl")
include("fd_plots.jl")
@echo densities()
@echo diagonal_operators()


text:
@testset "Spherical" begin

    @testset "simplest system" begin
        PolyImpl.@polyvar x
        H = StraightLineHomotopy([(x - 4.3im) * (x + (2.1 - 4im))], [(x - 2.0) * (x - (2.5+ 4.0im))])

        res = solve(H, [[-2.1 + 4.0im], [4.3im]], Spherical(), endgame_start=0.0)

        @test res[1].returncode == :Success
        @test res[2].returncode == :Success

        @test real(res[1].solution) ≈ [2.0]
        @test real(res[2].solution) ≈ [2.5]
        @test imag(res[2].solution) ≈ [4.0]
    end

end


text:
type SceneRecord
    scenes::Vector{Scene}
    timestep::Float64
    nscenes::Int # number of active scenes

    function SceneRecord(max_n_scenes::Int, timestep::Float64, max_n_vehicles::Int=500)
        scenes = Array{Scene}(max_n_scenes)
        for i in 1 : length(scenes)
            scenes[i] = Scene(max_n_vehicles)
        end
        new(scenes, timestep, 0)
    end
end

Base.show(io::IO, rec::SceneRecord) = print(io, "SceneRecord(nscenes=", rec.nscenes, ")")

Base.length(rec::SceneRecord) = rec.nscenes
function Base.deepcopy(rec::SceneRecord)
    retval = SceneRecord(length(rec.scenes), rec.timestep, length(rec.scenes[1].vehicles))
    for i in 1 : rec.nscenes
        copy!(retval.scenes[i], rec.scenes[i])
    end
    retval
end

record_length(rec::SceneRecord) = length(rec.scenes)
pastframe_inbounds(rec::SceneRecord, pastframe::Int) = 1 ≤ 1-pastframe ≤ rec.nscenes
get_scene(rec::SceneRecord, pastframe::Int) = rec.scenes[1 - pastframe]
get_elapsed_time(rec::SceneRecord, pastframe::Int) = (1-pastframe)*rec.timestep
function get_elapsed_time(
    rec::SceneRecord,
    pastframe_farthest_back::Int,
    pastframe_most_recent::Int,
    )

    (pastframe_most_recent - pastframe_farthest_back)*rec.timestep
end

function Base.empty!(rec::SceneRecord)
    rec.nscenes = 0
    rec
end

Base.getindex(rec::SceneRecord, vehicle_index::Int, pastframe::Int) = get_scene(rec, pastframe)[vehicle_index]
Base.getindex(rec::SceneRecord, vehicle_index::Int) = Base.getindex(rec, vehicle_index, 0)

function Base.setindex!(rec::SceneRecord, veh::Vehicle, vehicle_index::Int, pastframe::Int)

    veh_internal = rec[vehicle_index, pastframe]
    copy!(veh_internal, veh)

    rec
end
Base.setindex!(rec::SceneRecord, veh::Vehicle, vehicle_index::Int) = Base.setindex!(rec, veh, vehicle_index, 0)

function get_vehiclestate(rec::SceneRecord, id::Int, pastframe::Int=0)
    scene = get_scene(rec, pastframe)
    get_vehiclestate(scene, id)
end
function get_vehicle!(veh::Vehicle, rec::SceneRecord, id::Int, pastframe::Int=0)
    scene = get_scene(rec, pastframe)
    get_vehicle!(veh, scene, id)
end
get_vehicle(rec::SceneRecord, id::Int, pastframe::Int=0) = get_vehicle!(rec, id, pastframe)

function iscarinframe(rec::SceneRecord, id::Int, pastframe::Int=0)
    scene = get_scene(rec, pastframe)
    iscarinframe(scene, id)
end
function get_index_of_first_vehicle_with_id(rec::SceneRecord, id::Int, pastframe::Int=0)
    scene = get_scene(rec, pastframe)
    get_index_of_first_vehicle_with_id(scene, id)
end

function push_back_records!(rec::SceneRecord)
    for i in min(rec.nscenes+1,length(rec.scenes)) : -1 : 2
        copy!(rec.scenes[i], rec.scenes[i-1])
    end
    rec
end
function Base.insert!(rec::SceneRecord, scene::Scene, pastframe::Int=0)
    scene_internal = get_scene(rec, pastframe)
    copy!(scene_internal, scene)
    rec
end
function update!(rec::SceneRecord, scene::Scene)
    push_back_records!(rec)
    insert!(rec, scene, 0)
    rec.nscenes = min(rec.nscenes+1, record_length(rec))
    rec
end

function Base.get!(scene::Scene, rec::SceneRecord, pastframe::Int=0)
    scene_internal = get_scene(rec, pastframe)
    copy!(scene, scene_internal)
    scene
end


text:
using CalibrationErrors, Distances

const x1, y1 = 1, 3
const x2, y2 = [1, -1], [2, 1]

@testset "Squared exponential kernel" begin
    @test_throws ArgumentError SquaredExponentialKernel(-1)

    @test SquaredExponentialKernel() === SquaredExponentialKernel(1)
    @test SquaredExponentialKernel() === SquaredExponentialKernel(Euclidean())
    @test SquaredExponentialKernel() === SquaredExponentialKernel(1, Euclidean())

    for γ in (1, 2)
        kernel = SquaredExponentialKernel(γ, TotalVariation())
        @test kernel(x1, y1) ≈ exp(- γ)
        @test kernel(x2, y2) ≈ exp(- 2.25 * γ)

        kernel2 = SquaredExponentialKernel(γ, Cityblock())
        @test kernel2(x1, y1) ≈ exp(- 4 * γ)
        @test kernel2(x2, y2) ≈ exp(- 9 * γ)

        kernel3 = SquaredExponentialKernel(γ, Euclidean())
        @test kernel3(x1, y1) ≈ exp(- 4 * γ)
        @test kernel3(x2, y2) ≈ exp(- 5 * γ)
    end
end

@testset "Exponential kernel" begin
    @test_throws ArgumentError ExponentialKernel(-1)

    @test ExponentialKernel() === ExponentialKernel(1)
    @test ExponentialKernel() === ExponentialKernel(Euclidean())
    @test ExponentialKernel() === ExponentialKernel(1, Euclidean())

    for γ in (1, 2)
        kernel = ExponentialKernel(γ, TotalVariation())
        @test kernel(x1, y1) ≈ exp(- γ)
        @test kernel(x2, y2) ≈ exp(- 1.5 * γ)

        kernel2 = ExponentialKernel(γ, Cityblock())
        @test kernel2(x1, y1) ≈ exp(- 2 * γ)
        @test kernel2(x2, y2) ≈ exp(- 3 * γ)

        kernel3 = ExponentialKernel(γ, Euclidean())
        @test kernel3(x1, y1) ≈ exp(- 2 * γ)
        @test kernel3(x2, y2) ≈ exp(- sqrt(5) * γ)
    end
end


text:
# Pretty slow solution. Both part 1 and 2 take around 20 seconds
using DataStructures
function solve(game::Vector{Vector{Char}})
    game = [line[2:end-1] for line in game[2:end-1]]
    sy, sx = length(game), length(game[1])
    depth = sy - 1
    nswe = [(-1, 0), (1, 0), (0, 1), (0, -1)]
    adjacent(y, x) = [(y+ya, x+xa) for (ya, xa)=nswe if (y+ya)∈1:sy && (x+xa)∈1:sx]
    is_legal(yx) = occursin(game[yx[1]][yx[2]], ".ABCD")
    adjacent_legal = Dict()
    legal_coords = []
    target_lookup = Dict(a => ("."^i * a^(depth-i) for i=1:depth) for a in "ABCD")
    allowed_lookup = [Set(), Set((1, xa) for xa in (1, 2, 4, 6, 8, 10, 11))]

    for y=1:sy, x=1:sx
        if is_legal((y, x))
            adjacent_legal[(y, x)] = filter(is_legal, adjacent(y, x))
            (game[y][x] != '.' || game[y+1][x] == '#') && push!(legal_coords, (y, x))
        end
    end
    state_queue = PriorityQueue(game => 0)
    states_visited = Set{UInt64}()
    while !isempty(state_queue)
        state, cost = dequeue_pair!(state_queue)
        push!(states_visited, hash(state))
        if all(state[ys][xs] == ('A'+(xs-3)÷2) for (ys, xs) in legal_coords[end-(4*depth-1):end])
            return cost
        end
        for (y, x) in legal_coords
            if (a = state[y][x]) != '.'
                target = target_lookup[a]
                allowed = allowed_lookup[(y != 1) + 1]
                target_x = 3 + 2(a - 'A')
                slot = join(state[2+d][target_x] for d=0:depth-1)
                slot_pos = slot ∈ target ? (1 + sum(collect(slot) .== '.'), target_x) : (0, 0)
                isempty(allowed) && slot_pos == (0, 0) && continue
                visited = Set{Tuple{Int64, Int64}}()
                queue = Queue{Tuple{Int64, Int64, Int64}}()
                enqueue!(queue, (cost, y, x))
                while !isempty(queue)
                    new_cost, yc, xc = dequeue!(queue)
                    if slot_pos == (yc, xc) || (yc, xc) ∈ allowed
                        new_state = deepcopy(state)
                        new_state[yc][xc] = a
                        new_state[y][x] = '.'
                        if hash(new_state) ∉ states_visited
                            state_queue[new_state] = min(new_cost, get(state_queue, new_state, typemax(Int)))
                        end
                    end
                    push!(visited, (yc, xc))
                    for (ya, xa) in adjacent_legal[(yc, xc)]
                        if state[ya][xa] == '.' && (ya, xa) ∉ visited
                            adj_cost = new_cost + (10^(a - 'A'))
                            push!(visited, (ya, xa))
                            enqueue!(queue, (adj_cost, ya, xa))
                        end
                    end
                end
            end
        end
    end
end

function main()
    game = readlines() .|> collect
    println(solve(game))

    insert!(game, length(game)-1, collect("  #D#C#B#A#  "))
    insert!(game, length(game)-1, collect("  #D#B#A#C#  "))
    println(solve(game))
end
main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
#This file contains an Example Lesson plan for Julia tutor.

prompt_1 = ""*
"Welcome to your first day of class. Today we'll be doing a classic excercise. We will be learning "*
"how to print text to our console. Along the way we'll cover some ways to handle strings in Julia."*
"\n"*
"Strings are defined by double quotes \" \" in julia. For example, we could define a variable, "*
"called console_output to contain the text \"hello world\" like the following,\n"*
"console_output = \"hello world\"\n"*
"\n"*
"Okay now it's your turn,\n"
request_1 = "Define a string titled \'console_output\' to contain the text \"hello world\""
evaluator_1 = ExpressionEvaluator(:console_output, "hello world")
hint_1 = """Please ensure your solution is cases sensitive.
IE: "HeLlO WoRlD" is not the same as "hello world". """

#These four items compose a single lesson instance.
lesson_1 = Lesson( prompt_1, request_1, evaluator_1, hint_1 )

prompt_2 = ""*
"Congratulations you have instantiated your first variable in Julia! How might we display this"*
"text to console? There are two functions that let you do this: print() and println(). Print " *
"will display the String or variable as it is, but println will add a new line/carriage return"
" to the end of it.\n"*
"\n"*
"Now that our console_output variable is stored in memory, we can print it as follows:\n"*
"print(console_output)\n"
request_2 = "Print your \'console_output\' variable to the console using the `print()` function"
evaluator_2 = CommandEvaluator("hello world")
#Notice: we did not provide any hints for this excercise. Not nice, but either they get it,
#or they don't.
lesson_2 = Lesson( prompt_2, request_2, evaluator_2 )

#To create the lesson_plan (which is REQUIRED) we make an array of Lesson's.
lesson_plan = [ lesson_1, lesson_2 ]


text:
using Distributions
using GaussianMixtureTest
using Test

mu_true = [-2.0858,-1.4879]
wi_true = [0.0828,0.9172]
sigmas_true = [0.6735,0.2931]

m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
x = rand(m, 1000);

wi, mu, sigmas, ml = gmm(x, 2)
mhat = MixtureModel(map((u, v) -> Normal(u, v), mu, sigmas), wi)
@test loglikelihood(mhat, x) >= loglikelihood(m, x)

asymptoticdistribution(x, wi_true, mu_true, sigmas_true);
T, p = kstest(x, 2)

@test T >= 0
# write your own tests here


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

counts = Cint[mod(i,2) + 1 for i in 0:(size-1)]
check = collect(Iterators.flatten([fill(r, counts[r+1]) for r = 0:size-1]))

for T in Base.uniontypes(MPI.MPIDatatype)
    A = ArrayType{T}(fill(T(rank), counts[rank+1]))
    B = MPI.Allgatherv(A, counts, comm)
    @test B isa ArrayType{T}
    @test B == ArrayType{T}(check)

    # Test passing the output buffer
    B = ArrayType{T}(undef, sum(counts))
    MPI.Allgatherv!(A, B, counts, comm)
    @test B == ArrayType{T}(check)

    # Test assertion when output size is too small
    B = ArrayType{T}(undef, sum(counts)-1)
    @test_throws AssertionError MPI.Allgatherv!(A, B, counts, comm)

    # Test explicit MPI_IN_PLACE
    B = ArrayType(fill(T(rank), sum(counts)))
    MPI.Allgatherv!(MPI.IN_PLACE, B, counts, comm)
    @test B == ArrayType{T}(check)
end

MPI.Finalize()
@test MPI.Finalized()


text:
import POSIXClock

typeof(Pkg.installed("Gadfly"))==Void &&
  error("Please install Gadfly package to run this example.")
import Gadfly; pl=Gadfly

const interval = 200000 # 200us
function run()
  nItr = div(60000000000,interval)
  # pre-allocate samples array
  differ = Vector{Int64}(nItr+1)
  # pre-allocate timespec
  t = POSIXClock.timespec(0,0)
  actual_t = POSIXClock.timespec(0,0)

  # Lock future memory allocations, disable GC
  ccall(:mlockall, Cint, (Cint,), 2) # CAUTION: will crash Julia on future memory allocations
  gc_enable(false)

  n = 0
  POSIXClock.gettime!(t,POSIXClock.CLOCK_MONOTONIC)
  @time while n <= nItr
    POSIXClock.nanosleep!(t,interval)
    POSIXClock.gettime!(actual_t,POSIXClock.CLOCK_MONOTONIC)
    n+=1
    differ[n]=actual_t-t
  end

  # Unlock memory allocations, enable GC
  ccall(:munlockall, Cint, ())
  gc_enable(true)

  return differ
end

const d=run()
histplot=pl.plot(
  x=d/1000, pl.Geom.histogram(),
  pl.Guide.xlabel("Latency (us)"),
  pl.Guide.ylabel("Histogram counts"))


text:
module Outliers

using StatsBase

export outliers

include("outliers.jl")

end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: redshift_data
using AWS.Compat
using AWS.UUIDs
"""
    CancelStatement()

Cancels a running query. To be canceled, a query must be running.

# Required Parameters
- `Id`: The identifier of the SQL statement to cancel. This value is a universally unique identifier (UUID) generated by Amazon Redshift Data API. This identifier is returned by ExecuteStatment and ListStatements.

"""

cancel_statement(Id; aws_config::AWSConfig=global_aws_config()) = redshift_data("CancelStatement", Dict{String, Any}("Id"=>Id); aws_config=aws_config)
cancel_statement(Id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("CancelStatement", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Id"=>Id), args)); aws_config=aws_config)

"""
    DescribeStatement()

Describes the details about a specific instance when a query was run by the Amazon Redshift Data API. The information includes when the query started, when it finished, the query status, the number of rows returned, and the SQL statement.

# Required Parameters
- `Id`: The identifier of the SQL statement to describe. This value is a universally unique identifier (UUID) generated by Amazon Redshift Data API. This identifier is returned by ExecuteStatment and ListStatements.

"""

describe_statement(Id; aws_config::AWSConfig=global_aws_config()) = redshift_data("DescribeStatement", Dict{String, Any}("Id"=>Id); aws_config=aws_config)
describe_statement(Id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("DescribeStatement", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Id"=>Id), args)); aws_config=aws_config)

"""
    DescribeTable()

Describes the detailed information about a table from metadata in the cluster. The information includes its columns. A token is returned to page through the column list. Depending on the authorization method, use one of the following combinations of request parameters:    AWS Secrets Manager - specify the Amazon Resource Name (ARN) of the secret and the cluster identifier that matches the cluster in the secret.    Temporary credentials - specify the cluster identifier, the database name, and the database user name. Permission to call the redshift:GetClusterCredentials operation is required to use this method.

# Required Parameters
- `ClusterIdentifier`: The cluster identifier. This parameter is required when authenticating using either AWS Secrets Manager or temporary credentials.

# Optional Parameters
- `Database`: The name of the database. This parameter is required when authenticating using temporary credentials.
- `DbUser`: The database user name. This parameter is required when authenticating using temporary credentials.
- `MaxResults`: The maximum number of tables to return in the response. If more tables exist than fit in one response, then NextToken is returned to page through the results.
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
- `Schema`: The schema that contains the table. If no schema is specified, then matching tables for all schemas are returned.
- `SecretArn`: The name or ARN of the secret that enables access to the database. This parameter is required when authenticating using AWS Secrets Manager.
- `Table`: The table name. If no table is specified, then all tables for all matching schemas are returned. If no table and no schema is specified, then all tables for all schemas in the database are returned
"""

describe_table(ClusterIdentifier; aws_config::AWSConfig=global_aws_config()) = redshift_data("DescribeTable", Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier); aws_config=aws_config)
describe_table(ClusterIdentifier, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("DescribeTable", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier), args)); aws_config=aws_config)

"""
    ExecuteStatement()

Runs an SQL statement, which can be data manipulation language (DML) or data definition language (DDL). This statement must be a single SQL statement. Depending on the authorization method, use one of the following combinations of request parameters:    AWS Secrets Manager - specify the Amazon Resource Name (ARN) of the secret and the cluster identifier that matches the cluster in the secret.    Temporary credentials - specify the cluster identifier, the database name, and the database user name. Permission to call the redshift:GetClusterCredentials operation is required to use this method.

# Required Parameters
- `ClusterIdentifier`: The cluster identifier. This parameter is required when authenticating using either AWS Secrets Manager or temporary credentials.
- `Sql`: The SQL statement text to run.

# Optional Parameters
- `Database`: The name of the database. This parameter is required when authenticating using temporary credentials.
- `DbUser`: The database user name. This parameter is required when authenticating using temporary credentials.
- `SecretArn`: The name or ARN of the secret that enables access to the database. This parameter is required when authenticating using AWS Secrets Manager.
- `StatementName`: The name of the SQL statement. You can name the SQL statement when you create it to identify the query.
- `WithEvent`: A value that indicates whether to send an event to the Amazon EventBridge event bus after the SQL statement runs.
"""

execute_statement(ClusterIdentifier, Sql; aws_config::AWSConfig=global_aws_config()) = redshift_data("ExecuteStatement", Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Sql"=>Sql); aws_config=aws_config)
execute_statement(ClusterIdentifier, Sql, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("ExecuteStatement", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Sql"=>Sql), args)); aws_config=aws_config)

"""
    GetStatementResult()

Fetches the temporarily cached result of an SQL statement. A token is returned to page through the statement results.

# Required Parameters
- `Id`: The identifier of the SQL statement whose results are to be fetched. This value is a universally unique identifier (UUID) generated by Amazon Redshift Data API. This identifier is returned by ExecuteStatment and ListStatements.

# Optional Parameters
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
"""

get_statement_result(Id; aws_config::AWSConfig=global_aws_config()) = redshift_data("GetStatementResult", Dict{String, Any}("Id"=>Id); aws_config=aws_config)
get_statement_result(Id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("GetStatementResult", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Id"=>Id), args)); aws_config=aws_config)

"""
    ListDatabases()

List the databases in a cluster. A token is returned to page through the database list. Depending on the authorization method, use one of the following combinations of request parameters:    AWS Secrets Manager - specify the Amazon Resource Name (ARN) of the secret and the cluster identifier that matches the cluster in the secret.    Temporary credentials - specify the cluster identifier, the database name, and the database user name. Permission to call the redshift:GetClusterCredentials operation is required to use this method.

# Required Parameters
- `ClusterIdentifier`: The cluster identifier. This parameter is required when authenticating using either AWS Secrets Manager or temporary credentials.

# Optional Parameters
- `Database`: The name of the database. This parameter is required when authenticating using temporary credentials.
- `DbUser`: The database user name. This parameter is required when authenticating using temporary credentials.
- `MaxResults`: The maximum number of databases to return in the response. If more databases exist than fit in one response, then NextToken is returned to page through the results.
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
- `SecretArn`: The name or ARN of the secret that enables access to the database. This parameter is required when authenticating using AWS Secrets Manager.
"""

list_databases(ClusterIdentifier; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListDatabases", Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier); aws_config=aws_config)
list_databases(ClusterIdentifier, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListDatabases", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier), args)); aws_config=aws_config)

"""
    ListSchemas()

Lists the schemas in a database. A token is returned to page through the schema list. Depending on the authorization method, use one of the following combinations of request parameters:    AWS Secrets Manager - specify the Amazon Resource Name (ARN) of the secret and the cluster identifier that matches the cluster in the secret.    Temporary credentials - specify the cluster identifier, the database name, and the database user name. Permission to call the redshift:GetClusterCredentials operation is required to use this method.

# Required Parameters
- `ClusterIdentifier`: The cluster identifier. This parameter is required when authenticating using either AWS Secrets Manager or temporary credentials.
- `Database`: The name of the database. This parameter is required when authenticating using temporary credentials.

# Optional Parameters
- `DbUser`: The database user name. This parameter is required when authenticating using temporary credentials.
- `MaxResults`: The maximum number of schemas to return in the response. If more schemas exist than fit in one response, then NextToken is returned to page through the results.
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
- `SchemaPattern`: A pattern to filter results by schema name. Within a schema pattern, \"%\" means match any substring of 0 or more characters and \"_\" means match any one character. Only schema name entries matching the search pattern are returned.
- `SecretArn`: The name or ARN of the secret that enables access to the database. This parameter is required when authenticating using AWS Secrets Manager.
"""

list_schemas(ClusterIdentifier, Database; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListSchemas", Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Database"=>Database); aws_config=aws_config)
list_schemas(ClusterIdentifier, Database, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListSchemas", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Database"=>Database), args)); aws_config=aws_config)

"""
    ListStatements()

List of SQL statements. By default, only finished statements are shown. A token is returned to page through the statement list.

# Optional Parameters
- `MaxResults`: The maximum number of SQL statements to return in the response. If more SQL statements exist than fit in one response, then NextToken is returned to page through the results.
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
- `StatementName`: The name of the SQL statement specified as input to ExecuteStatement to identify the query. You can list multiple statements by providing a prefix that matches the beginning of the statement name. For example, to list myStatement1, myStatement2, myStatement3, and so on, then provide the a value of myStatement. Data API does a case-sensitive match of SQL statement names to the prefix value you provide.
- `Status`: The status of the SQL statement to list. Status values are defined as follows:    ABORTED - The query run was stopped by the user.    ALL - A status value that includes all query statuses. This value can be used to filter results.    FAILED - The query run failed.    FINISHED - The query has finished running.    PICKED - The query has been chosen to be run.    STARTED - The query run has started.    SUBMITTED - The query was submitted, but not yet processed.
"""

list_statements(; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListStatements"; aws_config=aws_config)
list_statements(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListStatements", args; aws_config=aws_config)

"""
    ListTables()

List the tables in a database. If neither SchemaPattern nor TablePattern are specified, then all tables in the database are returned. A token is returned to page through the table list. Depending on the authorization method, use one of the following combinations of request parameters:    AWS Secrets Manager - specify the Amazon Resource Name (ARN) of the secret and the cluster identifier that matches the cluster in the secret.    Temporary credentials - specify the cluster identifier, the database name, and the database user name. Permission to call the redshift:GetClusterCredentials operation is required to use this method.

# Required Parameters
- `ClusterIdentifier`: The cluster identifier. This parameter is required when authenticating using either AWS Secrets Manager or temporary credentials.
- `Database`: The name of the database. This parameter is required when authenticating using temporary credentials.

# Optional Parameters
- `DbUser`: The database user name. This parameter is required when authenticating using temporary credentials.
- `MaxResults`: The maximum number of tables to return in the response. If more tables exist than fit in one response, then NextToken is returned to page through the results.
- `NextToken`: A value that indicates the starting point for the next set of response records in a subsequent request. If a value is returned in a response, you can retrieve the next set of records by providing this returned NextToken value in the next NextToken parameter and retrying the command. If the NextToken field is empty, all response records have been retrieved for the request.
- `SchemaPattern`: A pattern to filter results by schema name. Within a schema pattern, \"%\" means match any substring of 0 or more characters and \"_\" means match any one character. Only schema name entries matching the search pattern are returned. If SchemaPattern is not specified, then all tables that match TablePattern are returned. If neither SchemaPattern or TablePattern are specified, then all tables are returned.
- `SecretArn`: The name or ARN of the secret that enables access to the database. This parameter is required when authenticating using AWS Secrets Manager.
- `TablePattern`: A pattern to filter results by table name. Within a table pattern, \"%\" means match any substring of 0 or more characters and \"_\" means match any one character. Only table name entries matching the search pattern are returned. If TablePattern is not specified, then all tables that match SchemaPatternare returned. If neither SchemaPattern or TablePattern are specified, then all tables are returned.
"""

list_tables(ClusterIdentifier, Database; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListTables", Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Database"=>Database); aws_config=aws_config)
list_tables(ClusterIdentifier, Database, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = redshift_data("ListTables", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClusterIdentifier"=>ClusterIdentifier, "Database"=>Database), args)); aws_config=aws_config)


text:
module DesignStructureMatrix

using Gadfly
using LinearAlgebra

export plotDSM
export toReachableMatrix
export OrderReachable
export Sequencing
export Clustering

## plot DSM by using Gadfly
function plotDSM(DSM,label)

    cDSM = copy(DSM);
    cDSM[diagind(cDSM)] .= -1

     spy(cDSM, Scale.y_discrete(labels = i->label[i]), Scale.x_discrete(labels = i->label[i]),
        Guide.ylabel(nothing), Guide.xlabel(nothing, orientation=:vertical), Guide.xticks(orientation=:vertical),
        Scale.color_continuous(colormap=Scale.lab_gradient( "darkgrey","aliceblue", "navy")),
        Theme(minor_label_font_size=16pt, key_position=:none, bar_spacing= 0.5mm, panel_fill="darkgrey"))

end

## convert DSM to ReachableMatrix
function toReachableMatrix(DSM)

    I = zeros(Int64,size(DSM)) + Diagonal(ones(Int64,size(DSM))); #identity matrix
    R1 = zeros(Int64,size(DSM))
    calcmax = 1000

    for i  = 1 : calcmax

        R1 = (DSM+I)^i
        (Rindex) =findall(x-> 1<x , R1);
        R1[Rindex] = ones(size(Rindex));

        R2 = (DSM+I)^(i+1)
        (Rindex) =findall(x-> 1<x , R2);
        R2[Rindex] = ones(size(Rindex));

        if R1 == R2
            break
        end

    end

    return R1
end

##order Reachable matrix for seaqencing
function OrderReachable(Reachable)

    original = copy(Reachable);
    (DSMsize,) = size(Reachable);
    valuesDSM = collect(1:DSMsize);
    count = 0;
    level = zeros(Int,DSMsize);

    while count < DSMsize

        R = Reachable;
        A = Reachable';
        RA = R .* A;
        (RAsize,) = size(RA);
        elements = collect(1:RAsize);

        for i = 1 : RAsize

            if RA[i,:] == R[i,:];
                    count += 1;
                    level[count] = valuesDSM[i];
                    elements[i] = 0;
            end
        end

        elements = filter(x-> x > 0, elements);
        Reachable = Reachable[elements , elements];
        valuesDSM = valuesDSM[elements];

    end

    orderedReachable = original[level,level]

    return orderedReachable , level

end

## Sequencing DSM
function Sequencing(DSM,label)
    cDSM = copy(DSM);
    clabel = copy(label);

    Reachable = toReachableMatrix(cDSM);

    (OrderedReachable,level) = OrderReachable(Reachable);

    SequencedDSM = cDSM[level,level];
    Sequencedlabel = clabel[level];

    return SequencedDSM, Sequencedlabel

end

##Clustering DSM
function Clustering(DSM,label)

    powcc = 1;
    powbid = 0;
    powdep = 1;
    rand_accept = 30;
    rand_bid = 30;
    times = 5;
    stable_limit = 2;

    original_DSM = copy(DSM);
    original_label = copy(label);
    (DSMsize,) = size(original_DSM);

    Clustermatrix = zeros(Int64,DSMsize,DSMsize) + Diagonal(ones(Int64,DSMsize,DSMsize)) ;

    cDSM = copy(original_DSM);

    function ClusterBid(Clustermatrix,cDSM,DSMsize,powdep,powbid,element)

        bestCluster = 0;
        bestBid = -1;

        for i = 1 : DSMsize

            Clusterlist= findall(x-> x == 1 , Clustermatrix[i,:]);
            (Clustersize,) = size(Clusterlist);
            bid = 0;

            for j = 1 : Clustersize

                bid += cDSM[element,Clusterlist[j]];

            end

            bid = bid ^ powdep / Clustersize ^ powbid;

            if bid > bestBid

                 bestBid = bid;
                bestCluster = i;
            end

        end

        return bestCluster

    end

    function TCC(cDSM,DSMsize,Clustermatrix,powcc)

        intraCost=0;
        extraCost=0;

        for i = 1: DSMsize

            for j = 1 : DSMsize

                Clusterofi = findall(x-> x == 1 , Clustermatrix[:,i]);
                Clusterofj = findall(x-> x == 1 , Clustermatrix[:,j]);
                (Clustersizeofi,) = size(Clustermatrix[Clusterofi,:]);

                cost = cDSM[i,j] + cDSM[j,i];

                if Clusterofi == Clusterofj

                    intraCost += cost*Clustersizeofi^powcc;

                else

                    extraCost += cost*DSMsize^powcc;
                end

            end

        end

        totalCost = intraCost + extraCost;

        return totalCost
    end

    TCCost = TCC(cDSM,DSMsize,Clustermatrix,powcc)

    costlist = zeros(DSMsize*times,1);

    for i = 1 : DSMsize * times;

        preClustermatrix = Clustermatrix;
        element = rand(1:DSMsize);
        bestCluster = ClusterBid(Clustermatrix,cDSM,DSMsize,powdep,powbid,element);
        Clustermatrix[:,element] = zeros(DSMsize,1);
        Clustermatrix[bestCluster,element] = 1;
        newCost = TCC(cDSM,DSMsize,Clustermatrix,powcc);
        if newCost <= maximum([TCCost,rand_accept])
            TCCost = newCost;
        else
            Clustermatrix = preClustermatrix;
        end

        costlist[i] = newCost
    end

    orderCluster = [sum(Clustermatrix,dims=2) Clustermatrix]
    orderedCluster = sortslices(orderCluster,dims = 1, rev=true)
    reorderedCluster = orderedCluster[:,2:DSMsize+1]

    Order = [];

    for i = 1 : DSMsize
        Order = [Order ; findall(x-> x == 1 , reorderedCluster[i,:])]
    end

    ## reorder DSM
    I =zeros(DSMsize,DSMsize) + Diagonal(ones(DSMsize,DSMsize))
    P =zeros(DSMsize,DSMsize) + Diagonal(ones(DSMsize,DSMsize))
    Clusteredlabel = copy(original_label)
    copylabel= copy(original_label)

    ## Permutation matrix
    for i = 1 : DSMsize
        P[i,:] = I[Order[i],:]
        Clusteredlabel[i] = copylabel[Order[i]]
    end

    ## order DSM
    ClusteredDSM = P*original_DSM*P'

    return ClusteredDSM, Clusteredlabel;

end

end # module


text:
# A struct that contain all the data required for converging ep
struct EPModel{T<:Real}
    scalefact::T
    updatealg!::Function
    epfields::EPFields{T}
    epmat::AbstractEPMat
    alpha::T
    beta_vec::SparseVector{T, Int}
    stat::Dict # Just data about the ep run time
end

function EPModel(
        S::AbstractArray{T,2}, b::AbstractArray{T}, lb::AbstractArray{T}, ub::AbstractArray{T};
        alpha::Real=Inf,                                        # inverse temperature
        beta_vec::AbstractVector{T} = spzeros(T, size(S, 2)),   # maxent inverse temperature vector
        solution::Union{EPOut{T}, Nothing} = nothing,           # start from a solution
        expval = nothing                                        # fix posterior probability experimental values for std and mean
    ) where {T<:Real}

    # Some checks
    M, N = size(S)
    M > N && @warn("M = $M ≥ N = $N")
    any(lb .> ub) && error("lower bound fluxes > upper bound fluxes. Consider swapping lower and upper bounds")

    # The scalefactor is just the maximum absolute bound (lb or ub).
    scalefact = get_scalefactor(lb, ub)

    # Create EPFields. If a solution is not given, the EPfields will be fresh
    epfields::EPFields = isnothing(solution) ? EPFields(N, expval, eltype(S)) :
        deepcopy(solution.sol) # preserve the original solution!

    # making a local copy to rescale
    lb, ub, b = copy.([lb, ub, b])

    #=
    Scale down μ, s, av, va of epfields and ub, lb and Y using the
    previous computed scalefactor.
    If epfields is fresh, it only will have effect on ub, lb and Y
    =#
    scaleepfield!(inv(scalefact), epfields, ub, lb, b) # scaling fields to [0,1]

    epmat = (alpha < Inf) ? EPMat(S, b, lb, ub, alpha) : EPMatT0(S, b, lb, ub)

    # One iteration of EP
    updatealg! = isinf(alpha) ? eponesweepT0! : eponesweep!

    beta_vec = prepare_beta_vec(epmat, beta_vec)

    return EPModel{T}(scalefact, updatealg!, epfields, epmat, alpha, beta_vec, Dict())

end

EPModel(metnet::MetNet; kwargs...) = EPModel(metnet.S, metnet.b, metnet.lb, metnet.ub; kwargs...)


text:
using PartialRejectionSampling
using Random

using Plots
using GraphPlot, Colors

using LightGraphs
const LG = LightGraphs

using Cairo, Compose

function plot(
    graph::LG.AbstractGraph,
    dims::Vector{Int}=zeros(Int, 2),
    path="";
    kwargs...
)
    if isempty(path)
        if any(dims .== 0)
            p = GraphPlot.gplot(graph;
                    nodelabel=LG.vertices(graph),
                    kwargs...
                    )
            return p
        else
            @assert prod(dims) == LG.nv(graph)
            pos = collect(Iterators.product(1:dims[2], 1:dims[1]))[:]
            locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)

            p = GraphPlot.gplot(graph,
                                locs_x,
                                reverse(locs_y);
                                nodelabel=LG.vertices(graph),
                                kwargs...
                )
            return p
        end
    else
        if any(dims .== 0)
            Compose.draw(
                PDF(path, 16cm, 16cm),
                GraphPlot.gplot(graph;
                                nodelabel=LG.vertices(graph),
                                kwargs...
                )
            )
        else
            @assert prod(dims) == LG.nv(graph)
            pos = collect(Iterators.product(1:dims[2], 1:dims[1]))[:]
            locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)

            Compose.draw(
                PDF(path, 16cm, 16cm),
                GraphPlot.gplot(graph,
                                locs_x,
                                reverse(locs_y);
                                nodelabel=LG.vertices(graph),
                                kwargs...
                )
            )
        end
    end
end

function plot(
    pp::PRS.Ising,
    state,
    dims=zeros(2),
    path="";
    kwargs...
)
    c_nodes = ifelse.(
            state .== 1,
            Colors.colorant"gray",
            Colors.colorant"white")

    plot(pp.graph, dims, path; nodefillc=c_nodes, kwargs...)
end

function plot(
    pp::PRS.HardCoreGraph,
    state,
    dims=zeros(2),
    path="";
    kwargs...
)
    c_nodes = [Colors.colorant"white" for _ in 1:LG.nv(pp.graph)]
    c_nodes[state] .= Colors.colorant"grey"

    plot(pp.graph, dims, path; nodefillc=c_nodes, kwargs...)
end

function plot(
    pp::PRS.SinkFreeGraph,
    sample,
    dims=zeros(2),
    path="";
    kwargs...
)
    c_nodes, c_edges = color_sinks(sample)

    plot(sample, dims, path; nodefillc=c_nodes, edgestrokec=c_edges, kwargs...)
end

function plot(
    pp::PRS.RootedSpanningForest,
    sample,
    dims=zeros(2),
    path="";
    kwargs...
)
    c_nodes, c_edges = color_cycles(sample)
    c_nodes[collect(pp.roots)] .= Colors.colorant"orange"

    plot(sample, dims, path; nodefillc=c_nodes, edgestrokec=c_edges)
end

function color_cycles(
    graph::LG.SimpleDiGraph{T},
    c_node=Colors.colorant"turquoise",
    c_node_cycle=Colors.colorant"red",
    c_edge=Colors.colorant"lightgray",
    c_edge_cycle=Colors.colorant"red"
) where {T}

    edge_map = PRS.edgemap(graph)
    edge_idx = T[]
    nodes = T[]
    for cycle in LG.simplecycles(graph)
        union!(nodes, cycle)
        for (x, y) in zip(cycle, circshift(cycle, -1))
            push!(edge_idx, edge_map[LG.Edge(x, y)])
        end
    end

    c_nodes = fill(c_node, LG.nv(graph))
    c_nodes[nodes] .= c_node_cycle

    c_edges = fill(c_edge, LG.ne(graph))
    c_edges[edge_idx] .= c_edge_cycle

    return c_nodes, c_edges
end

function color_sinks(
    graph::LG.SimpleDiGraph{T},
    c_node=Colors.colorant"turquoise",
    c_node_sink=Colors.colorant"red",
    c_edge=Colors.colorant"lightgray",
    c_edge_sink=Colors.colorant"red"
) where {T}

    nodes = T[]
    edge_idx = T[]
    edge_map = PRS.edgemap(graph)
    for v in PRS.sink_nodes(graph)
        push!(nodes, v)
        for w in LG.inneighbors(graph, v)
            push!(edge_idx, edge_map[LG.Edge(w, v)])
        end
    end

    c_nodes = fill(c_node, LG.nv(graph))
    c_nodes[nodes] .= c_node_sink

    c_edges = fill(c_edge, LG.ne(graph))
    c_edges[edge_idx] .= c_edge_sink

    return c_nodes, c_edges
end

function hard_core_graph_pedagogy(pp, dims, rng)

    path(i) = joinpath(@__DIR__,
                       "output",
                       "hard_core_graph",
                       join([lpad(i, 3, "0"), ".pdf"]))

    c_node_normal = Colors.colorant"turquoise";

    c_node_unccupied = Colors.colorant"white";
    c_node_occupied = Colors.colorant"gray";

    c_node_update = Colors.colorant"orange";
    c_node_bad = Colors.colorant"red";

    i = 0
    c_nodes = fill(c_node_normal, LG.nv(pp.graph))
    plot(pp.graph, dims, path(i); nodefillc=c_nodes)

    proba = pp.β / (one(pp.β) + pp.β)

    adj = LG.adjacency_matrix(pp.graph)
    occupied = Random.randsubseq(rng, LG.vertices(pp.graph), proba)

    while true

        i += 1
        c_nodes = fill(c_node_unccupied, LG.nv(pp.graph))
        c_nodes[occupied] .= c_node_occupied
        plot(pp.graph, dims, path(i); nodefillc=c_nodes)

        # Check if occupied vertices form an independent set
        sub_graph = LG.SimpleGraph(adj[occupied, occupied])
        LG.ne(sub_graph) == 0 && break

        bad = []
        independent, resample = [], Set()
        for cc in LG.connected_components(sub_graph)
            if length(cc) == 1  # Identify current independent vertices
                append!(independent, occupied[cc])
            else  # Construct the resampling set of vertices
                union!(resample, occupied[cc])
                union!(bad, occupied[cc])
                for v in cc
                    union!(resample, LG.neighbors(pp.graph, occupied[v]))
                end
            end
        end

        i += 1
        c_nodes[bad] .= c_node_bad
        plot(pp.graph, dims, path(i); nodefillc=c_nodes)

        Random.randsubseq!(rng, occupied, collect(resample), proba)
        append!(occupied, independent)

        i += 1
        c_nodes[collect(resample)] .= c_node_update;
        plot(pp.graph, dims, path(i); nodefillc=c_nodes)
    end
end

function rooted_spanning_tree_padagogy(pp, dims, rng)

    c_node_roots = Colors.colorant"lightgreen";
    c_node_normal = Colors.colorant"turquoise";
    c_node_update = Colors.colorant"orange";
    c_node_bad = Colors.colorant"red";

    c_edge_normal = Colors.colorant"lightgray";
    c_edge_update = Colors.colorant"orange";
    c_edge_bad = Colors.colorant"red";

    edge_map = PRS.edgemap(pp.graph)

    path(i) = joinpath(@__DIR__,
                       "output",
                       "rooted_spanning_tree",
                       join([lpad(i, 3, "0"), ".pdf"]))

    i = 0
    c_nodes = fill(c_node_normal, LG.nv(pp.graph))
    c_nodes[roots] .= c_node_roots
    plot(pp.graph, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edge_normal)

    i += 1
    g = PRS.random_neighbor_assignment(rng, pp.graph, roots)
    plot(pp.graph, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edge_normal)

    while true
        i += 1
        c_nodes, c_edges = color_cycles(pp.graph, c_node_normal, c_node_bad, c_edge_normal, c_edge_bad)
        c_nodes[roots] .= c_node_roots
        plot(pp.graph, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edges)

        vertices_in_cycles = Set(Iterators.flatten(LG.simplecycles(pp.graph)))
        isempty(vertices_in_cycles) && break

        # Resample the successor of vertices involved in cycles
        edges = []
        c_nodes = fill(c_node_normal, LG.nv(pp.graph))
        c_nodes[roots] .= c_node_roots
        c_edges = fill(c_edge_normal, LG.ne(pp.graph))

        for v in vertices_in_cycles
            c_nodes[v] = c_node_update
            # Remove current edge (v, w)
            w = LG.neighbors(pp.graph, v)[1]
            LG.rem_edge!(pp.graph, v, w)
            # Resample the successor w of v and add edge (v, w)

            w = rand(rng, LG.neighbors(pp.graph, v))
            LG.add_edge!(pp.graph, v, w)

            push!(edges, LG.Edge(v, w))
        end

        i += 1
        edge_map = PRS.edgemap(pp.graph)
        c_edges[[edge_map[e] for e in edges]] .= c_edge_update
        plot(pp.graph, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edges)
    end
end

function sink_free_graph_pedagogy(pp, dims, rng)
    c_node_normal = Colors.colorant"turquoise";
    c_node_update = Colors.colorant"orange";
    c_node_bad = Colors.colorant"red";

    c_edge_normal = Colors.colorant"lightgray";
    c_edge_update = Colors.colorant"orange";
    c_edge_bad = Colors.colorant"red";

    edge_map = PRS.edgemap(pp.graph)

    path(i) = joinpath(@__DIR__,
                       "output",
                       "sink_free_graph",
                       join([lpad(i, 3, "0"), ".pdf"]))

    i = 0
    plot(pp.graph, dims, path(i); nodefillc=c_node_normal, edgestrokec=c_edge_normal)

    i += 1
    g = PRS.random_edge_orientation(rng, pp.graph)
    plot(g, dims, path(i); nodefillc=c_node_normal, edgestrokec=c_edge_normal)

    while true

        i+=1
        c_nodes, c_edges = color_sinks(g, c_node_normal, c_node_bad, c_edge_normal, c_edge_bad)
        plot(g, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edges)

        sinks = PRS.sink_nodes(g)
        isempty(sinks) && break

        # Resample i.e. flip edges forming a sink uniformly at random (Bernoulli(0.5))
        c_nodes = fill(c_node_normal, LG.nv(g))
        c_edges = fill(c_edge_normal, LG.ne(g))

        edges = []
        for v in sinks
            c_nodes[v] = c_node_update
            inneigh_v = copy(LG.inneighbors(g, v))
            for w in inneigh_v
                if rand(rng) < 0.5
                    LG.rem_edge!(g, w, v)
                    LG.add_edge!(g, v, w)
                    push!(edges, LG.Edge(v, w))
                else
                    push!(edges, LG.Edge(w, v))
                end
            end
            edge_map = PRS.edgemap(g)
            c_edges[[edge_map[e] for e in edges]] .= c_edge_update
        end
        i += 1
        plot(g, dims, path(i); nodefillc=c_nodes, edgestrokec=c_edges)
    end
end


text:
using DelimitedFiles

# This is a function that reads in the Hutchinson operator and corresponding
#     probabilities and outputs a randomly selected transform
# This works by choosing a random number and then iterating through all
#     probabilities until it finds an appropriate bin
function select_array(hutchinson_op, probabilities)

    # random number to be binned
    rnd = rand()

    # This checks to see if a random number is in a bin, if not, that
    #     probability is subtracted from the random number and we check the
    #     next bin in the list
    for i = 1:length(probabilities)
        if (rnd < probabilities[i])
            return hutchinson_op[i]
        end
        rnd -= probabilities[i]
    end
end

# This is a general function to simulate a chaos game
# n is the number of iterations
# initial_location is the starting point of the chaos game
# hutchinson_op is the set of functions to iterate through
# probabilities is the set of probabilities corresponding to the likelihood
#     of choosing their corresponding function in hutchinson_op
function chaos_game(n::Int, initial_location, hutchinson_op, probabilities)

    # Initializing the output array and the initial point
    output_points = zeros(n,2)

    # extending point to 3D for affine transform
    point = [initial_location[1], initial_location[2], 1]

    for i = 1:n
        output_points[i,:] .= point[1:2]
        point = select_array(hutchinson_op, probabilities)*point
    end

    return output_points

end

barnsley_hutchinson = [[0.0 0.0 0.0;
                        0.0 0.16 0.0;
                        0.0 0.0 1.0],
                       [0.85 0.04 0.0;
                        -0.04 0.85 1.60;
                        0.0 0.0 1.0],
                       [0.20 -0.26 0.0;
                        0.23 0.22 1.60;
                        0.0 0.0 1.0],
                       [-0.15 0.28 0.0;
                        0.26 0.24 0.44;
                        0.0 0.0 1.0]]

barnsley_probabilities = [0.01, 0.85, 0.07, 0.07]
output_points = chaos_game(10000, [0,0],
                           barnsley_hutchinson, barnsley_probabilities)
writedlm("out.dat", output_points)


text:
#=
   Gender height.

   Identify a person's sex by height.

   Cf ~/cplint/gender_height.pl
      ~/blog/gender_height.blog
      ~/webppl/gender_height.wppl

   Note: The observations of height and gender seems to be much better
         using Turings standard "parameter observation" than using
         Dirac etc.

=#

# using Memoization
using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)
include("jl_utils.jl")

@model function gender_height(gender=missing,height=missing)
# @model function gender_height()

    male = 1
    female = 2
    gender ~ Categorical([0.5,0.5])

    # From https://en.wikipedia.org/wiki/List_of_average_human_height_worldwide
    # Here are the values of Sweden. I'm not sure about the variance of these heights...
    height ~ gender == male ? Normal(181.5,sqrt(50)) : Normal(166.8,sqrt(50))

    # height ~ gender == male ? Normal(178,17.7) : Normal(163,17.3)

end

function run_model(gender=missing,height=missing)

    model = gender_height(gender,height)

    num_chns = 4

    # chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)

    chns = sample(model, MH(), 10_000)
    # chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)

    # chns = sample(model, PG(20), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, PG(20), 1_000)

    # chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, IS(), 10_000)

    # chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
    # chns = sample(model, SMC(1000), 10_000)

    # chns = sample(model, NUTS(1000,0.65), 1_000)
    # chns = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)
    # chns = sample(model, Gibbs(MH(:gender),NUTS(10,0.65,:height)), 1_000)
    # chns = sample(model, Gibbs(MH(:gender),HMC(0.1,5,:height)), 1_000)
    # chns = sample(model, Gibbs(PG(10,:gender),HMC(0.1,5,:height)), 1_000)
    # chns = sample(model, Gibbs(MH(:gender),NUTS(1_000,0.65,:height)), 1_000)

    if ismissing(height)
        display(chns)
        # display(plot(chns))
    end
    show_var_dist_pct(chns, :gender)
    show_var_dist_pct(chns, :height, 20)
end

male = 1
female = 2
println("Male is coded as 1, female as 2.")
for height in [150.0,160.0,170.0,175.0,180.0, 185.0,190.0,200.0]
    println("\ngender:missing height:$height")
    run_model(missing,height)
    println()
end

println("\nmale, missing")
run_model(male,missing)

println("\nfemale, missing")
run_model(female,missing)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
## Transition cost object for continuous time approximation

struct TransitionCostContinuous{T}
    t::Vector{T}
    I_g::Vector{T}
    I_g2::Vector{T}
    I_tg::Vector{T}
end

function TransitionCostContinuous{T}(g, t::AbstractVector, tol=1e-3) where {T}
    N = length(t)
    I_g =  fill(zero(T), N)
    I_g2 = fill(zero(T), N)
    I_tg = fill(zero(T), N)
    for i=2:N
        I_g[i] = I_g[i-1] + quadgk(g, t[i-1], t[i], rtol=tol)[1]
        I_g2[i] = I_g2[i-1] + quadgk(τ -> g(τ)^2, t[i-1], t[i], rtol=tol)[1]
        I_tg[i] = I_tg[i-1] + quadgk(τ -> τ*g(τ), t[i-1], t[i], rtol=tol)[1]
    end
    TransitionCostContinuous{T}(t, I_g, I_g2, I_tg)
end

# getindex returns a quadratic form H([x1, x2]) that represents the
# transition cost from i to ip
function getindex(tc::TransitionCostContinuous, i::Integer, ip::Integer)
    t, I_g, I_g2, I_tg = tc.t, tc.I_g, tc.I_g2, tc.I_tg

    P = @SMatrix [(t[ip]-t[i])/3 (t[ip]-t[i])/6; (t[ip]-t[i])/6 (t[ip]-t[i])/3]

    q = -2* 1/(t[ip]-t[i]) *
    @SVector [-(I_tg[ip] - I_tg[i]) + t[ip]*(I_g[ip] - I_g[i]),
    (I_tg[ip] - I_tg[i]) - t[i]*(I_g[ip] - I_g[i])]

    r =  I_g2[ip] - I_g2[i]

    return QuadraticForm(P, q, r)
end

Base.size(tc::TransitionCostContinuous) = (length(tc.t)-1, length(tc.t))
Base.size(tc::TransitionCostContinuous, i::Integer) = size(tc)[i]

## Transition cost object for continuous time approximation

# Test non-uniform sampling
# non-uniform sampling needs to be handled

struct TransitionCostDiscrete{T,TimeType}
    t::TimeType
    P_matrices::Vector{SMatrix{2,2,T,4}}
    G1::Vector{T}
    G2::Vector{T}
    G3::Vector{T}
end

function TransitionCostDiscrete{T}(g::AbstractArray, t::TimeType=1:length(g)) where {T, TimeType<:AbstractArray{<:Integer}}
    N = length(g)

    # Find sums of g, k*g, and g^2
    G1 = zeros(T, N)
    G2 = zeros(T, N)
    G3 = zeros(T, N)
    # The sums corresponding to transitioning from [1, ip)
    # i.e. not including the cost at ip
    for k=2:N
        G1[k] = G1[k-1] + g[k-1]
        G2[k] = G2[k-1] + (k-1)*g[k-1]
        G3[k] = G3[k-1] + g[k-1]^2
    end

    # The P-matrices only depend on the distance d=ip-i
    P_matrices  = Vector{SMatrix{2,2,T,4}}(undef, N-1)
    P_matrices[1] = @SMatrix [1.0 0; 0 0]
    for d=2:N-1
        off_diag_elems = sum([k*(d - k) for k=0:d-1])
        P_matrices[d] = @SMatrix [P_matrices[d-1][1,1] + d^2    off_diag_elems;
        off_diag_elems            P_matrices[d-1][1,1]]
    end

    P_matrices = P_matrices ./ (1.0:N-1).^2 # FIXME: Some problem here
    TransitionCostDiscrete{T,TimeType}(t, P_matrices, G1, G2, G3)
end

# getindex returns a quadratic form H([x1, x2]) that represents the
# transition cost from i to ip
function getindex(tc::TransitionCostDiscrete, i::Integer, ip::Integer)
    t_i = tc.t[i]
    t_ip = tc.t[ip]
    q = -2* 1/(t_ip-t_i) *
    @SVector [-(tc.G2[t_ip] - tc.G2[t_i]) + t_ip*(tc.G1[t_ip] - tc.G1[t_i]),
    (tc.G2[t_ip] - tc.G2[t_i]) - t_i*(tc.G1[t_ip] - tc.G1[t_i])]

    r =  tc.G3[t_ip] - tc.G3[t_i]

    return QuadraticForm(tc.P_matrices[t_ip-t_i], q, r)
end

Base.size(tc::TransitionCostDiscrete) = (length(tc.t)-1, length(tc.t))
Base.size(tc::TransitionCostDiscrete, i::Integer) = size(tc)[i]


text:
"""
    institution_name(x)
    institution_name!(x, val)

Name of the institution in charge of the equipment that produced the composite instances.
"""
@defprop InstitutionName{:institution_name}::String

"""
    institutional_department_name(x)
    institutional_department_name!(x, val)

The department in the institution in charge of the equipment that produced the composite instances.
"""
@defprop InstitutionalDepartmentName{:institutional_department_name}::String

"""
    institution_address(x)
    institution_address!(x, val)

The address of the institution in charge of the equipment that produced the composite instances.
"""
@defprop InstitutionAddress{:institution_address}::String

"""
    InstitutionInformation

Metadata structure for general MRI sequence information.

## Supported Properties

$(description_list(institution_name,institutional_department_name,institution_address))

## Examples
```jldoctest
julia> using NeuroCore

julia> using NeuroCore.NeuroMetadata

julia> m = InstitutionInformation("a", "b", "c")
InstitutionInformation("a", "b", "c")

julia> m.institution_name
"a"

julia> m.institution_address
"c"

julia> m.institutional_department_name
"b"
```
"""
struct InstitutionInformation
    name::String
    department::String
    address::String
end

@properties InstitutionInformation begin
    institution_name(self) => :name
    institutional_department_name(self) => :department
    institution_address(self) => :address
end


text:
using Test
using NQCDynamics
using NQCDynamics: DynamicsMethods, Calculators
using NQCDynamics.DynamicsMethods.EhrenfestMethods
using OrdinaryDiffEq

@test Ehrenfest{Float64}(2) isa Ehrenfest
atoms = Atoms(:H)

@testset "Ehrenfest" begin
    sim = Simulation{Ehrenfest}(atoms, NQCModels.DoubleWell())

    r = zeros(size(sim))
    v = randn(size(sim))
    u = DynamicsVariables(sim, v, r, PureState(1, Adiabatic()))
    du = zero(u)

    @test DynamicsUtils.get_quantum_subsystem(u) ≈ Complex.([1 0; 0 0])

    DynamicsMethods.motion!(du, u, sim, 0.0)

    Calculators.evaluate_nonadiabatic_coupling!(sim.calculator, r)
    @test sim.calculator.nonadiabatic_coupling ≈ -sim.calculator.nonadiabatic_coupling'

    problem = ODEProblem(DynamicsMethods.motion!, u, (0.0, 1.0), sim)
    integrator = init(problem, Tsit5())

    σ = DynamicsUtils.get_quantum_subsystem(u)
    dσ = zero(σ)
    DynamicsUtils.set_quantum_derivative!(dσ, v, σ, sim)

    @testset "get_diabatic_population" begin
        population = Estimators.diabatic_population(sim, u)
        @test population ≈ [0.5, 0.5]
    end

    sim1 = Simulation{Ehrenfest}(atoms, NQCModels.TullyModelOne())
    r = fill(-5.0, size(sim))
    v1 = fill(8.9, size(sim)) ./ sim1.atoms.masses[1]
    z1 = DynamicsVariables(sim, v1, r, PureState(1, Adiabatic()))
    solution1 = run_trajectory(z1, (0.0, 2500.0), sim1; output=(:population), reltol=1e-6)

    @testset "run_trajectory" begin
        atoms = Atoms(2000)
        sim = Simulation{Ehrenfest}(atoms, NQCModels.TullyModelTwo())
        v = hcat(100 / 2000)
        r = hcat(-10.0)
        u = DynamicsVariables(sim, v, r, PureState(1, Adiabatic()))
        solution = run_trajectory(u, (0.0, 500.0), sim, output=(:hamiltonian), reltol=1e-6)
        @test solution.hamiltonian[1] ≈ solution.hamiltonian[end] rtol=1e-2
    end
end

@testset "Ehrenfest RPMD" begin
    atoms = Atoms(2000)
    sim = RingPolymerSimulation{Ehrenfest}(atoms, NQCModels.TullyModelTwo(), 10)
    v = fill(100 / 2000, 1, 1, 10)
    r = fill(-10.0, 1, 1, 10)
    u = DynamicsVariables(sim, v, r, PureState(1, Adiabatic()))
    solution = run_trajectory(u, (0.0, 500.0), sim, output=(:hamiltonian), dt=1)
    @test solution.hamiltonian[1] ≈ solution.hamiltonian[end] rtol=1e-2
end