text:
#################################################################################
#coeffHist=map(i->1/sum(pdf.(Geometric(.2),rand(Geometric(.2),20))),1:10)
#compare Array elements
#return { x : x∈A ⩓ x∈B }|indices of such x
#δ==Int64(2) → return { x : x∈B ⩓ x∉A }|indices of such x
# function compArr(A,B,δ=Int64(1))
# 	if δ==Int64(2)
# 		return B[findall(x->x∉A,B)],findall(x->x∉A,B)
# 	end
# 	return A[findall(x->x∈B,A)],findall(x->x∈B,A)
# end
using PDMats
using LinearAlgebra
using Distributions

function relabel(X)
	tmp=deepcopy(X)
	uni_tmp=unique(tmp)
	len_tmp=length(uni_tmp)
	map(i->tmp[tmp.==uni_tmp[i]].=10*i,1:len_tmp)
	return tmp
end

function countUnique(X)
	tmp1=unique(X)
	tmp2=map(i->sum(X.==tmp1[i]),1:length(tmp1))
	return DataFrame(unique=tmp1,count=tmp2)
end

function cartIDX(IDX_Cartesian_All::Array{Tuple{Int64,Int64},1},IDX::Int64)
	IDX_Bool=map(i->IDX_Cartesian_All[i]  == IDX_Cartesian_All[IDX],1:length(IDX_Cartesian_All))
	return IDX_Bool
end
function cartIDX(IDX_Cartesian_All::Array{Tuple{Int64,Int64},1},IDX_Cartesian::Tuple{Int64,Int64})
	IDX_Bool=map(i->IDX_Cartesian_All[i]  == IDX_Cartesian,1:length(IDX_Cartesian_All))
	return IDX_Bool
end

function heatSR(X::Array{Array{Float64,2},2})
	n=size(X,1)
	H=map(k->map(i->length(X[k,i]),1:n),1:n)
	H=reshape(vcat(H...),(n,n))
	return H
end
function dim(A)
	d1=length(A[:,1])
	d2=length(A[1,:])
	return d1,d2
end
function eye(d)
	E=zeros(Float64,d,d)
	map(i->E[i,i]=1,1:dim(E)[1])
	return E
end
function ctrans(X::Array{Int64,2})
	convert(Array{Float64,2},transpose(X))
end
function ctrans(X::Array{Float64,2})
	convert(Array{Float64,2},transpose(X))
end
function ctrans(X::Array{Array{Float64,2},1})
	convert(Array{Array,1},map(i->ctrans(𝛍[1][i])[:],collect(1:length(𝛍[1]))))
end
function dim2clean(MU_star::Array{Array,1})
	MU_star = dim2adj(MU_star)
	N=maximum(size(MU_star))
	𝛅=Int64(3)
	for i=1:N
		if isnan(MU_star[i])
			for j=1:𝛅
				MU_star[i,j]=-10.0
			end
		end
	end
	MU_star = dim2adj(MU_star)
	MU_star = MU_star[MU_star .> dim2adj([100.0,100.0,100.0])]
end
function dim2clean(MU_star::Array{Array,1},𝛆::Array{Float64,1})
	MU_star = dim2adj(MU_star)
	N=maximum(size(MU_star))
	𝛅=Int64(3)
	for i=1:N
		if isnan(MU_star[i])
			for j=1:𝛅
				MU_star[i,j]=-10.0
			end
		end
	end
	MU_star = dim2adj(MU_star)
	MU_star = MU_star[MU_star .>= dim2adj(𝛆)]
end
function dim2crop(X::Array{Float64,2},xmin::Float64,xmax::Float64,ymin::Float64,ymax::Float64)
	Xog = X
	X=X[X[:,1].>=xmin,:]
	X=X[X[:,1].<=xmax,:]
	X=X[X[:,2].>=ymin,:]
	X=X[X[:,2].<=ymax,:]
	val=(Xog[:,1] .∈ [X[:,1]]) .& (Xog[:,2] .∈ [X[:,2]])
	return[X,val]
end
function dim2adj(X::Array{Float64,2})
	J = size(X,1)
	MU_0 = Array{Array,1}(undef,J)
	for j=1:J
		MU_0[j] = X[j,:]
	end
	return MU_0
end
function dim2adj(X::Array{Float64,1})
	MU_0 = Array{Array,1}(undef,1)
	MU_0[1] = X[:,1]
	return MU_0
end
function dim2adj(X::Array{Array,1},δ=Int64(2))
	if δ != Int64(2)
		Y_raw = Array{Float64,2}(undef,maximum(size(X)),size(X[1],1))
		for i=1:maximum(size(X))
			for j=1:size(X[1],1)
				Y_raw[i,j] = X[i][j]
			end
		end
		return Y_raw
	end
	Y_raw = Array{Float64,2}(undef,maximum(size(X)),2)
	for i=1:maximum(size(X))
		Y_raw[i,1] = X[i][1]
		Y_raw[i,2] = X[i][2]
	end
	return Y_raw
end
####Generate MixtureModel off given data
function dim2mix(X::Array{Array,1},S=eye(size(X,2)))
	MM = MixtureModel(map(x->MvNormal(x,S),X))
	return MM
end
function dim2mix(X::Array{Float64,2},S=eye(size(X,2)))
	X = dim2adj(X)
	MM = MixtureModel(map(x->MvNormal(x,S),X))
	return MM
end
function dim2mix(X::Array{Float64,2},w::Array{Float64,1},S=eye(size(X,2)))
	#print("[dim2mix] X=",X,"\n")
	X = dim2adj(X)
	#print("[dim2mix..dim2adj(X)] X=",X,"\n")
	#print("[dim2mix] S=",S,"\n")
	MM = MixtureModel(map(x->MvNormal(x,S),X))
	return MM
end
function dim2mix(X::Array{Array,1},Σ::Array{Array{Float64,2}})
	MM = MixtureModel(map(i->MvNormal(X[i],Σ[i]),1:length(X)))
	return MM
end
function dim2mix(X::Array{Array,1},Σ::Array{Float64,2},p::Array{Float64,1})
	MM = MixtureModel(map(x->MvNormal(x,Σ),X),p)
	return MM
end
function dim2mix(X::Array{Float64,2},Σ::Array{Float64,2})
	X = dim2adj(X)
	MM = MixtureModel(map(x->MvNormal(x,Σ),X))
	return MM
end
####Eval MixtureModel (obsolete with kde update?)
####Eval MixtureModel on an input-space (for plotting)
# function dim2pdf(MM::MixtureModel{Multivariate,Continuous,MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}},Ω::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
# 	OMEGA=0*Matrix{Float64}(undef,length(Ω),length(Ω))
# 	for j=1:length(Ω)
# 		for i=1:length(Ω)
# 				OMEGA[i,j]=pdf(MM,[Ω[i],Ω[j]])
# 		end
# 	end
# 	return ç
# end
# function dim2pdf(MM::MixtureModel{Multivariate,Continuous,MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}},Ω1::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Ω2::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}})
# 	OMEGA=0*Matrix{Float64}(undef,length(Ω1),length(Ω2))
# 	for j=1:length(Ω1)
# 		for i=1:length(Ω2)
# 				OMEGA[i,j]=pdf(MM,[Ω1[i],Ω2[j]])
# 		end
# 	end
# 	return OMEGA
# end
####Evaluate bounds of data and generate approriate range
####in bounds for package KernelDensity
function dim2bounds(X::Array{Float64,2})
	M = 0*Matrix{Float64}(undef,2,2)
	M[1,1] = minimum(X[:,1])
	M[2,1] = maximum(X[:,1])
	M[1,2] = minimum(X[:,2])
	M[2,2] = maximum(X[:,2])
	return M
end
function dim2range(strt::Float64,stp::Float64,len::Int64)
	Omega = range(strt,stop=stp,length=Int64(len))
	return Omega
end
####Generate appropriate range in bounds for data X_raw in dimension d
####The range has to be generated for each dimension
function dim2range(X_raw::Array{Float64,2},d::Int64)
	if size(X_raw,1)<size(X_raw,2)
		X_raw=ctrans(X_raw)
	end
	len = Int64(floor(sqrt(maximum(size(X_raw)))))
	M = dim2bounds(X_raw)
	Omega = range(M[1,d],stop=M[2,d],length=Int64(len))
	return Omega
end
function dim2range(X_raw::Array{Float64,2})
	if size(X_raw,1)<size(X_raw,2)
		X_raw=ctrans(X_raw)
	end
	len = Int64(floor(sqrt(maximum(size(X_raw)))))
	M = dim2bounds(X_raw)
	Omega1 = range(M[1,1],stop=M[2,1],length=Int64(len))
	Omega2 = range(M[1,2],stop=M[2,2],length=Int64(len))
	delta = [Omega1[end]-Omega1[1];Omega2[end]-Omega2[1]]
	return [Omega1,Omega2,delta]
end
##Likelihood & Eval
# function dim2l(MvN::MvNormal{Float64,PDMat{Float64,Array{Float64,2}},Array{Float64,1}},y::Array{Float64,1},p=1.0)
# 	val=pdf(p*MvN,y)
# end
# function dim2L(MM::MixtureModel{Multivariate,Continuous,MvNormal{Float64,PDMat{Float64,Array{Float64,2}},Array{Float64,1}}},Y)
# 	κ=length(MM.components)
# 	N=size(Y,1)
# 	𝕷=prod([sum([MM.prior.p[k]*pdf(MM.components[k],Y[i]) for k in 1:κ]) for i in 1:N])
# 	return 𝕷
# end
# function dim2logL(MM::MixtureModel{Multivariate,Continuous,MvNormal{Float64,PDMat{Float64,Array{Float64,2}},Array{Float64,1}}},Y)
# 	κ=length(MM.components)
# 	#print("dim2log κ=",κ,"\n")
# 	N=size(Y,1)
# 	#print("dim2log N=",N,"\n")
# 	𝕷=sum([sum([(-1/2.0)*(log(det(MM.components[k].Σ.mat))+(Y[i]-MM.components[k].μ)'*(MM.components[k].Σ.mat)^-1*(Y[i]-MM.components[k].μ)+2*log(2*π)) for k in 1:κ]) for i in 1:N])
# 	#print("dim2log𝕷 =",𝕷,"\n")
# 	return 𝕷
# end
#MLE geom dis given FrNm
function geomMLE(FrNm::Array{Int64,1})
	β=collect(1e-8:1e-2:1.0)
	βpost=[]
	#Likelihood for β parameter of Geometric Distribution, based on Framenmbers of Observations
	FrNm = FrRef(FrNm)
	for j=1:length(β)
		push!(βpost,sum([FrNm[i]*log(1-β[j])+log(β[j])] for i in 1:length(FrNm))[1])
	end
	#MLE for β Parameter
	βhat=β[collect(1:length(β))[βpost.==maximum(βpost)]][1]
	#Geometric Distribution with MLE for Parameter
	G=Geometric(βhat)
	return G
end
function geomMLE(FrNm::Int64)
	G=Geometric(0.999)
	return G
end
function normalMLE(FrNm::Array{Int64,1})
	μMLE=1/length(FrNm) * sum(FrNm)
	σ=1.0
	nMLE=Normal(μMLE,σ)
	# σMLE=sqrt(1/length(FrNm) * sum( (FrNm .- μMLE).^2 ) )
	# if sum(FrNm .== μMLE) == length(FrNm) σMLE = 1.0 end #avoid 0 value for σMLE (occurs for all FrNm==μMLE)
	#nMLE=Normal(μMLE,σMLE)
	return nMLE
end
##ref FrNm of obs to first frame
function FrRef(FrNm::Array{Int64,1})
	FrNm = (sort(FrNm) .- minimum(FrNm)) .+1
end
function FrRef(FrNm::Int64)
	return FrNm
end
#compute data based on cluster corresponding assignment in AcollElem
# function Feval(AcollElem::Array{Int64,1})#,𝔽sim::Array{Int64,1})
# 	uniqueLengthElem=length(unique(AcollElem))
# 	uniqueAcollElem=unique(AcollElem)
# 	#get framenumbers for each cluster
# 	𝔽evalElem=Array{Array,1}(undef,uniqueLengthElem)
# 	𝔽evalElem=map(i->𝔽sim[AcollElem .== unique(AcollElem)[i]],1:uniqueLengthElem)
# 	#eval indicies of observations (to evaluate locations in 𝕐sim)
# 	𝔽evalElemIdx=Array{Array,1}(undef,uniqueLengthElem)
# 	𝔽evalElemIdx=map(i->findall(x->x∈𝔽evalElem[i],𝔽sim),1:uniqueLengthElem)
# 	#get observations corresponding to clusters
# 	𝕐evalElem=Array{Array,1}(undef,uniqueLengthElem)
# 	𝕐evalElem=map(i->𝕐sim[𝔽evalElemIdx[i],:],1:uniqueLengthElem)
# 	#FrmNm f.e. Cl ID||Indices of obs f.e. Cl ID||Loc of obs f.e. Cl ID||number of unique Cl||ID of unique Cl
# 	return 𝔽evalElem,𝔽evalElemIdx,𝕐evalElem,uniqueLengthElem,uniqueAcollElem
# end
function Feval2(AcollElem::Array{Int64,1},SEG_FrNm)#,𝔽sim::Array{Int64,1})
	uniqueLengthElem=length(unique(AcollElem))
	uniqueAcollElem=unique(AcollElem)
	#get framenumbers for each cluster
	𝔽evalElem=Array{Array,1}(undef,uniqueLengthElem)
	𝔽evalElem=map(i->SEG_FrNm[AcollElem .== unique(AcollElem)[i]],1:uniqueLengthElem)
	#eval indicies of observations (to evaluate locations in 𝕐sim)
	𝔽evalElemIdx=Array{Array,1}(undef,uniqueLengthElem)
	𝔽evalElemIdx=map(i->findall(x->x∈𝔽evalElem[i],SEG_FrNm),1:uniqueLengthElem)
	#get observations corresponding to clusters
	𝕐evalElem=Array{Array,1}(undef,uniqueLengthElem)
	𝕐evalElem=map(i->𝕐sim[𝔽evalElemIdx[i],:],1:uniqueLengthElem)
	#FrmNm f.e. Cl ID||Indices of obs f.e. Cl ID||Loc of obs f.e. Cl ID||number of unique Cl||ID of unique Cl
	return 𝔽evalElem,𝔽evalElemIdx,𝕐evalElem,uniqueLengthElem,uniqueAcollElem
