text:
"""
```
classicBP :: Function
```
This function is an important solver of this package. This solver implements the algorithm given in

Liberti, L., Lavor, C., & Maculan, N. (2008). A branch‐and‐prune algorithm for the molecular distance geometry problem. International Transactions in Operational Research, 15(1), 1-17.

The main difference of our implementation is that the input data (`NMRType`) can store informations of preprocessing function and as consequence we can optimize the search tree of this algorithm.
"""
function classicBP(NMRdata :: NMRType,
		   ε :: Float64,
		   virtual_ε :: Float64,
		  allmol :: Bool)
	# defining closure
	function classicBP_closure(l :: Int64,
				   pos::Int64,
				   mol :: MoleculeType,
				   C :: Array{Float64,2})
		if l == 1
			# first atom
			mol.atoms[1].element = NMRdata.info[1,:].nzval[1].atom1
			mol.atoms[1].x = 0.0
			mol.atoms[1].y = 0.0
			mol.atoms[1].z = 0.0
			#second atom
			mol.atoms[2].element = NMRdata.info[2,:].nzval[1].atom1
			mol.atoms[2].x = -NMRdata.info[1,2].dist
			mol.atoms[2].y = 0.0
			mol.atoms[2].z = 0.0
			# tird atom
			D12 = NMRdata.info[1,2].dist
			D13 = NMRdata.info[1,3].dist
			D23 = NMRdata.info[2,3].dist
			D14 = 0.0
			D24 = 0.0
			D34 = 0.0
			cθ,sθ = bondangle(D12,D13,D23)
			nop_node += [3,6,1,1]
			cω,sω = (0.0,0.0)
			mol.atoms[3].element = NMRdata.info[3,:].nzval[1].atom1
			mol.atoms[3].x = -D12+D23*cθ
			mol.atoms[3].y = D23*sθ
			mol.atoms[3].z = 0.0
			nop_node += [1,2,0,0]
			C = zeros(4,4)
			C[1,4] = mol.atoms[3].x
			C[2,4] = mol.atoms[3].y
			C[3,4] = mol.atoms[3].z
			C[1,1] = cθ
			C[1,2] = sθ
			C[2,1] = sθ
			C[2,2] = -cθ
			C[3,3] = -1.0
			C[4,4] = 1.0
			l = 4 # branching starts at atom 4
			pos = 4 # position in virtual path
		end
		B = zeros(4,4)
		λ = 1
		ρ = 1
		C_before = zeros(4,4)
		copyto!(C_before,C)
		keep = true
		while keep
			try
				D14 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos]].dist
			catch
				D14 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D24 = NMRdata.info[NMRdata.virtual_path[pos-2],NMRdata.virtual_path[pos]].dist
			catch
				D24 = sqrt((mol.atoms[NMRdata.virtual_path[pos-2]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-2]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-2]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D34 = NMRdata.info[NMRdata.virtual_path[pos-1],NMRdata.virtual_path[pos]].dist
			catch
				D34 = sqrt((mol.atoms[NMRdata.virtual_path[pos-1]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-1]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-1]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D12 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos-2]].dist
			catch
				D12 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos-2]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos-2]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos-2]].z)^2)
			end
			try
				D13 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos-1]].dist
			catch
				D13 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos-1]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos-1]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos-1]].z)^2)
			end
			try
				D23 = NMRdata.info[NMRdata.virtual_path[pos-2],NMRdata.virtual_path[pos-1]].dist
			catch
				D23 = sqrt((mol.atoms[NMRdata.virtual_path[pos-2]].x - mol.atoms[NMRdata.virtual_path[pos-1]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-2]].y - mol.atoms[NMRdata.virtual_path[pos-1]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-2]].z - mol.atoms[NMRdata.virtual_path[pos-1]].z)^2)
			end
			cθ,sθ = bondangle(D23,D24,D34)
			cω,sω = badtorsionangle(D12,D13,D14,D23,D24,D34)
			@debug "l value = $(l) and NMRdatavalue = $(NMRdata.virtual_path[pos]) in position $(pos)"
			B = torsionmatrix(cθ,sθ,cω,sω,D34,B,true)
			if l==NMRdata.virtual_path[pos]
				nop_node += [0,7,0,0] #torsion matrix
				nop_node += [3,6,1,1] # bond angle
				nop_node += [10,20,4,2] # bad torsion angle
				C = prodmatrix(C_before,B)
				nop_node += [24,33,0,0]
				keep = false
			else
				nop_node += [0,7,0,0] #torsion matrix
				nop_node += [3,6,1,1] # bond angle
				nop_node += [10,20,4,2] # bad torsion angle

				nop_vpath += [0,7,0,0] # torsion matrix
				nop_vpath += [3,6,1,1] # bond angle
				nop_vpath += [10,20,4,2] # bad torsion angle

				cpx = mol.atoms[NMRdata.virtual_path[pos]].x
				cpy = mol.atoms[NMRdata.virtual_path[pos]].y
				cpz = mol.atoms[NMRdata.virtual_path[pos]].z

				Virtual_Torsion = prodmatrix(C_before,B)
				nop_vpath += [24,33,0,0]
				nop_node += [23,33,0,0]

				if sqrt((Virtual_Torsion[1,4]- cpx)^2+(Virtual_Torsion[2,4]- cpy)^2+(Virtual_Torsion[3,4]- cpz)^2)> virtual_ε
					B = torsionmatrix(cθ,sθ,cω,sω,D34,B,false)
					C_before = prodmatrix(C_before,B)
					nop_vpath += [24,33,0,0]
					nop_node += [24,33,0,0]

				else
					copyto!(C_before,Virtual_Torsion)
				end
				@debug "virtual atom position  " C_before[1,4],C_before[2,4],C_before[3,4]
				pos = pos+1
			end
		end

		mol.atoms[l].element = NMRdata.info[l,:].nzval[1].atom1
		mol.atoms[l].x = C[1,4]
		mol.atoms[l].y = C[2,4]
		mol.atoms[l].z = C[3,4]
		count = [0,0,0,0]
		λ , count  = pruningtest(mol,l,NMRdata,ε,count)
		nop_ddf += count
		@debug "C at level $(l) right side " C
		if λ == 1
			if l<n
				@debug "Partial solution by right side at level $(l)"  mol
				n_branch +=1
				classicBP_closure(l+1,pos+1,mol,C)
			else
				nsol=nsol+1
				storage_mol[nsol] = copy(mol)
				@debug "Rank n was reached, a solution was found "
				return 0
			end
		else
			n_prune += 1
		end
		if allmol==false && nsol>0
			@debug "number of solutions"  nsol
			#@info "LDE = " LDE(mol,D,n,nad)
			@goto exit
		end
		B = torsionmatrix(cθ,sθ,cω,sω,D34,B,false)
		#nop_node += [0,0,0,0]
		C = prodmatrix(C_before,B)# tenho que otimizar este calculo
		nop_node += [24,33,0,0]
		mol.atoms[l].x = C[1,4]
		mol.atoms[l].y = C[2,4]
		mol.atoms[l].z = C[3,4]
		count = [0,0,0,0]
		ρ ,count = pruningtest(mol,l,NMRdata,ε,count) #preciso modificar
		nop_ddf += count
		@debug "C at level $(l) left side " C
		if ρ == 1
			if l<n
				@debug "Partial solution by left side at level $(l)" mol
				n_branch += 1
				classicBP_closure(l+1,pos+1,mol,C)
			else
				nsol = nsol+1
				storage_mol[nsol] = copy(mol)
				@debug "Rank n was reached, a solution was found "
				return 0
			end
		else
			n_prune += 1
		end

		@label exit
		return 0
	end # closure
	n = NMRdata.dim
	mol = MoleculeType(Vector{AtomType}(undef,n),0.0)
	for i=1:n
		mol.atoms[i] = AtomType(0.0,0.0,0.0)
	end
	C = zeros(4,4)
	nsol = 0
	storage_mol = Dict{Int64,MoleculeType}()
	n_prune = 0
	n_branch = 0
	#count_nop = [+-,*,/,√]
	nop_node = [0,0,0,0]
	nop_ddf = [0,0,0,0]
	nop_vpath = [0,0,0,0]
	classicBP_closure(1,1,mol,C)
	return nsol, storage_mol,Counter(nop_node,nop_vpath,nop_ddf,n_branch,n_prune)

end #solver classicBP

