text:
# This julia file builds on the Initial.julia
# It is based on the "How to efficiently select an arbitrary Clifford group element"
# Robert Koenig and John Smolin
# arXivv:1406.2170v1

# Copyright Robin Harper 2015-2020

# The premise here is that in order to generate a random Clifford, we can be sure of
# proper (Haar) randomness if either we generate all the cliffords for a certain quibit size
# and select one randomly OR if we can have a 1-1 mapping between the integers and the Cliffords
# and then we can just randomly select an integer. The latter is the one adopted in Koenig's paper.
# getNumberOfCliffords(n) returns the number of cliffords with n qubits.

# This ties into the Aaronsen/Gottesman paper becauase the generated cliffords are generated
# in such a way it is possible to create the tableau that would result IF the start myVector
# |000000> had been acted on by the chosen clifford unitary.
# We can then (decompose) the tableau using the algorithms specified in the Aaronson/Gottesman paper
# (see Initial.jl) to work out what combination of one-qubit (phase and hadmard) and two quibit
# (cnot) gates would create the unitary in the first place.

# The initial part of this file is a port of the python code in the Koenig/Smolin paper
# followed by an implementation of the decomposition algorithm in the Aaronson/Gottesman paper.
# There is also code to:
# -  Draw a decomposed circuit
# -  Display the kets of a |00..0> vector transformed by the circuit
# -  Form a 2^n 2^n complex matrix representing the "raw circuit" for use outside this formalism
# -  There is brute force method for finding the minimum set of gates to construct a given unitary
#    (even with only 3qubits this takes a tediously long time to run)
# -  There is the start of rationalise method, just now it eliminates chains of 4 phase gates or 2 hadamard gates.
#     Maybe the way to go here is to eliminate recognised patterns. [REF] uses a reverse tree building approach
#     which is interesting.
#
# Use is made of the ability of Julia to manipulate its own commands (a sort of hybrid functional/imperative)
# approach. With decomposition we are using arrays stored with the Tableau
# The first "commands" is a text represenetation of the commands needed to rebuild a decomposed end state.
# the second "executeCommands" contains the actual Julia instructions to rebuild the state.
# They are only populated when the appropriate decompositions have been made.
# This is better documented in the appropriate methods.

# The generation of the Clifford requires an element from the symplectic group
#include("Initial.jl")
#using Main.CHP

# Scratch variable in the global spce of Symplectic to be used with eval.
_t = Tableau(3)

# This just shows how big the groups get, Julia overflow is going to be a problem

function getNumberOfCliffords(n::BigInt)::BigInt
		return 2^(n^2+2*n)*prod([4^x-1 for x =1:n])
end

"""
    getNumberOfCliffords(n::Integer)::BigInt

returns the number of Cliffords for n qubits.
"""
function getNumberOfCliffords(n::Integer)::BigInt
		return getNumberOfCliffords(BigInt(n))
end

"""
    getNumberOfSymplecticCliffords(n::Integer)::BigInt

returns the number of Cliffords for n qubits, modulo the signs on each row
of the tableau. Note

getNumberOfCliffords(n) == getNumberOfSymplecticCliffords(n)*getNumerOfBitStringsCliffords(n)
"""
function getNumberOfSymplecticCliffords(n::BigInt)::BigInt
	return 2^(n^2)*prod([4^x-1 for x =1:n])
end

function getNumberOfSymplecticCliffords(n::Integer)::BigInt
	return getNumberOfSymplecticCliffords(BigInt(n))
end

"""
    getNumberOfBitStringsCliffords(n)

returns the number of bitstrings representing a +/- on each row for the Tableau.
Note:

getNumberOfCliffords(n) == getNumberOfSymplecticCliffords(n)*getNumerOfBitStringsCliffords(n)
"""
function getNumberOfBitStringsCliffords(n)
	return 2^(2*n)
end

"""
    cliffordToTableau(qubits,clifford,signs)

Converts an integer into a Clifford. You need to specify the number of qubits in the
tableau, the clifford number and the bit signs corresponding to the rows of the Tableau.
The commands used to create the Clifford can be found in Tableau.commands and Tableau.executeCommands.?
It is fine to specify numbers that are too large (they wrap), but the helper functions
```Julia
getNumberOfSymplecticCliffords(n::BigInt)::BigInt
getNumberOfBitStringsCliffords(n)
```
Will show how many cliffords there are for a certain number of qubits (Warning this is a very large number for n > 2)
"""
function cliffordToTableau(qubits,clifford,signs)
    t = Tableau(qubits)
	t.state = stabiliseSymp(symplectic(clifford,qubits),signs)
	decomposeState(t)
	return t
end

"""
    tableauToClifford(t::Tableau)

Returns the 'symplectic' number that corresponds to the clifford embedded in the Tableau.
Note that this is modulo the bit signs. I.e. the tableaus corresponding to e.g.
cliffordToTableau(1,3,3) and cliffordToTableau(1,3,4) will both return 3.
"""
function tableauToClifford(t::Tableau)
	symplecticinverse(decomposeStateToSymplectic(t.state))
end

"""
    decomposeState(tableau::Tableau,rationalise=true)

Decomposes the state stabilised by the Tableau using the Aaronsen/Gottesman method

That is it will decompose the "arbitrary" state into a series of CHP (cnot, hadamard and phase) gates.
This is unlikely to be the most concise decomposition. Unless rationalise is set to
false there will be some naive trimming of gates, basically self eliminating gates (e.g. 4 phase gates)
are removed. The commands (and the execCommands) in the tableau are wiped and recreated. Following the
decomposition the state is re-created using the new commands.
"""
function decomposeState(tableau::Tableau,rationalise=true)
	global _t
	tableau.commands = []
	tableau.executeCommands = []
	#addCommand(tableau,"output(state)",Expr(:call,:show,:_t))
    ss1=tableau.state
	while getState(ss1) < 11
		nextStep(tableau)
	end
	j=tableau.qubits
	addCommand(tableau,"initialise($j)",Expr(:call,:initialise,:_t))
	reverseCommands(tableau)
	removeRedundancy(tableau)
	_t = tableau
	tracking = _t.trackCommands
	_t.trackCommands = false
	for i in tableau.executeCommands
		eval(i)
	end
	_t.trackCommands = tracking

end

function drawCircuit(t::Tableau)
	commands = t.commands
    qasmCommands=[]
    labels=[]
	currentDir = pwd()
	for i = 1:size(commands,1)
		m = match(r"initialise\((.*)\)",commands[i])
		if (m!=nothing)
			for idx=1:t.qubits
    			#push!(qasmCommands,"qubit q$i")
                push!(labels,"q$idx")
			end
    	else
    		m=match(r"hadamard\((.*)\)",commands[i])
    		if (m!=nothing)
    			push!(qasmCommands,("H","q$(m.captures[1])"))
    		else
    			m=match(r"phase\((.*)\)",commands[i])
    			if (m!=nothing)
                    push!(qasmCommands,("P","q$(m.captures[1])"))
    			else
    				m=match(r"cnot\((.*),(.*)\)",commands[i])
    				if (m!=nothing)
    					push!(qasmCommands,("CNOT","q$(m.captures[2])","q$(m.captures[1])"))
    				end
    			end
    		end
    	end
    end
    return (qasmCommands,labels)
end

function qiskitCircuit(t::Tableau,name = "Circuit")
	commands = t.commands
    qasmCommands=["#Pass in the circuit and qubit register","def create$name(circs,qreg):"]

    labels=[]
	currentDir = pwd()
	for i = 1:size(commands,1)
		m = match(r"initialise\((.*)\)",commands[i])
		if (m!=nothing)
			for idx=1:Meta.parse(m.captures[1])
    			#push!(qasmCommands,"qubit q$i")
                push!(labels,"q$idx")
			end
    	else
    		m=match(r"hadamard\((.*)\)",commands[i])
    		if (m!=nothing)
    			push!(qasmCommands,"    circs.h(qreg[$(m.captures[1])])")
    		else
    			m=match(r"phase\((.*)\)",commands[i])
    			if (m!=nothing)
                    push!(qasmCommands,"    circs.s(qreg[$(m.captures[1])])")
    			else
    				m=match(r"cnot\((.*),(.*)\)",commands[i])
    				if (m!=nothing)
    					push!(qasmCommands,"    circs.cx(qreg[$(m.captures[1])],qreg[$(m.captures[2])])")
    				end
    			end
    		end
    	end
    end
	return join(qasmCommands,"\n")
end

###################### Functions used to create Clifford tabelau.

#takes square arrays and places them in a larger array, as follows
#   m1 0
#   0 m2
function directsum(m1,m2)
  n1=size(m1)[1]
  n2=size(m2)[1]
  out = zeros(n1+n2,n1+n2)
  out[1:n1,1:n1]=m1
  out[(n1+1):end,(n1+1):end]=m2
  return out
end

# Helper function used in decomposing the clifford
# builds up an array of the commands used to call the hadamard/phase/cnot.
# We are using two globals here for the sake of simplicity.
# execute commands is the Julia code required to call the relevant functions.

function addCommand(t::Tableau,comToAdd,comm)
	t.commands = append!(t.commands,[comToAdd])
	t.executeCommands=append!(t.executeCommands,[comm]);
end

# Returns the symplectic inner product of two vectors.

function inner(v,w)
	t=0
	for i in 1:(size(v)[1]>>1)
		 t+=v[2*i-1]*w[2*i]
		 t+=w[2*i-1]*v[2*i]
    end
    return t%2
 end

 function transvection(k,v)
 	if (k==0) return mod.(v,2)
 	end
 	return mod.((v+inner(k,v)*k),2)
 end

 function int2bits(i,n)
 	out = zeros(n)
 	for j in 1:n
 		out[j]=i&1
 		i>>=1
 	end
 	return out
 end

 function findtransvection(x,y)
 	out = zeros(UInt8,2,size(x)[1])
 	if (x==y)
 		return out
 	end
 	if inner(x,y) ==1
 		out[1,:]=mod.((x+y),2)
 		return out
 	end
 	z=zeros(size(x)[1])
 	for i in 1:(size(x)[1]>>1)
 		ii=2*i
 		if ((x[ii-1] + x[ii]) != 0) && ((y[ii-1]+y[ii])!=0)
 			# found the pair
 			z[ii-1]=(x[ii-1]+y[ii-1])%2
 			z[ii]=(x[ii]+y[ii])%2
 			if (z[ii-1]+z[ii])==0 #the same
 				z[ii]=1
 				if (x[ii-1] != x[ii])
 					z[ii-1]=1
 				end
 			end
 			out[1,:]=mod.((x+z),2)
 			out[2,:]=mod.((y+z),2)
 			return out
 		end
 	end
 	# ok we didn't find a pair
 	# need places where x is 00 and y does not
 	# and vice versa.
 	for i in 1:(size(x)[1]>>1)
 		ii=2*i
 		if ((x[ii-1]+x[ii]) !=0) && ((y[ii-1]+y[ii]) == 0)
 			if x[ii-1]==x[ii]
 				z[ii]=1
 			else
 				z[ii]=x[ii-1]
 				z[ii-1]=x[ii]
 			end
 			break
 		end
 	end
 	for i in 1:(size(x)[1]>>1)
 		ii=2*i
 		if ((x[ii-1]+x[ii]) ==0) && ((y[ii-1]+y[ii]) !=0 )
 			if y[ii-1]==y[ii]
 				z[ii]=1
 			else
 				z[ii]=y[ii-1]
 				z[ii-1]=y[ii]
 			end
 			break
 		end
 	end
 	# (x+z)%2 doesn't work on arrays, so mod function instead.
 	out[1,:]=mod.((x+z),2)
 	out[2,:]=mod.((y+z),2)
 	return out
 end

