text:
# https://www.cp.eng.chula.ac.th/~attawith/class/linint
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, DataStructures

"""
	presorted(V, EV [,preview=false] )

Generate the presorting needed by the scan-line algorithm (Bentley-Ottmann, 1979).

When two or more 2D (left to right oriented) lines have the same start point,
they are ordered according to forward values of the y-coord of the end line point.
This guarantess that each other line having the same first vertex is inserted in the right
scan-line data structure interval, i.e. between the two lines closely above and below it.

# Example

```
julia> lines = [[Float64[1,2],Float64[4,3]],
       [Float64[1,2],Float64[4,2.999]]]
2-element Array{Array{Array{Float64,1},1},1}:
 [[1.0, 2.0], [3.0, 3.0]]
 [[1.0, 2.0], [4.0, 2.999]]

julia> isless0 = (x,y) -> [x[1][1],x[1][2],x[2][2]] < [y[1][1],y[1][2],y[2][2]]
#66 (generic function with 1 method)

julia> sort(lines, lt=isless0)
2-element Array{Array{Array{Float64,1},1},1}:
 [[1.0, 2.0], [4.0, 2.999]]
 [[1.0, 2.0], [3.0, 3.0]]
```
"""
function presorted(V,EV, preview=false)
	lines = [[V[:,u],V[:,v]] for (u,v) in EV]
	orientedlines = Array{Array{Float64,1},1}[]
	for (v1,v2) in lines
		if v1 > v2
			v1,v2 = v2,v1
		end
		push!(orientedlines, [v1,v2])
	end
	isless0 = (x,y) -> [x[1][1],x[1][2],x[2][2]] < [y[1][1],y[1][2],y[2][2]]
	inputlines = sort(orientedlines, lt=isless0)
	newlines = map(cat,inputlines)
	W,EW  = lines2lar(newlines)
	if preview
		Plasm.view(Plasm.numbering(.1)((W,[[[k] for k=1:size(W,2)], EW ])))
	end
	vsorted = sort([[W[:,k],k] for k=1:size(W,2)])
	V = hcat([vsorted[k][1] for k=1:length(vsorted)]...)
	edgedict = OrderedDict([(oldv,newv) for (newv,(vert,oldv)) in enumerate(vsorted)])
	EV = [[edgedict[v1],edgedict[v2]] for (v1,v2) in EW]
	return V,EV
end

function lines2lar(lines)
	vertdict = OrderedDict{Array{Float64,1}, Int64}()
	EV = Array{Int64,1}[]
	idx = 0
	for h=1:length(lines)
		x1,y1,x2,y2 = cat(lines[h])

		if ! haskey(vertdict, [x1,y1])
			idx += 1
			vertdict[[x1,y1]] = idx
		end
		if ! haskey(vertdict, [x2,y2])
			idx += 1
			vertdict[[x2,y2]] = idx
		end
		v1,v2 = vertdict[[x1,y1]],vertdict[[x2,y2]]
		push!(EV, [v1,v2])
	end
	V = hcat(collect(keys(vertdict))...)
	return V,EV
end

function cuboids(n,scale=1.; grid=[1,1])
	assembly = []
	for k=1:n
		corner = rand(Float64, 2)
		dims = rand(Float64, 2)
		if grid==[1,1]
			V,(_,EV,_) = Lar.cuboid(corner,true,corner+dims)
		else
			V,(_,EV,_) = Lar.cuboidGrid(grid, true)
			g = Lar.Struct([ Lar.t(corner...), Lar.s(dims...), (V,EV) ])
			V,EV = Lar.struct2lar(g)
		end
		center = (corner + corner+dims)/2
		angle = rand(Float64)*2*pi
		obj = Lar.Struct([ Lar.t(center...), Lar.r(angle),
				Lar.s(scale,scale), Lar.t(-center...), (V,EV) ])
		push!(assembly, obj)
	end
	Lar.struct2lar(Lar.Struct(assembly))
end

function presortedlines(V,EV)
	lines = [[V[:,u],V[:,v]] for (u,v) in EV]
	orientedlines = Array{Array{Float64,1},1}[]
	for (v1,v2) in lines
		if v1 > v2
			v1,v2 = v2,v1
		end
		push!(orientedlines, [v1,v2])
	end
	outlines = sort(orientedlines)
end

"""
	intersection(V::Lar.Points, EV::Lar.Cells)
		(e1::Int, e2::Int)::Union{Nothing, Array}

Compute the intersection point of two line segments `line1` and `line2`,
of extreme point indices `EV[e1]` and `EV[e2]`, with given coordinates
in `V` array.

# Example

```julia
julia> line1 = V[:,EV[14]]
2×2 Array{Float64,2}:
 0.0379479  0.168838
 0.90315    1.0

julia> line2 = V[:,EV[23]]
2×2 Array{Float64,2}:
 0.0307485  0.312532
 0.93606    0.957881

julia> intersection(line1,line2)
2-element Array{Float64,1}:
 0.08846469268137241
 0.9405292234774614

julia> intersection(line1,line2)
2-element Array{Float64,1}:
 0.08846469268137241
 0.9405292234774614

julia> intersection(line1,-1*line2) # returns a `Nothing` value

julia>
```
"""
function intersection(line1,line2)::Union{Nothing, Array}
	x1,y1,x2,y2 = vcat(line1...)
	x3,y3,x4,y4 = vcat(line2...)

	# intersect lines e1,e2
	# 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
		p = [x1,y1] + α * ([x2,y2] - [x1,y1])
		return p
	end
end

"""
	computesegposition(SL,h,xe,ye)

Compute the current position of a neighbourhood line with respect to a new segment.

`SL` (the sweep line) is a SortedMultiDict; `h` is the semitoken of a segment on SL;
`(xe, ye)` are the coordinates of vertex of a new segment lying on `SL`, w.r.t.
compare the segment's `h` position.  The function returns a pair of values for `segA`
and `segB`, one of which is the empty array.
"""
function computesegposition(SL,h,xe,ye)
	segA = Array{Float64,1}[];
	segB = Array{Float64,1}[]
	u,v,type,edgeid = deref_value((SL,h))
	Dx = xe - u[1]
	Dy = Dx * (v[2]-u[2])/(v[1]-u[1])  ##### NO vertical edges allowed !!!
	y = u[2] + Dy
	if y > ye
		push!(segA, u, v, [edgeid])
	elseif y < ye
		push!(segB, u, v, [edgeid])
	elseif y == ye
		nothing::Nothing  # by now ... ;-)
	end
	return segA,segB
end

"""
	selectsegmentneighbors(SL::SortedMultiDict,i)

Compute `segA` and/or `segB`, i.e. the one/two segments adjacent to `segE`,
the segment of `E` event, on the `SL` sweepline. `segA` should be *above* `segE`;
`segB` should be *below* `segE`.
"""
function selectsegmentneighbor(SL,i)
	seghA, seghB = Array{Float64,1}[], Array{Float64,1}[]
	segkA, segkB = Array{Float64,1}[], Array{Float64,1}[]

	h = regress((SL,i))
	if h ≠ beforestartsemitoken(SL)
		seghA,seghB = computesegposition(SL,h,xe,ye)
		k = regress((SL,h))
		if k ≠ beforestartsemitoken(SL)
			segkA,segkB = computesegposition(SL,k,xe,ye)
		end
	end
	segA = union(seghA, segkA); segB = union(seghB, segkB);
	return segA, segB
end

"""
	swapsegments(SL,I,segA,segB)

Swap the positions of segments on the sweep line `SL` at the intersection point `I`.
"""
function swapsegments(SL,I,segA,segB)
	println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
	@show segA
	@show segB
	@show I
	# cancel both segA and segB from SL
	# ---------------------------------
	# get segment keys in SL
	a = segA[1]
	b = segB[1]
	# semitokens of segA and segB in SL
	a_st = searchsortedfirst(SL,a)
	b_st = searchsortedfirst(SL,b)
	if a_st ≠ pastendsemitoken(SL) && a_st ≠ pastendsemitoken(SL)
		# recover values of old segments from SL
		k_a,v_a = deref((SL, a_st))
		k_b,v_b = deref((SL, b_st))
		# delete old segments from SL
		delete!((SL, a_st))
		delete!((SL, b_st)) # delete old segment
		# swap the two new segments
		Dy = 0.0001
		newsegA = (I, v_a[2], v_a[3], v_a[4])
		newsegB = (I, v_b[2], v_b[3], v_b[4])
		if k_a[1] <= k_b[1]
			newsegA = (I+[0,Dy], v_a[2], v_a[3], v_a[4])
			keyA = (I+[0,Dy],"int",v_a[4])
		elseif k_a[1] > k_b[1]
			newsegB = (I+[0,Dy], v_b[2], v_b[3], v_b[4])
			keyB = (I+[0,Dy],"int",v_b[4])
		end
		# insert the new segments in SL
		insert!(SL, keyA, newsegA)
		insert!(SL, keyB, newsegB)
	end
	println("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@")
	return (keyA, newsegA),(keyB, newsegB)
end

function sweepline(V,EV)
	# event generation and ordering
	segments = presortedlines(V,EV)
	evpairs = [[(v1,v2,"start",k), (v2,v1,"end",k)] for (k,(v1,v2)) in enumerate(segments)]
	isless0 = (x,y) -> [x[1][1],x[1][2],x[2][2]] < [y[1][1],y[1][2],y[2][2]]
	events = sort(cat(evpairs),lt=isless0)
	# Initialize event queue ξ = all segment endpoints; Sort ξ by increasing x and y
	pqkeys = [(e[1],e[3],e[4]) for e in events]
	pqvalues = events
	ξ = PriorityQueue(zip(pqkeys,pqvalues))
	# Proper ordering object for the `SortedMultiDic` used by sweep line SL
	# Initialize sweep line SL to be empty
	SL = SortedMultiDict{Any,Any}()
	# Initialized output intersection list Λ to be empty
	Λ = Array{Array{Float64,1},1}[]

	while length(ξ) ≠ 0 # (ξ is nonempty)
@show peek(ξ)   # the next value of event (k => v) from ξ
		E = peek(ξ)[2]   # the next value of event (k => v) from ξ
		(v1, v2, nodetype, edgeid) = E
		# segE = E's segment
		e, segE = v1, (v1, v2, nodetype, edgeid) # key, value
		if E[3] == "start" # (E is a left endpoint)
			# Add segE to SL
			i = insert!(SL, e, segE)  # SortedMultiDict, key, value -> semitoken
			global xe,ye = e
			if first(SL) !== last(SL) # more than one segment in SL
				# compute segA and/or segB
				segA, segB = selectsegmentneighbor(SL,i)
				if segA ≠ []
				# segA = the segment above segE in SL
					a = Int(segA[3][1]); segA = segA[1:2] # edgeid of segA
					I = intersection(segE,segA)
					# (if Intersect( segE with segA) exists)
					if typeof(I) ≠ Nothing
						# Insert I into ξ
						#(keyE,newsegE),(keyA,newsegA) = swapsegments(SL,I,segE,segA)
						key = (I,"int",edgeid); val = (I,I,"int",a)
						enqueue!(ξ, key,val)
						#enqueue!(ξ, keyA,newsegA)
					end
				end
				if segB ≠ []
				# segB = the segment below segE in SL
					b = Int(segB[3][1]); segB = segB[1:2] # edgeid of segB
					I = intersection(segE,segB)
					# (if Intersect( segE with segB) exists)
					if typeof(I) ≠ Nothing
						# Insert I into ξ
						#(keyE,newsegE),(keyB,newsegB) = swapsegments(SL,I,segE,segB)
						key = (I,"int",edgeid); val = (I,I,"int",b)
						enqueue!(ξ, key,val)
						#enqueue!(ξ, keyB,newsegB)
					end
				end
			end
		elseif E[3] == "end" # (E is a right endpoint)
			# segE = E's segment
			# e = key(segE); st1 = corresponding semitoken in SL

			# compute E semitoken
			key = E[1]
			E_st = searchsortedlast(SL,key)

			# compute segA and/or segB
			segA, segB = selectsegmentneighbor(SL,E_st)
			# segA = the segment above segE in SL
			if segA ≠ []
				a = Int(segA[3][1]) # edgeid of segA
				segA = segA[1:2]
			end
			# segB = the segment below segE in SL
			if segB ≠ []
				b = Int(segB[3][1])
				segB = segB[1:2] # edgeid of segB
			end
			# Remove segE from SL
			delete!((SL,E_st))
			# (I = Intersect( segA with segB) exists)
			if segA ≠ [] && segB ≠ []
				I = intersection(segA, segB)
				if typeof(I) ≠ Nothing
					# (I is not in ξ already)
					# no problem anyway (in case I is overwritten)
					# Insert I into ξ
					#(keyA,newsegA),(keyB,newsegB) = swapsegments(SL,I,segA,segB)
					key = (I,"int",a); val = (I,I,"int",b)
					enqueue!(ξ, key,val)
					#enqueue!(ξ, keyB,newsegB)
				end
			end
		else # E is an intersection event
		@assert E[3] == "int"
			# Add E to the output list Λ
			push!(Λ, [E[1],E[2]])
			# segE1 above segE2 be E's intersecting segments in SL
			##i = insert!(SL, e, segE) ### ????

			# compute E semitoken
			key = E[1]
			E_st = searchsortedlast(SL,key)

			segE1, segE2 = selectsegmentneighbor(SL,E_st)
			if segE1 ≠ []
				h = Int(segE1[3][1]) # index of segment segE1
				v1 = V[:,EV[h][1]] # key of segE1 in SL
				hst1,hst2 = searchequalrange(SL::SortedMultiDict, v1)
				# extreme semitokens of segE1 in SL
				# segB = the segment below segE1 in SL
				segE1, segB = selectsegmentneighbor(SL,hst1)
				if segE1 ≠ [] && segB ≠ []
					I = intersection(segE1,segB)
					# (I = Intersect(segE1 with segB) exists)
					if typeof(I) ≠ Nothing
						# (if I is not in ξ already) no problem in case
						# Insert I into ξ
						(keyE1,newsegE1),(keyB,newsegB) = swapsegments(SL,I,segE1,segB)
						key = (I,"int",edgeid); val = (I,I,"int",b) # edgeid ??
						enqueue!(ξ, keyE1,newsegE1)
						enqueue!(ξ, keyB,newsegB)
					end
				end
			end
			if segE2 ≠ []
				k = Int(segE2[3][1])  # index of segment segE2
				v1 = V[:,EV[k][1]]  # key of segE2 in SL
				kst1,kst2 = searchequalrange(SL::SortedMultiDict, v1)
				# extreme semitokens of segE2 in SL
				# Swap their positions so that segE2 is now above segE1
				# ??? HOW TO DO IT ???
				# segA = the segment above segE2 in SL
				segA, segE2 = selectsegmentneighbor(SL,kst1)
				if segA ≠ [] && segE2 ≠ []  ## debug => false ==> don't compute last intersection
					I = intersection(segA, segE2)
					# (I = Intersect(segE2 with segA) exists)
					if typeof(I) ≠ Nothing
						# (if I is not in ξ already) no problem in case
						# Insert I into ξ
						(keyA,newsegA),(keyE2,newsegE2) = swapsegments(SL,I,segA,segE2)
						##key = (I,"int",edgeid); val = (I,I,"int",a) # edgeid ??
						enqueue!(ξ, keyA,newsegA)
						enqueue!(ξ, keyE2,newsegE2)
					end
				end
			end
			##delete!((SL, i)) ### ????
		end
		# remove E from ξ
		dequeue!(ξ)
	end
	return Λ
