text:
include("Graph.jl")

function CanonicalOrdering(G, Emb)
	processedNeighbor = zeros(Int, G.n)
	status = ones(Int, G.n)
	# status : 0 -> Processed
	#          1 -> Unprocessed
	#          2 -> One Neighbor Processed
	#          10 -> Ready To Be Processed
	Q = zeros(Int, 2*G.n)
	Order = Array{Int, 1}()
	front = 1
	rear = 0

	target(eid, u) = begin
		return (G.E[eid][2] != u) ? G.E[eid][2] : G.E[eid][3]
	end
	PushQ(x) = begin
		rear += 1
		Q[rear] = x
	end
	PopQ() = begin
		tmp = Q[front]
		front += 1
		return tmp
	end

	firstVertex = 1
	secondVertex = target(Emb[firstVertex][1], firstVertex)

	PushQ(firstVertex)
	status[firstVertex] = 10
	PushQ(secondVertex)
	status[secondVertex] = 10

	while front <= rear
		u = PopQ()
		if ((status[u] != 10) && (u != secondVertex))
			continue
		end
		for i = 1 : size(Emb[u], 1)
			priorEdge = (i == 1) ? Emb[u][size(Emb[u], 1)] : Emb[u][i-1]
			nextEdge = (i == size(Emb[u], 1)) ? Emb[u][1] : Emb[u][i+1]
			v = target(Emb[u][i], u)
			priorVertex = target(priorEdge, u)
			nextVertex = target(nextEdge, u)

			if status[v] == 1
				status[v] = 2
				processedNeighbor[v] = u
			elseif status[v] == 2
				x = processedNeighbor[v]
				if (((nextVertex == x) && !((firstVertex == u) && (secondVertex == x))) || ((priorVertex == x) && !((firstVertex == x) && (secondVertex == u))))
					status[v] = 10
				else
					status[v] = 10+1
				end
			elseif status[v] > 2
				processedBefore = (status[priorVertex] == 0)
				processedAfter = (status[nextVertex] == 0)
				if ((processedBefore == false) && (processedAfter == false))
					status[v] += 1
				end
				if ((processedBefore == true) && (processedAfter == true))
					status[v] -= 1
				end
			end
			if status[v] == 10
				PushQ(v)
			end
		end
		status[u] = 0
		push!(Order, u)
	end

	return Order
end

function PolylineDrawing(G, m, Order, Emb)
	edges = Array{Tuple{Int, Int}, 1}()
	A = zeros(Int32, G.n, G.n)
	newID = zeros(Int, G.n)
	foreach(i -> newID[Order[i]] = i, 1 : G.n)
	emb = Array{Array{Int, 1}, 1}(undef, G.n)
	for i = 1 : G.n
		emb[i] = []
		for j = 1 : size(Emb[Order[i]], 1)
			push!(emb[i], Emb[Order[i]][j])
		end
	end
	for (ID, u, v, w) in G.E
		nu = newID[u]
		nv = newID[v]
		A[nu, nv] = ID
		A[nv, nu] = ID
		(nu > nv) ? push!(edges, (nv, nu)) : push!(edges, (nu, nv))
	end
	InList = Array{Array{Int, 1}, 1}(undef, G.n)
	foreach(i -> InList[i] = [], 1 : G.n)
	OutList = Array{Array{Int, 1}, 1}(undef, G.n)
	foreach(i -> OutList[i] = [], 1 : G.n)

	target(u, eid) = (edges[eid][1] == u) ? edges[eid][2] : edges[eid][1]
	next(u, p) = (p == 1) ? size(emb[u], 1) : p-1
	prior(u, p) = (p == size(emb[u], 1)) ? 1 : p+1

	# Step 1 : calculate in and out list
	## calculate OutList[1]
	pos = 1
	while target(1, emb[1][pos]) != 2
		pos += 1
	end
	pos = next(1, pos)
	for i = 1 : size(emb[1], 1)
		if emb[1][pos] <= m
			push!(OutList[1], emb[1][pos])
		end
		pos = next(1, pos)
	end
	## calculate InList[2] and OutList[2]
	pos = 1
	while target(2, emb[2][pos]) != 1
		pos += 1
	end
	(emb[2][pos] <= m) ? push!(InList[2], emb[2][pos]) : nothing
	pos = next(2, pos)
	for i = 1 : size(emb[2], 1)-1
		if emb[2][pos] <= m
			push!(OutList[2], emb[2][pos])
		end
		pos = next(2, pos)
	end
	## calculate InList and OutList for node 3-n
	PRED = zeros(Int, G.n)
	SUCC = zeros(Int, G.n)
	SUCC[1] = 2
	PRED[2] = 1

	Val = zeros(Int, G.n)
	LL = zeros(Int, G.n)
	RR = zeros(Int, G.n)
	faceStart = 1
	sp = 2
	Val[1], Val[2] = 1, 2
	LL[1], LL[2] = 0, 1
	RR[1], RR[2] = 2, 0

	addBetween(lv, rv, x) = begin
		sp += 1
		LL[sp], Val[sp], RR[sp] = lv, x, rv
		RR[lv], LL[rv] = sp, sp
	end

	for k = 3 : G.n
		w = Array{Int, 1}()
		p = faceStart
		lv = -1
		rv = -1
		while p != 0
			if A[k, Val[p]] != 0
				push!(w, Val[p])
				(lv == -1) ? (lv = p) : nothing
				(lv != -1) ? (rv = p) : nothing
			end
			p = RR[p]
		end
		wn = size(w, 1)
		PRED[k] = w[1]
		SUCC[k] = w[wn]
		p = 1
		while target(k, emb[k][p]) != w[1]
			p += 1
		end
		q = next(k, p)
		while target(k, emb[k][q]) != w[wn]
			(emb[k][q] <= m) ? push!(OutList[k], emb[k][q]) : nothing
			q = next(k, q)
		end
		while true
			(emb[k][p] <= m) ? push!(InList[k], emb[k][p]) : nothing
			if target(k, emb[k][p]) == w[wn]
				break
			end
			p = prior(k, p)
		end
		addBetween(lv, rv, k)
	end

	indeg(x) = size(InList[x], 1)
	outdeg(x) = size(OutList[x], 1)

	# Step 2 : calculate outpoint and inpoint for each edge
	edgeinx = zeros(Int, m)
	edgeiny = zeros(Int, m)
	edgeoutx = zeros(Int, m)
	edgeouty = zeros(Int, m)
	outl = zeros(Int, G.n)
	outr = zeros(Int, G.n)
	inl = zeros(Int, G.n)
	inr = zeros(Int, G.n)

	isLinked(u, v) = ((u != 0) && (v != 0) && (A[u, v] <= m))

	for v = 1 : G.n
		## arrange outpoint
		outvPLUS = div(outdeg(v), 2)
		outvMINUS = div(outdeg(v)-1, 2)
		deltal = 1
		deltar = 1
		if indeg(v) >= 2
			outl[v], outr[v] = outvMINUS, outvPLUS
		elseif indeg(v) == 1
			if isLinked(v, SUCC[v])
				outl[v], outr[v] = outvPLUS, outvMINUS
				deltal = 0
			elseif isLinked(v, PRED[v])
				outl[v], outr[v] = outvMINUS, outvPLUS
				deltar = 0
			else
				outl[v], outr[v] = outvMINUS, outvPLUS
			end
		else
			outl[v], outr[v] = outvMINUS, outvPLUS
			deltal, deltar = 0, 0
		end
		tx, ty = -outl[v], deltal
		if outdeg(v) > 0
			for i = 1 : outl[v]
				edgeoutx[OutList[v][i]] = tx
				edgeouty[OutList[v][i]] = ty
				tx += 1
				ty += 1
			end
			edgeoutx[OutList[v][outl[v]+1]] = 0
			edgeouty[OutList[v][outl[v]+1]] = max(outl[v]+deltal-1, outr[v]+deltar-1, 0)
			tx, ty = 1, outr[v]+deltar-1
			for i = outl[v]+2 : outdeg(v)
				edgeoutx[OutList[v][i]] = tx
				edgeouty[OutList[v][i]] = ty
				tx += 1
				ty -= 1
			end
		end

		## arrange inpoint
		inl[v] = max(0, div(indeg(v)-3, 2))
		inr[v] = max(0, div(indeg(v)-2, 2))
		if indeg(v) >= 3
			edgeinx[InList[v][1]] = -inl[v]
			edgeiny[InList[v][1]] = 0
			tx, ty = -inl[v], -1
			for i = 2 : inl[v]+1
				edgeinx[InList[v][i]] = tx
				edgeiny[InList[v][i]] = ty
				tx += 1
				ty -= 1
			end
			edgeinx[InList[v][inl[v]+2]] = 0
			edgeiny[InList[v][inl[v]+2]] = -inr[v]
			tx, ty = 1, -inr[v]
			for i = inl[v]+3 : indeg(v)-1
				edgeinx[InList[v][i]] = tx
				edgeiny[InList[v][i]] = ty
				tx += 1
				ty += 1
			end
			edgeinx[InList[v][indeg(v)]] = inr[v]
			edgeiny[InList[v][indeg(v)]] = 0
		end
	end

	# Step 3 : arrange all nodes
	fill!(Val, 0)
	fill!(LL, 0)
	fill!(RR, 0)
	Father = zeros(Int, G.n)
	faceStart = 1
	sp = 2
	Val[1], Val[2] = 1, 2
	LL[1], LL[2] = 0, 1
	RR[1], RR[2] = 2, 0

	X = zeros(Int, G.n)
	Y = zeros(Int, G.n)
	X[1] = outl[1]
	X[2] = outr[1] + outl[2] + 1

	for k = 3 : G.n
		w = Array{Int, 1}()
		p = faceStart
		lv = -1
		rv = -1
		while p != 0
			if A[k, Val[p]] != 0
				push!(w, Val[p])
				(lv == -1) ? (lv = p) : nothing
				(lv != -1) ? (rv = p) : nothing
			end
			p = RR[p]
		end
		wn = size(w, 1)
		maxy = Y[w[1]]
		foreach(i -> maxy = max(maxy, Y[w[i]]), 2 : wn)
		Y[k] = inr[k] + maxy + 1
		sumx = 0
		for i = 2 : wn
			sumx += X[w[i]]
			X[w[i]] = sumx
		end
		## calculate nextLeftK and nextRightK
		nextLeftK, nextRightK = 0, 0
		pp = -1
		if outdeg(PRED[k]) > 0
			for i = 1 : outdeg(PRED[k])
				if target(PRED[k], OutList[PRED[k]][i]) == k
					pp = OutList[PRED[k]][i]
					break
				end
			end
			nextRightK = (pp == -1) ? (outr[PRED[k]]+1) : (edgeoutx[pp])
		end
		pp = -1
		if outdeg(SUCC[k]) > 0
			for i = 1 : outdeg(SUCC[k])
				if target(SUCC[k], OutList[SUCC[k]][i]) == k
					pp = OutList[SUCC[k]][i]
					break
				end
			end
			nextLeftK = (pp == -1) ? (-outl[SUCC[k]]-1) : (edgeoutx[pp])
		end
		dxl = nextRightK
		dxr = nextLeftK

		if indeg(k) >= 3
			dxt = edgeoutx[InList[k][inl[k]+2]]
			ct = target(k, InList[k][inl[k]+2])
			X[k] = max(X[ct] + dxt, dxl + outl[k])
			dta = X[k] - (X[ct] + dxt)
			X[SUCC[k]] = max(X[SUCC[k]] + dta - X[k], outr[k] - dxr)
			posT = 1
			while w[posT] != ct
				posT += 1
			end
			foreach(i -> X[w[i]] -= X[k], 2 : posT-1)
			foreach(i -> X[w[i]] += (dta - X[k]), posT : wn-1)
			foreach(i -> Father[w[i]] = k, 2 : wn-1)
		else
			X[k] = outl[k] + dxl
			(indeg(k) == 2 && isLinked(k, PRED[k]) == false) ? (X[k] = max(X[k], X[target(k, InList[k][1])] + edgeoutx[InList[k][1]])) : nothing
			if (indeg(k) == 2) && (!isLinked(k, SUCC[k])) && (X[k] > (X[target(k, InList[k][2])] + edgeoutx[InList[k][2]]))
				tmpdx = X[k] - (X[target(k, InList[k][2])] + edgeoutx[InList[k][2]])
				tmppt = 1
				while (w[tmppt] != target(k, InList[k][2]))
					tmppt += 1
				end
				foreach(i -> X[w[i]] += tmpdx, tmppt : wn)
			end
			X[SUCC[k]] = max(outr[k]-dxr, X[SUCC[k]]-X[k])
			foreach(i -> X[w[i]] -= X[k], 2 : wn-1)
			foreach(i -> Father[w[i]] = k, 2 : wn-1)
		end
		addBetween(lv, rv, k)
	end

	X2 = zeros(Int, G.n)
	finalw = Array{Int, 1}()
	p2 = faceStart
	while p2 != 0
		push!(finalw, Val[p2])
		p2 = RR[p2]
	end
	visited = zeros(Bool, G.n)
	sumx2 = 0
	for i = 1 : size(finalw, 1)
		sumx2 += X[finalw[i]]
		X2[finalw[i]] = sumx2
		visited[finalw[i]] = true
	end

	getX2(u) = begin
		if visited[u] == false
			X2[u] = X[u] + getX2(Father[u])
			visited[u] = true
		end
		return X2[u]
	end

	foreach(i -> X2[i] = getX2(i), 1 : G.n)

	# Step 4 : calculate path
	realX = zeros(Int, G.n)
	realY = zeros(Int, G.n)
	PATH = Array{Array{Tuple{Int, Int}, 1}, 1}(undef, m)
	deg = zeros(Int, G.n)
	foreach(i -> deg[i] = indeg(i) + outdeg(i), 1 : G.n)

	cross(x1, y1, x2, y2) = (x1*y2 - x2*y1)

	addToPATH(eid, x, y) = begin
		tn = size(PATH[eid], 1)
		if tn == 0
			push!(PATH[eid], (x, y))
		elseif (x != PATH[eid][tn][1]) || (y != PATH[eid][tn][2])
			if tn > 1
				if cross(PATH[eid][tn-1][1] - PATH[eid][tn][1], PATH[eid][tn-1][2] - PATH[eid][tn][2], x - PATH[eid][tn][1], y - PATH[eid][tn][2]) == 0
					PATH[eid][tn] = (x, y)
				else
					push!(PATH[eid], (x, y))
				end
			else
				push!(PATH[eid], (x, y))
			end
		end
	end

	for i = 1 : m
		(u, v) = edges[i]
		PATH[i] = []
		addToPATH(i, X2[u], Y[u])
		tx1, ty1 = X2[u]+edgeoutx[i], Y[u]+edgeouty[i]
		addToPATH(i, tx1, ty1)
		tx2, ty2 = X2[v]+edgeinx[i], Y[v]+edgeiny[i]
		((tx1 != tx2) && (ty1 != ty2)) ? addToPATH(i, tx1, ty2) : nothing
		addToPATH(i, tx2, ty2)
		addToPATH(i, X2[v], Y[v])
	end

	## dealing with one degree node
	hv = zeros(Bool, G.n)
	Snode = Array{Int, 1}()
	Sedge = Array{Int, 1}()
	Spath = Array{Tuple{Int, Int}, 1}()

	getEdgeID(u, pred) = begin
		for x in InList[u]
			if target(u, x) != pred
				return x, target(u, x)
			end
		end
		for x in OutList[u]
			if target(u, x) != pred
				return x, target(u, x)
			end
		end
		return nothing
	end

	addToSpath(u, eu) = begin
		if u < target(u, eu)
			for i = 1 : size(PATH[eu], 1)
				push!(Spath, PATH[eu][i])
			end
		else
			for i = size(PATH[eu], 1) : -1 : 1
				push!(Spath, PATH[eu][i])
			end
		end
	end

	getDir(x0, y0, x1, y1) = begin
		if (x0 == x1) && (y0 > y1)
			return 0, -1
		elseif (x0 == x1) && (y0 < y1)
			return 0, 1
		elseif (x0 > x1) && (y0 == y1)
			return -1, 0
		else
			return 1, 0
		end
	end

	DEC() = begin
		top = size(Spath, 1)
		while (top > 1) && (Spath[top] == Spath[top-1])
			pop!(Spath)
			top -= 1
		end
		if top > 1
			dx, dy = getDir(Spath[top][1], Spath[top][2], Spath[top-1][1], Spath[top-1][2])
			Spath[top] = (Spath[top][1]+dx, Spath[top][2]+dy)
		end
	end

	for i = 1 : G.n
		if (!hv[i]) && (deg[i] == 1)
			while size(Snode, 1) > 0
				pop!(Snode)
			end
			while size(Sedge, 1) > 0
				pop!(Sedge)
			end
			while size(Spath, 1) > 0
				pop!(Spath)
			end
			u = i
			eu, v = getEdgeID(u, -1)
			hv[v] = true
			push!(Snode, u)
			push!(Sedge, eu)
			addToSpath(u, eu)
			while deg[v] == 2
				predu = last(Snode)
				u = v
				eu, v = getEdgeID(u, predu)
				hv[v] = true
				push!(Snode, u)
				push!(Sedge, eu)
				addToSpath(u, eu)
			end
			tp = size(Spath, 1)
			# arrange the first node
			u = last(Snode)
			eu = last(Sedge)
			PATH[eu] = []
			push!(PATH[eu], last(Spath))
			if u < target(u, eu)
				if (edgeinx[eu] == 0) && (edgeiny[eu] == 0)
					DEC()
				else
					pop!(Spath)
					push!(Spath, (X2[target(u, eu)]+edgeinx[eu], Y[target(u, eu)]+edgeiny[eu]))
				end
			else
				if (edgeoutx[eu] == 0) && (edgeouty[eu] == 0)
					DEC()
				else
					pop!(Spath)
					push!(Spath, (X2[target(u, eu)]+edgeoutx[eu], Y[target(u, eu)]+edgeouty[eu]))
				end
			end
			push!(PATH[eu], last(Spath))
			X2[u] = last(Spath)[1]
			Y[u] = last(Spath)[2]
			pop!(Snode)
			pop!(Sedge)
			# arrange rest of the nodes
			while size(Snode, 1) > 0
				u = last(Snode)
				eu = last(Sedge)
				PATH[eu] = []
				push!(PATH[eu], last(Spath))
				DEC()
				push!(PATH[eu], last(Spath))
				X2[u] = last(Spath)[1]
				Y[u] = last(Spath)[2]
				pop!(Snode)
				pop!(Sedge)
			end
		end
	end

	## dealing with two degree node
	getTwoPaths(u) = begin
		ids = Vector{Int}()
		foreach(x -> push!(ids, x), InList[u])
		foreach(x -> push!(ids, x), OutList[u])
		for i in ids
			(PATH[i][1] != (X2[u], Y[u])) ? reverse!(PATH[i]) : nothing
		end
		return ids;
	end

	isBend(p1, p2, p3) = !(((p1[1] == p2[1]) && (p1[1] == p3[1])) || ((p1[2] == p2[2]) && (p1[2] == p3[2])))

	closePoint(u, v, eid) = begin
		tp = (u < v) ? (edgeoutx[eid], edgeouty[eid]) : (edgeinx[eid], edgeiny[eid])
		if (tp != (0, 0))
			return X2[u]+tp[1], Y[u]+tp[2]
		end
		dx, dy = getDir(X2[u], Y[u], X2[v], Y[v])
		return X2[u]+dx, Y[u]+dy
	end

	for i = 1 : G.n
		if (!hv[i]) && (deg[i] == 2)
			ids = getTwoPaths(i)
			if isBend(PATH[ids[1]][2], PATH[ids[1]][1], PATH[ids[2]][2])
				continue
			end
			if (size(PATH[ids[1]], 1) > 2)
				PATH[ids[1]] = PATH[ids[1]][2 : size(PATH[ids[1]], 1)]
				X2[i], Y[i] = PATH[ids[1]][1]
				PATH[ids[2]][1] = PATH[ids[1]][1]
			else
				X2[i], Y[i] = closePoint(target(i, ids[1]), i, ids[1])
				PATH[ids[1]][1] = (X2[i], Y[i])
				PATH[ids[2]][1] = (X2[i], Y[i])
			end
		end
	end

	# delete unused rows and columns
	Xs = Array{Int, 1}()
	Ys = Array{Int, 1}()
	NewX = Dict{Int, Int}()
	NewY = Dict{Int, Int}()
	for xx in X2
		push!(Xs, xx)
	end
	for yy in Y
		push!(Ys, yy)
	end
	for i = 1 : size(PATH, 1)
		for (xx, yy) in PATH[i]
			push!(Xs, xx)
			push!(Ys, yy)
		end
	end
	sort!(Xs)
	sort!(Ys)
	totalX = 1
	totalY = 1
	for xx in Xs
		if !haskey(NewX, xx)
			NewX[xx] = totalX
			totalX += 1
		end
	end
	for yy in Ys
		if !haskey(NewY, yy)
			NewY[yy] = totalY
			totalY += 1
		end
	end
	realPATH = Array{Array{Tuple{Int, Int}, 1}, 1}(undef, m)
	for i = 1 : size(PATH, 1)
		realPATH[i] = []
		for (xx, yy) in PATH[i]
			push!(realPATH[i], (NewX[xx], NewY[yy]))
		end
	end

	for i = 1 : G.n
		realX[Order[i]] = NewX[X2[i]]
		realY[Order[i]] = NewY[Y[i]]
	end
	return realX, realY, realPATH, totalX-1, totalY-1