"""
```
quaternionBP :: Function
```
This function defines a new solver. The implementation follows ideas describing in:

Fidalgo, F. Using Quaternion Geometric Algebra for efficient rotations in Branch and Prune Algorithtm to solve the Discretizable Molecular Distance Geometry Problem. In: Proceedings of AGACSE 2018, Campinas-SP, Brazil.
"""
#qbondangle = [4,5,2,2]
#qtorsionangle = [11,19]
function quaternionBP(NMRdata :: NMRType,
		      ε :: Float64,
		      virtual_ε :: Float64,
		     allmol :: Bool)
	# defining closure
	function quaternionBP_closure(l :: Int64,
				      pos::Int64,
				      mol :: MoleculeType,
				      Q :: Quaternion)
		if l == 1
			# first atom
			mol.atoms[1].element = NMRdata.info[1,:].nzval[1].atom1
			mol.atoms[1].x = 0.0
			mol.atoms[1].y = 0.0
			mol.atoms[1].z = 0.0
			#second atom
			mol.atoms[2].element = NMRdata.info[2,:].nzval[1].atom1
			mol.atoms[2].x = -NMRdata.info[1,2].dist
			mol.atoms[2].y = 0.0
			mol.atoms[2].z = 0.0
			#Q_before = Quaternion(0.0,0.0,-1.0,0.0)
			# tird atom
			D12 = NMRdata.info[1,2].dist
			D13 = NMRdata.info[1,3].dist
			D23 = NMRdata.info[2,3].dist
			D14 = 0.0
			D24 = 0.0
			D34 = 0.0
			cθ,sθ = qbondangle(D12,D13,D23)
			nop_node += [4,5,1,2]
			#Q = qprod(Q_before,Quaternion(sθ,0.0,0.0,cθ))
			#nop_node += [4,8,0,0]
			#qmol = rotopt(Q,D23)
			#nop_node += [4,10,0,0]
			Q = Quaternion(0.0,-cθ,-sθ,0.0)
			d = 2.0*D23
			qmol = Quaternion(0.0,d*(cθ^2-0.5),d*(cθ*sθ),0.0)
			nop_node = [1,4,0,0]
			mol.atoms[3].element = NMRdata.info[3,:].nzval[1].atom1
			mol.atoms[3].x = qmol.v1 + mol.atoms[2].x
			mol.atoms[3].y = qmol.v2 + mol.atoms[2].y
			mol.atoms[3].z = qmol.v3 + mol.atoms[2].z
			nop_node += [3,0,0,0]
			l = Int64(4) # branching starts at atom 4
			pos = Int64(4) # position in virtual path
		end
		lastpos = 1
		λ = 1
		ρ = 1
		a = 0.0
		b = 0.0
		c = 0.0
		d = 0.0
		Q_before = Q
		keep = true
		while keep
			try
				D14 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos]].dist
			catch
				D14 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D24 = NMRdata.info[NMRdata.virtual_path[pos-2],NMRdata.virtual_path[pos]].dist
			catch
				D24 = sqrt((mol.atoms[NMRdata.virtual_path[pos-2]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-2]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-2]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D34 = NMRdata.info[NMRdata.virtual_path[pos-1],NMRdata.virtual_path[pos]].dist
			catch
				D34 = sqrt((mol.atoms[NMRdata.virtual_path[pos-1]].x - mol.atoms[NMRdata.virtual_path[pos]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-1]].y - mol.atoms[NMRdata.virtual_path[pos]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-1]].z - mol.atoms[NMRdata.virtual_path[pos]].z)^2)
			end
			try
				D12 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos-2]].dist
			catch
				D12 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos-2]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos-2]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos-2]].z)^2)
			end
			try
				D13 = NMRdata.info[NMRdata.virtual_path[pos-3],NMRdata.virtual_path[pos-1]].dist
			catch
				D13 = sqrt((mol.atoms[NMRdata.virtual_path[pos-3]].x - mol.atoms[NMRdata.virtual_path[pos-1]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-3]].y - mol.atoms[NMRdata.virtual_path[pos-1]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-3]].z - mol.atoms[NMRdata.virtual_path[pos-1]].z)^2)
			end
			try
				D23 = NMRdata.info[NMRdata.virtual_path[pos-2],NMRdata.virtual_path[pos-1]].dist
			catch
				D23 = sqrt((mol.atoms[NMRdata.virtual_path[pos-2]].x - mol.atoms[NMRdata.virtual_path[pos-1]].x)^2 + (mol.atoms[NMRdata.virtual_path[pos-2]].y - mol.atoms[NMRdata.virtual_path[pos-1]].y)^2+ (mol.atoms[NMRdata.virtual_path[pos-2]].z - mol.atoms[NMRdata.virtual_path[pos-1]].z)^2)
			end
			cθ,sθ = qbondangle(D23,D24,D34)
			nop_node += [4,5,1,2]
			cω,sω = qtorsionangle(D12,D13,D14,D23,D24,D34)
			nop_node += [11,16,1,3]
			a = sθ*cω
			b = sθ*sω
			c = -cθ*sω
			d = cθ*cω
			nop_node += [0,4,0,0]
			@debug "l value = $(l) and NMRdatavalue = $(NMRdata.virtual_path[pos]) in position $(pos)"
			if l==NMRdata.virtual_path[pos]
				Q = qprod(Q_before,Quaternion(a,b,c,d))
				keep = false
				lastpos = pos
				nop_node += [12,16,0,0]
			else
				nop_vpath += [19,21,2,5] #bond+torsion+product
				Q_virtual = qprod(Q_before,Quaternion(a,b,c,d))
				nop_vpath += [12,16,0,0]
				nop_node += [12,16,0,0]
				qmol = rotopt(Q_virtual,D34)
				nop_vpath += [4,10,0,0]
				nop_node += [4,10,0,0]
				vx = qmol.v1 + mol.atoms[NMRdata.virtual_path[pos-1]].x-mol.atoms[NMRdata.virtual_path[pos]].x
				vy = qmol.v2 + mol.atoms[NMRdata.virtual_path[pos-1]].y-mol.atoms[NMRdata.virtual_path[pos]].y
				vz = qmol.v3 + mol.atoms[NMRdata.virtual_path[pos-1]].z-mol.atoms[NMRdata.virtual_path[pos]].z
				nop_node += [3,0,0,0]
				nop_vpath += [3,0,0,0]
				if sqrt(vx^2+vy^2+vz^2)> virtual_ε
					Q_before = qprod(Q_before,Quaternion(a,-b,-c,d))
					nop_node += [12,16,0,0]
					nop_vpath +=[12,16,0,0]
				else
					Q_before = Q_virtual
				end
				#qmol_aux = rot(Q_before,D34) #? is 34?
				#vx = qmol_aux.v1 + mol.atoms[NMRdata.virtual_path[pos-1]].x
				#vy = qmol_aux.v2 + mol.atoms[NMRdata.virtual_path[pos-1]].y
				#vz = qmol_aux.v3 + mol.atoms[NMRdata.virtual_path[pos-1]].z
				#@debug "repeated atom = $(vx),$(vy),$(vz)"
				pos = pos+1
			end
		end
		qmol = rotopt(Q,D34)
		nop_node += [4,10,0,0]
		mol.atoms[l].element = NMRdata.info[l,:].nzval[1].atom1
		mol.atoms[l].x = qmol.v1 + mol.atoms[NMRdata.virtual_path[lastpos-1]].x
		mol.atoms[l].y = qmol.v2 + mol.atoms[NMRdata.virtual_path[lastpos-1]].y
		mol.atoms[l].z = qmol.v3 + mol.atoms[NMRdata.virtual_path[lastpos-1]].z
		nop_node += [3,0,0,0]
		@debug "candidate atom by right side at level $(l) = $(mol.atoms[l].x),$(mol.atoms[l].y),$(mol.atoms[l].z)"
		count = [0,0,0,0]
		λ , count  = pruningtest(mol,l,NMRdata,ε,count)
		nop_ddf += count
		if λ == 1
			if l<n
				@debug "Partial solution by right side at level $(l)"  mol
				n_branch +=1
				quaternionBP_closure(l+1,pos+1,mol,Q)
			else
				nsol=nsol+1
				storage_mol[nsol] = copy(mol)
				@debug "Rank n was reached, a solution was found "
				return 0
			end
		else
			n_prune += 1
			# println("count prune = $(count_prun)")
		end
		if allmol==false && nsol>0
			@debug "number of solutions"  nsol
			@goto exit
		end
		Q = qprod(Q_before,Quaternion(a,-b,-c,d))
		nop_node += [12,16,0,0]
		qmol = rotopt(Q,D34)
		nop_node += [4,10,0,0]
		mol.atoms[l].x = qmol.v1 + mol.atoms[NMRdata.virtual_path[lastpos-1]].x
		mol.atoms[l].y = qmol.v2 + mol.atoms[NMRdata.virtual_path[lastpos-1]].y
		mol.atoms[l].z = qmol.v3 + mol.atoms[NMRdata.virtual_path[lastpos-1]].z
		nop_node += [3,0,0,0]
		@debug "candidate atom by left side at level $(l) = $(mol.atoms[l].x),$(mol.atoms[l].y),$(mol.atoms[l].z)"
		count = [0,0,0,0]
		ρ,count  = pruningtest(mol,l,NMRdata,ε,count)
		nop_ddf += count
		if ρ == 1
			if l<n
				@debug "Partial solution by left side at level $(l)" mol
				n_branch += 1
				quaternionBP_closure(l+1,pos+1,mol,Q)
			else
				nsol = nsol+1
				storage_mol[nsol] = copy(mol)
				@debug "Rank n was reached, a solution was found "
				return 0
			end
		else
			n_prune +=1
			# println("count prune = $(count_prun)")
		end

		@label exit
		return 0

	end #closure

	n = NMRdata.dim
	mol = MoleculeType(Vector{AtomType}(undef,n),0.0)
	for i=1:n
		mol.atoms[i] = AtomType(0.0,0.0,0.0)
	end
	Q = Quaternion(0.0,0.0,0.0,0.0)
	nsol = Int64(0)
	storage_mol = Dict{Int64,MoleculeType}()
	n_prune = 0
	n_branch = 0
	#count_nop = [+-,*,/,√]
	nop_node = [0,0,0,0]
	nop_ddf = [0,0,0,0]
	nop_vpath = [0,0,0,0]
	quaternionBP_closure(1,1,mol,Q)
	#println(" *** number of main operations evaluated $(count_op)")
	return nsol, storage_mol,Counter(nop_node,nop_vpath,nop_ddf,n_branch,n_prune)

end #solver quaternionBP


text:
using FlowAtlas
FlowAtlas.run("../data/workspace.wsp";files="../data/*/*.fcs",hold=false)

ext: jl

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
mutable struct ImageImporter <: AbstractImporter
    isenabled::Bool
end

function isenabled(importer::ImageImporter)
    return importer.isenabled
end

function enable!(importer::ImageImporter)
    importer.isenabled = true
end

function disable!(importer::ImageImporter)
    importer.isenabled  = false
end


text:
@testset "consumer" begin
    kfksvr = get(ENV, "KAFKA_SERVER", nothing)
    tpar = get(ENV, "TOPIC_PARTITION", nothing)
    if isnothing(kfksvr) || isnothing(tpar)
        @warn "Skip test consumer, KAFKA_SERVER=$kfksvr, TOPIC_PARTITION=$tpar"
        return
    end

    topic, par = split(tpar, ':')
    topic = string(topic)
    par = parse(Int, par)
    @info "connect $kfksvr ..."
    c = KafkaConsumer(kfksvr, "rdkafka-test")
    parlist = [(topic, par)]
    @info "subscribe $topic:$par ..."
    #subscribe(c, parlist)
    assign(c, parlist)
    seek(c, (topic, par), 0, 10000)
    pollthis() = poll(String, String, c)
    pollfirst(n) = begin
        for i = 1:10
            msg = pollthis()
            isnothing(msg) || return msg
        end
    end

    n = 10
    msg = pollfirst(n)
    if isnothing(msg)
        @warn "No result with fetch $n times"
        return
    end
    seek(c, (topic, par), 0)
    @test pollfirst(n) == msg
end


text:
import Base.push!

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

"""
    monomial_compress(io_equation, ode)

Compresses an input-output equation for the rank computation
Input:
- `io_equation` - input-output equation
- `ode` - the corresponding ODE model

Output:
- pair (coeffs, terms) such that:
    - sum of `coeffs[i] * terms[i] = io_equation`
    - `coeffs` involve only parameters, `terms` involve only inputs and outputs
    - length of the representation is the smallest possible
"""
function monomial_compress(io_equation, ode::ODE)
    return monomial_compress(io_equation, ode.parameters)
end

function monomial_compress(io_equation, params::Array{<: MPolyElem, 1})
    other_vars = [v for v in gens(parent(io_equation)) if !(var_to_str(v) in map(var_to_str, params))]
	coeffdict = extract_coefficients(io_equation, other_vars)
	expvect = collect(keys(coeffdict))
	coeffs = collect(values(coeffdict))
	termlist = map(x->prod(other_vars.^x), expvect)

    echelon_form = Array{Any, 1}()
    for (c, p) in zip(coeffs, termlist)
        for i in 1:length(echelon_form)
            basis_c = echelon_form[i][1]
            coef = coeff(c, leading_monomial(basis_c)) // leading_coefficient(basis_c)
            if coef != 0
                c = c - coef * basis_c
                echelon_form[i][2] += coef * p
            end
        end
        if c != 0
            push!(echelon_form, [c, p])
        end
    end

    result = ([a[1] for a in echelon_form], [a[2] for a in echelon_form])
    #s = 0
    #for (a, b) in zip(result[1], result[2])
    #    s += parent_ring_change(a, parent(io_equation)) * parent_ring_change(b, parent(io_equation))
    #end
    #println("====================")
    #println(s - io_equation)

    return result
end

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

# A data structure for storing exponent vectors as a trie in order
# to reduce the number of multiplications when computing many monomials
# in `massive_eval`
mutable struct ExpVectTrie
    depth::Int
    subtries::Dict{Int, ExpVectTrie}

    function ExpVectTrie(depth::Int)
        return new(depth, Dict{Int, ExpVectTrie}())
    end
end

function Base.push!(t::ExpVectTrie, vect::Array{Int, 1})
    if t.depth == 0
        if length(vect) != 0
            throw(DomainError("Inserting too long vector"))
        end
        return
    end

    f = last(vect)
    if !(f in keys(t.subtries))
        t.subtries[f] = ExpVectTrie(t.depth - 1)
    end
    push!(t.subtries[f], vect[1:end - 1])
end

"""
    get_max_below(t, vect)

Input:
- `t` - a trie with exponent vectors
- `vect` - yet another exponent vector

Output:
- a pair `(d, v)` where `v` is a vector in the trie which is componenwise ≤ `vect` and the difference `d` is as small as possible
"""
function get_max_below(t::ExpVectTrie, vect::Array{Int, 1})
    if t.depth == 0
        return (0, [])
    end
    min_diff, best_below = nothing, nothing
    for k in keys(t.subtries)
        last_diff = last(vect) - k
        if last_diff >= 0
            cur_diff, cur_below = get_max_below(t.subtries[k], vect[1:end - 1])
            if !isnothing(cur_diff)
                if isnothing(min_diff) || min_diff > cur_diff + last_diff
                    min_diff = cur_diff + last_diff
                    push!(cur_below, k)
                    best_below = cur_below
                end
            end
        end
    end
    return (min_diff, best_below)
end

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

"""
    massive_eval(polys, eval_dict)

Input:
- `polys` - a list of polynomials
- `eval_dict` - dictionary from variables to the values. Missing values are treated as zeroes

Output:
- a list of values of the polynomials

Evaluates a list of polynomails at a point. Assumes that multiplications are relatively expensive
(like in truncated power series) so all the monomials are precomputed first and the values of monomials
of lower degree are cached and used to compute the values of the monomials of higher degree
"""
function massive_eval(polys, eval_dict)
    R = parent(first(values(eval_dict)))
    point = [get(eval_dict, v, zero(R)) for v in gens(parent(first(polys)))]
    n = length(point)

    monomials = Set()
    for p in polys
        for exp in exponent_vectors(p)
            push!(monomials, exp)
        end
    end

    cache = Dict()
    cache[ [0 for i in 1:n] ] = one(R)
    cached_monoms = ExpVectTrie(n)
    push!(cached_monoms, [0 for _ in 1:n])

    for i in 1:n
        var_exp = [(i != j) ? 0 : 1 for j in 1:n]
        cache[var_exp] = point[i]
        push!(cached_monoms, var_exp)
    end

    for exp in sort!(collect(monomials), by=sum)
        if !(exp in keys(cache))
            monom_val = one(R)
            computed = [0 for i in 1:n]
            while sum(exp) > 0
                _, below = get_max_below(cached_monoms, exp)
                monom_val = monom_val * cache[below]
                exp = exp .- below
                computed = computed .+ below
                cache[computed] = monom_val
                push!(cached_monoms, computed)
            end
        end
    end

    results = []
    for p in polys
        res = zero(R)
        for (exp, coef) in zip(exponent_vectors(p), coefficients(p))
            res += coef * cache[exp]
        end
        push!(results, res)
    end
    return results
end

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

