text:
typealias Arrayd{dim} Array{Array{Float64,dim},1}

# generate the data:  X
tmpx = [rand(25), randn(50)/10 + .8]
tmpy = tmpx + rand(size(tmpx)) .* 1.1
N = length(tmpx)
X = Array{Float64,1}[[tmpx[i]-mean(tmpx), tmpy[i]-mean(tmpy)] for i=1:N]
kappa     = Array{Float64,1}[X[i]       for i in 1:N ]
eta_coeff = deepcopy(kappa)
phix      = X
Dphix     = Array{Float64,2}[eye(2) for i in 1:N]
sigma = 0.1
epsilon = 0.02
dim = 2
#---------------------------------------
function r(d,sigma)
	 exp(-.5*(d/sigma).^2)
end
function rp(d,sigma)
	 -d ./ (sigma^2 * exp(d.^2./(2*sigma^2)))
end
function rp_div_d(d,sigma)
	 -1.0./(sigma^2*exp(d.^2/(2*sigma^2)))
end
function rpp(d,sigma)
	 d.^2./(sigma^4*exp(d.^2/(2*sigma^2))) - 1./(sigma^2*exp(d.^2/(2*sigma^2)))
end
"TODO: test everything below this line"
function rppp(d,sigma)
	 (3*d)./(sigma^4*exp(d.^2/(2*sigma^2))) - d.^3./(sigma^6*exp(d.^2/(2*sigma^2)))
end
function R{T<:Number}(x::Array{T,1},y::Array{T,1},sigma)
	r(norm(x-y),sigma)
end
function gradR{T<:Number}(x::Array{T,1},y::Array{T,1},sigma)
	v=x-y
	n=norm(v)
	rp_div_d(n,sigma).*v
end
function outer{T<:Number}(u::Array{T,1},v::Array{T,1})
	length(u) == 1 ? u[1]*v[1] : u*transpose(v)
end
function g1g2R{T<:Number}(x::Array{T,1},y::Array{T,1},sigma)
	v=x-y
	n=norm(v)
	eey = length(x) == 1 ? 1.0 : eye(length(x))
	G=rp_div_d(n,sigma)*eey
	if n!=0
		u=v/n
		G+=outer(u,-u) *(rpp(n,sigma)-rp_div_d(n,sigma))
	end
	G
end
function g1g1R{T<:Number}(x::Array{T,1},y::Array{T,1},sigma)
	 -1*g1g2R(x,y,sigma)
end
function prodc{dim}(pp::Float64, cellA::Arrayd{dim})
	Array{Float64,dim}[ pp*cellA[i] for i = 1:length(cellA) ]
end
function target(kappa::Array{Array{Float64,1},1})
	targsig = 0.5
	density = Float64[]
	gradH   = Array{Float64,1}[]
	for k = 1:length(kappa)
		push!(density, exp(-(kappa[k][1].^2 + kappa[k][2].^2) / (2 * targsig^2)  ) / (2 * pi * targsig^2) )
		push!(gradH,[-kappa[k][1]/(targsig^2) , -kappa[k][2] / (targsig^2)])
	end
	density, gradH
end