end


text:
# helper functions for spherical harmonics spectra

"""
    SpectralSphericalMesh

Struct with mesh information.
"""
mutable struct SpectralSphericalMesh
    # grid info
    num_fourier::Int64
    num_spherical::Int64
    nλ::Int64
    nθ::Int64
    nd::Int64
    Δλ::Float64
    qwg::Array{Float64, 3}
    qnm::Array{Float64, 3}   # n,m coordinates
    wave_numbers::Array{Int64, 2}

    # variables
    var_grid::Array{Float64, 3}
    var_fourier::Array{ComplexF64, 3}
    var_spherical::Array{ComplexF64, 4}
    var_spectrum::Array{Float64, 3}

end
function SpectralSphericalMesh(nθ::Int64, nd::Int64)
    nλ = 2nθ
    Δλ = 2π / nλ

    num_fourier = Int64(floor((2 * nθ - 1) / 3)) # number of truncated zonal wavenumbers (m): minimum truncation given nθ - e.g.: nlat = 32 -> T21 (can change manually for more a severe truncation)
    num_spherical = Int64(num_fourier + 1) # number of total wavenumbers (n)

    radius = Float64(6371000)
    wave_numbers = compute_wave_numbers(num_fourier, num_spherical)

    qwg = zeros(Float64, num_fourier + 1, num_spherical + 1, nθ)
    qnm = zeros(Float64, num_fourier + 1, num_spherical + 2, nθ)

    var_fourier = zeros(Complex{Float64}, nλ, nθ, nd)
    var_grid = zeros(Float64, nλ, nθ, nd)
    nθ_half = div(nθ, 2)
    var_spherical =
        zeros(Complex{Float64}, num_fourier + 1, num_spherical + 1, nd, nθ_half)
    var_spectrum = zeros(Float64, num_fourier + 1, num_spherical + 1, nd)

    SpectralSphericalMesh(
        num_fourier,
        num_spherical,
        nλ,
        nθ,
        nd,
        Δλ,
        qwg,
        qnm,
        wave_numbers,
        var_grid,
        var_fourier,
        var_spherical,
        var_spectrum,
    )
end

"""
    compute_legendre!(num_fourier, num_spherical, sinθ, nθ)

Normalized associated Legendre polynomials, P_{m,l} = qnm

# Arguments:
- num_fourier
- num_spherical
- sinθ
- nθ

# References:
- Ehrendorfer, M. (2011) Spectral Numerical Weather Prediction Models, Appendix B, Society for Industrial and Applied Mathematics
- Winch, D. (2007) Spherical harmonics, in Encyclopedia of Geomagnetism and Paleomagnetism, Eds Gubbins D. and Herrero-Bervera, E., Springer

# Details (using notation and Eq. references from Ehrendorfer, 2011):

    l=0,1...∞    and m = -l, -l+1, ... l-1, l

    P_{0,0} = 1, such that 1/4π ∫∫YYdS = δ (where Y = spherical harmonics, S = domain surface area)
    P_{m,m} = sqrt((2m+1)/2m) cosθ P_{m-1m-1}
    P_{m+1,m} = sqrt(2m+3) sinθ P_{m m}
    sqrt((l^2-m^2)/(4l^2-1))P_{l,m} = P_{l-1, m} -  sqrt(((l-1)^2-m^2)/(4(l-1)^2 - 1))P_{l-2,m}

    THe normalization assures that 1/2 ∫_{-1}^1 P_{l,m}(sinθ) P_{n,m}(sinθ) dsinθ = δ_{n,l}

    Julia index starts with 1, so qnm[m+1,l+1] = P_l^m
"""
function compute_legendre!(num_fourier, num_spherical, sinθ, nθ)
    qnm = zeros(Float64, num_fourier + 1, num_spherical + 2, nθ)

    cosθ = sqrt.(1 .- sinθ .^ 2)
    ε = zeros(Float64, num_fourier + 1, num_spherical + 2)

    qnm[1, 1, :] .= 1
    for m in 1:num_fourier
        qnm[m + 1, m + 1, :] = -sqrt((2m + 1) / (2m)) .* cosθ .* qnm[m, m, :] # Eq. B.20
        qnm[m, m + 1, :] = sqrt(2m + 1) * sinθ .* qnm[m, m, :] # Eq. B.22
    end
    qnm[num_fourier + 1, num_fourier + 2, :] =
        sqrt(2 * (num_fourier + 2)) * sinθ .*
        qnm[num_fourier + 1, num_fourier + 1, :]

    for m in 0:num_fourier
        for l in (m + 2):(num_spherical + 1)
            ε1 = sqrt(((l - 1)^2 - m^2) ./ (4 * (l - 1)^2 - 1))
            ε2 = sqrt((l^2 - m^2) ./ (4 * l^2 - 1))
            qnm[m + 1, l + 1, :] =
                (sinθ .* qnm[m + 1, l, :] - ε1 * qnm[m + 1, l - 1, :]) / ε2 # Eq. B.18
        end
    end

    return qnm[:, 1:(num_spherical + 1), :]
end

"""
    compute_gaussian!(n)

Compute sin(latitude) and the weight factors for Gaussian integration

# Arguments
- n: number of latitudes

# References
- Ehrendorfer, M., Spectral Numerical Weather Prediction Models, Appendix B, Society for Industrial and Applied Mathematics, 2011

# Details (following notation from Ehrendorfer, 2011):
    Pn(x) is an odd function
    solve half of the n roots and weightes of Pn(x) # n = 2n_half
    P_{-1}(x) = 0
    P_0(x) = 1
    P_1(x) = x
    nP_n(x) = (2n-1)xP_{n-1}(x) - (n-1)P_{n-2}(x)
    P'_n(x) = n/(x^2-1)(xP_{n}(x) - P_{n-1}(x))
    x -= P_n(x)/P'_{n}()
    Initial guess xi^{0} = cos(π(i-0.25)/(n+0.5))
    wi = 2/(1-xi^2)/P_n'(xi)^2
"""
function compute_gaussian!(n)
    itermax = 10000
    tol = 1.0e-15

    sinθ = zeros(Float64, n)
    wts = zeros(Float64, n)

    n_half = Int64(n / 2)
    for i in 1:n_half
        dp = 0.0
        z = cos(pi * (i - 0.25) / (n + 0.5))
        for iter in 1:itermax
            p2 = 0.0
            p1 = 1.0

            for j in 1:n
                p3 = p2 # Pj-2
                p2 = p1 # Pj-1
                p1 = ((2.0 * j - 1.0) * z * p2 - (j - 1.0) * p3) / j  #Pj
            end
            # P'_n
            dp = n * (z * p1 - p2) / (z * z - 1.0)
            z1 = z
            z = z1 - p1 / dp
            if (abs(z - z1) <= tol)
                break
            end
            if iter == itermax
                @error("Compute_Gaussian! does not converge!")
            end
        end

        sinθ[i], sinθ[n - i + 1], = -z, z
        wts[i] = wts[n - i + 1] = 2.0 / ((1.0 - z * z) * dp * dp)
    end

    return sinθ, wts
end