end

# EXAMPLE 0
# data generation
lines = [[[1,3],[10,5]],[[2,6],[11,2.5]],
[[3,4],[8,2.25]],[[5,1.25],[12,5]]]
V,EV = lines2lar(lines)
Plasm.view(Plasm.numbering(2.)((V,[[[k] for k=1:size(V,2)], EV])))
# data sorting
V = Plasm.normalize(V,flag=true)
W,EW = presorted(V,EV)
Plasm.view(Plasm.numbering(.25)((W,[[[k] for k=1:size(W,2)], EW ])))

#=
# EXAMPLE 1
# data generation
V,EV = cuboids(10, .5)
V = Plasm.normalize(V,flag=true)
model2d = V,EV
Plasm.view(Plasm.numbering(.1)((V,[[[k] for k=1:size(V,2)], EV])))
# data sorting
W,EW = presorted(V,EV)
Plasm.view(Plasm.numbering(.15)((W,[[[k] for k=1:size(W,2)], EW ])))

# EXAMPLE 2
# data generation
V,(_,EV,FV) = Lar.cuboidGrid([3,3],true)
Plasm.view(Plasm.numbering(1)((V,[[[k] for k=1:size(V,2)], EV ])))
# data sorting
W,EW = presorted(V,EV)
Plasm.view(Plasm.numbering(1)((W,[[[k] for k=1:size(W,2)], EW ])))
=#

#=
V = [0.435352 0.449902 0.658165 0.672715 0.76547 0.685644 0.542851 0.463025 0.438641 0.446501 0.406689 0.414549 0.145006 0.104172 0.218371 0.177538 0.129832 0.11199 0.151923 0.134081 0.206417 0.174285 0.208523 0.176391 0.988503 0.779401 0.904498 0.695396 0.0573599 0.232703 0.122291 0.297634 0.412384 0.42493 0.578685 0.59123 0.290235 0.0044308 0.285804 0.0; 0.290109 0.272658 0.475885 0.458435 0.274124 0.513271 0.199815 0.438962 0.813539 0.876549 0.817525 0.880535 0.947022 0.879789 0.902464 0.835231 0.416359 0.365524 0.408606 0.357771 0.345616 0.120288 0.345315 0.119988 0.286924 0.348144 0.0 0.0612204 0.416447 0.157502 0.460414 0.20147 0.759923 0.751232 1.0 0.99131 0.409458 0.428899 0.34432 0.363761]

EV = Array{Int64,1}[[1, 2], [3, 4], [1, 3], [2, 4], [5, 6], [7, 8], [5, 7], [6, 8], [9, 10], [11, 12], [9, 11], [10, 12], [13, 14], [15, 16], [13, 15], [14, 16], [17, 18], [19, 20], [17, 19], [18, 20], [21, 22], [23, 24], [21, 23], [22, 24], [25, 26], [27, 28], [25, 27], [26, 28], [29, 30], [31, 32], [29, 31], [30, 32], [33, 34], [35, 36], [33, 35], [34, 36], [37, 38], [39, 40], [37, 39], [38, 40]]

sweepline(V,EV)
=#

max_stars_repo_name: cvdlab-courses/cat

max_issues_repo_name: cvdlab-courses/cat

max_forks_repo_name: cvdlab-courses/cat


text:
export ActionTransformedEnv

struct ActionTransformedEnv{P,M,E<:AbstractEnv} <: AbstractEnvWrapper
    env::E
    action_mapping::M
    action_space_mapping::P
end

"""
    ActionTransformedEnv(env;action_space_mapping=identity, action_mapping=identity)

`action_space_mapping` will be applied to `action_space(env)` and
`legal_action_space(env)`. `action_mapping` will be applied to `action` before
feeding it into `env`.
"""
ActionTransformedEnv(env; action_mapping = identity, action_space_mapping = identity) =
    ActionTransformedEnv(env, action_mapping, action_space_mapping)

RLBase.action_space(env::ActionTransformedEnv, args...) =
    env.action_space_mapping(action_space(env.env), args...)

RLBase.legal_action_space(env::ActionTransformedEnv, args...) =
    env.action_space_mapping(legal_action_space(env.env), args...)

(env::ActionTransformedEnv)(action, args...; kwargs...) =
    env.env(env.action_mapping(action), args...; kwargs...)

ext: jl


text:
# test constructors
np = 10
nd = 6

vars = (x = 1:3, v = 4:6, z = 1:6)
params = (a = 1, b = 1.0)

pa = zeros(nd, np)
l1 = ParticleList(pa; variables = vars)
l2 = ParticleList(eltype(pa), np, nd; variables = vars)

@test l1.list == l2.list
@test l1.views == l2.views
@test l1.params == l2.params
@test l1.params === l2.params

# test assertions
@test_throws AssertionError ParticleList(pa; variables = (x = 1:3, v = 4:6), parameters = (a = 1, x = 2))
@test_throws AssertionError ParticleList(pa; variables = (x = 1:3, list = 4:6), parameters = (a = 1, b = 2))
@test_throws AssertionError ParticleList(pa; variables = (x = 1:3, v = 4:6), parameters = (a = 1, list = 2))

# test getproperty
pa = rand(nd, np)
pl = ParticleList(pa; variables = vars, parameters = params)

# test list matrix
@test pl.list === pa

# test views
@test pl.x == view(pl.list, 1:3, :)
@test pl.v == view(pl.list, 4:6, :)
@test pl.z == view(pl.list, 1:6, :)

# test particle vector
@test pl.particles == [Particle(view(pa, :, i); variables = vars, parameters = params) for i in axes(pa, 2)]

# test variable views
@test pl.variables.x == [view(pa, 1:3, i) for i in axes(pa,2)]
@test pl.variables.v == [view(pa, 4:6, i) for i in axes(pa,2)]
@test pl.variables.z == [view(pa, 1:6, i) for i in axes(pa,2)]

# test parameters
@test pl.a == params.a
@test pl.b == params.b

# test array interface
@test eltype(pl) == eltype(pa)
@test length(pl) == np
@test size(pl) == (nd, np)

# test getindex methods
@test pl[1] == pl.particles[1]
@test pl[1,1] == pl.list[1,1] == pa[1,1]
@test pl[1,:] == pl.list[1,:] == pa[1,:]
@test pl[:,1] == pl.list[:,1] == pa[:,1]

# test setindex methods
# .
# .
# .

# test iterate methods
@test eachindex(pl) == eachindex(pl.particles)

i = 0
for p in pl
    global i += 1
    @test p == pl[i]
end


text:
using Test
using NCDatasets

function example_file(i,array, fname = tempname();
    varname = "var")
    @debug "fname $fname"

    NCDataset(fname,"c") do ds
        # Dimensions

        ds.dim["lon"] = size(array,1)
        ds.dim["lat"] = size(array,2)
        ds.dim["time"] = Inf

        # Declare variables

        ncvar = defVar(ds,varname, Float64, ("lon", "lat", "time"))
        ncvar.attrib["field"] = "u-wind, scalar, series"
        ncvar.attrib["units"] = "meter second-1"
        ncvar.attrib["long_name"] = "surface u-wind component"
        ncvar.attrib["time"] = "time"
        ncvar.attrib["coordinates"] = "lon lat"

        nclat = defVar(ds,"lat", Float64, ("lat",))
        nclat.attrib["units"] = "degrees_north"

        nclon = defVar(ds,"lon", Float64, ("lon",))
        nclon.attrib["units"] = "degrees_east"
        nclon.attrib["modulo"] = 360.0

        nctime = defVar(ds,"time", Float64, ("time",))
        nctime.attrib["long_name"] = "surface wind time"
        nctime.attrib["field"] = "time, scalar, series"
        nctime.attrib["units"] = "days since 2000-01-01 00:00:00 GMT"

        # Global attributes

        ds.attrib["history"] = "foo"

        # Define variables

        g = defGroup(ds,"group")
        ncvarg = defVar(g,varname, Float64, ("lon", "lat", "time"))
        ncvarg.attrib["field"] = "u-wind, scalar, series"
        ncvarg.attrib["units"] = "meter second-1"
        ncvarg.attrib["long_name"] = "surface u-wind component"
        ncvarg.attrib["time"] = "time"
        ncvarg.attrib["coordinates"] = "lon lat"

        ncvar[:,:,1] = array
        ncvarg[:,:,1] = array.+1
        #nclon[:] = 1:size(array,1)
        #nclat[:] = 1:size(array,2)
        nctime[:] = i
    end
    return fname
end

A = [randn(2,3),randn(2,3),randn(2,3)]

C = cat(A...; dims = 3)
CA = CatArrays.CatArray(3,A...)

idx_global,idx_local,sz = CatArrays.idx_global_local_(CA,(1:1,1:1,1:1))

@inferred CatArrays.idx_global_local_(CA,(1:1,1:1,1:1))

@test CA[1:1,1:1,1:1] == C[1:1,1:1,1:1]
@test CA[1:1,1:1,1:2] == C[1:1,1:1,1:2]
@test CA[1:2,1:1,1:2] == C[1:2,1:1,1:2]

@test CA[:,:,1] == C[:,:,1]
@test CA[:,:,2] == C[:,:,2]
@test CA[:,2,:] == C[:,2,:]
@test CA[:,1:2:end,:] == C[:,1:2:end,:]
@test CA[1,1,1] == C[1,1,1]

CA[2,2,:] = [1.,2.,3.]
@test A[1][2,2] == 1.
@test A[2][2,2] == 2.
@test A[3][2,2] == 3.

A = [randn(2,3),randn(2,3),randn(2,3)]
C = cat(A...; dims = 3)
fnames = example_file.(1:3,A)

mfds = NCDataset(fnames, deferopen = false);
varname = "var"
var = variable(mfds,varname);
data = var[:,:,:]

@test C == var[:,:,:]

@test_throws BoundsError var[:,:,end+1]
@test_throws BoundsError mfds[varname].var[:,:,end+1]

@test mfds.attrib["history"] == "foo"
@test var.attrib["units"] == "meter second-1"

@test dimnames(var) == ("lon", "lat", "time")
# lon does not vary in time and thus there should be no aggregation
lon = variable(mfds,"lon");
@test lon.attrib["units"] == "degrees_east"
@test size(lon) == (size(data,1),)

var = mfds[varname]
@test C == var[:,:,:]
@test dimnames(var) == ("lon", "lat", "time")

@test mfds.dim["lon"] == size(C,1)
@test mfds.dim["lat"] == size(C,2)
@test mfds.dim["time"] == size(C,3)

@test mfds.dim["time"] == size(C,3)

# save a aggregated file
fname_merged = tempname()
write(fname_merged,mfds)

ds_merged = NCDataset(fname_merged)
@test mfds.dim["time"] == size(C,3)
close(ds_merged)

# show
buf = IOBuffer()
show(buf,mfds)
occursin("time = 3",String(take!(buf)))

close(mfds)

# write
mfds = NCDataset(fnames,"a",deferopen = false);
mfds[varname][2,2,:] = 1:length(fnames)

for n = 1:length(fnames)
    NCDataset(fnames[n]) do ds
        @test ds[varname][2,2,1] == n
    end
end

mfds.attrib["history"] = "foo2"
sync(mfds)

NCDataset(fnames[1]) do ds
    @test ds.attrib["history"] == "foo2"
end

@test_throws NCDatasets.NetCDFError NCDataset(fnames,"not-a-mode")

@test keys(mfds) == [varname, "lat", "lon", "time"]
@test keys(mfds.dim) == ["lon", "lat", "time"]
@test NCDatasets.groupname(mfds) == "/"
@test size(mfds[varname]) == (2, 3, 3)
@test size(mfds[varname].var) == (2, 3, 3)
@test name(mfds[varname].var) == varname
@test NCDatasets.groupname(mfds.group["group"]) == "group"

# create new dimension in all files
mfds.dim["newdim"] = 123;
sync(mfds);
NCDataset(fnames[1]) do ds
    @test ds.dim["newdim"] == 123
end
close(mfds)

# multi-file merge

ampl = rand(50,50)
vel = rand(50,50)

fnames = [example_file(1, vel; varname = "vel"),
          example_file(1, ampl; varname = "ampl")]

ds = NCDataset(fnames,aggdim = "", deferopen = false);

@test ds["ampl"][:,:,1] == ampl
@test ds["vel"][:,:,1] == vel

@test sort(keys(ds)) == ["ampl", "lat", "lon", "time", "vel"]

# save a merged file
fname_merged = tempname()
write(fname_merged,ds)
close(ds)

ds_merged = NCDataset(fname_merged)
@test sort(keys(ds_merged)) == ["ampl", "lat", "lon", "time", "vel"]
@test ds_merged["ampl"][:,:,1] == ampl
@test ds_merged["vel"][:,:,1] == vel
close(ds_merged)
rm(fname_merged)
nothing


text:
import MPI

export distribute, MPI_rank, MPI_size

if !MPI.Initialized()
    MPI.Init()
end

MPI_rank = MPI.Comm_rank(MPI.COMM_WORLD)
MPI_size = MPI.Comm_size(MPI.COMM_WORLD)
MPI_COMM_WORLD = MPI.COMM_WORLD

mpirank(comm::MPI.Comm) = MPI.Comm_rank(comm)
mpisize(comm::MPI.Comm) = MPI.Comm_size(comm)

