text:
#*******************************************************************************
#
#
#                Aplicação do split-operator em um sistema 1D
#
#
#*******************************************************************************

include("save_the_data.jl")

# constantes físicas
Ry = 13.6056917253e3; # meV
a0 = 0.0529177208319; # nm
h2m = (Ry*a0^2)/0.067; # massa efetiva GaAs
hbar = 0.65821188926; # meV.ps

# definição da rede espacial
L = 1000;
Nx = 2*L;
x_grid = linspace(-L/2, L/2, Nx+1)
x_grid = x_grid[1:(end-1)];
Δx = x_grid[2] - x_grid[1];

# definição da rede do "espaço-recíproco"
k = fftshift(linspace(-pi/Δx, pi/Δx, Nx+1)[1:end-1]);
Δk = k[2]-k[1];

# Condição inicial
σ = 15 ;
k₀_max = pi/(2*Δx);
frac_max = 0.01; # fração do valor máximo de "k" que será usado na cond. inicial
k₀ = frac_max * k₀_max;  # o valor de "k₀" não pode ultrapassar π/(2*Δx)
x₀ = -200;
Ψ₀ = exp(-(x_grid-x₀).^2/(2*σ^2) + 1im*k₀*x_grid); # Gaussian packet

# Hamiltoniano separado em "T + V"
T_k = h2m * k.^2;

α, β, γ = 10, 0.07, 0;
V_x = α./(e.^(-β*(x_grid-γ))+1);

function psi_dt(ψ,Δt)
  # espaço de posições 'x'
  psi = ψ;
  psi = [exp(-1im*V_x[i]*Δt/(2*hbar))*psi[i] for i=1:Nx ];
  # espaço de momentos 'k'
  fft!(psi);
  psi = [ exp(-1im*T_k[i]*Δt/hbar)*psi[i] for i=1:Nx ];
  # espaço de posições 'x'
  ifft!(psi);
  psi = [exp(-1im*V_x[i]*Δt/(2*hbar))*psi[i] for i=1:Nx ];
  return psi;
end

function psit(Ψ; t_final = 10.0, N_steps = 100)
  tspan = linspace(0, t_final, N_steps);
  delta_t = tspan[2]-tspan[1];

  for i = 1:N_steps
    psi_t = psi_dt(Ψ,delta_t);
    Ψ = psi_t;
  end
  Ψ
end

t_f = 2.5;
ψ_1 = psit(Ψ_0, t_final = t_f, N_steps = 10);
nome_dados = string("psi_t_",round(t_f,2))

nome_arquivo = "dados_split_operator.h5";
nome_grupo = "split_1D";
nome_subg = "potencial_degrau";

atr_cond_ini = Dict("largura_pacote" => round(σ,1),
                        "x_0" => round(x₀,1),
                        "frac_k_max" => round(frac_max,3));

atr_pot = Dict("altura" => round(α,1),
                        "beta" => round(β,2),
                        "posição" => round(γ,1));

save_them_all( nome_arquivo,
               nome_grupo,
               nome_subg,
               ψ_1,
               nome_dados,
               atr_cond_ini,
               atr_pot);

max_stars_repo_name: mhlmedeiros/Split-Operator

max_issues_repo_name: mhlmedeiros/Split-Operator

max_forks_repo_name: mhlmedeiros/Split-Operator


text:
using LinQu, LinearAlgebra, Test, Random

include("extension_test.jl")
include("interface_test.jl")
include("qstates_test.jl")
include("toffoli_test.jl")
include("block_test.jl")

ext: jl


text:
# Winston:
using Winston
# optionally call figure prior to plotting to set the size
figure(width=600, height=400)
# plot some data
pl = plot(cumsum(rand(500) - 0.5), "r", cumsum(rand(500) - 0.5), "b")
# display the plot (not done automatically!)
display(pl)
println("Press enter to continue: ")
readline(STDIN)
# save the current figure
savefig("winston.svg")

max_stars_repo_path: Module 1/Chapter10/Old/winston.jl

max_issues_repo_path: Module 1/Chapter10/Old/winston.jl

max_forks_repo_path: Module 1/Chapter10/Old/winston.jl


text:
## converting pointers to an appropriate unsigned ##

const C_NULL = box(Ptr{Void}, 0)

# pointer to integer
convert(::Type{Uint}, x::Ptr) = box(Uint, unbox(Ptr,x))
convert{T<:Integer}(::Type{T}, x::Ptr) = convert(T,unsigned(x))

# integer to pointer
convert{T}(::Type{Ptr{T}}, x::Integer) = box(Ptr{T},unbox(Uint,uint(x)))

# pointer to pointer
convert{T}(::Type{Ptr{T}}, p::Ptr{T}) = p
convert{T}(::Type{Ptr{T}}, p::Ptr) = box(Ptr{T}, unbox(Ptr,p))

# object to pointer
convert(::Type{Ptr{Uint8}}, x::Symbol) = ccall(:jl_symbol_name, Ptr{Uint8}, (Any,), x)
convert(::Type{Ptr{Int8}}, x::Symbol) = ccall(:jl_symbol_name, Ptr{Int8}, (Any,), x)
convert(::Type{Ptr{Uint8}}, s::ByteString) = convert(Ptr{Uint8}, s.data)
convert(::Type{Ptr{Int8}}, s::ByteString) = convert(Ptr{Int8}, s.data)

convert{T}(::Type{Ptr{T}}, a::Array{T}) = ccall(:jl_array_ptr, Ptr{T}, (Any,), a)
convert(::Type{Ptr{None}}, a::Array) = ccall(:jl_array_ptr, Ptr{None}, (Any,), a)

pointer{T}(::Type{T}, x::Uint) = convert(Ptr{T}, x)
pointer{T}(::Type{T}, x::Ptr) = convert(Ptr{T}, x)
# note: these definitions don't mean any AbstractArray is convertible to
# pointer. they just map the array element type to the pointer type for
# convenience in cases that work.
pointer{T}(x::AbstractArray{T}) = convert(Ptr{T},x)
pointer{T}(x::AbstractArray{T}, i::Integer) = convert(Ptr{T},x)+(i-1)*sizeof(T)

# unsafe pointer to array conversions
pointer_to_array(p, d::Integer, own=false) = pointer_to_array(p, (d,), own)
function pointer_to_array{T,N}(p::Ptr{T}, dims::NTuple{N,Int}, own::Bool=false)
    ccall(:jl_ptr_to_array, Array{T,N}, (Any, Ptr{T}, Any, Int32),
          Array{T,N}, p, dims, own)
end
function pointer_to_array{T,N}(p::Ptr{T}, dims::NTuple{N,Integer}, own::Bool=false)
    for d in dims
        if !(0 <= d <= typemax(Int))
            error("invalid Array dimensions")
        end
    end
    pointer_to_array(p, convert((Int...), dims), own)
end
unsafe_load(p::Ptr,i::Integer) = pointerref(p, int(i))
unsafe_load(p::Ptr) = unsafe_load(p, 1)
unsafe_store!(p::Ptr{Any}, x::ANY, i::Integer) = pointerset(p, x, int(i))
unsafe_store!{T}(p::Ptr{T}, x, i::Integer) = pointerset(p, convert(T,x), int(i))
unsafe_store!{T}(p::Ptr{T}, x) = pointerset(p, convert(T,x), 1)

# convert a raw Ptr to an object reference, and vice-versa
unsafe_pointer_to_objref(p::Ptr) = pointertoref(unbox(Ptr{Void},p))
pointer_from_objref(x::Any) = ccall(:jl_value_ptr, Ptr{Void}, (Any,), x)

integer(x::Ptr) = convert(Uint, x)
unsigned(x::Ptr) = convert(Uint, x)

@eval sizeof{T<:Ptr}(::Type{T}) = $(div(WORD_SIZE,8))
eltype{T}(::Ptr{T}) = T

## limited pointer arithmetic & comparison ##

==(x::Ptr, y::Ptr) = uint(x) == uint(y)
-(x::Ptr, y::Ptr) = uint(x) - uint(y)

+(x::Ptr, y::Integer) = oftype(x, uint(uint(x) + y))
-(x::Ptr, y::Integer) = oftype(x, uint(uint(x) - y))
+(x::Integer, y::Ptr) = y + x

zero{T}(::Type{Ptr{T}}) = convert(Ptr{T}, 0)
zero{T}(x::Ptr{T}) = convert(Ptr{T}, 0)
one{T}(::Type{Ptr{T}}) = convert(Ptr{T}, 1)
one{T}(x::Ptr{T}) = convert(Ptr{T}, 1)

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
################################################################################
# Drone Example
################################################################################
# using Algames
using StaticArrays
using LinearAlgebra
using MeshCat