"""
    wronskian(io_equations, ode)

Input:
- `io_equations` - a set of io-equations in the form of the `Dict` as returned by `find_ioequations`
- `ode` - the `ODE` object

Output:
- a list of wronskians evaluated at a point modulo prime

Computes the wronskians of io_equations
"""
function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where P <: MPolyElem
    @debug "Compressing monomials"
    termlists = [monomial_compress(ioeq, ode)[2] for ioeq in values(io_equations)]
    @debug "Matrix sizes $(map(length, termlists))"

    # estimating the required order of truncation
    ord = max(map(length, termlists)...) + length(ode.x_vars)

    # reducing everything modulo prime
    PRIME = 2^31 - 1
    F = Nemo.GF(PRIME)
    polyring_red, gens_red = Nemo.PolynomialRing(F, map(var_to_str, gens(parent(termlists[1][1]))))
    termlists = [map(p -> parent_ring_change(p, polyring_red), tlist) for tlist in termlists]
    ode_red = reduce_ode_mod_p(ode, PRIME)

    @debug "Computing power series solution up to order $ord"
	ps = power_series_solution(
        ode_red,
        Dict(p => rand(1:100) for p in ode_red.parameters),
        Dict(x => rand(1:100) for x in ode_red.x_vars),
        Dict(u => [rand(1:100) for i in 0:ord] for u in ode_red.u_vars),
        ord
    )
    @debug "Computing the derivatives of the solution"
    ps_ext = Dict{MPolyElem, Nemo.gfp_abs_series}()# Generic.AbsSeries}()
    for v in vcat(ode_red.y_vars, ode_red.u_vars)
        cur_ps = ps[v]
        for i in 0:length(ode_red.x_vars)
            ps_ext[str_to_var(var_to_str(v) * "_$i", polyring_red)] = cur_ps
            cur_ps = ps_diff(cur_ps)
        end
    end

    @debug "Constructing wronskians"
    result = []
    for (i, tlist) in enumerate(termlists)
        n = length(tlist)
        evaled = massive_eval(tlist, ps_ext)
        S = Nemo.MatrixSpace(F, n, n)
        W = zero(S)
        for (i, ps) in enumerate(evaled)
            for j in 1:n
                W[i, j] = coeff(ps, j - 1)
            end
        end
        push!(result, W)
    end

    return result
end


text:
"""
$(TYPEDEF)

This is the supertype defining a section of the string into consideration which
can be considered as a "block". For instance, `abc { ... } def` contains one
braces block corresponding to the substring `{ ... }`. All subtypes of
AbstractBlock must have a `ss` field corresponding to the substring associated
to the block. See also [`Token`](@ref), [`OCBlock`](@ref).
"""
abstract type AbstractBlock end

# convenience functions to locate the substring associated to a block
from(β::AbstractBlock) = from(β.ss)
to(β::AbstractBlock)   = to(β.ss)
str(β::AbstractBlock)  = str(β.ss)

"""
$(TYPEDEF)

A token `τ::Token` denotes a part of the source string indicating a region that
may need further processing. It is identified by a symbol `τ.name` (e.g.:
`:MATH_ALIGN_OPEN`). Tokens are typically used in this code to identify
delimiters of environments. For instance, `abc \$ ... \$ def` contains two
tokens `:MATH_A` associated with the `\$` sign. Together they delimit here an
inline math expression.
"""
struct Token <: AbstractBlock
    name::Symbol
    ss::SubString
    lno::Int  # for LRINDENT it's useful to store line number
end
Token(n, s) = Token(n, s, 0)

to(β::Token) = ifelse(β.name == :EOS, from(β), to(β.ss))

context(t::Token) = context(str(t), from(t))

"""
$(TYPEDEF)

A special block to wrap special characters like like a html entity and have it
left unaffected by the Markdown to HTML transformation so that it can be
inserted "as is" in the HTML.
"""
struct HTML_SPCH <: AbstractBlock
    ss::SubString
    r::String
end
HTML_SPCH(ss) = HTML_SPCH(ss, "")

"""
$(TYPEDEF)

Prototype for an open-close block (see [`OCBlock`](@ref)) with the symbol of
the opening token (e.g. `:MATH_A`) and a corresponding list of closing tokens
(e.g. `(:MATH_A,)`). See also their definitions in `parser/md_tokens.jl`.
"""
struct OCProto
    name::Symbol
    otok::Symbol
    ctok::NTuple{N, Symbol} where N
    nest::Bool
end
OCProto(a...; nestable=false) = OCProto(a..., nestable)

"""
$(TYPEDEF)

Open-Close block, blocks that are defined by an opening `Token` and a closing
`Token`, they may be nested. For instance braces block are formed of an opening
`{` and a closing `}` and could be nested.
"""
struct OCBlock <: AbstractBlock
    name::Symbol
    ocpair::Pair{Token,Token}
    ss::SubString
end

"""
$(SIGNATURES)

Shorthand constructor to instantiate an `OCBlock` inferring the associated
substring from the `ocpair` (since it's the substring in between the tokens).
"""
OCBlock(η::Symbol, ω::Pair{Token,Token}, nestable::Bool=false) =
    OCBlock(η, ω, subs(str(ω.first), from(ω.first), to(ω.second)))

"""
$(SIGNATURES)

Convenience function to retrieve the opening token of an `OCBlock`.
"""
otok(ocb::OCBlock)::Token = ocb.ocpair.first

"""
$(SIGNATURES)

Convenience function to retrieve the closing token of an `OCBlock`.
"""
ctok(ocb::OCBlock)::Token = ocb.ocpair.second

"""
$(SIGNATURES)

Convenience function to return the content of an open-close block (`OCBlock`),
for instance the content of a `{...}` block would be `...`.
"""
function content(ocb::OCBlock)::SubString
    s = str(ocb.ss) # this does not allocate
    cfrom = nextind(s, to(otok(ocb)))
    c = ctok(ocb)
    if c.name == :EOS
        cto = from(c)
    else
        cto = prevind(s, from(c))
    end
    return subs(s, cfrom, cto)
end

"""
$(SIGNATURES)

Convenience function to strip the content of a block.
"""
stent = strip ∘ content

#=
FUNCTIONS / CONSTANTS THAT HELP DEFINE TOKENS
=#

"""
EOS

Convenience symbol to mark the end of the string to parse (helps with corner
cases where a token ends a document without being followed by a space).
"""
const EOS = '\0'

"""
SPACE_CHARS

Convenience list of characters that would correspond to a `\\s` regex. see also
https://github.com/JuliaLang/julia/blob/master/base/strings/unicode.jl.
To this are added things that would definitely switch token:
* EOS (end of string)
* '}' (if in command)
"""
const SPACE_CHAR = (' ', '\n', '\t', '\v', '\f', '\r', '\u85', '\ua0',
                    EOS, '}')

"""
SPACER

Convenience list of characters corresponding to digits.
"""
const NUM_CHAR   = ('1', '2', '3', '4', '5', '6', '7', '8', '9', '0')

"""
$(SIGNATURES)

Forward lookup checking if a sequence of characters matches `refstring` and is
followed (or not followed if `isfollowed==false`) by a character out of a list
of characters (`follow`). It returns

1. a number of steps indicating the number of characters to check,
2. whether there is an offset or not (if it is required to check a following
character or not),
3. a function that can be applied on a sequence of character.

# Example

```julia-repl
julia> (s, b, f) = isexactly("aabbc", "d", isfollowed=false);
julia> f("aabbcd")
false
julia> f("aabbce")
true
julia> s
5
```
"""
function isexactly(refstring::AS, follow::NTuple{K,Char} where K = (),
                   isfollowed=true)::Tuple{Int,Bool,Function,Bool}
    # number of code units from the start character
    steps = prevind(refstring, lastindex(refstring))
    # no offset (don't check next character)
    isempty(follow) && return (steps, false, (s,_) -> (s == refstring), false)
    # will include next char for verification (--> offset of 1)
    steps  = lastindex(refstring)
    nochop = !isfollowed || EOS ∈ follow
    # verification function; we want either (false false or true true))
    # NOTE, chop because we will take the next char with it

    λ(s, at_eos) = at_eos ?
                    # if at_eos  then we don't get an extra char, compare as is
                    # if isfollowed, then check EOS in follow
                    s == refstring && !isfollowed || EOS ∈ follow :
                    # if not, get extra char, compare  with
                    chop(s) == refstring && !xor(isfollowed, s[end] ∈ follow)

    return (steps, true, λ, nochop)
end

"""
$(SIGNATURES)

Check whether `c` is a letter or is in a vector of character `ac`.
"""
α(c::Char, ac::NTuple{K,Char}=()) where K = isletter(c) || (c ∈ ac)

"""
$(SIGNATURES)

Check whether `c` is alpha numeric or in vector of character `ac`
"""
αη(c::Char, ac::NTuple{K,Char}=()) where K = α(c, tuple(ac..., NUM_CHAR...))

"""
$(SIGNATURES)

Syntactic sugar for the incremental look case for which `steps=0` as well as
the `offset`. This is a case where from a start character we lazily accept the
next sequence of characters stopping as soon as a character fails to verify
`λ(c)`. See also [`isexactly`](@ref).
"""
incrlook(λ::Function, validator=nothing) = (0, false, λ, validator)

"""
$(SIGNATURES)

In combination with `incrlook`, checks to see if we have something that looks
like a @@div describing the opening of a div block. Triggering char is a first
`@`.
"""
is_div_open(i::Int, c::Char) = (i == 1 && return c == '@'; return α(c, ('-','_', ',', NUM_CHAR...)))

"""
$(SIGNATURES)

In combination with `incrlook`, checks to see if we have something that looks
like a triple backtick followed by a valid combination of letter defining a
language. Triggering char is a first backtick.
"""
is_language() = incrlook(_is_language, _validate_language)

function _is_language(i::Int, c::Char)
    i < 3  && return c == '`'  # ` followed by `` forms the opening ```
    i == 3 && return α(c)      # must be a letter
    return α(c, ('-',))        # can be a letter or a hyphen, for instance ```objective-c
end

_validate_language(stack::AS) = !isnothing(match(r"^```[a-zA-Z]", stack))

"""
$(SIGNATURES)

See [`is_language`](@ref) but with 5 ticks.
"""
is_language2() = incrlook(_is_language2, _validate_language2)

function _is_language2(i::Int, c::Char)
    i < 5  && return c == '`'
    i == 5 && return α(c)
    return α(c, ('-',))
end

_validate_language2(stack::AS) = !isnothing(match(r"^`````[a-zA-Z]", stack))

"""
$(SIGNATURES)

In combination with `incrlook`, checks to see if we have something that looks
like a html entity. Note that there can be fake matches, so this will need to
be validated later on; if validated it will be treated as HTML; otherwise it
will be shown as markdown. Triggerin char is a `&`.
"""
is_html_entity(i::Int, c::Char) = αη(c, ('#',';'))

"""
$(SIGNATURES)

Check if it looks like an emoji indicator `:...` note that it does
not take the final `:` this is checked and added in `validate_emoji!`.
"""
is_emoji(i::Int, c::Char) = αη(c, ('+','_','-'))

"""
$(SIGNATURES)

Check if it looks like `\\[\\^[\\p{L}0-9]+\\]:`.
"""
function is_footnote(i::Int, c::Char)
    i == 1 && return c == '^'
    i == 2 && return αη(c)
    i > 2  && return αη(c, (']', ':'))
end

"""
$SIGNATURES

Check if it looks like `---+`.
"""
is_hr1(::Int, c::Char) = c == '-'

"""
$SIGNATURES

Check if it looks like `___+`.
"""
is_hr2(::Int, c::Char) = c == '_'

"""
$SIGNATURES

Check if it looks like `***+`.
"""
is_hr3(::Int, c::Char) = c == '*'

"""
TokenFinder

Convenience type to define tokens. The Tuple comes from the output of functions
such as [`isexactly`](@ref).
"""
const TokenFinder =
        Pair{
            Tuple{
                Int,
                Bool,
                Function,
                Union{Bool,Nothing,Function}},
            Symbol}