end
##	CORE FUNCTIONS
function mu_star(mu1::Array{Float64,1},mu2::Array{Float64,1},W::Array{Float64,1})
	mu_star = (W[1]*mu1+W[2]*mu2)/(W[1]+W[2])
	return mu_star
end
function s_star(mu1::Array{Float64,1},mu2::Array{Float64,1},S1::Array{Float64,2},S2::Array{Float64,2},W=NaN)#W::Array{Float64,1})
	if isnan(W[1]) W=[1.,1.]; print("!!!\nWEIGHT ISNAN\n!!!\n") end
	Σ_star = (W[1]*S1 + W[2]*S2)/(W[1]+W[2]) + (W[1]*W[2]*(mu2-mu1)*(mu2-mu1)')/((W[1]+W[2])^2)
	return Σ_star
end
function dim2KLD(mu1::Array{Float64,1},mu2::Array{Float64,1},S1::Array{Float64,2},S2::Array{Float64,2},W=NaN)#W::Array{Float64,1})
	S=[S1,S2]
	if isnan(W[1]) W=[1.,1.]; print("!!!\nWEIGHT ISNAN\n!!!\n") end
	Σ_star = s_star(mu1,mu2,S[1],S[2],[W[1],W[2]])
	𝛅=Int64(2)
	KLD = (sum([W[i]/2.0 * (tr(Σ_star^-1 * S[i]) + log(det(Σ_star * S[i]^-1))-𝛅)] for i in 1:2)[1] + (W[1]*W[2]/(2.0*(W[1]+W[2])^2) * (mu1-mu2)'*Σ_star^-1*(mu1-mu2)))
	if isnan(KLD) KLD=0.0; print("!!!\nKLD ISNAN\n!!!\n") end
	return KLD
end
#Extrema of KLD matrix
function dim2KLDex(KLDfull::Array{Float64,2})
	𝛕=1e-5
	EX = 0*Array{Int64,2}(undef,3,2)
	K = maximum(size(KLDfull))
	#column1 minimum KLD indices [vectorindex, matrixindices[j,k]]
	EX[2,1] = findall(KLDfull .== minimum(KLDfull[findall((KLDfull .> 𝛕))]))[1][1] #column index
	EX[3,1] = findall(KLDfull .== minimum(KLDfull[findall((KLDfull .> 𝛕))]))[1][2] #row index
	EX[1,1] = Int64(floor((EX[3,1]-1)*K+EX[2,1])) #single indexing
	EXmin = KLDfull[EX[1,1]]
	#column2 maximum KLD indices [vectorindex, matrixindices[j,k]]
	EX[2,2] = findall(KLDfull .== maximum(KLDfull[findall((KLDfull .> 𝛕))]))[1][1]
	EX[3,2] = findall(KLDfull .== maximum(KLDfull[findall((KLDfull .> 𝛕))]))[1][2]#Int64(EX[1,2]-K*(EX[3,2]-1))
	EX[1,2] = Int64(floor((EX[3,2]-1)*K+EX[2,2]))
	EXmean=mean(collect(KLDfull[:]))
	#print("Idx min: ",EX[1,1],"= ","[",EX[2,1],",",EX[3,1],"]","min value: ", EXmin,"(mean value: ",EXmean,")\n")
	return EX
end
##Weighted Kullback Leibler Divergence
function dim2WKLDfull(MU::Array{Array,1},S::Array{Array{Float64,2},1},W::Array{Float64,1},SEG_FrNm,A=NaN,Atmp=NaN,WKLD_modINPUT=1)
	#	WKLD_modINPUT=1 	...	wKLD
	#	WKLD_modINPUT=2	...	KLD
	#	WKLD_modINPUT=3	...	experimental
	#entries are merged during each collapse, therefore K represents the current number of clusters
	K=maximum(size(MU))
	#if no A or Atmp is submitted, it's set as NaN (though A and Atmp are by default submitted during dim2WKLDfull computation)
	if isnan(A[1]) A=collect(1:K) end
	if isnan(Atmp[1]) Atmp=collect(1:K) end
	# MU=dim2adj(dim2adj(MU))
	global MU_star = map(i->MU[i],1:length(MU)) #(clustered-)data based on 𝕐sim
	global A_star = map(i->A[i],1:length(A)) #cluster ID

	W_star = map(i->W[i],1:length(W)) #weights
	A_coll = map(i->Atmp[i],1:length(Atmp)) #A_coll stores the cluster IDs which are updated for every collapse
	𝛅 = Int64(3) #dimension
	# S_star=S
	S_star= repeat([1.0*eye(𝛅)],length(S)) #COV Matrices
	K = maximum(size(MU_star))
	#print("K=",K,"\n")
	#initialize "distance" matrix
	KLD = 1e-15*ones(K,K)
	FrNm = Feval2(A_coll, SEG_FrNm)[1]
	# FrNm = Feval2(A_coll,SEG_FrNm)[1] #evaluate frames f.e. clust ID in A_coll
	#compute normal distribution with μMLE and σMLE for each current cluster
	#NormDistClustID = map(i->normalMLE(FrNm[i]),1:length(FrNm))
	GeomDistClustID = map(i->geomMLE(FrRef(FrNm[i])),1:length(FrNm))
	weightParaKLD = 1e-15*ones(K,K)
	weightParaKLD_norm = 1e-15*ones(K,K)
	if WKLD_modINPUT==1	# WEIGHTED clustering
		#KL distances
		#compute KL-distances for each and every pair of points
		for k=1:K
			for j=1:k
				if j==k
					weightParaKLD[j,k] = (1/exp(median(pdf.(GeomDistClustID[j],FrRef(FrNm[k])))))
				else
					# weightParaKLD[k,j] = (1/exp(exp(sum(pdf.(GeomDistClustID[j],FrNm[k])))))	#k in j geclustert
					# # print("pdf.(GeomDistClustID[$(j)],FrNm[$(k)])=$(pdf.(GeomDistClustID[j],FrNm[k]))\n")
					# # print("K= $(K), weightParaKLD[$(j),$(k)]=$(weightParaKLD[j,k])\n")
					# weightParaKLD[j,k] = (1/exp(exp(sum(pdf.(GeomDistClustID[k],FrNm[j])))))
					# # print("K= $(K), weightParaKLD[$(k),$(j)]=$(weightParaKLD[k,j])\n")
					weightParaKLD[j,k] = (1/exp(median(pdf.(GeomDistClustID[j],FrRef(FrNm[k])))))
					weightParaKLD[k,j] = (1/exp(median(pdf.(GeomDistClustID[k],FrRef(FrNm[j])))))
				end
			end
			# weightParaKLD[k,:]=weightParaKLD[k,:].*(1/sum(weightParaKLD[k,:]))
		end
		for k=1:K
			for j=1:K
				if j==k
					KLD[j,k]=1e5
				else
					# KLD[j,k] = dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0]) 	# no weight
					KLD[j,k] = (1/exp(sum(pdf.(GeomDistClustID[k],FrNm[j])))) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])	# weight

					# KLD[j,k] = (weightParaKLD[j,k]/exp(sum(pdf.(GeomDistClustID[k],FrNm[j])))) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])	# experimental
					# KLD[j,k] = dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[weightParaKLD[k,j]*(1/sum(weightParaKLD[k,:])),weightParaKLD[j,k]*(1/sum(weightParaKLD[:,k]))])	# experimental
					# KLD[j,k] = weightParaKLD[j,k]*(1/sum(weightParaKLD[:,k])) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
					# weightParaKLD_norm[j,k] = weightParaKLD[j,k]*(1/sum(weightParaKLD[:,k]))
					# # print("weightParaKLD[$k,$j]*(1/sum(weightParaKLD[$k,:]))= $(weightParaKLD[k,j]*(1/sum(weightParaKLD[k,:])))\n")
					# # print("weightParaKLD[$j,$k]*(1/sum(weightParaKLD[:,$k]))= $(weightParaKLD[j,k]*(1/sum(weightParaKLD[:,k])))\n")
				end
			end
		end

		#dim2KLDex computes extrema and corresponding indices in the distance matrix
		KLD[map(i->isnan(KLD[i]),1:length(KLD[:]))].= maximum(KLD[map(i->!isnan(KLD[i]),1:length(KLD[:]))])
		EX = dim2KLDex(KLD)
		#min1,min2 is the pair of values with minimal KLD and get merged
		min1=EX[2,1]
		min2=EX[3,1]
		global MU_star[EX[2,1]] = mu_star(MU_star[EX[2,1]],MU_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
		S_star[EX[2,1]] = s_star(MU_star[EX[2,1]],MU_star[EX[3,1]],S_star[EX[2,1]],S_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
		W_star[EX[2,1]] += W_star[EX[3,1]]
		S_star[EX[3,1]] = NaN.*S_star[EX[3,1]]
		global MU_star[EX[3,1]] =NaN.*MU_star[EX[3,1]]
		W_star[EX[3,1]] = NaN.*W_star[EX[3,1]]
		A_coll[A_coll .== A_star[EX[3,1]]].= A_star[EX[2,1]]
		#global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
		global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
		global A_star=A_star[map(i->A_star[i]!=0,1:length(A_star))]
		S_star=S_star[map(i->!isnan(S_star[i][1]),1:length(S_star))]
		global MU_star=MU_star[map(i->!isnan(MU_star[i][1]),1:length(MU_star))]
		W_star=W_star[map(i->!isnan(W_star[i]),1:length(W_star))]

		S_star=repeat([1.0*eye(𝛅)],length(S_star))
		KLD[EX[1,1]]=0.0
		#print("END OF while loop dim2WKLDfull\n")
		#EDIT 01062021
		# MU_star=MU_star[map(i->!isnan(MU_star[i][1]),1:length(MU_star))]
		# W_star=W_star[map(i->!isnan(W_star[i]),1:length(W_star))]
		# A_star=A_star[map(i->A_star[i]!=0,1:length(A_star))]
		return MU_star,W_star,S_star,A_star,A_coll,KLD,weightParaKLD,GeomDistClustID
##
	elseif WKLD_modINPUT==2	# UNweighted clustering
		for k=1:K
			for j=1:K
				if j==k
					KLD[j,k] = 1e5
				else
					KLD[j,k] = dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
				end
			end
		end
		#dim2KLDex computes extrema and corresponding indices in the distance matrix
		KLD[map(i->isnan(KLD[i]),1:length(KLD[:]))].= maximum(KLD[map(i->!isnan(KLD[i]),1:length(KLD[:]))])
		EX = dim2KLDex(KLD)
		#min1,min2 is the pair of values with minimal KLD and get merged
		min1=EX[2,1]
		min2=EX[3,1]
		global MU_star[EX[2,1]] = mu_star(MU_star[EX[2,1]],MU_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
		S_star[EX[2,1]] = s_star(MU_star[EX[2,1]],MU_star[EX[3,1]],S_star[EX[2,1]],S_star[EX[3,1]],[W_star[EX[2,1]],W_star[EX[3,1]]])
		W_star[EX[2,1]] += W_star[EX[3,1]]
		S_star[EX[3,1]] = NaN.*S_star[EX[3,1]]
		global MU_star[EX[3,1]] =NaN.*MU_star[EX[3,1]]
		W_star[EX[3,1]] = NaN.*W_star[EX[3,1]]
		A_coll[A_coll .== A_star[EX[3,1]]].= A_star[EX[2,1]]
		#global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
		global A_star[EX[3,1]] = 0*A_star[EX[3,1]]
		global A_star=A_star[map(i->A_star[i]!=0,1:length(A_star))]
		S_star=S_star[map(i->!isnan(S_star[i][1]),1:length(S_star))]
		global MU_star=MU_star[map(i->!isnan(MU_star[i][1]),1:length(MU_star))]
		W_star=W_star[map(i->!isnan(W_star[i]),1:length(W_star))]
		#S_star=repeat([1.0*eye(𝛅)],length(S_star))
		# KLD[EX[1,1]]=0.0
		#print("END OF while loop dim2WKLDfull\n")
		#EDIT 01062021
		MU_star=MU_star[map(i->!isnan(MU_star[i][1]),1:length(MU_star))]
		W_star=W_star[map(i->!isnan(W_star[i]),1:length(W_star))]
		A_star=A_star[map(i->A_star[i]!=0,1:length(A_star))]
		return MU_star,W_star,S_star,A_star,A_coll,KLD,weightParaKLD,GeomDistClustID
##
	elseif WKLD_modINPUT==3
		for k=1:K
			for j=1:K
				if j==k
					KLD[j,k]=1e5
				else
					KLD[j,k] = (1/exp(sum(pdf.(GeomDistClustID[k],FrNm[j])))) * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
				end
			end
		end
		return KLD
	end
end
# ##dim2KLD
# function dim2KLDfull_UNWEIGHTED(MU::Array{Array,1},S::Array{Array{Float64,2},1})
# 	#entries are merged during each collapse, therefore K represents the current number of clusters
# 	K=maximum(size(MU))
# 	#𝛅 = Int64(3) #dimension
# 	MU_star = map(i->MU[i],1:length(MU)) #(clustered-)data based on 𝕐sim
# 	S_star=S
# 	#S_star= repeat([1.0*eye(𝛅)],length(S)) #COV Matrices
# 	#print("K=",K,"\n")
# 	#initialize "distance" matrix
# 	KLD = 1e-15*ones(K,K)
# 	for k=1:K
# 		for j=1:K
# 			if j==k
# 				KLD[j,k] = 0
# 			else
# 				KLD[j,k] = abs(dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0]))
# 			end
# 		end
# 	end
# 	return KLD
# end
# ##
# function dim2KLDfull_WEIGHTED(MU::Array{Array,1},S::Array{Array{Float64,2},1},SEG_FrNm)
# 	#this version of the function
# 	MU_star = map(i->MU[i],1:length(MU)) #(clustered-)data based on 𝕐sim
# 	S_star=S
# 	K = maximum(size(MU_star))
# 	KLD = 1e-15*ones(K,K)
# 	#GeomDistClustID = map(i->geomMLE(FrRef(SEG_FrNm[i])),1:length(SEG_FrNm))
# 	GeomDistClustID = geomMLE(FrRef(SEG_FrNm))
# 	#GeomDistClustID = map(i->geomMLE(FrRef(tmp3.Distance[8].frame_SR[i])),1:length(tmp3.Distance[8].frame_SR))
# 	weightParaKLD = ones(K,K)
# 	for k=1:K
# 		for j=1:K
# 			if j==k
# 				KLD[j,k] = 0
# 			else
# 				weightParaKLD[j,k] = (1/exp(sum(pdf.(GeomDistClustID,SEG_FrNm))))
# 				KLD[j,k] = weightParaKLD[j,k] * dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0])
# 				# weightParaKLD[j,k] = 1/(exp(abs(cdf(GeomDistClustID,FrRef(SEG_FrNm)[j])-cdf(GeomDistClustID,FrRef(SEG_FrNm)[k]))))
# 				# KLD[j,k] = abs(weightParaKLD[j,k]*dim2KLD(MU_star[k],MU_star[j],S_star[k],S_star[j],[1.0,1.0]))
# 			end
# 		end
# 	end
# 	return KLD,weightParaKLD
# end
##
############################################################################
#	global variables can get changed inside a loop, but keep it's previously
#	assigned global value.
#	symbols: 𝕐 = \bbY
############################################################################

export dim
export eye
export ctrans
export dim2clean
export dim2crop
export dim2adj
export dim2mix
export dim2bounds
export dim2range
export geomMLE
export normalMLE
export FrRef
export Feval
export mu_star
export s_star
export dim2KLD
export dim2KLDex
export dim2WKLDfull
#export dim2KLDfull_UNWEIGHTED
#export dim2KLDfull_WEIGHTED
export heatSR

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
#=
 Example 1: Track the pathways of a user-defined water mass.
 Steps: (a) define the water mass 1). by a rectangular surface patch
            dyed with passive tracer concentration of 1,
            or 2. load a pre-defined surface patch in d_all.mat.
        (b) propagate the dye with the matrix A, with the result
            being the fraction of water originating from the
            surface region.
 See Section 2b of Gebbie & Huybers 2010, esp. eqs. (15)-(17).
=#
using Revise
using TMI, PyPlot, PyCall

TMIversion = "modern_90x45x33_GH10_GH12"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)

#- define the surface patch by the bounding latitude and longitude.
latbox = [50,60]; # 50 N -> 60 N, for example.

# mutable due to wraparound: don't use an immutable tuple
lonbox = [-50,0]; # 50 W -> prime meridian

# do numerical analysis
c = trackpathways(Alu,latbox,lonbox,γ);

# do plotting (could be a function)
plotextent(latbox, lonbox)

# plot a section at 330 east longitude (i.e., 30 west)
lon_section = 330; # only works if exact
lims = 0:5:100
tlabel = "water-mass fraction [%]"
sectionplot(100c,lon_section,lims,titlelabel = tlabel)

ext: jl


text:
if VERSION < v"1.3-"
    using Base.PCRE
    Base.:(^)(r::Regex, i::Integer) = Regex(string("(?:", r.pattern, "){$i}"), r.compile_options, r.match_options)
    function Base.:(*)(r1::Union{Regex,AbstractString,AbstractChar}, rs::Union{Regex,AbstractString,AbstractChar}...)
        mask = PCRE.CASELESS | PCRE.MULTILINE | PCRE.DOTALL | PCRE.EXTENDED # imsx
        match_opts   = nothing # all args must agree on this
        compile_opts = nothing # all args must agree on this
        shared = mask
        for r in (r1, rs...)
            r isa Regex || continue
            if match_opts === nothing
                match_opts = r.match_options
                compile_opts = r.compile_options & ~mask
            else
                r.match_options == match_opts &&
                    r.compile_options & ~mask == compile_opts ||
                    throw(ArgumentError("cannot multiply regexes: incompatible options"))
            end
            shared &= r.compile_options
        end
        unshared = mask & ~shared
        Regex(string(wrap_string(r1, unshared), wrap_string.(rs, Ref(unshared))...), compile_opts | shared, match_opts)
    end

    Base.:(*)(r::Regex) = r # avoids wrapping r in a useless subpattern
    wrap_string(r::Regex, unshared::UInt32) = string("(?", regex_opts_str(r.compile_options & unshared), ':', r.pattern, ')')
# if s contains raw"\E", split '\' and 'E' within two distinct \Q...\E groups:
    wrap_string(s::AbstractString, ::UInt32) =  string("\\Q", replace(s, raw"\E" => raw"\\E\QE"), "\\E")
    wrap_string(s::AbstractChar, ::UInt32) = string("\\Q", s, "\\E")

    regex_opts_str(opts) = (isassigned(_regex_opts_str) ? _regex_opts_str[] : init_regex())[opts]

# UInt32 to String mapping for some compile options
    const _regex_opts_str = Ref{Base.ImmutableDict{UInt32,String}}()

    init_regex() = _regex_opts_str[] = foldl(0:15, init=Base.ImmutableDict{UInt32,String}()) do d, o
        opt = UInt32(0)
        str = ""
        if o & 1 != 0
            opt |= PCRE.CASELESS
            str *= 'i'
        end
        if o & 2 != 0
            opt |= PCRE.MULTILINE
            str *= 'm'
        end
        if o & 4 != 0
            opt |= PCRE.DOTALL
            str *= 's'
        end
        if o & 8 != 0
            opt |= PCRE.EXTENDED
            str *= 'x'
        end
        Base.ImmutableDict(d, opt => str)
    end
end


text:
## ----------------------------------------------------------------------------------
_unbig(x) = round(Float64(x); digits = 15)
function trunc_sample(μ::Number, σ::Number, lb::Number, ub::Number; xbins = 1000)

    local_max = clamp(μ, lb, ub)
    max_range = max(abs(local_max - lb), abs(local_max - ub))
    log_range = range(log10(max_range), log10(σ/xbins), length = div(xbins, 2))
    xs = [
        local_max .- 10.0.^(log_range);
        local_max;
        local_max .+ 10.0.^(reverse(log_range))
    ]
    Txs = [big(x) for x in xs if lb <= x <= ub]
    Tpdf = ϕ.(Txs, big(μ), big(σ))
    Z = sum(@view(Tpdf[1:end - 1]) .* diff(Txs))
    Tpdf .= Tpdf ./ Z
    _unbig.(Txs), _unbig.(Tpdf)
end


text:
#under development

using Flux
using Flux: @functor

abstract type AbstractTransformer end

struct PwFFN
    din::Dense
    dout::Dense
end

struct Transformer <: AbstractTransformer
    mh::MultiheadAttention
    mhn::LayerNorm
    pw::PwFFN
    pwn::LayerNorm
    drop::Dropout
    proj::Dense
end

@functor Transformer

"""
    Transformer(size::Int, head::Int, ps::Int;
                future::Bool = true, act = relu, pdrop = 0.1)
    Transformer(size::Int, head::Int, hs::Int, ps::Int;
                future::Bool = true, act = relu, pdrop = 0.1)

Transformer layer.

`size` is the input size. if `hs` is not specify, use `div(size, head)` as the hidden size of multi-head attention.
`ps` is the hidden size & `act` is the activation function of the positionwise feedforward layer.
When `future` is `false`, the k-th token can't see the j-th tokens where j > k. `pdrop` is the dropout rate.
"""
function Transformer(size::Int, head::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1)
    rem(size, head) != 0 && error("size not divisible by head")
    Transformer(size, head, div(size, head), ps;future=future, act=act, pdrop=pdrop)
end

Transformer(size::Int, head::Int, hs::Int, ps::Int; future::Bool = true, act = relu, pdrop = 0.1) = Transformer(
    MultiheadAttention(head, size, hs, size; future=future, pdrop=pdrop),
    LayerNorm(size),
    PwFFN(size, ps, act),
    LayerNorm(size),
    Dropout(pdrop),
    Dense(hs,hs)
)

function (t::Transformer)(x::AbstractArray{T, N}, mask=nothing) where {T, N}
    h=t.mh(x,x,x,mask=mask)
    h=t.pwn(h .+ t.proj(h,h))
    h=t.pwn(h .+ t.pw(h))
end


text:
N = 800000 # number of observations
M = 40000 # number of workers
O = 5000 # number of firms
id1 = rand(1:M, N)
id2 = [rand(max(1, div(x, 8)-10):min(O, div(x, 8)+10)) for x in id1]
x1 = 5 * cos.(id1) + 5 * sin.(id2) + randn(N)
x2 =  cos.(id1) +  sin.(id2) + randn(N)
y= 3 .* x1 .+ 5 .* x2 .+ cos.(id1) .+ cos.(id2).^2 .+ randn(N)
df = DataFrame(id1 = id1, id2 = id2, x1 = x1, x2 = x2, y = y)
@time reg(df, @formula(y ~ x1 + x2 + fe(id1) + fe(id2)))


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ 245de79e-260c-11ec-082f-bdef59a44c40
using ControlSystems

# ╔═╡ f5e08ada-5f12-45a2-ba56-62af6bbe0a0d
begin
	J = 2.0
	b = 0.04
	K = 1.0
	R = 0.08
	L = 1e-4

	# dc motor model
	s = tf("s")
	P = K/(s*((J*s + b)*(L*s + R) + K^2))

	# Create an array of closed loop systems for different values of kp
	Tcl = TransferFunction[kp*P/(1 + kp * P) for kp = [1, 5, 15]]

	# Closed loop step response
	stepplot(Tcl, label= ["kp = 1" "kp = 5" "kp = 15"])
end

# ╔═╡ 3f85cf01-bcd7-4395-84bf-2c5eb597f497

# ╔═╡ Cell order:
# ╠═245de79e-260c-11ec-082f-bdef59a44c40
# ╠═f5e08ada-5f12-45a2-ba56-62af6bbe0a0d
# ╠═3f85cf01-bcd7-4395-84bf-2c5eb597f497


text:
using LinearAlgebra
using AdFem
using DelimitedFiles

A = readdlm("fenics/A.txt")
A1 = readdlm("fenics/B.txt")

A = [A A1]
mesh = Mesh(8, 8, 1. /8, degree=1)
B = compute_interaction_matrix(mesh)

@show norm(A - B)

mesh = Mesh(8, 8, 1. / 8, degree=2)

A = readdlm("fenics/A2.txt")
A1 = readdlm("fenics/B2.txt")
E = Int64.(readdlm("fenics/edges.txt"))
Edof = get_edge_dof(E, mesh)
A = [A A1]
B = Array(compute_interaction_matrix(mesh))

B[:, mesh.nnode + 1:mesh.ndof] = B[:, mesh.nnode + 1:mesh.ndof][:, Edof]
B[:, mesh.ndof + mesh.nnode + 1:2mesh.ndof] = B[:, mesh.ndof + mesh.nnode + 1:2mesh.ndof][:, Edof]
@show norm(A - B)


text:
using Test
#using AbstractPlotting
import Calculus
import ForwardDiff
#using GLMakie; GLMakie.activate!()
#using MakieLayout
import NonlinearEquations
import PyPlot
import Random
import SparseArrays
using StaticArrays

ForwardDiff_gradient(x...) = ForwardDiff.gradient(x...)

function kr(psi, alpha, N)
	if psi < 0
		m = (N - 1) / N
		denom = 1 + abs(alpha * psi) ^ N
		numer = 1 - abs(alpha * psi) ^ (N - 1) * denom ^ (-m)
		return numer ^ 2 / denom ^ (m / 2)
	else
		return one(psi)
	end
end

kr(x::AbstractArray) = kr(x[1], x[2], x[3])

function Calculus.differentiate(x::Calculus.SymbolParameter{:kr}, args, wrt)
	chain_part = map(x->Calculus.simplify(Calculus.differentiate(x, wrt)), args)
	if chain_part == [0, 0, 0]
		return :(0)
	else
		return :(sum(ForwardDiff_gradient(kr, SA[$(args...)]) .* SA[$(chain_part...)]))
	end
end

function hm(x, y)
	return 2 / (1 / x + 1 / y)
end

hm(x::AbstractArray) = hm(x[1], x[2])

function Calculus.differentiate(x::Calculus.SymbolParameter{:hm}, args, wrt)
	chain_part = map(x->Calculus.simplify(Calculus.differentiate(x, wrt)), args)
	if chain_part == [0, 0]
		return :(0)
	else
		return :(sum(ForwardDiff_gradient(hm, SA[$(args...)]) .* SA[$(chain_part...)]))
	end
end

function Calculus.differentiate(x::Calculus.SymbolParameter{:abs}, args, wrt)
	if length(args) > 1
		error("Too many arguments passed to abs()")
	end
	arg = args[1]
	return :(ifelse($arg > 0, 1, -1) * $(Calculus.differentiate(arg, wrt)))
end

function effective_saturation(alpha::Number, psi::Number, N::Number)
	m = (N - 1) / N
	if psi < 0
		return (1 + abs(alpha * psi) ^ N) ^ (-m)
	else
		return 1
	end
end

function saturation(psi::Array, coords::Array)
	Se = similar(psi)
	for i = 1:length(psi)
		if inclay(coords[1, i], coords[2, i])
			_, alpha, N, residual_saturation = params[:clay]
		else
			_, alpha, N, residual_saturation = params[:claysilt]
		end
		Se[i] = residual_saturation + effective_saturation(alpha, psi[i], N) * (1 - residual_saturation)
	end
	return Se
end

function callback(psi, residuals, J, i, saturation=saturation(psi, coords))
	if mod(i, 10) == 0
		fig, axs = PyPlot.subplots(2, 2, figsize=(16, 9))
		psi = reshape(psi, ns[2], ns[1])
		saturation = reshape(saturation, ns[2], ns[1])
		residuals = reshape(residuals, ns[2], ns[1])
		img = axs[1].imshow(map(x->x > 0 ? x : 0, psi), origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[1])
		img = axs[2].imshow(saturation, origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[2])
		goodindices = filter(i->coords[2, i] == maxs[2], 1:size(coords, 2))
		img = axs[3].imshow(map(x->x < 0 ? x : 0, psi), origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[3])
		img = axs[4].imshow(residuals, origin="lower", interpolation="nearest", extent=[mins[1], maxs[1], mins[2], maxs[2]])
		fig.colorbar(img, ax=axs[4])
		display(fig)
		println()
		println()
		PyPlot.close(fig)
	end
end

#=
function animatetransient(soln)
	hend = soln(tspan[end])
	scene, layout = layoutscene(resolution=(600, 600))
	t_observable = Node(0.0)
	h_observable = lift(t->reshape(soln(t), ns[2], ns[1])', t_observable)
	saturation_observable = lift(t->reshape(saturation(soln(t), coords), ns[2], ns[1])', t_observable)
	ax1 = layout[1, 1] = LAxis(scene; title="Pressure Head")
	ax2 = layout[2, 1] = LAxis(scene; title="Saturation")
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	heatmap!(ax1, xs, ys, h_observable; colorrange=extrema(hend))
	heatmap!(ax2, xs, ys, saturation_observable, colorrange=(0, 1))
	for ax in [ax1, ax2]
		xlims!(ax, mins[1], maxs[1])
		ylims!(ax, mins[2], maxs[2])
	end
	layout[1, 2] = LColorbar(scene; limits=extrema(hend), width=30)
	layout[2, 2] = LColorbar(scene; limits=(0, 1), width=30)
	#timeplot = layout[2, :] = LAxis(scene; title="time", xgridvisible=false, ygridvisible=false, yticklabelsvisible=false, yticksvisible=false)
	#timeplot.height=20
	#lines!(timeplot, t_observable, [0, 1], linewidth=5)
	#ylims!(timeplot, 0, 1)
	#xlims!(timeplot, 0, tspan[end])
	#display(scene)
	record(scene, "animation_$(ns[1]).mp4", 0:0.2:62; framerate=60) do t
		t_observable[] = t
	end
end

function plottransient(soln)
	hend = soln(tspan[end])
	scene, layout = layoutscene(resolution=(1200, 600))
	t_slider = layout[3, 1:3] = LSlider(scene; range=range(0, tspan[end]; length=500), startvalue=0)
	t_observable = lift(t->[t, t], t_slider.attributes[:value])
	h_observable = lift(t->reshape(soln(t), ns[2], ns[1])', t_slider.attributes[:value])
	saturation_observable = lift(t->reshape(saturation(soln(t), coords), ns[2], ns[1])', t_slider.attributes[:value])
	ax1 = layout[1, 1] = LAxis(scene; title="Pressure Head")
	ax2 = layout[1, 3] = LAxis(scene; title="Saturation")
	xs = range(mins[1]; stop=maxs[1], length=ns[1])
	ys = range(mins[2]; stop=maxs[2], length=ns[2])
	heatmap!(ax1, xs, ys, h_observable; colorrange=extrema(hend))
	heatmap!(ax2, xs, ys, saturation_observable, colorrange=(0, 1))
	for ax in [ax1, ax2]
		xlims!(ax, mins[1], maxs[1])
		ylims!(ax, mins[2], maxs[2])
	end
	layout[1, 2] = LColorbar(scene; limits=extrema(hend), width=30)
	layout[1, 4] = LColorbar(scene; limits=(0, 1), width=30)
	timeplot = layout[2, :] = LAxis(scene; title="time", xgridvisible=false, ygridvisible=false, yticklabelsvisible=false, yticksvisible=false)
	timeplot.height=20
	lines!(timeplot, t_observable, [0, 1], linewidth=5)
	ylims!(timeplot, 0, 1)
	xlims!(timeplot, 0, tspan[end])
	display(scene)
end
=#

function i2i1i2(i)
	i2 = 1 + mod(i - 1, ns[2])
	i1 = 1 + div(i - i2, ns[2])
	return i1, i2
end

function updateneighborarrays!(a1, a2, pair, value)
	n1, n2 = pair
	if i2i1i2(n1)[1] == i2i1i2(n2)[1]#the neighbors have the same x coordinate, update a2
		i1, i2_1 = i2i1i2(n1)
		i2_2 = i2i1i2(n2)[2]
		@assert i2_1 < i2_2
		a2[i2_1, i1] = value
	else#the neighbors have the same y coordinate, update a1
		i1_1, i2 = i2i1i2(n1)
		i1_2 = i2i1i2(n2)[1]
		@assert i1_1 < i1_2
		a1[i2, i1_1] = value
	end
end

function breakupgradient(grad)
	grad1 = zeros(ns[2], ns[1] - 1)
	grad2 = zeros(ns[2] - 1, ns[1])
	for (v, pair) in zip(grad, neighbors)
		updateneighborarrays!(grad1, grad2, pair, v)
	end
	return grad1, grad2
end

function plotgradient(grad, obsnode; vmin=minimum(grad), vmax=maximum(grad))
	grad1, grad2 = breakupgradient(grad)
	fig, axs = PyPlot.subplots(1, 2, figsize=(16, 9))
	img = axs[1].imshow(grad1, extent=[mins[1], maxs[1], mins[2], maxs[2]], origin="lower", interpolation="nearest", vmin=vmin, vmax=vmax)
	fig.colorbar(img, ax=axs[1])
	img = axs[2].imshow(grad2, extent=[mins[1], maxs[1], mins[2], maxs[2]], origin="lower", interpolation="nearest", vmin=vmin, vmax=vmax)
	fig.colorbar(img, ax=axs[2])
	display(fig)
	println()
	println()
	PyPlot.close(fig)
end


text:
using Roots

# Unoptimized reference implementation
# NB: ϕ is scaled by √2π compared to standard definition
function SD_reference(α, X::AbstractArray{T}) where T
	# (n(n-1))⁻¹α⁻⁵∑ᵢ∑ⱼϕⁱᵛ(α⁻¹(Xᵢ-Xⱼ))
	n = length(X)
	C = 1/(n*(n-1)*α^5)
	Y = X./α
	s = zero(T)
	@inbounds for i=1:n
		for j=i+1:n
			s += 2*ϕ4(Y[i]-Y[j])
		end
	end
	s += n*ϕ4(0)
	C * s
end

# NB: ϕ is scaled by √2π compared to standard definition
function TD_reference(b, X::AbstractArray{T}) where T
	# -(n(n-1))⁻¹b⁻⁷∑ᵢ∑ⱼϕᵛⁱ(b⁻¹(Xᵢ-Xⱼ))
	n = length(X)
	C = -1/(n*(n-1)*b^7)
	Y = X./b
	s = zero(T)
	@inbounds for i=1:n
		for j=i+1:n
			s += 2*ϕ6(Y[i]-Y[j])
		end
	end
	s += n*ϕ6(0)
	C * s
end

function objective_reference(h, X, α2Constant)
	# [ R(K) / (nσ⁴ₖSD(α₂(h))) ]^(1/5) - h = 0
	# R(K) = 1/(2√π)
	# σ⁴ₖ = 1
	n = length(X)
	α2 = α2Constant*h^(5/7) # 1.357[SD(a)/TD(b)]^(1/7) * h^(5/7)
	1/(√2*n*h^5) - SD_reference(α2,X)
end

function _bwsj_reference(X; rtol=0.1)
	n = length(X)
	q25,q75 = quantile(X,(0.25,0.75))
	λ = min(q75-q25, std(X)*1.349) # Givivng the same scale estimate as e.g. the KernSmooth R package.
	a = 0.920λ*n^(-1/7)
	b = 0.912λ*n^(-1/9)
	α2Constant = 1.357*(SD_reference(a,X)/TD_reference(b,X))^(1/7)

	# Figure out lower and upper bounds for search.
	# The Bandwidth scales roughly with n^(-1/5). TODO: Use this fact!

	lower = 0.1*λ
	for i=1:11
		i==11 && error("Failed to find lower bound")
		objective_reference(lower,X,α2Constant)>0 && break
		lower/=10
	end
	upper = 1*λ
	for i=1:11
		i==11 && error("Failed to find upper bound")
		objective_reference(upper,X,α2Constant)<0 && break
		upper*=10
	end

	find_zero(h->objective_reference(h,X,α2Constant), (lower,upper), Bisection(), xrtol=rtol)
end

bwsj_reference(X; kwargs...) = _bwsj_reference(issorted(X) ? X : sort(X); kwargs...)

function _smooth_reference(x,y,c,xeval,w)
	w .= -c.*(x.-xeval).^2
	w .= exp.(w .- maximum(w)) # multiply all weights with a factor for improved precision
	y'w / sum(w)
end

function smooth_reference(x::AbstractVector{<:T}, y::AbstractVector{<:S}, bandwidth, xeval) where {T<:Real,S<:Real}
	@assert length(x)==length(y)
	c = 1 ./(2 .* bandwidth.^2)
	w = similar(x, promote_type(T,eltype(bandwidth),eltype(xeval)))
	_smooth_reference.(Ref(x),Ref(y),c,xeval,Ref(w))
end

function _density_reference(x,bandwidth,c,xeval,w)
	w .= exp.(-c.*(x.-xeval).^2)
	sum(w) / (sqrt(2π)*bandwidth*length(x))
end

function density_reference(x::AbstractVector{<:T}, bandwidth, xeval) where {T<:Real}
	c = 1 ./(2 .* bandwidth.^2)
	w = similar(x, promote_type(T,eltype(bandwidth),eltype(xeval)))
	_density_reference.(Ref(x),bandwidth,c,xeval,Ref(w))
end


text:
#=
using MATLAB
function CorNewton3(G)
    @assert issymmetric(G)
    n = size(G, 1)
    mat"
    [$X,$y] = CorNewton3($(Array(G)),ones($n,1),1:$n,1:$n,0.0);
    "
    return Symmetric(X), y
end
=#

function fronorm(A, work)
    lda, n = size(A)
    ccall(
        (:dlansy_64_, "libopenblas64_"),
        Cdouble,
        (
            Ref{UInt8},
            Ref{UInt8},
            Ref{Int},
            Ptr{Float64},
            Ref{Int},
            Ptr{Float64},
        ),
        'F',
        A.uplo,
        n,
        A.data,
        lda,
        work,
    )
end

function symmetrize!(A)
    n = size(A, 1)
    for j = 1:n
        for i = 1:j-1
            @inbounds A.data[j, i] = A.data[i, j]
        end
    end
    return nothing
end

function symdot(A, B)
    symmetrize!(A)
    symmetrize!(B)
    return dot(A.data, B.data)
end

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
# Deprecated methods
using NNlib: DenseConvDims
import NNlib: stride, padding, dilation, flipkernel, spatial_dims, kernel_size,
    conv!, ∇conv_filter!, ∇conv_data!,
    maxpool!, meanpool!, ∇maxpool!, ∇meanpool!, PoolDims

using CUDA.CUDNN: scalingParameter, CUDNN_CONVOLUTION, convdims,
                  cudnnConvolutionDescriptor, cudnnConvolutionBwdDataAlgoPerf,
                  cudnnConvolutionForward!, cudnnConvolutionBwdFilterAlgoPerf,
                  cudnnConvolutionBackwardData, cudnnConvolutionBackwardFilter,
                  cudnnConvolutionBackwardBias

const CUDNNFloat = Union{Float16,Float32,Float64}

# Since CUDNN does not support 1D convolution, Conv in Flux will give a CUDNNError if the size is 1-dimensional.
fix1d(x) = x
fix1d(x::DenseCuArray{T, 3}) where T = reshape(x, 1, size(x, 1), size(x, 2), size(x, 3))
fix1d(cdims::DenseConvDims{1,K,C_in,C_out,S,P,D,F}) where {K,C_in,C_out,S,P,D,F} =
    DenseConvDims{2,(1,K...),C_in,C_out,(1,S...),(0,0,P...),(1,D...),F}((1,cdims.I...))
fix1d(pdims::PoolDims{1,K,S,P,D}) where {K,S,P,D,F} =
    PoolDims{2,(1,K...),(1,S...),(0,0,P...),(1,D...)}((1,pdims.I...), pdims.C_in)

# Convolution

function cudnnConvolutionDescriptor(cdims::DenseConvDims, x::DenseCuArray{T}) where T
    cdims, x = fix1d(cdims), fix1d(x)
    mode=(NNlib.flipkernel(cdims) ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION)
    cudnnConvolutionDescriptor(convdims(nnlibPadding(cdims),size(x),0),
                               convdims(NNlib.stride(cdims),size(x),1),
                               convdims(NNlib.dilation(cdims),size(x),1),
                               mode,
                               cudnnDataType(T),
                               math_mode(),
                               CUDNN_DEFAULT_REORDER,
                               Cint(1))
end

function conv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}, cdims::DenseConvDims;
               alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d = cudnnConvolutionDescriptor(cdims, x)
    cudnnConvolutionForward!(y, w, x, d; alpha, beta, z=y)
end

function NNlib.conv_bias_act!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T},
                            cdims::DenseConvDims, bias::DenseCuArray{T}, σ=identity;
                            z::DenseCuArray{T}=y, alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    d = cudnnConvolutionDescriptor(cdims, x)
    # only relu and identity are supported by cudnnConvolutionForward!
    activation = (σ == NNlib.relu ? CUDNN_ACTIVATION_RELU : CUDNN_ACTIVATION_IDENTITY)
    cudnnConvolutionForward!(y, w, x, d; z, bias, activation, alpha, beta)
    if activation === CUDNN_ACTIVATION_IDENTITY && σ ∉ (nothing, identity)
        y = σ.(y)
    end
    return y
end

function ∇conv_data!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, w::DenseCuArray{T},
                     cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(dx), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(w)
    convDesc = cudnnConvolutionDescriptor(cdims, dx)
    p = cudnnConvolutionBwdDataAlgoPerf(wDesc, w, yDesc, dy, convDesc, xDesc, dx)
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardData(handle(), alpha, wDesc, w, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, xDesc, dx)
    end
    return dx
end

function ∇conv_filter!(dw::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                       cdims::DenseConvDims; alpha=1, beta=0, algo=-1) where T<:CUDNNFloat
    if cudnnversion() < v"6"
        all(x -> x == 1, dilation(cdims)) || error("Only dilation = 1 is supported in cuDNN version < 6")
    end
    if algo != -1
        @warn "The algo option has been deprecated, the fastest algo is computed automatically" maxlog=1
    end
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta);
    xDesc, yDesc, wDesc = cudnnTensorDescriptor(x), cudnnTensorDescriptor(dy), cudnnFilterDescriptor(dw)
    convDesc = cudnnConvolutionDescriptor(cdims, x)
    p = cudnnConvolutionBwdFilterAlgoPerf(xDesc, x, yDesc, dy, convDesc, wDesc, dw);
    with_workspace(p.memory) do workspace
        cudnnConvolutionBackwardFilter(handle(), alpha, xDesc, x, yDesc, dy, convDesc, p.algo, workspace, sizeof(workspace), beta, wDesc, dw);
    end
    return dw
end

function ∇conv_bias!(db::DenseCuArray{T}, dy::DenseCuArray{T}; alpha=1, beta=0) where T<:CUDNNFloat
    alpha,beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    bDesc, yDesc = cudnnTensorDescriptor.((db,dy))
    cudnnConvolutionBackwardBias(handle(), alpha, yDesc, dy, beta, bDesc, db)
    return db
end

# Compatibility shims until users upgrade to new NNlib format
function conv!(y::DenseCuArray{T}, x::DenseCuArray{T}, w::DenseCuArray{T}; pad=0, stride=1, flipkernel=0, dilation=1, kwargs...) where {T<:CUDNNFloat}
    cdims = DenseConvDims(x, w; padding=pad, stride=stride, flipkernel=(flipkernel!=0), dilation=dilation)
    return conv!(y, x, w, cdims; kwargs...)
end

function ∇conv_filter!(dw::DenseCuArray{T}, dy::DenseCuArray{T}, x::DenseCuArray{T}; pad=0, stride=1, flipkernel=0, dilation=1, kwargs...) where {T<:CUDNNFloat}
    cdims = DenseConvDims(x, dw; padding=pad, stride=stride, flipkernel=(flipkernel!=0), dilation=dilation)
    # NOTE!!!  This compat shim re-arranges the argument order!
    return ∇conv_filter!(dw, x, dy, cdims; kwargs...)
end

function cudnnConvolutionForward(y::DenseCuArray{T,N}, x::DenseCuArray{T,N}, w::DenseCuArray{T,N},
                                 cdims::DenseConvDims; algo=0, alpha=1, beta=0) where {T,N}
    # @warn "`cudnnConvolutionForward(y,x,w,c::DenseConvDims)` is deprecated, please use one of the methods in `@doc cudnnConvolutionForward!`." maxlog=1
    cudnnConvolutionForward!(y, w, x; alpha, beta, padding=nnlibPadding(cdims), stride=NNlib.stride(cdims), dilation=NNlib.dilation(cdims), mode=(NNlib.flipkernel(cdims) ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION))
end

function cudnnConvolutionBiasActivationForward(y::DenseCuArray{T,N}, x::DenseCuArray{T,N}, w::DenseCuArray{T,N}, z::DenseCuArray{T,N}, bias::DenseCuArray{T,N},
                                               cdims::DenseConvDims; algo=0, alpha1=1, alpha2=1,
                                               activationMode=CUDNN_ACTIVATION_RELU, activationCoeff=0.0, activationReluNanOpt=CUDNN_NOT_PROPAGATE_NAN) where {T,N}
    # @warn "`cudnnConvolutionBiasActivationForward` is deprecated, please use one of the methods in `@doc cudnnConvolutionForward!`." maxlog=1
    cudnnConvolutionForward!(y, w, x; bias, activation=activationMode, z, alpha=alpha1, beta=alpha2, padding=nnlibPadding(cdims), stride=NNlib.stride(cdims), dilation=NNlib.dilation(cdims), mode=(NNlib.flipkernel(cdims) ? CUDNN_CROSS_CORRELATION : CUDNN_CONVOLUTION))
end


text:
@testset "Testing Ensemble Kalman Filter..." begin

using StateSpace, LinearAlgebra

# Define process model: Ricker growth model, with two populations
r1 = 0.03 							# 10% intrinsic growth rate for pop 1
r2 = 0.06 							# 14% intrinsic growth rate for pop 1
K1 = 20000 							# carrying capacity for pop 1
K2 = 10000 							# carrying capacity for pop 2

# process model
function f(x)
	[x[1] .* exp(r1 * (1 - (x[1] ./ K1))),
	 x[2] .* exp(r2 * (1 - (x[2] ./ K2)))]
end

V = Matrix(Diagonal([3000.0, 5000.0])) 		# covariance of (additive) process noise

# Define observation model
g(x) = 0.01 .* x
W = Matrix(Diagonal([100.0, 250.0]))

mod = NonlinearGaussianSSM(f, V, g, W)
x0 = MvNormal([1000.0, 500.0], [100.0 0.0; 0.0 100.0])

x1 = predict(mod, x0)
# println(mean(x1))

y1 = observe(mod, x1)
u1 = update(mod, x1, mean(y1))
# println(mean(u1))

xx, yy = simulate(mod, 200, x0)
yy[1, 50] = NaN

filt = EnKF()
fs = filter(mod, yy, x0, filt)
fs_ensemble = filter(mod, yy, x0, filt, return_ensemble=true)
@test size(fs_ensemble) == (mod.nx, filt.nparticles, size(yy, 2))
y_new = fs.observations[:, end] + randn(size(mod.V(1), 1)) / 10
update!(mod, fs, y_new)

ss = smooth(mod, fs)
@test loglikelihood(fs) < loglikelihood(ss)

end


text:
# Autogenerated wrapper script for cilantro_jll for armv6l-linux-gnueabihf-cxx03
export libcilantro

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("cilantro")
JLLWrappers.@declare_library_product(libcilantro, "libcilantro.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcilantro,
        "lib/libcilantro.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Base: show
using SHA
using StanSample
using Parquet
using DataFrames
using VersionedFile

mutable struct VersionedModel
    name::String
    # Stan Model
    stan_model::SampleModel

    # Versioning Detail Fields
    hash::String
    # Storage Information
    versionfile::VersionedFile
    version::Int
end

function VersionedModel(name::AbstractString, model::AbstractString, version_store_file::AbstractString)
    stan_model = SampleModel(name, model)
    hash = bytes2hex(sha3_512(lowercase(strip(model))))
    version = 1

    if isfile(version_store_file)
        versionstoredf = DataFrame(read_parquet(version_store_file))
        modelfromfile = filter(:name => modelnames -> modelnames == name, versionstoredf)
        if length(modelfromfile) == 0
            df = DataFrame(Name = [name], Model=[model], Version=[version_number], Hash=[hash])
            append!(versionstoredf, df)
        else
            duplicatemodelfromfile = filter(:hash => modelhash -> modelhash == hash, modelfromfile)
            if length(duplicatemodelfromfile) > 0
                version_number = duplicatemodelfromfile[:hash, 1]
            else
                version_number = max(duplicatemodelfromfile[:Version]) + 1
                df = DataFrame(Name = [name], Model=[model], Version=[version_number], Hash=[hash])
                append!(versionstoredf, df)
            end
        end
    else
        versionstoredf = DataFrame(Name = [name], Model=[model], Version=[version_number], Hash=[hash])
    end
    VersionedModel(stan_model, hash, version_number)
end


text:
using ImageFiltering: kernelfactors, imfilter!, imfilter, centered

"""
    inpaint_nans!(flow::Flow)

Inpaint NaNs in `flow` using surrounding non-NaN values.

# References
Oliveira, Manuel & Bowen, Brian & McKenna, Richard & Chang, Yu-Sung. (2001). Fast Digital Image Inpainting.. 261-266.
[Here.](https://www.researchgate.net/publication/220903053_Fast_Digital_Image_Inpainting)

See also: [`showflow`](@ref), [`Flow`](@ref)
"""
function inpaint_nans!(flow::Flow)

    # make a mask where nans are 0 and the rest is 1
    mask = .!isnan.(real(flow))
    @views flow[.!mask] .= 0

    # init interpolation kernel
    ker = centered([0.5, 2, 0.5])
    kernf = kernelfactors((ker, ker))

    zero_in_mask = any(mask .== 0)

    # init tmp vars
    div_coef = similar(flow, Float16)
    estim_flow = similar(flow)

    # until there 0s left in mask
    while zero_in_mask

        imfilter!(div_coef, mask, kernf, "symmetric")
        imfilter!(estim_flow, mask .* flow, kernf, "symmetric")

        # identify the newly estimated
        newly_inpainted = ((div_coef .!= 0) .& (mask .== 0))

        # set newly inpainted
        @views flow[newly_inpainted] .= estim_flow[newly_inpainted] ./ div_coef[newly_inpainted]

        # update mask
        # putting .= 1 here is a badly Errored bug.
        @views mask[newly_inpainted] .= true

        # update while condition
        zero_in_mask = any(mask .== 0)
    end

    return nothing
end


text:
export interpolateFace
#
# Given variables q at elem_volume nodes, compute corresponding gradients
# TODO: a higher level function
#

function calcGradient(mesh::AbstractDGMesh{Tmsh},
                      sbp::AbstractOperator{Tsbp},
                      eqn::EllipticData{Tsol, Tres, Tdim},
                      q::AbstractArray{Tsol, 3},
                      q_grad::AbstractArray{Tsol, 4}) where {Tmsh, Tsol, Tres, Tdim, Tsbp}
  @assert(size(q, 3) == mesh.numEl)
  @assert(size(q, 2) == mesh.numNodesPerElement)
  @assert(size(q, 1) == mesh.numDofPerNode)

  @assert(size(q_grad, 4) == Tdim)
  @assert(size(q_grad, 3) == mesh.numEl)
  @assert(size(q_grad, 2) == mesh.numNodesPerElement)
  @assert(size(q_grad, 1) == mesh.numDofPerNode)

  numElems = mesh.numEl
  numNodes = mesh.numNodesPerElement
  numVars = mesh.numDofPerNode
  numNodesPerElement = mesh.numNodesPerElement
  dim = Tdim
  # dim = Tdim

  q_grad[:,:,:,:] = 0.0

  Dx = Array{Tsbp}(numNodesPerElement, numNodesPerElement, 2)
  for e=1:numElems
    # First compute Dx for this element
    calcDx(mesh, sbp, Int(e), Dx)
    # q_grad = (Dxq, Dyq), Dx = H^(-1)Qx
    for n=1:numNodesPerElement
      for ivar=1:numVars
        for d=1:dim
          for col=1:numNodesPerElement
            q_grad[ivar, n, e, d] += Dx[n,col,d]*q[ivar, col, e]
          end
        end
      end
    end
  end
end

function interpolateFace(mesh::AbstractDGMesh,
                         sbp,
                         eqn,
                         opts,
                         q::Abstract3DArray,
                         q_face::AbstractArray{Tsol, 4}) where Tsol
  # interpolate solution
  interiorfaceinterpolate!(mesh.sbpface, mesh.interfaces, q, q_face)
end

function interpolateFace(mesh::AbstractDGMesh,
                         sbp,
                         eqn::EllipticData{Tsol, Tres, Tdim},
                         opts,
                         grad::AbstractArray{Tsol, 4},
                         grad_face::AbstractArray{Tsol, 5}) where {Tsol, Tres, Tdim}
  # interpolate gradient of solution
  for d=1: Tdim
    dqdx = sview(grad, :, :, :, d)
    dqdx_face = sview(grad_face, :, :, :, :, d)
    interiorfaceinterpolate!(mesh.sbpface, mesh.interfaces, dqdx, dqdx_face)
  end
end

#
# Calculate flux at edge cubature points using face-based form.
#
function calcFaceFlux(mesh::AbstractDGMesh{Tmsh},
                             sbp::AbstractOperator,
                             eqn::EllipticData{Tsol, Tres, Tdim},
                             opts,
                             interfaces::AbstractArray{Interface,1},
                             xflux_face::AbstractArray{Tres, 3},
                             yflux_face::AbstractArray{Tres, 3},
                             flux_face::AbstractArray{Tres, 3}) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  penalty_method = opts["Flux_name"]
  p = opts["order"]
  Cip = opts["Cip"]
  dq = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  factor_shahbazi = Float64(p + 1.0)*Float64(p + Tdim)/(2.0*Tdim)
  sbpface = mesh.sbpface
  numFacesPerElem = 3
  nrm = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  nrm1 = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  area = Array{Tmsh}(mesh.numNodesPerFace)
  lambda_dqdxL = Array{Tsol}(Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  lambda_dqdxR = Array{Tsol}(Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  eigMaxL = Array{Tmsh}(mesh.numDofPerNode)
  eigMaxR = Array{Tmsh}(mesh.numDofPerNode)

  #
  # For scalar penalty, we use penalty;
  # for matrix penalty, we use Sat.
  #
  penalty      = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  penalty_sat0 = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace)
  penalty_sat  = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace)
  penalty_shahbazi = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace)

  Sat = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerFace)
  #
  # |R 0| |Λxx Λxy| |R 0|^T
  # |0 R| |Λyx Λyy| |0 R|
  #
  sbpface = mesh.sbpface
  R = sview(sbpface.interp, :,:)
  RLR_L = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerFace)
  RLR_R = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerFace)
  RLR = Array{Tmsh}(mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerFace)
  BsqrtRHinvRtBsqrt = Array{Tmsh}(mesh.numNodesPerFace, mesh.numNodesPerFace)
  HRBRH = Array{Tmsh}(mesh.numNodesPerElement, mesh.numNodesPerElement)
  area_sum = sview(eqn.area_sum, :)
  perm = zeros(Tmsh, sbp.numnodes, sbpface.stencilsize)
  Hinv = zeros(Tmsh, sbp.numnodes, sbp.numnodes)
  Bsqrt = zeros(Tmsh, sbpface.numnodes, sbpface.numnodes)
  for s = 1:sbpface.stencilsize
    perm[sbpface.perm[s, 1], s] = 1.0
  end
  for i = 1:sbp.numnodes
    Hinv[i,i] = 1.0/sbp.w[i]
  end
  for i = 1:sbpface.numnodes
    Bsqrt[i,i] = sqrt(sbpface.wface[i])
  end

  BsqrtRHinvRtBsqrt = Bsqrt*R.'*perm.'*Hinv*perm*R*Bsqrt
  sigma = eigmax(BsqrtRHinvRtBsqrt)
  relax_coef = 1.0
  if haskey(opts, "unstable_coef")
    relax_coef = opts["unstable_coef"]
  end
  for f = 1:nfaces    # loop over faces
    face = interfaces[f]
    eL = face.elementL
    eR = face.elementR
    fL = face.faceL
    fR = face.faceR
    permL = sview(sbpface.perm, :, fL)
    permR = sview(sbpface.perm, :, fR)
    stencilSize = sbpface.stencilsize

    # Compute geometric info on face
    nrm = ro_sview(mesh.nrm_face, :, :, f)
    for n=1:mesh.numNodesPerFace
      # area[n] = norm(view(nrm, :, n))
      area[n] = norm(ro_sview(nrm, :, n))
        # norm vector in physical domain without any scale, ie, |nrm|=1
      nrm1[:, n] = nrm[:, n]/area[n]
      nrm1[:, n] = nrm[:, n]/area[n]
    end

    # Compute the size of element and face, and then meas(face)/meas(elem)
    elem_volL = 0.0
    elem_volR = 0.0
    face_area = 0.0
    for n = 1:mesh.numNodesPerElement
      elem_volL += sbp.w[n]/mesh.jac[n, eL]
      elem_volR += sbp.w[n]/mesh.jac[n, eR]
      # elem_volL += eqn.w[n, eL]
      # elem_volR += eqn.w[n, eR]
    end

    for n = 1:mesh.numNodesPerFace
      face_area += sbpface.wface[n]*area[n]
    end

    if penalty_method == "Shahbazi" || penalty_method == "SAT0"
      heL = elem_volL/eqn.area_sum[eL]
      heR = elem_volR/eqn.area_sum[eR]
      he = min(heL, heR)
    end
    area_weightL = area_sum[eL]/face_area
    area_weightR = area_sum[eR]/face_area

    #
    # compute lambda_dqdx
    #

    #
    # face variables
    #
    lambdaL = sview(eqn.lambda_face, :, :, :, 1, :, f)
    lambdaR = sview(eqn.lambda_face, :, :, :, 2, :, f)
    dqdxL = sview(eqn.q_grad_face, :, 1, :, f, 1)
    dqdyL = sview(eqn.q_grad_face, :, 1, :, f, 2)
    dqdxR = sview(eqn.q_grad_face, :, 2, :, f, 1)
    dqdyR = sview(eqn.q_grad_face, :, 2, :, f, 2)

    for n = 1:mesh.numNodesPerFace
      for var = 1:mesh.numDofPerNode
        lambda_dqdxL[1, var, n] = lambdaL[1, 1, var, n]*dqdxL[var, n] + lambdaL[1, 2, var, n]*dqdyL[var, n]
        lambda_dqdxL[2, var, n] = lambdaL[2, 1, var, n]*dqdxL[var, n] + lambdaL[2, 2, var, n]*dqdyL[var, n]
        lambda_dqdxR[1, var, n] = lambdaR[1, 1, var, n]*dqdxR[var, n] + lambdaR[1, 2, var, n]*dqdyR[var, n]
        lambda_dqdxR[2, var, n] = lambdaR[2, 1, var, n]*dqdxR[var, n] + lambdaR[2, 2, var, n]*dqdyR[var, n]
      end
    end

    #
    # First compute penalty
    #
    if penalty_method == "SAT0"  ### SATs, Cip = -1
      #
      # Compute SAT matrix (penalty parameter)
      # First compute the eigenvalue of
      #
      #	|Λxx Λxy| |J 0|
      #	|Λyx Λyy| |0 J|
      #
      eigMaxL[:] = -1.0
      eigMaxR[:] = -1.0
      lambdaL = sview(eqn.lambda, :,:,:,:,eL)
      lambdaR = sview(eqn.lambda, :,:,:,:,eR)

      for dof = 1 : mesh.numDofPerNode
        for n = 1 : mesh.numNodesPerElement
          # left element
          b = real(lambdaL[1,1,dof,n] + lambdaL[2,2,dof,n])
          ac = real(lambdaL[1,1,dof,n] * lambdaL[2,2,dof,n] - lambdaL[1,2,dof,n] * lambdaL[2,1,dof,n])
          root = 0.5*(b + sqrt(b*b - 4*ac))*mesh.jac[n, eL]
          eigMaxL[dof] = max(eigMaxL[dof], root)

          # right element
          b = real(lambdaR[1,1,dof,n] + lambdaR[2,2,dof,n])
          ac = real(lambdaR[1,1,dof,n] * lambdaR[2,2,dof,n] - lambdaR[1,2,dof,n] * lambdaR[2,1,dof,n])
          root = 0.5*(b + sqrt(b*b - 4*ac))*mesh.jac[n, eR]
          eigMaxR[dof] = max(eigMaxR[dof], root)
        end

        for n = 1:mesh.numNodesPerFace
          penalty[dof, n] = area_weightL*eigMaxL[dof] + area_weightR*eigMaxR[dof]
          # penalty[dof, n] = 2.0*he*max(eigMaxL[dof], eigMaxR[dof])
          # penalty[dof, n] *= 0.125*sigma
          penalty[dof, n] *= 0.25*sigma*area[n]*area[n]
        end
      end

    elseif  penalty_method == "Shahbazi"
      #
      # Reference:
      # Shahbazi, Mavriplis, Multigrid algorithms for high order discontinuous
      # Galerkin discretization of the compressible Navier-Stokes equations,
      # JCP, volume 228 Issue 21, November, 2009
      #
      lambda_faceL = sview(eqn.lambda_face, :,:, :, 1, :, f)
      lambda_faceR = sview(eqn.lambda_face, :,:, :, 2, :, f)

      #
      # Compute n^TΛn
      #
      for n = 1:mesh.numNodesPerFace
        for dof = 1:mesh.numDofPerNode
          penalty[dof, n] =  nrm1[1, n]*nrm1[1, n]*(lambda_faceL[1, 1, dof, n] + lambda_faceR[1,1,dof, n])
          penalty[dof, n] += nrm1[1, n]*nrm1[2, n]*(lambda_faceL[1, 2, dof, n] + lambda_faceR[1,2,dof, n])
          penalty[dof, n] += nrm1[2, n]*nrm1[1, n]*(lambda_faceL[2, 1, dof, n] + lambda_faceR[2,1,dof, n])
          penalty[dof, n] += nrm1[2, n]*nrm1[2, n]*(lambda_faceL[2, 2, dof, n] + lambda_faceR[2,2,dof, n])
          penalty[dof, n] *= 0.5*factor_shahbazi*area[n]/he
        end
      end

    elseif penalty_method == "SAT"
      #
      lambdaL = sview(eqn.lambda, :,:,:,:,eL)
      lambdaR = sview(eqn.lambda, :,:,:,:,eR)
      Sat[:,:,:] = 0.0

      # RLR_L and RLR_R are reused for each (d1, d2)
      for d1 = 1:Tdim
        for d2 = 1:Tdim
          # Compute R Λ{d1,d2}/H R^T
          # RLR_L[:,:,:] = 0.0
          # RLR_L[:,:,:] = 0.0
          RLR[:,:,:] = 0.0
          for n1 = 1:mesh.numNodesPerFace
            for n2 = 1:mesh.numNodesPerFace
              for k = 1:stencilSize
                # RLR_L[:, n2, n1] += lambdaL[d1, d2, :, permL[k]] / eqn.w[permL[k], eL] * R[k, n2] * R[k, n1]
                # RLR_R[:, n2, n1] += lambdaR[d1, d2, :, permR[k]] / eqn.w[permR[k], eR] * R[k, n2] * R[k, n1]
                # RLR[:, n2, n1] += (lambdaL[d1, d2, :, permL[k]] / eqn.w[permL[k], eL]+  lambdaR[d1, d2, :, permL[k]] / eqn.w[permR[k], eL])* R[k, n2] * R[k, n1]
                RLR[:, n2, n1] += (lambdaL[d1, d2, :, permL[k]] / eqn.w[permL[k], eL]*area_weightL+  lambdaR[d1, d2, :, permR[k]] / eqn.w[permR[k], eR] *area_weightR)* R[k, n2] * R[k, n1]
              end
            end
          end

          # N_{d1}i} RLR N_{d2}
          for n1 = 1:mesh.numNodesPerFace
            for n2 = 1:mesh.numNodesPerFace
              Sat[:, n2, n1] += nrm1[d1, n2]*RLR[:, n2, n1]*nrm1[d2, n1]
            end
          end
        end
      end

      #
      # Left multiply by numFacesPerElem/4*B
      #
      for row = 1:mesh.numNodesPerFace
        for col = 1:mesh.numNodesPerFace
          # Sat[:, row, col] *= 0.25*numFacesPerElem*sbpface.wface[col]*area[col]
          Sat[:, row, col] *= 0.25*sbpface.wface[col]*area[col]
        end
      end
    else
      error("We should never get here")
      for n = 1:mesh.numNodesPerFace
        for dof = 1:mesh.numDofPerNode
          penalty[dof, n] = Cip*p*p/he
        end
      end

    end

    qL = sview(eqn.q_face, :, 1, :, f)
    qR = sview(eqn.q_face, :, 2, :, f)

    for n = 1:mesh.numNodesPerFace
      for dof = 1:mesh.numDofPerNode
        dq[dof, n] = qL[dof, n] - qR[dof, n]
      end
    end

    for n=1:mesh.numNodesPerFace
      # the other is `[[v]]⋅f(u+, u-) = (v+ - v-)g(u+, u-)⋅norm`.
      # vector f(u+, u-) is stored in xflux_face and yflux_face,
      # while scalar g(u+, u-) is stored in flux_face.

      # term-2 vector [[q]]
      xflux_face[:, n, f] = -0.5*dq[:,n]*nrm[1, n]
      yflux_face[:, n, f] = -0.5*dq[:,n]*nrm[2, n]

      # term-3, {{∇q}}⋅norm
      flux_face[:, n, f] = -0.5*((lambda_dqdxL[1, :, n]+ lambda_dqdxR[1, :, n])*nrm[1, n] + (lambda_dqdxL[2, :, n]+ lambda_dqdxR[2, :, n])*nrm[2, n])

    end
    #
    # term-4, penalty term
    #

    # area comes from dΓ = area*dΓ_{ξ}
    if penalty_method == "Shahbazi" || penalty_method == "SAT0"|| penalty_method == "Hartman"
      for n=1:mesh.numNodesPerFace
        for dof = 1:mesh.numDofPerNode
          flux_face[dof, n, f] += penalty[dof, n]*dq[dof, n] * relax_coef
        end
      end
    elseif penalty_method == "SAT"
      for n=1:mesh.numNodesPerFace
        for n1 = 1:mesh.numNodesPerFace
          for dof = 1:mesh.numDofPerNode
            flux_face[dof, n, f] += Sat[dof, n, n1]*dq[dof, n1]*area[n] * relax_coef
          end
        end
      end
    else
      error("We should never get here")
    end
  end # end of loop over all interfaces

  return nothing