function quadrotor_experiment(vis::Visualizer; p::Int=2, animate::Bool=false, S::Int=1)
	T = Float64
	# Define the dynamics of the system
	model  = QuadrotorGame(p=p) # game with 3 players with double integrator dynamics in 3D

	# Define the horizon of the problem
	N = 20 # N time steps
	dt = 0.40 # each step lasts 0.1 second
	probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

	# Define the objective of each player
	# We use a LQR cost
	Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3,0.3,1,1,1]; 5*[1,0.3,0.3,1,1,1]])) for i=1:p] # Quadratic state cost
	R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
	# Desrired state
	xf = [SVector{12,T}([1.0,-0.5,+0.3,0,0,0, 0.3,0,0, 0,0,0]),
		  SVector{12,T}([1.5,-0.7,-0.3,0,0,0, 0.3,0,0, 0,0,0]),
		  SVector{12,T}([0.8, 0.1,+0.5,0,0,0, 0.3,0,0, 0,0,0]),
	      SVector{12,T}([0.8,-0.1,-0.5,0,0,0, 0.3,0,0, 0,0,0]),
	      ]
	xf = xf[1:p]
	# Desired control
	uf = [- model.mass * model.gravity[end] / 4 / model.kf * ones(SVector{model.mi[i],T}) for i=1:p]
	# Objectives of the game
	game_obj = GameObjective(Q,R,xf,uf,N,model)
	radius = 0.5*ones(p)
	μ = 20.0*ones(p)
	add_collision_cost!(game_obj, radius, μ)

	# Define the constraints that each player must respect
	game_con = GameConstraintValues(probsize)
	# Add collision avoidance
	radius = 0.08
	add_spherical_collision_avoidance!(game_con, radius)
	# Add wall constraint
	room_walls = [
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00,  1.00, -1.00], [0.00,  0.00, -1.00]),
	    Wall3D([-3.00, -1.00,  1.00], [3.00, -1.00,  1.00], [3.00,  1.00,  1.00], [0.00,  0.00,  1.00]),
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
	    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
	    ]
	door_cylinders = [
		CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00, -1.00, -1.00], :y, 2.00, 0.95),
		CylinderWall([0.00, -1.00,  1.00], :y, 2.00, 0.95),
		]
	add_wall_constraint!(game_con, room_walls)
	add_wall_constraint!(game_con, door_cylinders)

	build_wall!(vis, room_walls, α=0.1, name=:Room)
	build_cylinder!(vis, door_cylinders, radius=radius, α=0.1, name=:Door1)

	# Define the initial state of the system
	x0 = [
		[-1.0, -0.4, -0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		[-0.8,  0.4,  0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		[-1.4,  0.2,  0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		[-1.6, -0.3, -0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
		]
	x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

	# Define the Options of the solver
	opts = Options()
	# Define the game problem
	prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

	# Solve the problem
	t_elapsed = @elapsed for i = 1:S
		newton_solve!(prob)
	end
	t_elapsed /= S

	# Visualize the Results
	plot!(prob.model, prob.pdtraj.pr)
	plot!(prob.stats)

	build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)
	build_xf!(vis, model, xf, α=1.0, name=:Xf)

	animate && (anim = visualize_robot!(vis, model, prob.pdtraj.pr))
	return prob, t_elapsed
end

function double_integrator_3D_experiment(vis::Visualizer; p::Int=2, animate::Bool=false, S::Int=1)
	T = Float64
	# Define the dynamics of the system
	model  = DoubleIntegratorGame(p=p, d=3) # game with 3 players with double integrator dynamics in 3D

	# Define the horizon of the problem
	N = 20 # N time steps
	dt = 0.40 # each step lasts 0.1 second
	probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

	# Define the objective of each player
	# We use a LQR cost
	Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3,0.3]; 5*[1,0.3,0.3]])) for i=1:p] # Quadratic state cost
	R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
	# Desrired state
	xf = [SVector{6,T}([1.0,-0.5,+0.3, 0.3,0,0]),
		  SVector{6,T}([1.5,-0.7,-0.3, 0.3,0,0]),
		  SVector{6,T}([0.8, 0.1,+0.5, 0.3,0,0]),
	      SVector{6,T}([0.8,-0.1,-0.5, 0.3,0,0]),
	      ]
	xf = xf[1:p]
	# Desired control
	uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
	# Objectives of the game
	game_obj = GameObjective(Q,R,xf,uf,N,model)
	radius = 0.5*ones(p)
	μ = 20.0*ones(p)
	add_collision_cost!(game_obj, radius, μ)

	# Define the constraints that each player must respect
	game_con = GameConstraintValues(probsize)
	# Add collision avoidance
	radius = 0.08
	add_spherical_collision_avoidance!(game_con, radius)
	# Add wall constraint
	room_walls = [
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00,  1.00, -1.00], [0.00,  0.00, -1.00]),
	    Wall3D([-3.00, -1.00,  1.00], [3.00, -1.00,  1.00], [3.00,  1.00,  1.00], [0.00,  0.00,  1.00]),
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
	    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
	    ]
	door_cylinders = [
		CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00, -1.00, -1.00], :y, 2.00, 0.95),
		CylinderWall([0.00, -1.00,  1.00], :y, 2.00, 0.95),
		]
	add_wall_constraint!(game_con, room_walls)
	add_wall_constraint!(game_con, door_cylinders)

	build_wall!(vis, room_walls, α=0.1, name=:Room)
	build_cylinder!(vis, door_cylinders, radius=radius, α=0.1, name=:Door1)

	# Define the initial state of the system
	x0 = [
		[-1.0, -0.4, -0.6, 0.0, 0.0, 0.0],
		[-0.8,  0.4,  0.6, 0.0, 0.0, 0.0],
		[-1.4,  0.2,  0.3, 0.0, 0.0, 0.0],
		[-1.6, -0.3, -0.3, 0.0, 0.0, 0.0],
		]
	x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

	# Define the Options of the solver
	opts = Options()
	# Define the game problem
	prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

	# Solve the problem
	t_elapsed = @elapsed for i = 1:S
		newton_solve!(prob)
	end
	t_elapsed /= S

	# Visualize the Results
	plot!(prob.model, prob.pdtraj.pr)
	plot!(prob.stats)

	build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)
	build_xf!(vis, model, xf, α=1.0, name=:Xf)

	animate && (anim = visualize_robot!(vis, model, prob.pdtraj.pr))
	return prob, t_elapsed
end

function double_integrator_2D_experiment(vis::Visualizer; p::Int=2, animate::Bool=false, S::Int=1)
	T = Float64
	# Define the dynamics of the system
	model  = DoubleIntegratorGame(p=p, d=2) # game with 3 players with double integrator dynamics in 3D

	# Define the horizon of the problem
	N = 20 # N time steps
	dt = 0.40 # each step lasts 0.1 second
	probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

	# Define the objective of each player
	# We use a LQR cost
	Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3]; 5*[1,0.3]])) for i=1:p] # Quadratic state cost
	R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
	# Desrired state
	xf = [SVector{4,T}([1.0,-0.5, 0.3,0]),
		  SVector{4,T}([1.5,-0.7, 0.3,0]),
		  SVector{4,T}([0.8, 0.1, 0.3,0]),
	      SVector{4,T}([0.8,-0.1, 0.3,0]),
	      ]
	xf = xf[1:p]
	# Desired control
	uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
	# Objectives of the game
	game_obj = GameObjective(Q,R,xf,uf,N,model)
	radius = 0.5*ones(p)
	μ = 20.0*ones(p)
	add_collision_cost!(game_obj, radius, μ)

	# Define the constraints that each player must respect
	game_con = GameConstraintValues(probsize)
	# Add collision avoidance
	radius = 0.08
	add_spherical_collision_avoidance!(game_con, radius)
	# Add wall constraint
	room_walls = [
	    Wall([-3.00, -1.00], [3.00, -1.00], [0.00, -1.00]),
	    Wall([-3.00,  1.00], [3.00,  1.00], [0.00,  1.00]),
	    ]
	xc = [ 0.00, 0.00]
	yc = [-1.00, 1.00]
	circle_radius = [0.95, 0.95]

	add_wall_constraint!(game_con, room_walls)
	add_circle_constraint!(game_con, xc, yc, circle_radius)

	room_walls_vis = [
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
	    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
	    ]
	door_cylinders_vis = [
		CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
		]
	build_wall!(vis, room_walls_vis, α=0.1, name=:Room)
	build_cylinder!(vis, door_cylinders_vis, radius=radius, α=0.1, name=:Door1)

	# Define the initial state of the system
	x0 = [
		[-1.0, -0.4, 0.0, 0.0],
		[-0.8,  0.4, 0.0, 0.0],
		[-1.4,  0.2, 0.0, 0.0],
		[-1.6, -0.3, 0.0, 0.0],
		]
	x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

	# Define the Options of the solver
	opts = Options()
	# Define the game problem
	prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

	# Solve the problem
	t_elapsed = @elapsed for i = 1:S
		newton_solve!(prob)
	end
	t_elapsed /= S

	# Visualize the Results
	plot!(prob.model, prob.pdtraj.pr)
	plot!(prob.stats)

	build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)
	build_xf!(vis, model, xf, α=1.0, name=:Xf)

	animate && (anim = visualize_robot!(vis, model, prob.pdtraj.pr))
	return prob, t_elapsed
end

function unicycle_experiment(vis::Visualizer; p::Int=2, animate::Bool=false, S::Int=1)
	T = Float64
	# Define the dynamics of the system
	model = UnicycleGame(p=p) # game with 3 players with double integrator dynamics in 3D

	# Define the horizon of the problem
	N = 20 # N time steps
	dt = 0.40 # each step lasts 0.1 second
	probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

	# Define the objective of each player
	# We use a LQR cost
	Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3]; 5*[0.1,1.0]])) for i=1:p] # Quadratic state cost
	R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
	# Desrired state
	xf = [SVector{4,T}([1.0,-0.5, 0.0, 0.3]),
		  SVector{4,T}([1.5,-0.7, 0.0, 0.3]),
		  SVector{4,T}([0.8, 0.1, 0.0, 0.3]),
	      SVector{4,T}([0.8,-0.1, 0.0, 0.3]),
	      ]
	xf = xf[1:p]
	# Desired control
	uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
	# Objectives of the game
	game_obj = GameObjective(Q,R,xf,uf,N,model)
	radius = 0.5*ones(p)
	μ = 20.0*ones(p)
	add_collision_cost!(game_obj, radius, μ)

	# Define the constraints that each player must respect
	game_con = GameConstraintValues(probsize)
	# Add collision avoidance
	radius = 0.08
	add_spherical_collision_avoidance!(game_con, radius)
	# Add wall constraint
	room_walls = [
	    Wall([-3.00, -1.00], [3.00, -1.00], [0.00, -1.00]),
	    Wall([-3.00,  1.00], [3.00,  1.00], [0.00,  1.00]),
	    ]
	xc = [ 0.00, 0.00]
	yc = [-1.00, 1.00]
	circle_radius = [0.95, 0.95]

	add_wall_constraint!(game_con, room_walls)
	add_circle_constraint!(game_con, xc, yc, circle_radius)

	room_walls_vis = [
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
	    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
	    ]
	door_cylinders_vis = [
		CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
		]
	build_wall!(vis, room_walls_vis, α=0.1, name=:Room)
	build_cylinder!(vis, door_cylinders_vis, radius=radius, α=0.1, name=:Door1)

	# Define the initial state of the system
	x0 = [
		[-1.0, -0.4, 0.0, 0.0],
		[-0.8,  0.4, 0.0, 0.0],
		[-1.4,  0.2, 0.0, 0.0],
		[-1.6, -0.3, 0.0, 0.0],
		]
	x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

	# Define the Options of the solver
	opts = Options()
	# Define the game problem
	prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

	# Solve the problem
	t_elapsed = @elapsed for i = 1:S
		newton_solve!(prob)
	end
	t_elapsed /= S

	# Visualize the Results
	plot!(prob.model, prob.pdtraj.pr)
	plot!(prob.stats)

	build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)
	build_xf!(vis, model, xf, α=1.0, name=:Xf)

	animate && (anim = visualize_robot!(vis, model, prob.pdtraj.pr))
	return prob, t_elapsed
end

