text:
using Plasm
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using NearestNeighbors

#	transform sigma and related faces in space_idx
#	function face_mapping(V, FV, sigma, err=LinearAlgebraicRepresentation.ERR )
#		sigmavs = FV[sigma]; i = 1
#		compute affinely independent triple
#		while (-err < det(V[:, sigmavs[i:i+2]]) < err) && (i < length(sigmavs)-2)
#			i += 1
#		end
#		sigma translation
#		T = Lar.t(-V[:,sigmavs[i+1]]...)
#		W = T * [V[:,sigmavs] ; ones(1,length(sigmavs))]
#		t_vs = W[1:3,:]
#		translated sigma rotation
#		r1,r2,r3 = t_vs[:,i], t_vs[:,i+2], cross(t_vs[:,i], t_vs[:,i+2])
#		R = eye(4)
#		R[1:3,1:3] = inv([r1 r2 r3])
#		mapping = R * T
#		return mapping
#	end

# transform sigma and related faces in space_idx
function face_mapping(V, FV, sigma, err=LinearAlgebraicRepresentation.ERR )
	vs = FV[sigma]; i = 1
	# compute affinely independent triple
	n = length(vs)
	succ(i) = i % n + 1
	a = normalize(V[:,vs[succ(i)]] - V[:,vs[i]])
	b = normalize(V[:,vs[succ(succ(i))]] - V[:,vs[i]])
	c = cross(a,b)
	while (-err < det([a b c]) < err)
		i += 1
		a = normalize(V[:,vs[succ(i)]] - V[:,vs[i]])
		b = normalize(V[:,vs[succ(succ(i))]] - V[:,vs[i]])
		c = cross(a,b)
	end
	# sigma translation
	T = Lar.t(-V[:,vs[succ(i)]]...)
	W = T * [V[:,vs] ; ones(1,length(vs))]
	t_vs = W[1:3,:]
	# translated sigma rotation
	r1,r2,r3 = t_vs[:,i], t_vs[:,succ(succ(i))], cross(t_vs[:,i], t_vs[:,succ(succ(i))])
	eye(n) = LinearAlgebra.Matrix{Float64}(I,n,n)
	R = eye(4)
	R[1:3,1:3] = inv([r1 r2 r3])
	mapping = R * T
	return mapping
end

function polyline(V::Lar.Points)::Lar.LAR
	return V,[[k,k+1] for k=1:size(V,2)-1]
end