end # end of function calcFaceFlux


text:
#####
##### jit/call interface with GPUCompiler
#####

struct MixtapeCompilerTarget <: AbstractCompilerTarget end
struct MixtapeCompilerParams <: AbstractCompilerParams
    ctx::Any
end

struct Entry{F,RT,TT}
    f::F
    specfunc::Ptr{Cvoid}
    func::Ptr{Cvoid}
end

const jit_compiled_cache = Dict{UInt,Any}()

function jit(ctx::CompilationContext, f::F, tt::TT=Tuple{}) where {F,TT<:Type}
    fspec = FunctionSpec(f, tt, false, nothing) #=name=#
    job = CompilerJob(MixtapeCompilerTarget(), fspec, MixtapeCompilerParams(ctx))
    return GPUCompiler.cached_compilation(jit_compiled_cache, job, _jit, _jitlink)
end

@doc(
"""
    jit(ctx::CompilationContext, f::F, tt::TT = Tuple{}) where {F, TT <: Type}

Compile and specialize a method instance for signature `Tuple{f, tt.parameters...}` with pipeline parametrized by `ctx`.

Returns a callable "thunk" `Entry{F, RT, TT}` where `RT` is the return type of the instance after inference.
""", jit)

function _jitlink(job::CompilerJob, (rt, llvm_mod, func_name, specfunc_name))
    fspec = job.source
    jitted_mod = compile!(orc[], llvm_mod)
    specfunc_addr = addressin(orc[], jitted_mod, specfunc_name)
    specfunc_ptr = pointer(specfunc_addr)
    func_addr = addressin(orc[], jitted_mod, func_name)
    func_ptr = pointer(func_addr)
    if specfunc_ptr === C_NULL || func_ptr === C_NULL
        @error "Compilation error" fspec specfunc_ptr func_ptr
    end
    return Entry{typeof(fspec.f),rt,fspec.tt}(fspec.f, specfunc_ptr, func_ptr)