function bicycle_experiment(vis::Visualizer; p::Int=2, animate::Bool=false, S::Int=1)
	T = Float64
	# Define the dynamics of the system
	model = BicycleGame(p=p) # game with 3 players with double integrator dynamics in 3D

	# Define the horizon of the problem
	N = 20 # N time steps
	dt = 0.40 # each step lasts 0.1 second
	probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

	# Define the objective of each player
	# We use a LQR cost
	Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3]; 5*[1.0, 0.1]])) for i=1:p] # Quadratic state cost
	R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
	# Desrired state
	xf = [SVector{4,T}([1.0,-0.5, 0.3, 0.0]),
		  SVector{4,T}([1.5,-0.7, 0.3, 0.0]),
		  SVector{4,T}([0.8, 0.1, 0.3, 0.0]),
	      SVector{4,T}([0.8,-0.1, 0.3, 0.0]),
	      ]
	xf = xf[1:p]
	# Desired control
	uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
	# Objectives of the game
	game_obj = GameObjective(Q,R,xf,uf,N,model)
	radius = 0.5*ones(p)
	μ = 20.0*ones(p)
	add_collision_cost!(game_obj, radius, μ)

	# Define the constraints that each player must respect
	game_con = GameConstraintValues(probsize)
	# Add collision avoidance
	radius = 0.08
	add_spherical_collision_avoidance!(game_con, radius)
	# Add wall constraint
	room_walls = [
	    Wall([-3.00, -1.00], [3.00, -1.00], [0.00, -1.00]),
	    Wall([-3.00,  1.00], [3.00,  1.00], [0.00,  1.00]),
	    ]
	xc = [ 0.00, 0.00]
	yc = [-1.00, 1.00]
	circle_radius = [0.95, 0.95]

	add_wall_constraint!(game_con, room_walls)
	add_circle_constraint!(game_con, xc, yc, circle_radius)

	room_walls_vis = [
	    Wall3D([-3.00, -1.00, -1.00], [3.00, -1.00, -1.00], [3.00, -1.00,  1.00], [0.00, -1.00,  0.00]),
	    Wall3D([-3.00,  1.00, -1.00], [3.00,  1.00, -1.00], [3.00,  1.00,  1.00], [0.00,  1.00,  0.00]),
	    ]
	door_cylinders_vis = [
		CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
		CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
		]
	build_wall!(vis, room_walls_vis, α=0.1, name=:Room)
	build_cylinder!(vis, door_cylinders_vis, radius=radius, α=0.1, name=:Door1)

	# Define the initial state of the system
	x0 = [
		[-1.0, -0.4, 0.0, 0.0],
		[-0.8,  0.4, 0.0, 0.0],
		[-1.4,  0.2, 0.0, 0.0],
		[-1.6, -0.3, 0.0, 0.0],
		]
	x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

	# Define the Options of the solver
	opts = Options()
	# Define the game problem
	prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

	# Solve the problem
	t_elapsed = @elapsed for i = 1:S
		newton_solve!(prob)
	end
	t_elapsed /= S

	# Visualize the Results
	plot!(prob.model, prob.pdtraj.pr)
	plot!(prob.stats)

	build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj)
	build_xf!(vis, model, xf, α=1.0, name=:Xf)

	animate && (anim = visualize_robot!(vis, model, prob.pdtraj.pr))
	return prob, t_elapsed
end

S = 100

# vis1  = Visualizer()
# open(vis1)
prob1 = []
t1 = []
for p = 1:4
	prob_i, t_i = quadrotor_experiment(vis1, p=p, animate=false, S=S)
	push!(prob1, prob_i)
	push!(t1, t_i)
end
@show t1
# t = [0.12, 0.45, 1.34, 3.32]

# vis2  = Visualizer()
# open(vis2)
prob2 = []
t2 = []
for p = 1:4
	prob_i, t_i = double_integrator_3D_experiment(vis2, p=p, animate=false, S=S)
	push!(prob2, prob_i)
	push!(t2, t_i)
end
@show t2
# t = [0.0159, 0.214, 0.9569, 0.9119]

# vis3  = Visualizer()
# open(vis3)
prob3 = []
t3 = []
for p = 1:4
	prob_i, t_i = bicycle_experiment(vis3, p=p, animate=false, S=S)
	push!(prob3, prob_i)
	push!(t3, t_i)
end
@show t3
# t = [0.029, 0.125, 0.26, 0.454]

# vis4  = Visualizer()
# open(vis4)
prob4 = []
t4 = []
for p = 1:4
	prob_i, t_i = unicycle_experiment(vis4, p=p, animate=false, S=S)
	push!(prob4, prob_i)
	push!(t4, t_i)
end
@show t4
# t = [0.029, 0.132, 0.276, 0.745]

# vis5  = Visualizer()
# open(vis5)
prob5 = []
t5 = []
for p = 1:4
	prob_i, t_i = double_integrator_2D_experiment(vis5, p=p, animate=false, S=S)
	push!(prob5, prob_i)
	push!(t5, t_i)
end
@show t5
# t = [0.0095, 0.0695, 0.1640, 0.6938]

@show scn.(t1)
@show scn.(t2)
@show scn.(t3)
@show scn.(t4)
@show scn.(t5)

import Algames.plot!
function plot!(model::AbstractGameModel, traj::Algames.Traj; plt=plot())
    plot!(plt, legend=false, aspect_ratio=:equal)
    N = length(traj)
    for i = 1:model.p
        xi = [Algames.state(traj[k])[model.pz[i][1]] for k=1:N]
        yi = [Algames.state(traj[k])[model.pz[i][2]] for k=1:N]
        plot!(xi, yi, label=false)
        scatter!(xi, yi)
    end
    display(plt)
    return nothing
end

import Algames.plot!

function plot!(stats::Statistics; plt=plot(), lw::T=5.0) where {T}
	plot!(plt,
		size=(500,500),
		layout=(1,1,))
    iter = stats.iter
    dyn = log.(10, max.(1e-10, [stats.dyn_vio[i].max for i=1:iter]))
    con = log.(10, max.(1e-10, [stats.con_vio[i].max for i=1:iter]))
    sta = log.(10, max.(1e-10, [stats.sta_vio[i].max for i=1:iter]))
    opt = log.(10, max.(1e-10, [stats.opt_vio[i].max for i=1:iter]))
	y_min = minimum([dyn; con; sta; opt])
	y_max = maximum([dyn; con; sta; opt])
	# Set up plot
	plot!(plt[1,1],
		legend=:bottomleft,
		xlabel="Outer Loop Iterations",
		ylabel="log(cons. vio.)",
		title="Constraint Violation")
	# Add curves
	plot!(plt[1,1], dyn, linewidth=lw, label="dyn", legend=:bottomleft)
	plot!(plt[1,1], con, linewidth=lw, label="con", legend=:bottomleft)
	plot!(plt[1,1], sta, linewidth=lw, label="sta", legend=:bottomleft)
	plot!(plt[1,1], opt, linewidth=lw, label="opt", legend=:bottomleft)
	# Add rectangles
	plot_epochs!(plt, y_min, y_max, stats.outer_iter)

    display(plt)
    return nothing
end

import Algames.plot_epochs!
function plot_epochs!(plt, y_min::T, y_max::T, epochs::Vector{Int}) where {T}
	rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
	i_start = 1
	i_end = -1
	for k = 1:epochs[end]
		i_end = findlast(x -> x==k, epochs )
		plot!(rectangle(i_end-i_start,y_max-y_min,i_start,y_min), opacity=.1, label=false)
		i_start = i_end + 1
	end
	return nothing
end


text:
function solve_tspd_RL_test(n_nodes, n_groups, n_samples)
    n_runs = 10
    obj_DPS = zeros(n_runs)
    obj_RL = zeros(n_runs)
    for i in 1:n_runs
        x = rand(n_nodes) .* 100
        y = rand(n_nodes) .* 100
        x[1] = rand() # depot
        y[1] = rand() # depot

        # x = [0.8172268241831585, 73, 61, 44, 74, 23, 9, 51, 15, 9, 56]
        # y = [0.6284331187597952, 52, 46, 63, 15, 23, 80, 65, 97, 74, 84]

        result1 = solve_tspd(x, y, 1.0, 0.5; n_groups=n_groups)
        obj_DPS[i] = result1.total_cost

        result2 = solve_tspd_RL(x, y; n_samples=n_samples)
        obj_RL[i] = result2.total_cost
    end
    pct_gap = (mean(obj_RL) - mean(obj_DPS)) / mean(obj_DPS)
    @show mean(obj_DPS), mean(obj_RL), pct_gap
    @assert pct_gap < 0.05
    return mean(obj_DPS), mean(obj_RL), pct_gap
end

function test_RL()
    sizes = [11, 15, 20, 50, 100]
    n_groups = [1, 1, 2, 5, 10]
    n_samples = 10
    results = zeros(length(sizes), 3)
    for i in 1:length(sizes)
        @info "Testing n = $(sizes[i]) / DPS: n_groups = $(n_groups[i]) /  RL: n_samples = $n_samples"
        results[i, :] .= solve_tspd_RL_test(sizes[i], n_groups[i], n_samples)
    end
end


text:
__precompile__()

module MOOSE

using ContMechTensors

using WriteVTK
using JuAFEM

using ForwardDiff
using ForwardDiff: Partials, Dual, value, partials, npartials, setindex, convert
export Dual

using Reexport

const have_petsc = "PETSC_DIR" in keys(ENV)

if have_petsc
    @reexport using MiniPETSc

    import MiniPETSc.assemble!
    import MiniPETSc.solve!
    import MiniPETSc.zero!
    import MiniPETSc.plusEquals!
    import MiniPETSc.zeroRows!
    import MiniPETSc.serializeToZero
else
    import DummyMPI
    MPI = DummyMPI
    export MPI
end

export dofs, processor_id

export Node, Element, Mesh

export buildSquare, readXDAMesh

export Variable, System, addVariable!, addKernel!, addBC!, initialize!

export Solver, solve!, JuliaDenseImplicitSolver, JuliaDenseNonlinearImplicitSolver

if have_petsc
    export PetscImplicitSolver, PetscNonlinearImplicitSolver
end

export Kernel, Diffusion, Convection

export boundaryIDs, DirichletBC

export output, VTKOutput

export invalid_processor_id

# The default value type that will be used by MOOSE.jl
value_type = Float64

# Constructor for making a Dual with a particular value and index
dualVariable{T}(value::T, index::Int64, num_partials::Int64) = Dual(value, Partials(ntuple(n -> n != index ? zero(Float64) : 1.0, num_partials)))

const invalid_processor_id = -1

include("utils/PerfLog.jl")

export PerfLog, startRootLog!, stopRootLog!, startLog, stopLog, clear!

main_perf_log = PerfLog()

export main_perf_log

function finalize()
    println(main_perf_log)

    if !have_petsc
        MPI.Finalize()
    end
end

function __init__()
    if !have_petsc
        MPI.Init()
    end

    atexit(finalize)
end

include("numerics/JuliaSupport.jl")
include("mesh/DofObject.jl")
include("mesh/Node.jl")
include("mesh/Element.jl")
include("mesh/BoundaryInfo.jl")
include("mesh/Mesh.jl")
include("mesh/Generation.jl")
include("mesh/XDAReader.jl")
include("partitioners/Partitioner.jl")
include("partitioners/SimplePartitioner.jl")

if have_petsc
    include("partitioners/MetisPartitioner.jl")