function sigmamodel(V, copEV, FV, copFE, sigma, space_idx)
	bigpi = space_idx[sigma]
	Q = Lar.face_mapping(V', FV, sigma)
	EV = [findnz(copEV[e,:])[1] for e=1:size(copEV,1)]
	FE = [findnz(copFE[e,:])[1] for e=1:size(copFE,1)]
	sigma_edges = EV[FE[sigma]]
	sigma_verts = union(sigma_edges...)
	sigma_vdict =  Dict(zip(sigma_verts, 1:length(sigma_verts)))
	sigma_lines = [[sigma_vdict[v] for v in edge] for edge in sigma_edges]
	# sigma face on z=0 plane
	S = V'[:,sigma_verts] # sigma vertices
	Z = (Q * [S; ones(1,size(S,2))])[1:3,:] # sigma mapped in z=0
	v,ev = Z,sigma_lines
	#Plasm.view(Plasm.numbering(0.5)((v,[[[k] for k=1:size(v,2)],ev])))
	return v,ev,Q
end

function sigma_intersect(V, copEV, copFE, sigma, Q, bigpi)
	EV = [findnz(copEV[e,:])[1] for e=1:size(copEV,1)]
	FE = [findnz(copFE[f,:])[1] for f=1:size(copFE,1)]
	sigma_edges = EV[FE[sigma]]
	sigma_verts = union(sigma_edges...)
	sigma_vdict =  Dict(zip(sigma_verts, 1:length(sigma_verts)))
	sigma_lines = [[sigma_vdict[v] for v in edge] for edge in sigma_edges]
	# sigma face on z=0 plane
	S = V[:,sigma_verts] # sigma vertices
	Z = (Q * [S; ones(1,size(S,2))])[1:3,:] # sigma mapped in z=0
	#Plasm.view(Z, sigma_lines)
	# initialization of line storage (to be intersected)
	linestore = [[Z[:,v1] Z[:,v2]] for (v1,v2) in sigma_lines]
	linenum = length(sigma_lines)
	# reindexing of ∏(σ) chain
	bigpi_edges = [EV[fe] for fe in FE[bigpi]]
	bigpi_verts = union(union(bigpi_edges...)...)
	bigpi_vdict =  Dict(zip(bigpi_verts, 1:length(bigpi_verts)))
	bigpi_lines = [[sort([bigpi_vdict[v] for v in edge]) for edge in faceedges]
		for faceedges in bigpi_edges]
	# bigpi trasformed in 3D according to Q mapping
	P = V[:,bigpi_verts] # bigpi vertices
	W = (Q * [P; ones(1,size(P,2))])[1:3,:] # bigpi mapped by Q
	#Plasm.view(W, union(bigpi_lines...))
	# filter on bigpi_lines that do not cross z=0
	filtered_edges = [[[v1,v2] for (v1,v2) in face_lines
		if (sign(W[3,v1]) * sign(W[3,v2])) <= 0]
			for face_lines in bigpi_lines]
	filtered_edges = [edge for edge in filtered_edges if edge!=[]]
	#Plasm.view(W, union(filtered_edges...))
	#Plasm.view(Plasm.numbering()((W,[[[k] for k=1:size(W,2)],union(filtered_edges...)])))
	# computation of ordered z=0 points by face
	facepoints = []
	for face in filtered_edges
		points = []
		for (v1,v2) in face
			#v1[3] + (v2[3]-v1[3])*t = 0 # compute parameter of intersection point with z=0
			if (W[3,v2] - W[3,v1]) != 0
				t = -W[3,v1] / (W[3,v2] - W[3,v1])  # z1 + (z2-z1)*t = 0
				point = W[:,v1] + (W[:,v2] - W[:,v1]) * t
				push!(points,point)
			end
		end
		if length(Set(points)) > 1
			push!(facepoints,Set(points))
		end
		# TODO: linearly order facepoints[face]
	end
	# compute z_lines
	c = collect
	z0_lines = [[[c(ps)[k] c(ps)[k+1]] for k=1:2:(length(ps)-1)] for ps in facepoints]

	# intersecting lines upload in linestore
	linestore = [[Z[:,v1] Z[:,v2]] for (v1,v2) in sigma_lines]
	if z0_lines != []
		linestore = append!(linestore, union(z0_lines...))
	end
	return linestore,linenum
end

function computeparams(linestore,linenum)
	m = length(linestore)
	params = [[] for i=1:m]
	for h=1:(m-1)
		(x1,y1,x2,y2) = reshape(linestore[h][1:2,:],1,4)
		for k=(h+1):m
			(x3,y3,x4,y4) = reshape(linestore[k][1:2,:],1,4)
			# intersect lines h,k
			# http://www.cs.swan.ac.uk/~cssimon/line_intersection.html
			det = (x4-x3)*(y1-y2)-(x1-x2)*(y4-y3)
			α,β = 99,99
			if det != 0.0
				a = 1/det
				b = [y1-y2 x2-x1; y3-y4 x4-x3]  # x1-x2 => x2-x1 bug in the source link !!
				c = [x1-x3; y1-y3]
				(β,α) = a * b * c
			else
				if (y1==y2) == (y3==y4) # segments collinear
					α = -x1/(x2-x1)
					β = -x3/(x4-x3)
				elseif (x1==x2) == (x3==x4)
					α = -y1/(y2-y1)
					β = -y3/(y4-y3)
				else
					 # segments parallel: no intersection
					 (β,α) = 999
				end
			end
			if 0<=α<=1 && 0<=β<=1
				push!(params[h], α)
				push!(params[k], β)
			end
		end
	end
	out = []
	for line in params
		push!(line, 0.0)
		push!(line, 1.0)
		line = sort(collect(Set(line)))
		push!(out, line)
	end
	return out
end

function fragface(V, EV, FV, FE, space_idx, sigma)
	bigpi = space_idx[sigma]
	Q = Lar.face_mapping(V', FV, sigma)
	linestore, linenum = Lar.sigma_intersect(V', EV, FE, sigma, Q, bigpi)
	lineparams = Lar.computeparams(linestore,linenum)
	pairs = collect(zip(lineparams,linestore))

	linepoints = []
	number = 0
	for (k,(params,line)) in enumerate(pairs)
		v1,v2 = line[1:2,1],line[1:2,2]
		if params==[0.0, 1.0]
			points = [v1,v2]
			if k<=linenum number += 1 end
		elseif length(params)>2
			points = [v1+α*(v2-v1) for α ∈ params]
			if k<=linenum number += length(params)-1 end
		end
		push!(linepoints,points)
	end
	#linepoints = collect(Set(linepoints))
	verts = []
	edges = Array{Array{Int},1}()
	offset = 0
	for (h,pointline) in enumerate(linepoints)
		append!(verts, pointline)
		append!(edges, [[k,k+1] for k=1:length(pointline)-1]+offset)
		offset = length(verts)
	end
	verts = hcat(verts...)
	edges = convert(Lar.Cells, edges)
	return verts, edges, number
end

function removevertices(verts, edges, err=LinearAlgebraicRepresentation.ERR)
	vertices = collect(Set(vcat(edges...)))
	vertexnum = length(vertices)
	vdict = Dict(zip(vertices,1:vertexnum))
	cells0D = cells0D = zeros(2,vertexnum)
	cells1D = Array{Int64,1}[]
	for (v1,v2) in edges
		push!(cells1D,[vdict[v1],vdict[v2]])
	end
	#vdict = Dict(zip(values(vdict),keys(vdict))) # inverted vdict
	for key in keys(vdict)
		cells0D[:,vdict[key]] = verts[:,key]
	end
	return cells0D, cells1D
end

function merge_vertices(rV::Lar.Points, rEV::Lar.ChainOp, rFE::Lar.ChainOp,
err=Lar.ERR)
	verts = rV'
	edges = [findnz(rEV[e,:])[1] for e=1:size(rEV,1)]
	faces = [findnz(rFE[f,:])[1] for f=1:size(rFE,1)]

	function make_quotients(verts, faces, edges, err=Lar.ERR)
		PRECISION = 5.0

		vertsnum = size(verts, 2)
		kdtree = NearestNeighbors.KDTree(verts)
		todelete = []
		newverts = zeros(Int, vertsnum)
		i = 1
		for vi in 1:vertsnum
			if !(vi in todelete)
				nearvs = NearestNeighbors.inrange(kdtree, verts[:,vi], err)
				newverts[nearvs] = i
				nearvs = setdiff(nearvs, vi)
				todelete = union(todelete, nearvs)
				i = i + 1
			end
		end
		nverts = verts[:,setdiff(collect(1:vertsnum), todelete)]
		vertdict = DataStructures.OrderedDict()
		for k=1:size(nverts,2)
			key = map(Lar.approxVal(PRECISION), nverts[:,k])
			vertdict[key] = k
		end

		nedges = Array{Array{Int},1}()
		for (v1,v2) in edges
			(w1,w2) = abs.((v1,v2))
			key1 = map(Lar.approxVal(PRECISION), verts[:,w1])
			key2 = map(Lar.approxVal(PRECISION), verts[:,w2])
			push!(nedges, [vertdict[key1], vertdict[key2]])
		end
		nedges = Set(map(sort,nedges))
		nedges = sort(collect(nedges), by=x->(x[1],x[2]))
		nedges = filter( x->x[1]!=x[2], nedges )
#		nedges = sort(nedges, by=x->(x[1],x[2]))
#		nedges = sort(collect(Set(nedges)), by=x->(x[1],x[2]))
#		nedges = filter( x->x[1]!=x[2], nedges )
#		nedges = convert(Array{Array{Int64,1},1}, nedges)

		nfaces = Array{Array{Int},1}()
		for (k,face) in enumerate(faces)
			nface = Int64[]
			for edge in face
				(v1,v2) = edges[edge]
				(w1,w2) = abs.((v1,v2))
				key1 = map(Lar.approxVal(PRECISION), verts[:,w1])
				key2 = map(Lar.approxVal(PRECISION), verts[:,w2])
				push!(nface, vertdict[key1])
				push!(nface, vertdict[key2])
			end
			if nface ≠ []
				nface = sort(collect(Set(nface)))
				push!(nfaces,nface)
			end
		end
		nfaces = collect(Set(nfaces))
		nfaces = sort(nfaces, lt=lexless)
		nfaces = convert(Lar.Cells, nfaces)
		return nverts, nfaces, nedges
	end

	V,FV,EV = make_quotients(verts, faces, edges, err)
	return V,FV,EV
end

function (verts, edges, err=LinearAlgebraicRepresentation.ERR)
	vertsnum = size(verts, 2)
	kdtree = KDTree(verts)
	todelete = []
	newverts = zeros(Int, vertsnum)
	i = 1
	for vi in 1:vertsnum
	   if !(vi in todelete)
		   nearvs = inrange(kdtree, verts[:,vi], err)
		   newverts[nearvs] = i
		   nearvs = setdiff(nearvs, vi)
		   todelete = union(todelete, nearvs)
		   i = i + 1
	   end
	end
	nverts = verts[:,setdiff(collect(1:vertsnum), todelete)]
	vertdict = DataStructures.OrderedDict()
	for k=1:size(nverts,2)
		key = map(Lar.approxVal(PRECISION), nverts[:,k])
		vertdict[key] = k
	end
	nedges = Array{Array{Int},1}()
	for (v1,v2) in edges
		(w1,w2) = abs.((v1,v2))
		key1 = map(Lar.approxVal(PRECISION), verts[:,w1])
		key2 = map(Lar.approxVal(PRECISION), verts[:,w2])
		push!(nedges, [vertdict[key1], vertdict[key2]])
	end
	nedges = Set(map(sort,nedges))
	nedges = sort(collect(nedges), by=x->(x[1],x[2]))
	nedges = filter( x->x[1]!=x[2], nedges )
	return nverts,nedges
end

function preprocessing(V,copEV,FV)
	EV = [findnz(copEV[e,:])[1] for e=1:size(copEV,1)]
	#copEV = Lar.coboundary_0(EV)
	#copFE = build_copFE(FV, EV)
	copFE = Lar.build_copFE(FV, EV)
	copEF = copFE'
	EF = [findnz(copEF[e,:])[1] for e=1:size(copEF,1)]
	FE = [findnz(copFE[e,:])[1] for e=1:size(copFE,1)]

	sp_idx = Lar.Arrangement.spatial_index(V, copEV, copFE)
	space_idx = Array{Int,1}[]
	for sigma = 1:length(FV)
		edges = FE[sigma];
		# remove faces incident on sigma edges
		facesofedges = Set(vcat([EF[edge] for edge in edges]...));
		push!(space_idx, setdiff(sp_idx[sigma], facesofedges))
	end
	return V, copEV, copFE, space_idx
end

function computefragments(V, EV, FV, FE, space_idx, sigma)
	#v,ev = Lar.sigmamodel(V, EV, FE, sigma, space_idx)
	#Plasm.view(Plasm.numbering(0.4)((v,[ [[k] for k=1:size(v,2)],ev ])))

	verts, edges, linenum = Lar.fragface(V, EV, FV, FE, space_idx, sigma) #OK
	v,ev,Q = Lar.sigmamodel(V, EV, FV, FE, sigma, space_idx) #OK
	classify = Lar.pointInPolygonClassification(v,ev)
	internal = []
	for (k,edge)  in enumerate(edges)
		v1,v2 = edge
		if k <= linenum
			push!(internal, edge)
		else
			vm = (verts[:,v1]+verts[:,v2])/2
			if (classify(vm)=="p_in")
				push!(internal, edge)
			end
		end
	end

	nverts, nedges = Lar.removevertices(verts, internal)
	w, ev = Lar.mergevertices(nverts, nedges)
	cop_ev = Lar.build_copEV(ev,true)
	cop_fe = Lar.Arrangement.minimal_2cycles(w',cop_ev)	#OK
	shell_num = Lar.Arrangement.get_external_cycle(w', cop_ev, cop_fe)
	for j=1:cop_fe.n
		cop_fe[shell_num,j]=0
	end
	cop_fe = dropzeros(cop_fe)
	faces = [findnz(cop_fe[e,:])[1] for e=1:size(cop_fe,1) if e≠shell_num]
	fv = [collect(Set(vcat([ev[e] for e in face]...))) for face in faces]
	vv = [[k] for k=1:size(w,2)]
	#Plasm.view(Plasm.numbering(0.4)((w, [vv, ev, fv])))

	w = [ w; zeros(size(w,2))'; ones(size(w,2))' ]
	v = (inv(Q) * w)[1:3,:]
	vv = [[k] for k=1:size(w,2)]
	#Plasm.view(Plasm.numbering(0.4)((v, [vv, ev, fv])))
	return v', Lar.coboundary_0(ev), cop_fe
end


text:
using Snopt
using Test
using SparseArrays

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

@testset "rosenbrock" begin

function rosenbrock(g, df, dg, x, deriv)
    f = (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2
    fail = false

    if deriv
        df[1] = -2*(1 - x[1]) + 200*(x[2] - x[1]^2)*-2*x[1]
        df[2] = 200*(x[2] - x[1]^2)
    end

    return f, fail
end

x0 = [4.0; 4.0]
lx = [-5.0; -5.0]
ux = [5.0; 5.0]
lg = []
ug = []
rows = []
cols = []

xopt, fopt, info, out = snopta(rosenbrock, x0, lx, ux, lg, ug, rows, cols)

@test isapprox(xopt[1], 1.0; atol=1e-6)
@test isapprox(xopt[2], 1.0; atol=1e-6)
@test isapprox(fopt, 0.0; atol=1e-8)
@test info == "Finished successfully: optimality conditions satisfied"

end  # rosenbrock test set

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

@testset "barnes" begin

function barnes(g, df, dg, x, deriv)

    a1 = 75.196
    a3 = 0.12694
    a5 = 1.0345e-5
    a7 = 0.030234
    a9 = 3.5256e-5
    a11 = 0.25645
    a13 = 1.3514e-5
    a15 = -5.2375e-6
    a17 = 7.0e-10
    a19 = -1.6638e-6
    a21 = 0.0005
    a2 = -3.8112
    a4 = -2.0567e-3
    a6 = -6.8306
    a8 = -1.28134e-3
    a10 = -2.266e-7
    a12 = -3.4604e-3
    a14 = -28.106
    a16 = -6.3e-8
    a18 = 3.4054e-4
    a20 = -2.8673

    x1 = x[1]
    x2 = x[2]
    y1 = x1*x2
    y2 = y1*x1
    y3 = x2^2
    y4 = x1^2

    # --- function value ---

    f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
        a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
        a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
        a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
        a19*y1*y3 + a20*exp(a21*y1)
    fail = false

    # --- constraints ---

    g[1] = 1 - y1/700.0
    g[2] = y4/25.0^2 - x2/5.0
    g[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

    return f, fail

end

x0 = [10.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)
rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]
options = Dict(
    "Derivative option" => 0
)

xopt, fopt, info, out = snopta(barnes, x0, lx, ux, lg, ug, rows, cols, options)
@test isapprox(xopt[1], 49.5263; atol=1e-4)
@test isapprox(xopt[2], 19.6228; atol=1e-4)
@test isapprox(fopt, -31.6368; atol=1e-2)
@test info == "Finished successfully: optimality conditions satisfied"

end  # barnes test set

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

@testset "barnesgrad" begin

function barnesgrad(g, df, dg, x, deriv)

    a1 = 75.196
    a3 = 0.12694
    a5 = 1.0345e-5
    a7 = 0.030234
    a9 = 3.5256e-5
    a11 = 0.25645
    a13 = 1.3514e-5
    a15 = -5.2375e-6
    a17 = 7.0e-10
    a19 = -1.6638e-6
    a21 = 0.0005
    a2 = -3.8112
    a4 = -2.0567e-3
    a6 = -6.8306
    a8 = -1.28134e-3
    a10 = -2.266e-7
    a12 = -3.4604e-3
    a14 = -28.106
    a16 = -6.3e-8
    a18 = 3.4054e-4
    a20 = -2.8673

    x1 = x[1]
    x2 = x[2]
    y1 = x1*x2
    y2 = y1*x1
    y3 = x2^2
    y4 = x1^2

    # --- function value ---

    f = a1 + a2*x1 + a3*y4 + a4*y4*x1 + a5*y4^2 +
        a6*x2 + a7*y1 + a8*x1*y1 + a9*y1*y4 + a10*y2*y4 +
        a11*y3 + a12*x2*y3 + a13*y3^2 + a14/(x2+1) +
        a15*y3*y4 + a16*y1*y4*x2 + a17*y1*y3*y4 + a18*x1*y3 +
        a19*y1*y3 + a20*exp(a21*y1)
    fail = false

    # --- constraints ---

    g[1] = 1 - y1/700.0
    g[2] = y4/25.0^2 - x2/5.0
    g[3] = (x1/500.0- 0.11) - (x2/50.0-1)^2

    if deriv
        # --- derivatives of f ---

        dy1 = x2
        dy2 = y1 + x1*dy1
        dy4 = 2*x1
        dfdx1 = a2 + a3*dy4 + a4*y4 + a4*x1*dy4 + a5*2*y4*dy4 +
            a7*dy1 + a8*y1 + a8*x1*dy1 + a9*y1*dy4 + a9*y4*dy1 + a10*y2*dy4 + a10*y4*dy2 +
            a15*y3*dy4 + a16*x2*y1*dy4 + a16*x2*y4*dy1 + a17*y3*y1*dy4 + a17*y3*y4*dy1 + a18*y3 +
            a19*y3*dy1 + a20*exp(a21*y1)*a21*dy1

        dy1 = x1
        dy2 = x1*dy1
        dy3 = 2*x2
        dfdx2 = a6 + a7*dy1 + a8*x1*dy1 + a9*y4*dy1 + a10*y4*dy2 +
            a11*dy3 + a12*x2*dy3 + a12*y3 + a13*2*y3*dy3 + a14*-1/(x2+1)^2 +
            a15*y4*dy3 + a16*y4*y1 + a16*y4*x2*dy1 + a17*y4*y1*dy3 + a17*y4*y3*dy1 + a18*x1*dy3 +
            a19*y3*dy1 + a19*y1*dy3 + a20*exp(a21*y1)*a21*dy1

        df[1] = dfdx1
        df[2] = dfdx2

        # --- derivatives of g ---

        dg[1] = -x2/700.0  # 1, 1
        dg[2] = 2*x1/25^2   # 2, 1
        dg[3] = 1.0/500  # 3, 1
        dg[4] = -x1/700.0  # 1, 2
        dg[5] = -1.0/5   # 2, 2
        dg[6] = -2*(x2/50.0-1)/50.0  # 3, 2
    end

    return f, fail
end

x0 = [10.0; 10.0]
lx = [0.0; 0.0]
ux = [65.0; 70.0]
lg = -Inf*ones(3)
ug = zeros(3)
rows = [1, 2, 3, 1, 2, 3]
cols = [1, 1, 1, 2, 2, 2]
options = Dict(
    "Verify level" => 1
)

xopt, fopt, info, out = snopta(barnesgrad, x0, lx, ux, lg, ug, rows, cols, options)

@test isapprox(xopt[1], 49.5263; atol=1e-4)
@test isapprox(xopt[2], 19.6228; atol=1e-4)
@test isapprox(fopt, -31.6368; atol=1e-4)
@test info == "Finished successfully: optimality conditions satisfied"

end # barnesgrad test set

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

@testset "sparsederiv" begin

function sparsegrad(g, df, dg, x, deriv)

    f = x[1]^2 - x[2]
    fail = false

    # --- constraints ---

    g[1] = x[2] - 2*x[1]
    g[2] = -x[2]

    # --- derivatives of f ---

    df[1] = 2*x[1]
    df[2] = -1.0

    # --- derivatives of g ---

    dg[1] = -2.0  # 1, 1
    dg[2] = 1.0  # 1, 2
    dg[3] = -1.0   # 2, 2

    return f, fail

end

x0 = [0.0; 0.0]
lx = [-10.0, -10.0]
ux = [10.0, 10.0]
lg = -Inf*ones(2)
ug = zeros(2)
rows = [1, 1, 2]
cols = [1, 2, 2]
options = Dict(
    "Verify level" => 1
)

xopt, fopt, info, out = snopta(sparsegrad, x0, lx, ux, lg, ug, rows, cols, options)

@test isapprox(xopt[1], 1.0; atol=1e-6)
@test isapprox(xopt[2], 2.0; atol=1e-6)
@test isapprox(fopt, -1.0; atol=1e-6)
@test info == "Finished successfully: optimality conditions satisfied"

function sparsegrad2(g, df, dg, x, deriv)

    f = x[1]^2 - x[2]
    fail = false

    # --- constraints ---

    g[1] = x[2] - 2*x[1]
    g[2] = 3 - x[2]

    # --- derivatives of f ---

    df[1] = 2*x[1]
    df[2] = -1.0

    # --- derivatives of g ---

    dg[1] = -2.0
    dg[2] = 1.0
    dg[3] = -1.0

    return f, fail

end

xopt, fopt, info, out = snopta(sparsegrad2, x0, lx, ux, lg, ug, rows, cols, options)

@test isapprox(xopt[1], 1.5; atol=1e-6)
@test isapprox(xopt[2], 3.0; atol=1e-6)
@test isapprox(fopt, -0.75; atol=1e-6)
@test info == "Finished successfully: optimality conditions satisfied"

end # sparse test set

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

@testset "matyas" begin
    function matyas(g, df, dg, x, deriv)
        f = 0.26 * (x[1]^2 + x[2]^2) - 0.48 * x[1] * x[2]
        fail = false
        return f, fail
    end
    x0 = [5.0; 7]
    lx = [-10.0; -10]
    ux = [10.0; 10]
    lg = []
    ug = []
    rows = []
    cols = []
    options = Dict(
        "Derivative option" => 0
    )
    xopt, fopt, info, out = snopta(matyas, x0, lx, ux, lg, ug, rows, cols, options)
    @test isapprox(xopt[1], 0.0; atol=1e-4)
    @test isapprox(xopt[2], 0.0; atol=1e-4)
    @test isapprox(fopt, 0.0; atol=1e-3)
    @test info == "Finished successfully: optimality conditions satisfied"
end

ext: jl


text:
"""
    dna2rna(s::AbstractString)
Given: A DNA string t

having length at most 1000 nt.

Return: The transcribed RNA string of t
"""
function dna2rna(s::AbstractString)
	return replace(s, 'T' => 'U')
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


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

"""
    Scale(; low=0.25, high=0.75)

The scale transform of `x` is the value `(x .- xl) ./ (xh .- xl))`
where `xl = quantile(x, low)` and `xh = quantile(x, high)`.
"""
struct Scale <: Colwise
  low::Float64
  high::Float64

  function Scale(low, high)
    @assert 0 ≤ low ≤ high ≤ 1 "invalid quantiles"
    new(low, high)
  end
end

Scale(; low=0.25, high=0.75) = Scale(low, high)

assertions(::Type{Scale}) = [assert_continuous]

isrevertible(::Type{Scale}) = true

function colcache(transform::Scale, x)
  levels = (transform.low, transform.high)
  xl, xh = quantile(x, levels)
  xl == xh && ((xl, xh) = (zero(xl), one(xh)))
  (xl=xl, xh=xh)
end

colapply(::Scale, x, c)  = @. (x - c.xl) / (c.xh - c.xl)

colrevert(::Scale, y, c) = @. (c.xh - c.xl) * y + c.xl

"""
    MinMax()

The transform that is equivalent to `Scale(low=0, high=1)`.
"""
MinMax() = Scale(low=0.0, high=1.0)

"""
    Interquartile()

The transform that is equivalent to `Scale(low=0.25, high=0.75)`.
"""
Interquartile() = Scale(low=0.25, high=0.75)


text:
# Flags used to configure the FPGA manager when writing an image
# These are from include/linux/fpga/fpga-mgr.h
@enum FPGAFlags begin
    # Reconfigure the entire FPGA
    FPGA_FULL_RECONFIG        = 0

    # Do partial reconfiguration of the FPGA, if supported
    FPGA_PARTIAL_RECONFIG     = 1

    # FPGA has been configured prior to Linux booting
    FPGA_EXTERNAL_CONFIG      = 2

    # The bitstream is encrypted
    FPGA_ENCRYPTED_BITSTREAM  = 4

    # The SPI bitstream bit order is LSB first
    FPGA_BITSTREAM_LSB_FIRST  = 8

    # The FPGA bitstream is compressed
    FPGA_COMPRESSED_BITSTREAM = 16
end

@enum FPGAStates begin
    # Can't determine the state
    FPGA_STATE_UNKNOWN

    # The FPGA is powered off
    FPGA_STATE_POWER_OFF

    # The FPGA is powered up
    FPGA_STATE_POWER_UP

    # The FPGA is in the reset state
    FPGA_STATE_RESET

    # There is a firmware request in progress
    FPGA_STATE_FIRMWARE_REQ

    # There was an error with the firmware request
    FPGA_STATE_FIRMWARE_REQ_ERR

    # Preparing the FPGA for prgramming
    FPGA_STATE_WRITE_INIT

    # There was an error preparing the FPGA for programming
    FPGA_STATE_WRITE_INIT_ERR

    # Writing the bitstream image to the FPGA
    FPGA_STATE_WRITE

    # There was an error writing the bitstream image to the FPGAStates
    FPGA_STATE_WRITE_ERR

    # The FPGA is doing post-writing operations
    FPGA_STATE_WRITE_COMPLETE

    # There was an error while performing post-writing tasks
    FPGA_STATE_WRITE_COMPLETE_ERR

    # The FPGA is programmed and operating
    FPGA_STATE_OPERATING
end

# These are the strings that will be returned by the `state` sys file
const StateFromStrings = Dict{String, FPGAStates}(
    "unknown" => FPGA_STATE_UNKNOWN,
    "power off" => FPGA_STATE_POWER_OFF,
    "power up" => FPGA_STATE_POWER_UP,
    "reset" => FPGA_STATE_RESET,
    "firmware request" => FPGA_STATE_FIRMWARE_REQ,
    "firmware request error" => FPGA_STATE_FIRMWARE_REQ_ERR,
    "write init" => FPGA_STATE_WRITE_INIT,
    "write init error" => FPGA_STATE_WRITE_INIT_ERR,
    "write" => FPGA_STATE_WRITE,
    "write error" => FPGA_STATE_WRITE_ERR,
    "write complete" => FPGA_STATE_WRITE_COMPLETE,
    "write complete error" => FPGA_STATE_WRITE_COMPLETE_ERR,
    "operating" => FPGA_STATE_OPERATING,
    )

# Errors that could happen during full/partial reconfiguration
# Technically, more than one could happen
@enum FPGAStatus begin
    FPGA_STATUS_OPERATION_ERR
    FPGA_STATUS_CRC_ERR
    FPGA_STATUS_INCOMPATIBLE_IMAGE_ERR
    FPGA_STATUS_IP_PROTOCOL_ERR
    FPGA_STATUS_FIFO_OVERFLOW_ERR
end

# These are the strings that will be returned by the `status` sys file
const StatusFromStrings = Dict{String, FPGAStatus}(
    "reconfig operation error" => FPGA_STATUS_OPERATION_ERR,
    "reconfig CRC error" => FPGA_STATUS_CRC_ERR,
    "reconfig incompatible image" => FPGA_STATUS_INCOMPATIBLE_IMAGE_ERR,
    "reconfig IP protocol error" => FPGA_STATUS_IP_PROTOCOL_ERR,
    "reconfig fifo overflow error" => FPGA_STATUS_FIFO_OVERFLOW_ERR
    )


text:
include("utils.jl")

function run(fd)
    ntests = 3
    nincr_small = 1000
    nincr_big = 10000000

    test_sleep(fd, ntests)
    test_for_sleep(fd, ntests)

    test_add(fd, ntests, nincr_small)
    test_add(fd, ntests, nincr_big)

    test_add_if(fd, ntests, nincr_small)
    test_add_if(fd, ntests, nincr_big)

    test_subtract(fd, ntests, nincr_small)
    test_subtract(fd, ntests, nincr_big)

    test_subtract_if(fd, ntests, nincr_small)
    test_subtract_if(fd, ntests, nincr_big)

    test_multiply(fd, ntests, nincr_small)
    test_multiply(fd, ntests, nincr_big)

    test_multiply_if(fd, ntests, nincr_small)
    test_multiply_if(fd, ntests, nincr_big)

    test_divide(fd, ntests, nincr_small)
    test_divide(fd, ntests, nincr_big)

    test_divide_if(fd, ntests, nincr_small)
    test_divide_if(fd, ntests, nincr_big)
end

file_name = "/results/basic_operations/julia_basic_operations.txt"
path = string(chop_suffix_from_path(pwd(), "/utils/src"), file_name)
open(path, "w") do fd
    run(fd)
end

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
module Delaunay

export DelaunayNet, tri_points

using PyCall

const delaunay2D = PyNULL()

function __init__()
	if haskey(ENV, "Delaunay2D")
		pushfirst!(PyVector(pyimport("sys")["path"]), ENV["Delaunay2D"])
		copy!(delaunay2D, pyimport("delaunay2D"))
	else
		throw("Need to set ENV[\"Delaunay2D\"] to the path where the delaunay2D folder can be found. For me that is: \"C:\\Users\\matthew.heath\\repos\\Delaunay.jl\\src\\\" and reload or call Delaunay.__init__() ")
	end
end

struct DelaunayNet
	xs
	ys
	net
	function DelaunayNet(xs, ys)
		if delaunay2D == PyNULL()
			__init__()
		end
		new(xs, ys, delaunay2d(xs, ys))
	end
end

function delaunay2d(xs, ys)
	dt = delaunay2D.Delaunay2D()
	for i in 1:length(xs)
		dt.addPoint((xs[i], ys[i]))
	end
	map(t->(t[1]+1, t[2]+1, t[3]+1), dt.exportTriangles())
end

function tri_points(tri, xs, ys)
	map(i->xs[tri[i]], 1:3), map(i->ys[tri[i]], 1:3)
end

#####
end


text:
using Lint
using Compat
using Test

messageset(msgs) = Set(x.code for x in msgs)

include("variables.jl")
include("expression_iterator.jl")
include("exprutils.jl")
include("statictype.jl")

@testset "Lint Messages" begin
    include("messages.jl")
end

@testset "AST Linting" begin
	include("basics.jl")
	include("array.jl")
	include("badvars.jl")
	include("bitopbool.jl")
	include("comprehensions.jl")
	include("curly.jl")
	include("deadcode.jl")
	include("deprecate.jl")
	include("dictkey.jl")
	include("doc.jl")
	include("dupexport.jl")
	include("forloop.jl")
	include("funcall.jl")
	include("globals.jl")
	include("ifstmt.jl")
	include("import.jl")
	include("lambda.jl")
	include("macro.jl")
	include("meta.jl")
	include("pragma.jl")
	include("range.jl")
	include("ref.jl")
	include("similarity.jl")
	include("strings.jl")
	include("style.jl")
	include("I340.jl")
	include("I341.jl")
	include("I343.jl")
	include("I481.jl")
	include("W361.jl")
	include("E100.jl")
	include("throw.jl")
	include("tuple.jl")
	include("type.jl")
	include("typecheck.jl")
	include("undeclare.jl")
	include("using.jl")
	include("versions.jl")
	include("stagedfuncs.jl")
	include("incomplete.jl")
	include("misuse.jl")
end

@testset "Lint File" begin
    path = joinpath(@__DIR__, "DEMOFILE.jl")
    @test !isempty(lintfile(path))
end

try include("linthelper.jl") end
include("bugs.jl")
include("lintself.jl")

try
@testset "Server" begin
    include("server.jl")
end
end


text:
# Functionality for computing a suitable sampling parameter of an approximation problem

function samplingparameter(ss::SamplingStyle, platform::Platform, plt_par;
        verbose=false, L = -1, smpl_par = -1, options...)
    # Was the parameter specified as an option?
    if L != -1
        smpl_par = L
    elseif smpl_par != -1
        smpl_par
    else
        # It wasn't. Try to deduce its value from the options given.
        smpl_par = deduce_samplingparameter(ss, platform, plt_par; verbose=verbose, options...)
    end
    verbose && println("Sampling parameter: using smpl_par = $smpl_par")
    return smpl_par
end

samplingparameter(ss::GenericSamplingStyle, platform::Platform, plt_par; options...) = nothing

samplingparameter(ss::WeightedDiscreteStyle, platform::Platform, plt_par; options...) =
    samplingparameter(unweighted(ss), platform, plt_par; options...)

samplingparameter(ss::DiscreteGramStyle, platform::Platform, plt_par; options...) =
    samplingparameter(OversamplingStyle(), platform, plt_par; options...)

deduce_samplingparameter(ss::SamplingStyle, platform::Platform, plt_par; options...) =
    deduce_samplingparameter(ss, dictionary(platform, plt_par); options...)

deduce_samplingparameter(ss::OversamplingStyle, platform::Platform, plt_par; options...) =
    deduce_oversampling_parameter(ss, platform, plt_par; options...)

function deduce_oversampling_parameter(ss::OversamplingStyle, dict::Dictionary;
        verbose = false, oversamplingfactor = 2, M = -1, options...)
    # In the absence of smpl_par, we deduce M and then find the best matching smpl_par
    # M is either supplied, or we compute it based on the (default) oversamplingfactor
    if M == -1
        M = round(Int, oversamplingfactor * length(dict))
    end
    verbose && println("Sampling parameter: oversamplingfactor=$oversamplingfactor, length=$(length(dict)), M=$M")
    smpl_par = match_and_correct_sampling_parameter(dict, M; samplingstyle=ss, options...)
    verbose && println("Sampling parameter: best match for M = $M is smpl_par = $smpl_par")
    smpl_par
end

function deduce_oversampling_parameter(ss::OversamplingStyle, platform::Platform, plt_par;
        verbose = false, oversamplingfactor = 2, M = -1, options...)
    dict = dictionary(platform, plt_par)
    if M == -1
        M = round(Int, oversamplingfactor * length(dict))
    end
    verbose && println("Sampling parameter: oversamplingfactor=$oversamplingfactor, length=$(length(dict)), M=$M")
    smpl_par = match_and_correct_sampling_parameter(platform, plt_par, M; samplingstyle=ss, options...)
    verbose && println("Sampling parameter: best match for M = $M is smpl_par = $smpl_par")
    smpl_par
end

deduce_samplingparameter(ss::InterpolationStyle, dict::Dictionary; options...) =
    match_and_correct_sampling_parameter(dict, length(dict); samplingstyle=ss, options...)

# TODO: implement this one better (more general)
deduce_samplingparameter(::GramStyle, dict::Dictionary; options...) = length(dict)

deduce_samplingparameter(::RectangularGramStyle, dict::Dictionary; projectiondict, options...) =
    length(projectiondict)

deduce_samplingparameter(::GridStyle, dict::Dictionary; options...) = platformparameter(dict)

## Product sampling

oversamplingfactor_to_tuple(of::Real, dim::Int) = ntuple(t->of^(1/dim), Val(dim))
oversamplingfactor_to_tuple(of::Tuple, dim::Int) = (@assert length(of)==dim; of)

deduce_samplingparameter(ss::ProductSamplingStyle, platform::ProductPlatform, plt_par; options...) =
    map((style,pl,ppar)->samplingparameter(style, pl, ppar; options...),
        factors(ss), factors(platform), productparameter(platform, plt_par))

function deduce_samplingparameter(ss::ProductSamplingStyle{NTuple{N,OversamplingStyle}},
        platform::ProductPlatform, plt_par;
        oversamplingfactor=2, M = -1, verbose=false, options...) where {N}
    dict = dictionary(platform, plt_par)
    if !(M == -1)
        if M isa Integer
            oversamplingfactor = M / length(dict)
        else
            @assert length(M) == N
            oversamplingfactor = M ./ map(length, factors(dict))
        end
    end
    oversamplingfactor = oversamplingfactor_to_tuple(oversamplingfactor, N)
    smpl_par = tuple(map((pl,ppar,of)->deduce_samplingparameter(OversamplingStyle(), pl, ppar;
        oversamplingfactor=of, verbose=verbose, options...),
            factors(platform), productparameter(platform, plt_par), oversamplingfactor)...)
    verbose && println("Sampling parameter: best match for M = $M is smpl_par = $smpl_par")
    smpl_par
end


text:
function aggregate(obj::AggregatableStatistic, arg0::Iterable)
    return jcall(obj, "aggregate", void, (Iterable,), arg0)
end

function aggregate(obj::AggregatableStatistic, arg0::Vector{Object})
    return jcall(obj, "aggregate", void, (Vector{Object},), arg0)
end

function aggregate(obj::FourthMoment, arg0::SecondMoment)
    return jcall(obj, "aggregate", void, (SecondMoment,), arg0)
end

function clear(obj::FourthMoment)
    return jcall(obj, "clear", void, ())
end

function copy(obj::FourthMoment)
    return jcall(obj, "copy", FourthMoment, ())
end

function get_result(obj::FourthMoment)
    return jcall(obj, "getResult", jdouble, ())
end

function increment(obj::FourthMoment, arg0::jdouble)
    return jcall(obj, "increment", void, (jdouble,), arg0)
end


text:
struct NewtonSolver{T, AT, FT, TJ, TL, TS <: LinesearchState} <: AbstractNewtonSolver{T,AT}
    @newton_solver_variables

    function NewtonSolver{T,AT,FT,TJ,TL,TS}(x, y, F!, Jparams, linear_solver, linesearch, cache, config) where {T,AT,FT,TJ,TL,TS}
        status = NonlinearSolverStatus{T}(length(x))

        new(x, y, F!, Jparams, linear_solver, linesearch, cache, config, status)
    end
end

function NewtonSolver(x::AbstractVector{T}, y::AbstractVector{T}, F!; J! = nothing, linesearch = Static(), config = Options()) where {T}
    n = length(y)
    Jparams = JacobianParameters{T}(J!, F!, n)
    cache = NewtonSolverCache(x, y)
    linear_solver = LinearSolver(y)
    ls = LinesearchState(linesearch, linesearch_objective(F!, Jparams, cache))
    NewtonSolver{T, typeof(x), typeof(F!), typeof(Jparams), typeof(linear_solver), typeof(ls)}(x, y, F!, Jparams, linear_solver, ls, cache, config)
end

function solver_step!(s::NewtonSolver{T}) where {T}
    # shortcuts
    x = s.x
    y = s.y
    δ = s.cache.δx

    # compute Jacobian
    compute_jacobian!(s)

    # factorize linear solver
    factorize!(s.linear, s.cache.J)

    # copy previous solution
    s.cache.x₀ .= x

    # solve J δx = -f(x)
    rmul!(y, -1)
    ldiv!(δ, s.linear, y)

    # apply line search
    α = s.linesearch()
    x .+= α .* δ

    # compute residual
    s.F!(y, x)
    residual!(status(s), x, y)
end


text:
using JLD
using LinearAlgebra

# Setting number of threads to one
LinearAlgebra.BLAS.set_num_threads(1)

include("./SortSparse.jl")
include("./KoLesky.jl")
include("./CovFuncs.jl")

N = 1000000
x = rand(2,N)

ρList = [1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0]
timings = zeros(size(ρList))
nonz = zeros(size(ρList))
timingsAgg = zeros(size(ρList))
nonzAgg = zeros(size(ρList))
timingsBessel = zeros(size(ρList))
nonzBessel = zeros(size(ρList))
timingsBesselAgg = zeros(size(ρList))
nonzBesselAgg = zeros(size(ρList))

# Timings
for (ρInd, ρ) in enumerate(ρList)
  global timings
  global nonz

  covfunc!(r) = matern32!(r, 0.1)
  colptr, rowval, P, revP, distances =  sortSparseRev(x, ρ, 1)
  xLoc = x[:, P]
  skeletons = construct_skeletons(colptr, rowval, distances, 1.0)
  # Making sure everything is compiled
  if ρInd == 1
    @elapsed assembleL(skeletons, xLoc, covfunc!)
  end
  timings[ρInd] = @elapsed L = assembleL(skeletons, xLoc, covfunc!)
  nonz[ρInd] = nnz(L)
end

# Timings with aggregation
for (ρInd, ρ) in enumerate(ρList)
  global timingsAgg
  global nonzAgg
  covfunc!(r) = matern32!(r, 0.1)
  colptr, rowval, P, revP, distances =  sortSparseRev(x, ρ, 1)
  xLoc = x[:, P]
  skeletons = construct_skeletons(colptr, rowval, distances, 1.5)
  # Making sure everything is compiled
  if ρInd == 1
    @elapsed assembleL(skeletons, xLoc, covfunc!)
  end
  timingsAgg[ρInd] = @elapsed L = assembleL(skeletons, xLoc, covfunc!)
  nonzAgg[ρInd] = nnz(L)
end

# Timings with Bessel function
for (ρInd, ρ) in enumerate(ρList)
  global timingsBessel
  global nonzBessel
  function covfunc!(r)
    r .= matern.(r, 0.1, 1.5)
  end
  colptr, rowval, P, revP, distances =  sortSparseRev(x, ρ, 1)
  xLoc = x[:, P]
  skeletons = construct_skeletons(colptr, rowval, distances, 1.0)
  # Making sure everything is compiled
  if ρInd == 1
    @elapsed assembleL(skeletons, xLoc, covfunc!)
  end
  timingsBessel[ρInd] = @elapsed L = assembleL(skeletons, xLoc, covfunc!)
  nonzBessel[ρInd] = nnz(L)
end

# Timings with Bessel function and aggregation
for (ρInd, ρ) in enumerate(ρList)
  global timingsBesselAgg
  global nonzBesselAgg
  function covfunc!(r)
    r .= matern.(r, 0.1, 1.5)
  end
  colptr, rowval, P, revP, distances =  sortSparseRev(x, ρ, 1)
  xLoc = x[:, P]
  skeletons = construct_skeletons(colptr, rowval, distances, 1.5)
  # Making sure everything is compiled
  if ρInd == 1
    @elapsed assembleL(skeletons, xLoc, covfunc!)
  end
  timingsBesselAgg[ρInd] = @elapsed L = assembleL(skeletons, xLoc, covfunc!)
  nonzBesselAgg[ρInd] = nnz(L)
end

save("./out/jld/TimingVaryRho.jld",
     "ρList", ρList,
     "timings", timings,
     "nonz", nonz,
     "timingsAgg", timingsAgg,
     "nonzAgg", nonzAgg,
     "timingsBessel", timingsBessel,
     "nonzBessel", nonzBessel,
     "timingsBesselAgg", timingsBesselAgg,
     "nonzBesselAgg", nonzBesselAgg)


text:
#qt.h
@getf libqt mgl_create_graph_qt(draw::DrawCallback, title::Ptr{Cchar}, par::Ptr{Void}, load::LoadCallback)::HMGL
@getf libqt mgl_qt_run()::Cint


text:
"""
$(TYPEDEF)

Fluid property for scalar conservation laws

# Fields

$(FIELDS)
"""
@kwdef mutable struct Scalar{TA,TB} <: AbstractProperty
    a::TA = 1.0
    μᵣ::TB = 1e-8
end

"""
$(TYPEDEF)

Radiation property for linear Boltzmann equation

# Fields

$(FIELDS)
"""
@kwdef mutable struct Radiation{T1,T2,T3,T4} <: AbstractProperty
    Kn::T1 = 1.0
    σs::T2 = 1.0
    σa::T2 = 0.0
    m::T3 = 1e-3
    np::T4 = 1000
end

function Radiation(_Kn::Union{Real,AV}, _ss::Union{Real,AV}, _sa::Union{Real,AV})
    Kn = deepcopy(_Kn)
    σs = deepcopy(_ss)
    σa = deepcopy(_sa)
    m = 1e-3
    np = 1000

    return Radiation{typeof(Kn),typeof(σs),typeof(m),typeof(np)}(Kn, σs, σa, m, np)
end

"""
$(TYPEDEF)

Gas property

# Fields

$(FIELDS)
"""
@kwdef mutable struct Gas{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11<:Integer,T12} <: AbstractGas
    Kn::T1 = 1e-2
    Ma::T2 = 0.0
    Pr::T3 = 1.0
    K::T4 = 2.0
    γ::T5 = 5 / 3
    ω::T6 = 0.81
    αᵣ::T7 = 1.0
    ωᵣ::T8 = 0.5
    μᵣ::T9 = ref_vhs_vis(Kn, αᵣ, ωᵣ)
    m::T10 = 1e-3
    np::T11 = 1000
    fsm::T12 = nothing
end

"""
$(TYPEDEF)

Diatomic gas property

# Fields

$(FIELDS)
"""
@kwdef struct DiatomicGas{TA,TI,TF} <: AbstractGas
    Kn::TA = 1e-2
    Ma::TA = 0.0
    Pr::TA = 1.0
    K::TI = 2.0
    Kr::TI = 2.0
    γ::TF = 7 / 5
    ω::TF = 0.74
    αᵣ::TF = 1.0
    ωᵣ::TF = 0.5
    μᵣ::TF = ref_vhs_vis(Kn, αᵣ, ωᵣ)
    T₀::TF = 89.1 / 273
    Z₀::TF = 18.1
    σ::TF = 1 / 1.55
    ω₁::TF = 0.2354
    ω₂::TF = 0.3049
end

"""
$(TYPEDEF)

Multi-component gas property

# Fields

$(FIELDS)
"""
@kwdef struct Mixture{A,B,C,D,E,F,G,H,I} <: AbstractGas
    Kn::A = 1e-2
    Ma::B = 0.0
    Pr::C = 1.0
    K::D = 2.0
    γ::E = 5 / 3
    mi::F = 1.0
    ni::G = 0.5
    me::H = 0.5
    ne::I = 0.5
end

"""
$(TYPEDEF)

1D plasma property

# Fields

$(FIELDS)
"""
@kwdef struct Plasma1D{A,B,C,TD,E,F,G,H,I,J,TK,L,M,N,O,P} <: AbstractPlasma
    Kn::A = 1e-2
    Ma::B = 0.0
    Pr::C = 1.0
    K::TD = 2.0
    γ::E = 5 / 3

    mi::F = 1.0
    ni::G = 0.5
    me::H = 0.5
    ne::I = 0.5
    lD::J = 0.01
    rL::TK = 0.01

    sol::L = 100
    χ::M = 1
    ν::N = 1
    Ap::O
    An::O
    D::P
end

# unified consideration of deterministic and stochastic conditions
# {mi, ni, me, ne, sol, χ, ν} keep the same as before
function Plasma1D(
    Kn::Array,
    Ma::Union{Real,AA},
    Pr::Union{Real,AA},
    K::Union{Real,AA},
    γ::Union{Real,AA},
    mi::Real,
    ni::Real,
    me::Real,
    ne::Real,
    lD::Union{Real,AA},
    rL::Union{Real,AA},
    sol::Real,
    χ::Real,
    ν::Real,
)
    # A^+
    Ap = Array{Float64}(undef, 8, 8)
    Ap[1, 1] = (sol * χ) / 2.0
    Ap[7, 1] = χ / 2.0
    Ap[2, 2] = sol / 2.0
    Ap[6, 2] = 0.5
    Ap[3, 3] = sol / 2.0
    Ap[5, 3] = -1.0 / 2.0
    Ap[4, 4] = (sol * ν) / 2.0
    Ap[8, 4] = (sol^2 * ν) / 2.0
    Ap[3, 5] = -sol^2 / 2.0
    Ap[5, 5] = sol / 2.0
    Ap[2, 6] = sol^2 / 2.0
    Ap[6, 6] = sol / 2.0
    Ap[1, 7] = (sol^2 * χ) / 2.0
    Ap[7, 7] = (sol * χ) / 2.0
    Ap[4, 8] = ν / 2.0
    Ap[8, 8] = (sol * ν) / 2.0

    # A^-
    An = Array{Float64}(undef, 8, 8)
    An[1, 1] = -(sol * χ) / 2.0
    An[7, 1] = χ / 2.0
    An[2, 2] = -sol / 2.0
    An[6, 2] = 1.0 / 2.0
    An[3, 3] = -sol / 2.0
    An[5, 3] = -1.0 / 2.0
    An[4, 4] = -(sol * ν) / 2.0
    An[8, 4] = (sol^2 * ν) / 2.0
    An[3, 5] = -sol^2 / 2.0
    An[5, 5] = -sol / 2.0
    An[2, 6] = sol^2 / 2.0
    An[6, 6] = -sol / 2.0
    An[1, 7] = (sol^2 * χ) / 2.0
    An[7, 7] = -(sol * χ) / 2.0
    An[4, 8] = ν / 2.0
    An[8, 8] = -(sol * ν) / 2.0

    # eigenvalues of A
    D = [sol, sol, sol * χ, sol * ν, -sol, -sol, -sol * χ, -sol * ν]

    # inner constructor method
    return Plasma1D{
        typeof(Kn),
        typeof(Ma),
        typeof(Pr),
        typeof(K),
        typeof(γ),
        typeof(mi),
        typeof(ni),
        typeof(me),
        typeof(ne),
        typeof(lD),
        typeof(rL),
        typeof(sol),
        typeof(χ),
        typeof(ν),
        typeof(Ap),
        typeof(D),
    }(
        Kn,
        Ma,
        Pr,
        K,
        γ,
        mi,
        ni,
        me,
        ne,
        lD,
        rL,
        sol,
        χ,
        ν,
        Ap,
        An,
        D,
    )
end

"""
$(TYPEDEF)

2D plasma property

# Fields

$(FIELDS)
"""
@kwdef struct Plasma2D{A,B,C,D,E,F,G,H,I,J,TK,L,M,N,O,P} <: AbstractPlasma
    Kn::A = 1e-2
    Ma::B = 0.0
    Pr::C = 1.0
    K::D = 2
    γ::E = 5 / 3

    mi::F = 1.0
    ni::G = 0.5
    me::H = 0.5
    ne::I = 0.5
    lD::J = 0.01
    rL::TK = 0.01

    sol::L = 100
    χ::M = 1
    ν::N = 1
    A1p::O
    A1n::O
    A2p::O
    A2n::O
    D1::P
    D2::P
end

# unified consideration of deterministic and stochastic conditions
# {mi, ni, me, ne, sol, χ, ν} keep the same as before
function Plasma2D(
    Kn::Array,
    Ma::Union{Real,AA},
    Pr::Union{Real,AA},
    K::Union{Real,AA},
    γ::Union{Real,AA},
    mi::Real,
    ni::Real,
    me::Real,
    ne::Real,
    lD::Union{Real,AA},
    rL::Union{Real,AA},
    sol::Real,
    χ::Real,
    ν::Real,
)
    # A₁+
    A1p = zeros(8, 8)
    A1p[1, 1] = (sol * χ) / 2.0
    A1p[7, 1] = χ / 2.0
    A1p[2, 2] = sol / 2.0
    A1p[6, 2] = 0.5
    A1p[3, 3] = sol / 2.0
    A1p[5, 3] = -1.0 / 2.0
    A1p[4, 4] = (sol * ν) / 2.0
    A1p[8, 4] = (sol^2 * ν) / 2.0
    A1p[3, 5] = -sol^2 / 2.0
    A1p[5, 5] = sol / 2.0
    A1p[2, 6] = sol^2 / 2.0
    A1p[6, 6] = sol / 2.0
    A1p[1, 7] = (sol^2 * χ) / 2.0
    A1p[7, 7] = (sol * χ) / 2.0
    A1p[4, 8] = ν / 2.0
    A1p[8, 8] = (sol * ν) / 2.0

    # A₁-
    A1n = zeros(8, 8)
    A1n[1, 1] = -(sol * χ) / 2.0
    A1n[7, 1] = χ / 2.0
    A1n[2, 2] = -sol / 2.0
    A1n[6, 2] = 1.0 / 2.0
    A1n[3, 3] = -sol / 2.0
    A1n[5, 3] = -1.0 / 2.0
    A1n[4, 4] = -(sol * ν) / 2.0
    A1n[8, 4] = (sol^2 * ν) / 2.0
    A1n[3, 5] = -sol^2 / 2.0
    A1n[5, 5] = -sol / 2.0
    A1n[2, 6] = sol^2 / 2.0
    A1n[6, 6] = -sol / 2.0
    A1n[1, 7] = (sol^2 * χ) / 2.0
    A1n[7, 7] = -(sol * χ) / 2.0
    A1n[4, 8] = ν / 2.0
    A1n[8, 8] = -(sol * ν) / 2.0

    D1 = [sol, sol, sol * χ, sol * ν, -sol, -sol, -sol * χ, -sol * ν]

    # A₂+
    A2p = zeros(8, 8)
    A2p[1, 1] = sol / 2.0
    A2p[6, 1] = -1.0 / 2.0
    A2p[2, 2] = (sol * χ) / 2.0
    A2p[7, 2] = χ / 2.0
    A2p[3, 3] = sol / 2.0
    A2p[4, 3] = 1.0 / 2.0
    A2p[3, 4] = (sol^2) / 2.0
    A2p[4, 4] = sol / 2.0
    A2p[5, 5] = (sol * ν) / 2.0
    A2p[8, 5] = (sol^2 * ν) / 2.0
    A2p[1, 6] = -sol^2 / 2.0
    A2p[6, 6] = sol / 2.0
    A2p[2, 7] = (sol^2 * χ) / 2.0
    A2p[7, 7] = (sol * χ) / 2.0
    A2p[5, 8] = ν / 2.0
    A2p[8, 8] = (sol * ν) / 2.0

    # A₂-
    A2n = zeros(8, 8)
    A2n[1, 1] = -sol / 2.0
    A2n[6, 1] = -1.0 / 2.0
    A2n[2, 2] = -(sol * χ) / 2.0
    A2n[7, 2] = χ / 2.0
    A2n[3, 3] = -sol / 2.0
    A2n[4, 3] = 1.0 / 2.0
    A2n[3, 4] = (sol^2) / 2.0
    A2n[4, 4] = -sol / 2.0
    A2n[5, 5] = -(sol * ν) / 2.0
    A2n[8, 5] = (sol^2 * ν) / 2.0
    A2n[1, 6] = -sol^2 / 2.0
    A2n[6, 6] = -sol / 2.0
    A2n[2, 7] = (sol^2 * χ) / 2.0
    A2n[7, 7] = -(sol * χ) / 2.0
    A2n[5, 8] = ν / 2.0
    A2n[8, 8] = -(sol * ν) / 2.0

    D2 = [sol, sol, sol * χ, sol * ν, -sol, -sol, -sol * χ, -sol * ν]

    # inner constructor method
    return Plasma2D{
        typeof(Kn),
        typeof(Ma),
        typeof(Pr),
        typeof(K),
        typeof(γ),
        typeof(mi),
        typeof(ni),
        typeof(me),
        typeof(ne),
        typeof(lD),
        typeof(rL),
        typeof(sol),
        typeof(χ),
        typeof(ν),
        typeof(A1p),
        typeof(D1),
    }(
        Kn,
        Ma,
        Pr,
        K,
        γ,
        mi,
        ni,
        me,
        ne,
        lD,
        rL,
        sol,
        χ,
        ν,
        A1p,
        A1n,
        A2p,
        A2n,
        D1,
        D2,
    )
end


text:
module Demo
    using Dates
    using Plots
    #using LaTeXStrings

    import ..System_parameters:Params_set,parameterloading
    import ..LTK_universe:Universe
    import  ..MD:construct_MD_parameters
    import ..Measurements:Measurement_set,measurements
    import ..Verbose_print:println_verbose1,println_verbose2,Verbose_1
    import ..Heatbath:heatbath!

    system = Dict()
    actions = Dict()
    md = Dict()
    cg = Dict()
    wilson = Dict()
    staggered = Dict()
    measurement = Dict()
    # - - parameters - - - - - - - - - - -
# - - parameters - - - - - - - - - - -
    system["saveU_dir"] = ""
    system["verboselevel"] = 1
    #system["L"] = (6, 6, 6, 6)
    system["L"] = (4, 4, 4, 4)
    system["Nwing"] = 1
    system["Nsteps"] = 10000
    system["quench"] = true
    system["logfile"] = "Heatbath_L06060606_beta8.0_quenched.txt"
    system["initial"] = "hot"
    system["Dirac_operator"] = nothing
    system["log_dir"] = "./logs"
    system["Nthermalization"] = 10
    system["update_method"] = "Heatbath"
    system["randomseed"] = 111
    system["NC"] = 3
    system["saveU_every"] = 1
    system["BoundaryCondition"] = [1, 1, 1, -1]
    system["saveU_format"] = nothing
    #system["β"] = 6.93015
    system["β"] = 6.75850661032668353
    actions["use_autogeneratedstaples"] = false
    actions["couplingcoeff"] = Any[]
    actions["couplinglist"] = Any[]
    md["Δτ"] = 0.05
    md["SextonWeingargten"] = false
    md["MDsteps"] = 20
    cg["eps"] = 1.0e-19
    cg["MaxCGstep"] = 3000
    wilson["Clover_coefficient"] = 0
    wilson["r"] = 1
    wilson["hop"] = 0
    staggered["Nf"] = 0
    staggered["mass"] = 0
    measurement["measurement_methods"] = Dict[Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Polyakov_loop"), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Plaquette")]
    measurement["measurement_dir"] = "Heatbath_L06060606_beta8.0_quenched"
    measurement["measurement_basedir"] = "./measurements"
    # - - - - - - - - - - - - - - - - - - -

    function demo()

        params_set = Params_set(system,actions,md,cg,wilson,staggered,measurement)

        parameters = parameterloading(params_set)
        univ = Universe(parameters)

        mdparams = construct_MD_parameters(parameters)

        meas = Measurement_set(univ,parameters.measuredir,measurement_methods=parameters.measurement_methods)
        run_demo!(parameters,univ,meas)

    end

    function run_demo!(parameters,univ,meas)
        plt1 = histogram([0],label=nothing) #plot1
        plt2 = plot([], [],label=nothing) #plot2
        ylabel!("Plaquette")
        xlabel!("MC time")
        plt3 = histogram([0],label=nothing) #plot3
        plt4 = plot([], [],label=nothing) #plot4
        ylabel!("|Polyakov loop|")
        xlabel!("MC time")
        plt5 = scatter([], [],label=nothing,title="Polyakov loop") #plot5 scatter
        ylabel!("Im")
        xlabel!("Re")
        plt6 = plot([], [],label=nothing) #plot6
        ylabel!("Arg(Polyakov loop)")
        xlabel!("MC time")
        plot(plt1, plt2, plt3, plt4, plt5, plt6, layout = 6)

        hist_plaq=[]
        hist_poly=[]
        hist_poly_θ=[]

        function plot_refresh!(plt1,plt2,plt3,plt4,plt5,plt6,hist_plaq,hist_poly,hist_poly_θ, plaq, poly, itrj)
            bins=round(Int,log(itrj)*4+1)
            if itrj < 500
                bins=round(Int,log(itrj)*1.5+1)
            elseif itrj < 1000
                bins=round(Int,log(itrj)*2.5+1)
            elseif itrj < 2000
                bins=round(Int,log(itrj)*3.5+1)
            end
            append!(hist_plaq, plaq)
            append!(hist_poly, abs(poly))
            # omit un-thermalized part
            if (10 < itrj < 1000) & (itrj%5==0)
                print("remove unthermalize part $(length(hist_plaq)) -> ")
                popfirst!(hist_plaq)
                popfirst!(hist_poly)
                println(" $(length(hist_plaq))")
            end
            #
            plt1 = histogram(hist_plaq,bins=bins,label=nothing) #plot1
            #plot!(plt1,title="SU(3), Quenched, L=6^4,")
            plot!(plt1,title="SU(3), Quenched, L=4^4,")
            xlabel!("Plaquette")
            plot!(plt2,title="Heatbath")
            plt3 = histogram(hist_poly,bins=bins,label=nothing) #plot3
            xlabel!("|Polyakov loop|")
            #
            if itrj<500
                push!(plt2, 1, itrj, plaq)
            elseif 500<itrj<1000
                if itrj%2==0
                push!(plt2, 1, itrj, plaq)
                end
            elseif 1000<itrj
                if itrj%5==0
                push!(plt2, 1, itrj, plaq)
                end
            end
            push!(plt4, 1, itrj, abs(poly))
            if itrj<500
                push!(plt5, 1, real(poly), imag(poly))
            end
            if 500<itrj<1000
                if itrj%10==0
                    push!(plt5, 1, real(poly), imag(poly))
                end
            end
            if 1000<itrj
                if itrj%50==0
                    push!(plt5, 1, real(poly), imag(poly))
                end
            end
            #
            push!(plt6, 1, itrj, angle(poly))
            #
            plot(plt1, plt2, plt3, plt4, plt5, plt6, layout = 6)
            gui()
        end

        @assert parameters.update_method == "Heatbath"
        verbose = Verbose_1()
        Nsteps = parameters.Nsteps
        numaccepts = 0
        plaq,poly = measurements(0,univ.U,univ,meas;verbose = verbose) # check consistency of preparation.

        for itrj=1:Nsteps
            @time heatbath!(univ)
            plaq,poly = measurements(itrj,univ.U,univ,meas;verbose = verbose)
            plot_refresh!(plt1,plt2,plt3,plt4,plt5,plt6,hist_plaq,hist_poly,hist_poly_θ, plaq, poly, itrj)

            println_verbose1(verbose,"-------------------------------------")
            #println("-------------------------------------")
            flush(stdout)
            flush(verbose)
        end
    end

end


text:
"""
    hopper
    	model inspired by Raibert hopper
            from "Dynamically Stable Legged Locomotion"
		q = (x, z, t, r)
			x - lateral position
			z - vertical position
			t - body orientation
			r - leg length
"""
struct Hopper{T} <: Model{T}
    # dimensions
    nq::Int # generalized coordinates
    nu::Int # controls
    nw::Int # parameters

    # parameters
    mass_body::T
    mass_foot::T
    inertia_body::T
    body_radius::T
    foot_radius::T
    leg_len_max::T
    leg_len_min::T
    friction_joint::T

    # environment
    friction_body_world::T
    friction_foot_world::T
	gravity::T
end

# kinematics
function kinematics_body(::Hopper, q)
	[q[1],
	 q[2]]
end

function kinematics_body_jacobian(::Hopper, q)
    [1.0 0.0 0.0 0.0;
	 0.0 1.0 0.0 0.0]
end

function kinematics_foot(::Hopper, q)
	[q[1] + q[4] * sin(q[3]),
	 q[2] - q[4] * cos(q[3])]
end

function kinematics_foot_jacobian(::Hopper, q)
    [1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
	 0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]))]
end

# Lagrangian
function lagrangian(model::Hopper, q, q̇)
    L = 0.0

    # body kinetic energy
    v_body = q̇[1:2]
    ω_body = q̇[3]
    L += 0.5 * model.mass_body * dot(v_body, v_body)
    L += 0.5 * model.inertia_body * ω_body^2.0

    # body potential energy
    z_body = q[2]
    L -= model.mass_body * model.gravity * z_body

    # foot kinetic energy
    v_foot = kinematics_foot_jacobian(model, q) * q̇
    L += 0.5 * model.mass_foot * dot(v_foot, v_foot)

    # foot potential energy
    z_foot =  kinematics_foot(model, q)[2]
    L -= model.mass_foot * model.gravity * z_foot

    return L
end

# signed distance
function signed_distance(model::Hopper, q)
    [
     q[2] - model.body_radius,                    # body-world distance
     kinematics_foot(model, q)[2] - model.foot_radius, # foot-world distance
     q[4] - model.leg_len_min,                    # leg length min.
     model.leg_len_max - q[4]                     # leg length max.
    ]
end

function signed_distance(model::Hopper, q, pt1, pt2)
    # error between robot kinematics and point on nearest surface
    e1 = kinematics_body(model, q) - pt1
    e2 = kinematics_foot(model, q) - pt2
    @show pt1
    @show pt2
    [
     sqrt(dot(e1, e1)) - model.body_radius,  # body-world distance
     sqrt(dot(e2, e2)) - model.foot_radius,  # foot-world distance
     q[4] - model.leg_len_min,             # leg length min.
     model.leg_len_max - q[4]              # leg length max.
    ]
end

# input matrix
function input_jacobian(::Hopper, q)
	transpose([
     0.0 -sin(q[3]);
     0.0 cos(q[3]);
     1.0 0.0;
     0.0 1.0
    ])
end

# contact Jacobian
function contact_jacobian(::Hopper, q)
    [1.0 0.0 0.0 0.0; # note that body torque due to friction is added in residual
     0.0 1.0 0.0 0.0;
     1.0 0.0 (q[4] * cos(q[3])) sin(q[3]);
     0.0 1.0 (q[4] * sin(q[3])) (-1.0 * cos(q[3]));
     0.0 0.0 0.0 1.0;
     0.0 0.0 0.0 -1.0]
end

# nominal configuratoin
function nominal_configuration(model::Hopper)
    return [0.0, 1.0, 0.0 * π, 0.5]
end

function friction_coefficients(model::Hopper)
    return [model.friction_body_world; model.friction_foot_world]
end

# dimensions
nq = 4
nu = 2
nw = 0

# parameters
mass_body = 1.0
mass_foot = 0.1
inertia_body = 0.1
body_radius = 0.1
foot_radius = 0.05
leg_len_max = 1.0
leg_len_min = 0.25
friction_joint = 0.1

# environment
friction_body_world = 0.5
friction_foot_world = 0.5
gravity = 9.81

# model
hopper = Hopper(nq, nu, nw,
            mass_body, mass_foot, inertia_body,
            body_radius, foot_radius,
            leg_len_max, leg_len_min,
            friction_joint,
            friction_body_world, friction_foot_world,
            gravity)

# contact kinematics
hopper_contact_kinematics = [q -> kinematics_body(hopper, q),
    q -> kinematics_foot(hopper, q)]

hopper_contact_kinematics_jacobians = [q -> kinematics_body_jacobian(hopper, q),
    q -> kinematics_foot_jacobian(hopper, q)]

name(::Hopper) = :hopper


text:
@testset "LaguerrePolynomial" begin
    @testset "LaguerrePolynomial" begin
        # test construction of LaguerrePolynomial
        orders = [0,2,3]
        coeffs = [1.0, 2, 0.5]
        L = LaguerrePolynomial(orders, coeffs)
        @test L.orders == orders
        @test all(isapprox.(coeffs, L.coeffs))
        # test single Laguerre Polynomial
        L = LaguerrePolynomial(2, 1 // 2)
        @test L.orders == [2]
        @test L.coeffs == [1 // 2]

        # Laguerre Polynomial L3 with coefficient 1.
        L = LaguerrePolynomial(3)
        @test L.orders == [3]
        @test L.coeffs == [1]
        # test throw
        @test_throws DimensionMismatch LaguerrePolynomial([1,2], [2,3,2])
        @test_throws DomainError LaguerrePolynomial(-2, 1)
        # test string representation
        L = LaguerrePolynomial([2,3], [1,2])
        st = Laguerre.lashow(L)
        @test st == "(1)L2 + (2)L3"
    end

    @testset "LaguerrePolynomial addition" begin
        L1 = LaguerrePolynomial(2, 1 // 2)
        L2 = LaguerrePolynomial(3, 1.0)
        L = L1 + L2
        @test typeof(L) == LaguerrePolynomial{typeof(1.0)}
        @test L.orders == [2, 3]
        @test all(isapprox.(L.coeffs, [1 / 2, 1.0]))
        #
        L1 = LaguerrePolynomial(2, 1)
        L2 = LaguerrePolynomial(2, 1 // 3)
        L = L1 + L2
        @test typeof(L) == LaguerrePolynomial{typeof(1 // 3)}
        @test L.orders == [2]
        @test L.coeffs == [1 + 1 // 3]

        # test subtraction
        L = L1 - L2
        @test L.orders == [2]
        @test L.coeffs == [1 - 1 // 3]

    end

    @testset "LaguerrePolynomial * number" begin
        orders = [0,2,3]
        coeffs = [1, 2, 2]
        L1 = LaguerrePolynomial(orders, coeffs)
        L2 = L1 * 3
        # multiplication by a number
        @test L2.orders == [0,2,3]
        @test L2.coeffs == 3 * coeffs
        # division by a number
        @test_throws ErrorException L1 / 0
        L3 = L1 / 2
        @test L3.orders == L1.orders
        @test L3.coeffs == L1.coeffs / 2
    end

    @testset "LaguerrePolynomials product" begin
        k = 5
        L2 = LaguerrePolynomial(2, 1.0)
        Lk = LaguerrePolynomial(k, 1.0)
        L = L2 * Lk
        # L2Lk = k(k-1)/2 L_{k-2} -2k(k-1)L_{k-1}
        #   +(3k^2-k)L_k -2k(k+1)L_{k+1}
        #   (k+1)(k+2)/2 L_{k+2}
        res = k * (k - 1) / 2 * LaguerrePolynomial(k - 2)  -
        2 * k * (k - 1) * LaguerrePolynomial(k - 1) +
        (3 * k^2 - k) * LaguerrePolynomial(k) -
        2 * k * (k + 1) * LaguerrePolynomial(k + 1) +
        (k + 1) * (k + 2) / 2 * LaguerrePolynomial(k + 2)
        @test L == res

        # L1Lk = k L_{k-1} -2k L_k +(k+1) L_{k+1}
        L1 = LaguerrePolynomial(1, 1.0)
        Lk = LaguerrePolynomial(k, 1.0)
        L = L1 * Lk
        res = k * LaguerrePolynomial(k - 1) - 2 * k * LaguerrePolynomial(k) + (k + 1) * LaguerrePolynomial(k + 1)
        @test L ==  res
        @test Lk * LaguerrePolynomial(0) == Lk
    end
end


text:
module BooksController

using Genie.Renderer, SearchLight, Books, Genie.Router, Genie.Requests
using ViewHelper

function billgatesbooks()
  html!(:books, :billgatesbooks, context = @__MODULE__, books = SearchLight.all(Book))
end

function new()
  html!(:books, :new, context = @__MODULE__, book = Book())
end

function create()
  cover_path = if haskey(filespayload(), "book_cover")
      path = joinpath("img", "covers", filespayload("book_cover").name)
      write(joinpath("public", path), IOBuffer(filespayload("book_cover").data))

      path
    else
      ""
  end

  Book( title = @params(:book_title),
        author = @params(:book_author),
        cover = cover_path) |> save && redirect_to(:get_bgbooks)
end

function edit()
  book = try
    SearchLight.find_one!!(Book, @params(:id))
  catch ex
    return error_404("Book ID $(@params(:id))")
  end

  html!(:books, :edit, book = book, context = @__MODULE__)
end

function update()
  # nothing
end

function index()
  @vars(:books, SearchLight.all(Book))
  html!("""

Bill's Gates top $(length(@vars(:books))) recommended books

  <%
    @foreach(@vars(:books)) do book
      partial("app/resources/books/views/book.jl.html", book = book)
    end
  %>

  """, context = @__MODULE__)
end

function flaxtest()
  books = SearchLight.all(Book)
  Flax.section([
    Flax.h1("Bill's Gates top $(length(books))")
    Flax.ul(
      @foreach(books) do book
        partial("app/resources/books/views/book.jl.html", context = @__MODULE__, book = book)
      end
    )
  ])
end

function foo()
  "mee"
end

### API

module API

using ..BooksController
using Genie.Renderer
using SearchLight, Books
using JSON

function billgatesbooks()
  json!(:books, :billgatesbooks, books = SearchLight.all(Book), context = @__MODULE__)
end

function jsontest()
  respond(JSON.json(SearchLight.all(Book)), "application/json")
end

end # API

end # BooksController

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
immutable CollectionShrinker <: AbstractShrinker end

collection{K,V}(::Type{Dict{K,V}}, itr=Pair{K,V}[]) = Dict(itr)
collection{E}(::Type{Vector{E}}, itr=E[]) = collect(itr)
collection{E}(::Type{Set{E}}, itr=E[]) = Set(itr)

function Base.UTF8String(itr)
    buf = IOBuffer()
    for el in itr
        write(buf, el)
    end
    s = takebuf_string(buf)
    return s
end

collection(::Type{UTF8String}, itr="") = UTF8String(itr)

typealias AbstractLengthStrategy{D <: AbstractUniform} Union{JustStrategy, NumberStrategy{D, Int}}

# ----- Vectors -----
immutable VectorStrategy{E, L <: AbstractLengthStrategy, S} <: AbstractStrategy
    element_strategy::E
    length_strategy::L
    shrinker::S
end

Base.minimum(strategy::VectorStrategy) = minimum(strategy.length_strategy)
returntype{E, L, S}(::Type{VectorStrategy{E, L, S}}) = Vector{returntype(E)}

function vectors(element_strategy::AbstractStrategy,
    length_strategy::AbstractStrategy = signeds(0, 10);
    shrinker=CollectionShrinker())

    VectorStrategy(
        element_strategy,
        length_strategy,
        shrinker)
end

function arbitrary(strategy::VectorStrategy)
    n = arbitrary(strategy.length_strategy)
    [arbitrary(strategy.element_strategy) for i=1:n]
end

# ----- Dicts -----
immutable DictStrategy{K, V, L <: AbstractLengthStrategy, S} <: AbstractStrategy
    key_strategy::K
    value_strategy::V
    length_strategy::L
    shrinker::S
end

Base.minimum(strategy::VectorStrategy) = minimum(strategy.length_strategy)
returntype{K, V, L, S}(::Type{DictStrategy{K, V, L, S}}) = Dict{returntype(K), returntype(V)}

function dicts(
    key_strategy::AbstractStrategy,
    value_strategy::AbstractStrategy,
    length_strategy::AbstractStrategy;
    shrinker=CollectionShrinker())

    DictStrategy(key_strategy, value_strategy, length_strategy, shrinker)
end

function arbitrary{K, V, L, S}(strategy::DictStrategy{K, V, L, S})
    key_strategy = strategy.key_strategy
    value_strategy = strategy.value_strategy
    length_strategy = strategy.length_strategy

    n = arbitrary(length_strategy)
    dict = collection(returntype(typeof(strategy)))
    sizehint!(dict, n)

    i = 1
    n_duplicates = 0
    max_duplicates = max(div(n, 2), 10)
    while (i <= n)
        k = arbitrary(key_strategy)
        if !haskey(dict, k)
            v = arbitrary(value_strategy)
            dict[k] = v
            i += 1
        else
            n_duplicates += 1
        end

        if n_duplicates > max_duplicates
            error("Not enough unique keys. $max_duplicates found")
        end
    end
    return dict
end

# ----- Sets -----
immutable SetStrategy{E, L <: AbstractLengthStrategy, S} <: AbstractStrategy
    element_strategy::E
    length_strategy::L
    shrinker::S
end

function sets(
    element_strategy::AbstractStrategy,
    length_strategy::AbstractStrategy;
    shrinker = CollectionShrinker())

    SetStrategy(element_strategy, length_strategy, shrinker)
end

Base.minimum(strategy::VectorStrategy) = minimum(strategy.length_strategy)
returntype{E, L, S}(::Type{SetStrategy{E, L, S}}) = Set{returntype(E)}

function arbitrary{E, L, S}(strategy::SetStrategy{E, L, S})
    element_strategy = strategy.element_strategy
    length_strategy = strategy.length_strategy

    n = arbitrary(length_strategy)
    set = collection(returntype(typeof(strategy)))
    sizehint!(set, n)

    i = 1
    n_duplicates = 0
    max_duplicates = max(div(n, 2), 10)
    while i <= n
        v = arbitrary(element_strategy)
        if !(v in set)
            push!(set, v)
            i += 1
        else
            n_duplicates += 1
        end
        if n_duplicates > max_duplicates
            error("Not enough unique values. $max_duplicates found")
        end
    end
    return set
end

# ----- Strings -----
immutable UTF8StringStrategy{L <: AbstractLengthStrategy, S <: AbstractShrinker} <: AbstractStrategy
    dist::DiscreteUniform
    length_strategy::L
    specialvalues::Vector{Char}
    p::Bernoulli
    shrinker::S
end

_specialvalues(::Type{Char}) = Char['\0', ' ', '\n' , '\t']
function utf8strings(
    min::Int=0, max::Int=65535,
    length_strategy::AbstractLengthStrategy=signeds(0, 10);
    specialvalues = _specialvalues(Char),
    p=Bernoulli(0.1),
    shrinker=CollectionShrinker())

    UTF8StringStrategy(DiscreteUniform(min, max), length_strategy, specialvalues, p, shrinker)
end

Base.minimum(strategy::UTF8StringStrategy) = minimum(strategy.length_strategy)
returntype{L, S}(::Type{UTF8StringStrategy{L, S}}) = UTF8String

function arbitrary(strategy::UTF8StringStrategy)
    n = arbitrary(strategy.length_strategy)
    buf = IOBuffer()
    for i=1:n
        if Bool(rand(strategy.p))
            write(buf, rand(strategy.specialvalues))
        else
            write(buf, Char(rand(strategy.dist)))
        end
    end
    s = takebuf_string(buf)
    return s
end

# ----- Tuples -----
immutable TupleShrinker <: AbstractShrinker end

immutable TupleStrategy{TE} <: AbstractStrategy
    element_strategies::TE
end

function tuples(element_strategies::AbstractStrategy...)
    TupleStrategy(element_strategies)
end

function arbitrary{TE}(strategy::TupleStrategy{TE})
    map(arbitrary, strategy.element_strategies)
end

shrink(x, strategy::TupleStrategy) = Task(() -> nothing)

# ----- Types -----
immutable TypeStrategy{TE} <: AbstractStrategy
    constructor
    element_strategies::TE
end

function record(constructor, element_strategies)
    TypeStrategy(constructor, element_strategies)
end

function arbitrary{TE}(strategy::TypeStrategy{TE})
    strategy.constructor(map(arbitrary, strategy.element_strategies)...)
end

# ----- General -----
typealias VectorBisectionStrategy{E, L}     VectorStrategy{E, L, CollectionShrinker}
typealias DictBisectionStrategy{K, V, L}    DictStrategy{K, V, L, CollectionShrinker}
typealias SetBisectionStrategy{E, L}        SetStrategy{E, L, CollectionShrinker}
typealias UTF8StringBisectionStrategy{L}    UTF8StringStrategy{L, CollectionShrinker}

typealias AbstractBisectionStrategy Union{
    VectorBisectionStrategy,
    DictBisectionStrategy,
    SetBisectionStrategy,
    UTF8StringBisectionStrategy}

shrink(x, strategy::AbstractBisectionStrategy) = Task() do
    min_collection_size = minimum(strategy)
    T = returntype(typeof(strategy))
    n = length(x)
    if n > min_collection_size
        new_size = max(div(n, 2), min_collection_size)
        can_bisect = div(n, 2) >= min_collection_size

        itr = take(x, new_size)
        produce(collection(T, itr))

        if can_bisect
            itr = drop(x, new_size)
            produce(collection(T, itr))
        else
            itr = drop(x, n - new_size)
            produce(collection(T, itr))
        end
    end
end

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
using Base.Test
using FixedPointNumbers
using Compat
import Compat.String

function test_op{F,T}(fun::F, ::Type{T}, fx, fy, fxf, fyf, tol)
    # Make sure that the result is representable
    (typemin(T) <= fun(fxf, fyf) <= typemax(T)) || return nothing
    @assert abs(fun(fx, fy) - convert(T, fun(fxf, fyf))) <= tol
    @assert abs(convert(Float64, fun(fx, fy)) - fun(fxf, fyf)) <= tol
end

function test_fixed{T}(::Type{T}, f)
    values = [-10:0.01:10; -180:.01:-160; 160:.01:180]
    tol = 2.0^-f
    for x in values
        # Ignore values outside the representable range
        # typemin <, otherwise for -(-0.5)  > typemax
        if !(typemin(T) < x <= typemax(T))
            continue
        end
        isinteger(x) && @show x
        fx = convert(T,x)
        @test convert(T,convert(Float64, fx)) == fx
        @test convert(T,convert(Float64, -fx)) == -fx
        @test convert(Float64, -fx) == -convert(Float64, fx)

        fxf = convert(Float64, fx)

        rx = convert(Rational{BigInt},fx)
        @assert isequal(fx,rx) == isequal(hash(fx),hash(rx))

        for y in values
            if !(typemin(T) < y <= typemax(T))
                continue
            end

            fy = convert(T,y)
            fyf = convert(Float64, fy)

            @assert fx==fy || x!=y
            @assert fx<fy  || x>=y
            @assert fx<=fy || x>y

            test_op(+, T, fx, fy, fxf, fyf, tol)
            test_op(-, T, fx, fy, fxf, fyf, tol)
            test_op(*, T, fx, fy, fxf, fyf, tol)
            fy != 0 && test_op(/, T, fx, fy, fxf, fyf, tol)

            @assert isequal(fx,fy) == isequal(hash(fx),hash(fy))
        end
    end
end

@test isapprox(convert(Fixed{Int8,7}, 0.8), 0.797, atol=0.001)
@test isapprox(convert(Fixed{Int8,7}, 0.9), 0.898, atol=0.001)
@test_throws InexactError convert(Fixed{Int8, 7}, 0.999)
@test_throws InexactError convert(Fixed{Int8, 7}, 1.0)
@test_throws InexactError convert(Fixed{Int8, 7}, 1)
@test_throws InexactError convert(Fixed{Int8, 7}, 2)
@test_throws InexactError convert(Fixed{Int8, 7}, 128)

for (TI, f) in [(Int8, 8), (Int16, 8), (Int16, 10), (Int32, 16)]
    T = Fixed{TI,f}
    println("  Testing $T")
    test_fixed(T, f)
end

T = Fixed{Int8,7}
for i = -1.0:0.1:typemax(T)
    @test i % T === T(i)
end
@test ( 1.5 % T).i == round(Int,  1.5*128) % Int8
@test (-0.3 % T).i == round(Int, -0.3*128) % Int8

T = Fixed{Int16,9}
for i = -64.0:0.1:typemax(T)
    @test i % T === T(i)
end
@test ( 65.2 % T).i == round(Int,  65.2*512) % Int16
@test (-67.2 % T).i == round(Int, -67.2*512) % Int16

for T in [Fixed{Int8,7}, Fixed{Int16,8}, Fixed{Int16,10}]
    testapprox(T)  # defined in ufixed.jl
end

# reductions
F8 = Fixed{Int8,8}
a = F8[0.498, 0.1]
acmp = Float64(a[1]) + Float64(a[2])
@test sum(a) == acmp
@test sum(a, 1) == [acmp]

F6 = Fixed{Int8,6}
a = F6[1.2, 1.4]
acmp = Float64(a[1])*Float64(a[2])
@test prod(a) == acmp
@test prod(a, 1) == [acmp]

x = Fixed{Int8,8}(0.3)
for T in (Float16, Float32, Float64, BigFloat)
    y = convert(T, x)
    @test isa(y, T)
end

# Floating-point conversions
@test isa(float(one(Fixed{Int8,6})),   Float32)
@test isa(float(one(Fixed{Int32,18})), Float64)
@test isa(float(one(Fixed{Int32,25})), Float64)

# Show
x = Fixed{Int32,5}(0.25)
iob = IOBuffer()
show(iob, x)
str = String(take!(iob))
@test str == "0.25Q26f5"
@test eval(parse(str)) == x

for T in (Fixed{Int8,8}, Fixed{Int16,8}, Fixed{Int16,10}, Fixed{Int32,16})
    a = rand(T)
    @test isa(a, T)
    a = rand(T, (3, 5))
    @test isa(a, Array{T,2})
    @test size(a) == (3,5)
end


text:
using JLD2
using DataStructures
using StaticArrays
include("automata.jl")
include("nrm_cpu.jl")

function instantiate_model()

    propRules = [x -> x.rates[1]*x.reactants[1]                         # mh1 -> ph1
                x -> x.rates[5]*x.reactants[4]                          # ph1 -> 0
                x -> x.rates[30]*x.reactants[4]*((x.reactants[4]-1)/2)  # ph1 + ph1 -> ph11
                x -> x.rates[31]*x.reactants[7]                         # ph11 -> ph1 + ph1
                x -> x.rates[34]*x.reactants[4]*x.reactants[5]          # ph1 + ph7 -> ph17
                x -> x.rates[35]*x.reactants[9]                         # ph17 -> ph1 + ph7
                x -> x.rates[32]*x.reactants[4]*x.reactants[6]          # ph1 + ph6 -> ph16
                x -> x.rates[33]*x.reactants[8]                         # ph16 -> ph1 + ph6
                x -> x.rates[3]*x.reactants[2]                          # mh7 -> ph7
                x -> x.rates[7]*x.reactants[5]                          # ph7 -> 0
                x -> x.rates[40]*x.reactants[5]*(x.reactants[5]-1)/2    #ph7 + ph7 -> ph77
                x -> x.rates[41]*x.reactants[12]                        # ph77 -> ph7 + ph7
                x -> x.rates[38]*x.reactants[5]*x.reactants[6]          # ph7 + ph6 -> ph76
                x -> x.rates[39]*x.reactants[11]                        # ph76 -> ph7 + ph6
                x -> x.rates[2]*x.reactants[3]                          # mh6 -> ph6
                x -> x.rates[6]*x.reactants[6]                          # ph6 -> 0
                x -> x.rates[36]*x.reactants[6]*(x.reactants[6]-1)/2    # ph6 + ph6 -> ph66
                x -> x.rates[37]*x.reactants[10]                        # ph66 -> ph6 + ph6
                x -> x.rates[17]*x.reactants[7]                         # ph11 -> 0
                x -> x.rates[19]*x.reactants[9]                         # ph17 -> 0
                x -> x.rates[18]*x.reactants[8]                         # ph16 -> 0
                x -> x.rates[22]*x.reactants[12]                        # ph77 -> 0
                x -> x.rates[21]*x.reactants[11]                        # ph76 -> 0
                x -> x.rates[20]*x.reactants[10]                        # ph66 -> 0
                x -> x.rates[4]*x.reactants[13]                         # md -> pd
                x -> x.rates[8]*x.reactants[14]                         # pd -> 0
                x -> x.rates[9]*((1+(x.reactants[14]/x.rates[44]))/(1 + ((x.reactants[14]/x.rates[44])) + (x.reactants[7]/x.rates[42])^2 + (x.reactants[11]/x.rates[43])^2)) # fh1
                x -> x.rates[13]*x.reactants[1]                         # mh1 -> 0
                x -> x.rates[11]*((1+(x.reactants[14]/x.rates[44]))/(1 + ((x.reactants[14]/x.rates[44])) + (x.reactants[7]/x.rates[42])^2 + (x.reactants[11]/x.rates[43])^2)) # fh7
                x-> x.rates[15]*x.reactants[2]                          # mh7 -> 0
                x -> x.rates[10]                                        # 0 -> mh6
                x -> x.rates[14]*x.reactants[3]                         # mh6 -> 0
                x -> x.rates[12]*(1/(1 + (x.reactants[7]/x.rates[42])^2 + (x.reactants[11]/x.rates[43])^2)) # fd
                x -> x.rates[16]*x.reactants[13]]                       # md -> 0

    stoichiometry = [x -> insert!(x.pq_delayed[1], x.pos[1], (x.t[1] + x.rates[26]))    # mh1 -> ph1
                x -> x.reactants .+= @SMatrix [0 0 0 -1 0 0 0 0 0 0 0 0 0 0]    # ph1 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 -2 0 0 1 0 0 0 0 0 0 0]    # ph1 + ph1 -> ph11
                x -> x.reactants .+= @SMatrix [0 0 0 2 0 0 -1 0 0 0 0 0 0 0]    # ph11 -> ph1 + ph1
                x -> x.reactants .+= @SMatrix [0 0 0 -1 -1 0 0 0 1 0 0 0 0 0]   # ph1 + ph7 -> ph17
                x -> x.reactants .+= @SMatrix [0 0 0 1 1 0 0 0 -1 0 0 0 0 0]    # ph17 -> ph1 + ph7
                x -> x.reactants .+= @SMatrix [0 0 0 -1 0 -1 0 1 0 0 0 0 0 0]   # ph1 + ph6 -> ph16
                x -> x.reactants .+= @SMatrix [0 0 0 1 0 1 0 -1 0 0 0 0 0 0]    # ph16 -> ph1 + ph6
                x -> insert!(x.pq_delayed[2], x.pos[2], (x.t[1] + x.rates[28]))         # mh7 -> ph7
                x -> x.reactants .+= @SMatrix [0 0 0 0 -1 0 0 0 0 0 0 0 0 0]    # ph7 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 -2 0 0 0 0 0 0 1 0 0]    # ph7 + ph7 -> ph77
                x -> x.reactants .+= @SMatrix [0 0 0 0 2 0 0 0 0 0 0 -1 0 0]    # ph77 -> ph7 + ph7
                x -> x.reactants .+= @SMatrix [0 0 0 0 -1 -1 0 0 0 0 1 0 0 0]   # ph7 + ph6 -> ph76
                x -> x.reactants .+= @SMatrix [0 0 0 0 1 1 0 0 0 0 -1 0 0 0]    # ph76 -> ph7 + ph6
                x -> insert!(x.pq_delayed[3], x.pos[3], (x.t[1] + x.rates[27]))         # mh6 -> ph6
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 -1 0 0 0 0 0 0 0 0]    # ph6 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 -2 0 0 0 1 0 0 0 0]    # ph6 + ph6 -> ph66
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 2 0 0 0 -1 0 0 0 0]    # ph66 -> ph6 + ph6
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 -1 0 0 0 0 0 0 0]    # ph11 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 -1 0 0 0 0 0]    # ph17 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 -1 0 0 0 0 0 0]    # ph16 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 0 -1 0 0]    # ph77 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 -1 0 0 0]    # ph76 -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 -1 0 0 0 0]    # ph66 -> 0
                x -> insert!(x.pq_delayed[4], x.pos[4], (x.t[1] + x.rates[29]))         # md -> pd
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 0 0 0 -1]    # pd -> 0
                x -> insert!(x.pq_delayed[5], x.pos[5], (x.t[1] + x.rates[23]))         # 0 -> mh1
                x -> x.reactants .+= @SMatrix [-1 0 0 0 0 0 0 0 0 0 0 0 0 0]    # mh1 -> 0
                x -> insert!(x.pq_delayed[6], x.pos[6], (x.t[1] + x.rates[24]))         # 0 -> mh7
                x -> x.reactants .+= @SMatrix [0 -1 0 0 0 0 0 0 0 0 0 0 0 0]    # mh7 -> 0
                x -> insert!(x.pq_delayed[7], x.pos[7], (x.t[1] + 0))                   # 0 -> mh6
                x -> x.reactants .+= @SMatrix [0 0 -1 0 0 0 0 0 0 0 0 0 0 0]    # mh6 -> 0
                x -> insert!(x.pq_delayed[8], x.pos[8], (x.t[1] + x.rates[25]))         # 0 -> md
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 0 0 -1 0]    # md -> 0
                x -> x.reactants .+= @SMatrix [0 0 0 1 0 0 0 0 0 0 0 0 0 0]    # mh1 -> ph1
                x -> x.reactants .+= @SMatrix [0 0 0 0 1 0 0 0 0 0 0 0 0 0]    # mh7 -> ph7
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 1 0 0 0 0 0 0 0 0]    # mh6 -> ph6
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 0 0 0 1]    # md -> pd
                x -> x.reactants .+= @SMatrix [1 0 0 0 0 0 0 0 0 0 0 0 0 0]     # 0 -> mh1
                x -> x.reactants .+= @SMatrix [0 1 0 0 0 0 0 0 0 0 0 0 0 0]     # 0 -> mh7
                x -> x.reactants .+= @SMatrix [0 0 1 0 0 0 0 0 0 0 0 0 0 0]     # 0 -> mh6
                x -> x.reactants .+= @SMatrix [0 0 0 0 0 0 0 0 0 0 0 0 1 0]]    # 0 -> md
    model = Model(Vector{Cell}(), Vector{Float64}(),propRules,stoichiometry)
    return model
end

function simulate!(cell::Cell, m::Model, tmax::Float64, T_map::Dict{Integer, Integer}, D_map::Dict{Integer, Integer}, dgraph::Dict{Int64,Array{Int64}})
    potential_reactions = @SVector [1,9,15,25,27,29,31,33]
    while (cell.t[1] <= tmax)
        dt, event = jump(cell, cell.t[1])
        cell.t[1] += dt
        m.stoichiometry[event](cell)
        if event in potential_reactions
            update_ndelayed!(cell, T_map[event])
        end
        if (event == 39 || event == 28)
            store!(cell)
        end
        update_tk!(cell,dt)
        if event <= 34
            cell.Pk[event] += log(1/rand())
            calculate_propensity!(cell, m.propRules, dgraph[event])
        else
            calculate_propensity!(cell, m.propRules, dgraph[D_map[event]])
        end
        update_tjump!(cell)
    end
end

function fill_model!(ncell::Int, model::Model, r::Array{Float64})
    for i in 1:ncell
        push!(model.cells, Cell(reactants=[0 0 0 0 0 100 0 0 0 0 0 0 0 0],rates=r))
    end
end

function nrm(ncell=1, tmax=0.05,
    rates = [37.1,47.874,16.837,43.535,0.29825,0.25089,0.29916,0.30915,57.265,44.116,47.455,59.318,0.33807,0.19722,0.38719,0.18795,0.26387,0.28755,0.2728,0.30433,0.27177,0.28934,9.8992,8.7802,8.915,1.6025,0.91232,1.6195,10.967,0.017267,0.24785,0.0277,0.10982,0.0011008,0.24551,0.012698,0.28028,0.021359,0.080509,0.0054418,0.13634,711.49,280.18,511.8])
    m = instantiate_model()
    fill_model!(ncell, m, rates)
    initialize_propensities!(m)

    # Load the dependency graph.
    @load "src/dependency_graph.jld2" dgraph

    T_map = Dict{Integer, Integer}(1 => 1, 9 => 2, 15 => 3, 25 => 4, 27 => 5, 29 => 6, 31 => 7, 33 => 8)
    D_map = Dict{Integer, Integer}(35 => 1, 36 => 9, 37 => 15, 38 => 25, 39 => 27, 40 => 29, 41 => 31, 42 => 33)
    #nmax = 500000
    Threads.@threads for cell in m.cells
        simulate!(cell, m, convert(Float64,tmax), T_map, D_map, dgraph)
    end
    return m
end

max_stars_repo_name: syedtaz/project-cait

max_issues_repo_name: syedtaz/project-cait

max_forks_repo_name: syedtaz/project-cait


text:
## Low level gemm! call with pointers
## Borrowed from Knet.jl, adapted for compile-time constants

using LinearAlgebra
using LinearAlgebra.BLAS: libblas, BlasInt, @blasfunc

"""
    gemm!()

Low-level gemm!() call with pointers, borrowed from Knet.jl

Calculates `C = alpha*op(A)*op(B) + beta*C`, where:
  - `transA` and `transB` set `op(X)` to be either `identity()` or `transpose()`
  - alpha and beta are scalars
  - op(A) is an (M, K) matrix
  - op(B) is a (K, N) matrix
  - C is an (M, N) matrix.
"""
gemm!

# These are the datatypes we have fast GEMM for
gemm_datatype_mappings = (
    (:dgemm_, Float64),
    (:sgemm_, Float32),
    (:zgemm_, ComplexF64),
    (:cgemm_, ComplexF32),
)
for (gemm, elt) in gemm_datatype_mappings
    @eval begin
        @inline function gemm!(transA::Val, transB::Val,
                               M::Int, N::Int, K::Int,
                               alpha::$(elt), A::Ptr{$elt}, B::Ptr{$elt},
                               beta::$(elt), C::Ptr{$elt})
            # Convert our compile-time transpose marker to a char for BLAS
            convtrans(V::Val{false}) = 'N'
            convtrans(V::Val{true})  = 'T'

            if transA == Val(false)
                lda = M
            else
                lda = K
            end
            if transB == Val(false)
                ldb = K
            else
                ldb = N
            end
            ldc = M
            ccall((@blasfunc($(gemm)), libblas), Nothing,
                  (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ref{BlasInt}, Ref{$elt}, Ptr{$elt}, Ref{BlasInt},
                   Ptr{$elt}, Ref{BlasInt}, Ref{$elt}, Ptr{$elt},
                   Ref{BlasInt}),
                  convtrans(transA), convtrans(transB), M, N, K,
                  alpha, A, lda, B, ldb, beta, C, ldc)
        end
    end
end


text:
module Day16
export parse_input, part_1, part_2

struct Instruction
	opcode::Int
	A::Int
	B::Int
	C::Int
end

function Base.parse(::Type{Instruction}, input::AbstractString)::Instruction
    opcode, A, B, C = map(x -> parse(Int, x), split(input))
    return Instruction(opcode, A, B, C)
end

struct Example
	before::Array{Int, 1}
	instruction::Instruction
	after::Array{Int, 1}
end

function Base.parse(::Type{Example}, input::AbstractString)::Example
	before, instruction, after = split(input, "\n")
    return Example(
        parse_register(before),
        parse(Instruction, instruction),
        parse_register(after),
    )
end

function parse_register(register::AbstractString)
    m = match(r"\[(\d+), (\d+), (\d+), (\d+)\]", register)
	return map(x -> parse(Int, x), m.captures)
end

function parse_input(input::AbstractString)
	part1, part2 = split(input, "\n\n\n")

	example_strings = split(part1, "\n\n")

    examples = map(x -> parse(Example, x), example_strings)

    test_program = map(x-> parse(Instruction, x), split(strip(part2), "\n"))

    return examples, test_program
end

function part_1(examples)
    n_geq_3 = 0
    for example in examples
        n_fitting = 0
        for operation! in operations
            if test_operation(operation!, example)
                n_fitting += 1
            end
        end

        if n_fitting >= 3
            n_geq_3 += 1
        end
    end

    return n_geq_3
end

function part_2(examples, test_program)
    opcodes = derive_opcode_meaning(examples)
    return execute_program(test_program, opcodes)[1]
end

function execute_program(program, opcodes)
    register = zeros(Int, 4)
    for instruction in program
        opcodes[instruction.opcode](register, instruction)
    end
    return register
end

function test_operation(operation!, example)
    register = copy(example.before)
    operation!(register, example.instruction)
    return register == example.after
end

function derive_opcode_meaning(examples)
    opcodes = Dict(i => copy(operations) for i in 0:15)

    for example in examples
        opcode = example.instruction.opcode
        filter!(o -> test_operation(o, example), opcodes[opcode])
    end

    unknown = opcodes
    known = Dict()

    while !all(length(o) == 1 for o in values(opcodes))
        for (opcode, possible) in unknown
            if length(possible) == 1
                operation =  possible[1]
                known[opcode] = operation
                pop!(unknown, opcode)

                for (opcode, possible) in unknown
                    filter!(x -> x != operation, possible)
                end
            end
        end
    end

    return known
end

function addr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] + register[inst.B + 1]
end

function addi!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] + inst.B
end

function mulr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] * register[inst.B + 1]
end

function muli!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] * inst.B
end

function banr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] & register[inst.B + 1]
end

function bani!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] & inst.B
end

function borr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] | register[inst.B + 1]
end

function bori!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] | inst.B
end

function setr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1]
end

function seti!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = inst.A
end

function gtir!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] =  inst.A > register[inst.B + 1]
end

function gtri!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] > inst.B
end

function gtrr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] > register[inst.B + 1]
end

function eqir!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = inst.A == register[inst.B + 1]
end

function eqri!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] == inst.B
end

function eqrr!(register::Array{Int}, inst::Instruction)
    register[inst.C + 1] = register[inst.A + 1] == register[inst.B + 1]
end

operations = [addr!, addi!, mulr!, muli!, banr!, bani!, borr!, bori!, setr!, seti!, gtir!, gtri!, gtrr!, eqir!, eqri!, eqrr!]

end