end

function _jit(job::CompilerJob)
    rt, llvm_specfunc, llvm_func, llvm_mod = codegen(job.params.ctx, job.source.f, job.source.tt)
    specfunc_name = LLVM.name(llvm_specfunc)
    func_name = LLVM.name(llvm_func)
    linkage!(llvm_func, LLVM.API.LLVMExternalLinkage)
    linkage!(llvm_specfunc, LLVM.API.LLVMExternalLinkage)
    run_pipeline!(llvm_mod)
    return (rt, llvm_mod, func_name, specfunc_name)
end

#####
##### Experimental: support for AOT compilation
#####

# TODO: When https://github.com/JuliaGPU/GPUCompiler.jl/compare/jps/static-compile is merged -- refactor.

module MixtapeRuntime
    # dummy methods
    signal_exception() = return
    malloc(sz) = C_NULL
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return

    # for validation
    sin(x) = Base.sin(x)
end

import GPUCompiler: runtime_module
GPUCompiler.runtime_module(::CompilerJob{<:Any, MixtapeCompilerParams}) = MixtapeRuntime

const linker = Sys.isunix() ? "ld.lld" : Sys.isapple() ? "ld64.lld" : "lld-link"

function mixtape_job(ctx, @nospecialize(func), tt;
               cpu::String = (LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUName()),
               features::String=(LLVM.version() < v"8") ? "" : unsafe_string(LLVM.API.LLVMGetHostCPUFeatures()),
               name = GPUCompiler.safe_name(repr(func)))
    source = FunctionSpec(func, tt, false)
    target = GPUCompiler.NativeCompilerTarget(cpu=cpu, features=features)
    params = MixtapeCompilerParams(ctx)
    return CompilerJob(target, source, params)