"""
    trans_grid_to_spherical!(mesh::SpectralSphericalMesh, pfield::Array{Float64,2})

Transforms a variable on a Gaussian grid (pfield[nλ, nθ]) into the spherical harmonics domain (var_spherical2d[num_fourier+1, num_spherical+1])

Here λ = longitude, θ = latitude, η = sinθ, m = zonal wavenumber, n = total wavenumber:
var_spherical2d = F_{m,n}    # Output variable in spectral space (Complex{Float64}[num_fourier+1, num_spherical+1])
qwg = P_{m,n}(η)w(η)         # Weighted Legendre polynomials (Float64[num_fourier+1, num_spherical+1, nθ])
var_fourier2d = g_{m, θ}     # Untruncated Fourier transformation (Complex{Float64} [nλ, nθ])
pfield = F(λ, η)             # Input variable on Gaussian grid Float64[nλ, nθ]

# Arguments
- mesh: struct with mesh information
- pfield: variable on Gaussian grid to be transformed

# References
- Ehrendorfer, M., Spectral Numerical Weather Prediction Models, Appendix B, Society for Industrial and Applied Mathematics, 2011
- [Wiin1967](@cite)
"""
function trans_grid_to_spherical!(
    mesh::SpectralSphericalMesh,
    pfield::Array{Float64, 2},
)

    num_fourier, num_spherical = mesh.num_fourier, mesh.num_spherical
    var_fourier2d, var_spherical2d =
        mesh.var_fourier[:, :, 1] * 0.0, mesh.var_spherical[:, :, 1, :] * 0.0
    nλ, nθ, nd = mesh.nλ, mesh.nθ, mesh.nd

    # Retrieve weighted Legendre polynomials
    qwg = mesh.qwg # qwg[m,n,nθ]

    # Fourier transformation
    for j in 1:nθ
        var_fourier2d[:, j] = fft(pfield[:, j], 1) / nλ
    end

    # Complete spherical harmonic transformation
    @assert(nθ % 2 == 0)
    nθ_half = div(nθ, 2)
    for m in 1:(num_fourier + 1)
        for n in m:num_spherical
            var_fourier2d_t = transpose(var_fourier2d[m, :])  # truncates var_fourier(nlon, nhlat) to (nfourier,nlat)
            if (n - m) % 2 == 0
                var_spherical2d[m, n, :] .=
                    (
                        var_fourier2d_t[1:nθ_half] .+
                        var_fourier2d_t[nθ:-1:(nθ_half + 1)]
                    ) .* qwg[m, n, 1:nθ_half] ./ 2.0
            else
                var_spherical2d[m, n, :] .=
                    (
                        var_fourier2d_t[1:nθ_half] .-
                        var_fourier2d_t[nθ:-1:(nθ_half + 1)]
                    ) .* qwg[m, n, 1:nθ_half] ./ 2.0
            end
        end
    end

    return var_spherical2d
end

function compute_wave_numbers(num_fourier::Int64, num_spherical::Int64)
    """
    See wave_numers[i,j] saves the wave number of this basis
    """
    wave_numbers = zeros(Int64, num_fourier + 1, num_spherical + 1)

    for m in 0:num_fourier
        for n in m:num_spherical
            wave_numbers[m + 1, n + 1] = n
        end
    end

    return wave_numbers

end

ext: jl


text:
ENV["GKSwstype"] = "nul"  # For StatsPlots
println("pid: ", getpid())

# Read command line args.
if length(ARGS) > 1
  resultsdir = ARGS[1]
  awsbucket = ARGS[2]
  snum = parse(Int, ARGS[3])
  println("ARGS: ", ARGS)
else
  resultsdir = "results/test/"
  awsbucket = nothing
  snum = 1
end
flush(stdout)

# Create results/images dir if needed
imgdir = "$(resultsdir)/img"
mkpath(imgdir)

import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))

using CytofDensityEstimation
const CDE = CytofDensityEstimation
using LaTeXStrings
using Distributions
using StatsPlots
using HypothesisTests
import Random
using RCall
using BSON
import CytofDensityEstimation.Model.default_ygrid

include("scenarios.jl")

plotsize = (400, 400)

# Simulate data.
Random.seed!(5);
simdata = scenarios(snum)
yC, yT = simdata[:yC], simdata[:yT]

# Plot true data density.
function plot_trude_data_density()
  plot(default_ygrid(), pdf.(simdata[:mmC], default_ygrid()), lw=3,
       label=L"y_C", color=:blue)
  plot!(default_ygrid(), pdf.(simdata[:mmT], default_ygrid()), lw=3,
        label=L"y_T", color=:red)
end
plot_trude_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-true-density.pdf"))
closeall()

# Plot observed data density.
function plot_observed_data_density()
  density(yC[isfinite.(yC)],  lw=3, label=L"y_C", color=:blue)
  density!(yT[isfinite.(yT)], lw=3, label=L"y_T", color=:red)
end
plot_observed_data_density()
plot!(size=plotsize)
savefig(joinpath(imgdir, "data-kde.pdf"))
closeall()

# Define data, prior, and initial state.
Random.seed!(7)
K = 3
data = CDE.Model.Data(yC, yT)
prior_mu = let
  yfinite = [data.yC_finite; data.yT_finite]
  Normal(mean(yfinite), std(yfinite))
end
prior = CDE.Model.Prior(K, mu=prior_mu, nu=LogNormal(1.6, 0.4), p=Beta(100, 100),
                        omega=InverseGamma(.1, .1), psi=Normal(-1, 1))
state = CDE.Model.State(data, prior)
tuners = CDE.Model.Tuners(K, 0.1)

println("Priors:")
for fn in fieldnames(CDE.Model.Prior)
  println(fn, " => ", getfield(prior, fn))
end

# Parameters to fix
# - psi is near truth when mu, omega, nu, and eta are fixed, and n → ∞.
# - nu can be recovered when mu, omega, psi, and eta are fixed, and n → ∞.
# - mu is easy to recover when nu, omega, psi, and eta are fixed.
# - omega is near truth when mu, nu, psi, and eta are fixed; it requires a
#   large amount of data to estimate.
# - eta can be easily recovered when n → ∞.
fix = Symbol[]
# fix = Symbol[:omega]
# fix = Symbol[:psi]
# fix = Symbol[:nu, :psi]
# fix = Symbol[:mu, :omega, :nu, :eta]
# state.mu[1:3] .= simdata[:loc]
# state.nu[1:3] .= simdata[:df]
# state.psi[1:3] .= CDE.Util.toaltskew.(simdata[:scale], simdata[:skew])
# state.omega[1:3] .= CDE.Util.toaltscale.(simdata[:scale], simdata[:skew]) .^ 2
# state.etaC .= 0; state.etaC[:1:3] .= simdata[:etaC]
# state.etaT .= 0; state.etaT[:1:3] .= simdata[:etaT]
monitors = CDE.Model.default_monitors()

# Warmup.
# _tuners = deepcopy(tuners)
# state.beta = 1
# @time _, state, _ = CDE.Model.fit(
#     state, data, prior, _tuners, nsamps=[1], nburn=200, thin=1,
#     fix=[fix; [:beta]], monitors=monitors,
#     rep_beta_flipped=50)

# Run chain.
@time chain, laststate, summarystats = CDE.Model.fit(
    state, data, prior, tuners, nsamps=[10000], nburn=2000, thin=1,
    rep_aux=10,
    fix=fix, monitors=monitors, rep_beta_flipped=50)

# Save results
BSON.bson("$(resultsdir)/results.bson",
          Dict(:chain=>chain, :laststate=>laststate,
               :summarystats=>summarystats, :simdata=>simdata))

# Load via:
# using BSON, CytofDensityEstimation
out = BSON.load("$(resultsdir)/results.bson")

# Print KS statistic.
# ks_fit = ApproximateTwoSampleKSTest(yC, yT)
# println(ks_fit)
@rimport stats as rstats
ks_fit = rstats.ks_test(out[:simdata][:yC], out[:simdata][:yT])
println(ks_fit["p.value"][1])
flush(stdout)

# Print summary statistics.
CDE.Model.printsummary(out[:chain], out[:summarystats])
flush(stdout)

# Plot results.
@time CDE.Model.plotpostsummary(out[:chain], out[:summarystats], yC, yT,
                                imgdir, simdata=out[:simdata], bw_postpred=.2,
                                xlims_=(-6, 6))
flush(stdout)

if awsbucket != nothing
  CDE.Util.s3sync(from=resultsdir, to=awsbucket, tags=`--exclude '*.nfs'`)
end

println("Done!")
flush(stdout)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
using ArgParse
import ArgParse: parse_item

using Formatting
using Statistics

using Plots
# pyplot(reuse=true, size=(900,300))
pyplot()
theme(:ggplot2)

import TimeSeries
import Dates
import FracFin

function ArgParse.parse_item(::Type{AbstractVector{Int}}, S::AbstractString)
    # return ParseAbstractVector{Int}(S) #
    return ParseAbstractVector(Int, S)
end

function ArgParse.parse_item(::Type{NTuple{N, Int}}, S::AbstractString) where N
    return ParseNTuple(Int, S)
end

# function ParseAbstractVector{T}(S::AbstractString) where {T<:Real}
#     A = if something(findfirst(isequal(':'), S), 0) != 0  # step range
#         N = [parse(T, s) for s in split(S, ':')]
#         # collect(N[1]:N[2]:N[3])
#         N[1]:N[2]:N[3]
#     elseif S[1] == '[' && S[end] == ']'  # vector
#         [parse(T, s) for s in split(S[2:end-1], ',')]
#     else
#         error("Invalid input: $S")
#     end
#     # all((A .% 2) .== 0) || error("Invalid input: $S")

#     return A
# end

function ParseAbstractVector(T::Type{P}, S::AbstractString) where P<:Real
    A = if something(findfirst(isequal(':'), S), 0) != 0  # step range
        N = [parse(Int, s) for s in split(S, ':')]
        # collect(N[1]:N[2]:N[3])
        s = (length(N) == 2) ? 1 : N[2]
        StepRange(N[1], s, N[end])
    elseif S[1] == '[' && S[end] == ']'  # vector
        [parse(Int, s) for s in split(S[2:end-1], ',')]
    else
        error("Invalid input: $S")
    end
    # all((A .% 2) .== 0) || error("Invalid input: $S")

    return A
end

function ParseNTuple(T::Type{P}, S::AbstractString) where {P<:Real}
    return (parse(T, s) for s in split(S, ','))
#     A = if S[1] == '(' && S[end] == ')'  # tuple
#         (parse(T, s) for s in split(S[2:end-1], ','))
#     else
#         error("Invalid input: $S")
#     end
#     return A
end

function prepare_fGn(X0::AbstractVector{T}, dlag::Int) where {T<:Real}
    dX = fill(NaN, size(X0))  # fGn estimator is not friendly with NaN
    dX[1+dlag:end] = X0[1+dlag:end]-X0[1:end-dlag]
    return dX
end

function prepare_bspline(X0::AbstractVector{T}, sclrng::AbstractVector{Int}, vm::Int) where {T<:Real}
    W0, Mask = FracFin.cwt_bspline(X0, sclrng, vm, :left)
    t0 = findall(all(Mask, dims=2)[:])[1]
    # t1 = findall(all(Mask, dims=2)[:])[end]
    W0[1:t0-1,:] .= NaN
    return transpose(W0)
end

function equalize_daynight(sdata)
    edata = [sdata[1]]
    for n=2:length(sdata)
        d0 = TimeSeries.values(edata[n-1])
        d1 = TimeSeries.values(sdata[n])
        t1 = TimeSeries.timestamp(sdata[n])
        push!(edata, TimeSeries.TimeArray(t1, d1.-(d1[1]-d0[end]), TimeSeries.colnames(sdata[1])))
    end
    return edata
end

function parse_commandline()
    settings = ArgParseSettings("Apply rolling window estimator of fBm.")

    @add_arg_table settings begin
        "--wsize"
        help = "size of rolling window"
        arg_type = Int
        default = 180
        "--ssize"
        help = "size of sub window"
        arg_type = Int
        default = 60
        "--dlen"
        help = "length of decorrelation"
        arg_type = Int
        default = 5
        "--pov"
        help = "period of innovation"
        arg_type = Int
        default = 5
        "--ipts"
        help = "number of points per day (for intraday data only), 0 for no down-sampling"
        arg_type = Int
        default = 0
        "--idta"
        help = " starting time of truncation (for intraday data only)"
        arg_type = String
        default = "09:05"
        "--idtb"
        help = " ending time of truncation (for intraday data only)"
        arg_type = String
        default = "17:25"
        "--eql"
        help = "Equalize the opening and close values of two consecutive days (for ipts=0 only)"
        action = :store_true
        "--tfmt"
        help = "time format for parsing csv file"
        arg_type = String
        default = "yyyy-mm-dd HH:MM:SS"
        "--ncol"
        help = "index of the column to be analyzed, if the input file contains multiple columns"
        arg_type = Int
        default = 1
        "--verbose", "-v"
        help = "print message"
        action = :store_true
        "infile"
        help = "input file (with suffix '.csv' or '.txt') containing a time series"
        required = true
        "outdir"
        help = "output directory"
        required = true
        "powlaw"
        action = :command        # adds a command which will be read from an argument
        help = "power law method: 'ssize' and 'dlen' are ignored."
        "fGn-MLE"
        action = :command        # adds a command which will be read from an argument
        help = "fGn-MLE"
        "fWn-bspline-MLE"
        action = :command
        help = "fWn-bspline-MLE"
        "bspline-scalogram"
        action = :command
        help = "bspline-scalogram method: 'ssize' and 'dlen' are ignored."
    end

    @add_arg_table settings["powlaw"] begin    # add command arg_table: same as usual, but invoked on s["cmd"]
        "--pow"
        help = "compute the 'pow'-th moment "
        arg_type = Float64
        default = 1.
        "--dlags"
        help = "time lags for finite difference"
        arg_type = AbstractVector{Int}
        default = 2:10
    end

    settings["powlaw"].description = "power law description"  # this is how settings are tweaked for commands
    settings["powlaw"].commands_are_required = true # this makes the sub-commands optional
    settings["powlaw"].autofix_names = true # this uses dashes in long options, underscores in auto-generated dest_names

    @add_arg_table settings["fGn-MLE"] begin    # add command arg_table: same as usual, but invoked on s["cmd"]
        "--dlag"
        help = "time lag for finite difference"
        arg_type = Int
        default = 1
    end

    settings["fGn-MLE"].description = "fGn description"  # this is how settings are tweaked for commands
    settings["fGn-MLE"].commands_are_required = true # this makes the sub-commands optional
    settings["fGn-MLE"].autofix_names = true # this uses dashes in long options, underscores in auto-generated dest_names

    @add_arg_table settings["fWn-bspline-MLE"] begin    # add command arg_table: same as usual, but invoked on s["cmd"]
        "--vm"
        help = "vanishing moments of the wavelet"
        arg_type = Int
        default = 2
        "--sclidx"
        help = "indices of integer scales (scale must be an even number)"
        arg_type = AbstractVector{Int}
        default = 2:10
    end

    settings["fWn-bspline-MLE"].description = "fWn description"  # this is how settings are tweaked for commands
    settings["fWn-bspline-MLE"].commands_are_required = true # this makes the sub-commands optional

    @add_arg_table settings["bspline-scalogram"] begin    # add command arg_table: same as usual, but invoked on s["cmd"]
        "--vm"
        help = "vanishing moments of the wavelet"
        arg_type = Int
        default = 2
        "--sclidx"
        help = "indices of even integer scales (scale must be an even number)"
        arg_type = AbstractVector{Int}
        default = 2:10
        "--pow"
        help = "compute the 'pow'-th moment "
        arg_type = Float64
        default = 1.
    end

    settings["bspline-scalogram"].description = "fWn description"  # this is how settings are tweaked for commands
    settings["bspline-scalogram"].commands_are_required = true # this makes the sub-commands optional

    return parse_args(settings)