"""
$(SIGNATURES)

Go through a text left to right, one (valid) char at the time and keep track of
sequences of chars that match specific tokens. The list of tokens found is
returned.

**Arguments**

* `str`:   the initial text
* `dictn`: dictionary of possible tokens (multiple characters long)
* `dict1`: dictionaro of possible tokens (single character)
"""
function find_tokens(str::AS,
                     dictn::AbstractDict{Char,Vector{TokenFinder}},
                     dict1::AbstractDict{Char,Symbol})::Vector{Token}
    isempty(str) && return Token[]
    # storage to keep track of the tokens found
    tokens    = Vector{Token}()
    head_idx  = 1                # valid string index
    end_idx   = lastindex(str)   # outer bound
    while head_idx <= end_idx
        head_char = str[head_idx]
        # 1. is it one of the single-char token?
        if haskey(dict1, head_char)
            tok = Token(dict1[head_char], subs(str, head_idx))
            push!(tokens, tok)

        # 2. is it one of the multi-char token?
    elseif haskey(dictn, head_char)
            for ((steps, offset, λ, ν), case) ∈ dictn[head_char]
                #=
                ↪ steps = length of the lookahead, 0 if incremental (greedy)
                ↪ offset = if we need to check one character 'too much'
                (e.g. this is the case if we want to check that something
                is followed by a space)
                ↪ λ = the checker function
                    * for a fixed lookahead, it returns true if the segment
                    (head_idx → head_idx + steps) matches a condition
                    * for an incremental lookahead, it returns true if chars
                    given meet a condition, chars after head_idx are fed while
                    the condition holds.
                ↪ ν = the (optional) validator function in the case of a greedy
                    lookahead to check whether the sequence is valid

                Either way, we push to the 'memory' the exact span (without the
                potential offset) of the token and a symbol indicating what it
                is then we move the head at the end of the token (note that
                it's pushed by 1 again after the if-else-end to start again).
                =#
                if steps > 0 # exact match of a given fixed pattern
                    tail_idx = nextind(str, head_idx, steps)
                    # is there space for the fixed pattern?
                    at_eos = false
                    if ν && tail_idx == nextind(str, end_idx)
                        tail_idx = end_idx
                        at_eos = true
                    end
                    tail_idx > end_idx && continue
                    # consider the substring and verify whether it matches
                    cand_seq = subs(str, head_idx, tail_idx)
                    if λ(cand_seq, at_eos)
                        # if offset==True --> looked at 1 extra char (lookahead)
                        back_one = offset & !at_eos
                        head_idx = prevind(str, tail_idx, back_one)
                        tok = Token(case, chop(cand_seq, tail=back_one))
                        push!(tokens, tok)
                        # token identified, no need to check other cases.
                        break
                    end
                else # rule-based match: greedy catch until fail
                    nchars     = 1
                    tail_idx   = head_idx
                    probe_idx  = nextind(str, head_idx)
                    probe_idx > end_idx && continue
                    probe_char = str[probe_idx]
                    while λ(nchars, probe_char)
                        tail_idx   = probe_idx
                        probe_idx  = nextind(str, probe_idx)
                        probe_idx > end_idx && break
                        probe_char = str[probe_idx]
                        nchars   += 1
                    end
                    if tail_idx > head_idx
                        cand_seq = subs(str, head_idx, tail_idx)
                        # check if the validator is happy otherwise skip
                        isnothing(ν) || ν(cand_seq) || continue
                        # if it's happy push the token & move after the match
                        tok = Token(case, cand_seq)
                        push!(tokens, tok)
                        head_idx = tail_idx
                    end
                end
            end
        end
        head_idx  = nextind(str, head_idx)
    end
    # finally push the end token
    eos = Token(:EOS, subs(str, end_idx))
    push!(tokens, eos)
    return tokens
end


text:
"""
Unit tests for Node type.

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

# Standard library
using Test

# External packages
using ZMQ

# Boid.jl package
using Boid

# --- Test fixtures

# Concrete subtypes of abstract types
include("fixtures/TestChannelData.jl")
include("fixtures/TestControlLogicCore.jl")
include("fixtures/TestProcessingCore.jl")

# Constants
control_url = construct_ipc_url(".", "control-unit-test.zmq")
control_ipc_path = control_url[7:end]

output_url = construct_ipc_url(".", "output-unit-test.zmq")
output_ipc_path = output_url[7:end]

input_urls = Vector{String}()
input_ipc_paths = Vector{String}()
for i in 1:3
    local input_url = construct_ipc_url(".", "output-unit-test-$i.zmq")
    push!(input_urls, input_url)
    push!(input_ipc_paths, input_url[7:end])
end

# --- Set up/tear down methods

function _tearDown()
    # Remove control IPC
    rm(control_ipc_path, force=true)

    # Remove output IPC
    rm(output_ipc_path, force=true)

    # Remove input IPCs
    for ipc_path in input_ipc_paths
        rm(ipc_path, force=true)
    end
end

# --- Tests

@testset "Node: constructor" begin
    # --- Preparations

    id = "node-1"

    processing_core = TestProcessingCore()

    control_unit_params = Dict(
        "logic_core"=>TestControlLogicCore(),
        "url"=>control_url,
        "copy_logic_core"=>false)

    output_channel_params = Dict(
        "url"=>output_url,
        "data_type"=>TestChannelData)

    input_channel_params = Vector{Dict}()
    local input_url
    for input_url in input_urls
        push!(input_channel_params,
              Dict(
                "url"=>input_url,
                "data_type"=>TestChannelData))
    end

    # --- Tests

    # ------- Normal usage

    node = Node(id, processing_core,
                control_unit_params,
                output_channel_params,
                input_channel_params)

    @test node.control_unit isa ControlUnit
    @test node.output_channel isa OutputChannel
    @test node.input_channels isa Vector{InputChannel}
    @test node.input_data isa Vector
    @test length(node.input_data) == length(input_channel_params)
    @test !is_running(node)

    # TODO: add tests

    # ------ Error cases

    # Invalid `control_unit_params`
    required_params = ["logic_core", "url"]
    for param in required_params
        invalid_control_unit_params = copy(control_unit_params)
        delete!(invalid_control_unit_params, param)
        expected_message =
            "`control_unit_params` is missing required parameter: `$param`"
        @test_throws(ArgumentError(expected_message),
                     Node(id, processing_core,
                          invalid_control_unit_params,
                          output_channel_params,
                          input_channel_params))
    end

    # Invalid `output_channel_params`
    required_params = ["data_type", "url"]
    for param in required_params
        invalid_output_channel_params = copy(output_channel_params)
        delete!(invalid_output_channel_params, param)
        expected_message =
            "`output_channel_params` is missing required parameter: `$param`"
        @test_throws(ArgumentError(expected_message),
                     Node(id, processing_core,
                          control_unit_params,
                          invalid_output_channel_params,
                          input_channel_params))
    end

    # Invalid `input_channel_params`
    required_params = ["data_type", "url"]
    for i in 1:length(input_channel_params)
        for param in required_params
            invalid_input_channel_params = deepcopy(input_channel_params)
            delete!(invalid_input_channel_params[i], param)
            expected_message =
            "`input_channel_params[$i]` is missing required parameter: `$param`"
            @test_throws(ArgumentError(expected_message),
                         Node(id, processing_core,
                              control_unit_params,
                              output_channel_params,
                              invalid_input_channel_params))
        end
    end

    # --- Clean up

    _tearDown()
end

@testset "Node: run()" begin

    # --- Preparations

    id = "node-1"

    processing_core = TestProcessingCore()

    control_unit_params = Dict(
        "logic_core"=>TestControlLogicCore(),
        "url"=>control_url,
        "copy_logic_core"=>false)

    output_channel_params = Dict(
        "url"=>output_url,
        "data_type"=>TestChannelData)

    input_channel_params = Vector{Dict}()
    local input_url
    for input_url in input_urls
        push!(input_channel_params,
              Dict(
                "url"=>input_url,
                "data_type"=>TestChannelData))
    end

    # Construct Node
    node = Node(id, processing_core,
                control_unit_params,
                output_channel_params,
                input_channel_params)

    # Construct control socket
    control_socket = Socket(REQ)
    connect(control_socket, control_url)

    # --- Tests
#=
    # Start node
    run_node_task = @task run(node)
    schedule(run_node_task)

#=
    # Wait for node to finish initializing
    while !is_running(node)
        sleep(0.1)
    end

    # Check that node is in 'running' state
    @test is_running(node)

    # Check that input channels are all in the 'listening' state
    for channel in node.input_channels
        while !is_listening(channel)
            sleep(0.1)
        end
        @test channel.state.is_listening
    end

    # Send STOP signal
    send(control_socket, STOP)
    response = recv(control_socket, String)
    @test response == "SUCCESS"

    # Wait for node to finish running
    wait(run_node_task)
    @test !is_running(node)

=#
    # --- Clean up
=#
    _tearDown()
end


text:
function lag(X::Matrix{S}, lags::Int; pad::Bool = false, T_by_n::Bool = true,
             drop_obs::Int = 0) where {S<:Real}
    if pad
        keep_lags = max(lags - drop_obs, 0)
        start = 1 + max(drop_obs - lags, 0)
        if T_by_n
            out = fill!(Matrix{S}(undef, size(X, 1) - drop_obs, size(X, 2)), NaN)
            out[1 + keep_lags:end, :] = X[start:size(X, 1) - lags, :]
        else
            out = fill!(Matrix{S}(undef, size(X, 1), size(X, 2) - drop_obs), NaN)
            out[:, 1 + keep_lags:end] = X[:, start:size(X, 2) - lags]
        end
    else
        out = T_by_n ? X[1 + drop_obs:size(X, 1) - lags, :] :
            X[:, 1 + drop_obs:size(X, 2) - lags]
    end
    return out
end

function lag_data(data::Matrix{S}, lags::Int; use_intercept::Bool = true,
                  pad::Bool = false, padding::Matrix{S} = Matrix{S}(undef, 0, 0)) where {S<:Real}
    # Assumes data is nobs x T
    nobs, T = size(data)
    data = Matrix(data')

    # Construct XX matrix of covariates
    add_constant = use_intercept ? 1 : 0
    XX = fill!(Matrix{S}(undef, pad ? T : T - lags, lags * nobs + add_constant), NaN)
    if use_intercept
        XX[:, 1] .= one(S) # XX is T x n_regressors
    end

    for i = 1:lags
        XX[:, add_constant + (i - 1) * nobs + 1:add_constant + i * nobs] =
            lag(data, i; pad = pad, T_by_n = true, drop_obs = pad ? 0 : lags - i)
    end

    if pad && !isempty(padding)
        XX[1:lags, :] = padding
    end

    return XX
end

function compute_var_population_moments(data::Matrix{S}, lags::Int;
                                 use_intercept::Bool = false) where {S<:Real}
    # Compute population moments of sample data
    YY = convert(Matrix{S}, data[:, 1 + lags:end]')
    XX = lag_data(data, lags; use_intercept = use_intercept) # Construct XX matrix of covariates
    YYYY = YY' * YY
    XXYY = XX' * YY
    XXXX = XX' * XX

    return YYYY, XXYY, XXXX
end

"""
```
draw_stationary_VAR(YYYYC, XXYYC, XXXXC, T̄, n_obs, lags; standard_orientation = true)
draw_stationary_VAR(YYYYC, XXYYC, XXXXC, T̄; standard_orientation = true)
```
draws β and Σ from the distribution p(β, Σ | Y, θ) implied by the
population moments (or covariances) YYYYC, XXYYC, and XXXXC for a
VAR with parameters θ.

For example, if these population moments are generated by
a DSGE-VAR, then θ are the structural parameters of the DSGE
and the weight λ placed on the cross-restrictions implied
by the DSGE. The population moments would represent the
moments of the sample data and dummy observables generated
to implement the DSGE prior.

Given these moments, we compute the maximum-likelihood
estimates of β and Σ using OLS. Denote these estimates
by Β and S. Then we generate draws from p(β, Σ | Y, θ)
using the fact that
```
Σ | Y, θ ∼ ℐ𝒲 (T̄ × S, T̄ - (1 + lags * n_obs), n_obs),
β | Y, Σ,θ ∼ 𝒩 (B, Σ ⊗ XXXXC⁻¹).
```
Finally, we check that these draws generate a stationary
state space system. If they do not, then we keep drawing
until we obtain a pair of draws (β, Σ) that are stationary.

### Inputs
* `YYYYC::Matrix{<:Real}`: covariance of observables
* `XXYYC::Matrix{<:Real}`: covariance of observables with their lags
* `XXXXC::Matrix{<:Real}`: covariance of the lags of the observables
* `T̄::Int`: total number of time periods of observations, including sample observables
    from actual data and any dummy observables generated to implement priors.
* `n_obs::Int`: number of distinct observables
* `lags::Int`: number of lags in the VAR

### Keywords
* `standard_orientation::Bool`: if true, the draw of `β` has
    dimensions `(n_obs * lags) x n_obs`. Otherwise, it has the transposed dimensions.
* All other keywords are used for testing purposes.
"""
function draw_stationary_VAR(YYYYC::Matrix{S}, XXYYC::Matrix{S}, XXXXC::Matrix{S},
                             T̄::Int, n_obs::Int, lags::Int; standard_orientation::Bool = true,
                             testing::Bool = false,
                             test_Σ_draw_shock::Matrix{S} = Matrix{S}(undef, 0, 0),
                             test_β_draw_shock::Vector{S} = Vector{S}(undef, 0)) where {S<:Real}

    # Set up
    k = 1 + lags * n_obs
    inv_XXXXC = inv(XXXXC)
    β = inv_XXXXC * XXYYC
    inv_Σ_mul_T̄ = inv(YYYYC - XXYYC' * β)
    inv_Σ_mul_T̄ += inv_Σ_mul_T̄' # force to be positive definite
    inv_Σ_mul_T̄ ./= 2.
    cholmat = cholesky(inv_Σ_mul_T̄).L
    β = vec(β)

    if testing # just do one draw each
        # Draw from marginal posterior of Σ (based on DSGE-VAR)
        z = cholmat * test_Σ_draw_shock
        Σ_draw = inv(z * z')

        # Draw from the conditional posterior of β (based on DSGE-VAR)
        vc      = kron(Σ_draw, inv_XXXXC)
        vc      += vc'
        vc      ./= 2.
        β_draw  = convert(Matrix{S}, reshape(β + cholesky(vc).L * test_β_draw_shock, k, n_obs)')
    else
        stationary = false
        while !stationary
            # Draw from marginal posterior of Σ (based on DSGE-VAR)
            z = cholmat * randn(n_obs, T̄ - k)
            Σ_draw = inv(z * z')

            # Draw from the conditional posterior of β (based on DSGE-VAR)
            vc       = kron(Σ_draw, inv_XXXXC)
            vc       += vc'
            vc       ./= 2.
            β_draw   = convert(Matrix{S}, reshape(β + cholesky(vc).L * randn(n_obs * k), k, n_obs)')
            β_to_TTT = vcat(β_draw[:, (1+1):k],
                            hcat(Matrix{S}(I, n_obs * (lags - 1), n_obs * (lags - 1)),
                                 zeros(S, n_obs * (lags - 1), n_obs)))
            if maximum(abs.(eigen(β_to_TTT).values)) < 1
                stationary = true
            end
        end
    end

    if standard_orientation
        β_draw = convert(Matrix{S}, β_draw')
    end

    return β_draw, Σ_draw
end

function draw_stationary_VAR(YYYYC::Matrix{S}, XXYYC::Matrix{S}, XXXXC::Matrix{S},
                             T̄::Int; standard_orientation::Bool = true, testing::Bool = false,
                             test_Σ_draw_shock::Matrix{S} = Matrix{S}(undef, 0, 0),
                             test_β_draw_shock::Vector{S} = Vector{S}(undef, 0)) where {S<:Real}
    # Infer n_obs and lags from population moments
    n_obs = size(YYYYC, 1)
    lags = (size(XXXXC, 1) - 1) / n_obs

    return draw_stationary_VAR(YYYYC, XXYYC, XXXXC, T̄, n_obs, lags;
                               standard_orientation = standard_orientation, testing = testing,
                               test_Σ_draw_shock = test_Σ_draw_shock, test_β_draw_shock = test_β_draw_shock)
end

"""
```
draw_VECM(YYYYC, XXYYC, XXXXC, T̄, n_obs, lags, n_coint; standard_orientation = true)
```
draws β and Σ from the distribution p(β, Σ | Y, θ) implied by the
population moments (or covariances) YYYYC, XXYYC, and XXXXC for a
VECM with parameters θ.

For example, if these population moments are generated by
a DSGE-VECM, then θ are the structural parameters of the DSGE
and the weight λ placed on the cross-restrictions implied
by the DSGE. The population moments would represent the
moments of the sample data and dummy observables generated
to implement the DSGE prior.

Given these moments, we compute the maximum-likelihood
estimates of β and Σ using OLS. Denote these estimates
by Β and S. Then we generate draws from p(β, Σ | Y, θ)
using the fact that
```
Σ | Y, θ ∼ ℐ𝒲 (T̄ × S, T̄ - (1 + lags * n_obs), n_obs),
β | Y, Σ,θ ∼ 𝒩 (B, Σ ⊗ XXXXC⁻¹).
```

### Inputs
* `YYYYC::Matrix{<:Real}`: covariance of observables
* `XXYYC::Matrix{<:Real}`: covariance of observables with their lags
* `XXXXC::Matrix{<:Real}`: covariance of the lags of the observables
* `T̄::Int`: total number of time periods of observations, including sample observables
    from actual data and any dummy observables generated to implement priors.
* `n_obs::Int`: number of distinct observables
* `lags::Int`: number of lags in the VECM
* `n_coint::Int`: number of distinct cointegrating relationships

### Keywords
* `standard_orientation::Bool`: if true, the draw of `β` has
    dimensions `(n_obs * lags) x n_obs`. Otherwise, it has the transposed dimensions.
* All other keywords are used for testing purposes.
"""
function draw_VECM(YYYYC::Matrix{S}, XXYYC::Matrix{S}, XXXXC::Matrix{S},
                   T̄::Int, n_obs::Int, lags::Int, n_coint::Int; standard_orientation::Bool = true,
                   testing::Bool = false,
                   test_Σ_draw_shock::Matrix{S} = Matrix{S}(undef, 0, 0),
                   test_β_draw_shock::Vector{S} = Vector{S}(undef, 0)) where {S<:Real}

    # Set up
    k = 1 + lags * n_obs + n_coint
    inv_XXXXC = inv(XXXXC)
    β = inv_XXXXC * XXYYC
    inv_Σ_mul_T̄ = inv(YYYYC - XXYYC' * β)
    inv_Σ_mul_T̄ += inv_Σ_mul_T̄' # force to be positive definite
    inv_Σ_mul_T̄ ./= 2.
    cholmat = cholesky(inv_Σ_mul_T̄).L
    β = vec(β)

    if testing # just do one draw each
        # Draw from marginal posterior of Σ (based on DSGE-VECM)
        z = cholmat * test_Σ_draw_shock
        Σ_draw = inv(z * z')

        # Draw from the conditional posterior of β (based on DSGE-VECM)
        vc     = kron(Σ_draw, inv_XXXXC)
        vc    += vc'
        vc   ./= 2.
        β_draw = reshape(β + cholesky(vc).L * test_β_draw_shock, k, n_obs)
    else
        # Draw from marginal posterior of Σ (based on DSGE-VECM)
        z = cholmat * randn(n_obs, T̄ - k)
        Σ_draw = inv(z * z')

        # Draw from the conditional posterior of β (based on DSGE-VECM)
        vc     = kron(Σ_draw, inv_XXXXC)
        vc    += vc'
        vc    /= 2.
        β_draw = reshape(β + cholesky(vc).L * randn(n_obs * k), k, n_obs)
    end

    if !standard_orientation
        β_draw = convert(Matrix{S}, β_draw')
    end

    return β_draw, Σ_draw
end


text:
using Statistics
## Shootouts

mutable struct Shootout
  setups::Vector{Dict{Symbol,Any}}
  times#::Vector{Float64}
  errors#::Vector{uType}
  effs#::Vector{Float64} # Efficiencies
  effratios#::Matrix{uEltype}
  solutions
  names::Vector{String}
  N::Int
  bestidx::Int
  winner::String
end

mutable struct ShootoutSet
  shootouts::Vector{Shootout}
  probs#::Vector{DEProblem}
  probaux#::Vector{Dict{Symbol,Any}}
  N::Int
  winners::Vector{String}
end

function ode_shootout(args...;kwargs...)
  @warn("ode_shootout is deprecated. Use ShootOut instead")
  ShootOut(args...;kwargs...)
end

function Shootout(prob,setups;appxsol=nothing,names=nothing,error_estimate=:final,numruns=20,seconds=2,kwargs...)
  N = length(setups)
  @assert names === nothing || length(setups) == length(names)
  errors = Vector{Float64}(undef,N)
  solutions = Vector{Any}(undef,N)
  effs = Vector{Float64}(undef,N)
  times = Vector{Float64}(undef,N)
  effratios = Matrix{Float64}(undef,N,N)
  timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
  dense_errors = error_estimate ∈ DENSE_ERRORS
  if names === nothing
    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
  end
  for i in eachindex(setups)
    sol = solve(prob,setups[i][:alg];timeseries_errors=timeseries_errors,
    dense_errors = dense_errors,kwargs...,setups[i]...)

    if :prob_choice ∈ keys(setups[i])
      cur_appxsol = appxsol[setups[i][:prob_choice]]
    elseif prob isa AbstractArray
      cur_appxsol = appxsol[1]
    else
      cur_appxsol = appxsol
    end

    if cur_appxsol !== nothing
      errsol = appxtrue(sol,cur_appxsol)
      errors[i] = errsol.errors[error_estimate]
      solutions[i] = errsol
    else
      errors[i] = sol.errors[error_estimate]
      solutions[i] = sol
    end

    if haskey(setups[i], :prob_choice)
      _prob = prob[setups[i][:prob_choice]]
    elseif prob isa AbstractArray
      _prob = prob[1]
    else
      _prob = prob
    end

    benchmark_f = let _prob=_prob,alg=setups[i][:alg],sol=sol,kwargs=kwargs
      () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
                           timeseries_errors = false,
                           dense_errors = false, kwargs...)
    end
    benchmark_f() # pre-compile

    b_t = benchmark_f()
    if b_t > seconds
      times[i] = b_t
    else
      times[i] = mapreduce(i -> benchmark_f(), min, 2:numruns; init = b_t)
    end

    effs[i] = 1/(errors[i]*times[i])
  end
  for j in 1:N, i in 1:N
    effratios[i,j] = effs[i]/effs[j]
  end
  bestidx = findall((y)->y==maximum(effs),effs)[1]; winner = names[bestidx]
  return Shootout(setups,times,errors,effs,effratios,solutions,names,N,bestidx,winner)
end

function ode_shootoutset(args...;kwargs...)
  @warn("ode_shootoutset is deprecated. Use ShootoutSet instead")
  ShootoutSet(args...;kwargs...)
end

function ShootoutSet(probs,setups;probaux=nothing,
                     names=nothing,print_names=false,kwargs...)
  N = length(probs)
  @assert names === nothing || length(setups) == length(names)
  shootouts = Vector{Shootout}(undef,N)
  winners = Vector{String}(undef,N)
  if names === nothing
    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
  end
  if probaux === nothing
    probaux = Vector{Dict{Symbol,Any}}(undef,N)
    for i in 1:N
      probaux[i] = Dict{Symbol,Any}()
    end
  end
  for i in eachindex(probs)
    print_names && println(names[i])
    shootouts[i] = Shootout(probs[i],setups;names=names,kwargs...,probaux[i]...)
    winners[i] = shootouts[i].winner
  end
  return ShootoutSet(shootouts,probs,probaux,N,winners)
end

Base.length(shoot::Shootout) = shoot.N
Base.size(shoot::Shootout) = length(shoot)
Base.getindex(shoot::Shootout,i::Int) = shoot.effs[i]
Base.getindex(shoot::Shootout,::Colon) = shoot.effs
Base.firstindex(shoot::Shootout) = 1
Base.lastindex(shoot::Shootout) = lastindex(shoot.effs)

function Base.show(io::IO, shoot::Shootout)
  println(io,"Winner: $(shoot.winner)")
  println(io,"EffRatios: $(shoot.effratios[shoot.bestidx,:])")
end

Base.length(set::ShootoutSet) = set.N
Base.size(set::ShootoutSet) = length(set)
Base.getindex(set::ShootoutSet,i::Int) = set.shootouts[i]
Base.getindex(set::ShootoutSet,::Colon) = set.shootouts
Base.show(io::IO, set::ShootoutSet) = print(io,"ShootoutSet of $(set.N) shootouts ")
Base.firstindex(shoot::ShootoutSet) = 1
Base.lastindex(shoot::ShootoutSet) = lastindex(shoot.shootouts)

## WorkPrecisions

mutable struct WorkPrecision
  prob
  abstols
  reltols
  errors
  times
  name
  N::Int
end

mutable struct WorkPrecisionSet
  wps::Vector{WorkPrecision}
  N::Int
  abstols
  reltols
  prob
  setups
  names
  error_estimate
  numruns
end

function WorkPrecision(prob,alg,abstols,reltols,dts=nothing;
                       name=nothing,appxsol=nothing,error_estimate=:final,numruns=20,seconds=2,kwargs...)
  N = length(abstols)
  errors = Vector{Float64}(undef,N)
  times = Vector{Float64}(undef,N)
  if name === nothing
    name = "WP-Alg"
  end

  if haskey(kwargs, :prob_choice)
    _prob = prob[kwargs[:prob_choice]]
  elseif prob isa AbstractArray
    _prob = prob[1]
  else
    _prob = prob
  end

  let _prob = _prob
    timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
    dense_errors = error_estimate ∈ DENSE_ERRORS
    for i in 1:N
      if dts === nothing
        sol = solve(_prob,alg;kwargs...,abstol=abstols[i],
        reltol=reltols[i],timeseries_errors=timeseries_errors,
        dense_errors = dense_errors)
      else
        sol = solve(_prob,alg;kwargs...,abstol=abstols[i],
        reltol=reltols[i],dt=dts[i],timeseries_errors=timeseries_errors,
        dense_errors = dense_errors)
      end

      if haskey(kwargs, :prob_choice)
        cur_appxsol = appxsol[kwargs[:prob_choice]]
      elseif prob isa AbstractArray
        cur_appxsol = appxsol[1]
      else
        cur_appxsol = appxsol
      end

      if cur_appxsol !== nothing
        errsol = appxtrue(sol,cur_appxsol)
        errors[i] = mean(errsol.errors[error_estimate])
      else
        errors[i] = mean(sol.errors[error_estimate])
      end

      benchmark_f = let dts=dts,_prob=_prob,alg=alg,sol=sol,abstols=abstols,reltols=reltols,kwargs=kwargs
        if dts === nothing
          if typeof(_prob) <: DAEProblem
            () -> @elapsed solve(_prob, alg, sol.u, sol.t;
                                abstol = abstols[i],
                                reltol = reltols[i],
                                timeseries_errors = false,
                                dense_errors = false, kwargs...)
          else
            () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
                                abstol = abstols[i],
                                reltol = reltols[i],
                                timeseries_errors = false,
                                dense_errors = false, kwargs...)
          end
        else
          if typeof(_prob) <: DAEProblem
            () -> @elapsed solve(_prob, alg, sol.u, sol.t;
                                abstol = abstols[i],
                                reltol = reltols[i],
                                dt = dts[i],
                                timeseries_errors = false,
                                dense_errors = false, kwargs...)
          else
            () -> @elapsed solve(_prob, alg, sol.u, sol.t, sol.k;
                                abstol = abstols[i],
                                reltol = reltols[i],
                                dt = dts[i],
                                timeseries_errors = false,
                                dense_errors = false, kwargs...)
          end
        end
      end
      benchmark_f() # pre-compile

      b_t = benchmark_f()
      if b_t > seconds
        times[i] = b_t
      else
        times[i] = mapreduce(i -> benchmark_f(), min, 2:numruns; init = b_t)
      end
    end
  end
  return WorkPrecision(prob,abstols,reltols,errors,times,name,N)
end

function WorkPrecisionSet(prob,
                          abstols,reltols,setups;
                          print_names=false,names=nothing,appxsol=nothing,
                          error_estimate=:final,
                          test_dt=nothing,kwargs...)
  N = length(setups)
  @assert names === nothing || length(setups) == length(names)
  wps = Vector{WorkPrecision}(undef,N)
  if names === nothing
    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
  end
  for i in 1:N
    print_names && println(names[i])
    _abstols = get(setups[i],:abstols,abstols)
    _reltols = get(setups[i],:reltols,reltols)
        _dts = get(setups[i],:dts,nothing)
    filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[i])

    wps[i] = WorkPrecision(prob,setups[i][:alg],_abstols,_reltols,_dts;
                               appxsol=appxsol,
                               error_estimate=error_estimate,
                               name=names[i],kwargs...,filtered_setup...)
  end
  return WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,nothing)
end

@def error_calculation begin
  if !DiffEqBase.has_analytic(prob.f)
    t = prob.tspan[1]:test_dt:prob.tspan[2]
    brownian_values = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
    brownian_values2 = cumsum([[zeros(size(prob.u0))];[sqrt(test_dt)*randn(size(prob.u0)) for i in 1:length(t)-1]])
    np = NoiseGrid(t,brownian_values,brownian_values2)
    _prob = remake(prob,noise=np);
    true_sol = solve(_prob,appxsol_setup[:alg];kwargs...,appxsol_setup...)
  else
    _prob = prob
  end

  # Get a cache
  _abstols = get(setups[1],:abstols,abstols)
  _reltols = get(setups[1],:reltols,reltols)
  _dts = get(setups[1],:dts,zeros(length(_abstols)))
  filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[1])

  sol = solve(_prob,setups[1][:alg];
        kwargs...,filtered_setup...,abstol=_abstols[1],
        reltol=_reltols[1],dt=_dts[1],
        timeseries_errors=false,
        dense_errors = false)

  for j in 1:M, k in 1:N
    _abstols = get(setups[k],:abstols,abstols)
    _reltols = get(setups[k],:reltols,reltols)
    _dts = get(setups[k],:dts,zeros(length(_abstols)))
    filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[k])

    sol = solve(_prob,setups[k][:alg];
          kwargs...,filtered_setup...,abstol=_abstols[j],
          reltol=_reltols[j],dt=_dts[j],
          timeseries_errors=timeseries_errors,
          dense_errors = dense_errors)
    DiffEqBase.has_analytic(prob.f) ? err_sol = sol : err_sol = appxtrue(sol,true_sol)
    tmp_solutions[i,j,k] = err_sol
  end
end

function WorkPrecisionSet(prob::AbstractRODEProblem,abstols,reltols,setups,test_dt=nothing;
                          numruns=20,numruns_error = 20,
                          print_names=false,names=nothing,appxsol_setup=nothing,
                          error_estimate=:final,parallel_type = :none,
                          kwargs...)

  @assert names === nothing || length(setups) == length(names)
  timeseries_errors = DiffEqBase.has_analytic(prob.f) && error_estimate ∈ TIMESERIES_ERRORS
  weak_timeseries_errors = error_estimate ∈ WEAK_TIMESERIES_ERRORS
  weak_dense_errors = error_estimate ∈ WEAK_DENSE_ERRORS
  dense_errors = DiffEqBase.has_analytic(prob.f) && error_estimate ∈ DENSE_ERRORS
  N = length(setups); M = length(abstols)
  times = Array{Float64}(undef,M,N)
  tmp_solutions = Array{Any}(undef,numruns_error,M,N)
  if names === nothing
    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
  end
  time_tmp = Vector{Float64}(undef,numruns)

  # First calculate all of the errors
  if parallel_type == :threads
    Threads.@threads for i in 1:numruns_error
      @error_calculation
    end
  elseif parallel_type == :none
    for i in 1:numruns_error
      @info "Error calculation: $i/$numruns_error"
      @error_calculation
    end
  end

  _solutions_k = [[EnsembleSolution(tmp_solutions[:,j,k],0.0,true) for j in 1:M] for k in 1:N]
  solutions = [[DiffEqBase.calculate_ensemble_errors(sim;weak_timeseries_errors=weak_timeseries_errors,weak_dense_errors=weak_dense_errors) for sim in sol_k] for sol_k in _solutions_k]
  if error_estimate ∈ WEAK_ERRORS
    errors = [[solutions[j][i].weak_errors[error_estimate] for i in 1:M] for j in 1:N]
  else
    errors = [[solutions[j][i].error_means[error_estimate] for i in 1:M] for j in 1:N]
  end

  local _sol

  # Now time it
  for k in 1:N
    # precompile
    GC.gc()
    _abstols = get(setups[k],:abstols,abstols)
    _reltols = get(setups[k],:reltols,reltols)
    _dts = get(setups[k],:dts,zeros(length(_abstols)))
    filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[k])

    _sol = solve(prob,setups[k][:alg];
          kwargs...,filtered_setup...,abstol=_abstols[1],
          reltol=_reltols[1],dt=_dts[1],
          timeseries_errors=false,
          dense_errors = false)
    x = isempty(_sol.t) ? 0 : round(Int,mean(_sol.t) - sum(_sol.t)/length(_sol.t))
    GC.gc()
    for j in 1:M
      for i in 1:numruns
        time_tmp[i] = @elapsed sol = solve(prob,setups[k][:alg];
                        kwargs...,filtered_setup...,abstol=_abstols[j],
                        reltol=_reltols[j],dt=_dts[j],
                        timeseries_errors=false,
                        dense_errors = false)
      end
      times[j,k] = mean(time_tmp) + x
      GC.gc()
    end
  end

  wps = [WorkPrecision(prob,abstols,reltols,errors[i],times[:,i],names[i],N) for i in 1:N]
  WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,numruns_error)
end

function WorkPrecisionSet(prob::AbstractEnsembleProblem,abstols,reltols,setups,test_dt=nothing;
                          numruns=5, trajectories=1000,
                          print_names=false,names=nothing,appxsol_setup=nothing,expected_value=nothing,
                          error_estimate=:weak_final,ensemblealg=EnsembleThreads(),
                          kwargs...)

  @assert names === nothing || length(setups) == length(names)

  weak_timeseries_errors = error_estimate ∈ WEAK_TIMESERIES_ERRORS
  weak_dense_errors = error_estimate ∈ WEAK_DENSE_ERRORS

  N = length(setups); M = length(abstols)
  times = Array{Float64}(undef,M,N)
  solutions = Array{Any}(undef,M,N)
  if names === nothing
    names = [string(nameof(typeof(setup[:alg]))) for setup in setups]
  end
  time_tmp = Vector{Float64}(undef,numruns)

  # First calculate all of the errors
  for k in 1:N
    _abstols = get(setups[k],:abstols,abstols)
    _reltols = get(setups[k],:reltols,reltols)
    _dts = get(setups[k],:dts,zeros(length(_abstols)))
    filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[k])

    for j in 1:M
      sol = solve(prob,setups[k][:alg],ensemblealg;
        filtered_setup...,
        abstol=_abstols[j],
        reltol=_reltols[j],
        dt=_dts[j],
        timeseries_errors=false,
        dense_errors = false,
        trajectories=Int(trajectories),kwargs...)
      solutions[j,k] = sol
    end
    @info "$(setups[k][:alg]) ($k/$N)"
  end

  if error_estimate ∈ WEAK_ERRORS
    if expected_value != nothing
      if error_estimate == :weak_final
        errors = [[LinearAlgebra.norm(Statistics.mean(solutions[i,j].u .- expected_value))
          for i in 1:M] for j in 1:N]
      elseif error_estimate == :weak_l2
        errors = [[LinearAlgebra.norm(Statistics.mean(solutions[i,j] .- expected_value))
          for i in 1:M] for j in 1:N]
      else
        error("Error estimate $error_estimate is not implemented yet.")
      end
    else
      sol = solve(prob,appxsol_setup[:alg],ensemblealg;kwargs...,appxsol_setup...,
          timeseries_errors=false,dense_errors = false,trajectories=Int(trajectories))
      errors = [[LinearAlgebra.norm(Statistics.mean(solutions[i,j].u .- sol.u))
        for i in 1:M] for j in 1:N]
    end
  else
    error("use RODEProblem instead of EnsembleProblem for strong errors.")
  end

  local _sol

  # Now time it
  for k in 1:N
    # precompile
    GC.gc()
    _abstols = get(setups[k],:abstols,abstols)
    _reltols = get(setups[k],:reltols,reltols)
    _dts = get(setups[k],:dts,zeros(length(_abstols)))
    filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[k])

    _sol = solve(prob,setups[k][:alg],ensemblealg;
      filtered_setup...,
      abstol=_abstols[1],
      reltol=_reltols[1],
      dt=_dts[1],
      timeseries_errors=false,
      dense_errors = false,
      trajectories=Int(trajectories),kwargs...)
    #x = isempty(_sol.t) ? 0 : round(Int,mean(_sol.t) - sum(_sol.t)/length(_sol.t))
    GC.gc()
    for j in 1:M
      for i in 1:numruns
        time_tmp[i] = @elapsed sol = solve(prob,setups[k][:alg],ensemblealg;
                      filtered_setup...,
                      abstol=_abstols[j],
                      reltol=_reltols[j],
                      dt=_dts[j],
                      timeseries_errors=false,
                      dense_errors = false,
                      trajectories=Int(trajectories),kwargs...)
      end
      times[j,k] = mean(time_tmp) #+ x
      GC.gc()
    end
  end

  wps = [WorkPrecision(prob,abstols,reltols,errors[i],times[:,i],names[i],N) for i in 1:N]
  WorkPrecisionSet(wps,N,abstols,reltols,prob,setups,names,error_estimate,Int(trajectories))
end

function get_sample_errors(prob::AbstractRODEProblem,setup,test_dt=nothing;
                          appxsol_setup=nothing,
                          numruns,error_estimate=:final,
                          sample_error_runs = Int(1e7),
                          solution_runs,
                          parallel_type = :none,kwargs...)

  maxnumruns = findmax(numruns)[1]

  tmp_solutions_full = map(1:solution_runs) do i
    @info "Solution Run: $i"
    # Use the WorkPrecision stuff to calculate the errors
    tmp_solutions = Array{Any}(undef,maxnumruns,1,1)
    setups = [setup]
    abstols = [1e-2] # Standard default
    reltols = [1e-2] # Standard default
    M = 1; N = 1
    timeseries_errors = false; dense_errors = false
    if parallel_type == :threads
      Threads.@threads for i in 1:maxnumruns
        @error_calculation
      end
    elseif parallel_type == :none
      for i in 1:maxnumruns
        @error_calculation
      end
    end
    tmp_solutions = vec(tmp_solutions)
  end

  if DiffEqBase.has_analytic(prob.f)
    analytical_mean_end = mean(1:sample_error_runs) do i
      _dt = prob.tspan[2] - prob.tspan[1]
      if typeof(prob.u0) <: Number
        W = sqrt(_dt)*randn()
      else
        W = sqrt(_dt)*randn(size(prob.u0))
      end
      prob.f.analytic(prob.u0,prob.p,prob.tspan[2],W)
    end
  else
    # Use the mean of the means as the analytical mean
    analytical_mean_end = mean(mean(tmp_solutions[i].u[end] for i in 1:length(tmp_solutions)) for tmp_solutions in tmp_solutions_full)
  end

  if numruns isa Number
    mean_solution_ends = [mean([tmp_solutions[i].u[end] for i in 1:maxnumruns]) for tmp_solutions in tmp_solutions_full]
    return sample_error = 1.96std(norm(mean_sol_end - analytical_mean_end) for mean_sol_end in mean_solution_ends)/sqrt(numruns)
  else
    map(1:length(numruns)) do i
      mean_solution_ends = [mean([tmp_solutions[i].u[end] for i in 1:numruns[i]]) for tmp_solutions in tmp_solutions_full]
      sample_error = 1.96std(norm(mean_sol_end - analytical_mean_end) for mean_sol_end in mean_solution_ends)/sqrt(numruns[i])
    end
  end
end

Base.length(wp::WorkPrecision) = wp.N
Base.size(wp::WorkPrecision) = length(wp)
Base.getindex(wp::WorkPrecision,i::Int) = wp.times[i]
Base.getindex(wp::WorkPrecision,::Colon) = wp.times
Base.firstindex(wp::WorkPrecision) = 1
Base.lastindex(wp::WorkPrecision) = lastindex(wp.times)

function Base.show(io::IO, wp::WorkPrecision)
  println(io,"Name: $(wp.name)")
  println(io,"Times: $(wp.times)")
  println(io,"Errors: $(wp.errors)")
end

Base.length(wp_set::WorkPrecisionSet) = wp_set.N
Base.size(wp_set::WorkPrecisionSet) = length(wp_set)
Base.getindex(wp_set::WorkPrecisionSet,i::Int) = wp_set.wps[i]
Base.getindex(wp_set::WorkPrecisionSet,::Colon) = wp_set.wps
Base.firstindex(wp_set::WorkPrecisionSet) = 1
Base.lastindex(wp_set::WorkPrecisionSet)  = lastindex(wp_set.wps)

Base.show(io::IO, wp_set::WorkPrecisionSet) = println(io,"WorkPrecisionSet of $(wp_set.N) wps")


text:
import GradDescent: Optimizer

struct ModelOptimizer{M, O <: Optimizer}
    model::Type{M}
    optimizers::Dict{Symbol, O}
end

function ModelOptimizer(::Type{M}, opt::Optimizer) where M
    optimizers = Dict{Symbol, Optimizer}()
    for fld in fieldnames(M)
        optimizers[fld] = deepcopy(opt)
    end
    return ModelOptimizer(M, optimizers)
end

Base.show(io::IO, opt::ModelOptimizer) = print(io, "ModelOptimizer($(opt.model))")

ModelOptimizer(m, opt::Optimizer) = ModelOptimizer(typeof(m), opt)

function update_params!(m_opt::ModelOptimizer, m, grad)
    @assert(typeof(m) == typeof(grad), "Wrong gradient type: model is $(typeof(m)), " *
            "but gradient is $(typeof(grad))")
    for fld in fieldnames(m)
        # TODO: this will break on scalar parameters
        theta = getfield(m, fld)
        theta .-= update(m_opt.optimizers[fld], getfield(grad, fld))
        setfield!(m, fld, theta)
    end
end


text:
# This file assumes that the directories are organized as follows:
#
# ./SomeDirectory/RRTMGP.jl/
# ./SomeDirectory/rte-rrtmgp/
#

include("ConfigRRTMGP.jl")

cd(joinpath(root_dir,"rte-rrtmgp","build")) do
  rm("Makefile.libs"; force=true)
  rm("Makefile.conf"; force=true)
  rm("Makefile.rules"; force=true)
end

cd(joinpath(root_dir,"rte-rrtmgp","build")) do
  run(`make clean`)
end

cd(joinpath(root_dir,"rte-rrtmgp")) do
  rm(".DS_store"; force=true)
end
cd(joinpath(root_dir,"rte-rrtmgp","examples")) do
  rm(".DS_store"; force=true)
end

cd(rte_rrtmgp_rfmip_clear_sky) do
  conds_file = joinpath(".", "multiple_input4MIPs_radiation_RFMIP_UColorado-RFMIP-1-2_none.nc")
  rm(conds_file; force=true)
  rm(".DS_store"; force=true)
  rm("generate-output-file-templates.py"; force=true)
  rm("rld_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force=true)
  rm("rlu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force=true)
  rm("rrtmgp_rfmip_lw"; force=true)
  rm("rrtmgp_rfmip_sw"; force=true)
  rm("rsd_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force=true)
  rm("rsu_Efx_RTE-RRTMGP-181204_rad-irf_r1i1p1f1_gn.nc"; force=true)
  run(`make clean`)
end

cd(rte_rrtmgp_rfmip_all_sky) do
  run(`make clean`)
end


text:
@doc doc"""
    LogEuclideanMetric <: Metric

The LogEuclidean Metric consists of the Euclidean metric applied to all elements after mapping them
into the Lie Algebra, i.e. performing a matrix logarithm beforehand.
"""
struct LogEuclideanMetric <: RiemannianMetric end

@doc doc"""
    distance(M::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric}, x, y)