end

function aot(ctx::CompilationContext, f::F, tt::TT=Tuple{};
        path = tempname()) where {F, TT <: Type}
    open(path, "w") do io
        job = mixtape_job(ctx, f, tt)
        rt, _, _, mod = codegen(job.params.ctx, job.source.f, job.source.tt)
        GPUCompiler.finish_module!(job, mod)
        tm = GPUCompiler.llvm_machine(job.target)
        LLVM.emit(tm, mod, LLVM.API.LLVMObjectFile, path)
    end
    return (path, name)
end

#####
##### Call interface
#####

@generated function (entry::Entry{F,RT,TT})(args...) where {F,RT,TT}
    expr = quote
        args = Any[args...]
        ccall(entry.func, Any, (Any, Ptr{Any}, Int32), entry.f, args, length(args))
    end
    return expr
end

macro load_call_interface()
    expr = quote
        function cached_call(entry::Mixtape.Entry{F,RT,TT}, args...) where {F,RT,TT}

            # TODO: Fast ABI.
            #ccall(entry.func, Any, (Any, $(nargs...), ), entry.f, $(_args...))

            # Slow ABI. Requires an array allocation.
            expr = quote
                vargs = Any[args...]
                ccall($(entry.func), Any, (Any, Ptr{Any}, Int32), $(entry.f), vargs,
                      $(length(args)))
            end
            return expr
        end

        @generated function call(ctx::Mixtape.CompilationContext, f::F, args...) where {F}
            TT = Tuple{args...}
            entry = Mixtape.jit(ctx(), F.instance, TT)
            return cached_call(entry, args...)
        end
    end
    return esc(expr)