end

function main()
    parsed_args = parse_commandline()
    # println("Parsed args:")
    # for (arg,val) in parsed_args
    #     println("  $arg  =>  $val")
    # end

    # recover global parameters and options
    infile = parsed_args["infile"]
    outdir0 = parsed_args["outdir"]

    i1 = something(findlast(isequal('/'), infile), 0)
    i2 = something(findlast(isequal('.'), infile), length(infile)+1)
    sname = infile[i1+1:i2-1]  # filename without extension

    cmd = parsed_args["%COMMAND%"]  # name of command
    wsize = parsed_args["wsize"]  # size of rolling window
    ssize = parsed_args["ssize"]  # size of sub window
    dlen = parsed_args["dlen"]  # length of decorrelation
    pov = parsed_args["pov"]  # period of innovation
    ipts = parsed_args["ipts"]  # number of pts per day
    eql = parsed_args["eql"]
    idta = parsed_args["idta"]
    idtb = parsed_args["idtb"]
    tfmt = parsed_args["tfmt"]  # time format
    ncol = parsed_args["ncol"]  # index of column
    verbose = parsed_args["verbose"]  # print messages

    # recover command options
    prepare::Function = X->()  # function for data preparation
    estim::Function = X->()
    trans::Function = X->vec(X)
    outstr::String = ""

    if cmd  == "powlaw"
        ssize = wsize
        dlen = 1
        pow = parsed_args[cmd]["pow"]  # time-lags of finite difference
        dlags = parsed_args[cmd]["dlags"]  # time-lags of finite difference
        outstr = format("wsize[{}]_ipts[{}]_pov[{}]_pow[{}]_dlags[{}]", wsize, ipts, pov, pow, dlags)
        prepare = X -> X
        estim = X -> FracFin.powlaw_estim(X, dlags, pow)
    elseif cmd  == "fGn-MLE"
        dlag = parsed_args[cmd]["dlag"]  # time-lag of finite difference
        outstr = format("wsize[{}]_ssize[{}]_dlen[{}]_ipts[{}]_pov[{}]_dlag[{}]", wsize, ssize, dlen, ipts, pov, dlag)
        prepare = X -> prepare_fGn(X, dlag)
        estim = X -> FracFin.fGn_MLE_estim(X, dlag; method=:optim, ε=1e-2)
    elseif cmd == "fWn-bspline-MLE"
        sclrng = 2 * parsed_args[cmd]["sclidx"]  # range of scales
        vm = parsed_args[cmd]["vm"]  # vanishing moments
        outstr = format("wsize[{}]_ssize[{}]_dlen[{}]_ipts[{}]_pov[{}]_vm[{}]_sclrng[{}]", wsize, ssize, dlen, ipts, pov, vm, sclrng)
        prepare = X -> prepare_bspline(X, sclrng, vm)
        estim = X -> FracFin.fWn_bspline_MLE_estim(X, sclrng, vm; method=:optim, ε=1e-2)
    elseif cmd == "bspline-scalogram"
        ssize = wsize
        dlen = 1
        sclrng = 2 * parsed_args[cmd]["sclidx"]  # range of scales
        pow = parsed_args[cmd]["pow"]  # time-lags of finite difference
        vm = parsed_args[cmd]["vm"]  # vanishing moments
        outstr = format("wsize[{}]_ipts[{}]_pov[{}]_pow[{}]_vm[{}]_sclrng[{}]", wsize, ipts, pov, pow, vm, sclrng)
        prepare = X -> prepare_bspline(X, sclrng, vm)
        estim = X -> FracFin.bspline_scalogram_estim(X, sclrng, vm; p=pow, mode=:left)
        # trans = X -> var(X, dims=2)
        # trans = X -> vec(FracFin.robustvar(X, dims=2))
        # trans = X -> vec(mean(abs.(X .- mean(X, dims=2)).^pow, dims=2))
        trans = X -> vec(mean(abs.(X).^pow, dims=2))
    else
        error("Unknown command")
    end

    # load data
    if verbose
        printfmtln("Loading file {}...", infile)
    end
    # data0 = CSV.read(infile)
    toto = TimeSeries.readtimearray(infile, format=tfmt, delim=',')
    cname = TimeSeries.colnames(toto)[ncol]  # name of the column
    unit = Dates.Minute(minimum(diff(TimeSeries.timestamp(toto))))
    data0 = toto[cname]

    if Int(unit/Dates.Minute(1)) < 1440  # intraday data
        tdta, tdtb = Dates.Time(idta), Dates.Time(idtb)
        dha, dhb = Dates.Hour(tdta), Dates.Hour(tdtb)
        dma, dmb = Dates.Minute(tdta), Dates.Minute(tdtb)
        daytime = (dha+dma, dhb+dmb)

        Xm, Tm, Xt, Tt = [], [], [], []
        sdata0 = []

        if ipts>0
            wh = (dhb-dha) ÷ ipts + dmb-dma # length of intraday window
            nh = Int((dhb-dha)/Dates.Hour(1)) # number of hours in intraday data
            dh = nh ÷ ipts
            windows = [(w=daytime[1]+Dates.Hour(n); (w,w+wh)) for n=0:dh:(nh-1)]

            # splitted data of identical length
            Xms = []
            Tms = []
            for (wa,wb) in windows
                sdata0 = FracFin.window_split_timearray(data0, Dates.Hour(24), (wa,wb), fillmode=:fb, endpoint=false)
                push!(Xms, log.(hcat([TimeSeries.values(v) for v in sdata0]...)))
                # println(length(sdata0))
                # println(size(Xms[end]))
                push!(Tms, hcat([TimeSeries.timestamp(v) for v in sdata0]...))
            end

            Xm = reshape(vcat(Xms...), (size(Xms[1],1), :))
            Tm = reshape(vcat(Tms...), (size(Tms[1],1), :))
            Xt = mean(Xm, dims=1)[:]
            Tt = [Dates.DateTime(t) for t in Tm[1,:]]
        else
            sdata0 = FracFin.window_split_timearray(data0, Dates.Hour(24), daytime, fillmode=:fb, endpoint=false)
            sdata1 = eql ? sdata0 : equalize_daynight(sdata0)
            data1 = vcat(sdata1...)
            Tt = TimeSeries.timestamp(data1)
            Xm = log.(reshape(TimeSeries.values(data1), 1, :))
            Xt = log.(TimeSeries.values(vcat(sdata0...)))
            # Xm = log.(reshape(vcat(TimeSeries.values.(sdata0)...), 1, :))
            # Tt = vcat(TimeSeries.timestamp.(sdata0)...)
            # Xt = Xm[1, :]
        end
        # println(size(Xm))
        # println(size(Tm))
        # println(size(Xt))
        # println(size(Tt))

        Res = []
        for r = 1:size(Xm,1)
            Xi = prepare(Xm[r,:])  # input to rolling estimator
            res = FracFin.rolling_estim(estim, Xi, (wsize,ssize,dlen), pov, trans; mode=:causal)
            push!(Res, res)
        end

        Hm = hcat([[x[2][1][1] for x in res] for res in Res]...)'
        σm = hcat([[x[2][1][2] for x in res] for res in Res]...)'
        tidx = [x[1] for x in Res[1]]

        # embedding into a longer time series
        Ht = fill(NaN, length(Tt)); Ht[tidx] = mean(Hm, dims=1)[:]
        σt = fill(NaN, length(Tt)); σt[tidx] = mean(σm, dims=1)[:]
        Hs = fill(NaN, length(Tt)); Hs[tidx] = std(Hm, dims=1)[:]
        σs = fill(NaN, length(Tt)); σs[tidx] = std(σm, dims=1)[:]

        At = TimeSeries.TimeArray(Tt, [Xt Ht Hs σt σs], Symbol.(["Log_Price", "Hurst", "Hurst_std", "σ", "σ_std"]))

        # make the output folder
        outdir = format("{}/{}/{}/{}/{}/", outdir0, sname, cname, cmd, outstr)
        try
            mkpath(outdir)
        catch
        end

        TimeSeries.writetimearray(At, format("{}/estimate.csv",outdir))
        # At = TimeArray(T0[tidx], [X0[tidx] Ht σt Hs σs], ["Log_Price", "Hurst", "Hurst_std" "σ" "σ_std"])
        # TimeSeries.writetimearray(At, format("{}/estimate-ts.csv",outdir))

        fig1 = plot(Tt, Xt, ylabel="Log Price", label="")
        fig2 = plot(Tt, Ht, ribbon=(2*Hs, 2*Hs), ylim=[0,1], shape=:circle, ylabel="Hurst", label="")
        hline!(fig2, [0.5], width=[3], color=[:red], label="")
        fig3 = plot(Tt, σt, ribbon=(2*σs, 2*σs), shape=:circle, ylabel="σ", label="")
        plot(fig1, fig2, fig3, layout=(3,1), size=(1200,900), legend=true)
        savefig(format("{}/estimate.pdf",outdir))

        if ipts==0
            sdata1 = FracFin.window_split_timearray(At, Dates.Hour(24); fillmode=:o, endpoint=true)
            for Y0 in sdata1
                T0 = TimeSeries.timestamp(Y0)
                fig1 = plot(Y0[:Log_Price], ylabel="Log Price", label="")
                fig2 = plot(Y0[:Hurst], shape=:circle, ylabel="Hurst", label="")
                fig3 = plot(Y0[:σ], shape=:circle, ylabel="σ", label="")
                fig4 = plot(Y0[:Hurst], ylim=[0,1], shape=:circle, ylabel="Hurst", label="")
                hline!(fig4, [0.5], width=[3], color=[:red], label="")
                plot(fig1, fig2, fig3, fig4, layout=(4,1), size=(1200,900), xticks=T0[1]:Dates.Hour(1):T0[end], legend=true)
                savefig(format("{}/{}.pdf",outdir, Dates.Date(T0[1])))
                # TimeSeries.writetimearray(Yt, format("{}/csv/{}.csv",outdir, day))
            end
        end

        if verbose
            printfmtln("Outputs saved in {}", outdir)
        end
    else

    end
end

main()

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
using ITensors
using OptimKit
using Random
using Zygote
using Test

include("utils/circuit.jl")