# The steps are as set out on page 5 of the paper.
# We use bigInts because we can (and if n was greater than about 64 we would need to).
function symplectic(i,n)
	nn=2*n
	s=BigInt(2)^nn-1
	k=(i%s)+1
	i=div(i,s)
	# step 2
	# There is probably a better way to do this, but it works.
	f1=[x == '0' ? UInt8(0) : UInt8(1) for x  in reverse(string(k,base=2,pad=nn))]
	#step 3
	e1=zeros(UInt8,nn)
	#define the first basis
	e1[1]=1
	T=findtransvection(e1,f1)
	# step 4
	divisor = BigInt(2)^(nn-1)
	value = i%divisor
	bits = [x == '0' ? UInt8(0) : UInt8(1) for x  in reverse(string(value,base=2,pad=nn-1))]

	#step 5
	# Note that in Julia we need to expressly copy
	eprime=copy(e1)
	for j in  3:nn
		eprime[j]=bits[j-1]
	end
	h0=transvection(T[1:1,:]',eprime)
	h0=transvection(T[2:2,:]',h0)

	#step 6
	if bits[1]==1
		f1*=0
	end

	#step 7
	id2=[1 0;0 1]

	if n!=1
		g=directsum(id2,symplectic(i>>(nn-1),n-1))
	else
		g=id2
	end
	for j in 1:nn
		g[j,:]=transvection(T[1:1,:]',g[j:j,:]')
		g[j,:]=transvection(T[2:2,:]',g[j:j,:]')
		g[j,:]=transvection(h0,g[j:j,:]')
		g[j,:]=transvection(f1,g[j:j,:]')
	end
	return g
end

function numberofcosets(n)
	x = 2^(2*n-1)*(2^(2*n)-1)
	return x
end

function bits2int(b,nn)
	output = BigInt(0)
	tmp = BigInt(1)
	for j in 1:nn
		if b[j]==1
			output = output + tmp
		end
		tmp = tmp*2
	end
	return output
end

function symplecticinverse(gn)
	n = round(Int,size(gn)[1]/2)
	nn=2*n

	# step 1

	v = gn[1,:]
	w = gn[2,:]

	# step 2

	e1 = zeros(nn)
	e1[1] = 1
	T = findtransvection(v,e1)

	# step 3
	tw = copy(w)
	tw = transvection(T[1,:],tw)
	tw = transvection(T[2,:],tw)
	b = tw[1]
	h0 = zeros(nn)
	h0[1]=1
	h0[2]=0
	for j in 3:nn
		h0[j]=tw[j]
	end

	# step 4

	bb = zeros(nn-1)
	bb[1]=b
	for j in 3:nn
		bb[j-1] = tw[j]
	end

	zv = bits2int(v,nn)-1
	zw = bits2int(bb,nn-1)

	cvw = zw*(2^(2*n)-1)+zv

	if n == 1
		return cvw
	end

	gprime = copy(gn)
	if b == 0
		for j in 1:nn
			gprime[j,:] = transvection(T[2,:],transvection(T[1,:],gn[j,:]))
			gprime[j,:] = transvection(h0,gprime[j,:])
			gprime[j,:] = transvection(e1,gprime[j,:])
		end
	else
		for j in 1:nn
			gprime[j,:]=transvection(T[2,:],transvection(T[1,:],gn[j,:]))
			gprime[j,:]=transvection(h0,gprime[j,:])

		end
	end
	gnew = gprime[3:nn,3:nn]
	gnidx = symplecticinverse(gnew)*numberofcosets(n)+cvw
	return gnidx
end

#Splits up the symplectic into the alpha,beta,gamma and delta arrays that
#specify the action of the clifford unitary on the X and Z Paulis respectively

function parseSymplectic(symp)
	#note that here we have (so far) ignored that for a,b,c,d we need to multiply an r and s
	#then I am guessing we just use rj as bit string on rhs. ie r1 to rn
	# so with one bit, we have and additional 2*2 * 6
	# then with two bits its 2^2 * 2^2 = 16
	s2 = round(Int,(size(symp)[1]/2))
	a = zeros(s2,s2)
	b = zeros(s2,s2)
	c = zeros(s2,s2)
	d = zeros(s2,s2)
	for i=1:s2
		for j = 1:s2
			a[i,j]=symp[(i-1)*2+1,(j-1)*2+1]
			b[i,j]=symp[(i-1)*2+1,j*2]
			c[i,j]=symp[i*2,(j-1)*2+1]
			d[i,j]=symp[i*2,j*2]
		end
	end
	return (a,b,c,d)
end

# Another way of doing this (incomplete - just an aide memoir)
# reshape(collect(Iterators.flatten(collect(zip(a,b)))),3,6)

function decomposeStateToαβγδ(state)
    n = size(state)[1]
    boundary = round(Int,n/2)
    return (state[1:boundary,1:boundary],state[1:boundary,boundary+1:end],
        state[boundary+1:end,1:boundary],state[boundary+1:end,boundary+1:end])
end

function decomposeStateToSymplectic(state)
    n = round(Int,size(state)[1]/2)
    finalD = n*2
    sympy=zeros(finalD,finalD)
    (α,β,γ,δ)= decomposeStateToαβγδ(state)
    for i=1:n
        for j = 1:n
            sympy[(i-1)*2+1,(j-1)*2+1]= α[i,j]
            sympy[(i-1)*2+1,j*2]=β[i,j]
            sympy[i*2,(j-1)*2+1]= γ[i,j]
            sympy[i*2,j*2]= δ[i,j]
        end
    end
    return sympy
end

# Takes the symplectic, parses it and uses it to create an Aaronson/Gottesman tableau
#The bits specify the 'bit' pattern that controls the sign
#For example in the 1 qubit there are only 6 unique 'symplectic' patterns
#but each of the two rows in the "tableau" can have bit signs of 0 or 1, leading
# to 6*4 different cliffords.
#this is what the bits is.

function stabiliseSymp(symp,bits)
	(a,b,c,d)=parseSymplectic(symp)

	n=size(a)[1]
	top = hcat(a,b,[ (bits >> (x-1)) %2 for x=1:n])
	for co=1:n
		bits = div(bits,2)
	end
  	bottom = hcat(c,d,[ (bits >> (x-1)) %2 for x=1:n])
 	state = convert(Array{Int32,2},vcat(top,bottom))
 end

# Takes tableau state and works out what state we are in the the 11 step deomposition
# where we are taking an arbitrary tableau state and "decomposing" it with basic gates
# back to the initial |00...00> state.
function getState(state)
	as = state[:,1:end]
	# calculate the states, backwards, so we get the "most refined" first
	#define some useful comparison matrices
	# The tableau is split up as follows:
	# A | B
	# ______
	# C | D
	#
	n=div(size(state,1),2) # half the dimension of this 2n x (2n+1) matrix
	state_i = Matrix{Int32}(I,n,n)
	state_z = zeros(Int32,n,n)
	state_r = zeros(Int32,n,1)
	if (as==vcat(hcat(state_i,state_z,state_r),hcat(state_z,state_i,state_r)))
		return 11
	end
	#split up into four matrices
	A=as[1:n,1:n]
	B=as[1:n,n+1:2*n]
	C=as[n+1:2*n,1:n]
	D=as[n+1:2*n,n+1:2*n]
	R=as[1:end,end]
	R2=as[n+1:end,end:end]
	if (B==state_z && C==state_z && R2==state_r)
		return 10
	end
	if (C==state_z && A==B  && R2==state_r)
		return 9
	end
	if (A==state_i && C==state_z && D==state_i && R2==state_r)
		return 7
	end
	if (C==state_i && D==state_z && as[n+1:2*n,end:end] == state_r)
		return 6
	end
	if (D==state_z )
		return 5
	end
	if (C==D)
		return 4
	end
	if (C==state_i)
		return 2
	end
	if (rank(C) == n)
		return 1
	end
	return 0
end

# The next series of vectors each implement one of the manipulations specified in Aaronson/Gottesman
# Section VI Canonical form.

function getFullRank(tableau::Tableau)
	n=tableau.qubits # half the dimension of this 2n x (2n+1) matrix
	svec = tableau.state
	for ii=1:n
		for i=1:n
			r1 = rank(svec[n+1:2*n,1:n])
			if r1 == n
				return
			end
			hadamard(tableau,i,false)
			r2 = rank(svec[n+1:2*n,1:n])
			if (r2 == n )
				#println("hadamard(",i,")")
				#addCommand(tableau,"hadamard($i)",Expr(:call,:hadamard,:_t,i))
				return
			end
			if (r2 > r1)
				#println("hadamard(",i,")",Expr(:call,:hadamard,:svec,i))
				#addCommand(tableau,"hadamard($i)",Expr(:call,:hadamard,:_t,i))
				continue
			else
				hadamard(tableau,i,false)
			end
		end
		hadamard(tableau,ii,false)
	end
end

function getfirstOne(myVector)
	#println(myVector)
	for i=1:length(myVector)
		if myVector[i]==1
			#println("return ",i)
			return i
		end
	end
	#print("******* THIS MAY BE AN ERROR BUT GET FIRST ONE COULDNT FIND A 1********\n")
	#print("$myVector\n")
	return 0
end

function makeCtheI(tableau::Tableau,offset=-1)
	n=tableau.qubits # half the dimension of this 2n x (2n+1) matrix
	svec = tableau.state
	sanity = 100;
	if offset==-1 offset=n # no offset means we assume its the "C" we are making Int32
	end
	allDone = false;
	while !allDone
		for i=1:n
			#print("Checking for $i\n")
			allDone = identifytheRow(tableau,offset,i)
			if allDone == false
				break
			end
		end
		sanity = sanity -1;
		if (sanity < 1)
			print("Stuck in a loop making C, I\n")
			exit()
		end
	end
	# now we need double phases to make r_n+1..r_2n to be zero
end

function zapPhase(tableau::Tableau,offset=-1)
	n=tableau.qubits # half the dimension of this 2n x (2n+1) matrix
	svec = tableau.state
	if offset==-1 offset=n # no offset means we assume its the "C" we are making Int32
	end
	for i=1:n
		if (svec[offset+i,end] ==1)
			phase(tableau,i,false)
			#addCommand(t,"phase($i)",Expr(:call,:phase,:_t,i))
			phase(tableau,i,false)
			#addCommand(t,"phase($i)",Expr(:call,:phase,:_t,i))
		end
	end
end

# What I mean is take the passed in row, and use cnots to make it an identity element.
function identifytheRow(tableau::Tableau,offset,i)
	n=tableau.qubits # half the dimension of this 2n x (2n+1) state matrix
	#Make sure the diagonal element is 1, we can only use after the element to cnot with
	#println("Checking Diagonal ",i,",",i)
	svec = tableau.state
	if svec[offset+i,i] ==0 # need to get this diagonal equal to 1
		#println("It wasn't one in the diagonal we have offset ",offset, " i ", i)
		t = getfirstOne(svec[offset+i,i+1:n])
		#println(svec[offset+i,i:n])
		#println("Got a t back of ",t)
		if t==0
			#println("Couldnt make this the identity - redoing rank\n")
			getFullRank(tableau)
			return false
		end
		t=t+i
		#println("cnot(",t,",",i,")")
		cnot(tableau,t,i,false)
		#addCommand(tableau,"cnot($t,$i)",Expr(:call,:cnot,:_t,t,i))
		#println(svec[offset+i,1:n])
	end

	for j=1:i-1
		#println("Check ",i,",",j)
		if svec[offset+i,j]!= 0
			#println("Wasnt zero")
			#addCommand(tableau,"cnot($i,$j)",Expr(:call,:cnot,:_t,i,j))
			cnot(tableau,i,j,false)
			#println(svec[offset+i,1:n])
		end
	end
	#println("Up to diagonal cleared")
	#println(svec[offset+1,1:n])
	# so the diagonal element is now 1, make the rest zero
	for j=i+1:n
		#println("Checking to make zero",i,",",j)
		if svec[offset+i,j]==1
			#println("It was one")
			#println("cnot(",i,",",j,")")
			cnot(tableau,i,j,false)
			#addCommand(tableau,"cnot($i,$j)",Expr(:call,:cnot,:_t,i,j))
			#println(svec)
		end
	end
	return true
end

function diagonaliseD(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		if svec[n+i,n+i]== 0
			phase(tableau,i,false)
			#println("phase(",i,")");
			#addCommand("phase($i)",Expr(:call,:phase,:_t,i))
		end
	end
end

function diagonaliseB(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		if svec[i,n+i]== 0
			phase(tableau,i,false)
			#println("phase(",i,")");
			#addCommand(tableau,"phase($i)",Expr(:call,:phase,:_t,i))
		end
	end
end

function cmdm(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		for j=1:n
			if svec[n+i,j] == svec[n+i,n+j]
				continue
			end
			if svec[n+i,j] == 0
				cnot(tableau,j,i,false)
				#println("cnot(",j,",",i,")")
				#addCommand("cnot($j,$i)",Expr(:call,:cnot,:_t,j,i))
			else
				phase(tableau,j,false)
				#println("phase(",j,")")
				#addCommand("phase($j)",Expr(:call,:phase,:_t,i))
			end
		end
	end
end
function ambm(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		for j=1:n
			if svec[i,j] == svec[i,n+j]
				continue
			end
			if svec[i,j] == 0
				cnot(tableau,j,i,false)
				#println("cnot(",j,",",i,")")
				#addCommand(tableau,"cnot($j,$i)",Expr(:call,:cnot,:_t,j,i))
			else
				phase(tableau,j,false)
				#println("phase(",j,")")
				#addCommand(tableau,"phase($j)",Expr(:call,:phase,:_t,j))
			end
		end
	end
end

function zapD(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		for j=1:n
			if svec[n+i,n+j]==1
				phase(tableau,j,false)
				#println("phase(",j,")")
				#addCommand("phase($j)",Expr(:call,:phase,:_t,j))
			end
		end
	end
end

function zapB(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		for j=1:n
			if svec[i,n+j]==1
				phase(tableau,j,false)
				#println("phase(",j,")")
				#addCommand(tableau,"phase($j)",Expr(:call,:phase,:_t,j))
			end
		end
	end
end

function hadamardHard(tableau::Tableau)
	n=tableau.qubits
	svec = tableau.state
	for i=1:n
		hadamard(tableau,i,false)
		#println("hadamard(",i,")")
		#addCommand(tableau,"hadamard($i)",Expr(:call,:hadamard,:_t,i))
	end
end

function makeItR0(tableau::Tableau,offset=-1)
	n=tableau.qubits
	svec = tableau.state
	if offset<0
		offset=n
	end
	for i=1:n
		if svec[offset+i,2*n+1]== 1
			phase(tableau,i,false)
			phase(tableau,i,false)
			#println("phase(",i,")")
			#println("phase(",i,")")
			#addCommand(tableau,"phase($i)",Expr(:call,:phase,:_t,i))
			#addCommand(tableau,"phase($i)",Expr(:call,:phase,:_t,i))
		end
	end
end

function nextStep(tableau::Tableau)
	n=tableau.qubits # half the dimension of this 2n x (2n+1) matrix
	svec = tableau.state

	currentState = getState(svec)
	#print(currentState)
	if (currentState == 0)
		getFullRank(tableau)
	elseif currentState == 1
		makeCtheI(tableau)
	elseif currentState < 4
		diagonaliseD(tableau)
		cmdm(tableau)
	elseif currentState == 4
		zapD(tableau)
		makeItR0(tableau)
	elseif currentState == 5
		makeCtheI(tableau)
		zapPhase(tableau)
	elseif currentState == 6
		hadamardHard(tableau)
	elseif currentState < 9
		diagonaliseB(tableau)
		ambm(tableau)
	elseif currentState == 9
		zapB(tableau)
		makeItR0(tableau,0)
	elseif currentState == 10
		makeCtheI(tableau,0) # 0 = actually its A
		zapPhase(tableau,0)
	elseif currentState == 11
		println("Done")
	end
	return
end

# maxGates is the number of gates we are allowed (an integer offset to an array)
# gates is the "stack" <- an array of gates we are applying.
# increment index 1, if its greater then maxGate, increment the one above it.
function incrementGate!(gates,maxGate,offset=1)
	if (size(gates,1) < offset )
		push!(gates,1)
		return
	end
	currentValue = gates[offset]
	currentValue = currentValue +1
	if (currentValue > maxGate)
		gates[offset]=1
		incrementGate!(gates,maxGate,offset+1)
	else
		gates[offset]=currentValue
	end
end

function bruteForceBreadthFirst(clifford)
	global svec
	count = 0
	n=div(size(clifford,1),2) # half the dimension of this 2n x (2n+1) matrix
	if (n > 4)
		println("This might take some time!, there are ",n*(n+1) + n + n, " different gates , we might need ", n*n/log(n), " of them so I *might* have to check ", n*n/log(n)^(n*(n+1) + n + n), "combos !")
		println("Ctrl-C could be your friend")
	end
	gates = (Expr)[]
	svec = []
	# If we have a n quibit system then there are a total of
	# n*(n+1) possible cnots, n phase gates and n hadamards we can apply
	for i = 1:n
	  for j = 1:n
	  	if i!=j
	  		push!(gates,Expr(:call,:cnot,:svec,i,j,false))
	  	end
	  end
	  push!(gates,draw)
	  push!(gates,Expr(:call,:phase,:svec,i,false))
	end
	# so we now have an array of all possible gates we can apply to the system
	# we want to search through the gates, we "should" need at most n^2 of them.
	# The process will be to have a stack of gates to apply, check them out
	# if we match getState(11) i.e. they reverted the clifford to the initial state
	# we have found the gates needed
	# otherwise we increment our "bottom" gate, and if we have gone though them, recursively increment the
	# gate above it (incrementGate! does this for us)
	state = 0
	maxGates = size(gates,1);
	gatesToApply=(Int32)[]
	if getState(clifford) == 11
		println("Very funny you don't need to do anything")
		return
	end
	currentAt = 1
	while (state!=11 && size(gatesToApply,1) < n*n*n*n*n*n*n*n*n*n)
		svec = copy(clifford)
		incrementGate!(gatesToApply,maxGates)
		count += 1
		if (size(gatesToApply,1) > currentAt)
			print("Trying: ")
				for i=1:size(gatesToApply,1)
				print(gatesToApply[i]," ")
			end
			println("")
			currentAt+=1
		end
		for index=1:size(gatesToApply,1)
			#println("Going to apply ",gates[gatesToApply[index]])
			eval(gates[gatesToApply[index]])
		end
		state = getState(svec)
	end

	if (state == 11)
		println("Found it after (only) $count permuations")
		println("Gates applied")
		reverse!(gatesToApply)
		for i = 1:size(gatesToApply,1)
			println(gates[gatesToApply[i]])
		end
	else
		println("Didn't find it!")
	end
	return gatesToApply
end

#coding for current bits
# 1 = phase, 2 = hadamard
# 2 + x  = cnot control (with target at x)
# 2 + n + x = cnot target (with control at x)

function checkGate(tableau::Tableau,index,currentBits,n,toDelete)
	commands = tableau.commands
	executeCommands = tableau.executeCommands

	checking = commands[index]
	m = match(r"Tableau\((.*)\)",checking)
	if (m!=nothing)
    	return
    end
    m=match(r"hadamard\((.*)\)",checking)
    if (m!=nothing)
    	bit = parse(Int,(m.captures[1]))
    	if size(currentBits[bit],1) > 0
    		if currentBits[bit][end][1] == 2
    			# two hadamards make a nothing
    			(been,gone) = pop!(currentBits[bit])
    			push!(toDelete,gone)
    			push!(toDelete,index)
			else
				push!(currentBits[bit],(2,index))
			end
		else
				push!(currentBits[bit],(2,index))
		end
		return
	end
	m=match(r"phase\((.*)\)",checking)
    if (m!=nothing) # Its a phase, we need 4 of these
    	bit = parse(Int,(m.captures[1]))
    	if size(currentBits[bit],1) > 2
    		if currentBits[bit][end][1] == 1 && currentBits[bit][end-1][1] == 1 && currentBits[bit][end-2][1] == 1
    			# four phases make a nothing.
    			(b,g) = pop!(currentBits[bit])
    			push!(toDelete,g)
    			(b,g) = pop!(currentBits[bit])
    			push!(toDelete,g)
    			(b,g) = pop!(currentBits[bit])
    			push!(toDelete,g)
    			push!(toDelete,index)
			else
				push!(currentBits[bit],(1,index))
			end
		else
				push!(currentBits[bit],(1,index))
		end
		return
	end
	m=match(r"cnot\((.*),(.*)\)",checking)
    if (m!=nothing)
    		cbit = parse(Int,(m.captures[1]))
    		tbit = parse(Int,(m.captures[2]))
    		# so just now we are going to catch two cnots in a row
    		# needs to match both bits.
    		cbitNo = 2+tbit
    		tbitNo = 2+n+cbit
    		if (size(currentBits[cbit],1) > 0 && size(currentBits[tbit],1) > 0 && currentBits[cbit][end][1] == cbitNo &&
    			currentBits[tbit][end][1] == tbitNo)
    			(b,g)=pop!(currentBits[cbit])
    			push!(toDelete,g)
    			(b,g)=pop!(currentBits[tbit])
    			# dont need to push the deletion of this gate twice.
    			push!(toDelete,index)
    		else
    			push!(currentBits[cbit],(cbitNo,index))
    			push!(currentBits[tbit],(tbitNo,index))
    		end
    end
end

function removeRedundancy(tableau::Tableau)
	n = tableau.qubits
	executeCommands = tableau.executeCommands
	commands = tableau.commands
	toDelete = (Int32)[]
	bitsOf = (Array{Tuple{Int32,Int32}})[]
	for i=1:n
    	push!(bitsOf,[])
	end
	for i=1:size(commands,1)
		checkGate(tableau,i,bitsOf,n,toDelete)
	end
	deleteat!(commands,sort!(toDelete))
	deleteat!(executeCommands,toDelete)
end

##I am going to use the symplectic method to generate all the cliffords
# - Noting that we have 24, therefore there are 6 'symplectics' and 2*2=4 different bit patterns
# - Decompose, internally generates the stabilised regime and then works out what pattern of Hadamards and Phase gates are necessary to generate this, these are stored in "commands". [Here we are in the one Qubit regime - so there are no CNOTS]
# - I use a simple loop to execute the commands using a phase gate and hadamard gate. If we were to repeat this > 1 qubit, we would need to blow up the phase gates etc (for example, if we were in 3 qubit regime and it was phase(2) i.e. phase on the second qubit, we would multiply by $I\otimes P \otimes I$
# - Better still we could define our gates (here sphase, and shadamard) to have errors and get an erroneous clifford

"""
    generateRawCliffords(;phaseGate=[1 0;0 im],hadmardGate=1/sqrt(2)*[1 1;1 -1])::Array{Complex{Float64},2}[]

returns the single qubit cliffords, generated using the passed in  phase and hadamard gates (perfect gates default)
"""
function generateRawCliffords(;phaseGate=[1 0;0 im],hadamardGate=1/sqrt(2)*[1 1;1 -1])
	straightClifford = Array{Complex{Float64},2}[]
	state=[1 0;0 1]
	for i=1:6
    	for j=1:4
        	tab = cliffordToTableau(1,i,j)
        	for t in tab.commands
           		m = match(r"initialise\((.*)\)",t)
            	if (m!=nothing)
                	state=[1 0;0 1]
            	else
                	m=match(r"phase\((.*)\)",t)
                	if (m!=nothing)
           	        	state=state*phaseGate
            	    else
                    	m=match(r"hadamard\((.*)\)",t)
                    	if (m!=nothing)
                        	state=state*hadamardGate
                    	end
                	end
            	end
        	end
        	push!(straightClifford,round.(state,digits=10))
    	end
    end
    return straightClifford
end

"""
	drawCircuit(commands)
	Takes the command string and returns their qasm equivalent.
	Can then be passed into e.g. qiskit or something similar to draw them.
"""
function drawCircuit(commands)
	qasmCommands=[]
    labels=[]
	currentDir = pwd()
	for i = 1:size(commands,1)
		m = match(r"setup\((.*)\)",commands[i])
		if (m!=nothing)
			for idx=1:Meta.parse(m.captures[1])
    			#push!(qasmCommands,"qubit q$i")
                push!(labels,"q$idx")
			end
    	else
    		m=match(r"hadamard\((.*)\)",commands[i])
    		if (m!=nothing)
    			push!(qasmCommands,("H","q$(m.captures[1])"))
    		else
    			m=match(r"phase\((.*)\)",commands[i])
    			if (m!=nothing)
                    push!(qasmCommands,("P","q$(m.captures[1])"))
    			else
    				m=match(r"cnot\((.*),(.*)\)",commands[i])
    				if (m!=nothing)
    					push!(qasmCommands,("CNOT","q$(m.captures[2])","q$(m.captures[1])"))
    				end
    			end
    		end
    	end
    end
    return (qasmCommands,labels)
end

"""
	makeFromCommand(command)

	Takes the commands in a Tableau and uses them to build a matrix representing the Tableau (in the computational basis)
	Uses the original initialise command to determine the size of the system.
	Note: for many qubits this gets quite big (scales as \$2^n\$)
"""
function makeFromCommand(command)
	up=[1;0]
	down=[0;1]
	control=[1 0; 0 0]
	target= [0 0;0 1]
	id=[1 0;0 1]
	x=[0 1;1 0]
	maxBits=1
    currentState = [1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]
    si = [1 0;0 1]
    phaseG = [1 0;0 im]
    sphase1=kron([1 0;0 im],si)
    sphase2=kron(si,[1 0;0 im])
    shadamard=1/sqrt(2)*[1 1;1 -1]
    cnot12 = [1 0 0 0;0 1 0 0;0 0 0 1;0 0 1 0]
    cnot21 = [1 0 0 0;0 0 0 1;0 0 1 0;0 1 0 0]
    for t in command
        m = match(r"initialise\((.*)\)",t)
        if (m!=nothing)
        	bit =parse(Int,(m.captures[1]))
        	maxBits=bit
            state = id;
            for i=2:bit
            	state = kron(state,id);
            end
            currentState = state
        else
           m=match(r"phase\((.*)\)",t)
           if (m!=nothing)
                bit = parse(Int,(m.captures[1]))
                tphase = phaseG
                for i=1:bit-1
                	tphase = kron(id,tphase)
				end
				for i=bit+1:maxBits
					tphase = kron(tphase,id)
				end
				currentState = currentState*tphase
            else
                m=match(r"hadamard\((.*)\)",t)
                if (m!=nothing)
                    bit = parse(Int,(m.captures[1]))
                    thad = shadamard
                    for i=1:bit-1
                    	thad=kron(id,thad)
                    end
                    for i=bit+1:maxBits
                    	thad=kron(thad,id)
                    end
                    currentState=currentState*thad
                else
                    m=match(r"cnot\((.*),(.*)\)",t)
                    if (m!=nothing)
                        cbit = parse(Int,(m.captures[1]))
                        tbit = parse(Int,(m.captures[2]))

                        # just doing this each way, its easier for my mind.
                        if cbit < tbit
                        	startC=cbit
                        	endC=tbit
                        	l1 = id
                        	l2  = x
                        	for i=cbit+1:tbit-1
                        		l1=kron(id,l1)
                        		l2=kron(id,l2)
                        	end
                        	cbitbit = kron(control,l1)+kron(target,l2)
                        else
                        	startC=tbit
                        	endC=cbit
                        	l1 = id
                        	l2  = x
                        	for i=tbit+1:cbit-1
                        		l1=kron(l1,id)
                        		l2=kron(l2,id)
                        	end
                        	cbitbit = kron(l1,control)+kron(l2,target)
                        end
                        for i = 1:startC-1
                        	cbitbit=kron(id,cbitbit)
                        end
                        for i=endC+1:maxBits
                        	cbitbit=kron(cbitbit,id)
                        end
                        currentState = currentState*cbitbit
                    end
                end
            end
        end
    end
    return  currentState
end

"""
	makeFromCommand(t::Tableau)

	Uses the commands stored in a Tableau and uses them to build a matrix representing the Tableau (in the computational basis)
	Uses the original initialise command to determine the size of the system.
	Note: for many qubits this gets quite big (scales as \$2^n\$)
"""
function makeFromCommand(t::Tableau)
	return makeFromCommand(t.commands)
end

"""
	makeInverseFromCommand(t::Tableau)

	Uses the commands stored in a Tableau and uses them to build a matrix representing the inverse of the Tableau (in the computational basis)
	Uses the original initialise command to determine the size of the system.
	This might be useful for instance if you have a Tableau represnting Pauli measurements and you wanted the Clifford that maps them to the computational basis.
	Note: for many qubits this gets quite big (scales as \$2^n\$)
"""
function makeInverseFromCommand(t::Tableau)
	return makeFromCommand(reverse(t.commands[2:end]))
end

function reverseCommands(t::Tableau)
	len = length(t.commands)
	@assert length(t.executeCommands)==len
	newCommands = []
	newExecute  = []
	for i = len:-1:1
		m=match(r"phase\((.*)\)",t.commands[i])
		if (m!=nothing)
			push!(newCommands,t.commands[i])
			push!(newCommands,t.commands[i])
			push!(newCommands,t.commands[i])
			push!(newExecute,t.executeCommands[i])
			push!(newExecute,t.executeCommands[i])
			push!(newExecute,t.executeCommands[i])
        else
        	push!(newCommands,t.commands[i])
			push!(newExecute,t.executeCommands[i])
        end
    end
	t.commands = newCommands
	t.executeCommands = newExecute
    return
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -6.148474362391325 <= q <= 10.677081718106185)

                     add_NL_constraint(m, :(softplus(-0.2518902526786948 + 0.9847866884384731*softplus(0.2752793536861313 + -0.1568397657479923*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.7649756572019379*softplus(0.9918790053549298 + 0.85981905523253*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.10087195663512549*softplus(-0.8202416682813327 + -0.17766510212211006*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + softplus(0.638435553115452 + 0.09453389081965424*softplus(0.2752793536861313 + -0.1568397657479923*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.07927014075141203*softplus(0.9918790053549298 + 0.85981905523253*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.2624391124261729*softplus(-0.8202416682813327 + -0.17766510212211006*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) + softplus(0.5273095272255199 + 0.4025366346817978*softplus(0.2752793536861313 + -0.1568397657479923*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1062303426643516*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + 0.5642457669318222*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + -0.585290488444234*softplus(0.9918790053549298 + 0.85981905523253*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + 0.9639515148457134*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.08816238542180388*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5]))) + 0.1647823489958915*softplus(-0.8202416682813327 + -0.17766510212211006*softplus(0.860587138591963 + -0.8175452953351838*$(x[1]) + -0.2911974841119127*$(x[2]) + -0.03381280038289569*$(x[3]) + -0.021805451939043596*$(x[4]) + 0.5056820468199099*$(x[5])) + -0.1793020696087071*softplus(0.40215662414311426 + 0.47796901525165936*$(x[1]) + -0.8905174133817706*$(x[2]) + 0.3815447257070974*$(x[3]) + -0.7235559998197827*$(x[4]) + -0.43821768660130234*$(x[5])) + -0.41892665263312834*softplus(-0.371366893645773 + 0.03574445051584929*$(x[1]) + 0.07916257211879385*$(x[2]) + 0.5096612381208421*$(x[3]) + 0.1651890232922466*$(x[4]) + -0.4691295251335359*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@kwdef mutable struct DSubHandler <: AbstractHandler
    id = 0
end

@register(DSubHandler, -, [DBit, DBit] => [Number])

function (handler::DSubHandler)(netlist::Netlist,
                                inputs::Vector{Variable},
                                outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic subtract does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_sub$(handler.id)
        determ_sub #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) d_sub$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_sub$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

@kwdef mutable struct DSubAHandler <: AbstractHandler
    id = 0
end

@register(DSubAHandler, -, [DBit, Number] => [Number])

function (handler::DSubAHandler)(netlist::Netlist,
                                 inputs::Vector{Variable},
                                 outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic subtract does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_suba_fxp$(handler.id)
        determ_suba_fxp #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) d_suba_fxp$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_suba_fxp$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end

@kwdef mutable struct DSubBHandler <: AbstractHandler
    id = 0
end

@register(DSubBHandler, -, [Number, DBit] => [Number])

function (handler::DSubBHandler)(netlist::Netlist,
                                 inputs::Vector{Variable},
                                 outputs::Vector{Variable})
    # check size
    if getsize(netlist, getname(inputs[1])) != (1, 1) && getsize(netlist, getname(inputs[1])) != (1, 1)
        error("Deterministic subtract does not support matrix hardware generation yet.")
    end

    outstring = """
        $stdcomment
        // BEGIN d_subb_fxp$(handler.id)
        determ_subb_fxp #(
                .BIT_WIDTH(BIT_WIDTH),
                .INT_WIDTH(INT_WIDTH)
            ) d_subb_fxp$(handler.id) (
                .a($(inputs[1].name)),
                .b($(inputs[2].name)),
                .y($(outputs[1].name))
            );
        // END d_subb_fxp$(handler.id)
        \n"""

    handler.id += 1

    return outstring
end


text:
function buildstatetestgraph()
	g = Graph()

	australia = add!(g, Vertex("Country", @compat Dict{String,Any}("name"=>"Australia","population"=>23235800)))

	northernterritory = add!(g, Vertex("Territory", @compat Dict{String,Any}("name"=>"Northern Territory","population"=>241800)))
	tasmania = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"Tasmania","population"=>513400)))
	westernaustralia = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"Western Australia","population"=>2535700)))
	queensland = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"Queensland","population"=>4676400)))
	southaustralia = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"South Australia","population"=>1674700)))
	newsouthwales = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"New South Wales","population"=>7439200)))
	victoria = add!(g, Vertex("State", @compat Dict{String,Any}("name"=>"Victoria","population"=>5768600)))
	australiancapitalterritory = add!(g, Vertex("Territory", @compat Dict{String,Any}("name"=>"Australian Capital Territory","population"=>382900)))

	add!(g, Edge("IsIn", northernterritory, australia))
	add!(g, Edge("IsIn", tasmania, australia))
	add!(g, Edge("IsIn", westernaustralia, australia))
	add!(g, Edge("IsIn", queensland, australia))
	add!(g, Edge("IsIn", southaustralia, australia))
	add!(g, Edge("IsIn", newsouthwales, australia))
	add!(g, Edge("IsIn", victoria, australia))
	add!(g, Edge("IsIn", australiancapitalterritory, australia))

	add!(g, Edge("IsAdjacent", northernterritory, westernaustralia, @compat Dict{String,Any}("direction"=>"West")))
	add!(g, Edge("IsAdjacent", northernterritory, southaustralia, @compat Dict{String,Any}("direction"=>"South")))
	add!(g, Edge("IsAdjacent", northernterritory, queensland, @compat Dict{String,Any}("direction"=>"East")))

	add!(g, Edge("IsAdjacent", westernaustralia, northernterritory, @compat Dict{String,Any}("direction"=>"East")))
	add!(g, Edge("IsAdjacent", westernaustralia, southaustralia, @compat Dict{String,Any}("direction"=>"East")))

	add!(g, Edge("IsAdjacent", southaustralia, northernterritory, @compat Dict{String,Any}("direction"=>"North")))
	add!(g, Edge("IsAdjacent", southaustralia, westernaustralia, @compat Dict{String,Any}("direction"=>"West")))
	add!(g, Edge("IsAdjacent", southaustralia, victoria, @compat Dict{String,Any}("direction"=>"East")))
	add!(g, Edge("IsAdjacent", southaustralia, newsouthwales, @compat Dict{String,Any}("direction"=>"East")))
	add!(g, Edge("IsAdjacent", southaustralia, queensland, @compat Dict{String,Any}("direction"=>"East")))

	add!(g, Edge("IsAdjacent", newsouthwales, victoria, @compat Dict{String,Any}("direction"=>"South")))
	add!(g, Edge("IsAdjacent", newsouthwales, australiancapitalterritory))
	add!(g, Edge("IsAdjacent", newsouthwales, queensland, @compat Dict{String,Any}("direction"=>"North")))
	add!(g, Edge("IsAdjacent", newsouthwales, southaustralia, @compat Dict{String,Any}("direction"=>"West")))

	add!(g, Edge("IsAdjacent", victoria, newsouthwales, @compat Dict{String,Any}("direction"=>"North")))
	add!(g, Edge("IsAdjacent", victoria, southaustralia, @compat Dict{String,Any}("direction"=>"West")))

	add!(g, Edge("IsAdjacent", queensland, newsouthwales, @compat Dict{String,Any}("direction"=>"South")))
	add!(g, Edge("IsAdjacent", queensland, northernterritory, @compat Dict{String,Any}("direction"=>"West")))
	add!(g, Edge("IsAdjacent", queensland, southaustralia, @compat Dict{String,Any}("direction"=>"West")))

	sydney = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Sydney", "population"=>4667283)))
	perth = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Perth","population"=>1897548)))
	brisbane = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Brisbane", "population"=>2189878)))
	hobart = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Hobart", "population"=>282099)))
	canberra = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Canberra", "population"=>382900)))
	darwin = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Darwin", "population"=>131678)))
	melbourne = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Melbourne","population"=>4246345)))
	adelaide = add!(g, Vertex("City", @compat Dict{String,Any}("name"=>"Adelaide","population"=>1277174)))

	add!(g, Edge("IsCapitalOf", sydney, newsouthwales))
	add!(g, Edge("IsCapitalOf", perth, westernaustralia))
	add!(g, Edge("IsCapitalOf", brisbane, queensland))
	add!(g, Edge("IsCapitalOf", hobart, tasmania))
	add!(g, Edge("IsCapitalOf", darwin, northernterritory))
	add!(g, Edge("IsCapitalOf", melbourne, victoria))
	add!(g, Edge("IsCapitalOf", canberra, australiancapitalterritory))
	add!(g, Edge("IsCapitalOf", canberra, australia))
	add!(g, Edge("IsCapitalOf", adelaide, southaustralia))

	return g
end


text:
# Main module file

__precompile__()

"Thermophysical properties of water"
module WaterData

# Include all other parts of the module
include("config.jl")
include("constants.jl")
include("util.jl")
include("regions.jl")
include("phaseboundaries.jl")
include("eos.jl")
include("heatcapacity.jl")
include("thermalexpansivity.jl")
include("io.jl")

include("testhelpers.jl")

end # module WaterData


text:
using Dates: default_format

function Sink(fullpath::Union{AbstractString, IO};
              delim::Char=',',
              quotechar::Char='"',
              escapechar::Char='\\',
              missingstring::AbstractString="",
              dateformat=nothing,
              header::Bool=true,
              colnames::Vector{String}=String[],
              append::Bool=false,
              quotefields::Bool=false)
    delim = delim % UInt8; quotechar = quotechar % UInt8; escapechar = escapechar % UInt8
    dateformat = isa(dateformat, AbstractString) ? Dates.DateFormat(dateformat) : dateformat
    io = IOBuffer()
    options = CSV.Options(delim=delim, quotechar=quotechar, escapechar=escapechar, missingstring=missingstring, dateformat=dateformat)
    !append && header && !isempty(colnames) && writeheaders(io, colnames, options, Val{quotefields})
    return Sink(options, io, fullpath, position(io), !append && header && !isempty(colnames), colnames, length(colnames), append, Val{quotefields})
end

quoted(::Type{Val{true}},  val, q, e, d) =  string(q, replace(val, q => string(e, q)), q)
quoted(::Type{Val{false}}, val, q, e, d) = (q in val || d in val) ? string(q, replace(val, q => string(e, q)), q) : val

function writeheaders(io::IOBuffer, h::Vector{String}, options, quotefields)
    cols = length(h)
    q = Char(options.quotechar); e = Char(options.escapechar); d = Char(options.delim)
    for col = 1:cols
        Base.write(io, quoted(quotefields, h[col], q, e, d), ifelse(col == cols, NEWLINE, options.delim))
    end
    return nothing
end

# DataStreams interface
Data.streamtypes(::Type{CSV.Sink}) = [Data.Field]
Data.weakrefstrings(::Type{CSV.Sink}) = true

# Constructors
function Sink(sch::Data.Schema, T, append, file::Union{AbstractString, IO}; reference::Vector{UInt8}=UInt8[], kwargs...)
    sink = Sink(file; append=append, colnames=Data.header(sch), kwargs...)
    return sink
end

function Sink(sink, sch::Data.Schema, T, append; reference::Vector{UInt8}=UInt8[])
    sink.append = append
    sink.cols = size(sch, 2)
    !sink.header && !append && writeheaders(sink.io, Data.header(sch), sink.options, sink.quotefields)
    return sink
end

function Data.streamto!(sink::Sink, ::Type{Data.Field}, val, row, col::Int)
    q = Char(sink.options.quotechar); e = Char(sink.options.escapechar); d = Char(sink.options.delim)
    Base.write(sink.io, quoted(sink.quotefields, string(val), q, e, d), ifelse(col == sink.cols, NEWLINE, d))
    return nothing
end

function Data.streamto!(sink::Sink, ::Type{Data.Field}, val::Dates.TimeType, row, col::Int)
    v = Dates.format(val, sink.options.dateformat === nothing ? default_format(typeof(val)) : sink.options.dateformat)
    Base.write(sink.io, v, ifelse(col == sink.cols, NEWLINE, sink.options.delim))
    return nothing
end

const EMPTY_UINT8_ARRAY = UInt8[]
function Data.streamto!(sink::Sink, ::Type{Data.Field}, val::Missing, row, col::Int)
    Base.write(sink.io, sink.options.missingcheck ? sink.options.missingstring : EMPTY_UINT8_ARRAY, ifelse(col == sink.cols, NEWLINE, sink.options.delim))
    return nothing
end

function Data.close!(sink::CSV.Sink)
    io = isa(sink.fullpath, AbstractString) ? open(sink.fullpath, sink.append ? "a" : "w") : sink.fullpath
    Base.write(io, take!(sink.io))
    isa(sink.fullpath, AbstractString) && close(io)
    return sink
end

"""
`CSV.write(file_or_io::Union{AbstractString,IO}, source::Type{T}, args...; kwargs...)` => `CSV.Sink`

`CSV.write(file_or_io::Union{AbstractString,IO}, source::Data.Source; kwargs...)` => `CSV.Sink`

write a `Data.Source` out to a `file_or_io`.

Positional Arguments:

* `file_or_io`; can be a file name (string) or other `IO` instance
* `source` can be the *type* of `Data.Source`, plus any required `args...`, or an already constructed `Data.Source` can be passsed in directly (2nd method)

Keyword Arguments:

* `delim::Union{Char,UInt8}`; how fields in the file will be delimited; default is `UInt8(',')`
* `quotechar::Union{Char,UInt8}`; the character that indicates a quoted field that may contain the `delim` or newlines; default is `UInt8('"')`
* `escapechar::Union{Char,UInt8}`; the character that escapes a `quotechar` in a quoted field; default is `UInt8('\\')`
* `missingstring::String`; the ascii string that indicates how missing values will be represented in the dataset; default is the empty string `""`
* `dateformat`; how dates/datetimes will be represented in the dataset; default is ISO-8601 `yyyy-mm-ddTHH:MM:SS.s`
* `header::Bool`; whether to write out the column names from `source`
* `colnames::Vector{String}`; a vector of string column names to be used when writing the header row
* `append::Bool`; start writing data at the end of `io`; by default, `io` will be reset to the beginning or overwritten before writing
* `transforms::Dict{Union{String,Int},Function}`; a Dict of transforms to apply to values as they are parsed. Note that a column can be specified by either number or column name.

A few example invocations include:
```julia
# write out a DataFrame `df` to a file name "out.csv" with all defaults, including comma as delimiter
CSV.write("out.csv", df)

# write out a DataFrame, this time as a tab-delimited file
CSV.write("out.csv", df; delim='\t')

# write out a DataFrame, with missing values represented by the string "NA"
CSV.write("out.csv", df; missingstring="NA")

# write out a "header-less" file, with actual data starting on row 1
CSV.write("out.csv", df; header=false)

# write out a DataFrame `df` twice to a file, the resulting file with have twice the # of rows as the DataFrame
# note the usage of the keyword argument `append=true` in the 2nd call
CSV.write("out.csv", df)
CSV.write("out.csv", df; append=true)

# write a DataFrame out to an IOBuffer instead of a file
io = IOBuffer
CSV.write(io, df)

# write the result of an SQLite query out to a comma-delimited file
db = SQLite.DB()
sqlite_source = SQLite.Source(db, "select * from sqlite_table")
CSV.write("sqlite_table.csv", sqlite_source)
```
"""
function write end

function write(file::Union{AbstractString, IO}, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...) where {T}
    sink = Data.stream!(T(args...), CSV.Sink, file; append=append, transforms=transforms, kwargs...)
    return Data.close!(sink)
end
function write(file::Union{AbstractString, IO}, source; append::Bool=false, transforms::Dict=Dict{Int,Function}(), kwargs...)
    sink = Data.stream!(source, CSV.Sink, file; append=append, transforms=transforms, kwargs...)
    return Data.close!(sink)
end

write(sink::Sink, ::Type{T}, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(T(args...), sink; append=append, transforms=transforms); return Data.close!(sink))
write(sink::Sink, source; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))


text:
export diageye

import LinearAlgebra: Diagonal

xT_A_x(x::AbstractVector, A::AbstractMatrix) = dot(x, A * x)

diageye(::Type{T}, n::Int) where { T <: Real } = Matrix{T}(I, n, n)
diageye(n::Int)                                = diageye(Float64, n)


text:
using Pkg
using CSV, DataFrames
using Plots, Gadfly, Cairo
theme(:vibrant); # COlor schema
#https://docs.juliaplots.org/latest/generated/plotthemes/
using Statistics, Distributions
using Images, FileIO

str_cd = "C:\\Users\\justi\\Documents\\GitHub\\umich-robinson-plays\\eda\\"
str_cd_dfs = "C:\\Users\\justi\\Documents\\GitHub\\umich-robinson-plays\\data\\"

# ===================================================================================================================
# ===================================================================================================================
# ===================================================================================================================
# ========================================== Duncan Robinson ========================================================
# ===================================================================================================================
# ===================================================================================================================
# ===================================================================================================================
df = CSV.read(str_cd_dfs * "data_duncan_robinson_plays.csv", DataFrame);
# Transform  x - coordinate (long court) symmetric view
df[:, "event_coord_x"] = abs.(df[:, "event_coord_x"] .- 564);
df[df[:, :three_point_shot] .== true, :three_point_shot];

# Exploratory insights
#combine(groupby(df, [:type]), :type => length => :count)
#combine(groupby(df, [:game_id]), :type => length => :count)

df_sub = df[(df.type .== "fieldgoal"), :];
df_sub = df[df[:, :three_point_shot] .== true, :];

# Figure out three-point attempts

# Start with all fieldgoals
df_sub_3pt = df_sub[(df_sub.type .== "fieldgoal"), :]
# COllect three-point attemps and scores
#df_sub_3pt = df_sub[contains.(df_sub[:, "event_description"], "three"), :];
df_sub_3pt[:, "pointmade"] .= 0;
#df_sub_3pt[contains.(df_sub_3pt[:, "event_description"], "makes"), "pointmade"] .= 1;
df_sub_3pt[.!ismissing.(df_sub_3pt.points_scored .== 3), "pointmade"] .= 1;

# Full court
xtiksall = [-120, 525, 1140]
labelsxall = Dict(zip(xtiksall, ["baseline", "midcourt", "baseline"]))
# Half court
xtiks = [100, 525]
labelsx = Dict(zip(xtiks, ["midcourt", "baseline"]))

ytiks = [-50, 275, 580]
labelsy = Dict(zip(ytiks, ["right", "center", "left"]))

df_plt_3ptyes = df_sub_3pt[(df_sub_3pt.pointmade .== 1), ["event_coord_x", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptyes[:, 1], y=df_plt_3ptyes[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                      Coord.cartesian(xmin=xtiks[1], xmax=xtiks[2], ymin=ytiks[1], ymax=ytiks[3]),
                      Guide.xticks(ticks = xtiks), Guide.xlabel("Half-court (standardized)"),
                      Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                      Scale.x_continuous(labels = x -> labelsx[x]),
                      Scale.y_continuous(labels = x -> labelsy[x]),
                      #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                      Theme(background_color = nothing))
img = PNG(str_cd * "plt_3pt_scores_robinson.png", 18cm, 18cm)
draw(img, plt_3ptattempt)

df_plt_3ptatmt = df_sub_3pt[:, ["event_coord_x", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptatmt[:, 1], y=df_plt_3ptatmt[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                              Coord.cartesian(xmin=xtiks[1], xmax=xtiks[2], ymin=ytiks[1], ymax=ytiks[3]),
                              Guide.xticks(ticks = xtiks), Guide.xlabel("Half-court (standardized)"),
                              Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                              Scale.x_continuous(labels = x -> labelsx[x]),
                              Scale.y_continuous(labels = x -> labelsy[x]),
                              #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                              Theme(background_color = nothing))
img = PNG(str_cd * "plt_3pt_attempt_robinson.png", 18cm, 16cm)
draw(img, plt_3ptattempt)

# ===================================================================================================================
# ===================================================================================================================
# ===================================================================================================================
# ========================================== NCAA Statistics ========================================================
# ===================================================================================================================
# ===================================================================================================================
# ===================================================================================================================
df_all = CSV.read(str_cd_dfs * "data_all_plays.csv", DataFrame);

df_all[:, "event_coord_x_half"] = abs.(df_all[:, "event_coord_x"] .- 564);
# Data cleaning process

df_all[:, "event_coord_x_adj"] .= round.(Int, df_all.event_coord_x./10)*10;
df_all[:, "event_coord_y_adj"] .= round.(Int, df_all.event_coord_y./10)*10;
# Mark if points are made
df_all[:, "pointmade"] .= 0;
df_all[.!ismissing.(df_all.points_scored .== 3), "pointmade"] .= 1;

df_all_adj = combine(groupby(df_all, ["event_coord_x_adj", "event_coord_y_adj"]), :three_point_shot => sum => :temp_num_attempt);
# Evaluate what coordinates are giving min value issues
describe(df_all_adj, :min, :max, :q25, :median, cols = :temp_num_attempt);
# Select those aggregated coordinate that does not have much of attempts
df_row_remove = df_all_adj[df_all_adj[:, :temp_num_attempt] .>= 7, :];

df_all_filtered = rightjoin(df_all, df_row_remove, on = [:event_coord_x_adj, :event_coord_y_adj]);
df_all_filteredout = leftjoin(df_all, df_row_remove, on = [:event_coord_x_adj, :event_coord_y_adj]);

df_sub_all = df_all_filtered[df_all_filtered[:, :three_point_shot] .== true, :];

#Explanatory plot
df_plt_3ptyes = df_sub_all[(df_sub_all.pointmade .== 1), ["event_coord_x", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptyes[:, 1], y=df_plt_3ptyes[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                      Coord.cartesian(xmin=xtiksall[1], xmax=xtiksall[3], ymin=ytiks[1], ymax=ytiks[3]),
                      Guide.xticks(ticks = xtiksall), Guide.xlabel("Full court"),
                      Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                      Scale.x_continuous(labels = x -> labelsxall[x]),
                      Scale.y_continuous(labels = x -> labelsy[x]),
                      #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                      Theme(background_color = nothing))
img = PNG(str_cd * "plt_3pt_scores_all_fullcourt_ncaa.png", 29cm, 15cm)
draw(img, plt_3ptattempt)

#Alternate option to be save with image

img_path = str_cd*"basketball_halfcourt.png"
img = load(img_path)

df_plt_3ptatmt = df_sub_all[:, ["event_coord_x_half", "event_coord_y"]];
plt_img = Plots.plot([0, 535], [-20, 620], reverse(img, dims = 1), yflip = false, aspect_ratio = :none, showaxis = false)
histogram2d!(df_plt_3ptatmt[:, 1], df_plt_3ptatmt[:, 2], bins = 65, showaxis = false, ticks = false, normed = false, aspect_ratio = 1)
annotate!([(-10,570,text("Three-pooint attempts", 13, :left, :top, :blue))])
savefig(plt_img, str_cd * "shot_attempts.png") # save the fig referenced by plot_ref as fn

df_plt_3ptyes = df_sub_all[(df_sub_all.pointmade .== 1), ["event_coord_x_half", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptyes[:, 1], y=df_plt_3ptyes[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                      Coord.cartesian(xmin=xtiks[1], xmax=xtiks[2], ymin=ytiks[1], ymax=ytiks[3]),
                      Guide.xticks(ticks = xtiks), Guide.xlabel("Half-court (standardized)"),
                      Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                      Scale.x_continuous(labels = x -> labelsx[x]),
                      Scale.y_continuous(labels = x -> labelsy[x]),
                      #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                      Theme(background_color = nothing))
img = PNG(str_cd * "plt_3pt_scores_all_ncaa.png", 18cm, 18cm)
draw(img, plt_3ptattempt)

df_plt_3ptatmt = df_sub_all[:, ["event_coord_x_half", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptatmt[:, 1], y=df_plt_3ptatmt[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                              Coord.cartesian(xmin=xtiks[1], xmax=xtiks[2], ymin=ytiks[1], ymax=ytiks[3]),
                              Guide.xticks(ticks = xtiks), Guide.xlabel("Half-court (standardized)"),
                              Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                              Scale.x_continuous(labels = x -> labelsx[x]),
                              Scale.y_continuous(labels = x -> labelsy[x]),
                              #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                              Theme(background_color = nothing))
img = PNG(str_cd * "plt_3pt_attempt_all_ncaa.png", 18cm, 16cm)
draw(img, plt_3ptattempt)

#=
# THose shots that are made across the half court
condition_1 = (((df_all.team_basket .== "left") .& (df_all.event_coord_x .> 564)) .|
                               ((df_all.team_basket .== "right") .& (df_all.event_coord_x .< 564)))

suspected_hailmarry =  df_all[condition_1, :]

suspected_hailmarry[:, "pointmade"] .= 0;
suspected_hailmarry[.!ismissing.(suspected_hailmarry.points_scored .== 3), "pointmade"] .= 1;

# Start with all fieldgoals
suspected_hailmarry = suspected_hailmarry[(suspected_hailmarry.type .== "fieldgoal"), :]

df_plt_3ptyes = suspected_hailmarry[(suspected_hailmarry.pointmade .== 1), ["event_coord_x", "event_coord_y"]];
plt_3ptattempt = Gadfly.plot(x=df_plt_3ptyes[:, 1], y=df_plt_3ptyes[:, 2], Geom.hexbin(xbincount=20, ybincount=20),
                      Coord.cartesian(xmin=xtiksall[1], xmax=xtiksall[3], ymin=ytiks[1], ymax=ytiks[3]),
                      Guide.xticks(ticks = xtiksall), Guide.xlabel("Full court"),
                      Guide.yticks(ticks = ytiks), Guide.ylabel(""),
                      Scale.x_continuous(labels = x -> labelsxall[x]),
                      Scale.y_continuous(labels = x -> labelsy[x]),
                      #Guide.xticks(ticks=nothing, label = false), Guide.yticks(ticks=nothing),
                      Theme(background_color = nothing))

df_all[((((df_all.team_basket .== "left") .& (df_all.event_coord_x .> 564)) .|
        ((df_all.team_basket .== "right") .& (df_all.event_coord_x .< 564))) .&
        df_all.three_point_shot .== true) .|
        (abs.(df_all.event_coord_x .- 525) .< 100), :event_coord_x] .= 525
=#

max_stars_repo_name: justinjoliver/umich-robinson-plays

max_issues_repo_name: justinjoliver/umich-robinson-plays

max_forks_repo_name: justinjoliver/umich-robinson-plays


text:
module RecursiveFiltering

using KernelDensityEstimate, TransformUtils

export
  SysState,
  filterpropagate!,
  filterupdate!

include("KalmanFilter.jl")
include("KDEFilter.jl")

end


text:
#!/usr/bin/env julia
using GoTMSupport
using Base.Test

gopath = ENV["GOPATH"]

if isempty(gopath)
    error("You must set the GOPATH environment variable for this test to work")
end

choices_plist = completion_choices_plist("example.go", 37, 7, gopath * "/bin/gocode")
print(choices_plist)

# If the test works you should see something like this:
#
# ({display = "Add(u Vector2D) Vector2D";match = Add;insert = "(${1:u})$0";image = func;}, {display = "Dot(u Vector2D) float64";match = Dot;insert = "(${1:u})$0";image = func;}, {display = "Normal() Vector2D";match = Normal;insert = "()$0";image = func;}, {display = "X float64";match = X;insert = "$0";image = var;}, {display = "Y float64";match = Y;insert = "$0";image = var;})


text:
#####
##### 2D modes on triangles
#####

"""
    simplex_2D(a, b, i, j)

Evaluate 2D PKDO basis phi_ij at points (a,b) on the Duffy domain [-1,1]^2
"""
function simplex_2D(a, b, i, j)
    h₁ = jacobiP(a, 0, 0, i)
    h₂ = jacobiP(b, 2i+1, 0, j)
    return @. sqrt(2.0)*h₁*h₂*(1-b)^i # WARNING: (1-b)^i can blow up if N is too large (> 25)
end

"""
    grad_simplex_2D(a, b, id, jd)

Evalute the partial derivative w.r.t. (r,s) of 2D PKDO polynomials with
indices (id,jd) at (a,b)
"""

function grad_simplex_2D(a, b, id, jd)
    fa = jacobiP(a, 0, 0, id)
    gb = jacobiP(b, 2*id+1, 0, jd)
    dfa = grad_jacobiP(a, 0, 0, id)
    dgb = grad_jacobiP(b, 2*id+1, 0, jd)

    dmodedr = @. dfa*gb
    if id > 0
        dmodedr = dmodedr.*((0.5*(1 .- b)).^(id-1))
    end

    dmodeds = @. dfa*(gb*(0.5*(1+a)))
    if id > 0
        dmodeds = dmodeds.*((0.5*(1 .- b)).^(id-1))
    end

    tmp = @. dgb*((0.5*(1-b))^id)
    if id > 0
        tmp = tmp-0.5*id*gb.*((0.5*(1 .- b)).^(id-1))
    end
    dmodeds = dmodeds+fa.*tmp

    dmodedr = 2^(id+0.5)*dmodedr
    dmodeds = 2^(id+0.5)*dmodeds
    return dmodedr, dmodeds
end

"""
    rstoab(r, s, tol = 1e-12)

Converts from reference bi-unit right triangle coordinate (r,s) to polynomial
basis evaluation coordinates (a,b) on the domain [-1,1]^2
"""
function rstoab(r, s, tol = 1e-12)
    a = zeros(length(r))
    for n = 1:length(r)
        if abs(s[n] - 1) > tol
            a[n] = 2*(1+r[n])/(1-s[n])-1
        else
            a[n] = -1
        end
    end
    return a, s
end

"""
    quad_nodes_tri(N)

Returns quadrature nodes (from Gimbutas and Xiao 2010) which exactly integrate degree N polynomials
"""

function quad_nodes_tri(N)

    if N < 28
        rsw::Matrix{Float64} = readdlm(string(@__DIR__,"/QuadratureData/quad_nodes_tri_N", N, ".txt"),' ', Float64, '\n')
        r = rsw[:,1]
        s = rsw[:,2]
        w = rsw[:,3]
    else
        cubN = convert(Int,ceil((N+1)/2))
        r,s,w = stroud_quad_nodes(Tri(),cubN)
    end

    return r[:], s[:], w[:]
end

function stroud_quad_nodes(elem::Tri,N)
    cubA,cubWA = gauss_quad(0,0,N)
    cubB,cubWB = gauss_quad(1,0,N)

    cubA,cubB = vec.(meshgrid(cubA,cubB))

    r = @. 0.5*(1+cubA)*(1-cubB)-1
    s = cubB
    w = 0.5*cubWB*(cubWA')
    return vec.((r,s,w))
end

function equi_nodes(elem::Tri,N)
    Np = convert(Int,(N+1)*(N+2)/2)

    r = zeros(Np)
    s = zeros(Np)

    r1D = LinRange(-1,1,N+1)
    sk = 1
    for i = 0:N
        for j = 0:N-i
            r[sk] = r1D[i+1]
            s[sk] = r1D[j+1]
            sk += 1
        end
    end

    return r[:], s[:]
end

function quad_nodes(elem::Tri,N)
    return quad_nodes_tri(2*N)
end

function basis(elem::Tri,N,r,s)
    Np = convert(Int,(N+1)*(N+2)/2)
    V2D,V2Dr,V2Ds = ntuple(x->zeros(length(r), Np),3)
    a, b = rstoab(r, s)
    sk = 1
    for i = 0:N
        for j = 0:N-i
            V2D[:,sk] = simplex_2D(a,b,i,j)
            V2Dr[:,sk], V2Ds[:,sk] = grad_simplex_2D(a, b, i, j)
            sk = sk+1
        end
    end
    return V2D,V2Dr,V2Ds
end


text:
# Deal with multiple-source nodes.
function build_sync(s::Species, b::FastBreeder)

  # For each source used by this breeder, build the preceding set of
  # synchronisation operations and calculate the final state of the
  # representation graph for each, before calculating the synchronisation
  # operations required to produce a complete representation graph for
  # each.
  for child in sources(b)
    rep_graph = build_sync(RepresentationGraph(s), child)
    sync_path = species.full_sync!(rep_graph)
    build_sync(child.eigen, b.eigen, sync_path)
  end

  # Return the built breeder.
  return b
end

build_sync(g::RepresentationGraph, n::SelectionBreederSource) =
  species.initialize!(species.clone(g))

function build_sync(g::RepresentationGraph, n::VariationBreederSource)
  rep_graph::RepresentationGraph = g

  # For each source used by this breeding operator, calculate the shortest chain of
  # synchronisation operations required to generate the stage used by this operator.
  for child in sources(n)
    rep_graph = build_sync(g, child)
    sync_path = species.partial_sync!(rep_graph, n.stage_name)

    # If it's multi-source?
    # Mark all nodes except the edited node as dirty; DODGY.

    # Fetch the eigentypes of the two breeding nodes and implement the synchronisation
    # operation between them.
    build_sync(child.eigen, n.eigen, sync_path)
  end

  # Return the resulting representation graph.
  return rep_graph
end

function build_sync(from::Type, to::Type, ops::Vector{Tuple{AbstractString, AbstractString}})
  body = "function sync{I <: Individual}(from::Type{$(from)}, to::Type{$(to)}, s::Species{I}, inds::Vector{I})\n" *
    join(["convert!(s, \"$(s_from)\", \"$(s_to)\", inds)" for (s_from, s_to) in ops], "\n") *
    "\nreturn inds" *
    "\nend"

  # Dynamically evaluate the function definition
  eval(breeder, Base.parse(body))
end


text:
pages = Any[
        "Home" => "index.md",
        "Tutorials" => Any[
            "tutorials/using_catalyst.md",
            "tutorials/dsl.md",
            "tutorials/reaction_systems.md",
            "tutorials/basic_examples.md",
            "tutorials/compositional_modeling.md",
            "tutorials/symbolic_stoich.md",
            "tutorials/bifurcation_diagram.md",
            "tutorials/parameter_estimation.md",
            "tutorials/generating_reactions_programmatically.md"
        ],
        "FAQs" => "faqs.md",
        "API" => Any[
            "api/catalyst_api.md"
        ]
    ]


text:
"""
A simple example of creating a persistent session and logging into a web form. HTTP.jl does not have a distinct session object like requests.session() or rvest::html_session() but rather uses the `cookies` flag along with standard functions
"""
using HTTP

#dummy site, any credentials work
url = "http://quotes.toscrape.com/login"
session = HTTP.get(url; cookies = true)

credentials = Dict(
    "Username" => "username",
    "Password" => "password")

response = HTTP.post(url, credentials)


text:
using TimeSeriesIO

using Base.Test

include("test_dataframe.jl")
include("test_timearray.jl")


text:
struct Quadratic{BC<:BoundaryCondition} <: DegreeBC{2}
    bc::BC
end
(deg::Quadratic)(gt::GridType) = Quadratic(deg.bc(gt))

"""
Assuming uniform knots with spacing 1, the `i`th piece of quadratic spline
implemented here is defined as follows:

    y_i(x) = cm p(x-i) + c q(x) + cp p(1-(x-i))

where

    p(δx) = (δx - 1)^2 / 2
    q(δx) = 3/4 - δx^2

and the values for `cX` for `X ∈ {m,_,p}` are the pre-filtered coefficients.

For future reference, this expands to the following polynomial:

    y_i(x) = cm * 1/2 * (x-i-1)^2 + c * (3/4 - x + i)^2 + cp * 1/2 * (x-i)^2

When we derive boundary conditions we will use derivatives `y_1'(x-1)` and
`y_1''(x-1)`
"""
Quadratic

function positions(deg::Quadratic, ax, x)
    xm = roundbounds(x, ax)
    δx = x - xm
    expand_index(deg, fast_trunc(Int, xm), ax, δx), δx
end

value_weights(::Quadratic, δx) = (
    sqr(δx - SimpleRatio(1,2))/2,
    SimpleRatio(3,4) - sqr(δx),
    sqr(δx + SimpleRatio(1,2))/2)

gradient_weights(::Quadratic, δx) = (
    δx - SimpleRatio(1,2),
    -2 * δx,
    δx + SimpleRatio(1,2))

hessian_weights(::Quadratic, δx) = (oneunit(δx), -2*oneunit(δx), oneunit(δx))

expand_index(::Quadratic, xi::Number, ax::AbstractUnitRange, δx) = xi-1  # uses WeightedAdjIndex
# Others use WeightedArbIndex
expand_index(::Quadratic{<:Periodic}, xi::Number, ax::AbstractUnitRange, δx) =
    (modrange(xi-1, ax), modrange(xi, ax), modrange(xi+1, ax))
expand_index(::Quadratic{BC}, xi::Number, ax::AbstractUnitRange, δx) where BC<:Union{InPlace,InPlaceQ} =
    (max(xi-1, first(ax)), xi, min(xi+1, last(ax)))

padded_axis(ax::AbstractUnitRange, ::BSpline{<:Quadratic}) = first(ax)-1:last(ax)+1
padded_axis(ax::AbstractUnitRange, ::BSpline{Quadratic{BC}}) where BC<:Union{Periodic,InPlace,InPlaceQ} = ax

# # Due to padding we can extend the bounds
# lbound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = first(ax) - 0.5
# ubound(ax, ::BSpline{Quadratic{BC}}, ::OnGrid) where BC = last(ax) + 0.5

function inner_system_diags(::Type{T}, n::Int, ::Quadratic) where {T}
    du = fill(convert(T, SimpleRatio(1,8)), n-1)
    d = fill(convert(T, SimpleRatio(3,4)), n)
    dl = copy(du)
    (dl,d,du)
end

"""
`Quadratic{Flat}` `OnCell` and `Quadratic{Reflect}` `OnCell` amounts to setting
`y_1'(x) = 0` at x=1/2. Applying this condition yields

    -cm + c = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{BC}) where {T,TC,BC<:Union{Flat{OnCell},Reflect{OnCell}}}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = -1
    du[1] = dl[end] = 1
    lut!(dl, d, du), zeros(TC, n)
end

function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{InPlace{OnCell}}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = convert(T, SimpleRatio(7,8))
    lut!(dl, d, du), zeros(TC, n)
end

# InPlaceQ continues the quadratic at 2 all the way down to 1 (rather than 1.5)
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{InPlaceQ{OnCell}}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = SimpleRatio(9,8)
    dl[end] = du[1] = SimpleRatio(-1,4)
    # Woodbury correction to add 1/8 for row 1, col 3 and row n, col n-2
    rowspec = spzeros(T, n, 2)
    colspec = spzeros(T, 2, n)
    valspec = zeros(T, 2, 2)
    valspec[1,1] = valspec[2,2] = SimpleRatio(1,8)
    rowspec[1,1] = rowspec[end,2] = 1
    colspec[1,3] = colspec[2,end-2] = 1
    Woodbury(lut!(dl, d, du), rowspec, valspec, colspec), zeros(TC, n)
end

"""
`Quadratic{Flat}` `OnGrid` and `Quadratic{Reflect}` `OnGrid` amount to setting
`y_1'(x) = 0` at `x=1`. Applying this condition yields

    -cm + cp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{BC}) where {T,TC,BC<:Union{Flat{OnGrid},Reflect{OnGrid}}}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = -1
    du[1] = dl[end] = 0

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, oneunit(T)),
                                  (n, n-2, oneunit(T))
                                 )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Quadratic{Line}` `OnGrid` and `Quadratic{Line}` `OnCell` amount to setting
`y_1''(x) = 0` at `x=1` and `x=1/2` respectively. Since `y_i''(x)` is independent
of `x` for a quadratic b-spline, these both yield

    1 cm -2 c + 1 cp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{<:Line}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = 1
    du[1] = dl[end] = -2

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, 3, oneunit(T)),
                                  (n, n-2, oneunit(T)),
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Quadratic{Free}` `OnGrid` and `Quadratic{Free}` `OnCell` amount to requiring
an extra continuous derivative at the second-to-last cell boundary; this means
that `y_1''(3/2) = y_2''(3/2)`, yielding

    1 cm -3 c + 3 cp - cpp = 0
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{<:Free}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)
    d[1] = d[end] = 1
    du[1] = dl[end] = -3

    specs = WoodburyMatrices.sparse_factors(T, n,
                                    (1, 3, 3),
                                    (1, 4, -1),
                                    (n, n-2, 3),
                                    (n, n-3, -1))

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end

"""
`Quadratic{Periodic}` `OnGrid` and `Quadratic{Periodic}` `OnCell` close the system
by looking at the coefficients themselves as periodic, yielding

    c0 = c(N+1)

where `N` is the number of data points.
"""
function prefiltering_system(::Type{T}, ::Type{TC}, n::Int, degree::Quadratic{<:Periodic}) where {T,TC}
    dl,d,du = inner_system_diags(T,n,degree)

    specs = WoodburyMatrices.sparse_factors(T, n,
                                  (1, n, du[1]),
                                  (n, 1, dl[end])
                                  )

    Woodbury(lut!(dl, d, du), specs...), zeros(TC, n)
end


text:
#=
A unit fraction contains 1 in the numerator. The decimal representation of the unit fractions with denominators 2 to 10 are given:

1/2	= 	0.5
1/3	= 	0.(3)
1/4	= 	0.25
1/5	= 	0.2
1/6	= 	0.1(6)
1/7	= 	0.(142857)
1/8	= 	0.125
1/9	= 	0.(1)
1/10	= 	0.1
Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be seen that 1/7 has a 6-digit recurring cycle.

Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its decimal fraction part.
=#

# Multiplicative order
# http://mathworld.wolfram.com/MultiplicativeOrder.html
function multiplicative_order(a, n)
  if gcd(a,n) > 1
    return 0
  else
    order = 1
    mod_exp = a
    while mod_exp != 1
      order +=1
      mod_exp = (mod_exp * a) % n
    end
    return order
  end
end

function calc()
  largest = 0
  max_i = 0
  for i in 2:1000
    a = multiplicative_order(10, i)
    if a > largest
      largest = a
      max_i = i
    end
  end
  max_i
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
"""
Return leadfield for target location
"""
function calculate_specific_leadfield2(ldf::Leadfield, target::Coordinate; kwargs...)

    calculate_specific_leadfield2(ldf, find_location(ldf, target); kwargs...)
end

"""
Return leadfield for target location index
"""
function calculate_specific_leadfield2(ldf::Leadfield, target::Int; kwargs...)

    calculate_specific_leadfield2(ldf.L, target; kwargs...)
end

"""
Return leadfield for target location with null placed at supress location
"""
function calculate_specific_leadfield2(ldf::Leadfield, target::Union{Coordinate, Int}, supress::Union{Coordinate, Int}; bilateral::Bool=false, kwargs...)

    if bilateral & isa(supress, Coordinate)
        supress.x = -1 * supress.x
    elseif bilateral & isa(supress, Int)
        supress = find_location(ldf, Talairach(-1 * ldf.x[supress], ldf.y[supress], ldf.z[supress]))
    end

    target_ldf = calculate_specific_leadfield2(ldf, target; kwargs...)
    supress_ldf = calculate_specific_leadfield2(ldf, supress; kwargs...)

    hcat(target_ldf, supress_ldf)
end

"""
Return leadfield for target location with nulls placed in supress region
"""
function calculate_specific_leadfield2(ldf::Leadfield, target::Union{Coordinate, Int}, supress::Union{Coordinate, Int}, region_size::Real; bilateral::Bool=false, kwargs...)

    if bilateral & isa(supress, Coordinate)
        supress.x = -1 * supress.x
    elseif bilateral & isa(supress, Int)
        supress = find_location(ldf, Talairach(-1 * ldf.x[supress], ldf.y[supress], ldf.z[supress]))
    end

    target_ldf = calculate_specific_leadfield2(ldf, target; kwargs...)
    supress_ldf = calculate_specific_leadfield_region(ldf, supress, region_size; kwargs...)

    hcat(target_ldf, supress_ldf )
end

function calculate_specific_leadfield2(ldf::Leadfield, target::Union{Coordinate, Int}, supress::Array{EEG.Talairach,1}, region_size::Real;
    bilateral::Bool=false, keep_vecs::Int=3, kwargs...)

    target_ldf = calculate_specific_leadfield2(ldf, target; kwargs...)

    for s in supress

        supress_ldf = calculate_specific_leadfield_region(ldf, s, region_size; kwargs...)

        supress_ldf = reduce_leadfield(supress_ldf, keep_vecs, keep_dims = 0; kwargs...)

        target_ldf = hcat(target_ldf, supress_ldf)
    end

    return target_ldf
end

"""
Return leadfield for location at given index
"""
function calculate_specific_leadfield2{A <: AbstractFloat}(ldf::Array{A, 3}, target::Int; kwargs...)

    squeeze(ldf[target,:,:], 1)'
end

"""
Return leadfield for given region given index
"""
function calculate_specific_leadfield_region(ldf::Leadfield, target::Int, region::Real; bilateral::Bool=false, kwargs...)

    supress_ldf = calculate_specific_leadfield_region(ldf, Talairach(ldf.x[target], ldf.y[target], ldf.z[target]), region; kwargs...)
end

"""
Return leadfield for given region given coordinaates
"""
function calculate_specific_leadfield_region(ldf::Leadfield, target::Coordinate, region::Real; bilateral::Bool=false, kwargs...)

    to_supress = falses(size(ldf.x))

    if bilateral
        target.x = -1 * target.x
    end

    for loc in 1:length(ldf.x)

        if euclidean(target, [ldf.x[loc], ldf.y[loc], ldf.z[loc]]) < region
            to_supress[loc] = true
        end
    end
    Ls = ldf.L[to_supress, :, :]
    Ls = reshape(permutedims(Ls, [3, 2, 1]), size(ldf.L, 3), size(Ls, 1) * 3)
end

"""
Reduce leadfield dimensionality
"""
function reduce_leadfield{F <: AbstractFloat}(ldf::Array{F, 2}, keep_vecs; keep_dims = 3, kwargs...)
    hcat(ldf[:, 1:keep_dims], svdfact(ldf[:, keep_dims+1:end]).U[:, 1:keep_vecs])

end


text:
######### logistic regression on 1000 obs x 10 var  ###########

using MCMC

# generate a random dataset
srand(1)
const n = 1000
const nbeta = 10 # number of covariates, including intercept

const X = [ones(n) randn((n, nbeta-1))]  # covariates

const beta0 = randn((nbeta,))
const Y = rand(n) .< ( 1 ./ (1. + exp(- X * beta0))) # logistic response

# define model
ex = quote
	vars ~ Normal(0, 1.0)  # Normal prior, std 1.0 for predictors
	prob = 1 / (1. + exp(- X * vars))
	Y ~ Bernoulli(prob)
end

m = model(ex, vars=zeros(nbeta), gradient=true)

# run random walk metropolis (10000 steps, 1000 for burnin)
mcchain01 = run(m * RWM(0.05) * SerialMC(1000:10000))

describe(mcchain01)

# run Hamiltonian Monte-Carlo (10000 steps, 1000 for burnin, 2 inner steps, 0.1 inner step size)
mcchain02 = run(m * HMC(2, 0.1) * SerialMC(1000:10000))

acceptance(mcchain02)

# run NUTS HMC (10000 steps, 1000 for burnin)
mcchain03 = run(m * NUTS() * SerialMC(1000:10000))

var(mcchain03)


text:
module StippleUI

import Genie
import Stipple

const DEFAULT_WRAPPER = Genie.Renderer.Html.template

#===#

function theme() :: String
  Genie.Router.route("/css/stipple/quasar.min.css") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "css", "quasar.min.css"), String),
      :css) |> Genie.Renderer.respond
  end

  Stipple.Elements.stylesheet("/css/stipple/quasar.min.css")
end

#===#

function deps() :: String
  Genie.Router.route("/js/stipple/quasar.umd.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "quasar.umd.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Renderer.Html.script(src="/js/stipple/quasar.umd.min.js")
end

#===#

include("API.jl")
include("Badge.jl")
include("Banner.jl")
include("BigNumber.jl")
include("Button.jl")
include("Card.jl")
include("Checkbox.jl")
include("Chip.jl")
include("Dashboard.jl")
include("Dialog.jl")
include("Form.jl")
include("FormInput.jl")
include("Heading.jl")
include("Icon.jl")
include("List.jl")
include("Radio.jl")
include("Range.jl")
include("Select.jl")
include("Separator.jl")
include("Space.jl")
include("Table.jl")
include("Toggle.jl")

function __init__()
  push!(Stipple.Layout.THEMES, theme)
  push!(Stipple.DEPS, deps)
end

end


text:
using ForwardDiff
using LinearAlgebra
include("backtracking_line_search.jl")

function newton_algo(epsilon,objective,x)
    g = x -> ForwardDiff.gradient(objective, x); # g = ∇f
    gradient = g(x)
    hessian = ForwardDiff.hessian(objective, x)
    while sum(abs,gradient) > epsilon
        gradient = g(x)
        hessian = ForwardDiff.hessian(objective, x)
        x = x - (gradient' * inv(hessian))'
    end
    return x
end

max_stars_repo_name: zhuhs-q/basic-optimal-method

max_issues_repo_name: Jacob-jinxy/basic-optimal-method

max_forks_repo_name: Jacob-jinxy/basic-optimal-method


text:
f(x, y) = x * sqrt(y)
theoric(t) = (t ^ 2 + 4.0) ^ 2 / 16.0

rk4(f) = (t, y, δt) ->  # 1st (result) lambda
         ((δy1) ->      # 2nd lambda
         ((δy2) ->      # 3rd lambda
         ((δy3) ->      # 4th lambda
         ((δy4) -> ( δy1 + 2δy2 + 2δy3 + δy4 ) / 6 # 5th and deepest lambda: calc y_{n+1}
         )(δt * f(t + δt, y + δy3))         # calc δy₄
         )(δt * f(t + δt / 2, y + δy2 / 2)) # calc δy₃
         )(δt * f(t + δt / 2, y + δy1 / 2)) # calc δy₂
         )(δt * f(t, y))                    # calc δy₁

δy = rk4(f)
t₀, δt, tmax = 0.0, 0.1, 10.0
y₀ = 1.0

t, y = t₀, y₀
while t ≤ tmax
    if t ≈ round(t) @printf("y(%4.1f) = %10.6f\terror: %12.6e\n", t, y, abs(y - theoric(t))) end
    y += δy(t, y, δt)
    t += δt
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function save(f::Stream{format"STL_ASCII"}, mesh::AbstractMesh)
    io      = stream(f)
    vts     = decompose(Point{3, Float32}, mesh)
    fcs     = decompose(Face{3, Cuint, -1}, mesh)
    normals = decompose(Normal{3, Float32}, mesh)

    nV = length(vts)
    nF = length(fcs)

    # write the header
    write(io,"solid vcg\n")
    # write the data
    for i = 1:nF
        f = fcs[i]
        n = normals[f][1] # TODO: properly compute normal(f)
        v1, v2, v3 = vts[f]
        @printf io "  facet normal %e %e %e\n" n[1] n[2] n[3]
        write(io,"    outer loop\n")

        @printf io "      vertex  %e %e %e\n" v1[1] v1[2] v1[3]
        @printf io "      vertex  %e %e %e\n" v2[1] v2[2] v2[3]
        @printf io "      vertex  %e %e %e\n" v3[1] v3[2] v3[3]

        write(io,"    endloop\n")
        write(io,"  endfacet\n")
    end
    write(io,"endsolid vcg\n")
end

@compat show(io::IO, ::MIME"model/stl+ascii", mesh::AbstractMesh) = save(io, mesh)

function save(f::Stream{format"STL_BINARY"}, mesh::AbstractMesh)
    io      = stream(f)
    vts     = decompose(Point{3, Float32}, mesh)
    fcs     = decompose(Face{3, Cuint, -1}, mesh)
    normals = decompose(Normal{3, Float32}, mesh)
    nF = length(fcs)
    # Implementation made according to https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
    for i in 1:80 # write empty header
        write(io,0x00)
    end

    write(io, UInt32(nF)) # write triangle count
    for i = 1:nF
        f = fcs[i]
        n = normals[f][1] # TODO: properly compute normal(f)
        v1, v2, v3 = vts[f]
        for j=1:3; write(io, n[j]); end # write normal

        for v in [v1, v2, v3]
            for j = 1:3
                write(io, v[j]) # write vertex coordinates
            end
        end
        write(io,0x0000) # write 16bit empty bit
    end
end

function load(fs::Stream{format"STL_BINARY"}, MeshType=GLNormalMesh)
    #Binary STL
    #https://en.wikipedia.org/wiki/STL_%28file_format%29#Binary_STL
    io = stream(fs)
    read(io, 80) # throw out header

    triangle_count = read(io, UInt32)
    FaceType    = facetype(MeshType)
    VertexType  = vertextype(MeshType)
    NormalType  = normaltype(MeshType)

    faces       = Array(FaceType,   triangle_count)
    vertices    = Array(VertexType, triangle_count*3)
    normals     = Array(NormalType, triangle_count*3)
    i = 0
    while !eof(io)
        faces[i+1]      = Face{3, Int, -1}(i*3, i*3+1, i*3+2)
        normals[i*3+1]  = NormalType(read(io, Float32), read(io, Float32), read(io, Float32))
        normals[i*3+2]  = normals[i*3+1] # hurts, but we need per vertex normals
        normals[i*3+3]  = normals[i*3+1]
        vertices[i*3+1] = VertexType(read(io, Float32), read(io, Float32), read(io, Float32))
        vertices[i*3+2] = VertexType(read(io, Float32), read(io, Float32), read(io, Float32))
        vertices[i*3+3] = VertexType(read(io, Float32), read(io, Float32), read(io, Float32))

        skip(io, 2) # throwout 16bit attribute
        i += 1
    end
    return MeshType(vertices=vertices, faces=faces, normals=normals)
end

function load(fs::Stream{format"STL_ASCII"}, MeshType=GLNormalMesh)
    #ASCII STL
    #https://en.wikipedia.org/wiki/STL_%28file_format%29#ASCII_STL
    io = stream(fs)

    FaceType   = facetype(MeshType)
    VertexType = vertextype(MeshType)
    vs         = VertexType[]
    fs         = FaceType[]

    topology   = false
    vert_count = 0
    vert_idx   = [0,0,0]

    while !eof(io)
        line = split(lowercase(readline(io)))
        if !isempty(line) && line[1] == "facet"
            #normal = NormalType(line[3:5])
            readline(io) # Throw away outerloop
            for i=1:3
                vertex = VertexType(split(readline(io))[2:4])
                if topology
                    idx = findfirst(vertices(mesh), vertex)
                end
                if topology && idx != 0
                    vert_idx[i] = idx
                else
                    push!(vs, vertex)
                    vert_count += 1
                    vert_idx[i] = vert_count
                end
            end
            readline(io) # throwout endloop
            readline(io) # throwout endfacet
            push!(fs, Face{3, Int, 0}(vert_idx...))
        end
    end
    return MeshType(vs, fs)
end