end

@doc(
"""
    @load_call_interface()

A macro which expands to load a generated function `call` into the scope of the calling module. This generated function can be applied to signature argument types `Tuple{ctx<:CompilationContext, f<:Function, args...}`. `call` then creates a new instance of `ctx` and calls `Mixtape.jit` -- it then caches a `ccall` which calls a function pointer to the GPUCompiler-compiled LLVM module.

The `call` interface currently uses a slow ABI `ccall` -- which costs an array allocation for each toplevel `call`. This allocation is required to construct a `Vector{Any}` for the arguments and pass a pointer to it over the line, where the call unboxes each argument.
""", :(@Mixtape.load_call_interface))


text:
module TestBaselet
using Test

@testset "$file" for file in sort([
    file for file in readdir(@__DIR__) if match(r"^test_.*\.jl$", file) !== nothing
])
    # Skip inference test on Julia 1.0.  It may work in 1.1 or 1.2 but
    # they are note tested.
    if VERSION < v"1.3-" && file == "test_inference.jl"
        @info "Skip $file for Julia $VERSION"
        continue
    elseif (
        lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true" &&
        file == "test_inference.jl"
    )
        @info "Skip $file on PkgEval."
        continue
    end
    include(file)
end

end  # module


text:
# A more precise example of $\text H_2$ with 4 qubits hamiltonian.