end

include("mesh/MeshInitialization.jl")
include("systems/Variable.jl")
include("kernels/Kernel.jl")
include("bcs/BoundaryCondition.jl")
include("systems/System.jl")
include("solvers/Solver.jl")
include("solvers/Assembly.jl")
include("outputs/Output.jl")

include("solvers/JuliaDenseImplicitSolver.jl")
include("solvers/JuliaDenseNonlinearImplicitSolver.jl")

if have_petsc
    include("solvers/PetscImplicitSolver.jl")
    include("solvers/PetscNonlinearImplicitSolver.jl")
end

include("kernels/Diffusion.jl")
include("kernels/Convection.jl")

include("bcs/DirichletBC.jl")

include("outputs/VTKOutput.jl")

end


text:
###############################################################################
#
# Functions to display time in @error @info, and @warn
#
###############################################################################
"""
A function to covert information to time + information

    terror(info::String)

Display the error info with date and time before the given
- `info` Infomation to display

---
Examples
```julia
# add a time stamp before the message
@error terror("This is an error with time stamp at the beginning!");
```
"""
function terror(info::String)
    return "$(format(now(),"yyyy-mm-dd HH:MM:SS"))\n       $(info)"
end

"""
A function to covert information to time + information

    tinfo(info::String)

Display the error info with date and time before the given
- `info` Infomation to display

---
Examples
```julia
# add a time stamp before the message
@info tinfo("This is an info with time stamp at the beginning!");
```
"""
function tinfo(info::String)
    return "$(format(now(),"yyyy-mm-dd HH:MM:SS"))\n      $(info)"
end

"""
A function to covert information to time + information

    twarn(info::String)

Display the error info with date and time before the given
- `info` Infomation to display

---
Examples
```julia
# add a time stamp before the message
@warn twarn("This is an info with time stamp at the beginning!");
```
"""
function twarn(info::String)
    return "$(format(now(),"yyyy-mm-dd HH:MM:SS"))\n         $(info)"
end


text:
function genData(m::RandomSearch,nCons::Int,nProdTot::Int,nProdCons::Int; fname::String = "" )
	Random.seed!(m.seed)

	# Consumer index into full data
	seli = [ 1+ii*nProdCons:(ii+1)*nProdCons for ii = 0:nCons-1]

	# Generate products
	pid,chars = genProducts(m::Model,seli,nProdTot,nProdCons)

	# Products' positions (same for all consumers)
	pos = fill(collect(0:nProdCons-1),nCons)

	# Gather in table for convenience
	d = (seli = seli,pid = pid,chars = chars)

	# Get reservation value of continuing search
	cfun = getCfun(m.cfun)
	μ  = sum(mean(m.dChars) .* m.β[1:end-m.outDummy]) + mean(m.dV)
	σ2  = sum(var(m.dChars) .* m.β[1:end-m.outDummy] .^2) + var(m.dV)
	zs = calcZs(0.0,Normal(μ,sqrt(σ2)),cfun(m.cs))[1]

	xβv = chars * m.β .+ rand(m.dV,size(chars,1))
	u   = xβv .+ rand(m.dE,size(d.chars,1))
	path,purch = genSearchPaths(m,d,nProdCons,zs,u)

	path = [path[path[:,i].>0,i] for i =1:nCons]

	# Gather again data
	d = (seli = seli,pid = pid,chars = chars,pos = pos,path = path, purch = purch, u=u , zs = zs)

	if !isempty(fname)
		@info "Saved data in bson format to $fname."
		@save fname d
	end

	return d
end

###################################################
function genSearchPaths(m::RandomSearch,d,maxProd,zs,u)
	n = length(d.seli)

	# Same for all consumers
	path = zeros(Int64,maxProd,n)
	purch = fill(0,n)

  @inbounds Threads.@threads for k = 1:Threads.nthreads()
		threadRange = getrange(n)
		locPathi 	= zeros(Int64,maxProd)
		for ii in threadRange # Loop over consumers allocated to specific thread
			# Select consumer
			r = d.seli[ii]
			nProd = length(d.seli[ii])
			# Reset locals
			locPathi .= 0
			# generate search paths (updates locPathi/locSearchi and returns purchi)
			purchi = genSearchPathi!(m,locPathi,zs,u[r])
			# Fill in values
			path[:,ii] = locPathi
			purch[ii] = purchi
		end

	end
	return path,purch
end

###################################################
function genSearchPathi!(m::RandomSearch,path,zs,u)
# NOTE: ASSUMES TOP DOWN (which is same as random if Xs are random)
	# Initial status
	cc_s = 1
	ind_p = 1 # index of maximum utility (initially outside option)
	u_max = u[ind_p] # max utility
	# Generate search paths
	while true
		if cc_s == size(u,1) # Searched all
			purch = findmax(u)[2]
			return purch
		elseif zs > u_max # search
			cc_s  += 1
			path[cc_s] = cc_s
			u_max = if u[cc_s] > u_max
						ind_p = cc_s
						u[cc_s]
					else
						u_max
					end
		else # Purchase and end search
			purch = ind_p
			return purch
			break # end with purchase
		end
	end
end

###################################################
function negLogLik(m::RandomSearch,θ::Vector{T},d,cfun; e = [], v = [] ) where T<:Real

	# Extract Parameters
	β = θ[1:end-1]
	cs = θ[end]

	# Draw random numbers if not provided
	if isempty(e)
		Random.seed!(m.seed+23)
		e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	end
	if isempty(v)
		Random.seed!(m.seed+3598)
		v = rand(m.dV,size(d.chars,1),m.options.nDraws)
	end

	# Form utility
	xβv = d.chars * β .+ v
	u = xβv .+ e

	# Search values
	μ  = sum(mean(m.dChars) .* m.β[1:end-m.outDummy]) + mean(m.dV)
	σ2  = sum(var(m.dChars) .* m.β[1:end-m.outDummy] .^2) + var(m.dV)
	zs = calcZs(0.0,Normal(μ,sqrt(σ2)),cfun(m.cs))[1]

	# Loop through consumers
	logLik = fill(0.0,Threads.nthreads())

	@inbounds Threads.@threads for k = 1:Threads.nthreads()
		locSumLL = 0.0
		locV = zeros(Float64,3,size(e,2))
		for ii in getrange(length(d.seli))
			locV .= 0.0
			fillV_RS!(locV,d.path[ii],d.purch[ii],u,zs,d.seli[ii],m.smo)
			locSumLL += log(mean(1.0 ./ (1.0 .+ sum(locV,dims=1))))
		end
		logLik[k] = locSumLL
	end
	sum_logLik = abs(sum(logLik)) == Inf ? -1e10 : sum(logLik)
	return -sum_logLik
end

##
# ind_p is index of purchased
function fillV_RS!(v,path,ind_p,u,zs,r,smo)
	nSearch = length(path)
	# Continuation, requires u <= reservation value for those where continued
	if nSearch > 0
		for ee = 1:size(u,2)
			# um = u[r[1],ee]
			# v[1,ee] += exp(-smo*(zs - um))
			# for ll in path[1:end-1]
			# 	um = max(um,u[r[ll],ee])
			# 	v[1,ee] += exp(-smo*(zs - um))
			# end
			u0 = u[r[1],ee]
			um = nSearch > 1 ? max(u0,findmax_range(u,r[path[1:end-1]],ee)[1],u0) : u0
			v[1,ee] += exp(-smo*(zs - um))
		end
	end

	# Stopping, requires max util of searched > reservation value
	if nSearch < length(r)-1
		for ee = 1:size(u,2)
			u0 = u[r[1],ee]
			um = nSearch > 0 ? max(u0,findmax_range(u,r[path],ee)[1]) : u0
			v[2,ee] += exp(-smo*(um-zs))
		end
	end

	# Purchase requires being maximum utility
	if nSearch > 0
		for ee = 1:size(u,2)
			v[3,ee] += exp(-smo*(u[r[ind_p],ee] - u[r[1],ee]))
			for nn in path
				v[3,ee] += exp(-smo*(u[r[ind_p],ee] - u[r[nn],ee]))
			end
			v[3,ee] -= 1 # adjustment since also includes purchased above
		end
	end
  	return
end

###################################################
function estimateModel(m::RandomSearch,d; startvals = [] , pars_true = []  )

	# Preliminaries
	cfun = 	getCfun(m.cfun)

	nBeta = size(d.chars,2)

	if isempty(pars_true)
		pars_true = fill(NaN,nBeta+1)
	end

	if isempty(startvals)
		startvals = rand(Uniform(),nBeta+1)
	end

	# Random draws, pre-allocate
	Random.seed!(m.seed)
	e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	v = rand(m.dV,size(d.chars,1),m.options.nDraws)

	############################################################################
	# Otimization

	# Define objective function
	obj_fun(θ,p) = negLogLik(m,θ,d,cfun;v=v,e=e)
	obj = OptimizationFunction(obj_fun,m.options.autodiff)
	prob = OptimizationProblem(obj,startvals)

	# Run optimization
	res = solve(prob,m.options.algo;
				show_trace = true , show_every=m.options.showEvery,
				maxiters = m.options.maxIter,
				x_tol = m.options.x_tol,
				f_tol = m.options.f_tol,
				g_tol = m.options.g_tol)

	pars_opt = res.minimizer
	fval_opt = res.minimum

	# Calculate obj. fun at start and true val
	startVal = obj_fun(startvals,[])
	trueVal = isnan(pars_true[1]) ? NaN : obj_fun(pars_true,[])

	# Extract results
	β = pars_opt[1:end-1]
	cs = pars_opt[end]
	finalVal = fval_opt

	infotxt = "
	####################################
	# Results of Estimation:
	Estimation model:	$(typeof(m))
	Final Objective Value: 		$finalVal
	Start Objective Value:		$startVal
	True Objective Value:		$trueVal
	Parameters:
		beta  = $(round.(β,digits=2))
		beta0 = $(round.(pars_true[1:end-1],digits=2))
		cs = $(round(cs,digits=2)) | cs0 = $(round(pars_true[end],digits=2))
	#################################### "

	@info infotxt

	return β,cs,finalVal

end