#### test these ##########
function transd_dlDphix_dt{dim}(eta_coeff,  kappa,  phix, dlDphix::Arrayd{dim}, sigma)
	n_knots = length(eta_coeff)
	n_phis  = length(dlDphix)
	returnval = dim == 1 ? Array{Float64,1}[[0.0] for i in 1:n_phis] : Array{Float64,dim}[zeros(dim,dim) for i in 1:n_phis]
	for col = 1:dim, i = 1:n_phis, j = 1:n_knots
		returnval[i][:,col] += gradR(phix[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlDphix[i][:,col]
	end
	prodc(-1.0, returnval)
end
function transd_dlphix_dt{dim}(eta_coeff, kappa, phix, dlphix, Dphix, dlDphix::Arrayd{dim}, sigma)
	n_knots  = length(eta_coeff)
	n_phis  = length(dlphix)
	returnval = Array{Float64,1}[zeros(dim) for i in 1:n_phis]
	for i = 1:n_phis, j = 1:n_knots
		returnval[i] +=  gradR(phix[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlphix[i]
		for col = 1:dim
			returnval[i] +=  g1g1R(phix[i], kappa[j],sigma) * Dphix[i][:,col] * transpose(eta_coeff[j]) * dlDphix[i][:,col]
		end
	end
	prodc(-1.0, returnval)
end
function transd_dlkappa_dt{dim}(eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix::Arrayd{dim}, sigma)
	n_knots  = length(eta_coeff)
	n_phis   = length(dlphix)
	returnval = Array{Float64,1}[zeros(dim) for i in 1:n_knots]
	for i = 1:n_knots, j = 1:n_knots
		returnval[i] +=  -dot(eta_coeff[i],eta_coeff[j]) .* (g1g1R(kappa[i], kappa[j], sigma) * dleta_coeff[i])
		returnval[i] +=  -dot(eta_coeff[i],eta_coeff[j]) .* (g1g2R(kappa[j], kappa[i], sigma) * dleta_coeff[j])
		returnval[i] +=  -1.0 * gradR(kappa[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlkappa[i]
		returnval[i] +=         gradR(kappa[j], kappa[i],sigma) * transpose(eta_coeff[i]) * dlkappa[j]
	end
	for i = 1:n_knots, j = 1:n_phis
		returnval[i] +=  -1.0*gradR(phix[j], kappa[i],sigma) * transpose(eta_coeff[i]) * dlphix[j]
		for col = 1:dim
			returnval[i] +=  g1g2R(phix[j], kappa[i],sigma) * Dphix[j][:,col] * transpose(eta_coeff[i]) * dlDphix[j][:,col]
		end
	end
	prodc(-1.0, returnval)
end
function transd_dleta_dt{dim}(eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix::Arrayd{dim}, sigma)
	n_knots = length(kappa)
	n_phis  = length(phix)
	returnval = Array{Float64,1}[ zeros(dim) for i in 1:n_knots]
	for i = 1:n_knots, j = 1:n_knots
		returnval[i] +=  -eta_coeff[j] * transpose(gradR(kappa[i], kappa[j],sigma)) * dleta_coeff[i]
		returnval[i] +=  -eta_coeff[j] * transpose(gradR(kappa[j], kappa[i],sigma)) * dleta_coeff[j]
		returnval[i] +=  R(kappa[j], kappa[i], sigma) * dlkappa[j]
	end
	for i = 1:n_knots, j = 1:n_phis
		returnval[i] +=  R(phix[j], kappa[i], sigma) * dlphix[j]
		for col = 1:dim
			returnval[i] += dot(Dphix[j][:,col], gradR(phix[j], kappa[i],sigma)) * dlDphix[j][:,col]
		end
	end
	prodc(-1.0, returnval)
end
#### against these ...test these ##########
# transDphix   = Array{Float64,2}[zeros(dim,dim) for i in 1:N]
function transd_dlDphix_dt!{dim}(transDphix::Arrayd{dim}, eta_coeff::Arrayd{1}, kappa::Arrayd{1}, phix::Arrayd{1}, dlDphix::Arrayd{dim}, sigma, epsilon)
	for i = 1:length(transDphix)
		transDphix[i] = zero(transDphix[1])
		for j = 1:length(eta_coeff)
			for col = 1:dim
				transDphix[i][:,col] =  transDphix[i][:,col] - epsilon * gradR(phix[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlDphix[i][:,col]
			end
		end
	end
end
# transdlphix   = Array{Float64,1}[zeros(dim) for i in 1:N]
function transd_dlphix_dt!{oneOrTwo}(transdlphix::Arrayd{1}, eta_coeff::Arrayd{1}, kappa::Arrayd{1}, phix::Arrayd{1}, dlphix::Arrayd{1}, Dphix::Arrayd{oneOrTwo}, dlDphix::Arrayd{oneOrTwo}, sigma, epsilon)
	dim = length(kappa[1])
	for i = 1:length(transdlphix)
		transdlphix[i] = zero(transdlphix[1])
		for j = 1:length(kappa)
			transdlphix[i] += epsilon * gradR(phix[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlphix[i]
			for col = 1:dim
				transdlphix[i] += epsilon * g1g1R(phix[i], kappa[j],sigma) * Dphix[i][:,col] * transpose(eta_coeff[j]) * dlDphix[i][:,col]
			end
		end
	end
end
# transdlKappa   = Array{Float64,1}[zeros(dim) for i in 1:N]
function transd_dlkappa_dt!{dim}(transdlKappa::Arrayd{1}, eta_coeff::Arrayd{1}, dleta_coeff::Arrayd{1}, kappa::Arrayd{1}, dlkappa::Arrayd{1}, phix::Arrayd{1}, dlphix::Arrayd{1}, Dphix::Arrayd{dim}, dlDphix::Arrayd{dim}, sigma, epsilon)
	for i = 1:length(transdlKappa)
		transdlKappa[1] = zero(transdlKappa[1])
		for j = 1:length(kappa)
			transdlKappa[i] += epsilon * dot(eta_coeff[i],eta_coeff[j]) .* (g1g1R(kappa[i], kappa[j], sigma) * dleta_coeff[i])
			transdlKappa[i] += epsilon * dot(eta_coeff[i],eta_coeff[j]) .* (g1g2R(kappa[j], kappa[i], sigma) * dleta_coeff[j])
			transdlKappa[i] += epsilon * gradR(kappa[i], kappa[j],sigma) * transpose(eta_coeff[j]) * dlkappa[i]
			transdlKappa[i] -= epsilon * gradR(kappa[j], kappa[i],sigma) * transpose(eta_coeff[i]) * dlkappa[j]
		end
	end
	for i = 1:length(transdlKappa)
		for j = 1:length(phix)
			transdlKappa[i] += epsilon * gradR(phix[j], kappa[i],sigma) * transpose(eta_coeff[i]) * dlphix[j]
			for col = 1:dim
				transdlKappa[i] -=  epsilon * g1g2R(phix[j], kappa[i],sigma) * Dphix[j][:,col] * transpose(eta_coeff[i]) * dlDphix[j][:,col]
			end
		end
	end
end
# transdlEta   = Array{Float64,1}[zeros(dim) for i in 1:N]
function transd_dleta_dt!{dim}(transdlEta::Arrayd{1}, eta_coeff::Arrayd{1}, dleta_coeff::Arrayd{1}, kappa::Arrayd{1}, dlkappa::Arrayd{1}, phix::Arrayd{1}, dlphix::Arrayd{1}, Dphix::Arrayd{dim}, dlDphix::Arrayd{dim}, sigma, epsilon)
	for i = 1:length(kappa)
		transdlEta[i] = zero(transdlEta[i])
		for j = 1:length(kappa)
			transdlEta[i] +=  epsilon * eta_coeff[j] * transpose(gradR(kappa[i], kappa[j],sigma)) * dleta_coeff[i]
			transdlEta[i] +=  epsilon * eta_coeff[j] * transpose(gradR(kappa[j], kappa[i],sigma)) * dleta_coeff[j]
			transdlEta[i] -=  epsilon * R(kappa[j], kappa[i], sigma) * dlkappa[j]
		end
	end
	for i = 1:length(kappa), j = 1:length(phix)
		transdlEta[i] -= epsilon * R(phix[j], kappa[i], sigma) * dlphix[j]
		for col = 1:dim
			transdlEta[i] -= epsilon * dot(Dphix[j][:,col], gradR(phix[j], kappa[i],sigma)) * dlDphix[j][:,col]
		end
	end
end

dlphix      = Array{Float64,1}[target(phix)[2][i]/length(phix) for i = 1:length(phix)]
dleta_coeff = Array{Float64,1}[zeros(dim) for i = 1:N]
dlkappa     = Array{Float64,1}[zeros(dim) for i = 1:N]
transdlphix  = Array{Float64,1}[zeros(dim) for i in 1:N]
transdlKappa = Array{Float64,1}[zeros(dim) for i in 1:N]
transdlEta   = Array{Float64,1}[zeros(dim) for i in 1:N]
transDphix   = Array{Float64,2}[zeros(dim,dim) for i in 1:N]
dlDphix = Array{Float64,2}[inv(Dphix[i])/length(phix) for i = 1:length(phix)]

@time transDphixtest = transd_dlDphix_dt(eta_coeff,  kappa,  phix, dlDphix, sigma)
@time transd_dlDphix_dt!(transDphix, eta_coeff, kappa, phix, dlDphix, sigma, epsilon)

@time transdlphixtest = transd_dlphix_dt(eta_coeff, kappa, phix, dlphix, Dphix, dlDphix, sigma)
@time transd_dlphix_dt!(transdlphix, eta_coeff, kappa, phix, dlphix, Dphix, dlDphix, sigma, epsilon)

@time transdlKappatest = transd_dlkappa_dt(eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix, sigma)
@time transd_dlkappa_dt!(transdlKappa, eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix, sigma, epsilon)

@time transdlEtatest = transd_dleta_dt(eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix, sigma)
@time transd_dleta_dt!(transdlEta, eta_coeff, dleta_coeff, kappa, dlkappa, phix, dlphix, Dphix, dlDphix, sigma, epsilon)

error("exiting the script now")

[transDphix prodc(epsilon,transDphixtest)] # same
[transdlphix prodc(epsilon, transdlphixtest)] # different
[transdlKappa prodc(epsilon, transdlKappatest)] #same
[transdlEta prodc(epsilon, transdlEtatest)]

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

# speed check for Array{Float64,2} vrs Array{Array{Float64,1},1}
using  PyCall
@pyimport matplotlib.pyplot as plt
include("src/flow_ode.jl")
include("src/rfuncs.jl")
include("src/targets.jl")

# typealias Arrayd{dim} Array{Array{Float64,dim},1}
tmpx = [rand(50), randn(75)/10.0 + .8]
tmpy = tmpx + rand(size(tmpx)) .* 1.1
N = length(tmpx)
X = Array{Float64,1}[[tmpx[i]-mean(tmpx), tmpy[i]-mean(tmpy)] for i=1:N]
kappa_cell     = Array{Float64,1}[X[i]       for i in 1:N ]
eta_coeff_cell = Array{Float64,1}[zero(kappa_cell[i]) for i in 1:N ]

kappa_splat     = Array(Float64,(2,N))
eta_coeff_splat     = Array(Float64,(2,N))
for i = 1:N
	kappa_splat[:,i] = X[i]
	eta_coeff_splat[:,i] = [0.0, 0.0]
end
function d_eta_dt1(eta_coeff::Array{Float64,2}, kappa::Array{Float64,2})
	sigma = 1.0
	dd = size(kappa,1)
	tmp = zeros(kappa)
	for i = 1:size(kappa,2)
	 	for j = 1:size(kappa,2)
	 		dtij = dot(eta_coeff[:,i],eta_coeff[:,j])
	 		gRij = gradR(kappa[:,i], kappa[:,j],sigma)
	 		for z = 1:size(kappa,1)
				tmp[z,i] -=  dtij*gRij[z]
			end
		end
	end
	tmp
end
function d_eta_dt2(eta_coeff::Array{Float64,2}, kappa::Array{Float64,2})
	sigma = 1.0
	dd = size(kappa,1)
	tmp = zeros(kappa)
	for i = 1:size(kappa,2)
	 	for j = 1:size(kappa,2)
	 		dtij = 0.0
	 		for z = 1:size(kappa,1)
	 			dtij += eta_coeff[z,i]*eta_coeff[z,j]
	 		end
	 		gRij = gradR([1.0, 2.0], [1.0, 2.0],sigma)
	 		for z = 1:size(kappa,1)
				tmp[z,i] -=  dtij*gRij[z]
			end
		end
	end
	tmp
end
function d_eta_dt(eta_coeff, kappa)
	sigma = 1.0
	dd = size(eta_coeff[1])
	tmp = Array{Float64,1}[zeros(dd) for i in 1:length(eta_coeff)]
	for i = 1:length(kappa)
	 	for j = 1:length(kappa)
			tmp[i] -= dot(eta_coeff[i],eta_coeff[j]) * gradR(kappa[i], kappa[j],sigma)
		end
	end
	tmp
end

d_eta_dt1(eta_coeff_splat, kappa_splat)
d_eta_dt2(eta_coeff_splat, kappa_splat)
d_eta_dt(eta_coeff_cell, kappa_cell)

@time d_eta_dt1(eta_coeff_splat, kappa_splat);
@time d_eta_dt2(eta_coeff_splat, kappa_splat);
@time d_eta_dt(eta_coeff_cell, kappa_cell);

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

# speed check for Array{Float64,2} vrs Array{Array{Float64,1},1}
using  PyCall
@pyimport matplotlib.pyplot as plt
include("src/flow_ode.jl")
include("src/rfuncs.jl")
include("src/targets.jl")
# typealias
# typealias Arrayd{dim} Array{Array{Float64,dim},1}
tmpx = [rand(50), randn(75)/10.0 + .8]
tmpy = tmpx + rand(size(tmpx)) .* 1.1
N = length(tmpx)
X = Array{Float64,1}[[tmpx[i]-mean(tmpx), tmpy[i]-mean(tmpy)] for i=1:N]
kappa_cell     = Array{Float64,1}[X[i]       for i in 1:N ]
eta_coeff_cell = Array{Float64,1}[zero(kappa_cell[i]) for i in 1:N ]

function d_eta_dt(eta_coeff::Arrayd{1}, kappa::Arrayd{1}, sigma, epsilon)
	dEtaDt   = Array{Float64,1}[zeros(2) for i in 1:length(eta_coeff)]
	for i = 1:length(dEtaDt)
		dEtaDt[i] = zero(kappa[1]) # you've got to start it at zero
		for j = 1:length(kappa)
			dEtaDt[i] -=  epsilon * dot(eta_coeff[i],eta_coeff[j]) * gradR(kappa[i], kappa[j], sigma)
		end
	end
	dEtaDt
end
function d_kappa_dt(eta_coeff::Arrayd{1}, kappa::Arrayd{1}, sigma, epsilon)
	dKappaDt = Array{Float64,1}[zeros(2) for i in 1:length(kappa)]
	for i = 1:length(dKappaDt)
		dKappaDt[i] = zero(kappa[1]) # you've got to start it at zero
		for j = 1:length(kappa)
			dKappaDt[i] += epsilon * eta_coeff[j] * R(kappa[i], kappa[j], sigma)
		end
	end
	dKappaDt
end
function d_phix_dt(eta_coeff::Arrayd{1}, kappa::Arrayd{1}, phix::Arrayd{1}, sigma, epsilon)
	dPhixDt  = Array{Float64,1}[zeros(2) for i in 1:length(phix)]
	for i = 1:length(dPhixDt)
		dPhixDt[i] = zero(phix[1])
		for j = 1:length(kappa)
			dPhixDt[i] += epsilon * eta_coeff[j]* R(phix[i], kappa[j], sigma)
		end
	end
	dPhixDt
end
function d_Dphix_dt(eta_coeff::Arrayd{1}, kappa::Arrayd{1}, phix::Arrayd{1}, Dphix::Arrayd{2}, sigma, epsilon)
    dDphixDt = Array{Float64,2}[zeros(2,2) for i in 1:length(phix)]
	for i = 1:length(dDphixDt)
		dDphixDt[i] = zero(Dphix[1]) # you've got to start it at zero
		for j = 1:length(kappa)
			dDphixDt[i] += epsilon * eta_coeff[j] * transpose(gradR(phix[i], kappa[j], sigma)) * Dphix[i]
		end
	end
	dDphixDt
end

sigma = 0.1
epsilon = 0.1
D =  Array{Float64,2}[ones(2,2) for i in 1:length(X)]
@time d_eta_dt(eta_coeff_cell, kappa_cell, sigma, epsilon);
@time d_kappa_dt(eta_coeff_cell, kappa_cell, sigma, epsilon);
@time d_phix_dt(eta_coeff_cell, kappa_cell, X, sigma, epsilon);
@time d_Dphix_dt(eta_coeff_cell, kappa_cell, X, D, sigma, epsilon);
@time td_Dphix_dt(eta_coeff_cell, kappa_cell, X, D, sigma, epsilon);

# from the  above code it is clear that the bottle neck is gradR not the looping

type Flow
	kappa::Array{Float64,2}
	eta_coeff::Array{Float64,2}
end

flowinstance1 = Flow(rand(2,2), rand(3,3))
flowinstance1.kappa # prints out kappa
flowinstance1.eta_coeff # prints out eta

# if you want a default for one of the parameters you can make these constructors
Flow(x)= Flow(x, eye(2,2))
Flow()= Flow(eye(2,2), eye(2,2))

flowinstance2 = Flow(rand(2,2))
flowinstance3 = Flow()

# Define + on Flow types
function +(a::Flow, b::Flow)
	Flow(a.kappa + b.kappa, a.eta_coeff + b.eta_coeff)
end

flowinstance2 + flowinstance3

# to try and make things faster you can define an immutable type which tells the complier that
# the fields will not be re-bound.

immutable ImFlow
	kappa::Array{Float64,2}
	eta_coeff::Array{Float64,2}
end

imflowinstance1 = ImFlow(rand(2,2), rand(3,3))

imflowinstance1.kappa = eye(2,2) # this gives an error since I'm trying to re-bind kappa
imflowinstance1.kappa[1] = 2.0 # this is ok since I'm just changing the entry kappa[1] and not the container kappa

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

#---------------------------------------
# kernel evals and derivatives...these are local to this module
#------------------------------------
r(d::Real,sigma) = exp(-0.5 * d * d / (sigma * sigma))
function rp_div_d(d::Real,sigma)
	s2 = sigma * sigma
	-exp(-0.5 * d * d / s2) / s2
end
function rpp(d::Real,sigma)
	rd = r(d,sigma)
	s2 = sigma * sigma
	d * d * rd / (s2 * s2) - rd / s2
end
function R{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
	r(norm(x-y),sigma)
end
function gradR{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
	v=x-y
	n=norm(v)
	v * rp_div_d(n,sigma)
end
function outer{T<:Real}(u::Array{T,1},v::Array{T,1})
	length(u) == 1 ? u[1]*v[1] : u*transpose(v)
end
function g1g2R{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
	v = x-y
	n = norm(v)
	u = v/n
	eey = length(x) == 1 ? 1.0 : eye(length(x))
	rpd = rp_div_d(n,sigma)
	if n == 0
		G = -rpp(n,sigma) * eey
		return G
	else
		G = -rpd * eey
		G += outer(u,-u) * (rpp(n,sigma) - rpd)
		return G
	end
end
function g1g1R{T<:Real}(x::Array{T,1},y::Array{T,1},sigma)
	 -g1g2R(x,y,sigma)
end

x=[1.0, 2.0]
y=[1.0, 2.2]
@time g1g2R(x, y, 2.0)
@time g1g2R2(x, y, 2.0)

@time g1g2R(x, y, 2.0)
@time g1g2R2(x, y, 2.0)

#---------------------
# I want to test the impact of  the type conversion
# ------------------------

include("src/rfuncs.jl")
typealias Array1 Array{Array{Float64,1},1} # for lists of vectors
typealias Array2 Array{Array{Float64,2},1} # for lists of jacobians
array1(dim, k) = Array{Float64,1}[zeros(dim) for i=1:k] # constructor for Array1
array2(dim, k) = Array{Float64,2}[zeros(dim, dim) for i=1:k] # constructor for Array2
array2eye(dim, k) = Array{Float64,2}[eye(dim) for i=1:k] # constructor for Array2
immutable Flow
	kappa::Array1
	eta_coeff::Array1
	phix::Array1
	Dphix::Array2
	dim::Int64
end
Flow(dim::Int64, k) = Flow(array1(dim, k), array1(dim, k), array1(dim, k), array2(dim, k), dim) # zero Flow constructor
for opt = (:+, :-, :*) # this defines +, - and * between Flow and Float64
	@eval begin
		function ($opt)(yin::Flow, a::Real)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = Flow(yin.dim, 0)
			for i = 1:nKap
				push!(yout.kappa, ($opt)(a,yin.kappa[i]))
				push!(yout.eta_coeff, ($opt)(a,yin.eta_coeff[i]))
			end
			for i = 1:nPhi
				push!(yout.phix, ($opt)(a,yin.phix[i]))
				push!(yout.Dphix, ($opt)(a,yin.Dphix[i]))
			end
			yout
		end # end function definition
	end # end the begin block
	@eval ($opt)(a::Real, yin::Flow) = ($opt)(yin::Flow, a::Float64)
end # end for loop over functions
# ytest1 = Flow(2, 3)
# ytest2 = ytest1 + 5.7
# ytest3 = ytest2 * 1.7 + 4
function d_forward_dt(yin::Flow, sigma)
	nKap = length(yin.kappa)
	nPhi = length(yin.phix)
	yout = Flow(yin.dim, nKap)
	# eta
	for i = 1:nKap, j = 1:nKap
		yout.eta_coeff[i] -= dot(yin.eta_coeff[i], yin.eta_coeff[j]) * gradR(yin.kappa[i], yin.kappa[j], sigma)
	end
	# kappa
	for i = 1:nKap, j = 1:nKap
		yout.kappa[i] += yin.eta_coeff[j] * R(yin.kappa[i], yin.kappa[j], sigma)
	end
	# phi
	for i = 1:nPhi, j = 1:nKap
		yout.phix[i] += yin.eta_coeff[j] * R(yin.phix[i], yin.kappa[j], sigma)
	end
	# Dphi
	for i = 1:nPhi, j = 1:nKap
		yout.Dphix[i] += yin.eta_coeff[j] * transpose(gradR(yin.phix[i], yin.kappa[j], sigma)) * yin.Dphix[i]
	end
	yout
end

#transpose
immutable TrFlow
	kappa::Array1
	eta_coeff::Array1
	phix::Array1
	Dphix::Array2
	dlkappa::Array1
	dleta_coeff::Array1
	dlphix::Array1
	dlDphix::Array2
	dim::Int64
end

TrFlow(dim::Int64, k) = TrFlow(array1(dim, k), array1(dim, k), array1(dim, k), array2(dim, k), array1(dim, k), array1(dim, k), array1(dim, k), array2(dim, k), dim) # zero Flow constructor

for opt = (:+, :-, :*) # this defines +, - and * between Flow and Float64
	@eval begin
		function ($opt)(yin::TrFlow, a::Real)
			nKap = length(yin.kappa)
			nPhi = length(yin.phix)
			yout = TrFlow(yin.dim, 0)
			for i = 1:nKap
				push!(yout.kappa, ($opt)(a,yin.kappa[i]))
				push!(yout.eta_coeff, ($opt)(a,yin.eta_coeff[i]))
				push!(yout.dlkappa, ($opt)(a,yin.dlkappa[i]))
				push!(yout.dleta_coeff, ($opt)(a,yin.dleta_coeff[i]))
			end
			for i = 1:nPhi
				push!(yout.phix, ($opt)(a,yin.phix[i]))
				push!(yout.Dphix, ($opt)(a,yin.Dphix[i]))
				push!(yout.dlphix, ($opt)(a,yin.dlphix[i]))
				push!(yout.dlDphix, ($opt)(a,yin.dlDphix[i]))
			end
			yout
		end # end function definition
	end # end the begin block
	@eval ($opt)(a::Real, yin::Flow) = ($opt)(yin::Flow, a::Float64)
end # end for loop over functions

function d_backward_dt(yin::TrFlow, sigma)
	nKap = length(yin.kappa)
	nPhi = length(yin.phix)
	yout = TrFlow(yin.dim, nKap)
	# eta
	for i = 1:nKap, j = 1:nKap
		yout.eta_coeff[i] += dot(yin.eta_coeff[i], yin.eta_coeff[j]) * gradR(yin.kappa[i], yin.kappa[j], sigma)
	end
	# kappa
	for i = 1:nKap, j = 1:nKap
		yout.kappa[i] -= yin.eta_coeff[j] * R(yin.kappa[i], yin.kappa[j], sigma)
	end
	# phi
	for i = 1:nPhi, j = 1:nKap
		yout.phix[i] -= yin.eta_coeff[j] * R(yin.phix[i], yin.kappa[j], sigma)
	end
	# Dphi
	for i = 1:nPhi, j = 1:nKap
		yout.Dphix[i] -= yin.eta_coeff[j] * transpose(gradR(yin.phix[i], yin.kappa[j], sigma)) * yin.Dphix[i]
	end
	# dlDphix
	for col = 1:yin.dim, i = 1:nPhi, j = 1:nKap
		yout.dlDphix[i][:,col] +=  gradR(yin.phix[i], yin.kappa[j], sigma) * transpose(yin.eta_coeff[j]) * yin.dlDphix[i][:,col]
	end
	# dlphix
	for i = 1:nPhi, j = 1:nKap
		yout.dlphix[i] += gradR(yin.phix[i], yin.kappa[j],sigma) * transpose(yin.eta_coeff[j]) * yin.dlphix[i]
		for col = 1:yin.dim
			yout.dlphix[i] +=  g1g1R(yin.phix[i], yin.kappa[j], sigma) * yin.Dphix[i][:,col] * transpose(yin.eta_coeff[j]) * yin.dlDphix[i][:,col]
		end
	end
	# dlkappa
	for i = 1:nKap, j = 1:nKap
		yout.dlkappa[i] -= dot(yin.eta_coeff[i], yin.eta_coeff[j]) .* (g1g1R(yin.kappa[i], yin.kappa[j], sigma) * yin.dleta_coeff[i])
		yout.dlkappa[i] -= dot(yin.eta_coeff[i], yin.eta_coeff[j]) .* (g1g2R(yin.kappa[j], yin.kappa[i], sigma) * yin.dleta_coeff[j])
		yout.dlkappa[i] += gradR(yin.kappa[i], yin.kappa[j],sigma) * transpose(yin.eta_coeff[j]) * yin.dlkappa[i]
		yout.dlkappa[i] -= gradR(yin.kappa[j], yin.kappa[i],sigma) * transpose(yin.eta_coeff[i]) * yin.dlkappa[j]
	end
	for i = 1:nKap, j = 1:nPhi
		yout.dlkappa[i] -= gradR(yin.phix[j], yin.kappa[i], sigma) * transpose(yin.eta_coeff[i]) * yin.dlphix[j]
		for col = 1:yin.dim
			yout.dlkappa[i] += g1g2R(yin.phix[j], yin.kappa[i],sigma) * yin.Dphix[j][:,col] * transpose(yin.eta_coeff[i]) * yin.dlDphix[j][:,col]
		end
	end
	# dleta_coeff
	for i = 1:nKap, j = 1:nKap
		yout.dleta_coeff[i] -= yin.eta_coeff[j] * transpose(gradR(yin.kappa[i], yin.kappa[j], sigma)) * yin.dleta_coeff[i]
		yout.dleta_coeff[i] -= yin.eta_coeff[j] * transpose(gradR(yin.kappa[j], yin.kappa[i], sigma)) * yin.dleta_coeff[j]
		yout.dleta_coeff[i] += R(yin.kappa[j], yin.kappa[i], sigma) * yin.dlkappa[j]
	end
	for i = 1:nKap, j = 1:nPhi
		yout.dleta_coeff[i] += R(yin.phix[j], yin.kappa[i], sigma) * yin.dlphix[j]
		for col = 1:yin.dim
			yout.dleta_coeff[i] += dot(yin.Dphix[j][:,col], gradR(yin.phix[j], yin.kappa[i],sigma)) * yin.dlDphix[j][:,col]
		end
	end
	yout
end

tmpx = [rand(50), randn(75)/10.0 + .8]; tmpy = tmpx + rand(size(tmpx)) .* 1.1; nKap = length(tmpx)
X = Array{Float64,1}[[tmpx[i]-mean(tmpx), tmpy[i]-mean(tmpy)] for i=1:nKap]

dim = length(X[1])
yFlow = Flow(deepcopy(X), array1(dim,nKap), deepcopy(X), array2(dim,nKap), dim)
yTrFlow = TrFlow(deepcopy(X), array1(dim,nKap), deepcopy(X), array2(dim,nKap), deepcopy(X), array1(dim,nKap), deepcopy(X), array2(dim,nKap), dim)

@time d_forward_dt(yFlow, 1.0);
@time d_backward_dt(yTrFlow, 1.0);

max_stars_repo_name: EthanAnderes/MleDensityWarp

max_issues_repo_name: EthanAnderes/MleDensityWarp

max_forks_repo_name: EthanAnderes/MleDensityWarp


text:
using JLD
include("Transformation.jl")
include("AbstractSystem.jl")
include("Tree.jl")
include("Evaluation.jl")
include("BackPropogation.jl")
include("Facility.jl")
import Base.push!

type Axiom
    tree1 :: Tree
    tree2 :: Tree
    index :: Dict
end

function Base.show(io :: IO, m :: Axiom)
    print(io, "Axiom[")
    print(io, m.tree1)
    print(io, ", ")
    print(io, m.tree2)
    print(io, "]")
end

function beautify(m :: Axiom)
    string("Axiom:\n", beautify(m.tree1, 1), "\n", beautify(m.tree2, 1), "\n")
end

Axiom(tree1 :: Tree, tree2 :: Tree) = Axiom(tree1, tree2, add!(index(tree1), index(tree2)))

Axiom(skeleton1, skeleton2) = Axiom(toTree(skeleton1), toTree(skeleton2))

function push!(index :: Dict, ops :: Array)
    n = length(ops)
    for i in 1:n
        for t in index[i]
            t.op = ops[i]
        end
    end
end

function push!(index :: Dict, ops :: Dict)
    for key in keys(ops)
        if haskey(index, key)
            ts = index[key]
            for t in ts
                t.op = ops[key]
            end
        end
    end
end

function push!(axiom :: Axiom, ops)
    push!(axiom.index, variables)
end

function init_axiom!(axiom :: Axiom, variables)
    push!(axiom.index, variables)
    init_tree!(axiom.tree1)
    init_tree!(axiom.tree2)
end

init_axioms! = distribute(init_axiom!)

## add some loss function to deal with degenerating problem?

function loss(a, b)
    b - a
end

function train_axiom!(axiom :: Axiom, variables, n = 1, randomize = identity)
    push!(axiom.index, variables)
    d1 :: Array{Float64, 1} = axiom.tree1.value[:d]
    d2 :: Array{Float64, 1} = axiom.tree2.value[:d]
    v1 :: Array{Float64, 1} = axiom.tree1.value[:value]
    v2 :: Array{Float64, 1} = axiom.tree2.value[:value]
    for i in 1:n
        randomize(variables)
        push!(axiom.index, variables)
        eval_tree!(axiom.tree2)
        eval_tree!(axiom.tree1)
        for j in 1:length(d1)
            d1[j] = loss(v1[j], v2[j])
            d2[j] = loss(v2[j], v1[j])
        end
        bp_tree!(axiom.tree1)
        bp_tree!(axiom.tree2)
    end
end

train_axioms! = distribute(train_axiom!)

## To prevent degeneration problem, we use anti-training to deal with the problem.
## The current method to deal with degeneration is anti-training: to push random evaluation results far away,
## but the defect now is that it is unstable in long term.
## One possible remedy is to anti-train only when the random evaluation results are "too near".

const DefaultThreshold = 0.1

function anti_train_axiom!(axiom :: Axiom, n = 1, randomize = identity, thres = DefaultThreshold)
    d1 :: Array{Float64, 1} = axiom.tree1.value[:d]
    d2 :: Array{Float64, 1} = axiom.tree2.value[:d]
    v1 :: Array{Float64, 1} = axiom.tree1.value[:value]
    v2 :: Array{Float64, 1} = axiom.tree2.value[:value]
    for i in 1:n
        randomize(axiom.tree1)
        randomize(axiom.tree2)
        eval_tree!(axiom.tree2)
        eval_tree!(axiom.tree1)
        if mean((v1 - v2) .^ 2) < thres
            for j in 1:length(d1)
                d1[j] = - loss(v1[j], v2[j])
                d2[j] = - loss(v2[j], v1[j])
            end
            bp_tree!(axiom.tree1)
            bp_tree!(axiom.tree2)
        end
    end
end

anti_train_axioms! = distribute(anti_train_axiom!)

include("DataStructs.jl")

UChar = Class("Char", 8)
Variable = Class("Variable", 100)
Stream = List(UChar, Variable, [:v, :u, :v1])

Prog = Class("Prog", 100)
Bindingc = Class("Binding", 100)
Binding = Pair(Variable, Prog, Bindingc, [:b, :v, :p])

Framec = Class("Frame", 100)
Frame = List(Bindingc, Framec, [:f1, :b, :f2])

Envc = Class("Env", 200)
Env = List(Framec, Envc, [:e1, :f, :e2])

lookup = DFunction("lookup", [Envc, Variable], Prog)
set = DFunction("set", [Envc, Variable, Prog], Envc)
add_binding = DFunction("add_binding", [Framec, Variable, Prog], Framec)
def = DFunction("def", [Envc, Variable, Prog], Envc)
extend = DFunction("extend", [Envc, Variable, Prog], Envc)

var = DFunction("var", [Variable], Prog)
definition = DFunction("definition", [Variable, Prog], Prog)
assignment = DFunction("assignment", [Variable, Prog], Prog)
procedure = DFunction("procedure", [Variable, Prog], Prog)
func_call = DFunction("func_call", [Prog, Prog], Prog)
Seq = List(Prog, Prog)

f_eval = DFunction("f_eval", [Prog, Envc], Prog)
s_eval = DFunction("eval", [Prog, Envc], Envc)

@load "infinity.jld"

Interpreter = Dict()
I_names = Dict()
I_names[:base] = [:UChar, :Variable, :Stream, :Prog, :Bindingc, :Binding, :Framec, :Frame, :Envc, :Env]
I_names[:env] = [:lookup, :set, :add_binding, :def, :extend]
I_names[:prog] = [:var, :definition, :assignment, :procedure, :func_call, :Seq]
I_names[:repl] = [:f_eval, :s_eval]
inject(I_names, Interpreter)
Interpreter

axiom_lookup1 = Axiom([lookup, [Env.cons, [Frame.cons, [Binding.pair, :v, :p], :f], :e], :v], :p)
axiom_lookup2 = Axiom([lookup, [Env.cons, [Frame.cons, [Binding.pair, :v1, :p1], :f2], :e2], :v2],
                      [lookup, [Env.cons, :f2, :e2], :v2])
axiom_lookup3 = Axiom([lookup, [Env.cons, Frame.empty, :e], :v], [lookup, :e, :v])

axiom_set1 = Axiom([set, [Env.cons, [Frame.cons, [Binding.pair, :v1, :p1], :f2], :e2], :v1, :p2],
                   [Env.cons, [Frame.cons, [Binding.pair, :v1, :p2], :f2], :e2])
axiom_set2 = Axiom([set, [Env.cons, [Frame.cons, [Binding.pair, :v1, :p1], :f2], :e2], :v2, :p2],
                   [Env.cons, [Frame.cons, [Binding.pair, :v1, :p1],
                                           [Env.first, [set, [Env.cons, :f2, :e2], :v2, :p2]]],
                              [Env.rest, [set, [Env.cons, :f2, :e2], :v2, :p2]]])
axiom_set3 = Axiom([set, [Env.cons, Frame.empty, :e], :v, :p], [Env.cons, Frame.empty, [set, :e, :v, :p]])

axiom_add_binding = Axiom([add_binding, :f, :v, :p], [Frame.cons, [Binding.pair, :v, :p], :f])
axiom_def = Axiom([def, :e, :v, :p], [Env.first!, :e, [add_binding, [Env.first, :e], :v, :p]])
axiom_extend = Axiom([extend, :e, :v, :p], [def, [Env.cons, Frame.empty, :e], :v, :p])

axiom_var_s = Axiom([s_eval, [var, :v], :e], :e)
axiom_var_f = Axiom([f_eval, [var, :v], :e], [lookup, :e, :v])
axiom_definition_s = Axiom([s_eval, [definition, :v, :p], :e], [def, [s_eval, :p, :e], :v, [f_eval, :p, :e]])
axiom_definition_f = Axiom([f_eval, [definition, :v, :p], :e], [f_eval, :p, :e])
axiom_assignment_s = Axiom([s_eval, [assignment, :v, :p], :e], [set, [s_eval, :p, :e], :v, [f_eval, :p, :e]])
axiom_assignment_f = Axiom([f_eval, [assignment, :v, :p], :e], [f_eval, :p, :e])
axiom_proc_s = Axiom([s_eval, [procedure, :v, :p], :e], :e)
axiom_proc_f = Axiom([f_eval, [procedure, :v, :p], :e], [procedure, :v, :p])
axiom_func_s1 = Axiom([s_eval, [func_call, :p1, :p2], :e],
                      [s_eval, [func_call, [f_eval, :p1, :e], [f_eval, :p2, :e]], :e])
axiom_func_f1 = Axiom([f_eval, [func_call, :p1, :p2], :e],
                      [f_eval, [func_call, [f_eval, :p1, :e], [f_eval, :p2, :e]], :e])
axiom_func_s2 = Axiom([s_eval, [func_call, [procedure, :v, :p1], :p2], :e],
                      [Env.rest, [s_eval, :p1, [extend, :e, :v, :p2]]])
axiom_func_f2 = Axiom([f_eval, [func_call, [procedure, :v, :p1], :p2], :e], [f_eval, :p1, [extend, :e, :v, :p2]])
axiom_seq_s = Axiom([s_eval, [Seq.cons, :p1, :p2], :e], [s_eval, :p2, [s_eval, :p1, :e]])
axiom_seq_f = Axiom([f_eval, [Seq.cons, :p1, :p2], :e], [f_eval, :p2, [s_eval, :p1, :e]])

axioms = Dict()
axioms_lookup = [axiom_lookup1, axiom_lookup2, axiom_lookup3]
axioms_set = [axiom_set1, axiom_set2, axiom_set3]
axioms[:base_env] = [Stream.axioms, Binding.axioms, Frame.axioms, Env.axioms]
axioms[:env] = [axioms_lookup, axioms_set, axiom_add_binding, axiom_def, axiom_extend]
axioms_v = [axiom_var_s, axiom_var_f]
axioms_d = [axiom_definition_s, axiom_definition_f]
axioms_a = [axiom_assignment_s, axiom_assignment_f]
axioms_p = [axiom_proc_s, axiom_proc_f]
axioms_f = [axiom_func_s1, axiom_func_f1, axiom_func_s2, axiom_func_f2]
axioms_s = [axiom_seq_s, axiom_seq_f]
axioms[:prog] = [axioms_v, axioms_d, axioms_a, axioms_p, axioms_f, axioms_s];

## The Repl facility, note that for the program structure parsing, we rely on the julia parser;
## and we also need to read the variable name from string.
include("Repl.jl")
r = REPL(Env.empty, Seq.empty)

function object_init(dict :: Dict, name, num :: Int64, class :: Class)
    dict[Symbol(name, num)] = Object(string(name, num), class)
end

function object_init(dict :: Dict, name, class :: Class)
    dict[Symbol(name)] = Object(string(name), class)
end

function object_init(dict :: Dict, name_dict :: Dict, num :: Int64)
    for k in keys(name_dict)
        object_init(dict, k, name_dict[k])
        for i in 1:num
            object_init(dict, k, i, name_dict[k])
        end
    end
end

function object_init(dict :: Dict, name_dict :: Dict)
    for k in keys(name_dict)
        object_init(dict, k, name_dict[k])
    end
end

type Memory
    class :: Class
    mclass :: Class

    encode :: DFunction
    decode :: DFunction

    axioms :: Array{Axiom, 1}
end

function Base.show(io :: IO, m :: Memory)
    print(io, string("Module.Memory(", m.class.class_name, "->", m.mclass.class_name, ")"))
end

function Memory(class :: Class, mclass :: Class, symbols = [:s, :m])
    encode = DFunction("encode", [class], mclass)
    decode = DFunction("decode", [mclass], class)

    sym = symbols[1]
    msym = symbols[2]
    axiom1 = Axiom([decode, [encode, sym]], sym)
    axiom2 = Axiom([encode, [decode, msym]], msym)
    axioms = [axiom1, axiom2]

    Memory(class, mclass, encode, decode, axioms)
end

function Memory(class :: Class, type_len :: Int64, symbols = [:s, :m])
    memory_name = string("Memory", "_", class.class_name)
    mclass = Class(memory_name, type_len)
    Memory(class, mclass, symbols)
end

## @time l = List(Sensor, 70)

naive_randomize = distribute(function(o :: Object) randn!(o.value) end)

function another_naive_randomize(tree :: Tree)
    foreach(another_naive_randomize, tree.subtrees)
    if typeof(tree.op) == Object
        randn!(tree.op.value)
    end
    tree
end

# UChar = Class("Char", 8)
# Variable = Class("Variable", 100)
# Prog = Class("Prog", 100)
# Bindingc = Class("Binding", 100)
# Binding = Pair(Variable, Prog, Bindingc, [:b, :v, :p])
# Framec = Class("Frame", 100)
# Frame = List(Bindingc, Framec, [:f1, :b, :f2])
# Envc = Class("Env", 200)
# Env = List(Framec, Envc, [:e1, :f, :e2])

dict = Dict()
ndict = Dict()
ndict[:u] = UChar
ndict[:v] = Variable
ndict[:p] = Prog
ndict[:b] = Bindingc
ndict[:f] = Framec
ndict[:e] = Envc
object_init(dict, ndict, 2)

init_axioms!(axioms, dict)
while true
    anti_train_axioms!(axioms, 1000, another_naive_randomize)
    train_axioms!(axioms, dict, 10000, naive_randomize)
    save("infinity.jld", Interpreter)
end

ext: jl

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

@testset "test tet4 stiffness matrix" begin
    el = Element(Tet4, [1, 2, 3, 4])
    el["youngs modulus"] = 96.0
    el["poissons ratio"] = 1/3
    x1 = [2.0, 3.0, 4.0]
    x2 = [6.0, 3.0, 2.0]
    x3 = [2.0, 5.0, 1.0]
    x4 = [4.0, 3.0, 6.0]
    u1 = u2 = u3 = u4 = zeros(3)
    el["geometry"] = Vector{Float64}[x1, x2, x3, x4]
    u = Vector{Float64}[u1, u2, u3, u4]
    pr = Problem(Elasticity, "tet4", 3)
    as = Assembly()
    assemble!(as, pr, el, 0.0)
    Kt = full(as.K)
    Kt_expected = [
         149.0   108.0   24.0   -1.0    6.0   12.0  -54.0   -48.0    0.0  -94.0   -66.0  -36.0
         108.0   344.0   54.0  -24.0  104.0   42.0  -24.0  -216.0  -12.0  -60.0  -232.0  -84.0
          24.0    54.0  113.0    0.0   30.0   35.0    0.0   -24.0  -54.0  -24.0   -60.0  -94.0
          -1.0   -24.0    0.0   29.0  -18.0  -12.0  -18.0    24.0    0.0  -10.0    18.0   12.0
           6.0   104.0   30.0  -18.0   44.0   18.0   12.0   -72.0  -12.0    0.0   -76.0  -36.0
          12.0    42.0   35.0  -12.0   18.0   29.0    0.0   -24.0  -18.0    0.0   -36.0  -46.0
         -54.0   -24.0    0.0  -18.0   12.0    0.0   36.0     0.0    0.0   36.0    12.0    0.0
         -48.0  -216.0  -24.0   24.0  -72.0  -24.0    0.0   144.0    0.0   24.0   144.0   48.0
           0.0   -12.0  -54.0    0.0  -12.0  -18.0    0.0     0.0   36.0    0.0    24.0   36.0
         -94.0   -60.0  -24.0  -10.0    0.0    0.0   36.0    24.0    0.0   68.0    36.0   24.0
         -66.0  -232.0  -60.0   18.0  -76.0  -36.0   12.0   144.0   24.0   36.0   164.0   72.0
         -36.0   -84.0  -94.0   12.0  -36.0  -46.0    0.0    48.0   36.0   24.0    72.0  104.0]
    if !isapprox(Kt, Kt_expected)
        info("Test failed")
        info("Kt_expected")
        dump(Kt_expected)
        info("Kt")
        dump(Kt)
    end
    @test isapprox(Kt, Kt_expected)
end


text:
s = """
function f(x)
    return x+1
    println(x)
end
"""
msgs = lintstr(s)
@test length(msgs) == 1
@test msgs[1].code == :W641
@test contains(msgs[1].message, "unreachable code after return")


text:
export DSSAIntegrator

mutable struct DSSAIntegrator{F,uType,tType,P,S,CB,SA,OPT,TS,chanType,chanS} <: DiffEqBase.DEIntegrator{SSAStepper,Nothing,uType,tType}
    f::F
    u::uType
    t::tType
    tprev::tType
    p::P
    sol::S
    i::Int
    tstop::tType
    cb::CB
    saveat::SA
    save_everystep::Bool
    save_end::Bool
    cur_saveat::Int
    opts::OPT
    tstops::TS
    tstops_idx::Int
    u_modified::Bool
    keep_stepping::Bool  # false if should terminate a simulation
    de_chan::chanType
    chan_sol::chanS
    delayjumpsets::DelayJumpSet
    save_delay_channel::Bool
end

mutable struct DSSASolution{uType,tType,chanS}
    u::uType
    t::tType
    channel::chanS
    odesol::DiffEqBase.ODESolution
end
function DSSASolution(ode_sol::DiffEqBase.ODESolution, channel)
    @unpack u, t = ode_sol
    DSSASolution{typeof(u),typeof(t),typeof(channel)}(u, t, channel, ode_sol)
end

function DiffEqBase.u_modified!(integrator::DSSAIntegrator, bool::Bool)
    integrator.u_modified = bool
end

function DiffEqBase.__solve(djump_prob::DelayJumpProblem,
    alg::SSAStepper;
    kwargs...)
    integrator = init(djump_prob, alg; kwargs...)
    solve!(integrator)
    if integrator.save_delay_channel
        DSSASolution(integrator.sol, integrator.chan_sol)
    else
        integrator.sol
    end
end

## Initiate the Problem

function DiffEqBase.__init(djump_prob::DelayJumpProblem,
    alg::SSAStepper;
    save_start=true,
    save_end=true,
    seed=nothing,
    alias_jump=Threads.threadid() == 1,
    saveat=nothing,
    callback=nothing,
    tstops=eltype(djump_prob.prob.tspan)[],
    numsteps_hint=100)
    if !(djump_prob.prob isa DiscreteProblem)
        error("SSAStepper only supports DiscreteProblems.")
    end
    @assert isempty(djump_prob.jump_callback.continuous_callbacks) # check if continuous_callbacks is empty
    ## Set up the cb which encodes the aggregator information
    if alias_jump
        cb = djump_prob.jump_callback.discrete_callbacks[end]
        if seed !== nothing
            Random.seed!(cb.condition.rng, seed)
        end
    else
        cb = deepcopy(djump_prob.jump_callback.discrete_callbacks[end])
        if seed === nothing
            Random.seed!(cb.condition.rng, seed_multiplier() * rand(UInt64))
        else
            Random.seed!(cb.condition.rng, seed)
        end
    end

    opts = (callback=CallbackSet(callback),)
    prob = djump_prob.prob
    de_chan0 = convert(Vector{Vector{typeof(prob.tspan[1])}}, djump_prob.de_chan0)
    if save_start
        t = [prob.tspan[1]]
        u = [copy(prob.u0)]
        chan_sol = [deepcopy(de_chan0)] # DelaySSA: build chan_sol
    else
        t = typeof(prob.tspan[1])[]
        u = typeof(prob.u0)[]
        chan_sol = typeof(de_chan0)[] # DelaySSA: build chan_sol
    end

    sol = DiffEqBase.build_solution(prob, alg, t, u, dense=false, calculate_error=false, destats=DiffEqBase.DEStats(0), interp=DiffEqBase.ConstantInterpolation(t, u))
    save_everystep = any(cb.save_positions)

    if typeof(saveat) <: Number
        _saveat = prob.tspan[1]:saveat:prob.tspan[2]
    else
        _saveat = saveat
    end

    if _saveat !== nothing && !isempty(_saveat) && _saveat[1] == prob.tspan[1]
        cur_saveat = 2
    else
        cur_saveat = 1
    end

    if _saveat !== nothing && !isempty(_saveat)
        sizehint!(u, length(_saveat) + 1)
        sizehint!(t, length(_saveat) + 1)
    elseif save_everystep
        sizehint!(u, numsteps_hint)
        sizehint!(t, numsteps_hint)
    else
        sizehint!(u, save_start + save_end)
        sizehint!(t, save_start + save_end)
    end

    integrator = DSSAIntegrator(prob.f, copy(prob.u0), prob.tspan[1], prob.tspan[1], prob.p, sol, 1, prob.tspan[1], cb, _saveat, save_everystep, save_end, cur_saveat, opts, tstops, 1, false, true, deepcopy(de_chan0), chan_sol, djump_prob.delayjumpsets, djump_prob.save_delay_channel)

    cb.initialize(cb, integrator.u, prob.tspan[1], integrator)
    DiffEqBase.initialize!(opts.callback, integrator.u, prob.tspan[1], integrator)
    integrator
end

function DiffEqBase.add_tstop!(integrator::DSSAIntegrator, tstop)
    if tstop > integrator.t
        future_tstops = @view integrator.tstops[integrator.tstops_idx:end]
        insert_index = integrator.tstops_idx + searchsortedfirst(future_tstops, tstop) - 1
        Base.insert!(integrator.tstops, insert_index, tstop)
    end
end

function DiffEqBase.solve!(integrator::DSSAIntegrator)

    end_time = integrator.sol.prob.tspan[2]
    while should_continue_solve(integrator) # It stops before adding a tstop over
        step!(integrator)
    end
    last_t = integrator.t
    integrator.t = end_time

    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        saveat_function_direct_method_test!(integrator, last_t)
    else
        saveat_function!(integrator, last_t)
    end

    if integrator.save_end && integrator.sol.t[end] != end_time
        saveat_end_function!(integrator, last_t)
    end
    DiffEqBase.finalize!(integrator.opts.callback, integrator.u, integrator.t, integrator)
end

function saveat_end_function!(integrator, prev_t)
    # save last t
    end_time = integrator.sol.prob.tspan[2]
    push!(integrator.sol.t, end_time)
    t_final_gap = end_time - prev_t

    # save last u
    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        update_delay_at_tstop_test!(integrator.cb.affect!, integrator, integrator.p, prev_t, t_final_gap)
    end
    push!(integrator.sol.u, copy(integrator.u))

    # save last de_chan
    if integrator.save_delay_channel
        last_chan = deepcopy(integrator.de_chan)
        if !(typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation)
            shift_delay_channel!(last_chan, t_final_gap)
            update_delay_channel!(last_chan)
        end
        push!(integrator.chan_sol, last_chan)
    end
end

"""
    function saveat_function!(integrator, prev_t)

Note that this function does not change u and de_chan, but only takes the snapshots of u and de_chan accordingly
"""
function saveat_function!(integrator, prev_t)
    if integrator.saveat !== nothing && !isempty(integrator.saveat)
        # Split to help prediction
        last_saved_t = prev_t
        while integrator.cur_saveat < length(integrator.saveat) && (integrator.saveat[integrator.cur_saveat] < integrator.t)
            tgap = integrator.saveat[integrator.cur_saveat] - last_saved_t
            push!(integrator.sol.t, integrator.saveat[integrator.cur_saveat])
            push!(integrator.sol.u, copy(integrator.u))
            if integrator.save_delay_channel
                prev_de_chan = deepcopy(integrator.de_chan)
                shift_delay_channel!(prev_de_chan, tgap)
                push!(integrator.chan_sol, prev_de_chan)
            end
            last_saved_t = integrator.saveat[integrator.cur_saveat]
            integrator.cur_saveat += 1
        end
    end
end

"""
    function saveat_function_direct_method_test!(integrator, prev_t)

Note that this function does not change u and de_chan for the next jump, but only takes the snapshots of u and de_chan accordingly. Because in next `execute_jumps`, shadow_integrator will overwrite integrator and when t = tend, the initial state will be recovered
"""
function saveat_function_direct_method_test!(integrator, prev_t)
    # Special to Direct method
    if integrator.saveat !== nothing && !isempty(integrator.saveat)
        last_saved_t = prev_t
        # save for last step and copied = false
        changed = false
        aggregator = integrator.cb.affect!
        if (integrator.cur_saveat < length(integrator.saveat)) && (integrator.saveat[integrator.cur_saveat] < integrator.t)
            shadow_integrator = aggregator.shadow_integrator
            shadow_integrator.u = copy(integrator.u)
            shadow_integrator.de_chan = deepcopy(integrator.de_chan)
            aggregator.copied = true
            changed = true
        end
        while integrator.cur_saveat < length(integrator.saveat) && (integrator.saveat[integrator.cur_saveat] < integrator.t)

            tgap = integrator.saveat[integrator.cur_saveat] - last_saved_t
            push!(integrator.sol.t, integrator.saveat[integrator.cur_saveat])

            update_delay_at_tstop_test!(aggregator, shadow_integrator, integrator.p, last_saved_t, tgap)
            push!(integrator.sol.u, copy(shadow_integrator.u))

            if integrator.save_delay_channel
                prev_de_chan = deepcopy(shadow_integrator.de_chan)
                push!(integrator.chan_sol, prev_de_chan)
            end

            last_saved_t = integrator.saveat[integrator.cur_saveat]
            integrator.cur_saveat += 1
        end
        if aggregator.copied && changed
            tend = integrator.t == integrator.sol.prob.tspan[2] ? integrator.sol.prob.tspan[2] : aggregator.next_jump_time
            last_tgap = tend - last_saved_t
            update_delay_at_tstop_test!(aggregator, shadow_integrator, integrator.p, last_saved_t, last_tgap)
        end
    end
end

# The Jump aggregators should not register the next jump through add_tstop! for SSAIntegrator
# such that we can achieve maximum performance
@inline function register_next_jump_time!(integrator::DSSAIntegrator, p::AbstractDSSAJumpAggregator, t)
    integrator.tstop = p.next_jump_time
    nothing
end

## Delay SSA modify
function DiffEqBase.step!(integrator::DSSAIntegrator)
    integrator.tprev = integrator.t
    next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)

    doaffect = false
    if !isempty(integrator.tstops) &&
       integrator.tstops_idx <= length(integrator.tstops) &&
       integrator.tstops[integrator.tstops_idx] < next_jump_time

        integrator.t = integrator.tstops[integrator.tstops_idx]
        integrator.tstops_idx += 1
    else
        integrator.t = integrator.tstop
        doaffect = true # delay effect until after saveat
    end

    if typeof(integrator.cb.affect!) <: DelayDirectJumpAggregation
        saveat_function_direct_method_test!(integrator, integrator.tprev)
    else
        saveat_function!(integrator, integrator.tprev)
    end

    # FP error means the new time may equal the old if the next jump time is
    # sufficiently small, hence we add this check to execute jumps until
    # this is no longer true.
    while integrator.t == integrator.tstop
        doaffect && integrator.cb.affect!(integrator)
    end

    if !(typeof(integrator.opts.callback.discrete_callbacks) <: Tuple{})
        discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!(integrator, integrator.opts.callback.discrete_callbacks...)
    else
        saved_in_cb = false
    end

    !saved_in_cb && savevalues!(integrator)

    nothing
end

function DiffEqBase.savevalues!(integrator::DSSAIntegrator, force=false)
    saved, savedexactly = false, false

    # No saveat in here since it would only use previous values,
    # so in the specific case of DSSAStepper it's already handled

    if integrator.save_everystep || force
        saved = true
        savedexactly = true
        push!(integrator.sol.t, integrator.t)
        push!(integrator.sol.u, copy(integrator.u))
        if integrator.save_delay_channel
            push!(integrator.chan_sol, deepcopy(integrator.de_chan))
        end
    end

    saved, savedexactly
end

function should_continue_solve(integrator::DSSAIntegrator)
    end_time = integrator.sol.prob.tspan[2]

    # we continue the solve if there is a tstop between now and end_time
    has_tstop = !isempty(integrator.tstops) &&
                integrator.tstops_idx <= length(integrator.tstops) &&
                integrator.tstops[integrator.tstops_idx] < end_time

    # we continue the solve if there will be a jump between now and end_time
    has_jump = integrator.t < integrator.tstop < end_time

    integrator.keep_stepping && (has_jump || has_tstop)
end

num_constant_rate_jumps(aggregator::AbstractDSSAJumpAggregator) = length(aggregator.rates)


text:
########################################
## File Name: nominal_policies_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/14
## Description: Test code for src/nominal_policies.jl
########################################

using LinearAlgebra

@testset "Nominal Policies Test" begin

# Localization Task
dtc = 0.01;
dto = 0.5;
dtexec = dto*0.8;
Q = Matrix(1.0I, 2, 2);
Cu = Matrix(1.0I, 2, 2);
transPos = TransModel_Pos();
ukfModel = UKFPosRange();
gradientPolicy = GradientMultiTargetLocalizationPolicy();

bPos = BelMvNormal(0.,zeros(2),100*Matrix(1.0I, 2, 2));
bVecPos = [bPos for ii = 1:100];
p_robot = PhysPos(0.,[5.,7.]);
sV = AugState(p_robot,bVecPos);
u_param_min_2D = [-2.,-2.];
u_param_max_2D = [2.,2.];

uPos = control_nominal(gradientPolicy,ukfModel,sV,u_param_min_2D,u_param_max_2D,dto,Q)
@test uPos.t == 0.0
@test isapprox(uPos.vel, [0.556133, 0.831093], rtol=1e-4)

# Manipulation Task
dtc = 0.02;
dto = 0.2;
dtexec = dto*0.8;
QMP = [Matrix(1.0I, 6, 6) zeros(6,5); zeros(5,6) 0.00*Matrix(1.0I, 5, 5)];
RMP = Matrix(1.0I, 9, 9);
CsMP = QMP;
CuMP = Matrix(1.0I, 3, 3);
pcontrolPolicy = ManipulatePositionControlPolicy();

xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.2);
uMP = MControl2D(0.,6.,7.,8.);
bMP = BelMvNormal(xMP.t,vec(xMP),diagm([10.,10.,1.,10.,10.,10.,1.,1.,10.,10.,0.1]));
u_param_min_3D = [-3.,-3.,-3.];
u_param_max_3D = [3.,3.,3.];
posGain = 1.0;
rotGain = 0.2;

uMP = control_nominal(pcontrolPolicy,bMP,u_param_min_3D,u_param_max_3D,posGain,rotGain)
@test uMP.t == 0.0
@test uMP.fx == -0.1
@test uMP.fy == -0.2
@test uMP.tr == 0.5283185307179586

end


text:
module PolyhedralOmega
using LinearAlgebra, IterTools, MultivariatePolynomials, TypedPolynomials, AbstractAlgebra, TaylorSeries
#using LazySet

include("SymbolicCone.jl")
include("./SmithNormalForm/src/SmithNormalForm.jl")

function macmahon(A::Matrix{Int64}, b::Vector{Int64})
    sizeA = size(A)
    Id = Matrix(1I, sizeA[2], sizeA[2])
    V = vcat(Id, A)
    q = append!(zeros(Float64, sizeA[2]), -b)
    #o = zeros(Bool, size(q,1))
    o = zeros(Bool, sizeA[2])
    return SymbolicCone(V, q, o)
end

function flip(C::SymbolicCone)
    V = C.V
    s = 1
    for i = 1:size(V, 2)
        for j = 1:size(V, 1)
            if (V[j, i] != 0)
                forward = (sign(V[j, i]) == 1)
                if (!forward)
                    s = s * (-1)
                    V[:, i] = (-1) * V[:, i]
                    C.o[i] = xor(C.o[i], true)
                end
                break
            end
        end
    end
    C.sign = s
    return C
end

function elimLastCoordinate(C::SymbolicCone)
    V = C.V
    n = size(V, 1)
    k = size(V, 2)
    q = C.q
    o = C.o
    res = []
    II = []

    Iplus = filter(x -> V[n, x] > 0, Vector(1:k)) #one of them should have zero eq or not?
    Iminus = filter(x -> V[n, x] < 0, Vector(1:k))
    #if (q[n] > 0 || (q[n] == 0 && size(Iplus, 1) <= size(Iminus, 1) && size(Iplus, 1) >= 1) )
    #if (q[n] == 0)
    #    println("Iequal")
    #    return [SymbolicCone(prim(V[1:n-1, :]), q[1:n-1], o[1:k])]
    #end
    if (q[n] < 0)
        #res = [SymbolicCone(prim(V[1:n-1, :]), q[1:n-1], o[1:k])]
        println("Iplus")
        II = Iplus
    else
        println("Iminus")
        II = Iminus
        res = [SymbolicCone(prim(V[1:n-1, :]), q[1:n-1], o[1:k])]
    end
    #Cprime if q[n] is zero
    println("II: ", II)

    sgnq = q[n]
    for j in II
        L = []
        for i in 1:k
            if i == j
                push!(L, (sgnq >= 0 ? -1 : 1)*V[:,j])
            else
                push!(L, (sgnq >= 0 ? 1 : -1)*( (V[n,i] * V[:,j]) + (-V[n,j]*V[:,i]) ))
            end
        end
        println("L: ", L)
        new_o = o[j] != false ? vcat(o[1:j-1], false, o[j+1:end]) : deepcopy(o)
        new_V = prim(hcat(L...))
        new_V = new_V[1:(size(new_V, 1)-1),:]

        new_q = map(+, q, (-q[n] / (V[n, j])) * V[:, j])
        new_q = new_q[1:n-1]
        new_cone = SymbolicCone(new_V, new_q, new_o)

        print("\nnew cone:")
        println(new_cone)

        flip_res = flip(new_cone)
        push!(res, flip_res)

        print("\nflip:")
        println(flip_res)
    end
    return res
end

function prim_v(v::Array{Int64})
    d = abs(gcd(v))
    if (d == 1)
        return v
    else
        return reduce(vcat, map(vi -> floor(Int, vi / d), v))
    end
end

function prim(V::Matrix{Int64})
    return Matrix{Int64}(transpose(reduce(hcat, map(i -> prim_v(V[i, :]), Vector(1:size(V, 1))))))
end

function eliminateCoordinates(C::SymbolicCone, k::Int)
    res = elimLastCoordinate(C)
    innerres = []

    println("first elc: ", res)
    for i = 1:(k-1)#becuse we call eliminate last coordinate before
        for j = 1:length(res)
            eres =  elimLastCoordinate(res[j])
            append!(innerres, eres)
            println("j: ", j, " eres: ", eres)
        end
        res = innerres
        innerres = []
    end

    return res
end

function enumerateFundamentalParallelePiped(C::SymbolicCone)
    SMFRes = SmithNormalForm.smith(C.V)
    S = SmithNormalForm.diagm(SMFRes)
    Uinv = inv(SMFRes.S)
    Winv = inv(SMFRes.T)
    dimension = size(C.V, 2) # num of rows
    ambientDimension = size(C.V, 1) # num of cols
    #println("C.V: ", C.V)
    #println("\nUinv: ", Uinv, "\nWinv: ", Winv)
    #println("size 1: ", size(S,1), "\nsize 2: ", size(S,2))

    diagonals = Int64[]
    for i in 1:dimension
        if(i <= size(S,1) && i <= size(S,2))
            #println("i: ", i)
            push!(diagonals, S[i,i])
        end
    end
    #println("Diagonals: ", diagonals)

    lastDiagonal = diagonals[end]

    # sprime = [Integer(sk / si) for si in s]
    primeDiagonals = Int64[]
    for d in diagonals
        push!(primeDiagonals, Int64(lastDiagonal/d))
    end
    #println("Prime diagonals: ", primeDiagonals)

    # qhat = Uinv * q
    apex = C.q
    #println("q: ", apex, "\nV: ", C.V)

    qhat = Uinv*apex
    #println("qhat: ", qhat)

    # Wprime
    Wprime = [Winv[j,i]*primeDiagonals[i] for i = 1:dimension, j = 1:dimension] #Winv * primeDiagonals
    #println("wprime: ", Wprime)

    tmpWprime = deepcopy(Wprime)

    #println("ftmpwprime: ", tmpWprime)
    for i in 1:dimension
        tmpWprime[dimension-i+1,:] = Wprime[i,:]
    end

    #println("tmpwprime: ", tmpWprime)
    #Wprime = tmpWprime

    # qtrans
    qtrans = [sum([-Wprime[j,i] * qhat[i] for i = 1:dimension]) for j = 1:dimension]
    #println("qtrans: ", qtrans)

    #qfrac
    qfrac = [qtrans[i] - floor.(Int, qtrans[i]) for i = 1:dimension]
    #println("qfrac: ", qfrac)

    #qint
    qint = [ floor.(Int, qi) for qi in qtrans ]
    #println("qint: ", qint)

    #qsummand
    qsummand = [Int64(qi) for qi in (lastDiagonal*apex + C.V*qfrac) ]
    #println("qsummand", qsummand)

    #println("dim: ", dimension, " adim: ", ambientDimension)
    #println("o: ", C.o, " qfrac: ", qfrac)
    #openness
    openness = [ (qfrac[j] == 0 ? C.o[j] : 0) for j in 1:ambientDimension]
    #println("openness: ", openness)

    #bigP
    #res1 = [[1:1:diagonals[i];] for i= 1:dimension]
    #println("res1: ", res1)

    # CartesianProduct( *[xrange(s[i]) for i in 1:k] )
    L = []
    P = []
    for v in IterTools.product([1:diagonals[i] for i in 1:dimension]...)
        push!(P, v)
        innerRes = []
        j = 1
        for qj in qint
            inner = 0
            i = 1
            for vi in v
                inner += Wprime[i,j] * vi
                i += 1
            end
            inner += qj
            inner = inner % lastDiagonal

            if inner == 0 && C.o[j]
                inner = lastDiagonal
            end
            append!(innerRes, inner)
            j += 1
        end
        #println("inneerres: ", innerRes)

        outerRes = []
        for l in 1:ambientDimension
            outer = 0
            j = 1
            for innerResi in innerRes
                outer += C.V[l,j] * innerResi
                j += 1
            end
            append!(outerRes, outer) # outerRes is an integral vector
        end
        #println("outerRes: ", outerRes)
        #push!(L, tuple(collect( ((ai + bi) / lastDiagonal) for (ai,bi) in collect(zip(outerRes, qsummand)) )))
        push!(L, collect( Int64((ai + bi) / lastDiagonal) for (ai,bi) in collect(zip(outerRes, qsummand)) ))
    end
    C.fp = L
end

function computeRationalFunction(C::SymbolicCone)
    numOfVars = size(C.q,1)
    VV = [string("x_",i) for i in 1:numOfVars]
    S, X = PolynomialRing(QQ, VV, ordering=:deglex)
    num = 0
    #println("VV: ", VV, "\nQQ: ", QQ, "\nS: ", S, "\nX: ", X, "\n")
    for p in C.fp
        tmp = 1
        for i in 1:length(p)
            tmp *=X[i]^p[i]
        end
        num += tmp
    end
    #println("Numerator:", num)
    den = 1
    for j in 1:size(C.V,2)
        tmp = 1
        for i in 1:size(C.V,1)
            tmp *=X[i]^C.V[:,j][i]
        end
        den *= 1 - tmp
    end
    #println("Denominator:", den)
    ratfun = num//den
    #println("ratfun:", ratfun)

    C.ratfun = ratfun
end

function solve(A::Matrix{Int64}, b::Vector{Int64})
    @polyvar RationalFunctions
    ratfun = 0
    C = macmahon(A, b)

    print("MacMahon cone: ")
    println(C)

    ListOfSymbolicCones = eliminateCoordinates(C, size(C.V, 1) - size(C.V, 2))
    #println("eliminate coordinates: ", ListOfSymbolicCones)
    println("Result Cones: ")
    for cone in ListOfSymbolicCones
        #enumerateFundamentalParallelePiped(cone)
        #computeRationalFunction(cone)
        #ratfun += cone.sign*cone.ratfun
        println(cone)
        #println("fp: ", cone.fp)
        #println("ratfun: ", cone.ratfun)
        #println()
    end
    return ratfun

end

end

max_stars_repo_name: kevseryolcu/PolyhedralOmega

max_issues_repo_name: kevseryolcu/PolyhedralOmega

max_forks_repo_name: kevseryolcu/PolyhedralOmega


text:
using PackmolInputCreator
using Test

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


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "RelicToolkit"
version = v"0.4.0-613-g613795f7-3"

# Collection of sources required to build RelicToolkit
sources = [
    "https://github.com/blenessy/relic.git" => "613795f7481bf7806ba09bf248ab36aec599f4ab",
]

# Bash recipe for building across all platforms
script = raw"""
if [ "${target#*-apple-}" != "$target" ]; then
    opsys=MACOSX
elif [ "${target%-mingw32}" != "$target" ]; then
    opsys=WINDOWS
else
    opsys=LINUX
fi
cd "$WORKSPACE/srcdir/relic"

cat >>src/relic_core.c <<EOF1
#include "relic_bn.h"
#include "relic_md.h"
#include "relic_fp.h"
#include "relic_pc.h"
const size_t JL_RLC_BN_SIZE = RLC_BN_SIZE;
const size_t JL_RLC_MD_LEN = RLC_MD_LEN;
const size_t JL_BN_ST_SIZE  = sizeof(bn_st);
const size_t JL_DIG_T_SIZE  = sizeof(dig_t);
const size_t JL_FP_ST_SIZE  = sizeof(fp_st);
const size_t JL_FP2_ST_SIZE = sizeof(fp2_st);
const size_t JL_FP3_ST_SIZE = sizeof(fp3_st);
const size_t JL_G1_ST_SIZE  = sizeof(g1_st);
const size_t JL_G2_ST_SIZE  = sizeof(g2_st);
EOF1

cat >build.sh <<EOF2
set -exo pipefail
mkdir "build_\$1"
cd "build_\$1"
cmake -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_TOOLCHAIN_FILE="/opt/$target/${target}_gcc.toolchain" \
  -DALLOC=AUTO \
  -DARITH=\$2 \
  -DBENCH=0 \
  -DCOLOR=off \
  -DCHECK=off \
  -DCMAKE_BUILD_TYPE=RELEASE \
  -DCMAKE_POSITION_INDEPENDENT_CODE=on \
  -DCOMP="-O3 -funroll-loops -fomit-frame-pointer" \
  -DDOCUM=off \
  -DSEED=udev \
  -DEP_PLAIN=off \
  -DEP_SUPER=off \
  -DFP_METHD="BASIC;COMBA;COMBA;MONTY;LOWER;SLIDE" \
  -DFP_PMERS=off \
  -DFP_PRIME=\$3 \
  -DFP_QNRES=on \
  -DAMALG=on \
  -DFPX_METHD="INTEG;INTEG;LAZYR" \
  -DMD_METHD=SH256 \
  -DMULTI=PTHREAD \
  -DOPSYS="${opsys:-MACOSX}" \
  -DPP_EXT=LAZYR \
  -DPP_METHD="LAZYR;OATEP" \
  -DQUIET=off \
  -DSHLIB=on \
  -DSTBIN=off \
  -DSTLIB=off \
  -DTESTS=0 \
  -DTIMER=CYCLE \
  -DVERBS=off \
  -DWITH="BN;DV;FP;FPX;EP;EPX;EC;PP;PC;MD" \
  -DWSIZE="${nbits:-64}" \
  ..
make
make install
# workaround mingw install bug
if [ "${target%-mingw32}" != "$target" ]; then
    rm -f "$prefix/lib/librelic.$dlext.a"
    cp ./lib/librelic.$dlext "$prefix/lib/"
fi
mv "$prefix/lib/librelic.$dlext" "$prefix/lib/librelic_\$1.$dlext"
EOF2
sh build.sh gmp_pbc_bls381 gmp 381
sh build.sh gmp_pbc_bn254 gmp 254
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    #Linux(:aarch64, libc=:glibc),
    #Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    #Linux(:powerpc64le, libc=:glibc),
    #Linux(:i686, libc=:musl),
    #Linux(:x86_64, libc=:musl),
    #Linux(:aarch64, libc=:musl),
    #Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    #FreeBSD(:x86_64),
    # The resulting lib is called librelic_*.dll.a need to investigate the problem
    Windows(:i686),
    Windows(:x86_64),
]
#platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "librelic_gmp_pbc_bls381", :librelic_gmp_pbc_bls381),
    LibraryProduct(prefix, "librelic_gmp_pbc_bn254", :librelic_gmp_pbc_bn254),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/RelicToolkitBuilder

max_issues_repo_name: blenessy/RelicToolkitBuilder

max_forks_repo_name: blenessy/RelicToolkitBuilder


text:
abstract type Dump end

struct NoDump<: Dump end

struct LAMMPSTrj <: Dump
    filename::String
    freq::Int64
    is_velocity::Bool
    is_force::Bool
    io::IOStream
    function LAMMPSTrj( filename::String,
                        freq::Int64;
                        is_velocity=false,
                        is_force=false)
        io=open(filename,"w")
        new(filename::String,
            freq::Int64,
            is_velocity::Bool,
            is_force::Bool,
            io::IOStream
        )
    end
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Deep learning with Flux")

#

#
# # Intro to Flux.jl

#-

# In the previous course, we learned how machine learning allows us to classify data as apples or bananas with a single neuron. However, some of those details are pretty fiddly! Fortunately, Julia has a powerful package that does much of the heavy lifting for us, called [`Flux.jl`](https://fluxml.github.io/).
#
# *Using `Flux` will make classifying data and images much easier!*

#-

# ## Using `Flux.jl`
#
# We can get started with `Flux.jl` via:

## using Pkg; Pkg.add(["Flux", "Plots"])
using Flux, Plots

# #### Helpful built-in functions
#
# When working we'll `Flux`, we'll make use of built-in functionality that we've had to create for ourselves in previous notebooks.
#
# For example, the sigmoid function, σ, that we have been using already lives within `Flux`:

#nb ?σ
#jl @doc σ

#-

plot(σ, -5, 5, label="\\sigma", xlabel="x", ylabel="\\sigma\\(x\\)")

# Importantly, `Flux` allows us to *automatically create neurons* with the **`Dense`** function. For example, in the last notebook, we were looking at a neuron with 2 inputs and 1 output:
#
#
#
#  We could create a neuron with two inputs and one output via

model = Dense(2, 1, σ)

# This `model` object comes with places to store weights and biases:

model.W

#-

model.b

#-

typeof(model.W)

#-

x = rand(2)
model(x)

#-

σ.(model.W*x + model.b)

# Unlike in previous notebooks, note that `W` is no longer a `Vector` (1D `Array`) and `b` is no longer a number! Both are now stored in so-called `TrackedArray`s and `W` is effectively being treated as a matrix with a single row. We'll see why below.

#-

# Other helpful built-in functionality includes ways to automatically calculate gradients and also the cost function that we've used in the previous course -
#
# $$L(w, b) = \sum_i \left[y_i - f(x_i, w, b) \right]^2$$
#
# If you normalize by dividing by the total number of elements, this becomes the "mean square error" function, which in `Flux` is named **`Flux.mse`**.

methods(Flux.mse)

# ### Bringing it all together
#
# Load the datasets that contain the features of the apple and banana images.

using CSV, DataFrames

apples = DataFrame(CSV.File(datapath("data/apples.dat"), delim='\t', allowmissing=:none, normalizenames=true))
bananas = DataFrame(CSV.File(datapath("data/bananas.dat"), delim='\t', allowmissing=:none, normalizenames=true));

#-

x_apples  = [ [row.red, row.green] for row in eachrow(apples)]
x_bananas = [ [row.red, row.green] for row in eachrow(bananas)];

# Concatenate the x (features) together to create a vector of all our datapoints, and create the corresponding vector of known labels:

xs = [x_apples; x_bananas]
ys = [fill(0, size(x_apples)); fill(1, size(x_bananas))];

#-

model = Dense(2, 1, σ)

# We can evaluate the model (currently initialized with random weights) to see what the output value is for a given input:

model(xs[1])

# And of course we can examine the current loss value for that datapoint:

loss = Flux.mse(model(xs[1]), ys[1])

#-

typeof(loss)

# ### Backpropagation

model.W

#-

model.W.grad

#-

using Flux.Tracker
back!(loss)

#-

model.W.grad

# Now we have all the tools necessary to build a simple gradient descent algorithm!

#-

# ### The easy way
#
# You don't want to manually write out gradient descent algorithms every time! Flux, of course, also brings in lots of optimizers that can do this all for you.

#nb ?SGD
#jl @doc SGD

#-

#nb ?Flux.train!
#jl @doc Flux.train!

# So we can simply define our loss function, an optimizer, and then call `train!`. That's basic machine learning with Flux.jl.

model = Dense(2, 1, σ)
L(x,y) = Flux.mse(model(x), y)
opt = SGD(params(model))
Flux.train!(L, zip(xs, ys), opt)

# ## Visualize the result

contour(0:.1:1, 0:.1:1, (x, y) -> model([x,y])[].data, fill=true)
scatter!(first.(x_apples), last.(x_apples), label="apples")
scatter!(first.(x_bananas), last.(x_bananas), label="bananas")
xlabel!("mean red value")
ylabel!("mean green value")

max_stars_repo_path: Courses/Deep learning with Flux/1300.Intro-to-Flux.jl.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
include("../src/coefexp_ab.jl")

using Printf

function _printnumstr( num::BigFloat)
    str = @sprintf("%1.200e", num)
    println("parse(")
    println("    BigFloat," )
    println("    \"$(str[1:70])\" *")
    println("    \"$(str[71:140])\" *")
    println("    \"$(str[141:end])\"")
    print(") ")
end

function print_for_test(
    order,
    epsilon::Rational{BigInt},
    n_tau, dt::Rational{BigInt},
    list_j,
    str_plus
)

    prec=precision(BigFloat)
    setprecision(1024)
    coef_tau = [collect(0:n_tau/2 - 1); collect(-n_tau/2:-1)] # the result must be real
    par = CoefExpABRational(order, float(epsilon), coef_tau, float(dt) )
    ordp1=order+1
    println("# CoefExpAB order=$order epsilon=$epsilon n_tau =$n_tau dt=$dt")
    println("# tab_res_coef is AB coefficient for each value from 1 to n_tau")
    println("# this file is generated by gendatacoefexp.jl file")
    println("function get_coef_ab_for_test$str_plus()")
    println("    tabres = zeros(Complex{BigFloat}, $n_tau, $ordp1, $ordp1)")
    for j in list_j
        println("    tabres[ :, :, $j] .= [")
        for i_ell=1:n_tau
            print("    ")
            res = view(par.tab_coef, i_ell, :, j)
            for i = 1:ordp1
                if i <= j
                    _printnumstr(real(res[i]))
                    print("+ im * ")
                    _printnumstr(imag(res[i]))
                else
                    print(" 0")
                end
            end
            println("")
        end
        println("]")
    end
    println("    return tabres, $list_j")
    println("end")
    setprecision(prec)
end

print_for_test(15, big"1"//10^10, 32, big"1"//10000, [1,5,10,16],"10")


text:
using JuMP
using CPLEX
#using GLPK

function createLSP(params, nodes, demands, costs)

	#model = Model(GLPK.Optimizer)
	model = Model(CPLEX.Optimizer)

	#= Paramètres (Constantes) du problème

		demands[i][t] : d_{i,t}= demande du revendeur i au pas de temps t | i in [1,n], t in [1,l]=T
		nodes[i]["h"] : h_i = coût de stockage d'un produit pour le revendeur/fournisseur i | i in [0,n]
		node[i]["L"] : L_i le stock maximal du revendeur/fournisseur i | i in [0,n]
		params["u"] : u = coût de production d'un produit pour le fournisseur
		params["f"] : f = coût de set up (pour chaque période où le fournisseur décide de produire, il paye un coût f)
		params["C"] : C = Capacité de production max du fournisseur (=M dans le sujet)

	=#

	#= Variables de décision

		p[t] : p_t  (Variable de production): quantité produite à la période t
		y[t] : y_t  (Variable de lancement): variable binaire qui vaut 1 si une production est lancée à la période t, et 0 sinon

		I[i,j] : I_{i,t} (Variable de stockage): quantité en stock à la fin de la période t pour i.

		q[i,j] : q_{i,t} (Variable d’approvisionnement): quantité produite pour le revendeur i à la période t.

	=#

	#Variables Réel positive
	@variable(model, p[1:params["l"]] >= 0)
	@variable(model, I[0:params["n"], 0:params["l"]] >= 0)
	@variable(model, q[1:params["n"], 1:params["l"]] >= 0)

	#Variable binaire
	@variable(model, y[1:params["l"]], Bin)

	for i in 0:params["n"]
		@constraint(model, I[i, 0] == nodes[i]["L0"])
	end

	###Fonction objectif

	obj = 0
	for t in 1:params["l"]

		#On ajoute le cout de production de touts les produits au temps t
		obj += params["u"] * p[t]

		#Si on produit au temps t, on ajoute le cout de setup
		obj += params["f"] * y[t]

		#??? On ne prend pas le cout de stockage de fournisseur ?
		for i in 1:params["n"]
			obj += nodes[i]["h"] * I[i, t]
		end

	end
	@objective(model, Min, obj)

	###Contraintes
	for t in 1:params["l"]

		@constraint(model, I[0, t-1] + p[t] == I[0, t] + sum(q[i, t] for i in 1:params["n"])) #contraintes 1

		for i in 1:params["n"]
			@constraint(model, I[i, t-1] + q[i, t] == demands[i, t] + I[i, t]) #contraintes 2
			@constraint(model, I[i, t-1] + q[i, t] <= nodes[i]["L"]) #contraintes 5
		end

		@constraint(model, p[t] <= params["C"] * y[t]) #constraintes 3
		@constraint(model, I[0, t-1] <= nodes[0]["L"]) #contraintes 4

	end

	return model

end


text:
module TestThermodynamics

using Base.Test
using SoilTracers: Constants
using SoilTracers.PhysChem: c2k, k2c

@time @testset "Basic thermodynamics" begin
    temp1_C = 0.
    temp2_C = 100.
    temp1_K = c2k(temp1_C)
    temp2_K = c2k(temp2_C)
    @test temp1_K === temp1_C + Constants.T_0
    @test temp2_K === temp2_C + Constants.T_0

    @test temp1_C === k2c(temp1_K)
    @test temp2_C === k2c(temp2_K)
end

end


text:
issquare(A::AbstractMatrix) = ==(size(A)...)

mat(a::Vector) = reshape(a,length(a),1)


text:
# List of available algorithms for AIBECS time stepping

@testset "Time steppers" begin
    δt = ustrip(1.0u"yr" |> u"s")
    Nt = 10000 # number of time steps taken
    nt = length(T_all)
    n = nt * nb
    @unpack xgeo = p
    x₀ = xgeo * ones(n)

    Implicit_methods = [
                        AIBECS.crank_nicolson_step,
                        AIBECS.crank_nicolson_step!,
                        AIBECS.euler_backward_step,
                        AIBECS.euler_backward_step!,
                       ]
    @testset "Implicit time steps" begin
        @testset "$(string(sf))" for sf in Implicit_methods
            x = copy(x₀)
            for i in 1:Nt
                x .= sf(x, p, δt, F, ∇ₓF)
            end
            @test norm(x) / norm(F(x,p)) > ustrip(upreferred(1u"kyr"))
        end
    end

    Explicit_methods = [
                        AIBECS.euler_forward_step,
                        AIBECS.euler_forward_step!,
                       ]
    @testset "Explicit time steps" begin
        @testset "$(string(sf))" for sf in Explicit_methods
            x = copy(x₀)
            for i in 1:100Nt              # more smaller steps
                x .= sf(x, p, 0.01δt, F) # for explicit scheme
            end
            @test norm(x) / norm(F(x,p)) > ustrip(upreferred(1u"kyr"))
        end
    end

    @testset "Crank-Nicolson leapfrog" begin
        x, xᵢ, xᵢ₋₁ = copy(x₀), copy(x₀), copy(x₀)
        for i in 1:Nt
            x .= AIBECS.crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, δt, T, ∇ₓL, NL)
            xᵢ₋₁ .= xᵢ
            xᵢ .= x
        end
        @test norm(x) / norm(F(x,p)) > ustrip(upreferred(1u"kyr"))
        A⁺, A⁻ = AIBECS.crank_nicolson_leapfrog_step_A⁺_and_A⁻(p, δt, T, ∇ₓL)
        x, xᵢ, xᵢ₋₁ = copy(x₀), copy(x₀), copy(x₀)
        for i in 1:Nt
            x .= AIBECS.crank_nicolson_leapfrog_step(xᵢ, xᵢ₋₁, p, A⁺, A⁻, NL)
            xᵢ₋₁ .= xᵢ
            xᵢ .= x
        end
        @test norm(x) / norm(F(x,p)) > ustrip(upreferred(1u"kyr"))
    end

end


text:
"""
    erdos_renyi!(g::AbstractGraph, n_steps::Integer)

Erdos-Renyi style graph evolution, adds edges randomly at each step

Arguments
* `g`      : An instance of type AbstractGraph
* `n_steps`: Number of edges to add to the AbstractGraph
Returns
* `g`, updates `g` in-place
"""
function erdos_renyi!(g::AbstractGraph, n_steps::Integer)
    for t = 1:n_steps
        edge = choose_edge(g)
        add_edge!(g, edge)
    end

    finalize_observables!(g)
    return g
end

"""
    bohman_frieze!(g::AbstractGraph, n_steps::Integer; K::Integer=2)

Achlioptas process, implementation of Bohman-Frieze bounded size rule

Arguments
* `g`      : An instance of type AbstractGraph
* `n_steps`: Number of edges to add to the AbstractGraph
Keyword Arguments
* `K`      : Bounded size of clusters upon which to determine edge acceptance (default = 2)
Returns
* `g`, updates `g` in-place
"""
function bohman_frieze!(g::AbstractGraph, n_steps::Integer; K::Integer = 2)
    for t = 1:n_steps
        edge₁ = choose_edge(g)
        edge₂ = choose_edge(g, edge₁)

        if length(g.clusters[g.cluster_ids[edge₁[1]]]) < K &&
           length(g.clusters[g.cluster_ids[edge₁[2]]]) < K
            add_edge!(g, edge₁)
        else
            add_edge!(g, edge₂)
        end
    end

    finalize_observables!(g)
    return g
end

"""
    product_rule!(g::AbstractGraph, n_steps::Integer)

Achlioptas process, implementation of the product rule

Arguments
* `g`      : An instance of type AbstractGraph
* `n_steps`: Number of edges to add to the AbstractGraph
Returns
* `g`, updates `g` in-place
"""
function product_rule!(g::AbstractGraph, n_steps::Integer)
    for t = 1:n_steps
        edge₁ = choose_edge(g)
        edge₂ = choose_edge(g, edge₁)

        if g.cluster_ids[edge₁[1]] == g.cluster_ids[edge₁[2]]
            add_edge!(g, edge₁)
        elseif g.cluster_ids[edge₂[1]] == g.cluster_ids[edge₂[2]]
            add_edge!(g, edge₂)
        elseif length(g.clusters[g.cluster_ids[edge₁[1]]]) *
               length(g.clusters[g.cluster_ids[edge₁[2]]]) <
               length(g.clusters[g.cluster_ids[edge₂[1]]]) *
               length(g.clusters[g.cluster_ids[edge₂[2]]])
            add_edge!(g, edge₁)
        else
            add_edge!(g, edge₂)
        end
    end

    finalize_observables!(g)
    return g
end

"""
    stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer)

Achlioptas process, a probability based model for accepting edges

Arguments
* `g`      : An instance of type AbstractGraph
* `n_steps`: Number of edges to add to the AbstractGraph
Returns
* `g`, updates `g` in-place
"""
function stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer)
    for t = 1:n_steps
        edge₁ = choose_edge(g)
        edge₂ = choose_edge(g, edge₁)

        if g.cluster_ids[edge₁[1]] == g.cluster_ids[edge₁[2]]
            add_edge!(g, edge₁)
        elseif g.cluster_ids[edge₂[1]] == g.cluster_ids[edge₂[2]]
            add_edge!(g, edge₂)
        else
            C₁ =
                length(g.clusters[g.cluster_ids[edge₁[1]]]) +
                length(g.clusters[g.cluster_ids[edge₁[2]]])
            C₂ =
                length(g.clusters[g.cluster_ids[edge₂[1]]]) +
                length(g.clusters[g.cluster_ids[edge₂[2]]])
            p = (1 / C₁) / (1 / C₁ + 1 / C₂)

            if p > rand(g.rng)
                add_edge!(g, edge₁)
            else
                add_edge!(g, edge₂)
            end
        end
    end

    finalize_observables!(g)
    return g
end

"""
    stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer, q::Integer)

q-edge Achlioptas process, a probability based model for accepting edges

Arguments
* `g`      : An instance of type AbstractGraph
* `n_steps`: Number of edges to add to the AbstractGraph
* `q`      : Number of edges to evaluate
Returns
* `g`, updates `g` in-place
"""
function stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer, q::Integer)
    for t = 1:n_steps
        edges = [choose_edge(g)]
        for i = 1:q-1
            push!(edges, choose_edge(g, edges))
        end

        edge_added = false
        for edge in edges
            if g.cluster_ids[edge[1]] == g.cluster_ids[edge[2]]
                add_edge!(g, edge)
                edge_added = true
                break
            end
        end
        if edge_added
            continue
        end

        C = [
            (
                length(g.clusters[g.cluster_ids[edge[1]]]) +
                length(g.clusters[g.cluster_ids[edge[2]]])
            ) for edge in edges
        ]
        p = (1 ./ C) ./ sum(1 ./ C)
        R = rand(g.rng)
        P = 0
        for i = 1:q
            P += p[i]
            if R < P
                add_edge!(g, edges[i])
                break
            end
        end
    end

    finalize_observables!(g)
    return g
end

"""
    stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer, t_data::Dict, savepath::String)

Customized version of `stochastic_edge_acceptance!(g::AbstractGraph, n_steps::Integer)` that
saves the cluster size distrubtions at `t₀` and `t₁` to a .csv in savepath

Arguments
* `g`       : An instance of type AbstractGraph
* `n_steps` : Number of edges to add to the AbstractGraph
* `t_data`  : Dict with `(N, seed)` as keys and `(t₀, t₁)` as values
* `savepath`: Path which to save the cluster size distributions in
Returns
* `g`, updates `g` in-place
"""
function stochastic_edge_acceptance!(
    g::AbstractGraph,
    n_steps::Integer,
    t_data::Dict,
    savepath::String,
)
    t₀, t₁ = t_data[(g.N, Int(g.rng.seed[1]))]
    for t = 1:n_steps
        edge₁ = choose_edge(g)
        edge₂ = choose_edge(g, edge₁)

        if g.cluster_ids[edge₁[1]] == g.cluster_ids[edge₁[2]]
            add_edge!(g, edge₁)
        elseif g.cluster_ids[edge₂[1]] == g.cluster_ids[edge₂[2]]
            add_edge!(g, edge₂)
        else
            C₁ =
                length(g.clusters[g.cluster_ids[edge₁[1]]]) +
                length(g.clusters[g.cluster_ids[edge₁[2]]])
            C₂ =
                length(g.clusters[g.cluster_ids[edge₂[1]]]) +
                length(g.clusters[g.cluster_ids[edge₂[2]]])
            p = (1 / C₁) / (1 / C₁ + 1 / C₂)

            if p > rand(g.rng)
                add_edge!(g, edge₁)
            else
                add_edge!(g, edge₂)
            end
        end

        if t == t₀
            savefile = joinpath(
                savepath,
                "t_0",
                "$(typeof(g))_stochastic_edge_acceptance_$(Int(log2(g.N)))_t_0_cluster_size_distribution_seed_$(Int(g.rng.seed[1])).csv",
            )
            save_cluster_size_dict(g, savefile)
        end
        if t == t₁
            savefile = joinpath(
                savepath,
                "t_1",
                "$(typeof(g))_stochastic_edge_acceptance_$(Int(log2(g.N)))_t_1_cluster_size_distribution_seed_$(Int(g.rng.seed[1])).csv",
            )
            save_cluster_size_dict(g, savefile)
        end
    end

    finalize_observables!(g)
    return g
end

"""
    save_cluster_size_dict(d, n, savefile)
Saves the input cluster size distribution dictionary `d` to `savefile`

Arguments
* `d`       : Cluster size distribution dictionary, i.e. `g.cluster_sizes`
* `N`       : Number of nodes in `g`, i.e. `g.N`
* `savefile`: Path which to save the cluster size distributions in
"""
function save_cluster_size_dict(g, savefile)
    open(savefile, "w") do f
        write(f, "# N = $(g.N)\n")
        write(f, "# seed = $(Int(g.rng.seed[1]))\n")
        write(f, "cluster_size,cluster_size_count\n")
        for (key, value) in g.cluster_sizes
            write(f, "$(key),$(value)\n")
        end
    end
end


text:
#
# errors.jl --
#
# Management of errors.
#
#-------------------------------------------------------------------------------
#
# This file is part of "AndorCameras.jl" released under the MIT license.
#
# Copyright (C) 2017-2019, Éric Thiébaut.
#

Base.showerror(io::IO, e::AndorError) =
    print(io, "Error in ", e.func, ": ", geterrormessage(e.code), " [",
          geterrorsymbol(e.code), "]")

const _ERR_REASONS = Dict{AT.STATUS,String}(
    AT.SUCCESS                      => "Function call has been successful",
    AT.ERR_NOTINITIALISED           => "Function called with an uninitialized handle",
    AT.ERR_NOTIMPLEMENTED           => "Feature has not been implemented for the chosen camera",
    AT.ERR_READONLY                 => "Feature is read only",
    AT.ERR_NOTREADABLE              => "Feature is currently not readable",
    AT.ERR_NOTWRITABLE              => "Feature is currently not writable/executable",
    AT.ERR_OUTOFRANGE               => "Value/index is outside the valid range",
    AT.ERR_INDEXNOTAVAILABLE        => "Index is currently not available",
    AT.ERR_INDEXNOTIMPLEMENTED      => "Index is not implemented for the chosen camera",
    AT.ERR_EXCEEDEDMAXSTRINGLENGTH  => "String value provided exceeds the maximum allowed length",
    AT.ERR_CONNECTION               => "Error connecting to or disconnecting from hardware",
    AT.ERR_NODATA                   => "No Internal Event or Internal Error",
    AT.ERR_INVALIDHANDLE            => "Invalid device handle passed to function ",
    AT.ERR_TIMEDOUT                 => "The function timed out while waiting for data arrive in output queue",
    AT.ERR_BUFFERFULL               => "The input queue has reached its capacity",
    AT.ERR_INVALIDSIZE              => "The size of a queued buffer did not match the frame size",
    AT.ERR_INVALIDALIGNMENT         => "A queued buffer was not aligned on an 8-byte boundary",
    AT.ERR_COMM                     => "An error has occurred while communicating with hardware",
    AT.ERR_STRINGNOTAVAILABLE       => "Index/string is not available",
    AT.ERR_STRINGNOTIMPLEMENTED     => "Index/string is not implemented for the chosen camera",
    AT.ERR_NULL_FEATURE             => "NULL feature name passed to function",
    AT.ERR_NULL_HANDLE              => "Null device handle passed to function",
    AT.ERR_NULL_IMPLEMENTED_VAR     => "Feature not implemented",
    AT.ERR_NULL_READABLE_VAR        => "Readable not set",
    AT.ERR_NULL_WRITABLE_VAR        => "Writable not set",
    AT.ERR_NULL_MINVALUE            => "NULL min value",
    AT.ERR_NULL_MAXVALUE            => "NULL max value",
    AT.ERR_NULL_VALUE               => "NULL value returned from function",
    AT.ERR_NULL_STRING              => "NULL string returned from function",
    AT.ERR_NULL_COUNT_VAR           => "NULL feature count",
    AT.ERR_NULL_ISAVAILABLE_VAR     => "Available not set",
    AT.ERR_NULL_MAXSTRINGLENGTH     => "Max string length is NULL",
    AT.ERR_NULL_EVCALLBACK          => "Event callback parameter is NULL",
    AT.ERR_NULL_QUEUE_PTR           => "Pointer to queue is NULL",
    AT.ERR_NULL_WAIT_PTR            => "Wait pointer is NULL",
    AT.ERR_NULL_PTRSIZE             => "Pointer size is NULL",
    AT.ERR_NOMEMORY                 => "No memory has been allocated for the current action",
    AT.ERR_DEVICEINUSE              => "Function failed to connect to a device because it is already being used",
    AT.ERR_HARDWARE_OVERFLOW        => "The software was not able to retrieve data from the card or camera fast enough to avoid the internal hardware buffer bursting")

geterrormessage(code::Integer) =
    get(_ERR_REASONS, code, "Unknown error code")

const _ERR_SYMBOLS = Dict{AT.STATUS,Symbol}(
    AT.SUCCESS                      => :AT_SUCCESS,
    AT.ERR_NOTINITIALISED           => :AT_ERR_NOTINITIALISED,
    AT.ERR_NOTIMPLEMENTED           => :AT_ERR_NOTIMPLEMENTED,
    AT.ERR_READONLY                 => :AT_ERR_READONLY,
    AT.ERR_NOTREADABLE              => :AT_ERR_NOTREADABLE,
    AT.ERR_NOTWRITABLE              => :AT_ERR_NOTWRITABLE,
    AT.ERR_OUTOFRANGE               => :AT_ERR_OUTOFRANGE,
    AT.ERR_INDEXNOTAVAILABLE        => :AT_ERR_INDEXNOTAVAILABLE,
    AT.ERR_INDEXNOTIMPLEMENTED      => :AT_ERR_INDEXNOTIMPLEMENTED,
    AT.ERR_EXCEEDEDMAXSTRINGLENGTH  => :AT_ERR_EXCEEDEDMAXSTRINGLENGTH,
    AT.ERR_CONNECTION               => :AT_ERR_CONNECTION,
    AT.ERR_NODATA                   => :AT_ERR_NODATA,
    AT.ERR_INVALIDHANDLE            => :AT_ERR_INVALIDHANDLE,
    AT.ERR_TIMEDOUT                 => :AT_ERR_TIMEDOUT,
    AT.ERR_BUFFERFULL               => :AT_ERR_BUFFERFULL,
    AT.ERR_INVALIDSIZE              => :AT_ERR_INVALIDSIZE,
    AT.ERR_INVALIDALIGNMENT         => :AT_ERR_INVALIDALIGNMENT,
    AT.ERR_COMM                     => :AT_ERR_COMM,
    AT.ERR_STRINGNOTAVAILABLE       => :AT_ERR_STRINGNOTAVAILABLE,
    AT.ERR_STRINGNOTIMPLEMENTED     => :AT_ERR_STRINGNOTIMPLEMENTED,
    AT.ERR_NULL_FEATURE             => :AT_ERR_NULL_FEATURE,
    AT.ERR_NULL_HANDLE              => :AT_ERR_NULL_HANDLE,
    AT.ERR_NULL_IMPLEMENTED_VAR     => :AT_ERR_NULL_IMPLEMENTED_VAR,
    AT.ERR_NULL_READABLE_VAR        => :AT_ERR_NULL_READABLE_VAR,
    AT.ERR_NULL_WRITABLE_VAR        => :AT_ERR_NULL_WRITABLE_VAR,
    AT.ERR_NULL_MINVALUE            => :AT_ERR_NULL_MINVALUE,
    AT.ERR_NULL_MAXVALUE            => :AT_ERR_NULL_MAXVALUE,
    AT.ERR_NULL_VALUE               => :AT_ERR_NULL_VALUE,
    AT.ERR_NULL_STRING              => :AT_ERR_NULL_STRING,
    AT.ERR_NULL_COUNT_VAR           => :AT_ERR_NULL_COUNT_VAR,
    AT.ERR_NULL_ISAVAILABLE_VAR     => :AT_ERR_NULL_ISAVAILABLE_VAR,
    AT.ERR_NULL_MAXSTRINGLENGTH     => :AT_ERR_NULL_MAXSTRINGLENGTH,
    AT.ERR_NULL_EVCALLBACK          => :AT_ERR_NULL_EVCALLBACK,
    AT.ERR_NULL_QUEUE_PTR           => :AT_ERR_NULL_QUEUE_PTR,
    AT.ERR_NULL_WAIT_PTR            => :AT_ERR_NULL_WAIT_PTR,
    AT.ERR_NULL_PTRSIZE             => :AT_ERR_NULL_PTRSIZE,
    AT.ERR_NOMEMORY                 => :AT_ERR_NOMEMORY,
    AT.ERR_DEVICEINUSE              => :AT_ERR_DEVICEINUSE,
    AT.ERR_HARDWARE_OVERFLOW        => :AT_ERR_HARDWARE_OVERFLOW)

geterrorsymbol(code::Integer) =
    get(_ERR_SYMBOLS, code, :AT_ERR_UNKNOWN)


text:
"""
    secant(f,p0,p1)

adalah fungsi untuk menyelesaikan `f(x)=0` menggunakan metode secant. Metode ini
memerlukan dua nilai tebakan awal `p0` dan `p1`.

# Examples
```jl
julia> f(x) = x^3-3*x+2;

julia> pk, flag, M = secant(f,-2.6,-2.4);

julia> pk
-2.0

julia> flag
0

julia> M
9×3 Array{Float64,2}:
 0.0  -2.6      NaN
 1.0  -2.4      NaN
 2.0  -2.1066     0.293401
 3.0  -2.02264    0.0839576
 4.0  -2.00151    0.0211303
 5.0  -2.00002    0.00148856
 6.0  -2.0        2.25138e-5
 7.0  -2.0        2.26855e-8
 8.0  -2.0        3.40616e-13
```
return solusi `p0` dengan `flag` bernilai 0 jika metode bisection berhasil menemukan
solusi dan gagal jika tidak 0. Serta, matriks `M` yang berisi catatan proses tiap
iterasi `[k, pk, error]`
"""
function secant(f,p0,p1)
    # Definisikan nilai toleransi, maksimum iterasi dan tebakan awal.
    delta = 10^-12;
    maxi = 100;
    M = [0 p0 NaN
         1 p1 NaN];
    p = [p0 p1];
    pk = p1;
    flag = 1;
    # Mulai langkah iterasi
    for k = 2:maxi
        # rumus metode secant
        pk=p[k]-f(p[k])*(p[k]-p[k-1])/(f(p[k])-f(p[k-1]));
        p = [p pk]
        # Hitung nilai galat mutlak dan relatif
        err=abs(p[k+1]-p[k]);
        rel=2*err/(abs(p[k]) + eps());
        M = [M; [k pk err]]
        # Kriteria penghentian iterasi jika galat memenuhi toleransi.
        if err<delta || rel<delta
            flag = 0;
            break
        end
    end
    return pk, flag, M
end


text:
# Trust-region
# ------------------------------------------------------------
@with_kw mutable struct TRData{T <: AbstractFloat}
    Q̃::T = 1e10
    Δ::T = 1.0
    cΔ::Int = 0
    incumbent::Int = 1
    major_iterations::Int = 0
    minor_iterations::Int = 0
end

@with_kw mutable struct TRParameters{T <: AbstractFloat}
    γ::T = 1e-4
    Δ::T = 1.0
    Δ̅::T = 1000.0
end

"""
    TrustRegion

Functor object for using trust-region regularization in an L-shaped algorithm. Create by supplying a [`TR`](@ref) object through `regularize ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

...
# Parameters
- `γ::T = 1e-4`: Relative tolerance for deciding if a minor iterate should be accepted as a new major iterate.
- `Δ::AbstractFloat = 1.0`: Initial size of ∞-norm trust-region.
- `Δ̅::AbstractFloat = 1000.0`: Maximum size of ∞-norm trust-region.
...
"""
struct TrustRegion{T <: AbstractFloat, A <: AbstractVector} <: AbstractRegularization
    data::TRData{T}
    parameters::TRParameters{T}

    ξ::A
    Q̃_history::A
    Δ_history::A
    incumbents::Vector{Int}

    function TrustRegion(ξ₀::AbstractVector; kw...)
        T = promote_type(eltype(ξ₀), Float32)
        ξ₀_ = convert(AbstractVector{T}, copy(ξ₀))
        A = typeof(ξ₀_)
        return new{T, A}(TRData{T}(), TRParameters{T}(;kw...), ξ₀_, A(), A(), Vector{Int}())
    end
end

function initialize_regularization!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    tr.data.Δ = tr.parameters.Δ
    set_trustregion!(lshaped, tr)
    return nothing
end

function log_regularization!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    @unpack Q̃,Δ,incumbent = tr.data
    push!(tr.Q̃_history, Q̃)
    push!(tr.Δ_history, Δ)
    push!(tr.incumbents, incumbent)
    return nothing
end

function log_regularization!(lshaped::AbstractLShapedSolver, t::Integer, tr::TrustRegion)
    @unpack Q̃,Δ,incumbent = tr.data
    tr.Q̃_history[t] = Q̃
    tr.Δ_history[t] = Δ
    tr.incumbents[t] = incumbent
    return nothing
end

function take_step!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack τ = lshaped.parameters
    @unpack Q̃ = tr.data
    @unpack γ = tr.parameters
    need_update = false
    t = timestamp(lshaped)
    Q̃t = incumbent_objective(lshaped, t, tr)
    if Q + τ <= Q̃ && (tr.data.major_iterations == 1 || Q <= Q̃t - γ*abs(Q̃t-θ))
        need_update = true
        enlarge_trustregion!(lshaped, tr)
        tr.data.cΔ = 0
        tr.ξ .= current_decision(lshaped)
        tr.data.Q̃ = Q
        tr.data.incumbent = timestamp(lshaped)
        tr.data.major_iterations += 1
    else
        need_update = reduce_trustregion!(lshaped, tr)
        tr.data.minor_iterations += 1
    end
    if need_update
        set_trustregion!(lshaped, tr)
    end
    return nothing
end

function process_cut!(lshaped::AbstractLShapedSolver, cut::HyperPlane{FeasibilityCut}, tr::TrustRegion)
    @unpack τ = lshaped.parameters
    if !satisfied(cut, tr.ξ, τ)
        A = [I cut.δQ; cut.δQ' 0*I]
        b = [zeros(length(tr.ξ)); -gap(cut, tr.ξ)]
        t = A\b
        tr.ξ .= tr.ξ + t[1:length(tr.ξ)]
        set_trustregion!(lshaped, tr)
    end
    return nothing
end

function set_trustregion!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    @unpack Δ = tr.data
    nt = nthetas(lshaped)
    l = max.(StochasticPrograms.get_stage_one(lshaped.stochasticprogram).colLower, tr.ξ .- Δ)
    append!(l, fill(-Inf,nt))
    u = min.(StochasticPrograms.get_stage_one(lshaped.stochasticprogram).colUpper, tr.ξ .+ Δ)
    append!(u, fill(Inf,nt))
    MPB.setvarLB!(lshaped.mastersolver.lqmodel, l)
    MPB.setvarUB!(lshaped.mastersolver.lqmodel, u)
    return nothing
end

function enlarge_trustregion!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack τ, = lshaped.parameters
    @unpack Δ = tr.data
    @unpack Δ̅ = tr.parameters
    t = timestamp(lshaped)
    Δ̃ = incumbent_trustregion(lshaped, t, tr)
    ξ = incumbent_decision(lshaped, t, tr)
    Q̃ = incumbent_objective(lshaped, t, tr)
    if Q̃ - Q >= 0.5*(Q̃-θ) && abs(norm(ξ-lshaped.x,Inf) - Δ̃) <= τ
        # Enlarge the trust-region radius
        tr.data.Δ = max(Δ, min(Δ̅, 2*Δ̃))
        return true
    else
        return false
    end
end

function reduce_trustregion!(lshaped::AbstractLShapedSolver, tr::TrustRegion)
    @unpack Q,θ = lshaped.data
    @unpack Q̃,Δ,cΔ = tr.data
    t = timestamp(lshaped)
    Δ̃ = incumbent_trustregion(lshaped, t, tr)
    Q̃ = incumbent_objective(lshaped, t, tr)
    ρ = min(1, Δ̃)*(Q-Q̃)/(Q̃-θ)
    if ρ > 0
        tr.data.cΔ += 1
    end
    if ρ > 3 || (cΔ >= 3 && 1 < ρ <= 3)
        # Reduce the trust-region radius
        tr.data.cΔ = 0
        tr.data.Δ = min(Δ, (1/min(ρ,4))*Δ̃)
        return true
    else
        return false
    end
end

# API
# ------------------------------------------------------------
"""
    TR

Factory object for [`TrustRegion`](@ref). Pass to `regularize ` in the `LShapedSolver` factory function. Equivalent factory calls: `TR`, `WithTR`, `TrustRegion`, `WithTrustRegion`. See ?TrustRegion for parameter descriptions.

"""
struct TR <: AbstractRegularizer
    parameters::Dict{Symbol,Any}
end
TR(; kw...) = TR(Dict{Symbol,Any}(kw))
WithTR(; kw...) = TR(Dict{Symbol,Any}(kw))
TrustRegion(; kw...) = TR(Dict{Symbol,Any}(kw))
WithTrustRegion(; kw...) = TR(Dict{Symbol,Any}(kw))

function (tr::TR)(x::AbstractVector)
    return TrustRegion(x; tr.parameters...)
end

function str(::TR)
    return "L-shaped using trust-region"
end


text:
using HDF5

out = ARGS[3]

sm = h5read(ARGS[1], "/eigenvalues")
lg = h5read(ARGS[2], "/eigenvalues")
egs = vcat(sm, lg)
min, max = extrema(egs)
open(out, "w") do io
    write(io,"condition_number,smallest,largest\n$(max/min),$min,$max")
end


text:
struct SecondaryLimbDark{LD1<:AbstractLimbDark,LD2<:AbstractLimbDark,T} <: AbstractLimbDark
    primary_driver::LD1
    secondary_driver::LD2
    brightness_ratio::T
end

@doc raw"""
    SecondaryLimbDark(primary::AbstractLimbDark,
                      secondary::AbstractLimbDark;
                      brightness_ratio=1)
    SecondaryLimbDark(u_p::AbstractVector, u_s=u_p; kwargs...)

Compose two limb darkening laws together to add a secondary eclipse. If vectors of coefficients are provided, laws will automatically be constructed using [`PolynomialLimbDark`](@ref). The surface brightness ratio is given in terms of the host; e.g., if the companion is half as bright as the host, the ratio would be 0.5.

!!! note "Interface"
    [`SecondaryLimbDark`](@ref) only works with an orbit, since the companion's reference frame needs to be calculated. This means you can't call it using an impact parameter like `ld(b, r)` directly.

# Mathematical form
```math
f(t, r) = \frac{2f_p(t, r) + \eta r^2 f_s(t', r')}{1 + f_p(t, r) + \eta r^2 f_s(t', r')}
```
where $f_p$ is to the primary flux, $f_s$ is to the secondary flux, and $\eta$ is the surface brightness ratio. $t'$ and $r'$ correspond to the time and radius ratio from the companion's reference frame.

# Examples
```jldoctest second
# equal size and limb darkening
r = 1.0
u = [0.4, 0.26]
# companion is 1/10 as bright
brightness_ratio = 0.1
ld = SecondaryLimbDark(u; brightness_ratio)
orbit = SimpleOrbit(period=2, duration=0.5)
fp = ld(orbit, 0, r) # primary egress
fs = ld(orbit, 1, r) # secondary egress

fp ≈ brightness_ratio * fs

# output
true
```
"""
function SecondaryLimbDark(driver1::AbstractLimbDark, driver2::AbstractLimbDark; brightness_ratio=1)
    return SecondaryLimbDark(driver1, driver2, brightness_ratio)
end

function SecondaryLimbDark(u1::AbstractVector, u2=u1; kwargs...)
    driver1 = PolynomialLimbDark(u1)
    driver2 = PolynomialLimbDark(u2)
    return SecondaryLimbDark(driver1, driver2; kwargs...)
end

function compute(ld::SecondaryLimbDark, orbit::AbstractOrbit, t, r; kwargs...)
    f1 = ld.primary_driver(orbit, t, r; kwargs...)
    orbit2 = Orbits.flip(orbit, r * orbit.r_star)
    f2 = ld.secondary_driver(orbit2, t, 1/r; kwargs...)
    flux_ratio = ld.brightness_ratio * r^2
    return (f1 + flux_ratio * f2) / (1 + flux_ratio)
end