using Yao
using YaoExtensions

# Here we we invoke [OpenFermion](https://github.com/quantumlib/OpenFermion) in Python
# to get the hamiltonian under [Jordan-Weigner transformation](https://en.wikipedia.org/wiki/Jordan%E2%80%93Wigner_transformation)
# (with 2 electrons and 4 sto-3g basis set orbitals).

using PyCall

py"""
import numpy as np
from openfermion.config import *
from openfermion.hamiltonians import MolecularData
from openfermion.transforms import jordan_wigner

def make_H2_mol(bond_len):

    ## Returns:
    ##     molecule(openfermion.hamiltonians.MolecularData):
    ## MolecularData for H2 with certain bond length

    ## load molecule
    atom_1 = 'H'
    atom_2 = 'H'
    basis = 'sto-3g'
    multiplicity = 1
    charge = 0
    coordinate_1 = (0.0, 0.0, 0.0)
    coordinate_2 = (0.0, 0.0, bond_len)
    geometry = [(atom_1, coordinate_1), (atom_2, coordinate_2)]
    molecule = MolecularData(geometry, basis, multiplicity,
    charge, description=str(bond_len))
    ## molecule = run_pyscf(molecule,run_scf=1,run_ccsd=1,run_fci=1)
    ## can calculate on arbitrary bond length if package openfermionpyscf is available,
    ## only existing molecule data can be used
    molecule.load()
    return molecule

def get_hamiltonian(bond_len):
    ## get the coefficients and pauli terms of Hamiltonian of H_2 with given bond length
    ## Returns:
    ##     (PyArray,PyArray)
    molecule = make_H2_mol(bond_len)
    qubit_hamiltonian = jordan_wigner(molecule.get_molecular_hamiltonian())
    qubit_hamiltonian.compress()
    return list(qubit_hamiltonian.terms.keys()),list(qubit_hamiltonian.terms.values())

"""

# Define the function to transfer OpenFermion hamiltonian to Yao circuit

function get_fermion_hamiltonian(n::Int,terms::Array,coefs::Vector)
    gates=Dict("Z"=>Z,"X"=>X,"Y"=>Y)
    to_pauli(t::Tuple{Int,String})=put(n,t[1]+1=>get(gates,t[2],ErrorException("Invalid")))
    return sum(coefs[2:end].*map(x->reduce(*,map(to_pauli,x)),terms[2:end]))
end

# Use VQE code given in the previous Yao example

using Flux: Optimise

function train!(circ, hamiltonian; optimizer, niter::Int=200,verbose::Bool=true)
     params = parameters(circ)
     dispatch!(circ, :random)
     for i=1:niter
         _, grad = expect'(hamiltonian, zero_state(nqubits(circ)) => circ)
         Optimise.update!(optimizer, params, grad)
         dispatch!(circ, params)
         if verbose
            println("Energy =
            $(expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ) |> real)")
         end
     end
     return expect(hamiltonian, zero_state(nqubits(hamiltonian)) |> circ)
end

# Train VQE to search ground energy on various bond lengths.

bond_lens = Array(0.2:0.1:1.5)#\AA
pes = Vector{Real}()

for l in bond_lens
    terms, coefs = py"get_hamiltonian"(l)
    e = coefs[1]

    h = get_fermion_hamiltonian(4,terms,coefs)
    c = variational_circuit(4)
    emin_vqe = train!(c, h; optimizer=Optimise.ADAM(0.1),verbose=false)
    println("Estimated minimum energy at bond length $(l): $(e+emin_vqe|>real)")
    push!(pes,e+emin_vqe|>real)

    using LinearAlgebra
    emin = eigvals(Matrix(mat(h)))[1]
    @assert isapprox(emin, emin_vqe, atol=1e-2)
end

# Plot the potential energy surface (PES)

using Plots

plot(bond_lens,pes, markershape = :circle, linestyle = :dash, label = "VQE",ylabel = "Ground State Energy (Hartree)", xlabel = "Bond Length (Å) ")

# Since OpenFermion only generates molecule data at discrete bond length with minimum step to be 0.1 Å,
#
# the VQE gives lowest energy -1.1361862609 Hartree at 0.7 Å.
#
# As a comparison, the ground energy at stable bond length given by FCI (Full Configuration Interaction) is -1.13730688  Hartree at 0.734832 Å.


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function domainmaps(a::Array{Float64,1}, b::Array{Float64,1})
    # a, b, c are three vectors describing the simulation box
    # a = [a1, a2, a3]
    # b = [b1, b2, b3]
    # c = [c1, c2, c3]

    R2B = [a b]; # transform the unit square to the simulation box (a, b)
    B2R = inv(R2B); # transform the simulation box (a, b) to the unit square

    return B2R, R2B

end

function domainmaps(a::Array{Float64,1}, b::Array{Float64,1}, c::Array{Float64,1})
    # a, b, c are three vectors describing the simulation box
    # a = [a1, a2, a3]
    # b = [b1, b2, b3]
    # c = [c1, c2, c3]

    R2B = [a b c]; # transform the unit cube to the simulation box (a, b, c)
    B2R = inv(R2B); # transform the simulation box (a, b, c) to the unit cube

    return B2R, R2B

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
export simulate

mutable struct OpenableCards
    by_rarity::Vector{Vector{Int}}
    by_set_rarity::Dict{Tuple{Symbol,Int},Vector{Int}}
    fresh_by_rarity::Vector{Vector{Int}}
    fresh_by_set_rarity::Dict{Tuple{Symbol,Int},Vector{Int}}
    sets::Vector{Symbol}
end

function OpenableCards(x::Vector{Vector{Int}}, y::Dict{Tuple{Symbol,Int},Vector{Int}}, sets::Vector{Symbol})
    OpenableCards(x,y,deepcopy(x),deepcopy(y),sets)
end

function update_openables(sets; parameters, db, collection)
    openables = OpenableCards(
            map(i -> setdiff( findall(x ->  x["rarity"] == i && in(x["set"], sets), db),
                              parameters.noopen_cards), 1:4),
            Dict(
                ((s,r), setdiff(indices_set_rarity(db,s,r), parameters.noopen_cards)) for s in sets, r in 1:4
            ), sets
        )

    if parameters.prevent_duplicates
        for i in eachindex(collection)
            if collection[i] >= 4
                remove_fully_collected!(i; openable_cards = openables)
            end
        end
    end

    openables
end

"""
    simulate(nreps::Int, deckinfos; kwargs...)

Main simulation function for a fixed deck order, producing a `SimOutput`.

# Arguments
- `nreps::Int`: number of iterations to run
- `deckinfos`: A vector produced by `deckinfo` representing a sequence of decks to acquire
- `parameters`: A `SimParameters` value for simulation parameters
- `db`: Database of card info
- `sets`: Sets to simulate.
"""
function simulate(nreps::Int, deckinfos; parameters = SimParameters(), db = card_db, sets = sets)
    collection = deepcopy(parameters.starter_cards)
    fresh_openable_cards = update_openables(sets; parameters = parameters, db = db, collection = collection)

    account = AccountState(parameters)

    output = SimOutput(nreps, length(deckinfos), sets, collection)

    if parameters.track_collection_progress
        output.collection_progress = Array{Union{Missing,Int}}(missing, nreps, parameters.max_track_progress_packs, length(collection))
    end

    if parameters.welcome_bundle
        account.bonus_packs[parameters.welcome_bundle_set] += 5
        set = next_set_welcome_bundle!(collection, account, first(deckinfos), sets, fresh_openable_cards)
        account.bonus_packs[set] += 12
    end

    freshaccount = deepcopy(account)

    fresh_fixed_pack_track = deepcopy(parameters.fixed_pack_track)

    pack_contents = [
        zeros(Int, 5),
        zeros(Int, 2),
        zeros(Int, 1),
        zeros(Int, 1)
    ]

    packs_opened = Symbol[]
    bonus_packs_opened = Symbol[]
    icrs_gained = Int[]
    pack_cards_gained = Int[]
    wc_crafts_gained = Int[]

    overall_icrs_gained = Int[]
    overall_pack_cards_gained = Int[]
    overall_wc_crafts_gained = Int[]

    before = NamedTuple{(:card, :needed, :rarity),Tuple{Int64,Float64,Int64}}[]
    after = NamedTuple{(:card, :needed, :rarity),Tuple{Int64,Float64,Int64}}[]

    for rep in 1:nreps
        pitytimer_parameters = rand(parameters.pitytimer_distribution)

        resize!(overall_icrs_gained, 0)
        resize!(overall_pack_cards_gained, 0)
        resize!(overall_wc_crafts_gained, 0)

        if parameters.fixed_starter
            collection = deepcopy(parameters.starter_cards)
        else
            collection = deepcopy(rand(startersequences))
        end
        openable_cards = deepcopy(fresh_openable_cards)

        if parameters.kaladesh_grant
            for i in eachindex(collection)
                collection[i] += kldgrant[i]
                if collection[i] > 4
                    collection[i] = 4
                end
            end
        end
        account = deepcopy(freshaccount)

        if parameters.prevent_duplicates
            for i in eachindex(collection)
                if collection[i] >= 4
                    remove_fully_collected!(i; openable_cards = openable_cards)
                end
            end
        end

        parameters.fixed_pack_track = deepcopy(fresh_fixed_pack_track)

        for deckindex in eachindex(deckinfos)
            resize!(packs_opened, 0)
            resize!(bonus_packs_opened, 0)
            resize!(icrs_gained, 0)
            resize!(pack_cards_gained, 0)
            resize!(wc_crafts_gained, 0)

            current_deck = deckinfos[deckindex]
            before = deck_status(collection = collection, current_deck = current_deck)

            while missing_cards(collection, current_deck)
                if account.vault_pct > 100.0
                        open_vault!(account, output; rep = rep)
                end

                if can_finish_deck(collection, account, current_deck)
                    spend_wildcards!(account = account, collection = collection,
                     current_deck =  current_deck, output = output, rep = rep, deckindex = deckindex,
                     logging_vector = wc_crafts_gained, parameters = parameters, openable_cards = openable_cards)

                    @assert !missing_cards(collection, current_deck)
                    @goto cleanup
                end

                if any(x -> x < Inf, parameters.wc_upgrade_rate)
                    convert_wildcards!(account, parameters)
                end

                set = parameters.nextset(collection, account, current_deck, sets, openable_cards)

                if account.bonus_packs[set] == 0
                    if haskey(parameters.openable_schedule, output.total_earned_packs[rep])
                        openable_cards = update_openables(parameters.openable_schedule[output.total_earned_packs[rep]]; parameters = parameters,
                            collection = collection, db = db)
                    end

                    icr_sequence = nextitem!(parameters.icrs_per_pack)
                    for dist in icr_sequence
                        icr_rarity = rand(dist)
                        process_icr!(collection = collection, openable_cards = openable_cards,
                            account = account, rarity = icr_rarity, parameters = parameters,
                            logging_vector = icrs_gained)
                    end

                    # TODO: Raredrafting simulation?

                    fixedpack = nextitem!(account.fixed_pack_track)

                    if !ismissing(fixedpack)
                        set = fixedpack
                    end

                    if !in(set, openable_cards.sets)
                        set = rand(openable_cards.sets)
                    end

                    push!(packs_opened, set)
                    output.total_earned_packs[rep] += 1
                else
                    account.bonus_packs[set] -= 1
                    push!(bonus_packs_opened, set)
                end

                account.vault_pct += parameters.pack_vault_pct

                timer_result = nextitem!(account.wc_timer)
                if timer_result > 0
                    account.wc_count[timer_result] += 1
                end

                open_pack!(
                  pack_contents;
                  set =  set,
                  collection = collection,
                  pitytimer_parameters = pitytimer_parameters,
                  account =  account,
                  openable_cards = openable_cards,
                  parameters = parameters,
                  logging_vector = pack_cards_gained,
                  output = output,
                  rep = rep
                )

                @label cleanup

                if parameters.track_collection_progress
                    output.collection_progress[rep,output.total_packs[rep],:] = deepcopy(collection)
                end

            end

            for packset in packs_opened
                output.packs_opened[packset][rep,deckindex] += 1
            end

            for packset in bonus_packs_opened
                output.bonus_packs_opened[packset][rep,deckindex] += 1
            end

            append!(overall_icrs_gained, icrs_gained)
            append!(overall_pack_cards_gained, pack_cards_gained)
            append!(overall_wc_crafts_gained, wc_crafts_gained)

            after = deck_status(collection = collection,current_deck =  current_deck)

            for i in eachindex(after)
                rarity = before[i].rarity
                card = before[i].card
                output.card_sources[rep,deckindex,rarity,1] += before[i].total
                output.card_sources[rep,deckindex,rarity,2] += before[i].needed
                output.card_sources[rep,deckindex,rarity,3] += after[i].needed
            end

            cards_gained_counter!(output, wc_crafts_gained, before; rep = rep,
             deckindex = deckindex, column = 4)
            cards_gained_counter!(output, pack_cards_gained, before; rep = rep,
             deckindex = deckindex, column = 5)
            cards_gained_counter!(output, icrs_gained, before; rep = rep,
             deckindex = deckindex, column = 6)
            cards_gained_counter!(output, overall_wc_crafts_gained, before; rep = rep,
             deckindex = deckindex, column = 7)
            cards_gained_counter!(output, overall_pack_cards_gained, before; rep = rep,
             deckindex = deckindex, column = 8)
            cards_gained_counter!(output, overall_icrs_gained, before; rep = rep,
              deckindex = deckindex, column = 9)
        end
        output.ending_collection[rep,:] = collection
    end

    maxpacks = maximum(output.total_packs)
    if length(output.collection_progress) > 0 && maxpacks < parameters.max_track_progress_packs
        output.collection_progress = output.collection_progress[:,1:maxpacks,:]
    end

    output
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
using StatsBase, Distributions, Plots; pyplot()