function distribute(size::Int, comm_size::Int, rank::Int)
    """
    Compute the first element and size for a given rank

    Parameters
    ----------
    size : Int
        Length of array
    comm_size : Int
        Number of MPI  processes
    rank: Int
        MPI rank

    Returns
    -------
    size : Int
        Size
    offset : Int
        Offsets for the MPI process
    """
    if comm_size > size
        if rank + 1 <= size
            return rank + 1, 1
        else
            return 1, 0
        end
    end

    base = div(size, comm_size)
    leftover = size % comm_size
    if rank+1 <= leftover
        size = base+1
        start = rank * (base+1) + 1
    else
        size = base
        start = (base+1) * leftover + (rank-leftover) * base + 1
    end
    return start, size
end

function Allreduce(data, mpi_sum, comm)
    if comm === nothing
        return data
    end
    return MPI.Allreduce(data, mpi_sum, comm)
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
@with_kw struct OptimizationState{T} <: OptimizerTrace
    it::Int64
    time::Float64 = 0.0
    f_x::Float64
    g_x::Float64
    norm_step::Float64 = -1.0
    norm_minsubgradient_tangent::Float64 = -1.0
    norm_minsubgradient_normal::Float64 = -1.0
    ncalls_f::Int64 = -1
    ncalls_g::Int64 = -1
    ncalls_∇f::Int64 = -1
    ncalls_proxg::Int64 = -1
    ncalls_gradₘF::Int64 = -1
    ncalls_HessₘF::Int64 = -1
    ncalls_retr::Int64 = -1
    niter_CG::Int64 = -1
    niter_manls::Int64 = -1
    additionalinfo::T = NamedTuple()
end

const OptimizationTrace{T} = Vector{OptimizationState{T}}


text:
using TheAlgorithms
using Test

using LinearAlgebra

@testset "TheAlgorithms" begin
    include("basic.jl")
    include("conversions.jl")
    include("cipher.jl")
    include("data_structures.jl")
    include("knapsack.jl")
    include("longest_increasing_subsequence.jl")
    include("machine_learning.jl")
    include("math.jl")
    include("matrix.jl")
    include("project-rosalind.jl")
    include("scheduling.jl")
    include("searches.jl")
    include("sorts.jl")
    include("statistics.jl")
    include("strings.jl")
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
# Wilkinson-type polynomial defined by its roots:
# wₙ(x) = (x-1)⋅(x-2)⋅ ⋯ ⋅ (x-n)

using Polynomials

# function generate_wilkinson(n)#, T=BigFloat)   # SLOW

#     p = poly(collect(1:n))
#     p_prime = polyder(p)

#     coeffs = p.a #[convert(T, coeff) for coeff in p.a]
#     # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a]

#     function f(x)
#         total = coeffs[1]
#         for i in 2:length(coeffs)
#             total += coeffs[i]*x^(i-1)
#         end
#         total
#     end

#     f
# end

function generate_wilkinson_horner(n)
    p = poly(collect(1:n))  # make polynomial with roots 1 to n

    coeffs = p.a

    expr = :($(coeffs[end]))

    for i in length(coeffs)-1:-1:1
        expr = :($expr*x + $(coeffs[i]))
    end

    expr
end

function subscriptify(n::Int)
    subscript_digits = [c for c in "₀₁₂₃₄₅₆₇₈₉"]
    dig = reverse(digits(n))
    join([subscript_digits[i+1] for i in dig])
end

# Generate Wilkinson functions W₃ etc.:
for n in 1:15
    fn_name = symbol(string("W", subscriptify(n)))
    expr = generate_wilkinson_horner(n)

    @eval $(fn_name)(x) = $(expr)
end