Compute the distance on the [`SymmetricPositiveDefinite`](@ref) manifold between
`x` and `y` as a [`MetricManifold`](@ref) with [`LogEuclideanMetric`](@ref).
The formula reads

```math
    d_{\mathcal P(n)}(x,y) = \lVert \Log x - \Log y \rVert_{\mathrm{F}}
```

where $\operatorname{Log}$ denotes the matrix logarithm and
$\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
"""
function distance(
    M::MetricManifold{SymmetricPositiveDefinite{N},LogEuclideanMetric},
    x,
    y,
) where {N}
    return norm(log(Symmetric(x)) - log(Symmetric(y)))
end


text:
using FundamentalsNumericalComputation

fncdir = dirname(pathof(FNC))

for root in ["intro","linsys","leastsq","nonlineqn","localapprox","ivp","appendix"]
	for mdfile in filter(s->endswith(s,".md"),readdir(root))
		str = readlines(joinpath(root,mdfile))
		itext = findfirst(occursin.("{proof:function}",str))
		offset = 0
		while !isnothing(itext)
			itext += offset
			funcname = match(r"`+{proof:function} (.*)",str[itext]).captures[1]
			@show root,mdfile,itext,funcname

			funcstr = ""

			for srcfile in filter(s->startswith(s,"chapter"),readdir(fncdir))
				src = readlines(joinpath(fncdir,srcfile))
				isrc = findfirst(occursin.("function $funcname",src))
				isnothing(isrc) && continue

				k = isrc + findfirst(occursin.("return",src[isrc+1:end]))
				endsrc = k + findfirst(startswith.(src[k+1:end],"end"))

				k = findlast(startswith.(src[1:isrc-1],"\"\"\""))
				startsrc = findlast(startswith.(src[1:k-1],"\"\"\""))

				funcstr = src[startsrc:endsrc]
				break
			end

			n = itext + findfirst(occursin.("{code-block}",str[itext+1:end]))
			offset = n + findfirst(startswith.(str[n+1:end],"```"))

			newstr = [":lineno-start: 1";funcstr]
			str = [str[1:n];newstr;str[offset:end]]
			offset = n + length(newstr) + 1

			itext = findfirst(occursin.("{proof:function}",str[offset+1:end]))
		end

		open(joinpath(root,"$mdfile"),"w") do io
			[ println(io,s) for s in str ];
		end
	end
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
# ConvDual, Duality
# Problem spec - input: Hyperrectangle with uniform radisu, outputL half space constraint

using NeuralVerification
using Test

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

at = @__DIR__

# Problem type - input:Hyperrectangle, output:Hyperrectangle
print("###### Group 3 - input:Hyperrectangle, output:Hpolytope (with 1 constraint) ######\n")
# Small MNIST Network

mnist_small = read_nnet("$at/../examples/networks/mnist_small.nnet")

# entry 23 in MNIST datset
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
#output_center = [-1311.1257826380004,4633.767704436501,-654.0718535670002,-1325.349417307,1175.2361184373997,-1897.8607293569007,-470.3405972940001,830.8337987382,-377.7467076115001,572.3674015264198]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

#output_low = output_center .- out_epsilon
#output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
#outputSet = Hyperrectangle(low=output_low, high=output_high)

A = Matrix(undef, 2, 1)
A = [1.0, -1, 0, 0, 0, 0, 0, 0, 0 ,0]'
b = [0.0]
outputSet = HPolytope(A, b)

########### Small ############
# Problem type - input:Hyperrectangle, output:Hpolytope (one constraint)

problem_hyperrect_halfspace_small = Problem(mnist_small, inputSet, outputSet)
print("\n\n################## Small ##################\n")

# Duality
print("\nDuality - Small")
optimizer = GLPKSolverMIP()
solver = Duality(optimizer)
#timed_result = @timed solve(solver, problem_hyperrect_halfspace_small)
t = timed_result[2]
out = timed_result[1]
#print(" - Time: $t s - Output: $out")

# ConvDual
print("\nConvDual - Small")
optimizer = GLPKSolverMIP()
solver = ConvDual()
timed_result = @timed solve(solver, problem_hyperrect_halfspace_small)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

########### Deep ############
print("\n\n################## Deep ##################\n")

mnist_deep = read_nnet("$at/../examples/networks/mnist_large.nnet")
problem_hyperrect_halfspace_deep = Problem(mnist_deep, inputSet, outputSet)

# Duality
print("\nDuality - Deep")
optimizer = GLPKSolverMIP()
solver = Duality(optimizer)
#timed_result = @timed solve(solver, problem_hyperrect_halfspace_deep)

# ConvDual
print("\nConvDual - Deep")
optimizer = GLPKSolverMIP()
solver = ConvDual()
timed_result = @timed solve(solver, problem_hyperrect_halfspace_deep)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")

########### Wide ############
print("\n\n################## Wide ##################\n")

mnist_wide = read_nnet("$at/../examples/networks/mnist-1-100.nnet")
problem_hyperrect_halfspace_wide = Problem(mnist_wide, inputSet, outputSet)

# Duality
print("\nDuality - Wide")
optimizer = GLPKSolverMIP()
solver = Duality(optimizer)
#timed_result = @timed solve(solver, problem_hyperrect_halfspace_deep)
t = timed_result[2]
out = timed_result[1]
#print(" - Time: $t s - Output: $out")

# ConvDual
print("\nConvDual - Wide")
optimizer = GLPKSolverMIP()
solver = ConvDual()
timed_result = @timed solve(solver, problem_hyperrect_halfspace_deep)
t = timed_result[2]
out = timed_result[1]
print(" - Time: $t s - Output: $out")


text:
push!(LOAD_PATH,"../src/")

using Documenter, LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using CAGD

makedocs(
	sitename = "CSG.jl ",
	assets = ["assets/lar.css", "assets/logo.png"],
	format = Documenter.HTML(
		prettyurls = get(ENV, "CI", nothing) == "true"
	),
	pages = [
		"Home" => "index.md",
		"Tutorial" => [
			"Model" => "model.md",
			"Boolean forms" => "boolean.md",
			"Examples" => "examples.md"
		],
		"Computational pipeline" => [
			"Spatial indices" => "spatial_indices.md",
			"Planar arrangement" => "planar_arrangement.md",
			"Splitting" => "splitting.md",
			"Congruence" => "congruence.md",
			"Topological Gift Wrapping" => "tgw.md",
			"Spatial_arrangement" => "spatial_arrangement.md",
			"Shells" => "shells.md"
		],
		"Utils" => "utils.md"
	]
)


text:
@testset "Lines2D" begin

	@testset "LineSegment" begin
		img = zeros(Gray{N0f8}, 10, 10)
		expected = copy(img)
		expected[diagind(expected)] .= oneunit(Gray{N0f8})
		res = @inferred draw(img, LineSegment(1,1,10,10))
		@test expected == res
		res = @inferred draw(img, LineSegment(1,1,10,10), bresenham)
		@test expected == res
		res = @inferred draw(img, LineSegment(1,1,10,10), oneunit(Gray), bresenham)
		@test expected == res
		res = @inferred draw(img, LineSegment(1,1,10,10), oneunit(Gray))
		@test expected == res
		res = @inferred draw(img, LineSegment(CartesianIndex(1, 1), CartesianIndex(10, 10)))
		@test expected == res
		res = @inferred draw(img, LineSegment(CartesianIndex(1, 1), CartesianIndex(10, 10)), bresenham)
		@test expected == res
		res = @inferred draw(img, LineSegment(CartesianIndex(1, 1), CartesianIndex(10, 10)), oneunit(Gray), bresenham)
		@test expected == res
		res = @inferred draw(img, LineSegment(CartesianIndex(1, 1), CartesianIndex(10, 10)), oneunit(Gray))
		@test expected == res
		expected[diagind(expected)] .= Gray{N0f8}(0.5)
		res = @inferred draw(img, LineSegment(1,1,10,10), Gray{N0f8}(0.5))
		@test expected == res
		img2 = zeros(RGB, 10, 10)
		expected = copy(img2)
		expected[diagind(expected)] .= RGB{N0f8}(0.2, 0.3, 0.4)
		res = @inferred draw(img2, LineSegment(1,1,10,10), RGB{N0f8}(0.2,0.3,0.4))
		@test expected == res
	end

	@testset "Bresenham" begin
		img = zeros(Gray{N0f8}, 10, 10)
		expected = copy(img)
		expected[diagind(expected)] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 1, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[diagind(expected)[1:5]] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 1, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[10, :] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 10, 1, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[5, :] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 5, 1, 5, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[5, 1:5] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 5, 1, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[:, 10] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 10, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[:, 5] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 5, 10, 5, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[1:5, 5] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 5, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[10:9:91] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 1, 10, 10, 1, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[10:9:55] .= oneunit(Gray{N0f8})
		res = @inferred bresenham(copy(img), 10, 1, 5, 6, oneunit(Gray{N0f8}))
		@test expected == res

		img2 = zeros(RGB, 10, 10)
		expected = copy(img2)
		expected[diagind(expected)] .= RGB{N0f8}(0.2, 0.3, 0.4)
		res = @inferred bresenham(copy(img2), 1, 1, 10, 10, RGB{N0f8}(0.2, 0.3, 0.4))
		@test expected == res
	end

	@testset "Xiaolin-Wu" begin
		img = zeros(Gray{N0f8}, 10, 10)
		expected = copy(img)
		expected[diagind(expected)] .= oneunit(Gray{N0f8})
		expected[1, 1] = Gray{N0f8}(0.5)
		expected[10, 10] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 1, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[diagind(expected)[1:5]] .= oneunit(Gray{N0f8})
		expected[1, 1] = Gray{N0f8}(0.5)
		expected[5, 5] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 1, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[10, :] .= oneunit(Gray{N0f8})
		expected[10, 1] = Gray{N0f8}(0.5)
		expected[10, 10] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 10, 1, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[5, :] .= oneunit(Gray{N0f8})
		expected[5, 1] = Gray{N0f8}(0.5)
		expected[5, 10] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 5, 1, 5, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[5, 1:5] .= oneunit(Gray{N0f8})
		expected[5, 1] = Gray{N0f8}(0.5)
		expected[5, 5] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 5, 1, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[:, 10] .= oneunit(Gray{N0f8})
		expected[1, 10] = Gray{N0f8}(0.5)
		expected[10, 10] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 10, 10, 10, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[:, 5] .= oneunit(Gray{N0f8})
		expected[1, 5] = Gray{N0f8}(0.5)
		expected[10, 5] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 5, 10, 5, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[1:5, 5] .= oneunit(Gray{N0f8})
		expected[1, 5] = Gray{N0f8}(0.5)
		expected[5, 5] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 5, 5, 5, oneunit(Gray{N0f8}))
		@test expected == res

		expected = copy(img)
		expected[10:9:91] .= oneunit(Gray{N0f8})
		expected[10, 1] = Gray{N0f8}(0.5)
		expected[1, 10] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 1, 10, 10, 1, oneunit(Gray{N0f8}))
		@test expected == res
		expected = copy(img)
		expected[10:9:55] .= oneunit(Gray{N0f8})
		expected[10, 1] = Gray{N0f8}(0.5)
		expected[5, 6] = Gray{N0f8}(0.5)
		res = @inferred xiaolin_wu(copy(img), 10, 1, 5, 6, oneunit(Gray{N0f8}))
		@test expected == res
	end

	@testset "LineNormal" begin
		img = zeros(Gray, 10, 10)
		expected = copy(img)
		expected[1:10, 1] .= 1
		@inferred draw!(img, LineNormal(1, 0))
		@test img == expected
		fill!(img,zero(Gray))
		res = @inferred draw(img, LineNormal(1, 0))
		@test expected == res
		expected[1:10, 1] .= Gray(0.5)
		res = @inferred draw(img, LineNormal(1, 0), Gray(0.5))
		@test expected == res
		fill!(expected,zero(Gray))
		expected[1:10,2] .= expected[1:10,7] .= expected[10,1:10] .= Gray(0.3)
		@inferred draw!(img, LineNormal(2, 0), Gray(0.3))
		@inferred draw!(img, LineNormal(7, 0), Gray(0.3))
		@inferred draw!(img, LineNormal(10, π/2), Gray(0.3))
		@test img == expected
		fill!(img,zero(Gray))
		fill!(expected,zero(Gray))
		for i in 1:10
			expected[i,10-i+1] = 1
		end
		res = @inferred draw(img, LineNormal(7.5, π/4))
		@test expected == res
		for i in 1:10
			expected[i,10-i+1] = 0.8
		end
		fill!(expected, Gray(0))
		for i in 1:10
			expected[i,i] = Gray(0.8)
		end
		res = @inferred draw(img, LineNormal((0,-π/4)), Gray(0.8))
		@test expected == res
		img = zeros(Gray,3,3)
		expected = copy(img)
		res = @inferred draw(img, LineNormal(4,0))
		@test expected == res
		res = @inferred draw(img, LineNormal(4,π/2))
		@test expected == res
		res = @inferred draw(img, LineNormal(5,π/4))
		@test expected == res
	end

	@testset "LineTwoPoints" begin
		img = zeros(Gray{N0f8}, 10, 10)
		expected = copy(img)
		expected[diagind(expected)] .= oneunit(Gray{N0f8})
		res = @inferred draw(img, LineTwoPoints(1,1,10,10))
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(1,1,10,10), bresenham)
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(1,1,10,10), oneunit(Gray), bresenham)
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(1,1,10,10), oneunit(Gray))
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(CartesianIndex(1, 1), CartesianIndex(10, 10)))
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(CartesianIndex(1, 1), CartesianIndex(10, 10)), bresenham)
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(CartesianIndex(1, 1), CartesianIndex(10, 10)), oneunit(Gray), bresenham)
		@test expected == res
		res = @inferred draw(img, LineTwoPoints(CartesianIndex(1, 1), CartesianIndex(10, 10)), oneunit(Gray))
		@test expected == res
		expected[diagind(expected)] .= Gray{N0f8}(0.5)
		res = @inferred draw(img, LineTwoPoints(1,1,10,10), Gray{N0f8}(0.5))
		@test expected == res
		img2 = zeros(RGB, 10, 10)
		expected = copy(img2)
		expected[diagind(expected)] .= RGB{N0f8}(0.2, 0.3, 0.4)
		res = @inferred draw(img2, LineTwoPoints(1,1,10,10), RGB{N0f8}(0.2,0.3,0.4))
		@test expected == res
		@test img2 == @inferred draw(img2, LineTwoPoints(1,11,10,11))
	end
end


text:
# Nominal Models
@testset "Nominal Models" begin
    data = CSV.read(joinpath(pkgdir(Econometrics), "data", "insure.csv"), DataFrame, select = [:insure, :age, :male, :nonwhite, :site]) |>
        dropmissing |>
        (data -> transform!(data, [:insure, :site] .=> categorical, renamecols = false))
    model = fit(
        EconometricModel,
        @formula(insure ~ age + male + nonwhite + site),
        data,
        contrasts = Dict(:insure => DummyCoding(base = "Uninsure")),
        )
    @test sprint(show, model) ==
          "Probability Model for Nominal Response\nCategories: Uninsure, Indemnity, Prepaid\nNumber of observations: 615\nNull Loglikelihood: -555.85\nLoglikelihood: -534.36\nR-squared: 0.0387\nLR Test: 42.99 ∼ χ²(10) ⟹ Pr > χ² = 0.0000\nFormula: insure ~ 1 + age + male + nonwhite + site\n───────────────────────────────────────────────────────────────────────────────────────────────────\n                                       PE         SE       t-value  Pr > |t|       2.50%     97.50%\n───────────────────────────────────────────────────────────────────────────────────────────────────\ninsure: Indemnity ~ (Intercept)   1.28694     0.59232     2.17271     0.0302   0.123689   2.4502\ninsure: Indemnity ~ age           0.00779612  0.0114418   0.681372    0.4959  -0.0146743  0.0302666\ninsure: Indemnity ~ male         -0.451848    0.367486   -1.22957     0.2193  -1.17355    0.269855\ninsure: Indemnity ~ nonwhite     -0.217059    0.425636   -0.509965    0.6103  -1.05296    0.618843\ninsure: Indemnity ~ site: 2       1.21152     0.470506    2.57493     0.0103   0.287497   2.13554\ninsure: Indemnity ~ site: 3       0.207813    0.366293    0.56734     0.5707  -0.511547   0.927172\ninsure: Prepaid ~ (Intercept)     1.55666     0.596327    2.61041     0.0093   0.385533   2.72778\ninsure: Prepaid ~ age            -0.00394887  0.0115993  -0.340439    0.7336  -0.0267287  0.018831\ninsure: Prepaid ~ male            0.109846    0.365187    0.300793    0.7637  -0.607343   0.827035\ninsure: Prepaid ~ nonwhite        0.757718    0.419575    1.80592     0.0714  -0.0662835  1.58172\ninsure: Prepaid ~ site: 2         1.32456     0.469789    2.81947     0.0050   0.401941   2.24717\ninsure: Prepaid ~ site: 3        -0.380175    0.372819   -1.01973     0.3083  -1.11235    0.352001\n───────────────────────────────────────────────────────────────────────────────────────────────────"
    β, V, σ = coef(model), vcov(model), stderror(model)
    @test β ≈ [
        1.286943,
        0.0077961,
        -0.4518496,
        -0.2170589,
        1.211563,
        0.2078123,
        1.556656,
        -0.0039489,
        0.1098438,
        0.7577178,
        1.324599,
        -0.3801756,
        ] rtol = 1e-3
    @test V ≈
          [
        0.35084518 -0.00590456 -0.02519209 -0.01837512 -0.08995794 -0.08264268 0.29928937 -0.00510107 -0.02257270 -0.01585492 -0.07392512 -0.06808314
        -0.00590456 0.00013092 -0.00046225 -0.00028066 0.00041004 0.00040470 -0.00509454 0.00011354 -0.00039387 -0.00024122 0.00033847 0.00035236
        -0.02519209 -0.00046225 0.13504644 0.01565749 0.01009559 0.00933967 -0.02288847 -0.00039497 0.11365141 0.01369872 0.00945927 0.00802593
        -0.01837512 -0.00590456 -0.02519209 0.18116609 0.01886597 -0.02821348 -0.01638629 -0.00023179 0.01314666 0.15068111 0.01664999 -0.02345786
        -0.08995794 -0.00590456 -0.02519209 -0.01837512 0.22138224 0.064050660 -0.07438947 0.00035751 0.00845363 0.01594691 0.19895500 0.05189742
        -0.08264268 0.00040470 0.00933967 -0.02821348 0.06405066 0.134170290 -0.06796097 0.00034828 0.00826890 -0.02353016 0.05182727 0.11060487
        0.29928937 -0.00509454 -0.02288847 -0.01638629 -0.07438947 -0.06796097 0.35560782 -0.00601966 -0.02683161 -0.02105510 -0.08991231 -0.08165509
        -0.00510107 0.00011354 -0.00039497 -0.00023179 0.00035751 0.00034828 -0.00601966 0.00013455 -0.00047223 -0.00028680 0.00039990 0.0004179
        -0.02257270 -0.00039387 0.11365141 0.01314666 0.00845363 0.00826890 -0.02683161 -0.00047223 0.13336253 0.01625685 0.01101870 0.00922333
        -0.01585492 -0.00024122 0.01369872 0.15068111 0.01594691 -0.02353016 -0.02105510 -0.00028680 0.01625685 0.17604391 0.02099458 -0.03008391
        -0.07392512 0.00033847 0.00945927 0.01664999 0.19895500 0.05182727 -0.08991231 0.00039990 0.01101870 0.02099458 0.22070774 0.06273376
        -0.06808314 0.00035236 0.00802593 -0.02345786 0.05189742 0.11060487 -0.08165509 0.00041790 0.00922333 -0.03008391 0.06273376 0.13899387
        ] |>
        Hermitian rtol = 1e-3
    @test σ ≈ [
        0.5923219,
        0.0114418,
        0.3674867,
        0.4256361,
        0.4705127,
        0.3662926,
        0.5963286,
        0.0115994,
        0.3651883,
        0.4195759,
        0.4697954,
        0.3728188,
        ] rtol = 1e-4
end


text:
using DrWatson
@quickactivate "GCN_HM_GRN-Integration"
DrWatson.greet()