binomialRV(n,p) = sum(rand(n) .< p)

p, n, N = 0.5, 10, 10^6

bDist = Binomial(n,p)
xGrid = 0:n
bPmf = [pdf(bDist,i) for i in xGrid]
data = [binomialRV(n,p) for _ in 1:N]
pmfEst = counts(data,0:n)/N

plot( xGrid, pmfEst,
	line=:stem, marker=:circle,
	c=:blue, ms=10, msw=0, lw=4, label="MC estimate")
plot!( xGrid, bPmf,
	 line=:stem, marker=:xcross, c=:red,
	 ms=6, msw=0, lw=2, label="PMF", xticks=(0:1:10),
	 ylims=(0,0.3), xlabel="x", ylabel="Probability")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using BetterBanner
using Test

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


text:
# Code for computing with domains.

############################
# The union of two domains
############################

# DD can be any collection type: a Tuple, or a Vector
# Eventually it will likely be forced to be an AbstractSet
#
"""
A `UnionDomain` represents the union of a set of domains.
"""
struct UnionDomain{DD,T} <: Domain{T}
    domains  ::  DD
end

"""
The `UnionDomain` constructor can be invoked in one of three ways:
- with a list of arguments: UnionDomain(d1, d2)
- with a single domain: UnionDomain(d::Domain)
- or with any iterable list of domains: UnionDomain(domains)
"""
UnionDomain(domains...) = UnionDomain(domains)
UnionDomain(d::Domain) = UnionDomain((d,))

UnionDomain(domains) = _UnionDomain(expand_union(promote_domains(domains)...))
_UnionDomain(domains) = UnionDomain{typeof(domains),eltype(first(domains))}(domains)

"Return all arguments in a tuple of domains, with each `UnionDomain` replaced by its elements."
expand_union() = ()
expand_union(domain, domains...) = (domain, expand_union(domains...)...)
expand_union(domain::UnionDomain, domains...) = (elements(domain)..., expand_union(domains...)...)

convert(::Type{Domain}, v::AbstractVector{<:Domain}) = UnionDomain(v)
convert(::Type{Domain}, v::AbstractSet{<:Domain}) = UnionDomain(v)

convert(::Type{Domain{T}}, d::UnionDomain{S}) where {S,T} =
    UnionDomain(convert_domain.(T, elements(d)))

Domain(v::AbstractVector) = convert(Domain, v)
Domain(v::AbstractSet) = convert(Domain, v)

elements(d::UnionDomain) = d.domains

hash(d::UnionDomain, h::UInt) = hash(Set(d.domains), h)

union(domains::Domain...) = UnionDomain(domains...)

# Catch one specific case
union(d1::Domain, d2::Domain)  = d1 == d2 ? d1 : UnionDomain(d1, d2)

# The union of domains corresponds to a logical OR of their characteristic functions
indomain(x, d::UnionDomain) = mapreduce(d->in(x, d), |, elements(d))
approx_indomain(x, d::UnionDomain, tol = default_tolerance(d)) = mapreduce(d->approx_in(x, d, tol), |, elements(d))

point_in_domain(d::UnionDomain) = point_in_domain(element(d,1))

==(a::UnionDomain, b::UnionDomain) = Set(elements(a)) == Set(elements(b))

isempty(d::UnionDomain) = all(isempty, d.domains)

function show(io::IO, d::UnionDomain)
    print(io, "a union of $(numelements(d)) domains:\n")
    for (i,e) in enumerate(elements(d))
        print(io, "\t$i.\t: ", e, "\n")
    end
end

## ualgebra

for op in (:+, :-, :*, :/)
    @eval begin
        $op(domain::UnionDomain, x::Number) = UnionDomain(broadcast($op, elements(domain), x))
        $op(x::Number, domain::UnionDomain) = UnionDomain(broadcast($op, x, elements(domain)))
    end
end

\(x::Number, domain::UnionDomain) = UnionDomain(broadcast(\, x, elements(domain)))

for (op, mop) in ((:minimum, :min), (:maximum, :max), (:infimum, :min), (:supremum, :max))
    @eval $op(d::UnionDomain) = mapreduce($op, $mop, elements(d))
end

setdiff(d1::UnionDomain, d2::UnionDomain) = UnionDomain(setdiff.(elements(d1), Ref(d2)))

function setdiff(d1::UnionDomain, d2::Domain)
    s = Set(elements(d1))
    # check if any element is in d1 and just remove
    s2 = Set(setdiff(s, tuple(d2)))
    s2 ≠ s && return UnionDomain(s2)

    UnionDomain(setdiff.(elements(d1), Ref(d2)))
end

function setdiff(d1::Domain, d2::UnionDomain)
    ret = d1
    for d in elements(d2)
        ret = setdiff(ret, d)
    end
    ret
end

# use \ as a synomym for setdiff, in the context of domains (though, generically,
# \ means left division in Julia)
\(d1::Domain, d2::Domain) = setdiff(d1, d2)

##############################
# The intersection of domains
##############################

"""
An `IntersectionDomain` represents the intersection of a set of domains.
"""
struct IntersectionDomain{DD,T} <: Domain{T}
    domains ::  DD
end

"""
The `IntersectionDomain` constructor can be invoked in one of three ways:
- with a list of arguments: IntersectionDomain(d1, d2)
- with a single domain: IntersectionDomain(d::Domain)
- or with any iterable list of domains: IntersectionDomain(domains)
"""
IntersectionDomain(domains...) = IntersectionDomain(domains)
IntersectionDomain(d::Domain) = IntersectionDomain((d,))

IntersectionDomain(domains) = _IntersectionDomain(expand_intersection(promote_domains(domains)...))
_IntersectionDomain(domains) = IntersectionDomain{typeof(domains),eltype(first(domains))}(domains)

"Return all arguments in a tuple of domains, with each `UnionDomain` replaced by its elements."
expand_intersection() = ()
expand_intersection(domain, domains...) = (domain, expand_intersection(domains...)...)
expand_intersection(domain::IntersectionDomain, domains...) = (elements(domain)..., expand_intersection(domains...)...)

elements(d::IntersectionDomain) = d.domains

intersect(d1::Domain, d2::Domain) = d1 == d2 ? d1 : IntersectionDomain(d1, d2)
function intersect(d1::UnionDomain, d2::UnionDomain)
    d1 == d2 && return d1
    union(intersect.(Ref(d1), elements(d2))...)
end
intersect(d1::UnionDomain, d2::Domain) = union(intersect.(d1.domains, Ref(d2))...)
intersect(d1::Domain, d2::UnionDomain) = union(intersect.(Ref(d1), d2.domains)...)

# The intersection of domains corresponds to a logical AND of their characteristic functions
indomain(x, d::IntersectionDomain) = mapreduce(d->in(x, d), &, elements(d))

(&)(d1::Domain, d2::Domain) = intersect(d1,d2)

function intersect(d1::ProductDomain, d2::ProductDomain)
    if numelements(d1) == numelements(d2)
        ProductDomain([intersect(element(d1,i), element(d2,i)) for i in 1:numelements(d1)]...)
    else
        IntersectionDomain(d1, d2)
    end
end

convert(::Type{Domain{T}}, d::IntersectionDomain{S}) where {S,T} =
    IntersectionDomain(convert_domain.(T, elements(d)))

function show(io::IO, d::IntersectionDomain)
    print(io, "the intersection of $(numelements(d)) domains:\n")
    print(io, "\t1.\t: ", element(d,1), "\n")
    print(io, "\t2.\t: ", element(d,2), "\n")
end

#########################################
### The difference between two domains
#########################################

struct DifferenceDomain{D1,D2,T} <: Domain{T}
    d1    ::  D1
    d2    ::  D2
end

DifferenceDomain(d1, d2) = _DifferenceDomain(promote_domain(d1, d2)...)
_DifferenceDomain(d1, d2) = DifferenceDomain{typeof(d1),typeof(d2),eltype(d1)}(d1, d2)

convert(::Type{Domain{T}}, d::DifferenceDomain{S}) where {S,T} =
    DifferenceDomain(convert_domain(T, d.d1), convert_domain(T, d.d2))

function setdiff(d1::Domain{T}, d2::Domain) where T
    d1 == d2 && return EmptySpace{T}()
    DifferenceDomain(d1, d2)
end

setdiff(d1::Domain, d2) = DifferenceDomain(d1, d2)
setdiff(d1, d2::Domain) = DifferenceDomain(d1, d2)

# The difference between two domains corresponds to a logical AND NOT of their characteristic functions
indomain(x, d::DifferenceDomain) = in(x, d.d1) && (~in(x, d.d2))

\(d1::Domain, d2) = setdiff(d1, d2)

function show(io::IO, d::DifferenceDomain)
    print(io, "the difference of 2 domains:\n")
    print(io, "\t1.\t: ", d.d1, "\n")
    print(io, "\t2.\t: ", d.d2, "\n")
end


text:
using Plots, LaTeXStrings, Statistics

function Deposition(;len, tot_time, time_steps)
    Time = ceil.(Int, exp.(2:(tot_time-2)/(time_steps):tot_time))
    surf = [0 for i=1:len]
    VarList = [0.0 for i=1:time_steps]
    for n in 2:time_steps+1
        randsurf = rand(1:len,(Time[n]-Time[n-1]))
        for i in randsurf
            MAX = FindMax(surf, i, len)
            surf[i] = MAX
        end
        VarList[n-1] = std(surf)
    end
    return VarList
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function FindMax(surface, index_,L_surf)
    i1 , i2 = sides(index_, L_surf)
    maxlen = max(surface[i1],surface[index_] + 1,surface[i2])
    return maxlen
end

function Linear_fit(;Time, VarList, time_steps)
    A = [hcat(log.(Time[1:time_steps])) reshape(ones(time_steps), time_steps, 1)]
    b = reshape(log.(VarList[1:time_steps]), time_steps, 1)
    line = (A \ b)
    return line
end

# function FindLine(tot_steps, Time, VarList)
#     retline = [0.0 0.0]
#     for i in 5:tot_steps
#         Paraline = Dict(
#                         :Time => Time,
#                         :VarList => VarList,
#                         :time_steps => i
#                             )
#         Line = Linear_fit(;Paraline...)
#         absl = abs(Line[1]*Time[i] + Line[2] - VarList[i])
#         standev = std(Line[1].*Time .+ Line[2] .- VarList)
#         if absl > 0.2 && standev > 1 && retline == [0.0 0.0]
#             retline = Line
#         end
#         if absl > 1 && standev > 3
#             return retline, i
#         end
#     end
# end

theme(:dark)
gr()

iternum = 1000
Parameters = Dict(:len => 200,
                    :tot_time => 13,
                        :time_steps => 50)
allVar = [ [0.0 for i in 1:Parameters[:time_steps]] for j = 1:iternum]
meanVar = [0.0 for i in 1:Parameters[:time_steps]]
vars = [0.0 for i in 1:Parameters[:time_steps]]
for i in 1:iternum
    VarList = Deposition(;Parameters...)
    allVar[i] = VarList
    meanVar += VarList
    print("\r$i")
end
meanVar /= iternum
for i in 1:Parameters[:time_steps]
    vars[i] = std(log.(hcat(allVar...))[i,:])
end

Time = exp.(0:(Parameters[:tot_time])/(Parameters[:time_steps]-1):Parameters[:tot_time])
scatter(log.(Time[3:end]), log.(meanVar[3:end]),
    # xlims = (1, Parameters[:tot_time]),
    xlabel= L"Log\ Time",
    ylabel= L"Log\ W_{(t)}",
    title= L"Log-Log\ Plot\ ~W_{(t)}-Time~\ (L = %$(Parameters[:len]))",
    label = L"Data\ point",
    yerror = vars,
    legend = 150)

Paraline = Dict(
                :Time => Time[3:end],
                :VarList => meanVar[3:end],
                :time_steps => 13)
Line = Linear_fit(;Paraline...)
X = log.(Time[3]:Time[30])
Y = X .* Line[1] .+ Line[2]

plot!(X,Y,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))")

savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\Ballistic-Deposition\\Deposition\\Ballistic-Deposition-Kardar-Parisi-Zhang\\Fig\\W-t(L=$(Parameters[:len])).png")

max_stars_repo_name: shahmari/Ballistic-Deposition

max_issues_repo_name: shahmari/Ballistic-Deposition

max_forks_repo_name: shahmari/Ballistic-Deposition


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestModel

using Base.Test
using MXNet

function test_feedforward()
  info("Model::FeedForward::constructor")
  let x = @mx.var x
    m = mx.FeedForward(x)
    @test m.arch === x
    @test length(m.ctx) == 1
  end

  info("Model::FeedForward::constructor::keyword context")
  let x = @mx.var x
    m = mx.FeedForward(x, context = mx.cpu())
    @test m.arch === x
    @test length(m.ctx) == 1
  end

  let x = @mx.var x
    m = mx.FeedForward(x, context = [mx.cpu(), mx.cpu(1)])
    @test m.arch === x
    @test length(m.ctx) == 2
  end
end

@testset "Model Test" begin
  test_feedforward()
end

end  # module TestModel

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
using DelimitedFiles
using Printf

function count_trees(sf, dh, dv)
    nrows = size(sf)[1]
    ncols = length(sf[1,1][:])

    pos = 1
    ntree = 0

    for i in 1+dv:dv:nrows
        # println(i)
        # println(sf[i,1][:])

        pos = pos + dh
        if pos > ncols
            pos = pos - ncols
        end

        if sf[i,1][pos] == '#'
            ntree += 1
        end
    end
    return ntree
end

slope = readdlm("./input.dat")

println("Part a:\t", count_trees(slope, 3, 1))

a = count_trees(slope, 1, 1)
b = count_trees(slope, 3, 1)
c = count_trees(slope, 5, 1)
d = count_trees(slope, 7, 1)
e = count_trees(slope, 1, 2)

println("===================")
println("R1, D1\t", a)
println("R3, D1\t", b)
println("R5, D1\t", c)
println("R7, D1\t", d)
println("R1, D2\t", e)
println("===================")
println("Product: ", a*b*c*d*e)