text:
function nullspace(A::Operator{T};tolerance=10eps(real(T)),maxlength=1_000_000) where T
   K=transpose_nullspace(qr(A'),tolerance,maxlength)
    # drop extra rows, and use QR to determine rank
    Q,R = qr(K,Val(true))
    ind=findfirst(r->abs(r)≤100tolerance,diag(R))
    Kret=ind==0 ? Q : Q[:,1:ind-1]
    Fun(Space(fill(domainspace(A),(1,ind-1))),vec(Kret'))
end

function transpose_nullspace(QR::QROperator,tolerance,maxlength)
    T=eltype(QR)
    resizedata!(QR,:,100)

    m=size(QR.H,1)
    K=zeros(100,m-1)
    K[1:m-1,1:m-1]=Matrix(I,m-1,m-1)

    for k=1:10
        v=QR.H[:,k]
        QQ=Matrix{T}(I,m,m)-2v*v'
        K[:,:]=K*QQ[1:end-1,2:end]
        K[k+m-1,:]=QQ[end,2:end]
    end
    k=11
    α=0.9

    while slnorm(K,floor(Int,k^α),:) >tolerance*k
        if k > maxlength
            @warn "Max length of $maxlength reached."
            break
        end

        if k ≥ QR.ncols
            resizedata!(QR,:,k+100)
        end

        v=QR.H[:,k]
        QQ=(Matrix(I,m,m)-2v*v')
        K[:,:]=K*QQ[1:end-1,2:end]

        if k+m-1 > size(K,1)
            K=pad(K,k+m+100,:)
        end
        K[k+m-1,:]=QQ[end,2:end]
        k+=1
    end

    K[1:floor(Int,k^α),:]
end

nullspace(A::AbstractArray{OO};kwds...) where {OO<:Operator} = nullspace(interlace(A);kwds...)


text:
include("../rgi-setup.jl")
using BSON

function load_one(region, dir, dirbase="output/rasters-fit_sigma")
    base = dirname(@__FILE__)
    region_str = VAWTools.int2str2(region)
    dir_ = joinpath(base, "..", dirbase, dir)
    loaded1 = BSON.load(dir_*"/$(region_str)_out1.bson");
    loaded2 = BSON.load(dir_*"/$(region_str)_out2.bson");
    out1 = loaded1[:out1];
    out2 = loaded2[:out2];
    glnrs1, glnrs2 = try
        (loaded1[:glnrs1],
         loaded1[:glnrs2])
    catch
        (nothing, nothing)
    end
    # there was a mistake in glnrs1 for region=3
    if region==3
        for (i,o) in enumerate(out1)
            if (o isa Dict) && (o[:rgi]=="RGI60-03.00452" || o[:rgi]=="RGI60-03.00410")
                out1[i] = ErrorException("no h-data")
            end
        end
    end
    return out1, out2, glnrs1, glnrs2
end

function load_all(dir, dirbase="output/rasters-fit_sigma", regions=[3,4,7,11,14])
    out1s, out2s = Dict(), Dict()
    gnr1, gnr2 = Dict(), Dict()
    for r in regions
        print("Loading region $r...")
        try
            o1,o2,g1,g2 = load_one(r, dir, dirbase)
            out1s[r] = o1
            out2s[r] = o2
            gnr1[r] = g1
            gnr2[r] = g2
            println("done.")
        catch e
            println("error: $e.")
        end
    end
    return out1s, out2s, gnr1, gnr2
end

nothing


text:
module ITensorImpurity
  using ITensors

  export anderson, apply_layers

  function ITensors.op(o::ITensors.SiteOp, s::AbstractVector{<:Index})
    return op(o.name, s, o.site)
  end

  ITensors.sites(𝒽::ITensors.MPOTerm) = [only(o.site) for o in 𝒽.ops]
  ITensors.names(𝒽::ITensors.MPOTerm) = [o.name for o in 𝒽.ops]

  function ITensors.op(𝒽::ITensors.MPOTerm, s::AbstractVector{<:Index})
    n⃗ = ITensors.sites(𝒽)
    p = sortperm(n⃗)
    opnames = ITensors.names(𝒽)
    s⃗ = s[n⃗][p]
    𝒽′ = ITensors.MPOTerm(𝒽.coef, [ITensors.SiteOp(opnames[n], p[n]) for n in eachindex(𝒽.ops)])
    return prod(MPO(ITensors.OpSum([𝒽′]), s⃗))
  end

  # Trotter layer
  function Base.exp(::Val{:trotter_order_1}, δt::Number, ℋ::OpSum, s::AbstractVector{<:Index})
    # TODO: Sort and merge terms on the same sites.
    circuit = Vector{ITensor}(undef, length(ℋ))
    for n in 1:length(ℋ)
      # TODO: Use `cis` in the imaginary `δt` case.
      circuit[n] = exp(δt * op(ℋ[n], s))
    end
    return circuit
  end

  function Base.exp(::Val{:trotter_order_2}, δt::Number, ℋ::OpSum, s::Vector{<:Index})
    expH1 = exp(Val{:trotter_order_1}(), δt / 2, ℋ, s)
    return vcat(expH1, reverse(expH1))
  end

  function contract_slice(t1::ITensor, t2::ITensor; slice_ind)
    @assert hasind(t1, slice_ind)
    @assert hasind(t2, slice_ind)
    i = slice_ind
    r = ITensor()
    for n in 1:dim(i)
      pn = onehot(i => n)
      r += (t1 * pn) * (t2 * pn)
    end
    return r
  end

  function contract_slice(::typeof(linkinds), tn)
    ∏tn = tn[1]
    for n in 2:length(tn)
      l = commonind(∏tn, tn[n])
      ∏tn = contract_slice(∏tn, tn[n]; slice_ind=l)
    end
    return ∏tn
  end

  const σ⃗ = ("↑", "↓")

  function anderson(N; Nᴸ, t, V, t′, U, Vᵍ)
    ℋ = OpSum()
    for n in 1:Nᴸ
      if n < Nᴸ
        for σ in σ⃗
          ℋ .-= t, "c†$σ", n, "c$σ", n + 1
          ℋ .-= t, "c†$σ", n + 1, "c$σ", n
        end
      end
      ℋ .-= V / 2, "ntot", n
    end
    for σ in σ⃗
      ℋ .-= t′, "c†$σ", Nᴸ, "c$σ", Nᴸ + 1
      ℋ .-= t′, "c†$σ", Nᴸ + 1, "c$σ", Nᴸ
      ℋ .-= t′, "c†$σ", Nᴸ + 1, "c$σ", Nᴸ + 2
      ℋ .-= t′, "c†$σ", Nᴸ + 2, "c$σ", Nᴸ + 1
    end
    ℋ .+= U, "n↑↓", Nᴸ + 1
    ℋ .+= Vᵍ, "ntot", Nᴸ + 1
    for n in (Nᴸ + 2):N
      if n < N
        for σ in σ⃗
          ℋ .-= t, "c†$σ", n, "c$σ", n + 1
          ℋ .-= t, "c†$σ", n + 1, "c$σ", n
        end
      end
      ℋ .+= V / 2, "ntot", n
    end
    return ℋ
  end

  function J(ψᵗ, t′; Nᴸ)
    J⃗ᴸᵗ = Dict()
    J⃗ᴿᵗ = Dict()
    for σ in σ⃗
      J⃗ᴸᵗ[σ] = im * t′ * correlation_matrix(ψᵗ, "c†$σ", "c$σ"; site_range=Nᴸ:(Nᴸ + 1))
      J⃗ᴿᵗ[σ] = im * t′ * correlation_matrix(ψᵗ, "c†$σ", "c$σ"; site_range=(Nᴸ + 1):(Nᴸ + 2))
    end
    Jᴸᵗ = sum(σ -> J⃗ᴸᵗ[σ][1, 2] - J⃗ᴸᵗ[σ][2, 1], σ⃗)
    Jᴿᵗ = sum(σ -> J⃗ᴿᵗ[σ][1, 2] - J⃗ᴿᵗ[σ][2, 1], σ⃗)
    return (Jᴸᵗ + Jᴿᵗ) / 2
  end

  function cc(ψ; Nᴸ, t′)
    ccs = Dict()
    for σ in σ⃗
      ccs[σ] = im * t′ * correlation_matrix(ψ, "c†$σ", "c$σ"; site_range=Nᴸ:(Nᴸ + 2))
    end
    return ccs
  end

  function apply_layers(U⃗, ψ0; f, kwargs...)
    nlayers = length(U⃗)
    f⃗ = Vector{Any}(undef, nlayers + 1)
    ψτ = ψ0
    f⃗[1] = f(ψτ)
    for n in 1:nlayers
      @show n, nlayers
      ψτ = apply(U⃗[n], ψτ; kwargs...)
      @show maxlinkdim(ψτ)
      @show fτ = f(ψτ)
      f⃗[n + 1] = fτ
    end
    return ψτ, f⃗
  end
end # module ITensorImpurity

max_stars_repo_name: mtfishman/ITensorImpurity

max_issues_repo_name: mtfishman/ITensorImpurity

max_forks_repo_name: mtfishman/ITensorImpurity


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Seven, The Treachery of Whales."

include("Utils.jl")

using DelimitedFiles
using Statistics
using .Utils

"Calculate the minimum amount of fuel used to align the `crabs`’ position."
function solve1(crabs :: Array{Int64}) :: Int64
    median_crab = crabs |> median |> floor

    return crabs .- median_crab .|> abs |> sum
end

#= Initially I only used the floor function and it worked properly, but as an
 = user in Reddit pointed out, floor() isn’t always the right rounding.
 =#

"""Calculate the minimum amount of fuel used to align the `crabs`’ position,
with an arithmetic series usage of fuel."""
function solve2(crabs :: Array{Int64}) :: Int64
    mean_crab = crabs |> mean |> floor |> Int64

    Σ₁ = crabs .- mean_crab .|> abs .|> ∑
    Σ₂ = crabs .- (mean_crab + 1) .|> abs .|> ∑

    return [Σ₁, Σ₂] .|> sum |> minimum
end

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = readdlm(first(ARGS), ',', Int)

    S₁, S₂ = input_data |> solve1, input_data |> solve2

    @assert S₁ == 343441
    @assert S₂ == 98925151

    S₁ |> println
    S₂ |> println
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    FlipSignBridge{T, S1, S2, F, G}

Bridge a `G`-in-`S1` constraint into an `F`-in-`S2` constraint by multiplying
the function by `-1` and taking the point reflection of the set across the
origin. The flipped `F`-in-`S` constraint is stored in the `constraint`
field by convention.
"""
abstract type FlipSignBridge{
    T,
    S1<:MOI.AbstractSet,
    S2<:MOI.AbstractSet,
    F<:MOI.AbstractFunction,
    G<:MOI.AbstractFunction,
} <: SetMapBridge{T,S2,S1,F,G} end

function MOI.delete(
    model::MOI.ModelLike,
    bridge::FlipSignBridge,
    i::MOIB.IndexInVector,
)
    func = MOI.get(model, MOI.ConstraintFunction(), bridge.constraint)
    idx = setdiff(1:MOI.output_dimension(func), i.value)
    new_func = MOIU.eachscalar(func)[idx]
    set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
    new_set = MOI.update_dimension(set, MOI.dimension(set) - 1)
    MOI.delete(model, bridge.constraint)
    bridge.constraint = MOI.add_constraint(model, new_func, new_set)
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::FlipSignBridge,
    change::MOI.ScalarCoefficientChange,
)
    MOI.modify(
        model,
        bridge.constraint,
        MOI.ScalarCoefficientChange(change.variable, -change.new_coefficient),
    )
    return
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::FlipSignBridge,
    change::MOI.MultirowChange{T},
) where {T}
    new_coefficients = Tuple{Int64,T}[
        (index, -coef) for (index, coef) in change.new_coefficients
    ]
    MOI.modify(
        model,
        bridge.constraint,
        MOI.MultirowChange(change.variable, new_coefficients),
    )
    return
end

"""
    GreaterToLessBridge{
        T,
        F<:MOI.AbstractScalarFunction,
        G<:MOI.AbstractScalarFunction
    } <: FlipSignBridge{T, MOI.GreaterThan{T}, MOI.LessThan{T}, F, G}

Transforms a `G`-in-`GreaterThan{T}` constraint into an `F`-in-`LessThan{T}`
constraint.
"""
struct GreaterToLessBridge{
    T,
    F<:MOI.AbstractScalarFunction,
    G<:MOI.AbstractScalarFunction,
} <: FlipSignBridge{T,MOI.GreaterThan{T},MOI.LessThan{T},F,G}
    constraint::MOI.ConstraintIndex{F,MOI.LessThan{T}}
end

const GreaterToLess{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToLessBridge{T},OT}

function MOIB.map_set(::Type{<:GreaterToLessBridge}, set::MOI.GreaterThan)
    return MOI.LessThan(-set.lower)
end

function MOIB.inverse_map_set(::Type{<:GreaterToLessBridge}, set::MOI.LessThan)
    return MOI.GreaterThan(-set.upper)
end

function concrete_bridge_type(
    ::Type{<:GreaterToLessBridge{T}},
    G::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.GreaterThan{T}},
) where {T}
    F = MOIU.promote_operation(-, T, G)
    return GreaterToLessBridge{T,F,G}
end

"""
    LessToGreaterBridge{
        T,
        F<:MOI.AbstractScalarFunction,
        G<:MOI.AbstractScalarFunction
    } <: FlipSignBridge{T, MOI.LessThan{T}, MOI.GreaterThan{T}, F, G}

Transforms a `G`-in-`LessThan{T}` constraint into an `F`-in-`GreaterThan{T}`
constraint.
"""
struct LessToGreaterBridge{
    T,
    F<:MOI.AbstractScalarFunction,
    G<:MOI.AbstractScalarFunction,
} <: FlipSignBridge{T,MOI.LessThan{T},MOI.GreaterThan{T},F,G}
    constraint::MOI.ConstraintIndex{F,MOI.GreaterThan{T}}
end

const LessToGreater{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToGreaterBridge{T},OT}

function MOIB.map_set(::Type{<:LessToGreaterBridge}, set::MOI.LessThan)
    return MOI.GreaterThan(-set.upper)
end

function MOIB.inverse_map_set(
    ::Type{<:LessToGreaterBridge},
    set::MOI.GreaterThan,
)
    return MOI.LessThan(-set.lower)
end

function concrete_bridge_type(
    ::Type{<:LessToGreaterBridge{T}},
    G::Type{<:MOI.AbstractScalarFunction},
    ::Type{MOI.LessThan{T}},
) where {T}
    F = MOIU.promote_operation(-, T, G)
    return LessToGreaterBridge{T,F,G}
end

"""
    NonnegToNonposBridge{
        T,
        F<:MOI.AbstractVectorFunction,
        G<:MOI.AbstractVectorFunction
    } <: FlipSignBridge{T, MOI.Nonnegatives, MOI.Nonpositives, F, G}

Transforms a `G`-in-`Nonnegatives` constraint into a `F`-in-`Nonpositives`
constraint.
"""
mutable struct NonnegToNonposBridge{
    T,
    F<:MOI.AbstractVectorFunction,
    G<:MOI.AbstractVectorFunction,
} <: FlipSignBridge{T,MOI.Nonnegatives,MOI.Nonpositives,F,G}
    constraint::MOI.ConstraintIndex{F,MOI.Nonpositives}
end

const NonnegToNonpos{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonnegToNonposBridge{T},OT}

function concrete_bridge_type(
    ::Type{<:NonnegToNonposBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.Nonnegatives},
) where {T}
    F = MOIU.promote_operation(-, T, G)
    return NonnegToNonposBridge{T,F,G}
end

"""
    NonposToNonnegBridge{
        T,
        F<:MOI.AbstractVectorFunction,
        G<:MOI.AbstractVectorFunction,
    } <: FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives, F, G}

Transforms a `G`-in-`Nonpositives` constraint into a `F`-in-`Nonnegatives`
constraint.
"""
mutable struct NonposToNonnegBridge{
    T,
    F<:MOI.AbstractVectorFunction,
    G<:MOI.AbstractVectorFunction,
} <: FlipSignBridge{T,MOI.Nonpositives,MOI.Nonnegatives,F,G}
    constraint::MOI.ConstraintIndex{F,MOI.Nonnegatives}
end

const NonposToNonneg{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonposToNonnegBridge{T},OT}

function MOIB.map_set(::Type{<:NonposToNonnegBridge}, set::MOI.Nonpositives)
    return MOI.Nonnegatives(set.dimension)
end

function MOIB.inverse_map_set(
    ::Type{<:NonposToNonnegBridge},
    set::MOI.Nonnegatives,
)
    return MOI.Nonpositives(set.dimension)
end

function concrete_bridge_type(
    ::Type{<:NonposToNonnegBridge{T}},
    G::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.Nonpositives},
) where {T}
    F = MOIU.promote_operation(-, T, G)
    return NonposToNonnegBridge{T,F,G}
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
#=
LDA.jl

Adham Beyki, odinay@gmail.com
27/10/2015
=#

##########################################
###### Latent Dirichlet Allocation ######
##########################################
immutable LDA{T} <: hMixtureModel
    component::T
    KK::Int
    aa::Float64

    LDA{T}(c::T, KK::Int, aa::Float64) = new(c, KK, aa)
end
LDA{T}(c::T, KK::Int, aa::Real) = LDA{typeof(c)}(c, KK, convert(Float64, aa))

function Base.show(io::IO, lda::LDA)
    println(io, "Latent Dirichlet Allocation model with $(lda.KK) $(typeof(lda.component)) components")
end

function storesample(zz::Vector{Vector{Int}}, sample_n::Int, filename::ASCIIString)

    n_groups = length(zz)
    n_group_j = [length(zz[jj]) for jj=1:n_groups]
    println(typeof(n_group_j))
    zz_flat = Int[]
    for jj = 1:n_groups
        append!(zz_flat, zz[jj])
    end

    if endswith(filename, "_")
        dummy_filename = string(filename, sample_n, ".h5")
    else
        dummy_filename = string(filename, "_", sample_n, ".h5")
    end

    println("storing on disk...")
    HDF5.h5open(dummy_filename, "w") do file
        HDF5.write(file, "n_groups", n_groups)
        HDF5.write(file, "n_group_j", n_group_j)
        HDF5.write(file, "zz_flat", zz_flat)
        HDF5.write(file, "sample_n", sample_n)
    end
end

function collapsed_gibbs_sampler{T1, T2}(
    lda::LDA{T1},
    xx::Vector{Vector{T2}},
    zz::Vector{Vector{Int}},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    store_every::Int=250, results_path="", filename::ASCIIString="LDA_results_")

    components = Array(typeof(lda.component), lda.KK)
    for kk = 1:lda.KK
        components[kk] = deepcopy(lda.component)
    end

    n_groups       = length(xx)
    n_group_j      = [length(zz[jj]) for jj = 1:n_groups]
    n_iterations   = n_burnins + (n_samples)*(n_lags+1)
    lda_aa         = fill(lda.aa, lda.KK)
    pp             = zeros(Float64, lda.KK)
    nn             = zeros(Int, n_groups, lda.KK)
    log_likelihood = 0.0

    # initializing the components
    tic()
    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
            log_likelihood += loglikelihood(components[kk], xx[jj][ii])
        end
    end
    elapsed_time = toq()

    # starting the MCMC chain
    for iteration = 1:n_iterations

        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, time=%.2f, likelihood=%.2f", iteration, lda.KK, elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()
        @inbounds for jj = randperm(n_groups)
            @inbounds for ii = randperm(n_group_j[jj])

                # 1
                # remove the datapoint
                kk = zz[jj][ii]
                delitem!(components[kk], xx[jj][ii])
                nn[jj, kk] -= 1

                # 2
                # sample zz
                @inbounds for kk = 1:lda.KK
                    pp[kk] = log(nn[jj, kk] + lda.aa) + logpredictive(components[kk], xx[jj][ii])
                end

                lognormalize!(pp)
                kk = sample(pp)

                # 3
                # add the datapoint to the newly sampled cluster
                zz[jj][ii] = kk
                additem!(components[kk], xx[jj][ii])
                nn[jj, kk] += 1
                log_likelihood += loglikelihood(components[kk], xx[jj][ii])
            end #ii
        end # jj
        elapsed_time = toq()

        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(zz, sample_n, normalized_filename)
            end
        end
    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(zz, sample_n, normalized_filename)
end

function posterior{T1, T2}(lda::LDA{T1}, xx::Vector{Vector{T2}}, zz::Vector{Vector{Int}})

    components = Array(typeof(lda.component), lda.KK)
    for kk = 1:lda.KK
        components[kk] = deepcopy(lda.component)
    end

    n_groups = length(xx)
    n_group_j = [length(zz[jj]) for jj = 1:n_groups]
    nn = zeros(Int, n_groups, lda.KK)

    for jj = 1:n_groups
        for ii = 1:n_group_j[jj]
            kk = zz[jj][ii]
            additem!(components[kk], xx[jj][ii])
            nn[jj, kk] += 1
        end
    end

    return([posterior(components[kk]) for kk =1:lda.KK], nn)
end


text:
using Plots
theme(:juno)
theme(:default)

push!(LOAD_PATH,"src")

using DalitzPlotAnalysis
using amplitudes_compass
using PWAHelper
using SDMHelper
using FittingPWALikelihood

mass_bin_name = "2320_2340"

for app in ["rd", "mc", "fu"]
    @eval $(Symbol("kinvar_"*app)) = joinpath("data","$(mass_bin_name)_t1_"*$app*".txt")
    @eval $(Symbol("basisfunc_"*app)) = joinpath("data","$(mass_bin_name)_t1_"*$app*".jld")
end

##### Converting data data #####

iffile("data/$(mass_bin_name)_rd.root") || error("Cannot find the file")
iffile("data/$(mass_bin_name)_mc.root") || error("Cannot find the file")
@time let tmin = "0.1", tmax = "0.112853"
    run(pipeline(`src/Extract data/$(mass_bin_name)_rd.root D $(tmin) $(tmax)`, stdout=kinvar_rd))
    run(pipeline(`src/Extract data/$(mass_bin_name)_mc.root M $(tmin) $(tmax)`, stdout=kinvar_mc))
    run(pipeline(`src/Extract data/$(mass_bin_name)_mc.root F $(tmin) $(tmax)`, stdout=kinvar_fu))
end

################################

@time precalculate_compass_basis(kinvar_rd, basisfunc_rd)
@time precalculate_compass_basis(kinvar_mc, basisfunc_mc)
@time precalculate_compass_basis(kinvar_fu, basisfunc_fu)

@time const PsiMC = read_precalc_basis(basisfunc_mc);

@time const sum_mat = [sum(PsiMC[e,i]'*PsiMC[e,j] for e in 1:size(PsiMC,1))
    for i=1:Nwaves, j=1:Nwaves] /size(PsiMC,1);
# 10.593515 seconds (101.00 k allocations: 4.752 MiB)

let sum_mat_n = [sum_mat[i,j]/sqrt(sum_mat[i,i]*sum_mat[j,j]) for i=1:Nwaves, j=1:Nwaves];
    heatmap(real(sum_mat_n))
end

const noϵ = [i==1 for i=1:size(wavesfile[:,6],1)]
const posϵ = [ϵ=="+" for ϵ in wavesfile[:,6]]
const negϵ = [ϵ=="-" for ϵ in wavesfile[:,6]]

const ModelBlocks = [noϵ, posϵ]
const Npar = size(get_parameter_map(ModelBlocks),2)

const PsiDT = read_precalc_basis(basisfunc_rd);

LLH, GRAD, LLH_and_GRAD!, HES = createLLHandGRAD(PsiDT, sum_mat, ModelBlocks);

test_t = rand(Npar)-0.5
@time @show LLH(test_t)
@time minpars = minimize(LLH, LLH_and_GRAD!; verbose=1, starting_pars=test_t)

minpars0 = vcat(readdlm("minpars_r1.txt")...);
@time minpars = minimize(LLH, LLH_and_GRAD!; verbose=1, starting_pars=minpars0)
writedlm("minpars_r1.txt", minpars)

@time minpars_tries = let Natt = 5
    tries = [];
    for i in 1:Natt
        _minpars = minimize(LLH, LLH_and_GRAD!; verbose=1, starting_pars=rand(Npar))
        _llh = LLH(_minpars)
        push!(tries,(_llh,_minpars))
    end
    tries
end

# three well pronounced minimas
@show [v[1] for v in minpars_tries]

minpars_first = minpars_tries[3][2]
writedlm("minpars_r2.txt", minpars_first)
minpars_second = minpars_tries[1][2]
writedlm("minpars_r2.txt", minpars_second)
minpars_third = minpars_tries[5][2]
writedlm("minpars_r3.txt", minpars_third)

#####################################################################################
#####################################################################################
@time hes_anal = HES(minpars)
inv_hes = inv(hes_anal)

diag_error = sqrt.(diag(inv_hes))
# diag(inv_hes);
let tog = [[minpars[i],diag_error[i]] for i in 1:length(minpars)]
    stog = sort(tog, by=x->abs(x[1]), rev=true)
    hcat_stog = hcat(stog...)
    plot(abs.(hcat_stog[1,:]), yerr = hcat_stog[2,:], ylim=(0,1))
end

#####################################################################################
#####################################################################################

@time weights = let mpars = minpars,
    Tmap = get_parameter_map(ModelBlocks),
    pblocks = make_pblock_masks(ModelBlocks)
        [cohsq(extnd(PsiMC[i,:],Tmap).*mpars,pblocks) for i in 1:size(PsiMC,1)];
end;

const MC1 = readdlm(kinvar_mc);
@time const MC3 = hcat([swap_kin_parameters(MC1[e,:]...) for e in 1:size(MC1,1)]...)';

histogram(sqrt.(vcat(MC1[:,2],MC3[:,2])), weights=vcat(weights,weights), bins=(linspace(0.3,2.2,100)))

const DT1 = readdlm(kinvar_rd);
@time const DT3 = hcat([swap_kin_parameters(DT1[e,:]...) for e in 1:size(DT1,1)]...)';

histogram2d(DT1[:,2], DT3[:,2], bins=linspace(0,6,100))

stephist(sqrt.(vcat(MC1[:,2],MC3[:,2])), weights=vcat(weights,weights)/size(MC1,1)*size(DT1,1),
    bins=linspace(0.3,2.2,100),lab="weighted MC", size=(800,500), xlab="M3pi (GeV)")
let v = stephist!(sqrt.(vcat(DT1[:,2],DT3[:,2])), bins=linspace(0.3,2.2,100),lab="data")
    yarr = v.subplots[1][2][:y]
    xarr = v.subplots[1][2][:x]
    @show length(xarr)
    x2 = [(xarr[2i]+xarr[2(i+1)])/2 for i in 1:(div(length(xarr),2)-1)];
    y2 = [yarr[2i] for i in 1:(div(length(yarr),2)-1)];
    plot!(x2, y2, err=sqrt.(y2), lab="", marker = (0, stroke(0.5, :gray)), l=nothing)
end
# savefig(joinpath("plots","data_with_errors.png"))
# savefig(joinpath("plots","data_with_errors.pdf"))

#####################################################################################
#####################################################################################

@time const PsiFU = read_precalc_basis(basisfunc_fu);

@time const BmatFU = [sum(PsiFU[e,i]'*PsiFU[e,j] for e in 1:size(PsiFU,1))
    for i=1:Nwaves, j=1:Nwaves] /size(PsiMC,1);
# BmatFU
SDM = size(PsiDT,1)*pars_to_SDM(minpars, BmatFU, ModelBlocks)

minpars./SDM_to_pars(SDM/size(PsiDT,1), BmatFU, ModelBlocks)

SDM_RD = let path_to_SDM = "/localhome/mikhasenko/cernbox/tmp/pwa_from_scratch_data"
    readdlm(joinpath(path_to_SDM,"0.100000-0.112853","sdm92.re")) +
    1im*readdlm(joinpath(path_to_SDM,"0.100000-0.112853","sdm92.im"))
end

SDM_RD_err = let path_to_SDM = "/localhome/mikhasenko/cernbox/tmp/pwa_from_scratch_data"
    readdlm(joinpath(path_to_SDM,"0.100000-0.112853","sdm92-err.re")) +
    1im*readdlm(joinpath(path_to_SDM,"0.100000-0.112853","sdm92-err.im"))
end

minpars_FHDR = SDM_to_pars(SDM_RD/size(PsiDT,1), BmatFU, ModelBlocks)

plot([real.(diag(SDM_RD)) real.(diag(SDM))], lab=["F.H.-D.R." "M.M."],
    xlab="# wave", title = "Diagonal of the SDM", size=(800,500))
savefig(joinpath("plots","sdm_results_$(mass_bin_name)_t1.png"))
savefig(joinpath("plots","sdm_results_$(mass_bin_name)_t1.pdf"))
plot(-real.(diag(SDM))./real.(diag(SDM_RD)) .+ 1, lab=["1 - M.M. / F.H.-D.R."],
    xlab="# wave", title = "Diagonal of the SDM", size=(1800,500),
    xticks=[i for i in 1:88])
savefig(joinpath("plots","sdm_diff_$(mass_bin_name)_t1.pdf"))

let rd=real.(diag(SDM_RD)), wn = wavenames
    sort([(wn[i],rd[i]) for i in 1:length(rd)], by=x->x[2],rev=true)
end
let
    plot(real.(SDM_RD[2,:]), frame=:origin)
    plot!(real.(SDM[2,:]))
end
plot(hcat(real.(SDM_RD[2,:]), real.(SDM[2,:])), lab=["F.H.-D.R." "M.M."],
    xlab="# wave", title = "Real part of the second row of SDM", size=(800,500))
savefig(joinpath("plots","sdm2_results.png"))
savefig(joinpath("plots","sdm2_results.pdf"))

sum(diag(SDM))
sum(diag(SDM_RD))

diag(SDM)

minpars_rd = SDM_to_pars(SDM_RD/size(PsiDT,1), BmatFU, ModelBlocks)
minpars_cv = SDM_to_pars(SDM   /size(PsiDT,1), BmatFU, ModelBlocks)
@show [minpars_cv minpars_rd]
minpars_rd2 = [copysign(v,minpars[i]) for (i,v) in enumerate(minpars_rd)]
minpars

arg(z) = atan2(imag(z), real(z))
SDM_with_sign = diag(SDM).*cis.([arg(SDM[2,i]) for i in 1:size(SDM,1)]);

gr()
let v = plot(size=(500,400))
    cutoff = 0.2*max(abs.(SDM_with_sign)...)
    for (i,v) in enumerate(SDM_with_sign)
        plot!([0, real(v)], [0, imag(v)], lab="", l=(1))
        (abs(v) > cutoff) && annotate!(
            [(0.75real(v), 0.75imag(v), text(wavenames[i][3:end],10))])
    end
    plot!()
end
savefig("/tmp/arrow_intensities.pdf")
writedlm("/tmp/SDM.2300.re",real(SDM))
writedlm("/tmp/SDM.2300.im",imag(SDM))

#####################################################################################
#####################################################################################

writedlm(joinpath("data","bootstrap_tmp","minpars_r1.txt"), minpars)
writedlm(joinpath("data","bootstrap_tmp","minpars_r2.txt"), minpars_second)
writedlm(joinpath("data","bootstrap_tmp","minpars_r3.txt"), minpars_third)
BootstrapResults = let Nb = 100
    # path to save
    path_to_tmp_res = joinpath("data","bootstrap_tmp")
    # path to save
    minpars0 = vcat(readdlm(joinpath(path_to_tmp_res,"minpars_r3.txt"))...);
    res = Matrix{Float64}(Nb,length(minpars0))
    llh = Vector{Float64}(Nb)
    @progress for b in 15:Nb
        Nd = size(PsiDT,1)
        const _PsiDT = hcat([PsiDT[rand(1:Nd),:] for i in 1:Nd]...).'
        # @show size(_PsiDT)
        _LLH, _GRAD, _LLH_and_GRAD!, _HES = createLLHandGRAD(_PsiDT, sum_mat, ModelBlocks);
        @time _minpars = minimize(_LLH, _LLH_and_GRAD!; verbose=1, starting_pars=minpars0)
        # @show minpars
        res[b,:] = _minpars
        llh[b] = _LLH(_minpars)
        writedlm(joinpath(path_to_tmp_res,"BootstrapResults-$(b).txt"), res[b,:])
        writedlm(joinpath(path_to_tmp_res,"llh-$(b).txt"), llh[b])
    end
    res
end
# BootstrapResults = let Nb = 100, path_to_tmp_res = joinpath("data","bootstrap_tmp")
#     res = Matrix{Float64}(Nb,Npar)
#     for b in 1:Nb
#         res[b,:] = readdlm(joinpath(path_to_tmp_res,"BootstrapResults-$(b).txt"))
#     end
#     res
# end
# writedlm(joinpath("data","bootstrap_tmp","bootstrap_r3.txt"), BootstrapResults)
# llh = let Nb = 100, path_to_tmp_res = joinpath("data","bootstrap_tmp")
#     res = Vector{Float64}(Nb)
#     for b in 1:Nb
#         res[b] = readdlm(joinpath(path_to_tmp_res,"llh-$(b).txt"))[1]
#     end
#     res
# end
# writedlm(joinpath("data","bootstrap_tmp","llh_r3.txt"), llh)
# histogram!(llh, bins=30, lab = "-LogLH")
# vline!([LLH(minpars_third)], lab = "main fit")
# savefig("plots/llh_three_minima.pdf")

vline!([minpars[1]])
histogram(BootstrapResults[:,1])
scatter(BootstrapResults[:,103], BootstrapResults[:,102])

btp_error = sqrt.([cov(BootstrapResults[:,i]) for i in 1:Npar])
btp_mean = [mean(BootstrapResults[:,i]) for i in 1:Npar]

plot([minpars btp_mean], lab=["main fit" "bootstrap mean"], xlab = "# parameter")

let m = max(diag_error...)
	for (i,v) in enumerate(diag_error)
		(v == m) && @show i
	end
end

let tog = [[minpars[i],diag_error[i]] for i in 1:length(minpars)]
    stog = sort(tog, by=x->abs(x[1]), rev=true)
    hcat_stog = hcat(stog...)
    plot(abs.(hcat_stog[1,:]), yerr = hcat_stog[2,:], ylim=(0,1),
    lab="", title="PWA results and the errors Hessian", xlab="# parameter", size=(800,500))
end
savefig(joinpath("plots","parameters_with_hessian_errors.png"))
savefig(joinpath("plots","parameters_with_hessian_errors.pdf"))

let tog = [[btp_mean[i], btp_error[i]] for i in 1:length(minpars)]
    stog = sort(tog, by=x->abs(x[1]), rev=true)
    hcat_stog = hcat(stog...)
    plot!(abs.(hcat_stog[1,:]), yerr = hcat_stog[2,:], ylim=(0,1),
    lab="", title="PWA results and the errors Bootstrap", xlab="# parameter", size=(800,500))
end
savefig(joinpath("plots","parameters_with_bootstrap_errors.png"))
savefig(joinpath("plots","parameters_with_bootstrap_errors.pdf"))

###
SDM_BootstrapResults = hcat(
	[real.(diag(size(PsiDT,1)*pars_to_SDM(BootstrapResults[i,:], BmatFU, ModelBlocks)))
		for i in 1:size(BootstrapResults,1)]...);

plot(SDM_BootstrapResults, lab="")

let v = SDM_BootstrapResults
	p = sort([v[i,:] for i in 1:size(v,1)], by=x->mean(x), rev=true)
	quantiles = [quantile(p[i],[0.16,0.84]) for i in 1:length(p)]
	plot([quantiles[i][1] for i in 1:length(p)], fillrange=[quantiles[i][2] for i in 1:length(p)], lab="16%-84% quantiles", c=:lightgrey)
	plot!([mean(p[i]) for i in 1:length(p)], lab="Bootstrap mean")
	pm = sort(real.(diag(SDM)), rev=true)
	plot!(pm, lab="main fit", size=(1000,400))
end
let v = real.(diag(SDM_RD)), ve = real.(diag(SDM_RD_err))
	vs = sort([[v[i],ve[i]] for i in 1:length(v)], by=x->x[1], rev=true);
	plot!([i[1] for i in vs], yerr = [i[2] for i in vs], lab = "Florain-Dima")
end
savefig(joinpath("plots","compare_everything.pdf"))

let i = 8
	histogram(SDM_BootstrapResults[i,:], bins=20)
	vline!([(real(SDM[i,i]) + sqrt(cov(SDM_BootstrapResults[i,:]))*[-1,0,1])...])
end

max(abs.(minpars_second./minpars_first)...)

vcat([1,3,4]...,[4])

let v = BootstrapResults.'
	p0 = sort([vcat(v[i,:]..., minpars_first[i], minpars_second[i]) for i in 1:size(v,1)], by=x->mean(x), rev=true)
	# p = [v[i,:] for i in 1:size(v,1)]
	p = [p0[i][1:end-2] for i in 1:length(p0)]
	quantiles = [quantile(p[i][1:end-2],[0.16,0.84]) for i in 1:length(p)]
	plot([quantiles[i][1] for i in 1:length(p)], fillrange=[quantiles[i][2] for i in 1:length(p)], lab="16%-84% quantiles", c=:lightgrey)
	plot!([mean(p[i]) for i in 1:length(p)], lab="Bootstrap mean", size=(2000,800))
	plot!([p0[i][end-1] for i in 1:length(p0)], lab="deepest min")
	plot!([p0[i][end] for i in 1:length(p0)], lab="second min")
end

savefig(joinpath("plots","minimias.pdf"))


text:
"""
MCTS solver with DPW

Fields:

    depth::Int64:
        Maximum rollout horizon and tree depth.
        default: 10

    exploration_constant::Float64:
        Specified how much the solver should explore.
        In the UCB equation, Q + c*sqrt(log(t/N)), c is the exploration constant.
        default: 1.0

    n_iterations::Int64
        Number of iterations during each action() call.
        default: 100

    max_time::Float64
        Maximum amount of CPU time spent iterating through simulations.
        default: Inf

    k_action::Float64
    alpha_action::Float64
    k_state::Float64
    alpha_state::Float64
        These constants control the double progressive widening. A new state
        or action will be added if the number of children is less than or equal to kN^alpha.
        defaults: k:10, alpha:0.5

    keep_tree::Bool
        If true, store the tree in the planner for reuse at the next timestep (and every time it is used in the future). There is a computational cost for maintaining the state dictionary necessary for this.
        default: false

    enable_action_pw::Bool
        If true, enable progressive widening on the action space; if false just use the whole action space.
        default: true

    check_repeat_state::Bool
    check_repeat_action::Bool
        When constructing the tree, check whether a state or action has been seen before (there is a computational cost to maintaining the dictionaries necessary for this)
        default: true

    tree_in_info::Bool:
        If true, return the tree in the info dict when action_info is called. False by default because it can use a lot of memory if histories are being saved.
        default: false

    rng::AbstractRNG:
        Random number generator

    estimate_value::Any (rollout policy)
        Function, object, or number used to estimate the value at the leaf nodes.
        If this is a function `f`, `f(mdp, s, depth)` will be called to estimate the value.
        If this is an object `o`, `estimate_value(o, mdp, s, depth)` will be called.
        If this is a number, the value will be set to that number.
        default: RolloutEstimator(RandomSolver(rng))

    init_Q::Any
        Function, object, or number used to set the initial Q(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_Q(o, mdp, s, a)` will be called.
        If this is a number, Q will always be set to that number.
        default: 0.0

    init_N::Any
        Function, object, or number used to set the initial N(s,a) value at a new node.
        If this is a function `f`, `f(mdp, s, a)` will be called to set the value.
        If this is an object `o`, `init_N(o, mdp, s, a)` will be called.
        If this is a number, N will always be set to that number.
        default: 0

    next_action::Any
        Function or object used to choose the next action to be considered for progressive widening.
        The next action is determined based on the MDP, the state, `s`, and the current `DPWStateNode`, `snode`.
        If this is a function `f`, `f(mdp, s, snode)` will be called to set the value.
        If this is an object `o`, `next_action(o, mdp, s, snode)` will be called.
        default: RandomActionGenerator(rng)

    default_action::Any
        Function, action, or Policy used to determine the action if POMCP fails with exception `ex`.
        If this is a Function `f`, `f(pomdp, belief, ex)` will be called.
        If this is a Policy `p`, `action(p, belief)` will be called.
        If it is an object `a`, `default_action(a, pomdp, belief, ex)` will be called, and if this method is not implemented, `a` will be returned directly.
        default: `ExceptionRethrow()`
"""
mutable struct DPWSolver <: AbstractMCTSSolver
    depth::Int
    exploration_constant::Float64
    n_iterations::Int
    max_time::Float64
    k_action::Float64
    alpha_action::Float64
    k_state::Float64
    alpha_state::Float64
    keep_tree::Bool
    enable_action_pw::Bool
    check_repeat_state::Bool
    check_repeat_action::Bool
    tree_in_info::Bool
    rng::AbstractRNG
    estimate_value::Any
    init_Q::Any
    init_N::Any
    next_action::Any
    default_action::Any
end

"""
    DPWSolver()

Use keyword arguments to specify values for the fields
"""
function DPWSolver(;depth::Int=10,
                    exploration_constant::Float64=1.0,
                    n_iterations::Int=100,
                    max_time::Float64=Inf,
                    k_action::Float64=10.0,
                    alpha_action::Float64=0.5,
                    k_state::Float64=10.0,
                    alpha_state::Float64=0.5,
                    keep_tree::Bool=false,
                    enable_action_pw::Bool=true,
                    check_repeat_state::Bool=true,
                    check_repeat_action::Bool=true,
                    tree_in_info::Bool=false,
                    rng::AbstractRNG=Random.GLOBAL_RNG,
                    estimate_value::Any = RolloutEstimator(RandomSolver(rng)),
                    init_Q::Any = 0.0,
                    init_N::Any = 0,
                    next_action::Any = RandomActionGenerator(rng),
                    default_action::Any = ExceptionRethrow()
                   )
    DPWSolver(depth, exploration_constant, n_iterations, max_time, k_action, alpha_action, k_state, alpha_state, keep_tree, enable_action_pw, check_repeat_state, check_repeat_action, tree_in_info, rng, estimate_value, init_Q, init_N, next_action, default_action)
end

#=
mutable struct StateActionStateNode
    N::Int
    R::Float64
    StateActionStateNode() = new(0,0)
end

mutable struct DPWStateActionNode{S}
    V::Dict{S,StateActionStateNode}
    N::Int
    Q::Float64
    DPWStateActionNode(N,Q) = new(Dict{S,StateActionStateNode}(), N, Q)
end

mutable struct DPWStateNode{S,A} <: AbstractStateNode
    A::Dict{A,DPWStateActionNode{S}}
    N::Int
    DPWStateNode{S,A}() where {S,A} = new(Dict{A,DPWStateActionNode{S}}(),0)
end
=#

mutable struct DPWTree{S,A}
    # for each state node
    total_n::Vector{Int}
    children::Vector{Vector{Int}}
    s_labels::Vector{S}
    s_lookup::Dict{S, Int}

    # for each state-action node
    n::Vector{Int}
    q::Vector{Float64}
    transitions::Vector{Vector{Tuple{Int,Float64}}}
    a_labels::Vector{A}
    a_lookup::Dict{Tuple{Int,A}, Int}

    # for tracking transitions
    n_a_children::Vector{Int}
    unique_transitions::Set{Tuple{Int,Int}}

    function DPWTree{S,A}(sz::Int=1000) where {S,A}
        sz = min(sz, 100_000)
        return new(sizehint!(Int[], sz),
                   sizehint!(Vector{Int}[], sz),
                   sizehint!(S[], sz),
                   Dict{S, Int}(),

                   sizehint!(Int[], sz),
                   sizehint!(Float64[], sz),
                   sizehint!(Vector{Tuple{Int,Float64}}[], sz),
                   sizehint!(A[], sz),
                   Dict{Tuple{Int,A}, Int}(),

                   sizehint!(Int[], sz),
                   Set{Tuple{Int,Int}}()
                  )
    end
end

function insert_state_node!(tree::DPWTree{S,A}, s::S, maintain_s_lookup=true) where {S,A}
    push!(tree.total_n, 0)
    push!(tree.children, Int[])
    push!(tree.s_labels, s)
    snode = length(tree.total_n)
    if maintain_s_lookup
        tree.s_lookup[s] = snode
    end
    return snode
end

function insert_action_node!(tree::DPWTree{S,A}, snode::Int, a::A, n0::Int, q0::Float64, maintain_a_lookup=true) where {S,A}
    push!(tree.n, n0)
    push!(tree.q, q0)
    push!(tree.a_labels, a)
    push!(tree.transitions, Vector{Tuple{Int,Float64}}[])
    sanode = length(tree.n)
    push!(tree.children[snode], sanode)
    push!(tree.n_a_children, 0)
    if maintain_a_lookup
        tree.a_lookup[(snode, a)] = sanode
    end
    return sanode
end

Base.isempty(tree::DPWTree) = isempty(tree.n) && isempty(tree.q)

struct DPWStateNode{S,A} <: AbstractStateNode
    tree::DPWTree{S,A}
    index::Int
end

children(n::DPWStateNode) = n.tree.children[n.index]
n_children(n::DPWStateNode) = length(children(n))
isroot(n::DPWStateNode) = n.index == 1

mutable struct DPWPlanner{P<:Union{MDP,POMDP}, S, A, SE, NA, RNG} <: AbstractMCTSPlanner{P}
    solver::DPWSolver
    mdp::P
    tree::Union{Nothing, DPWTree{S,A}}
    solved_estimate::SE
    next_action::NA
    rng::RNG
end

function DPWPlanner(solver::DPWSolver, mdp::P) where P<:Union{POMDP,MDP}
    se = convert_estimator(solver.estimate_value, solver, mdp)
    return DPWPlanner{P,
                      statetype(P),
                      actiontype(P),
                      typeof(se),
                      typeof(solver.next_action),
                      typeof(solver.rng)}(solver,
                                          mdp,
                                          nothing,
                                          se,
                                          solver.next_action,
                                          solver.rng
                     )
end

Random.seed!(p::DPWPlanner, seed) = Random.seed!(p.rng, seed)


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

# -------------
# COMPOSITIONS
# -------------

"""
    ilr(c)

Isometric log-ratio transformation of composition `c`.
"""
function ilr(c::Composition{D}) where {D}
  w = components(c) .+ eps()
  l = log.(w)
  T = eltype(w)
  x = MVector(ntuple(i->zero(T), D-1))
  for i in 1:D-1
    s = zero(T)
    sqrtinv = 1/sqrt(i*(i+1))
    for j in 1:i+1
      if j < i+1
        s += - sqrtinv * l[j]
      elseif j==i+1
        s += i * sqrtinv * l[j]
      end
    end
    x[i] = s
  end
  SVector(x)
end

"""
    ilrinv(x)

Inverse ilr transformation of coordinates `x`.
"""
function ilrinv(x::SVector{D}) where {D}
  z = MVector(ntuple(i->0.0, D+1))
  for i in 1:D+1
    s = 0.0
    for j in 1:D
      sqrtinv = 1/sqrt(j*(j+1))
      if i < j+1
        s += - sqrtinv * x[j]
      elseif i == j+1
        s += j * sqrtinv * x[j]
      end
    end
    z[i] = exp(s)
  end
  Composition(𝒞(z))
end

ilrinv(x::AbstractVector) = ilrinv(SVector{length(x)}(x))

# -------
# TABLES
# -------

"""
    ILR([refvar])

Isometric log-ratio transform following the
[TableTransforms.jl](https://github.com/JuliaML/TableTransforms.jl)
interface.

Optionally, specify the reference variable `refvar` for the ratios.
Default to the last column of the input table.
"""
struct ILR <: LogRatio
  refvar::Union{Symbol,Nothing}
end

ILR() = ILR(nothing)

refvar(transform::ILR, vars) =
  isnothing(transform.refvar) ? last(vars) : transform.refvar

newvars(::ILR, n) = collect(n)[begin:end-1]

oldvars(::ILR, vars, rvar) = [collect(vars); rvar]

applymatrix(::ILR, X) = mapslices(ilr ∘ Composition, X, dims=2)

function revertmatrix(::ILR, Y)
  D = size(Y, 2)
  f = components ∘ ilrinv ∘ SVector{D}
  mapslices(f, Y, dims=2)
end


text:
# START TYPE ExploreSolver
mutable struct ExploreSolver
    inst::Instance
    nb_test::Int          # Nombre total de voisins testés
    nb_move::Int          # nombre de voisins acceptés (améliorants ou non)
    nb_degrad::Int        # nombre de mouvements dégradants
    nb_improve::Int       # Nombre de mouvements améliorants

    cursol::Solution      # Solution courante
    bestsol::Solution     # meilleure Solution rencontrée

    function ExploreSolver(inst::Instance)
        this=new()
        this.inst = inst
        this.nb_test = 0
        this.nb_move = 0
        this.nb_degrad = 0
        this.nb_improve = 0

        this.cursol = Solution(inst)
        ln2("Solution correspondant à l'ordre de l'instance")
        ln2(to_s(this.cursol))

        initial_sort!(this.cursol, presort=:shuffle)
        ln2("Solution initiale envoyée au solver")
        ln2(to_s(this.cursol))

        this.bestsol = Solution(this.cursol)

        return this
    end
end

function solve(sv::ExploreSolver, itermax::Int, neighbour_operator!)
    ln2("BEGIN solve(ExploreSolver)")
    ### Début du glouton  ### FIXME

    # planes_lb = zeros(Int, inst.nb_planes)
    # planes_ub = zeros(Int, inst.nb_planes)
    planes_mean_bound = zeros(sv.inst.nb_planes)

    for (i, plane) in enumerate(sv.inst.planes)
        # planes_lb[i] = plane.lb
        # planes_ub[i] = plane.ub
        planes_mean_bound[i] = (plane.lb + plane.ub)
        # pas besoin de diviser par 2 car le tri sera le même
    end

    # VECT = low+planes_ub

    # ordre des avions selon la quantité (e_i+l_i)/2
    order = sortperm(planes_mean_bound)
    permu!(sv.cursol, 1:sv.inst.nb_planes, order)

    ### Fin du glouton  ###

    lg1("iter <nb_move>=<nb_improve>+<nb_degrade> => <bestcost>")

    current_costs = Vector{Float64}(undef, itermax)
    for iter in 1:itermax
        prevcost = sv.cursol.cost
        neighbour_operator!(sv.cursol)
        sv.nb_move += 1
        degrad = sv.cursol.cost - prevcost
        ln4("degrad=$(degrad)")
        if degrad < 0
            # Ce voisin est meilleur : on l'accepte
            lg3("+")
            sv.nb_improve += 1
            # mise a jour éventuelle de la meilleure solution
            if sv.cursol.cost < sv.bestsol.cost
                # La sauvegarde dans bestsol n'est utile que si on ne fait une descente pure
                copy!(sv.bestsol, sv.cursol)
                if lg1()
                    msg = string("\niter ", sv.nb_move, "=", sv.nb_improve, "+",
                                 sv.nb_degrad)
                    if lg2()
                        # affiche coût + ordre des avions
                        msg *= string(" => ", to_s(sv.bestsol))
                    else
                        # affiche seulement le coût
                        msg *= string(" => ", sv.bestsol.cost)
                    end
                    print(msg)
                end
            end
        else # degrad < 0
            # Ce voisin est plus mauvais : on l'accepte aussi (car exploration) !!
            sv.nb_degrad += 1
            if Args.get(:level) in 3:3
                print("-")
            end
            if lg4()
                msg = string("\n     ", sv.nb_move, ":", sv.nb_improve, "+/",
                             sv.nb_degrad, "- cursol=", to_s(sv.cursol) )
                print(msg)
            end
        end
        current_costs[iter] = sv.cursol.cost
    end
    ln2("\nEND solve(ExploreSolver)")
    return current_costs
end


text:
function make_figure8_path()
  X = Float64[]
  Y = Float64[]
  N = 100
  R = 40

  y = collect(range(0.0, stop=R, length=N))
  x = fill(0.0, length(y))
  append!(X, x)
  append!(Y, y)

  th = range(0.0, stop=(3/2 * pi), length=(3 * N))
  y = similar(th)
  x = similar(th)
  for i in 1:length(th)
    x[i] = R * cos(th[i]) - R
    y[i] = R * sin(th[i]) + R
  end
  append!(X, x)
  append!(Y, y)

  x = collect(range(-R, stop=R, length=(2 * N)))
  y = fill(0.0, length(x))
  append!(X, x)
  append!(Y, y)

  th = range(pi / 2, stop=(-pi), length=(3 * N))
  y = similar(th)
  x = similar(th)
  for i in 1:length(th)
    x[i] = R * cos(th[i]) + R
    y[i] = R * sin(th[i]) - R
  end
  append!(X, x)
  append!(Y, y)

  y = collect(range(-R, stop=0.0, length=N))
  x = fill(0.0, length(y))
  append!(X, x)
  append!(Y, y)

  #remove_duplicates(X, Y)
  (X, Y) = resample(X, Y, 300)

  # compute the first derivative and normalized
  dX = similar(X)
  dY = similar(Y)

  for i = 1:length(X)-1
    dX[i] = X[i + 1] - X[i]
    dY[i] = Y[i + 1] - Y[i]
  end
  dX[end] = X[end] - X[end-1]
  dY[end] = Y[end] - Y[end-1]

  # second derivative
  ddX = similar(dX)
  ddY = similar(dY)
  for i = 1:length(X)-1
    ddX[i] = dX[i + 1] - dX[i]
    ddY[i] = dY[i + 1] - dY[i]
  end
  ddX[end] = X[end] - X[end-1]
  ddY[end] = Y[end] - Y[end-1]

  Sx = similar(X)
  Sy = similar(X)
  Px = similar(X)
  Py = similar(X)
  for i in 1:length(X)
    len = sqrt(dX[i]^2 + dY[i]^2)
    Sx[i] = dX[i] / len
    Sy[i] = dY[i] / len

    Px[i] = -Sy[i]
    Py[i] = Sx[i]
  end

  # path and curvature
  olds = 0.0
  S = similar(X)
  R = similar(X)
  for i = 1:length(X)
    dr = [dX[i]; dY[i]]
    ddr = [ddX[i]; ddY[i]]

    k = cross([dr; 0], [ddr; 0])[3] / norm(dr, 2)^3
    R[i] = (abs(k) < 1e-3) ? Inf : 1 / k
    S[i] = olds
    olds += norm(dr, 2)
  end

  return Path(X, Y, dX, dY, ddX, ddY, Sx, Sy, Px, Py, R, S)
end

function remove_duplicates(x::Array{Float64}, y::Array{Float64})
  @assert length(x) == length(y)
  len = length(x)
  i = 1
  while i < len
    if x[i] == x[i + 1] && y[i] == y[i + 1]
      deleteat!(x, i)
      deleteat!(y, i)

      len = length(x)
    end
    i += 1
  end
end

function resample(x::Array{Float64}, y::Array{Float64}, n::Int)
  @assert length(x) == length(y) > 2

  xp = fill(0.0, n)
  yp = fill(0.0, n)

  s = cumsum(sqrt.(x.^2 + y.^2))
  sp = range(0.0, stop=s[end], length=n)

  xp[1] = x[1]
  yp[1] = y[1]
  for i in 2:(n-1)
    idx = binary_search(s, sp[i])
    idx = idx == length(x) ? idx - 1 : idx
    xp[i] = (x[idx+1] - x[idx])/(s[idx+1] - s[idx]) * (sp[i] - s[idx]) + x[idx]
    yp[i] = (y[idx+1] - y[idx])/(s[idx+1] - s[idx]) * (sp[i] - s[idx]) + y[idx]
  end
  xp[end] = x[end]
  yp[end] = y[end]

  return (xp, yp)
end


text:
# The file test.mseed comes from an older IRIS libmseed, found by anowacki
# It has a more complicated structure than the test.mseed file in more recent
# versions of libmseed, which reads with no issues
printstyled("  mini-SEED file read\n", color=:light_green)

test_sac_file     = string(path, "/SampleFiles/SAC/test_le.sac")
test_mseed_file   = string(path, "/SampleFiles/SEED/test.mseed")
test_mseed_pat    = string(path, "/SampleFiles/SEED/t*.mseed")
mseed_vals_file   = string(path, "/SampleFiles/SEED/test_mseed_vals.txt")

@test_throws ErrorException read_data("mseed", test_sac_file)

S = read_data("mseed", test_mseed_file, v=0)
@test isequal(S.id[1], "NL.HGN.00.BHZ")
@test ≈(S.fs[1], 40.0)
@test ≈(S.gain[1], 1.0)
@test isequal(string(u2d(S.t[1][1,2]*1.0e-6)), "2003-05-29T02:13:22.043")
@test ≈(S.x[1][1:5], [ 2787, 2776, 2774, 2780, 2783 ])

# Test breaks if memory-resident SeisIOBuf structure SEED is not reset
S1 = read_data("mseed", test_mseed_file, v=0)
if Sys.iswindows() == false
  S2 = read_data("mseed", test_mseed_pat, v=0)
  @test S == S1 == S2
end

mseed_vals = readdlm(mseed_vals_file, ',', comments=true, comment_char='#')
seedvals = Dict{String,Any}()
ntest = size(mseed_vals,1)
for i = 1:ntest
  seedvals[mseed_vals[i,1]] = Float32.(mseed_vals[i, 2:end])
end

if safe_isdir(path*"/SampleFiles/Restricted")
  printstyled("    mseed test files (time gaps, blockette types)\n", color=:light_green)
  redirect_stdout(out) do
    seed_support()
    files = ls(path*"/SampleFiles/Restricted/*mseed")
    for f in files
      println(stdout, "attempting to read ", f)
      S = SeisData()
      read_data!(S, "mseed", f, v=3)
      @test isempty(S) == false

      # Test that our encoders return the expected values
      (tmp, fname) = splitdir(f)
      if haskey(seedvals, fname)
        x = get(seedvals, :fname, Float32[])
        nx = lastindex(x)
        y = getindex(getfield(S, :x), 1)[1:nx]
        @test isapprox(x,y)
      end

      if occursin("text-encoded", f)
        @test haskey(S.misc[1], "seed_ascii") == true
        str = split(S.misc[1]["seed_ascii"][1], "\n", keepempty=false)
        @test occursin("Quanterra Packet Baler Model 14 Restart.", str[1])

      elseif occursin("detection.record",f )
        ev_rec = get(S.misc[1], "seed_event", "no record")[1]
        @test ev_rec == "2004,7,28,20,28,6,185,80.0,0.39999998,18.0,dilatation,1,3,2,1,4,0,2,0,Z_SPWWSS"

      elseif occursin("SHW.UW", f)
        @test size(S.t[1],1) >= 158
        @test size(S.t[2],1) >= 8
        @test string(u2d(S.t[1][1,2]*1.0e-6)) == "1980-03-22T20:45:18.349"
        @test isequal(S.id, String[ "UW.SHW..EHZ", "UW.SHW..SHZ" ])
        @test ≈(S.fs, Float64[104.085000, 52.038997])
        @test ≈(S.x[1][1:5], Float64[-68.0, -57.0, -71.0, -61.0, -52.0])
        fnames = ls(path*"/SampleFiles/Restricted/1980*SHZ.D.SAC")

        C = S[2]
        @test w_time(t_win(C.t, C.fs), C.fs) == C.t
        t = t_win(C.t, C.fs)[:,1]
        W = Array{DateTime,1}(undef, 0)
        for i=1:length(t)
          # push!(W, round(u2d(t[i]*1.0e-6), Second))
          push!(W, u2d(t[i]*1.0e-6))
        end

        Y = Array{DateTime,1}(undef,0)
        for f in fnames
          seis = read_data("sac", f)[1]
          push!(Y, u2d(seis.t[1,2]*1.0e-6))
        end
        #[round(u2d(i*1.0e-6), Second) for i in t]
        # println("W = ", string.(W))
        # println("Y = ", string.(Y))
        Δ = [abs(.001*(W[i]-Y[i]).value)*C.fs for i=1:length(Y)]
        @test maximum(Δ) < 1.0
      else
        @test isempty(S) == false
      end
    end
  end
else
  printstyled("  extended SEED tests skipped. (files not found; is this Appveyor?)\n", color=:green)
end


text:
"""
    parity_sign(P)

Given an array or tuple of integers representing
a permutation or a subset of a permutation,
compute the parity sign defined as -1 for a
permutation consisting of an odd number of swaps
and +1 for an even number of swaps. This
implementation uses an O(n^2) algorithm and is
intended for small permutations only.
"""
function parity_sign(P)::Int
  L = length(P)
  s = +1
  for i in 1:L, j in (i + 1):L
    s *= sign(P[j] - P[i])
  end
  return s
end

isfermionic(qv::QNVal) = (modulus(qv) < 0)

isfermionic(qn::QN) = any(isfermionic, qn)

has_fermionic_subspaces(i::Index) = false

function has_fermionic_subspaces(i::QNIndex)
  for b in 1:nblocks(i)
    isfermionic(qn(i, b)) && (return true)
  end
  return false
end

isfermionic(i::Index) = has_fermionic_subspaces(i)

has_fermionic_subspaces(is::IndexSet) = false

function has_fermionic_subspaces(is::Union{QNIndexSet,NTuple{N,QNIndex}}) where {N}
  for i in is, b in 1:nblocks(i)
    isfermionic(qn(i, b)) && (return true)
  end
  return false
end

has_fermionic_subspaces(T) = has_fermionic_subspaces(inds(T))

"""
    fparity(qn::QN)
    fparity(qn::IndexVal)

Compute the fermion parity (0 or 1) of a QN of IndexVal,
defined as the sum mod 2 of each of its fermionic
QNVals (QNVals with negative modulus).
"""
function fparity(qn::QN)
  p = 0
  for qv in qn
    if isfermionic(qv)
      p += val(qv)
    end
  end
  return mod(p, 2)
end

fparity(iv::Pair{<:Index}) = fparity(qn(iv))

Base.isodd(q::QN) = isodd(fparity(q))
Base.isodd(iv::Pair{<:Index}) = isodd(fparity(iv))

"""
    compute_permfactor(p,iv_or_qn::Vararg{T,N})

Given a permutation p and a set "s" of QNIndexVals or QNs,
if the subset of index vals which are fermion-parity
odd undergo an odd permutation (odd number of swaps)
according to p, then return -1. Otherwise return +1.
"""
function compute_permfactor(p, iv_or_qn...; range=1:length(iv_or_qn))::Int
  using_auto_fermion() || return 1
  N = length(iv_or_qn)
  oddp = @MVector zeros(Int, N)
  n = 0
  @inbounds for j in range
    if fparity(iv_or_qn[p[j]]) == 1
      n += 1
      oddp[n] = p[j]
    end
  end
  return parity_sign(oddp[1:n])
end

# Default implementation for non-QN IndexVals
NDTensors.permfactor(p, ivs...) = 1

function NDTensors.permfactor(p, ivs::Vararg{Pair{QNIndex},N}; kwargs...) where {N}
  using_auto_fermion() || return 1
  return compute_permfactor(p, ivs...; kwargs...)
end

function NDTensors.permfactor(
  perm, block::NDTensors.Block{N}, inds::Union{QNIndexSet,NTuple{N,QNIndex}}; kwargs...
) where {N}
  using_auto_fermion() || return 1
  qns = ntuple(n -> qn(inds[n], block[n]), N)
  return compute_permfactor(perm, qns...; kwargs...)
end

# Version of getperm which is type stable
# and works for Tuple or Vector inputs
function vec_getperm(s1, s2)
  N = length(s1)
  p = Vector{Int}(undef, N)
  for i in 1:N
    @inbounds p[i] = NDTensors._findfirst(==(@inbounds s1[i]), s2)
  end
  return p
end

@inline function NDTensors.compute_alpha(
  ElR,
  labelsR,
  blockR,
  input_indsR,
  labelsT1,
  blockT1,
  indsT1::NTuple{N1,QNIndex},
  labelsT2,
  blockT2,
  indsT2::NTuple{N2,QNIndex},
) where {N1,N2}
  if !using_auto_fermion()
    !has_fermionic_subspaces(indsT1) || !has_fermionic_subspaces(indsT2)
    return one(ElR)
  end

  # the "indsR" argument to compute_alpha from NDTensors
  # may be a tuple of QNIndex, so convert to an IndexSet
  indsR = IndexSet(input_indsR)

  nlabelsT1 = NDTensors.sort(labelsT1; rev=true)
  nlabelsT2 = NDTensors.sort(labelsT2)

  # Make orig_labelsR from the order of
  # indices that would result by just
  # taking the uncontracted indices of
  # T1 and T2 in their input order:
  NR = length(labelsR)
  orig_labelsR = zeros(Int, NR)
  u = 1
  for ls in (nlabelsT1, nlabelsT2), l in ls
    if l > 0
      orig_labelsR[u] = l
      u += 1
    end
  end

  permT1 = NDTensors.getperm(nlabelsT1, labelsT1)
  permT2 = NDTensors.getperm(nlabelsT2, labelsT2)
  permR = vec_getperm(labelsR, orig_labelsR)

  alpha1 = NDTensors.permfactor(permT1, blockT1, indsT1)
  alpha2 = NDTensors.permfactor(permT2, blockT2, indsT2)
  alphaR = NDTensors.permfactor(permR, blockR, indsR)

  alpha_arrows = one(ElR)
  for n in 1:length(indsT1)
    l = labelsT1[n]
    i = indsT1[n]
    qi = qn(i, blockT1[n])
    if l < 0 && dir(i) == Out && fparity(qi) == 1
      alpha_arrows *= -1
    end
  end

  α = alpha1 * alpha2 * alphaR * alpha_arrows

  return α
end

# Flip signs of selected blocks of T prior to
# it being multiplied by a combiner ITensor
# labelsR gives the ordering of indices after the product
function NDTensors.before_combiner_signs(
  T,
  labelsT_,
  indsT::NTuple{NT,QNIndex},
  C,
  labelsC_,
  indsC::NTuple{NC,QNIndex},
  labelsR,
  indsR::NTuple{NR,QNIndex},
) where {NC,NT,NR}
  if !using_auto_fermion() || !has_fermionic_subspaces(T)
    return T
  end

  T = copy(T)

  labelsC = [l for l in labelsC_]
  labelsT = [l for l in labelsT_]

  # number of uncombined indices
  Nuc = NC - 1

  ci = NDTensors.cinds(storage(C))[1]
  combining = (labelsC[ci] > 0)

  isconj = NDTensors.isconj(storage(C))

  if combining
    #println("Combining <<<<<<<<<<<<<<<<<<<<<<<<<<<")

    nlabelsT = Int[]

    if !isconj
      # Permute uncombined indices to front
      # in same order as indices passed to the
      # combiner constructor
      append!(nlabelsT, labelsC[2:end])
    else # isconj
      # If combiner is conjugated, put uncombined
      # indices in *opposite* order as on combiner
      append!(nlabelsT, reverse(labelsC[2:end]))
    end
    @assert all(l -> l < 0, nlabelsT)

    for l in labelsT
      if l > 0 #uncontracted
        append!(nlabelsT, l)
      end
    end
    @assert length(nlabelsT) == NT

    # Compute permutation that moves uncombined indices to front
    permT = vec_getperm(nlabelsT, labelsT)

    for blockT in keys(blockoffsets(T))
      # Compute sign from permuting uncombined indices to front:
      alphaT = NDTensors.permfactor(permT, blockT, indsT)

      neg_dir = !isconj ? In : Out
      alpha_arrows = 1
      alpha_mixed_arrow = 1
      C_dir = dir(indsC[1])
      for n in 1:length(indsT)
        i = indsT[n]
        qi = qn(i, blockT[n])
        if labelsT[n] < 0 && fparity(qi) == 1
          alpha_mixed_arrow *= (dir(i) != C_dir) ? -1 : +1
          alpha_arrows *= (dir(i) == neg_dir) ? -1 : +1
        end
      end

      fac = alphaT * alpha_arrows

      if isconj
        fac *= alpha_mixed_arrow
      end

      if fac != 1
        Tb = blockview(T, blockT)
        scale!(Tb, fac)
      end
    end # for blockT

  elseif !combining
    #
    # Uncombining ---------------------------
    #
    #println("Uncombining >>>>>>>>>>>>>>>>>>>>>>>>>>>")

    nc = findfirst(l -> l < 0, labelsT)
    nlabelsT = [labelsT[nc]]
    ic = indsT[nc]

    for l in labelsT
      (l > 0) && append!(nlabelsT, l)
    end

    # Compute sign for permuting combined index to front
    # (sign alphaT to be computed for each block below):
    permT = vec_getperm(nlabelsT, labelsT)

    #
    # Note: other permutation of labelsT which
    # relates to two treatments of isconj==true/false
    # in combining case above is handled as a
    # post-processing step in NDTensors.after_combiner_signs
    # implemented below
    #

    for blockT in keys(blockoffsets(T))
      alphaT = NDTensors.permfactor(permT, blockT, indsT)

      neg_dir = !isconj ? Out : In
      qic = qn(ic, blockT[nc])
      alpha_arrows = (fparity(qic) == 1 && dir(ic) == neg_dir) ? -1 : +1

      fac = alphaT * alpha_arrows

      if fac != 1
        Tb = blockview(T, blockT)
        scale!(Tb, fac)
      end
    end
  end

  return T
end #NDTensors.before_combiner_signs

function NDTensors.after_combiner_signs(
  R, labelsR, indsR::NTuple{NR,QNIndex}, C, labelsC, indsC::NTuple{NC,QNIndex}
) where {NC,NT,NR}
  ci = NDTensors.cinds(store(C))[1]
  combining = (labelsC[ci] > 0)
  combining && error("NDTensors.after_combiner_signs only for uncombining")

  if !using_auto_fermion() || !has_fermionic_subspaces(R)
    return R
  end

  R = copy(R)

  # number of uncombined indices
  Nuc = NC - 1

  isconj = NDTensors.isconj(store(C))

  if !combining
    if !isconj
      #println("!!! Doing uncombining post-processing step")
      rperm = ntuple(i -> (Nuc - i + 1), Nuc) # reverse permutation
      NDTensors.scale_blocks!(
        R, block -> NDTensors.permfactor(rperm, block, indsR; range=1:Nuc)
      )
    else
      #println("!!! Doing conjugate uncombining post-processing step")
      C_dir = dir(inds(C)[1])

      function mixed_arrow_sign(block)
        alpha_mixed_arrow = 1
        for n in 1:Nuc
          i = indsR[n]
          qi = qn(i, block[n])
          if dir(i) == C_dir && fparity(qi) == 1
            alpha_mixed_arrow *= -1
          end
        end
        return alpha_mixed_arrow
      end

      NDTensors.scale_blocks!(R, block -> mixed_arrow_sign(block))
    end
  end

  return R
end #NDTensors.after_combiner_signs


text:
function tcm_randomCompleteTree( N, depth )

	tree_left = Int64[]
	tree_right = Int64[]
	tree_product = Int64[]
	tree_isLeaf = Bool[]
	#tree_choice_prediction = Int64[]

	tree_unused_products = Array{ Array{Int64,1} }( 2^(depth+1) - 1)
	tree_include_set = Array{ Array{Int64,1} }( 2^(depth+1) - 1)

	tree_unused_products[1] = collect(1:(N-1))
	tree_include_set[1] = [N]

	counter = 1;
	num_nodes = 1;

	for d in 0:(depth-1)
		num_d_nodes = 2^d;

		#next_children_are_leaves = d == n_levels ? true : false;

		for i in 1:num_d_nodes
			push!(tree_left, counter + 1)
			push!(tree_right, counter + 2)

			temp = rand(tree_unused_products[num_nodes])

			push!(tree_product, temp)
			push!(tree_isLeaf, false)
			tree_unused_products[counter+1] = setdiff(tree_unused_products[num_nodes], temp)
			tree_unused_products[counter+2] = setdiff(tree_unused_products[num_nodes], temp)
			tree_include_set[counter+1] = [tree_include_set[num_nodes]; temp]
			tree_include_set[counter+2] = copy(tree_include_set[num_nodes])

			counter += 2
			num_nodes += 1
		end
	end

	num_d_nodes = 2^depth;
	for i in 1:num_d_nodes
		push!(tree_left, 0)
		push!(tree_right, 0)
		push!(tree_product, rand(tree_include_set[num_nodes]))
		push!(tree_isLeaf, true)
		counter += 2
		num_nodes += 1
	end

	return tree_left, tree_right, tree_product, tree_isLeaf
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
import Automa
import Automa.RegExp: @re_str
using BenchmarkTools
using Random: seed!

seed!(1234)
data = String(vcat([push!(rand(b"ACGTacgt", 59), UInt8('\n')) for _ in 1:1000]...))

VISUALIZE = false
function writesvg(name, machine)
    dot = joinpath(@__DIR__, "$(name).dot")
    svg = joinpath(@__DIR__, "$(name).svg")
    info("writing $(dot)")
    write(dot, Automa.machine2dot(machine))
    info("writing $(svg)")
    run(`dot -Tsvg -o $(svg) $(dot)`)
end

# Case 1
# ------

println("Case 1 ", raw"([A-z]*\r?\n)*")
match(data) = occursin(r"^(:?[A-z]*\r?\n)*$", data)
@assert match(data)
println("PCRE:                 ", @benchmark match(data))

machine = Automa.compile(re"([A-z]*\r?\n)*")
VISUALIZE && writesvg("case1", machine)
context = Automa.CodeGenContext(generator=:goto, checkbounds=false)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl:            ", @benchmark match(data))

context = Automa.CodeGenContext(generator=:goto, checkbounds=false, loopunroll=10)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl (unrolled): ", @benchmark match(data))

# Case 2
# ------

println()
println("Case 2 ", raw"([A-Za-z]*\r?\n)*")
match(data) = occursin(r"^(:?[A-Za-z]*\r?\n)*$", data)
@assert match(data)
println("PCRE:                 ", @benchmark match(data))

machine = Automa.compile(re"([A-Za-z]*\r?\n)*")
VISUALIZE && writesvg("case2", machine)
context = Automa.CodeGenContext(generator=:goto, checkbounds=false)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl:            ", @benchmark match(data))

context = Automa.CodeGenContext(generator=:goto, checkbounds=false, loopunroll=10)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl (unrolled): ", @benchmark match(data))

# Case 3
# ------

println()
println("Case 3 ", raw"([ACGTacgt]*\r?\n)*")
match(data) = occursin(r"^(:?[ACGTacgt]*\r?\n)*$", data)
@assert match(data)
println("PCRE:                 ", @benchmark match(data))

machine = Automa.compile(re"([ACGTacgt]*\r?\n)*")
VISUALIZE && writesvg("case3", machine)
context = Automa.CodeGenContext(generator=:goto, checkbounds=false)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl:            ", @benchmark match(data))

context = Automa.CodeGenContext(generator=:goto, checkbounds=false, loopunroll=10)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl (unrolled): ", @benchmark match(data))

# Case 4
# ------

println()
println("Case 4 ", raw"([A-Za-z\*-]*\r?\n)*")
match(data) = occursin(r"^(:?[A-Za-z\*-]*\r?\n)*$", data)
@assert match(data)
println("PCRE:                 ", @benchmark match(data))

machine = Automa.compile(re"([A-Za-z\*-]*\r?\n)*")
VISUALIZE && writesvg("case4", machine)
context = Automa.CodeGenContext(generator=:goto, checkbounds=false)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl:            ", @benchmark match(data))

context = Automa.CodeGenContext(generator=:goto, checkbounds=false, loopunroll=10)
@eval function match(data)
    $(Automa.generate_init_code(context, machine))
    p_end = p_eof = lastindex(data)
    $(Automa.generate_exec_code(context, machine))
    return cs == 0
end
@assert match(data)
println("Automa.jl (unrolled): ", @benchmark match(data))


text:
export @primitive

const PRIMITIVES = Any[]

"""
    generate_forward_stub(name::Symbol, op)

Generate forward stub which forward primitive circuit to instruction interfaces.
"""
function generate_forward_stub(name::Symbol, op)
    quoted_name = QuoteNode(name)
    stub = gensym(name)

    return quote
        function $stub(::$(Circuit){$quoted_name}, r::$(AbstractRegister), locs::$(Locations))
            $(YaoAPI).instruct!(r, $op, Tuple(locs))
            return
        end

        function $stub(::$(Circuit){$quoted_name}, r::$(TraceTape), locs::$(Locations))
            $(GlobalRef(Compiler, :trace!))(
                r,
                Expr(:quantum, :gate, GlobalRef(YaoLang, $quoted_name), locs),
            )
            return
        end

        function $stub(
            ::$(Circuit){$quoted_name},
            r::$(AbstractRegister),
            locs::$(Locations),
            ctrl_locs::$(CtrlLocations),
        )
            raw_ctrl_locs = Tuple(ctrl_locs.storage)
            ctrl_cfg = map(Int, (ctrl_locs.configs...,))
            $(YaoAPI).instruct!(r, $op, Tuple(locs), raw_ctrl_locs, ctrl_cfg)
            return
        end

        function $stub(
            ::$(Circuit){$quoted_name},
            r::$(TraceTape),
            locs::$(Locations),
            ctrl_locs::$(CtrlLocations),
        )
            $(GlobalRef(Compiler, :trace!))(
                r,
                Expr(:quantum, :ctrl, GlobalRef(YaoLang, $quoted_name), locs, ctrl_locs),
            )
            return
        end

        (::$PrimitiveCircuit{$quoted_name})() = $Circuit{$quoted_name}($stub)
        Core.@__doc__ const $name = $PrimitiveCircuit{$quoted_name}()
        push!($(GlobalRef(Compiler, :PRIMITIVES)), $name)
        $name
    end
end

function primitive_m(x::Symbol)
    generate_forward_stub(x, :(Val($(QuoteNode(x)))))
end

function primitive_m(ex::Expr)
    def = splitdef(ex; throw = false)
    def === nothing && return assign_statement(ex)

    haskey(def, :name) || throw(Meta.ParseError("Invalid Syntax: expect a function name"))
    name = def[:name]
    quoted_name = QuoteNode(name)
    stub = gensym(name)
    if haskey(def, :args)
        args = map(rm_annotations, def[:args])
    else
        args = ()
    end

    mat_stub_def = deepcopy(def)
    mat_stub = gensym(:mat)
    mat_stub_def[:name] = mat_stub

    primitive_def = deepcopy(def)
    primitive_def[:name] = :(::$(PrimitiveCircuit{name}))
    primitive_def[:body] = quote
        m = $(Expr(:call, mat_stub, args...))
        return Circuit{$quoted_name}($stub, (m, $(args...)))
    end

    circ = gensym(:circ)
    register = gensym(:register)
    locs = gensym(:locs)
    ctrl_locs = gensym(:ctrl_locs)
    matrix = gensym(:m)

    stub_def = Dict{Symbol,Any}()
    stub_def[:name] = stub
    stub_def[:args] =
        Any[:($circ::Circuit{$quoted_name}), :($register::$AbstractRegister), :($locs::$Locations)]
    stub_def[:body] = quote
        $matrix = $circ.free[1]
        YaoAPI.instruct!($register, $matrix, Tuple($locs))
        return
    end

    ctrl_stub_def = Dict{Symbol,Any}()
    ctrl_stub_def[:name] = stub
    ctrl_stub_def[:args] = Any[
        :($circ::Circuit{$quoted_name}),
        :($register::$AbstractRegister),
        :($locs::$Locations),
        :($ctrl_locs::$CtrlLocations),
    ]
    ctrl_stub_def[:body] = quote
        $matrix = $circ.free[1]
        # issue #10
        raw_ctrl_locs = ($(ctrl_locs).storage...,)
        ctrl_cfg = map(Int, ($(ctrl_locs).configs...,))
        YaoAPI.instruct!($register, $matrix, Tuple($locs), raw_ctrl_locs, ctrl_cfg)
        return
    end

    trace_stub_def = Dict{Symbol,Any}()
    trace_stub_def[:name] = stub
    trace_stub_def[:args] =
        Any[:($circ::Circuit{$quoted_name}), :($register::$TraceTape), :($locs::$Locations)]
    trace_stub_def[:body] = Expr(
        :call,
        GlobalRef(Compiler, :trace!),
        register,
        :(Expr(
            :quantum,
            :gate,
            Expr(:call, GlobalRef(YaoLang, $quoted_name), $circ.free[2:end]...),
            $locs,
        )),
    )

    trace_ctrl_stub_def = Dict{Symbol,Any}()
    trace_ctrl_stub_def[:name] = stub
    trace_ctrl_stub_def[:args] = Any[
        :($circ::Circuit{$quoted_name}),
        :($register::$TraceTape),
        :($locs::$Locations),
        :($ctrl_locs::$CtrlLocations),
    ]
    trace_ctrl_stub_def[:body] = Expr(
        :call,
        GlobalRef(Compiler, :trace!),
        register,
        :(Expr(
            :quantum,
            :ctrl,
            Expr(:call, GlobalRef(YaoLang, $quoted_name), $circ.free[2:end]...),
            $locs,
            $ctrl_locs,
        )),
    )

    quote
        $(combinedef(mat_stub_def))
        $(combinedef(stub_def))
        $(combinedef(ctrl_stub_def))
        $(combinedef(primitive_def))
        $(combinedef(trace_stub_def))
        $(combinedef(trace_ctrl_stub_def))
        Core.@__doc__ const $name = $(PrimitiveCircuit{name})()
        push!($(GlobalRef(Compiler, :PRIMITIVES)), $name)
        $name
    end
end

function assign_statement(ex::Expr)
    ex.head === :(=) ||
        throw(Meta.ParseError("Invalid Syntax, expect <primitive gate name> = <matrix expr>, got $ex"))
    ex.args[1] isa Symbol || throw(Meta.ParseError("Invalid Syntax, expect Symbol got $(ex.args[1])"))
    name = ex.args[1]
    matrix_const = gensym(:matrix_const)

    return quote
        const $matrix_const = $(esc(ex.args[2]))
        $(generate_forward_stub(name, matrix_const))
    end
end

"""
    @primitive ex

Define a primitive quantum instruction. `ex` can be a Symbol, if the corresponding instruction
interface of `YaoAPI.instruct!` is implemented. Or `ex` can be an assignment statement for constant
instructions. Or `ex` can be a function that returns corresponding matrix given a set of classical
parameters.

# Example

Since the instructions interface `YaoAPI.instruct!` of Pauli operators are defined, we can use

```julia
@primitive X
```

to declare a Pauli X primitive instruction.

Or we can also define a Hadamard primitive instruction via its matrix form

```julia
@primitive H = [1 1;1 -1]/sqrt(2)
```

For parameterized gates, such as phase shift gate, we can define it as

```julia
@primitive shift(θ::Real) = Diagonal([1.0, exp(im * θ)])
```
"""
macro primitive(ex)
    return esc(primitive_m(ex))
end