@testset "optimization" begin
  @testset "Energy minimization" begin
    N = 3
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓")
    H = prod(Hmpo)
    ψ₀ = prod(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::ITensor, ψ::ITensor)
      ψdag = dag(ψ)
      return (ψdag' * H * ψ)[] / (ψdag * ψ)[]
    end
    E(ψ::ITensor) = E(H, ψ)
    ∇E(ψ::ITensor) = E'(ψ)
    fg(ψ::ITensor) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(;
      c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0, verbosity=0
    )
    algorithm = LBFGS(3; maxiter=20, gradtol=1e-8, linesearch=linesearch)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    D, _ = eigen(H; ishermitian=true)
    @test E(H, ψ) < E(H, ψ₀)
    @test E(H, ψ) ≈ minimum(D)
  end

  @testset "Energy minimization (MPS)" begin
    N = 4
    χ = 4
    s = siteinds("S=1/2", N; conserve_qns=true)
    os = OpSum()
    for n in 1:(N - 1)
      os .+= 0.5, "S+", n, "S-", n + 1
      os .+= 0.5, "S-", n, "S+", n + 1
      os .+= "Sz", n, "Sz", n + 1
    end
    Hmpo = MPO(os, s)

    Random.seed!(1234)
    ψ₀mps = randomMPS(s, n -> isodd(n) ? "↑" : "↓"; linkdims=χ)

    H = ITensors.data(Hmpo)
    ψ₀ = ITensors.data(ψ₀mps)
    # The Rayleigh quotient to minimize
    function E(H::Vector{ITensor}, ψ::Vector{ITensor})
      N = length(ψ)
      ψdag = dag.(addtags.(ψ, "bra"; tags="Link"))
      ψ′dag = prime.(ψdag)
      e = ITensor(1.0)
      for n in 1:N
        e = e * ψ′dag[n] * H[n] * ψ[n]
      end
      norm = ITensor(1.0)
      for n in 1:N
        norm = norm * ψdag[n] * ψ[n]
      end
      return e[] / norm[]
    end
    E(ψ) = E(H, ψ)
    ∇E(ψ) = E'(ψ)
    fg(ψ) = (E(ψ), ∇E(ψ))
    linesearch = HagerZhangLineSearch(; c₁=0.1, c₂=0.9, ϵ=1e-6, θ=1 / 2, γ=2 / 3, ρ=5.0)
    algorithm = LBFGS(5; maxiter=50, gradtol=1e-4, linesearch=linesearch, verbosity=0)
    ψ, fψ, gψ, numfg, normgradhistory = optimize(fg, ψ₀, algorithm)
    sweeps = Sweeps(5)
    setmaxdim!(sweeps, χ)
    fψmps, ψmps = dmrg(Hmpo, ψ₀mps, sweeps; outputlevel=0)
    @test E(H, ψ) ≈ inner(ψmps', Hmpo, ψmps) / inner(ψmps, ψmps) rtol = 1e-2
  end

  @testset "State preparation (full state)" begin
    function Rylayer(N, θ⃗)
      return [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:N]
    end

    function CXlayer(N)
      return [("CX", (n, n + 1)) for n in 1:2:(N - 1)]
    end

    # The variational circuit we want to optimize
    function variational_circuit(θ⃗)
      N = length(θ⃗)
      return vcat(Rylayer(N, θ⃗), CXlayer(N))
    end

    N = 4
    Random.seed!(1234)
    θ⃗ = 2π .* rand(N)
    gates = variational_circuit(θ⃗)

    s = siteinds("Qubit", N)
    ψₘₚₛ = MPS(s, "0")
    ψ = prod(ψₘₚₛ)
    U = buildcircuit(gates, s)
    # Create the target state
    Uψ = apply(U, ψ)

    @test inner_circuit(Uψ, U, ψ) ≈ 1

    function loss(θ⃗)
      gates = variational_circuit(θ⃗)
      U = buildcircuit(gates, s)
      return -abs(inner_circuit(Uψ, U, ψ))^2
    end

    θ⃗₀ = randn!(copy(θ⃗))
    fg(x) = (loss(x), convert(Vector, loss'(x)))
    θ⃗ₒₚₜ, fₒₚₜ, gₒₚₜ, numfg, normgradhistory = optimize(fg, θ⃗₀, GradientDescent())
    @test loss(θ⃗ₒₚₜ) ≈ loss(θ⃗) rtol = 1e-2
  end

  @testset "State preparation (MPS)" begin
    for gate in ["Ry"]#="Rx", =#
      nsites = 4 # Number of sites
      nlayers = 2 # Layers of gates in the ansatz
      gradtol = 1e-3 # Tolerance for stopping gradient descent

      # A layer of the circuit we want to optimize
      function layer(nsites, θ⃗)
        gate_layer = [(gate, (n,), (θ=θ⃗[n],)) for n in 1:nsites]
        CX_layer = [("CX", (n, n + 1)) for n in 1:2:(nsites - 1)]
        return [gate_layer; CX_layer]
      end

      # The variational circuit we want to optimize
      function variational_circuit(nsites, nlayers, θ⃗)
        range = 1:nsites
        circuit = layer(nsites, θ⃗[range])
        for n in 1:(nlayers - 1)
          circuit = [circuit; layer(nsites, θ⃗[range .+ n * nsites])]
        end
        return circuit
      end

      Random.seed!(1234)

      θ⃗ᵗᵃʳᵍᵉᵗ = 2π * rand(nsites * nlayers)
      𝒰ᵗᵃʳᵍᵉᵗ = variational_circuit(nsites, nlayers, θ⃗ᵗᵃʳᵍᵉᵗ)

      s = siteinds("Qubit", nsites)
      Uᵗᵃʳᵍᵉᵗ = ops(𝒰ᵗᵃʳᵍᵉᵗ, s)

      ψ0 = MPS(s, "0")

      # Create the random target state
      ψᵗᵃʳᵍᵉᵗ = apply(Uᵗᵃʳᵍᵉᵗ, ψ0; cutoff=1e-8)

      #
      # The loss function, a function of the gate parameters
      # and implicitly depending on the target state:
      #
      # loss(θ⃗) = -|⟨θ⃗ᵗᵃʳᵍᵉᵗ|U(θ⃗)|0⟩|² = -|⟨θ⃗ᵗᵃʳᵍᵉᵗ|θ⃗⟩|²
      #
      function loss(θ⃗)
        nsites = length(ψ0)
        s = siteinds(ψ0)
        𝒰θ⃗ = variational_circuit(nsites, nlayers, θ⃗)
        Uθ⃗ = ops(𝒰θ⃗, s)
        ψθ⃗ = apply(Uθ⃗, ψ0; cutoff=1e-8)
        return -abs(inner(ψᵗᵃʳᵍᵉᵗ, ψθ⃗))^2
      end

      θ⃗₀ = randn!(copy(θ⃗ᵗᵃʳᵍᵉᵗ))

      @test loss(θ⃗₀) ≉ loss(θ⃗ᵗᵃʳᵍᵉᵗ)

      loss_∇loss(x) = (loss(x), convert(Vector, loss'(x)))
      @show gate
      algorithm = LBFGS(; gradtol=gradtol, verbosity=2)
      θ⃗ₒₚₜ, lossₒₚₜ, ∇lossₒₚₜ, numfg, normgradhistory = optimize(
        loss_∇loss, θ⃗₀, algorithm
      )

      @test loss(θ⃗ₒₚₜ) ≈ loss(θ⃗ᵗᵃʳᵍᵉᵗ) rtol = 1e-5
    end
  end

  @testset "VQE (MPS)" begin
    nsites = 4 # Number of sites
    nlayers = 2 # Layers of gates in the ansatz
    gradtol = 1e-3 # Tolerance for stopping gradient descent

    # The Hamiltonian we are minimizing
    function ising_hamiltonian(nsites; h)
      ℋ = OpSum()
      for j in 1:(nsites - 1)
        ℋ += -1, "Z", j, "Z", j + 1
      end
      for j in 1:nsites
        ℋ += h, "X", j
      end
      return ℋ
    end

    # A layer of the circuit we want to optimize
    function layer(nsites, θ⃗)
      RY_layer = [("Ry", (n,), (θ=θ⃗[n],)) for n in 1:nsites]
      CX_layer = [("CX", (n, n + 1)) for n in 1:2:(nsites - 1)]
      return [RY_layer; CX_layer]
    end

    # The variational circuit we want to optimize
    function variational_circuit(nsites, nlayers, θ⃗)
      range = 1:nsites
      circuit = layer(nsites, θ⃗[range])
      for n in 1:(nlayers - 1)
        circuit = [circuit; layer(nsites, θ⃗[range .+ n * nsites])]
      end
      return circuit
    end

    s = siteinds("Qubit", nsites)

    h = 1.3
    ℋ = ising_hamiltonian(nsites; h=h)
    H = MPO(ℋ, s)
    ψ0 = MPS(s, "0")

    #
    # The loss function, a function of the gate parameters
    # and implicitly depending on the Hamiltonian and state:
    #
    # loss(θ⃗) = ⟨0|U(θ⃗)† H U(θ⃗)|0⟩ = ⟨θ⃗|H|θ⃗⟩
    #
    function loss(θ⃗)
      nsites = length(ψ0)
      s = siteinds(ψ0)
      𝒰θ⃗ = variational_circuit(nsites, nlayers, θ⃗)
      Uθ⃗ = ops(𝒰θ⃗, s)
      ψθ⃗ = apply(Uθ⃗, ψ0; cutoff=1e-8)
      return inner(ψθ⃗', H, ψθ⃗; cutoff=1e-8)
    end

    Random.seed!(1234)
    θ⃗₀ = 2π * rand(nsites * nlayers)

    loss_∇loss(x) = (loss(x), convert(Vector, loss'(x)))
    algorithm = LBFGS(; gradtol=gradtol, verbosity=0)
    θ⃗ₒₚₜ, lossₒₚₜ, ∇lossₒₚₜ, numfg, normgradhistory = optimize(loss_∇loss, θ⃗₀, algorithm)

    sweeps = Sweeps(5)
    setmaxdim!(sweeps, 10)
    e_dmrg, ψ_dmrg = dmrg(H, ψ0, sweeps; outputlevel=0)

    @test loss(θ⃗ₒₚₜ) ≈ e_dmrg rtol = 1e-1
  end
end


text:
nfeatures, ntargets = 2, 3
τ = 100
X = rand(nfeatures, τ) # feature matrix: nfeatures x nobs
Y = rand(ntargets, τ) # target matrix: ntargets x nobs
model = Affine{Float64}(nfeatures, ntargets) # some learnable Transformation?
strategy = SGD(model) # the spec and state of a learning algorithm

# loop through batches checking for early stopping after each batch
function learn!(model, strategy, data::BatchIterator)
    pre_hook(model, strategy)
    for (i, batch) in enumerate(data)
        # update the params for this batch
        learn!(model, strategy, batch)

        iter_hook(model, strategy, i)
        finished(model, strategy, i) && break
    end
    post_hook(model, strategy)
end

function learn!(model, strategy, batch::Batch)
    θ = params(model)
    ∇ = grad(model)
    ∇avg = zeros(θ)
    scalar = 1 / length(batch)
    for (input,target) in batch
        # forward and backward passes for this datapoint
        transform!(model, target, input)
        grad!(model)

        # add to the total param change for this strategy/gradient
        for (i,j) in zip(eachindex(∇avg), eachindex(∇))
            ∇avg[i] += ∇[j] * scalar
        end
    end

    # update the params using the average gradient
    update!(θ, strategy, ∇avg)
end

# NOTE: a strategy would have to implement update! and finished... the rest is optional

# do a kfolds training

scores = []
for (train, validate) in kfolds(X, Y, k = 5, batchsize = 20)
    # here train and validate are type BatchIterator
    learn!(model, strategy, train)

    # add the score (accuracy? loss?)
    push!(scores, score(model, validate))
end
info("Accuracy: $(mean(scores)) (+/- $(std(scores)))")

max_stars_repo_name: JuliaML/StochasticOptimization

max_issues_repo_name: JuliaML/StochasticOptimization

max_forks_repo_name: JuliaML/StochasticOptimization


text:
using ExprTools
using Test

@testset "ExprTools.jl" begin
    include("function.jl")
    include("method.jl")
    include("def_tools.jl")
    include("type_utils.jl")
end


text:
using Clang
using Clang.LibClang
using Test

@testset "typedef" begin
    trans_unit = parse_header(joinpath(@__DIR__, "c", "typedef.h"))
    ctx = DefaultContext()
    push!(ctx.trans_units, trans_unit)

    # get root cursor
    root_cursor = getcursor(trans_unit)
    cursors = children(root_cursor)
    # case 1
    volatile_array = cursors[1]
    volatile_type = underlying_type(volatile_array)
    @test kind(volatile_type) != CXType_Unexposed
    @test spelling(volatile_array) == "foo"
    # case 2
    char_array = cursors[2]
    char_array_type = underlying_type(char_array)
    @test kind(char_array_type) != CXType_Unexposed
    @test spelling(char_array) == "arrayType"
    # case 3
    struct_alias = cursors[4]
    struct_alias_type = underlying_type(struct_alias)
    @test kind(struct_alias_type) != CXType_Unexposed
    @test spelling(struct_alias) == "alias"
    # case 4
    node = cursors[5]
    node_type = type(node)
    @test canonical(node) == cursors[5]
    @test canonical(node_type) |> typedecl == cursors[7]
    wrap!(ctx, cursors[5])
    expr = :(struct Node
                data::Cint
                nextptr::Ptr{Node}
            end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:Node].items[1] == expr

    node_typedef = cursors[6]
    node_typedef_type = underlying_type(node_typedef)
    @test kind(node_typedef_type) != CXType_Unexposed
    @test spelling(node_typedef) == "Node"
    # issue #213
    wrap!(ctx, cursors[9])
    expr = :(const bar = _bar)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:bar].items[1] == expr
    # issue #214
    wrap!(ctx, cursors[10])
    expr = :(const secretunion = Cvoid)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:secretunion].items[1] == expr
    wrap!(ctx, cursors[11])
    expr = :(const SECRET = secretunion)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:SECRET].items[1] == expr
end

@testset "typedef_anon" begin
    # parse file
    trans_unit = parse_header(joinpath(@__DIR__, "c", "typedef_anon.h"))
    ctx = DefaultContext()
    push!(ctx.trans_units, trans_unit)

    # get root cursor
    root_cursor = getcursor(trans_unit)
    cursors = children(root_cursor)
    # case 1
    anonymous_typedef_enum = cursors[1]
    @test name(anonymous_typedef_enum) == ""
    next_cursor = cursors[2]
    @test name(next_cursor) == "Fruit"
    @test kind(next_cursor) == CXCursor_TypedefDecl
    @test name(children(next_cursor)[1]) == ""
    @test is_typedef_anon(anonymous_typedef_enum, next_cursor) == true
    # case 2
    common_enum = cursors[3]
    @test name(common_enum) == "Banana"
    next_cursor = cursors[4]
    @test kind(next_cursor) == CXCursor_TypedefDecl
    @test name(children(next_cursor)[1]) != ""
    @test is_typedef_anon(common_enum, next_cursor) == false
    # case 3
    anonymous_typedef_struct = cursors[5]
    @test name(anonymous_typedef_struct) == ""
    next_cursor = cursors[6]
    @test name(next_cursor) == "Foo"
    @test kind(next_cursor) == CXCursor_TypedefDecl
    @test name(children(next_cursor)[1]) == ""
    @test is_typedef_anon(anonymous_typedef_struct, next_cursor) == true
    # case 4
    attribute_struct = cursors[7]
    @test name(attribute_struct) == ""
    next_cursor = cursors[8]
    @test name(next_cursor) == "char2"
    @test kind(next_cursor) == CXCursor_TypedefDecl
    @test name(children(next_cursor)[1]) == ""
    @test is_typedef_anon(attribute_struct, next_cursor) == true
    # case 5
    anonymous_typedef_union = cursors[9]
    @test name(anonymous_typedef_union) == ""
    next_cursor = cursors[10]
    @test name(next_cursor) == "handle"
    @test kind(next_cursor) == CXCursor_TypedefDecl
    @test name(children(next_cursor)[1]) == ""
    @test is_typedef_anon(anonymous_typedef_union, next_cursor) == true

    ctx.children = cursors
    ctx.children_index = 9
    wrap!(ctx, anonymous_typedef_union)
    expr = :(struct handle
                 ptr::Ptr{Cvoid}
             end)
    Base.remove_linenums!(expr)
    @test ctx.common_buffer[:handle].items[1] == expr
end

@testset "typedef_pointer" begin
    trans_unit = parse_header(joinpath(@__DIR__, "c", "typedef_pointer.h"))
    ctx = DefaultContext()
    push!(ctx.trans_units, trans_unit)

    # get root cursor
    root_cursor = getcursor(trans_unit)
    cursors = children(root_cursor)
    # case 1
    int_pointer = cursors[1]
    int_pointer_type = underlying_type(int_pointer)
    @test kind(int_pointer_type) != CXType_Unexposed
    @test spelling(int_pointer) == "intptr"
    # case 2
    function_pointer = cursors[2]
    function_pointer_type = underlying_type(function_pointer)
    @test kind(function_pointer_type) != CXType_Unexposed
    @test spelling(function_pointer) == "funcptr"
    # case 3
    FOO = cursors[4]
    @test clang2julia(underlying_type(FOO)) == :FOO

    FOOP = cursors[5]
    @test clang2julia(underlying_type(FOOP)) == :(Ptr{FOO})

    wrap!(ctx, cursors[6])
    expr = :(
        function test1(a::Cdouble,b::Cdouble,c::FOOP)
            ccall((:test1,libxxx),FOOP,(Cdouble,Cdouble,FOOP),a,b,c)
        end)
    Base.remove_linenums!(expr)
    @test ctx.api_buffer[1] == expr
end


text:
list = ["NAME: Siddhant Sharma","EMAIL: siddhaantsharma.ss@gmail.com", "SLACKUSER: @Siddhant", "BIOSTACK: Medicinal chemistry and cheminformatics","TWITTER: @sidd2508", "HAMMINGD: 5"]
 for i in list
       println(i)
       end

#to do hamming distances quickly in Julia, use https://github.com/JuliaStats/Distances.jl

max_stars_repo_name: ssiddhantsharma/team-greider

max_issues_repo_name: ssiddhantsharma/team-greider

max_forks_repo_name: ssiddhantsharma/team-greider


text:
export @MT

struct MTBatchSampler{A} <: AbstractIterativeSampler
    samplers::A
end

macro MT_call(fun, val)
    push!(fun.args, Expr(:kw, :par_type, val))
    return quote
        $(esc(fun))
    end
end

macro MT(fun)
    push!(fun.args, Expr(:kw, :par_type, :_tc))
    return quote
        # a global thread_cache
        #TODO make it gensymmed
        global _tc = ThreadsCache()
        s = Vector{Any}(undef, Threads.nthreads())
        Threads.@threads for i=1:Threads.nthreads()
            s[Threads.threadid()] = $(esc(fun))
        end

        # make it type stable
        T = typeof(first(s))
        sT = Vector{T}(undef, length(s))
        sT .= s
        MTBatchSampler(sT)
    end
end

function sample!(s::MTBatchSampler)
    N = length(s.samplers)
    data_1 = Vector{Any}()
    data_prec = Vector{Any}(undef, N)
    Threads.@threads for i=1:N
        L, prec = sample!(s.samplers[i])
        i == 1 && push!(data_1, L)
        data_prec[i] = prec
    end
    T = typeof(first(data_prec))
    data_prec = Vector{T}(data_prec)
    return data_1[1], data_prec
end

function compute_observables(s::MTBatchSampler)
    data = Vector{Any}()
    Threads.@threads for i=1:length(s.samplers)
        res = compute_observables(s.samplers[i])
        i == 1 && push!(data, res)
    end
    return data[1]
end

function add_observable!(s::MTBatchSampler, name, obs)
    Threads.@threads for i=1:length(s.samplers)
        add_observable!(s.samplers[i], name, obs)
    end
    return s
end

function precondition!(data::Vector, params, iter_n)
    parallel_preconditioner = false
    data1 = first(data)
    if params isa SR
        if first(data1.S) isa SrMat
            if first(data1.S).parallel_cache isa ThreadsCache
                parallel_preconditioner = true
            end
        end
    end

    if parallel_preconditioner
        _par_precondition!(data, params, iter_n)
    else
        precondition!(data1, params, iter_n)
    end
end

function _par_precondition!(data, params::SR, iter_n)
    res_data = Vector{Any}()
    Threads.@threads for i=1:length(data)
        res = _precondition!(data[i], params, iter_n)
        i == 1 && push!(res_data, res)
    end
    return res_data[1]
end


text:
export
    decimalperiod,
    ispandigital

#=
http://mathworld.wolfram.com/DecimalExpansion.html

"When a rational number m/n with (m,n)==1 is expanded, the period begins
after s terms and has length t, where s and t are the smallest numbers
satisfying 10^s=10^(s+t) (mod n)."
=#
decimalperiod(_, n) = begin
    n0 = one(n)
    for (p, k) in factor(n)
        p == 2 && continue
        p == 5 && continue
        n0 *= p^k
    end
    n0 == 1 ? zero(n) : multiplicativeorder(10, n0)
end

ispandigital(n::Integer) = digits(n) |> ispandigital
ispandigital(digits::Array{T,1}) where {T<:Number} = begin
    r = 0
    l = length(digits)
    for d in digits
        1 <= d <= l || return false
        r |= 1 << (d - 1)
    end
    r == 2^l - 1
end


text:
using DegeneracyHunter

# Create an optimization model 'm'
include("zero_flowrates.jl")

# Load MILP solver
using Gurobi

# Search for irreducible degenerate sets
sets, dd, ds = DegeneracyHunter.degeneracyHunter(m, GurobiSolver());

println("Done.")


text:
# Decision trees as found in DecisionTree Julia package.
module DecisionTreeWrapper

importall Orchestra.Types
importall Orchestra.Util

import DecisionTree
DT = DecisionTree

export PrunedTree,
       RandomForest,
       DecisionStumpAdaboost,
       fit!,
       transform!

# Pruned ID3 decision tree.
type PrunedTree <: Learner
  model
  options

  function PrunedTree(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_options => {
        # Merge leaves having >= purity_threshold combined purity.
        :purity_threshold => 1.0
      }
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(tree::PrunedTree, instances::Matrix, labels::Vector)
  impl_options = tree.options[:impl_options]
  tree.model = DT.build_tree(labels, instances)
  tree.model = DT.prune_tree(tree.model, impl_options[:purity_threshold])
end
function transform!(tree::PrunedTree, instances::Matrix)
  return DT.apply_tree(tree.model, instances)
end

# Random forest (C4.5).
type RandomForest <: Learner
  model
  options

  function RandomForest(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_options => {
        # Number of features to train on with trees.
        :num_subfeatures => nothing,
        # Number of trees in forest.
        :num_trees => 10,
        # Proportion of trainingset to be used for trees.
        :partial_sampling => 0.7
      }
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(forest::RandomForest, instances::Matrix, labels::Vector)
  # Set training-dependent options
  impl_options = forest.options[:impl_options]
  if impl_options[:num_subfeatures] == nothing
    num_subfeatures = size(instances, 2)
  else
    num_subfeatures = impl_options[:num_subfeatures]
  end
  # Build model
  forest.model = DT.build_forest(
    labels,
    instances,
    num_subfeatures,
    impl_options[:num_trees],
    impl_options[:partial_sampling]
  )
end

function transform!(forest::RandomForest, instances::Matrix)
  return DT.apply_forest(forest.model, instances)
end

# Adaboosted C4.5 decision stumps.
type DecisionStumpAdaboost <: Learner
  model
  options

  function DecisionStumpAdaboost(options=Dict())
    default_options = {
      # Output to train against
      # (:class).
      :output => :class,
      # Options specific to this implementation.
      :impl_options => {
        # Number of boosting iterations.
        :num_iterations => 7
      }
    }
    new(nothing, nested_dict_merge(default_options, options))
  end
end

function fit!(adaboost::DecisionStumpAdaboost,
  instances::Matrix, labels::Vector)

  # NOTE(svs14): Variable 'model' renamed to 'ensemble'.
  #              This differs to DecisionTree
  #              official documentation to avoid confusion in variable
  #              naming within Orchestra.
  ensemble, coefficients = DT.build_adaboost_stumps(
    labels, instances, adaboost.options[:impl_options][:num_iterations]
  )
  adaboost.model = {
    :ensemble => ensemble,
    :coefficients => coefficients
  }
end

function transform!(adaboost::DecisionStumpAdaboost, instances::Matrix)
  return DT.apply_adaboost_stumps(
    adaboost.model[:ensemble], adaboost.model[:coefficients], instances
  )
end

end # module


text:
import StatsBase: countmap
function raw_majority_vote(labels)
    votes = countmap(labels)
    votes = sort(collect(votes), by=x->x[2], rev=true)
    winner = votes[1][1]
    return winner
end

@assert raw_majority_vote(['a', 'b', 'c', 'b']) == 'b'

function majority_vote(labels)
    """Assumes that labels are ordered from nearest to farthest."""
    vote_counts = countmap(labels)
    vote_counts = sort(collect(vote_counts), by=x->x[2], rev=true)
    winner = vote_counts[1][1]
    if(length(vote_counts) == 1)
        return winner
    end

    winner_count= vote_counts[1][2]
    second_count = vote_counts[2][2]
    # num_winners = length([count
                    #    for count in vote_counts.values()
                    #    if count == winner_count])

    if (winner_count != second_count)
        return winner                     # unique winner, so return it
    else
        return majority_vote(labels[1:length(labels)-1]) # try again without the farthest
    end
end

# Tie, so look at first 4, then 'b'
@assert majority_vote(['a', 'b', 'c', 'b', 'a'])== 'b'

using Distances

struct LabeledPoint
    point
    label
    LabeledPoint(point, label) = new(point, label)
end

function knn_classify(k, labeled_points, new_point)

    by_distance = sort(labeled_points, by = x -> SqEuclidean()(x.point, new_point))

    k_nearest_labels = [lp.label for lp in by_distance[1:k]]

    return majority_vote(k_nearest_labels)
end

function parse_iris_row(line)
    """
    sepal_length, sepal_width, petal_length, petal_width, class
    """
    row = split(line, ",")
    measurements = [parse(Float64, value) for value in row[1:end-1]]
    # class is e.g. "Iris-virginica"; we just want "virginica"
    label = split(row[end], "-")[end]

    return LabeledPoint(measurements, label)
end

data = open(f->read(f,String), "iris.data")
irislines = []
for i in split(data, "\n")
    push!(irislines,i)
end

iris_data = [parse_iris_row(line) for line in irislines]
using Random
Random.seed!(12)
function split_data(origdata, prob)
    """Split data into fractions [prob, 1 - prob]"""
    data = copy(origdata)                    # Make a shallow copy
    data = shuffle(data)              # because shuffle modifies the list.
    cut = Int(floor(length(data) * prob))       # Use prob to find a cutoff
    return data[1:cut], data[cut+1:end]     # and split the shuffled list there.
end

iris_train, iris_test = split_data(iris_data, 0.70)
@assert length(iris_train) == 0.7*150
@assert length(iris_test) == 0.3 * 150

confusion_matrix = Dict()
let num_correct = 0

    for iris in iris_test
        if(iris.label == "virginica")
            # continue
        end
        predicted = knn_classify(5, iris_train, iris.point)
        actual = iris.label

        if predicted == actual
            num_correct += 1
        end

        confusion_matrix[(predicted, actual)] = get(confusion_matrix,(predicted, actual),0) + 1
    end

    pct_correct = num_correct / length(iris_test)

    print(pct_correct)
    print( confusion_matrix)

end


text:
using XScales
using Test

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


text:
# rotate quaternion
quaternion_rotate(q1::Quaternion,q2::Quaternion) = q2 * q1 / q2

# rotate vector
vector_rotate(v::AbstractVector,q::Quaternion) = Vmat(quaternion_rotate(Quaternion(v), q))
# vector_rotate(v::StaticVector,q::Quaternion) = q * v
∂vector_rotate∂q(p::AbstractVector, q::Quaternion) = VLmat(q) * Lmat(Quaternion(p)) * Tmat() + VRᵀmat(q) * Rmat(Quaternion(p))

# rotation matrix
rotation_matrix(q::Quaternion) = VRᵀmat(q) * LVᵀmat(q)
function ∂rotation_matrix∂q(q::Quaternion, p::AbstractVector; attjac::Bool=false)
    M = ∂VRᵀmat∂q(LVᵀmat(q) * p) + VRᵀmat(q) * ∂LVᵀmat∂q(p)
    attjac && (M *= LVᵀmat(q))
    return M
end
function ∂rotation_matrix_inv∂q(q::Quaternion, p::AbstractVector; attjac::Bool=false)
    M = ∂rotation_matrix∂q(inv(q), p) * Tmat()
    attjac && (M *= LVᵀmat(q))
    return M
end


text:
function interpp1(X,Y,V)
    itp = interpolate((X,), Y, Gridded(Linear()))
    sz = size(V)
    y=zeros(sz[1],sz[2])
    for i=1:sz[1]
        for j=1:sz[2]
            y[i,j]=itp(V[i,j])
        end
    end
    return y
end

max_stars_repo_name: umbertosaetti/J-GenHel

max_issues_repo_name: umbertosaetti/J-GenHel

max_forks_repo_name: umbertosaetti/J-GenHel


text:
struct Trip
    stop_times::Vector{StopTime}
    route::Int64
    service::Int64
    pattern::Int64
end


text:
function passthedoorman(word)
    m = match(r"(.)\1",word)
    m ≡ nothing ? 0 : 3(m.captures[1][1]-'a'+1)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using NPZ
using TomoForward
using Images
using Plots
using XfromProjections
using StaticArrays
using PyCall
using Logging
include("../../phantoms/simple_phantoms.jl")
include("../tv_primaldual_flow.jl")
include("../optical_flow.jl")
include("../../../src/analytic/gridrec.jl")

cwd = @__DIR__
cd(cwd)

like(x::T, y) where T = convert(T, y)

proj = npzread("data/projections.npy")
theta = npzread("data/theta.npy")

nangles, imageZ, detcount = size(proj)

H, W = 100,100

function radon_operator(H,W, detcount, angles)
    proj_geom = ProjGeom(1.0, detcount, angles)
    A = fp_op_parallel2d_line(proj_geom, H, W)
    return A
end

pr = Float64.(reshape(sum(proj,dims=2), (nangles, detcount)))

window_step = 2
nframes = Int(nangles/window_step)

frames = zeros(H,W,nframes)
As = map(t -> radon_operator(H,W,detcount,theta[t:t+1]),1:window_step:nangles)
bs = reshape(pr,(detcount*window_step,nframes))
niter=500
u0s = zeros(H,W,nframes)

w_tv = 0.3
w_flow  = 0.05

@info "Reconstructing using joint motion estimation and reconstruction"
c=10.0
us_flow = recon2d_tv_primaldual_flow(As, bs, u0s, 20, niter, w_tv, w_flow, c)

@info "Reconstruction using tv regularization frame by frame"
us_tv = zeros(H,W,nframes)
for t = 1:nframes
    A = As[t]
    p = bs[:,t]
    u0 = u0s[:,:,t]
    us_tv[:,:,t] .= recon2d_tv_primaldual!(us_tv[:,:,t], A, p, niter, w_tv, c)
end

@info "Preparing results in human readable format"
anim = @animate for t=1:nframes
    l = @layout [a b c]
    p1 = plot(Gray.(frames[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="Ground truth")
    p2 = plot(Gray.(us_flow[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="Flow")
    p3 = plot(Gray.(us_tv[:,:,t]), aspect_ratio=:equal, framestyle=:none, title="TV")
    plot(p1, p2, p3, layout = l)
end

cwd = @__DIR__
path = normpath(joinpath(@__DIR__, "result"))
!isdir(path) && mkdir(path)
cd(path)
gif(anim, "reconstruction_flow.gif", fps = 1)
cd(cwd)

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
import Random: GLOBAL_RNG
import Random
import SpecialFunctions: loggamma, gamma, beta
#import Random:rand

import Base.+
import Base.-
import Base.*
import Base./

+(dist::Distributions.Normal{Float64}, x::Float64) =
    Distributions.Normal(Distributions.location(dist) + x, Distributions.scale(dist))
-(dist::Distributions.Normal{Float64}, x::Float64) =
    Distributions.Normal(Distributions.location(dist) - x, Distributions.scale(dist))
*(dist::Distributions.Normal{Float64}, x::Float64) =
    Distributions.Normal(Distributions.location(dist) * x, Distributions.scale(dist) * x)
/(dist::Distributions.Normal{Float64}, x::Float64) =
    Distributions.Normal(Distributions.location(dist) / x, Distributions.scale(dist) / x)

+(dist::Distributions.Uniform{Float64}, x::Float64) =
    Distributions.Uniform(Distributions.location(dist) + x, Distributions.scale(dist))
-(dist::Distributions.Uniform{Float64}, x::Float64) =
    Distributions.Uniform(Distributions.location(dist) - x, Distributions.scale(dist))
*(dist::Distributions.Uniform{Float64}, x::Float64) =
    Distributions.Uniform(Distributions.location(dist), Distributions.scale(dist) * x)
/(dist::Distributions.Uniform{Float64}, x::Float64) =
    Distributions.Uniform(Distributions.location(dist), Distributions.scale(dist) / x)

*(dist::Distributions.LogNormal{Float64}, x::Float64) = Distributions.LogNormal(
    Distributions.location(dist) + log(x),
    Distributions.scale(dist),
)
/(dist::Distributions.LogNormal{Float64}, x::Float64) = Distributions.LogNormal(
    Distributions.location(dist) - log(x),
    Distributions.scale(dist),
)

# Truncation

function truncate_rand(
    distribution::Distributions.UnivariateDistribution,
    bounds::Vector{Float64};
    n::Int64 = 1,
)
    return map(_ -> min(max(rand(distribution), bounds[1]), bounds[2]), 1:n)
end

function truncate_rand(
    distribution::Distributions.MultivariateDistribution,
    bounds::Vector{Vector{Float64}};
    n::Int64 = 1,
)
    return map(y -> map(rand(distribution), bounds) do x, bound
        min(max(x, bound[1]), bound[2])
    end, 1:n)
end

include("polyagamma/polyagamma.jl")
include("truncated_normal/truncated_normal_2.jl")
include("gamma/gamma.jl")
include("truncated_gamma/truncated_gamma.jl")
include("truncated_inverse_gaussian/truncated_inverse_gaussian.jl")


text:
#
# Describes trajectories for all `AbstractOrbit`s.
#

"""
An alias for a `Vector` of `AbstractOrbits`.
"""
const Trajectory{T} = Vector{T} where T <: AbstractOrbit

Base.show(io::IO, traj::Trajectory{<:CartesianOrbit}) = println(io, "Cartesian Two-body trajectory with ", length(traj), " steps")
Base.show(io::IO, traj::Trajectory{<:KeplerianOrbit}) = println(io, "Keplerian Two-body trajectory with ", length(traj), " steps")
Base.show(io::IO, ::MIME"text/plain", traj::Trajectory) = show(io,traj)


text:
using Documenter
using PotentialDB, PotentialCalculation

makedocs(sitename="PotentialDB",
         pages=["Home"=>"index.md",
                "Usage"=>"usage.md"])

deploydocs(
      repo = "github.com/MatrixLabTools/PotentialDB.jl.git"
)


text:
srand(17)
N = 1000
z = 3
g =  random_regular_graph(N, z, seed=1)
Js = random_couplings(g)
σ, E = ground_state_τeo(g, Js, maxiters=10_000)
@test typeof(E) == Int
@test E == -1266

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
function evaluation_performed(obj::BatchLSObserver, arg0::jint, arg1::jint, arg2::Vector{Orbit}, arg3::ParameterDriversList, arg4::ParameterDriversList, arg5::ParameterDriversList, arg6::EstimationsProvider, arg7::LeastSquaresProblem_Evaluation)
    return jcall(obj, "evaluationPerformed", void, (jint, jint, Vector{Orbit}, ParameterDriversList, ParameterDriversList, ParameterDriversList, EstimationsProvider, LeastSquaresProblem_Evaluation), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end


text:
module src

# =================================
"""
    linear_function_interpolation( x0::Real, xsamples::Vector{T} where T <: Real, fsamples::Vector{T} where T <: Real )

linear interpolation on one-parameter function `f(x)` which is defined by a vector of sample `x` points
and a vector of sample `f(x)` corresponding to sample `x`;
`xsamples` should be increasingly sorted.
returns the function value evaluated at the given point `x0`.
if `x0<findmin(x)` or `x0>findmax(x)`, this function uses the linear expolation
of minimum/maximum sample `x` and the second minimum/maximum sample `x`.
"""
function linear_function_interpolation( x0::Real, xsamples::Vector{T} where T <: Real, fsamples::Vector{T} where T <: Real )
    issorted(xsamples) ? nothing : throw(ErrorException("xsamples should be pre-sorted by ascending"))
    local N = length(fsamples); local xmax = xsamples[end]; local xmin = xsamples[1]
    local x0_location = (x0 - xmin) / (xmax - xmin) * (N - 1) + 1  # locating x0
    local nearest_xsample_loc = floor(Int,x0_location)  # the nearest sample x point (on the left side) to interpolate x0
    local dist_from_nearest_xsample = x0_location - nearest_xsample_loc # the fractional distance between x0 and the nearest x point
    # NOTE: relationship: x[near] < x0 < x[near+1]

    # case: if touching the upper bound
    if (x0 > xmax) | (nearest_xsample_loc == N)  # to prevent BoundsError at index N
        res = fsamples[end]
        return res::Float64
    # case: if touching the lower bound
    elseif (x0 <= xmin) | (nearest_xsample_loc == 1)
        res = fsamples[1] - (1 - x0_location) * ( fsamples[2] - fsamples[1] )
        return res::Float64
    # case: normal case
    else
        res = (1 - dist_from_nearest_xsample) * fsamples[nearest_xsample_loc] +
            dist_from_nearest_xsample * fsamples[nearest_xsample_loc + 1]
        return res::Float64
    end # if

end # linear_function_interpolation

# =================================
"""
    utilfunc_1var(x::Real ; type::String = "logarithm", pars::Union{Tuple,NamedTuple} = ()::Tuple{} )

a dispatch of different kinds of utility functions with one argument (e.g. consumption).
supported utility functions & required parameters:
1. "logarithm" + () <- empty tuple <- ``u(x)=\\log(x)``
2. "CRRA" + (η = ,) <- relative risk aversion coefficient <- ``u(x)=\\frac{x^{1-\\eta}-1}{1-\\eta}``
3. "CARA" + (ρ = ,) <- absolute risk aversion coefficient <- ``u(x)=-e^{-\\rho x}``
4. "CD" + (θ = ,) <- share of consumption <- ``u(x)=x^{\\theta}``
"""
function utilfunc_1var(x::Real ; type::String = "logarithm", pars::Union{Tuple,NamedTuple} = ()::Tuple{} )
    if type == "logarithm"
        return u_log_1var(x)

    return nothing
end # utilfunc_1var

end # src
#

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
using Test

@testset "0 degree rotation is and identity" begin
    op = vertex_map(rotation(0), 10, 5)

    @test op((1, 1)) == (1, 1)
    @test op((10, 5)) == (10, 5)
    @test op((1, 5)) == (1, 5)
    @test op((10, 1)) == (10, 1)
    @test op((4, 3)) == (4, 3)
    @test op((2, 2)) == (2, 2)
end

@testset "90 degree rotation of square lattice correctly transforms vertices" begin
    op = vertex_map(rotation(90), 10, 5)

    @test op((1, 1)) == (10, 1)
    @test op((5, 10)) == (1, 5)
    @test op((1, 10)) == (1, 1)
    @test op((5, 1)) == (10, 5)
    @test op((4, 3)) == (8, 4)
    @test op((3, 4)) == (7, 3)
    @test op((2, 2)) == (9, 2)
end

@testset "180 degree rotation of square lattice correctly transforms vertices" begin
    op = vertex_map(rotation(180), 10, 5)

    @test op((1, 1)) == (10, 5)
    @test op((10, 5)) == (1, 1)
    @test op((1, 5)) == (10, 1)
    @test op((10, 1)) == (1, 5)
    @test op((4, 3)) == (7, 3)
end

@testset "270 degree rotation of square lattice correctly transforms vertices" begin
    op = vertex_map(rotation(270), 10, 5)

    @test op((1, 1)) == (1, 5)
    @test op((5, 1)) == (1, 1)
    @test op((5, 10)) == (10, 1)
    @test op((1, 10)) == (10, 5)
    @test op((4, 3)) == (3, 2)
    @test op((3, 4)) == (4, 3)
end

@testset "Reflection around x axis correctly transforms vertices" begin
    op = vertex_map(reflection(:x), 10, 5)

    @test op((1, 1)) == (10, 1)
    @test op((10, 5)) == (1, 5)
    @test op((1, 5)) == (10, 5)
    @test op((10, 1)) == (1, 1)
    @test op((4, 3)) == (7, 3)
    @test op((2, 2)) == (9, 2)
end

@testset "Reflection around y axis correctly transforms vertices" begin
    op = vertex_map(reflection(:y), 10, 5)

    @test op((1, 1)) == (1, 5)
    @test op((10, 5)) == (10, 1)
    @test op((1, 5)) == (1, 1)
    @test op((10, 1)) == (10, 5)
    @test op((4, 3)) == (4, 3)
    @test op((2, 2)) == (2, 4)
end

@testset "Reflection around diag correctly transforms vertices" begin
    op = vertex_map(reflection(:diag), 10, 5)

    @test op((1, 1)) == (1, 1)
    @test op((5, 10)) == (10, 5)
    @test op((1, 10)) == (10, 1)
    @test op((5, 1)) == (1, 5)
    @test op((4, 3)) == (3, 4)
    @test op((2, 2)) == (2, 2)
end

@testset "Reflection around antydiag correctly transforms vertices" begin
    op = vertex_map(reflection(:antydiag), 10, 5)

    @test op((1, 1)) == (10, 5)
    @test op((5, 10)) == (1, 1)
    @test op((1, 10)) == (1, 5)
    @test op((5, 1)) == (10, 1)
    @test op((4, 3)) == (8, 2)
    @test op((3, 4)) == (7, 3)
    @test op((2, 2)) == (9, 4)
end

for transform ∈ all_lattice_transformations
@testset "$(transform) of square lattice is a bijection" begin
    op = vertex_map(transform, 3, 3)
    all_points = [(i, j) for i ∈ 1:3, j ∈ 1:3]
    @test Set(op.(all_points)) == Set(all_points)
end
end

for transform ∈ (rotation.([90, 270])..., reflection.([:diag, :antydiag])...)
@testset "$(transform) of rectangular lattice is bijection flipping lattice dimensions" begin
    op = vertex_map(transform, 2, 3)
    original_points = [(i, j) for i ∈ 1:3, j ∈ 1:2]
    transformed_points = [(i, j) for i ∈ 1:2, j ∈ 1:3]
    @test Set(op.(original_points)) == Set(transformed_points)
end
end

for transform ∈ (rotation.([0, 180])..., reflection.([:x, :y])...)
@testset "$(transform) of rectangular lattice is bijection leaving lattice dimensions unchanged" begin
    op = vertex_map(transform, 7, 3)
    all_points = [(i, j) for i ∈ 1:7, j ∈ 1:3]
    @test Set(op.(all_points)) == Set(all_points)
end
end


text:
function mps(s)
# create minimum-phase spectrum sm from complex spectrum s
   return exp( fft( fold( ifft( log( clipdb(complex(s),-100) )))));
end


text:
### File for the run! function for Stacks
# Compressed Stack structure is efficient within this run method

### Functions to run a CompressedStack similarly than for a classic stack
function run!(stack::CompressedStack)
  while !eof(stack.input)
    elt = readinput(stack)
    while !isempty(stack) && stack.pop_condition(stack,elt)
      pop!(stack)
    end
    if stack.push_condition(stack, elt)
      push!(stack, elt)
    end
  end
end
# run! function specific to the reconstruction procedure of Compressed Stacks
function run!(stack::CompressedStack, limit::Int)
  while limit >= 0
    elt = readinput(stack)
    limit -= 1
    while !isempty(stack) && stack.pop_condition(stack,elt)
      pop!(stack)
    end
    if stack.push_condition(stack, elt)
      push!(stack, elt)
    end
  end
end

### Functions to run a CompressedStack with extra control on the pop condition
function runCountingPop!(stack::CompressedStack)
  while !eof(stack.input)
    elt = readinput(stack)
    popCounter=0
    while !isempty(stack) && stack.pop_condition(stack,elt,popCounter)
      pop!(stack)
      popCounter=popCounter+1
    end
    print(stack)
    println("LEAVING POP LOOOP!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ",elt, " empty ",isempty(stack))
    #print(stack)
    if stack.push_condition(stack, elt)
      println("pushing ",elt)
      push!(stack, elt)
      print(stack)
    end
  end
end
# run! function specific to the reconstruction procedure of Compressed Stacks
function runCountingPop!(stack::CompressedStack, limit::Int)
  println("CALLED RUN COUNTING POP!")
  while limit >= 0
    elt = readinput(stack)
    limit -= 1
    popCounter=0
    while !isempty(stack) && stack.pop_condition(stack,elt,popCounter)
      pop!(stack)
      popCounter=popCounter+1
    end
    if stack.push_condition(stack, elt)
      push!(stack, elt)
    end
  end
end

### For Normal Stacks
function run!(stack::NormalStack)
  while !eof(stack.input)
    elt = readinput(stack)
    while !isempty(stack) && stack.pop_condition(stack,elt)
      pop!(stack)
    end
    if stack.push_condition(stack, elt)
      push!(stack, elt)
    end
  end
end

### Comparison between both Stacks
function run!(cs::CompressedStack, ns::NormalStack)
  println("\n\n!!!Start of the comparative run!!!\n")
  while !eof(cs.input)
    if eof(ns.input) != eof(cs.input)
      println("eof_ns = $(eof(ns.input)); eof_cs = $(eof(cs.input))")
      println("pos_ns = $(position(ns.input)); pos_cs = $(position(cs.input))")
      print(cs)
      print(ns)
      error("The two stacks have different input reading: eof")
    elseif isempty(cs) != isempty(ns)
      error("Only one of the two stacks is empty")
    end
    elt = readinput(cs)
    if elt != readinput(ns)
      error("The content in the input file are read differently before a push!")
    elseif cs.pop_condition(cs,elt) != ns.pop_condition(ns,elt)
      error("The pop conditions of the stacks give different boolean values")
    end

    while !isempty(cs) && cs.pop_condition(cs,elt)
      if top(ns) != top(cs)
              print(cs)
              print(ns)
        println("top: cs=$(top(cs)), ns=$(top(ns))")
        error("The top element of the stacks are different after a pop!")
      end
      pop!(cs)
      pop!(ns)
    end

    if cs.push_condition(cs,elt) != ns.push_condition(ns,elt)
      error("The push conditions of the stacks give different boolean values")
    end
    if cs.push_condition(cs, elt)
      push!(cs, elt)
      push!(ns, elt)
    end
  end
  println("\n\n!!!End of the comparative run!!!\n")
  print(cs)
  print(ns)
end


text:
module Microclimate

using Unitful,
      NCDatasets,
      FieldMetadata,
      FieldDefaults

using Unitful: W, m, °C, hr, mol, K, s, J, Mg, g, kg, L, kPa, Pa

import FieldMetadata: default, units, limits, @units, @limits

export AbstractMicroclimate,
       AbstractLayeredMicroclimate,
       AbstractMicroclimGrid, MicroclimGrid,
       AbstractMicroclimPoint, MicroclimPoint,
       AbstractMicroclimInstant, MicroclimInstant,
       AbstractMicroclimControl, MicroclimControl,
       LinearLayerInterpolator, LinearLayerInterpolators

export radiation, snowdepth, airtemperature, relhumidity, windspeed,
       soiltemperature, soilwaterpotential, soilwatercontent

export mean_soiltemperature, mean_soilwaterpotential, mean_soilwatercontent

export weightedmean, layermax

export load_grid, load_point

export download_microclim

include("interpolation.jl")
include("types.jl")
include("netcdf.jl")
include("files.jl")

end # module


text:
"""
`SatResample` -- conventiently resample satellite data to common local grids.
"""
module SatResample

# Write your package code here.
using Proj4
using NearestNeighbors
using Statistics
using StaticArrays

include("sample_grids.jl")

"""
	resample(mx, my, data, xgrid, ygrid, maxdist, nneighbors)

returns unstructured ungridded data on mx,my coordinates onto
a regular 2D gridded data on xgrid, ygrid
"""
function resample(mx,my,data,xgrid,ygrid,maxdist,nneighbors)

	if isa(data[1],AbstractArray) #multiple data layers
		md=true
	else
		md=false
	end

	AA=[SVector(i) for i in Iterators.product(ygrid,xgrid)][:]
	BB=SVector.(zip(my,mx))
	kdt=KDTree(BB)#;leafsize=3,reorder=false)
	ind,dist=knn(kdt,AA,nneighbors,true)

	if md
		outimg= [fill(NaN,(length(ygrid),length(xgrid))) for i in 1:length(data)]
	else
		outimg=fill(NaN,(length(ygrid),length(xgrid)))
	end

	for (n,(i,d)) in enumerate(zip(ind,dist))
		c=0
		for di in d
			if di<maxdist
				c+=1
			else
				break
			end
		end
		if c!=0
			if md
				for ii = 1:length(data)
					outimg[ii][n]=mean(@views data[ii][i[1:c]])
				end
			else
				outimg[n]=mean(@views data[i[1:c]])
			end
		end
	end
	return outimg
end

"""
	Grid(proj, xlim, ylim, resolution)

Grid instance holding projection, border information and resolution
"""
struct Grid
	"The projection which is a Proj4 instance"
	proj::Proj4.Projection
	"A tuple holding minimum and maximum values in x direction in m"
	xlim::NTuple{2,Float64}
	"A tuple holding minimum and maximum values in x direction in m"
	ylim::NTuple{2,Float64}
	"the resolution in m"
	resolution::Float64
end

"""
	gen_ps_grid(clon, clat, width, height, resolution)

generates polar stereographic grid with center coordinates of `clat`, `clon`
and `width` and `height` and `resolution` in m. This is a special case with
the North pole always rotated upwards. For arbitrary rotated grid see [`gen_ps_grid`](@ref)
"""
function gen_ps_grid(clon,clat,width,height,resolution)
	proj=Projection("+proj=stere +lat_0=$clat +lon_0=$clon +ellps=WGS84")
	Grid(proj,(-width/2, width/2), (-height/2, height/2),resolution)
end

"""
	gen_ps_grid(clon, clat, width, height, resolution, rotation)

generates polar stereographic grid with center coordinates of `clat`, `clon`
and `width` and `height` and `resolution` in m rotated by `rotation` degrees`.
"""
function gen_ps_grid(clon,clat,width,height,resolution,rotation)
	proj=Projection("+proj=stere +lat_0=$clat +lat_ts=$clat +lon_0=$(rotation) +ellps=WGS84")
	x,y=lonlat2xy(Float64[clon clat],proj)
	nproj=Projection("+proj=stere +lat_0=$clat +lat_ts=$clat +lon_0=$rotation +x_0=$(-x) +y_0=$(-y) +ellps=WGS84")
	return Grid(nproj,(-width/2, width/2), (-height/2, height/2),resolution)
end

"""
	resample(mx, my, data, grid, maxdist, nneighbor)

resample `data` which is already in `grid` coordinates `mx` and `my`.
`maxdist` is the maximum distance in grid units to take into account
`nneighbor` is the number of neighbors to take into account for the resampling
"""
function resample(mx,my,data,grid::Grid,maxdist,nneighbor)
	xg,yg=genxy(grid)
	resample(mx,my,data,xg,yg,maxdist,nneighbor)
end

"""
	resample(lon, lat, data, grid, maxdist, nneighbor)

resample `data` which is `lat`, `lon` coordinates into the `grid`.
`maxdist` is the maximum distance in grid units to take into account
`nneighbor` is the number of neighbors to take into account for the resampling
"""
function llresample(lon,lat,data,grid::Grid,maxdist,nneighbor)
	xx,yy=xyfromlonlat(lon,lat,grid)
	resample(xx,yy,data,grid,maxdist,nneighbor)
end

"""
	xyfromlonlat(lon, lat, grid)

convert `lat`, `lon` into x and y coordinates of the `grid`
"""
function xyfromlonlat(lon,lat,grid::Grid)
	eachcol(lonlat2xy(Float64[lon lat],grid.proj))
end

"""
	genxy(grid)

generates grid coordinates `x` and `y` for a given grid instance
"""
function genxy(grid::Grid)
	xg=grid.xlim[1]+grid.resolution/2:grid.resolution:grid.xlim[2]
	yg=grid.ylim[1]+grid.resolution/2:grid.resolution:grid.ylim[2]
	return xg,yg
end

end


text:
import Base.==

const no_attr_uuid = string(Base.Random.uuid4())

getattr_safe(n :: HTMLText, a :: AbstractString) = no_attr_uuid
getattr_safe(n :: HTMLNode, a :: AbstractString) = get(n.attributes, a, no_attr_uuid)

tag_safe(n :: HTMLText) = no_attr_uuid
tag_safe(n :: HTMLNode) = tag(n)

abstract type GumboPred end

immutable IsTag <: GumboPred; tag; end
IsTag(t :: AbstractString) = IsTag(Symbol(t))
matchp(p :: IsTag, n :: HTMLNode) = tag_safe(n) == p.tag
==(a :: IsTag, b :: IsTag) = a.tag == b.tag

immutable HasId <: GumboPred; id; end
matchp(p :: HasId, n :: HTMLNode) = getattr_safe(n, "id") == p.id
==(a :: HasId, b :: HasId) = a.id == b.id

immutable HasClass <: GumboPred; regex :: Regex; end
HasClass(s :: AbstractString) = HasClass(Regex("\\b$s\\b"))
matchp(p :: HasClass, n :: HTMLNode) = ismatch(p.regex, getattr_safe(n, "class"))
==(a :: HasClass, b :: HasClass) = a.regex == b.regex

immutable MatchGroup <: GumboPred; ps :: Array{GumboPred}; end
function matchp(ps :: MatchGroup, n :: HTMLNode)
    for p=ps.ps
        if !matchp(p, n)
            return false
        end
    end
    return true
end
==(a :: MatchGroup, b :: MatchGroup) = a.ps == b.ps

function _css_part_to_matcher(desc)
    function matcher(e)
        if e[1] == '#'
            HasId(e[2:end])
        elseif e[1] == '.'
            HasClass(e[2:end])
        else
            IsTag(e)
        end
    end

    preds = [matcher(e) for e=matchall(r"(#[\w-]+)|(\w+)|(\.[\w-]+)", desc)]
    if length(preds) == 1
        preds[1]
    else
        MatchGroup(preds)
    end
end

function compile_css(css)
    ops = split(strip(css), r"\s+")
    parts = GumboPred[_css_part_to_matcher(e) for e=ops]

    if length(parts) == 1
        parts[1]
    else
        parts
    end
end

macro css_str(css)
    quote
        $(compile_css(css))
    end
end

import Base.select
function select(p :: GumboPred, n :: HTMLNode, out = HTMLNode[] :: Array{HTMLNode})
    if matchp(p, n)
        push!(out, n)
    end
    for e=children(n)
        select(p, e, out)
    end
    out
end

function select(p :: GumboPred, ns :: Array{HTMLNode}, out = HTMLNode[] :: Array{HTMLNode})
    for n=ns
        select(p, n, out)
    end
    out
end

function select(ps :: Array{GumboPred}, items :: Array{HTMLNode}, out = HTMLNode[] :: Array{HTMLNode})
    @assert length(ps) > 0

    let p=ps[1]
        next_items = HTMLNode[]
        for e=items
            select(p, e, next_items)
        end
        items = next_items
    end

    if length(ps) == 1
        return items
    end

    for p=ps[2:end]
        next_items = HTMLNode[]
        for e=items
            for c=children(e)
                select(p, c, next_items)
            end
        end
        items = next_items
    end
    items
end

select(ps :: Array{GumboPred}, n :: HTMLNode, out = HTMLNode[] :: Array{HTMLNode}) =
    select(ps, [n], out)

select(p :: GumboPred, d :: HTMLDocument, out = HTMLNode[] :: Array{HTMLNode}) =
    select(p, children(d.root), out)

select(ps :: Array{GumboPred}, d :: HTMLDocument, out = HTMLNode[] :: Array{HTMLNode}) =
    select(ps, children(d.root), out)


text:
#******************************************************************************#
#                                                                              #
#                                                         :::      ::::::::    #
#    decode.jl                                          :+:      :+:    :+:    #
#                                                     +:+ +:+         +:+      #
#    By: mngomane <mngomane@student.42.fr>          +#+  +:+       +#+         #
#                                                 +#+#+#+#+#+   +#+            #
#    Created: 2015/03/13 22:44:33 by mngomane          #+#    #+#              #
#    Updated: 2015/03/13 22:44:33 by mngomane         ###   ########.fr        #
#                                                                              #
#******************************************************************************#

args = length(ARGS)

if	args == 1
	index = 1;

	tab = split(ARGS[1], ".")
	len = length(tab)
	while index <= len
		print(char(int(tab[index])))
		index += 1
	end
	print("\n")
end

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
import DiffEqJump

struct SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}
    integrator::DiffEqJump.SSAIntegrator{F,uType,tType,P,S,CB,SA,OPT,TS}
end

Base.IteratorSize(::Type{SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}}) where {F,uType,tType,P,S,CB,SA,OPT,TS} = Base.SizeUnknown()
Base.eltype(::Type{SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}}) where {F,uType,tType,P,S,CB,SA,OPT,TS} = Tuple{uType,tType,Int}

function Base.iterate(iter::SSAIter)
    integrator = iter.integrator
    end_time = integrator.sol.prob.tspan[2]
    if integrator.t >= end_time
        return nothing
    end
    aggregator = integrator.cb.condition
    i = aggregator.next_jump

    next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)
    if !isempty(integrator.tstops) &&
       integrator.tstops_idx <= length(integrator.tstops) &&
       integrator.tstops[integrator.tstops_idx] < next_jump_time
        next_jump_time = integrator.tstops[integrator.tstops_idx]
        i = 0
    end
    t = min(next_jump_time, end_time)
    if t == end_time
        i = 0
    end
    return (integrator.u, t, i), ()
end

function Base.iterate(iter::SSAIter, state::Tuple{})
    integrator = iter.integrator
    if DiffEqJump.should_continue_solve(integrator)
        end_time = integrator.sol.prob.tspan[2]
        step!(integrator)
        aggregator = integrator.cb.condition
        i = aggregator.next_jump

        next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)
        if !isempty(integrator.tstops) &&
           integrator.tstops_idx <= length(integrator.tstops) &&
           integrator.tstops[integrator.tstops_idx] < next_jump_time
            next_jump_time = integrator.tstops[integrator.tstops_idx]
            i = 0
        end
        t = min(next_jump_time, end_time)
        if t == end_time
            i = 0
        end

        return (integrator.u, t, i), ()
    end
    nothing
end

SciMLBase.getsyms(iter::SSAIter) = SciMLBase.getsyms(iter.integrator.sol)
SciMLBase.getindepsym(iter::SSAIter) = SciMLBase.getindepsym(iter.integrator.sol)


text:
module SFGReservoirs

using DifferentialEquations
using SFGAnalysis
using StaticArrays
using Optim

include("sfspectra.jl")
include("reservoirs.jl")

export  model,
        State

end # module