######################################################
function calcWelfareDemand(m::RandomSearch,d; e = [] , v = [] )

	# Draw random numbers if not provided
	if isempty(e)
		Random.seed!(m.seed+43)
		e = rand(m.dE,size(d.chars,1),m.options.nDraws)
	end
	if isempty(v)
		Random.seed!(m.seed+323)
		v = rand(m.dV,size(d.chars,1),m.options.nDraws)
	end

	# Derived values
	nCons = length(d.seli)
	maxProd = maximum(vcat([length(d.seli[i]) for i = 1:length(d.seli)]))

	# Form utility
	xβv = d.chars * m.β .+ v
	u   = xβv .+ e

	# Search costs
	cfun = getCfun(m.cfun)

	μ  = sum(mean(m.dChars) .* m.β[1:end-m.outDummy]) + mean(m.dV)
	σ2  = sum(var(m.dChars) .* m.β[1:end-m.outDummy] .^2) + var(m.dV)
	zs = calcZs(0.0,Normal(μ,sqrt(σ2)),cfun(m.cs))[1]

	# Pre allocate
	sumU = zeros(Threads.nthreads())
	sumC = zeros(Threads.nthreads())
	demand = zeros(Int64,maxProd,Threads.nthreads())
	cs = cfun(m.cs) # Note, need to put here as cfun(m.cs) does not work in threaded
	@inbounds Threads.@threads for k = 1:Threads.nthreads()
		threadRange = getrange(length(d.seli))

		locPathi 	= zeros(Int64,maxProd)
		locU = 0.0
		locC = 0.0
		locD = zeros(Int64,maxProd)
		for ii in threadRange
			@views posi = d.pos[ii]
			for dd in 1:size(e,2)
				@views ud 	= u[d.seli[ii],dd]
				locPathi 	.= 0
				purchi 	= genSearchPathi!(m,locPathi,zs,ud)
				locU 	+= ud[purchi]
				for ss = 2:maxProd # NOTE: locPathi first entry is = 0 due to outside option
					if locPathi[ss] == 0
						break
					else
						locC += cs # note: cs = cfun(m.cs) from above
					end
				end
				locD[purchi] += 1
			end
		end
		sumU[k] = locU / size(e,2)
		sumC[k] = locC / size(e,2)
		demand[:,k] .= locD
	end
	return sum(sumU) - sum(sumC), sum(sumU), sum(sumC), sum(demand,dims=2) ./ size(e,2) ./ nCons
end


text:
module DisjunctiveProgramming

using JuMP, IntervalArithmetic, Symbolics, Suppressor

export add_disjunction!, add_proposition!, reformulate_disjunction
export @disjunction, @proposition
export choose!

include("constraint.jl")
include("logic.jl")
include("bounds.jl")
include("utils.jl")
include("bigm.jl")
include("hull.jl")
include("reformulate.jl")
include("macros.jl")

end # module


text:
#Bootstrap confidence band for univariate nonparametric regression
function bootstrapCB(xdata::RealVector, ydata::RealVector;
    xeval::RealVector=xdata, B::Int64=500, reg::Function=locallinear, lb::Real=-Inf, ub::Real=Inf, kernel::Function=gaussiankernel, h::Real=-Inf, alpha::Real=0.05)

    n=length(xdata)
    y_matrix=zeros(B, length(xeval))
    cb=zeros(2, length(xeval))

    xdata, xeval, kernel = boundit(xdata, xeval, kernel, lb, ub)
    if h <= 0
        h = bwreg(xdata, ydata, reg, kernel)
    end

    yhat=reg(xdata, ydata, xeval=xdata, h=h, kernel=kernel)
	#in the wild bootstrap setting, new error terms are generated from a distribution with 0 mean, equal second and third moment
	#usually work well unless the error terms are extremly small, then the cb may not cover prediction.
    e = ydata .- yhat
    coef=[(1-sqrt(5))/2,(1+sqrt(5))/2]

    for b in 1:B
        boot_coef=coef[rand(Categorical([(5+sqrt(5))/10, 1-(5+sqrt(5))/10]), n)]
        boot_e=boot_coef .* e
        boot_y=yhat .+ boot_e
        y_matrix[b, :]=reg(xdata, boot_y, xeval=xeval, h=h, kernel=kernel)
    end
    for i in 1:length(xeval)
        cb[:, i]=quantile(y_matrix[:, i], [alpha/2, 1-alpha/2])
    end
    cb
end

# #Bootstrap confidence band for multivariate nonparametric regression
# function BootstrapCB(B::Int64, xeval::Matrix{Float64},
#   xdata::Matrix{Float64}, ydata::Vector{Float64}, reg::Function=LP0, kernel::Function=GaussianKernel,
# h::Vector{Float64}=BandwidthLSCVReg(xdata,ydata,reg,kernel))

#   (n, p)=size(xdata)
#   (m,p1 )= size(xeval)
#   if p != p1
#     error("xeval should have same dimension as xdata")
#   end

#   y_matrix=zeros(B, m)
#   cb=zeros(2, m)

#   yhat=reg(xdata, xdata, ydata, kernel, h)
#   mhat=reg(xdata, xdata, ydata, kernel, h)
#   e = ydata .- yhat
#   coef=[-1,0,2]

#   for b in 1:B
# #     print(b, "\t")
#     boot_coef=coef[rand(Categorical([1/3, 1/2, 1/6]), n)]
#     boot_e=boot_coef .* e
#     boot_y=mhat .+ boot_e
#     y_matrix[b, :]=reg(xeval, xdata, boot_y,kernel, h)
#   end
#   for i in 1:m
#     cb[:, i]=quantile(y_matrix[:, i], [.975, .025])
#   end
#   cb
# end


text:
@testset "Box" begin

box = Box()

@test box.id == 1
@test GRAYTR.can_intersect(box)

# Create with transformation
box = Box(GRAYTR.scaling(2.0, 1.0, 1.0))
@test GRAYTR.obj_bounds(box) == GRAYTR.BoundingBox(Point3(-1, -1, -1), Point3(1, 1, 1))
@test GRAYTR.world_bounds(box) == GRAYTR.BoundingBox(Point3(-2, -1, -1), Point3(2, 1, 1))

box = Box(2, GRAYTR.scaling(2.0, 1.0, 1.0))
@test GRAYTR.obj_bounds(box) == GRAYTR.BoundingBox(Point3(-1, -1, -1), Point3(1, 1, 1))
@test box.id == 2

# Transform after creation
T = GRAYTR.scaling(1.0, 2.0, 1.0)
box = T(Box())
@test GRAYTR.obj_bounds(box) == GRAYTR.BoundingBox(Point3(-1, -1, -1), Point3(1, 1, 1))
@test GRAYTR.world_bounds(box) == GRAYTR.BoundingBox(Point3(-1, -2, -1), Point3(1, 2, 1))

# intersectP
D = Box()

R = GRAYTR.Ray(Point3(0.5, 0, 2), Vector3(0, 0, -1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.5, 0, -2), Vector3(0, 0, 1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.5, 0, 2), Vector3(0, 0, -1), 10, Inf, 1)
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.5, 0, -2), Vector3(0, 0, 1), 10, Inf, 1)
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 2, 2), Vector3(0, 0, -1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0, 2, -2), Vector3(0, 0, 1))
@test !GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.5, 0, -2), Vector3(0, 0, 1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.5, 0, 2), Vector3(0, 0, -1))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.0, 0, -0.2), Vector3(0, 1, 0))
@test GRAYTR.intersectP(R, D)

R = GRAYTR.Ray(Point3(0.0, 0, 0.0), Vector3(0, -1, 0))
@test GRAYTR.intersectP(R, D)

# intersect

R = GRAYTR.Ray(Point3(0.5, 0, 2), Vector3(0, 0, -1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0.5, 0, 1.0)
@test dg.n ≈ normalize(Normal3(0, 0, 1))
@test t ≈ 1.0

R = GRAYTR.Ray(Point3(0.5, 0, -3), Vector3(0, 0, 1))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0.5, 0, -1.0)
@test dg.n ≈ normalize(Normal3(0.0, 0, -1.0))
@test t ≈ 2.0

R = GRAYTR.Ray(Point3(2, 0, -0.5), Vector3(-1, 0, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(1.0, 0.0, -0.5)
@test dg.n ≈ normalize(Normal3(1.0, 0.0, 0.0))
@test t ≈ 1.0

# Ray along Y axis, hitting box

R = GRAYTR.Ray(Point3(0, 2, -0.5), Vector3(0, -1, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0.0, 1.0, -0.5)
@test dg.n ≈ normalize(Normal3(0.0, 1.0, 0.0))
@test t ≈ 1.0

R = GRAYTR.Ray(Point3(0, -2, -0.5), Vector3(0, 1, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg.p ≈ Point3(0.0, -1.0, -0.5)
@test dg.n ≈ normalize(Normal3(0.0, -1.0, 0.0))
@test t ≈ 1.0

# Ray along Y axis, not hitting box

R = GRAYTR.Ray(Point3(0, 2, -1.5), Vector3(0, -1, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

R = GRAYTR.Ray(Point3(0, -2, -1.5), Vector3(0, 1, 0))
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

# Box not inside ray limits

R = GRAYTR.Ray(Point3(0, 0, 2), Vector3(0, 0, -1), 10, Inf, 1)
dg, t, e = GRAYTR.shape_intersect(R, D)
@test dg == nothing
@test isnan(t)
@test isnan(e)

end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
#=
FloydWarshallTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/FloydWarshall.jl")
matrix =  GBMatrix([[0, 5, 4, 7] [3, 0, 7, 2] [4, 1, 0, 3] [1, 6, 5, 0]])
@test floydwarshall(matrix, 4) == [[8, 13, 12, 7] [3, 8, 7, 2] [4, 9, 8, 3] [1, 6, 5, 0]]


text:
using HDF5
using Base.Test
using Compat

@testset "extendtest" begin

    fn = joinpath(tempdir(),"test.h5")

    fid = h5open(fn, "w")
    g = g_create(fid, "shoe")
    d = d_create(g, "foo", datatype(Float64), ((10,20),(100,200)), "chunk", (1,1))
    #println("d is size current $(map(int,HDF5.get_dims(d)[1])) max $(map(int,HDF5.get_dims(d)[2]))")
    dims, max_dims = HDF5.get_dims(d)
    ui64(n) = UInt64(n)
    @test dims == (ui64(10),ui64(20))
    @test max_dims == (ui64(100),ui64(200))
    set_dims!(d, (100,150))
    dims, max_dims = HDF5.get_dims(d)
    @test dims == (ui64(100),ui64(150))
    @test max_dims == (ui64(100),ui64(200))
    d[1,1:5] = [1.1231,1.313,5.123,2.231,4.1231]
    set_dims!(d, (1,5))
    @test size(d) == (1,5)
    #println("d is size current $(map(int,HDF5.get_dims(d)[1])) max $(map(int,HDF5.get_dims(d)[2]))")
    b = d_create(fid, "b", Int, ((1000,),(-1,)), "chunk", (100,)) #-1 is equivalent to typemax(Hsize) as far as I can tell
    #println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")
    b[1:200] = ones(200)
    dims, max_dims = HDF5.get_dims(b)
    @test dims == (ui64(1000),)
    @test max_dims == (HDF5.MAXIMUM_DIM,)
    set_dims!(b, (10000,))
    dims, max_dims = HDF5.get_dims(b)
    @test dims == (ui64(10000),)
    @test max_dims == (HDF5.MAXIMUM_DIM,)
    #println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")
    # b[:] = [1:10000] # gave error no method endof(HDF5Dataset{PlainHDF5File},),
    # so I defined endof(dset::HDF5Dataset) = length(dset), and exported endof
    # but that didn't fix the error, despite the endof function working
    # d[1] produces error ERROR: Wrong number of indices supplied, should datasets support linear indexing?
    b[1:10000] = [1:10000;]
    #println(b[1:100])

    close(fid)

    fid = h5open(fn, "r")
    d_again = fid["shoe/foo"]
    dims, max_dims = HDF5.get_dims(d_again)
    @test dims==(ui64(1),ui64(5))
    @test max_dims == (ui64(100),ui64(200))
    @test (sum(d_again[1,1:5])-sum([1.1231,1.313,5.123,2.231,4.1231])) == 0
    #println("d is size current $(map(int,HDF5.get_dims(re_d)[1])) max $(map(int,HDF5.get_dims(re_d)[2]))")
    @test fid["b"][1:10000] == [1:10000;]
    b_again = fid["b"]
    dims, max_dims = HDF5.get_dims(b_again)
    @test dims == (ui64(10000),)
    @test max_dims == (HDF5.MAXIMUM_DIM,)
    #println("b is size current $(map(int,HDF5.get_dims(b)[1])) max $(map(int,HDF5.get_dims(b)[2]))")

    close(fid)

end


text:
function ConstantMagneticPotentialDifference(;name, V_m=1.0)
    val = V_m
    @named two_port_elementary = TwoPortElementary()
    @unpack port_p, port_n = two_port_elementary
    @parameters V_m
    @variables Phi(t)
    eqs = [
        V_m ~ port_p.V_m - port_n.V_m,
        Phi ~ port_p.Phi,
        0 ~ port_p.Phi + port_n.Phi,
    ]
    extend(ODESystem(eqs, t, [Phi], [V_m], systems=[port_p, port_n], defaults=Dict(V_m => val), name=name), two_port_elementary)
end

function ConstantMagneticFlux(;name, Phi=1.0)
    val = Phi
    @named two_port_elementary = TwoPortElementary()
    @unpack port_p, port_n = two_port_elementary
    @parameters Phi
    @variables V_m(t)
    eqs = [
        V_m ~ port_p.V_m - port_n.V_m,
        Phi ~ port_p.Phi,
        0 ~ port_p.Phi + port_n.Phi,
    ]
    extend(ODESystem(eqs, t, [V_m], [Phi], systems=[port_p, port_n], defaults=Dict(Phi => val), name=name), two_port_elementary)
end


text:
# Autogenerated wrapper script for libcleri_jll for armv6l-linux-gnueabihf
export libcleri

using PCRE2_jll
JLLWrappers.@generate_wrapper_header("libcleri")
JLLWrappers.@declare_library_product(libcleri, "libcleri.so.0")
function __init__()
    JLLWrappers.@generate_init_header(PCRE2_jll)
    JLLWrappers.@init_library_product(
        libcleri,
        "lib/libcleri.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module ARCTools

using Glob
using Parameters
using ReadableRegex

export Resources
export QSub, qsub
export JuliaCmd, JuliaJob

#### QSub

const Maybe{T} = Union{T, Nothing}

"""
https://confluence.it.ubc.ca/display/UARC/Running+Jobs
"""
@with_kw mutable struct Resources
    walltime::Maybe{String} = nothing
    select::Maybe{Int} = nothing
    ncpus::Maybe{Int} = nothing
    ngpus::Maybe{Int} = nothing
    mpiprocs::Maybe{Int} = nothing
    ompthreads::Maybe{Int} = ncpus
    mem::Maybe{Int} = nothing
    gpu_mem::Maybe{Int} = nothing
end

function resource_string(r::Resources, prop::Symbol)
    val = getfield(r, prop)
    units = (prop === :mem || prop === :gpu_mem) ? "gb" : ""
    sep = (prop === :walltime) ? "," : ":"
    val === nothing ? "" : "$(prop)=$(val)$(units)$(sep)"
end

function resource_list(r::Resources)
    r_list = ""
    for prop in [
            # Note: order matters here for :walltime and :select
            :walltime,
            :select,
            :ncpus,
            :ngpus,
            :mpiprocs,
            :ompthreads,
            :mem,
            :gpu_mem,
        ]
        r_list *= resource_string(r, prop)
    end

    # Return resource list, removing trailing separator
    @assert !isempty(r_list)
    r_list[begin:end-1]
end

#### QSub

@with_kw mutable struct QSub
    resources::Resources
    modules::Vector{String}
    account::String
    jobname::String
    stdout::String = ""
    stderr::String = ""
    interactive::Bool = false
    X11forwarding::Bool = interactive
    @assert !(X11forwarding && !interactive) "X11 forwarding requires interactive=true"
end

function pbs_header(q::QSub)
    pbs_directives = Any[
        "#PBS -l $(resource_list(q.resources))",
        "#PBS -A $(q.account)",
        "#PBS -N $(q.jobname)",
        "#PBS -o $(q.stdout)",
        "#PBS -e $(q.stderr)",
        "#PBS -j oe",
    ]
    if q.interactive
        push!(pbs_directives, "#PBS -I")
        push!(pbs_directives, "#PBS -q interactive_$(q.resources.ngpus !== nothing ? "gpu" : "cpu")")
    end
    if q.X11forwarding
        push!(pbs_directives, "#PBS -X")
    end

    """
    #!/bin/bash
    $(join(pbs_directives, "\n"))

    # Load modules
    $(join(["module load $mod" for mod in q.modules], "\n"))

    # Potentially useful environment variables exported by PBS
    #   PBS_JOBDIR      Pathname of job's staging and execution directory on the primary execution host
    #   PBS_JOBID       Job identifier given by PBS when the job is submitted
    #   PBS_JOBNAME     Job name specified by submitter
    #   PBS_NODEFILE    Name of file containing the list of nodes assigned to the job
    #   PBS_O_HOME      User's home directory. Value of HOME taken from user's submission environment
    #   PBS_O_PATH      User's PATH. Value of PATH taken from user's submission environment
    #   PBS_O_WORKDIR   Absolute path to directory where qsub is run. Value taken from user's submission environment
    #   TMPDIR          Pathname of job's scratch directory. Set when PBS assigns it
    """
end

#### JuliaCmd

@with_kw mutable struct JuliaCmd
    bindir::String
    threads::Int
    project::String
    script::String
    flags::Dict{String, Any} = Dict{String, Any}()
    args::Vector{String} = String[]
    env::Dict{String, Any} = Dict{String, Any}()
    secrets::Vector{String} = String[]
end

function try_find_julia_bindir()
    jldir = chomp(read(`which julia`, String)) # path of default julia
    jldir = chomp(read(`readlink -e $(jldir)`, String)) # follow symlinks
    bindir = dirname(jldir)
    @assert isdir(bindir)
    return bindir
end

function bash_commands(jl::JuliaCmd)
    envfile = logfile(JuliaJob(), jl.project, jl.script, "env.toml")
    flags = [isnothing(v) ? k : "$k=$v" for (k, v) in jl.flags]
    julia = join([joinpath(jl.bindir, "julia"); flags], " ")

    """
    # Ensure project directory path exists
    mkdir -p $(jl.project)
    cd $(jl.project)

    # Set environment variables
    $(join(["export $k=$v" for (k,v) in jl.env], "\n"))

    # Log environment variables
    $(julia) -e 'using TOML; TOML.print(ENV)' | sort -h $(
        isempty(jl.secrets) ?
            "> $(envfile)" :
            "| grep -Ev '^($(join(collect(jl.secrets), "|")))' > $(envfile)"
    )

    # Run julia script
    $(join([julia; jl.script; "--"; jl.args], " "))
    """
end

#### JuliaCmd + QSub

struct JuliaJob end

function Base.basename(::JuliaJob, script)
    @assert endswith(script, ".jl")
    basename(script)[begin:end-3]
end

function logfile(::JuliaJob, project, script, suffix)
    joinpath(
        mkpath(joinpath(project, "pbs")),
        basename(JuliaJob(), script) * "_" * suffix,
    )
end

function qsub(
        q::QSub,
        jl::JuliaCmd,
    )
    logfile_(suffix) = logfile(JuliaJob(), jl.project, jl.script, suffix)

    # Set standard output/standard error log files
    isempty(q.stdout) && (q.stdout = logfile_("stdout.txt"))
    isempty(q.stderr) && (q.stderr = logfile_("stderr.txt"))

    # Default environment
    get!(jl.env, "JULIA_BINDIR", jl.bindir)
    get!(jl.env, "JULIA_PROJECT", jl.project)
    get!(jl.env, "JULIA_NUM_THREADS", jl.threads)

    # Default flags
    get!(jl.flags, "--startup-file", "no")
    get!(jl.flags, "--history-file", "no")
    get!(jl.flags, "--optimize", nothing)
    get!(jl.flags, "--quiet", nothing)

    if q.interactive
        get!(jl.flags, "-i", nothing)
    end

    # Write pbs script to file
    open(logfile_("job.pbs"), "w+") do io
        print(io,
            """
            $(pbs_header(q))
            $(bash_commands(jl))
            """
        )
    end

    `qsub $(logfile_("job.pbs"))`
end

#### Copy project

function copy_project(;
        project,
        folders = [
            "src",
            "scripts",
        ],
        dest = mktempdir(),
    )
    files = joinpath.(project, [
        "Project.toml",
        "Manifest.toml",
        folders...,
    ])
    files = filter(ispath, files)

    for file in files
        cp(file, joinpath(dest, basename(file)); force=false, follow_symlinks=true)
    end

    return dest
end

end # module ARCTools


text:
# cartpole_sample.jl
using Gym

env = GymEnv("CartPole-v0")

action_space = env.action_space

# @show action_space.n

obs_space = env.observation_space

# @show obs_space.shape

reward = 0.0
episode_count = 10

for i=1:episode_count
    total = 0
    ob = reset!(env)
    render(env)#comment out this line if you do not want to visualize the environment
    while true
        action = sample(env.action_space)
        ob, reward, done, information = step!(env, action)
        total += reward
        render(env)#comment out this line if you do not want to visualize the environment
        done && break
    end
    println("episode $i total Rewards: $total")
end

close!(env)

max_stars_repo_name: JuliaTagBot/JuliaGymDemo

max_issues_repo_name: JuliaTagBot/JuliaGymDemo

max_forks_repo_name: JuliaTagBot/JuliaGymDemo


text:
# This file contains tests for input type

# A module which contains the code to
# perform tests on .Input module from Scats
module TestInput

using Formatting                 # Formatted strings
using Scats.Internal.Prec        # Precision module from Scats
using Scats: Internal.Input, API # API and .Input module from Scats
using Test                       # A package to perform tests

# Load the general utilities
include("utilities.jl")

# Print the header
println("\e[1;32mRUNNING\e[0m: input_test.jl")

"""
Read input data from a file
"""
macro input_read(file)
    return esc(:(s.Input.read!($file)))
end

"""
Generate an example of a file containing input data
"""
macro input_example(file)
    return esc(:(s.Input.example($file)))
end

"""
Reset the values of the input data
"""
macro input_reset()
    return esc(:(s.Input.reset!()))
end

"""
Test the values of input data
"""
macro test_values()
    return esc(quote
        @test s.Input.N == 230
        @test s.Input.Δt == 1
        @test s.Input.q == 0.01
        @test s.Input.t == [ t for t in 0:229 ]
        @test s.Input.x == [ x for x in 0:229 ]
    end)
end

"""
Create a temporary file, write exemplary
generator parameters in it and read them
"""
macro read_example_data()
    return esc(quote
        @file
        @input_example file
        @input_read file
    end)
end

"""
Test the end of file exception on an empty file
"""
macro test_eof_exception()
    return esc(quote
        try
            @input_read file
        catch e
            @test e isa Input.ScatsInputEOF
            @test sprint(showerror, e) == retrieve_messages(
                Input.ScatsInputEOF, file
            )
        end
    end)
end

@construct_exceptions "Input.ScatsInput"

# Test exceptions related to file status
@testset "Check file status" begin

    @instantiate

    # Test when a wrong path has been specified
    try
        @input_read "Wrong file path!"
    catch e
        @test e isa Input.ScatsInputNotAFile
        @test sprint(showerror, e) == retrieve_messages(
            Input.ScatsInputNotAFile, "Wrong file path!"
        )
    end

    @file
    @test_eof_exception

    # Test the end of file exception on files with different number of lines
    for i in 1:13

        # Generate valid data on input lines
        if !(i in range(2, 11, step = 3))
            println(io, "Line ", i)
        elseif i == 2
            println(io, 1)
        else
            println(io, RT(1))
        end

        # Update the file immediately
        flush(io)

        @test_eof_exception

    end

end

# Test different ways to read input data
@testset "Read `good` data" begin

    @instantiate

    @read_example_data

    @test_values
    @input_reset

    # Read the input data (another way)
    s.Input.read!(file)

    @test_values
    @input_reset

    # Read the input data (another way)
    s.Input(file)

    @test_values

end

# Test the creation of an example of input data
@testset "Check creation of an example" begin

    @instantiate
    @read_example_data
    @test_values

    # Create a temporary directory
    dir = mktempdir()

    # Test the exception being thrown when a directory has been passed
    try
        @input_example dir
    catch e
        @test e isa Input.ScatsInputIsADir
        @test sprint(showerror, e) == retrieve_messages(
            Input.ScatsInputIsADir, dir,
        )
    end

    # Write input data (another way)
    s.input_example(file)

    @input_read file
    @test_values

    # Test the exception being thrown when a directory has been passed
    try
        s.input_example(dir)
    catch e
        @test e isa Input.ScatsInputIsADir
        @test sprint(showerror, e) == retrieve_messages(
            Input.ScatsInputIsADir, dir,
        )
    end

end

# Test all exceptions related to invalid input
@testset "Read `bad` data" begin

    @instantiate

    # Construct an ordered list of exceptions
    exceptions = construct_exceptions(
        "WR", "N", "Δt", "q", "t", "x"
    )

    # Get an ordered list of expected messages
    messages = retrieve_messages(exceptions, "input")

    # Create a temporary file to contain valid parameters
    good_file, _ = mktemp()

    @input_example good_file

    # Corrupt a file on a specific line
    @inline function break_a_line!(ln::Int)

        @file

        open(good_file, "r") do good_io

            k = 0

            # Print lines of the `good` file into the `bad` file
            for line in eachline(good_io)

                k += 1

                # Corrupt a line
                if k == ln
                    line = "Hello there!\n"
                end

                # Print a line
                println(io, line)

            end

        end

        # Close the stream
        close(io)

        # Save the temporary file as `input`
        cp(file, "input", force = true)

    end

    # Test all cases of corrupted input
    for i in 1:5

        # Calculate the number of a line to corrupt
        n = 2 + (i - 1) * 3

        # Corrupt a specific line
        break_a_line!(n)

        # Test an exception
        try
            @input_read "input"
        catch e
            @test e isa exceptions[i]
            @test sprint(showerror, e) == string(messages[i])
        end

    end

    # Remove the `input` file
    isfile("input") && rm("input")

end

# Test writing input data to a file
@testset "Check writing" begin

    @instantiate

    # Specify expected values to be writed
    contents = [
        "230", " 1.000000000000000E+00", " 1.000000000000000E-02",
        join([ sprintf1(Prec.RF, s) for s in 0:229 ], " "^3),
        join([ sprintf1(Prec.RF, s) for s in 0:229 ], " "^3)
    ]

    # Create a temporary file to contain input data (from example)
    @file

    @input_example file
    @input_read file

    # Create another temporary file to contain input data (written)
    @file

    # Write input data
    s.write_input(file)

    # Read lines from the written file
    lines = readlines(file)

    # Test the values
    for i in 1:5
        n = 2 + (i - 1) * 3
        @test lines[n] == contents[i]
    end

    # Write input data (another way)
    s.Input.write(file)

    # Read lines from the written file
    lines = readlines(file)

    # Test the values
    for i in 1:5
        n = 2 + (i - 1) * 3
        @test lines[n] == contents[i]
    end

end

# Test values resetting
@testset "Check resetting" begin

    @instantiate
    @read_example_data

    @input_reset

    # Test values
    @test s.Input.N == 0
    @test s.Input.Δt == 0
    @test s.Input.q == 0
    @test s.Input.t == []
    @test s.Input.x == []

end

end


text:
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Elements
import ClimateMachine.Mesh.Elements: baryweights
using GaussQuadrature
using Base.Threads

# Depending on CliMa version
# old, should return a tuple of polynomial orders
# polynomialorders(::DiscontinuousSpectralElementGrid{T, dim, N}) where {T, dim, N} = Tuple([N for i in 1:dim])
# new, should return a tuple of polynomial orders
polynomialorders(::DiscontinuousSpectralElementGrid{T, dim, N}) where {T, dim, N} = N

# utils.jl
"""
function cellaverage(Q; M = nothing)

# Description
Compute the cell-average of Q given the mass matrix M.
Assumes that Q and M are the same size

# Arguments
`Q`: MPIStateArrays (array)

# Keyword Arguments
`M`: Mass Matrix (array)

# Return
The cell-average of Q
"""
function cellaverage(Q; M = nothing)
    if M==nothing
        return nothing
    end
    return (sum(M .* Q, dims = 1) ./ sum(M , dims = 1))[:]
end

"""
function coordinates(grid::DiscontinuousSpectralElementGrid)

# Description
Gets the (x,y,z) coordinates corresponding to the grid

# Arguments
- `grid`: DiscontinuousSpectralElementGrid

# Return
- `x, y, z`: views of x, y, z coordinates
"""
function coordinates(grid::DiscontinuousSpectralElementGrid)
    x = view(grid.vgeo, :, grid.x1id, :)   # x-direction
    y = view(grid.vgeo, :, grid.x2id, :)   # y-direction
    z = view(grid.vgeo, :, grid.x3id, :)   # z-direction
    return x, y, z
end

"""
function cellcenters(Q; M = nothing)

# Description
Get the cell-centers of every element in the grid

# Arguments
- `grid`: DiscontinuousSpectralElementGrid

# Return
- Tuple of cell-centers
"""
function cellcenters(grid::DiscontinuousSpectralElementGrid)
    x, y, z = coordinates(grid)
    M = view(grid.vgeo, :, grid.Mid, :)  # mass matrix
    xC = cellaverage(x, M = M)
    yC = cellaverage(y, M = M)
    zC = cellaverage(z, M = M)
    return xC[:], yC[:], zC[:]
end

# find_element.jl
# 3D version
function findelement(xC, yC, zC, location, p, lin)
    ex, ey, ez = size(lin)
    # i
    currentmin = ones(1)
    minind = ones(Int64, 1)
    currentmin[1] = abs.(xC[p[lin[1,1,1]]] .- location[1])
    for i in 2:ex
        current = abs.(xC[p[lin[i,1,1]]] .- location[1])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    i = minind[1]
    # j
    currentmin[1] = abs.(yC[p[lin[1,1,1]]] .- location[2])
    minind[1] = 1
    for i in 2:ey
        current = abs.(yC[p[lin[1,i,1]]] .- location[2])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    j = minind[1]
    # k
    currentmin[1] = abs.(zC[p[lin[1,1,1]]] .- location[3])
    minind[1] = 1
    for i in 2:ez
        current = abs.(zC[p[lin[1,1, i]]] .- location[3])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    k = minind[1]
    return p[lin[i,j,k]]
end

# 2D version
function findelement(xC, yC, location, p, lin)
    ex, ey = size(lin)
    # i
    currentmin = ones(1)
    minind = ones(Int64, 1)
    currentmin[1] = abs.(xC[p[lin[1,1]]] .- location[1])
    for i in 2:ex
        current = abs.(xC[p[lin[i,1]]] .- location[1])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    i = minind[1]
    # j
    currentmin[1] = abs.(yC[p[lin[1,1]]] .- location[2])
    minind[1] = 1
    for i in 2:ey
        current = abs.(yC[p[lin[1,i]]] .- location[2])
        if current < currentmin[1]
            currentmin[1] = current
            minind[1] = i
        end
    end
    j = minind[1]
    return p[lin[i,j]]
end

# gridhelper.jl

struct InterpolationHelper{S, T}
    points::S
    quadrature::S
    interpolation::S
    cartesianindex::T
end

function InterpolationHelper(g::DiscontinuousSpectralElementGrid)
    porders = polynomialorders(g)
    if length(porders) == 3
        npx, npy, npz = porders
        rx, wx = GaussQuadrature.legendre(npx+1, both)
        ωx = baryweights(rx)
        ry, wy = GaussQuadrature.legendre(npy+1, both)
        ωy = baryweights(ry)
        rz, wz = GaussQuadrature.legendre(npz+1, both)
        ωz = baryweights(rz)
        linlocal = reshape(collect(1:(npx+1)*(npy+1)*(npz+1)), (npx+1, npy+1, npz+1))
        return InterpolationHelper((rx, ry, rz), (wx, wy, wz), (ωx, ωy, ωz), linlocal)
    elseif length(porders) == 2
        npx, npy = porders
        rx, wx = GaussQuadrature.legendre(npx+1, both)
        ωx = baryweights(rx)
        ry, wy = GaussQuadrature.legendre(npy+1, both)
        ωy = baryweights(ry)
        linlocal = reshape(collect(1:(npx+1)*(npy+1)), (npx+1, npy+1))
        return InterpolationHelper((rx, ry), (wx, wy), (ωx, ωy), linlocal)
    else
        println("Not supported")
        return nothing
    end
    return nothing
end

struct ElementHelper{S, T, U, Q, V, W}
    cellcenters::S
    coordinates::T
    cartesiansizes::U
    polynomialorders::Q
    permutation::V
    cartesianindex::W
end

addup(xC, tol) = sum(abs.(xC[1] .- xC) .≤ tol)

# only valid for cartesian domains
function ElementHelper(g::DiscontinuousSpectralElementGrid)
    porders = polynomialorders(g)
    x, y, z = coordinates(g)
    xC, yC, zC = cellcenters(g)
    ne = size(x)[2]
    ex = round(Int64, ne / addup(xC, 10^4 * eps(maximum(abs.(x)))))
    ey = round(Int64, ne / addup(yC, 10^4 * eps(maximum(abs.(y)))))
    ez = round(Int64, ne / addup(zC, 10^4 * eps(maximum(abs.(z)))))

    check = ne == ex * ey * ez
    check ? true : error("improper counting")
    p = getperm(xC, yC, zC, ex, ey, ez)
    # should use dispatch ...
    if length(porders) == 3
        npx, npy, npz = porders
        lin = reshape(collect(1:length(xC)), (ex, ey, ez))
        return ElementHelper((xC, yC, zC), (x, y, z), (ex, ey, ez), porders, p, lin)
    elseif length(porders) == 2
        npx, npy = porders
        lin = reshape(collect(1:length(xC)), (ex, ey))
        check = ne == ex * ey
        check ? true : error("improper counting")
        return ElementHelper((xC, yC), (x, y), (ex, ey), porders, p, lin)
    else
        println("no constructor for polynomial order = ", porders)
        return nothing
    end
    return nothing
end

struct GridHelper{S,T,V}
    interpolation::S
    element::T
    grid::V
end

function GridHelper(g::DiscontinuousSpectralElementGrid)
    return GridHelper(InterpolationHelper(g), ElementHelper(g), g)
end

function getvalue(f, location, gridhelper::GridHelper)
    ih = gridhelper.interpolation
    eh = gridhelper.element
    porders = gridhelper.element.polynomialorders
    if length(porders) == 3
        npx, npy, npz = gridhelper.element.polynomialorders
        fl = reshape(f, (npx+1, npy+1, npz+1, prod(eh.cartesiansizes)))
        ip = getvalue(fl, eh.cellcenters..., location,
                eh.permutation, eh.cartesianindex, ih.cartesianindex,
                eh.coordinates..., ih.points..., ih.interpolation...)
        return ip
    elseif length(porders) == 2
        npx, npy = gridhelper.element.polynomialorders
        fl = reshape(f, (npx+1, npy+1, prod(eh.cartesiansizes)))
        ip = getvalue(fl, eh.cellcenters..., location,
                eh.permutation, eh.cartesianindex, ih.cartesianindex,
                eh.coordinates..., ih.points..., ih.interpolation...)
        return ip
    end
    return nothing
end

# lagrange_interpolation.jl
function checkgl(x, rx)
    for i in eachindex(rx)
        if abs(x - rx[i]) ≤ eps(rx[i])  return i end
    end
    return 0
end

function lagrange_eval(f, newx, newy, newz, rx, ry, rz, ωx, ωy, ωz)
    icheck = checkgl(newx, rx)
    jcheck = checkgl(newy, ry)
    kcheck = checkgl(newz, rz)
    numerator = zeros(1)
    denominator = zeros(1)
    for k in eachindex(rz)
        if kcheck ==0
            Δz = (newz .- rz[k])
            polez = ωz[k] ./ Δz
            kk = k
        else
            polez = 1.0
            k = eachindex(rz)[end]
            kk = kcheck
        end
        for j in eachindex(ry)
            if jcheck ==0
                Δy = (newy .- ry[j])
                poley = ωy[j] ./ Δy
                jj = j
            else
                poley = 1.0
                j = eachindex(ry)[end]
                jj = jcheck
            end
            for i in eachindex(rx)
                if icheck == 0
                    Δx = (newx .- rx[i])
                    polex = ωx[i] ./ Δx
                    ii = i
                else
                    polex = 1.0
                    i = eachindex(rx)[end]
                    ii = icheck
                end
                numerator[1] += f[ii,jj,kk] * polex * poley * polez
                denominator[1] += polex * poley * polez
            end
        end
    end
    return numerator[1] / denominator[1]
end

function lagrange_eval(f, newx, newy, rx, ry, ωx, ωy)
    icheck = checkgl(newx, rx)
    jcheck = checkgl(newy, ry)
    numerator = zeros(1)
    denominator = zeros(1)
    for j in eachindex(ry)
        if jcheck ==0
            Δy = (newy .- ry[j])
            poley = ωy[j] ./ Δy
            jj = j
        else
            poley = 1.0
            j = eachindex(ry)[end]
            jj = jcheck
        end
        for i in eachindex(rx)
            if icheck == 0
                Δx = (newx .- rx[i])
                polex = ωx[i] ./ Δx
                ii = i
            else
                polex = 1.0
                i = eachindex(rx)[end]
                ii = icheck
            end
            numerator[1] += f[ii,jj] * polex * poley
            denominator[1] += polex * poley
        end
    end
    return numerator[1] / denominator[1]
end

function lagrange_eval_nocheck(f, newx, newy, newz, rx, ry, rz, ωx, ωy, ωz)
    numerator = zeros(1)
    denominator = zeros(1)
    for k in eachindex(rz)
            Δz = (newz .- rz[k])
            polez = ωz[k] ./ Δz
        for j in eachindex(ry)
                Δy = (newy .- ry[j])
                poley = ωy[j] ./ Δy
            for i in eachindex(rx)
                    Δx = (newx .- rx[i])
                    polex = ωx[i] ./ Δx
                    numerator[1] += f[i,j,k] * polex * poley * polez
                    denominator[1] += polex * poley * polez
            end
        end
    end
    return numerator[1] / denominator[1]
end

function lagrange_eval_nocheck(f, newx, newy, rx, ry, ωx, ωy)
    numerator = zeros(1)
    denominator = zeros(1)
    for j in eachindex(ry)
            Δy = (newy .- ry[j])
            poley = ωy[j] ./ Δy
        for i in eachindex(rx)
                Δx = (newx .- rx[i])
                polex = ωx[i] ./ Δx
                numerator[1] += f[i,j] * polex * poley
                denominator[1] += polex * poley
        end
    end
    return numerator[1] / denominator[1]
end

function lagrange_eval_nocheck(f, newx, rx, ωx)
    numerator = zeros(1)
    denominator = zeros(1)
    for i in eachindex(rx)
            Δx = (newx .- rx[i])
            polex = ωx[i] ./ Δx
            numerator[1] += f[i] * polex * poley
            denominator[1] += polex * poley
    end
    return numerator[1] / denominator[1]
end

# 3D, only valid for rectangles
function getvalue(fl, xC, yC, zC, location, p, lin, linlocal, x, y, z, rx, ry, rz, ωx, ωy, ωz)
    e = findelement(xC, yC, zC, location, p, lin)
    # need bounds to rescale, only value for cartesian

    xmax = x[linlocal[length(rx),1,1], e]
    xmin = x[linlocal[1,1,1], e]
    ymax = y[linlocal[1,length(ry),1], e]
    ymin = y[linlocal[1,1,1], e]
    zmax = z[linlocal[1,1,length(rz)], e]
    zmin = z[linlocal[1,1,1], e]

    # rescale new point to [-1,1]³
    newx = 2 * (location[1] - xmin) / (xmax - xmin) - 1
    newy = 2 * (location[2] - ymin) / (ymax - ymin) - 1
    newz = 2 * (location[3] - zmin) / (zmax - zmin) - 1

    return lagrange_eval(view(fl,:,:,:,e), newx, newy, newz, rx, ry, rz, ωx, ωy, ωz)
end

# 2D
function getvalue(fl, xC, yC, location, p, lin, linlocal, x, y, rx, ry,  ωx, ωy)
    e = findelement(xC, yC, location, p, lin)
    # need bounds to rescale
    xmax = x[linlocal[length(rx),1,1], e]
    xmin = x[linlocal[1,1,1], e]
    ymax = y[linlocal[1,length(ry),1], e]
    ymin = y[linlocal[1,1,1], e]

    # rescale new point to [-1,1]²
    newx = 2 * (location[1] - xmin) / (xmax - xmin) - 1
    newy = 2 * (location[2] - ymin) / (ymax - ymin) - 1

    return lagrange_eval(view(fl,:,:,e), newx, newy, rx, ry, ωx, ωy)
end

# permutations.jl
function getperm(xC, yC, zC, ex, ey, ez)
    pz = sortperm(zC)
    tmpY = reshape(yC[pz], (ex*ey, ez))
    tmp_py = [sortperm(tmpY[:,i]) for i in 1:ez]
    py = zeros(Int64,length(pz))
    for i in eachindex(tmp_py)
        n = length(tmp_py[i])
        ii = (i-1) * n + 1
        py[ii : ii + n-1] .= tmp_py[i] .+ ii .- 1
    end
    tmpX = reshape(xC[pz][py], (ex, ey*ez))
    tmp_px = [sortperm(tmpX[:,i]) for i in 1:ey*ez]
    px = zeros(Int64,length(pz))
    for i in eachindex(tmp_px)
        n = length(tmp_px[i])
        ii = (i-1) * n + 1
        px[ii : ii + n-1] .= tmp_px[i] .+ ii .- 1
    end
    p = [pz[py[px[i]]] for i in eachindex(px)]
    return p
end


text:
module TestWDICountries

using Base.Test
using DataFrames
using WorldBankDataTd
using Dates

## the worldbank database shows some peculiarities
url =
    "http://api.worldbank.org/countries/afr?per_page=25000&format=json"

## missing values for aggregates are sometimes denoted by "NA"
json = [WorldBankDataTd.download_parse_json(url)[2]]
@test json[1]["incomeLevel"]["id"] == "NA"
@test json[1]["region"]["id"] == "NA"

## and sometimes denoted by ""
@test json[1]["lendingType"]["id"] == ""
@test json[1]["latitude"] == ""
@test json[1]["longitude"] == ""
@test json[1]["adminregion"]["id"] == ""
@test json[1]["adminregion"]["value"] == ""
@test json[1]["capitalCity"] == ""

## iso2c is "NA" for Namibia
url =
    "http://api.worldbank.org/countries/nam?per_page=25000&format=json"

json = [WorldBankDataTd.download_parse_json(url)[2]]
@test json[1]["iso2Code"] == "NA"

end