text:
# using Algames# using Algames
# using BenchmarkTools
# using LinearAlgebra
# using StaticArrays
# using TrajectoryOptimization
# const TO = TrajectoryOptimization
# const AG = Algames
#
# # Define the dynamics model of the game.
# struct InertialUnicycleGame{T} <: AbstractGameModel
#     n::Int  # Number of states
#     m::Int  # Number of controls
#     mp::T
# 	pu::Vector{Vector{Int}} # Indices of the each player's controls
# 	px::Vector{Vector{Int}} # Indices of the each player's x and y positions
#     p::Int  # Number of players
# end
# InertialUnicycleGame() = InertialUnicycleGame(
# 	12,
# 	6,
# 	1.0,
# 	[[1,2],[3,4],[5,6]],
# 	[[1,2],[5,6],[9,10]],
# 	3)
# Base.size(::InertialUnicycleGame) = 12,6,[[1,2],[3,4],[5,6]],3 # n,m,pu,p
#
# # Instantiate dynamics model
# model = InertialUnicycleGame()
# n,m,pu,p = size(model)
# T = Float64
# px = model.px
# function TO.dynamics(model::InertialUnicycleGame, x, u) # Non memory allocating dynamics
#     qd1 = @SVector [cos(x[3]), sin(x[3])]
#     qd1 *= x[4]
#     qd2 = @SVector [cos(x[7]), sin(x[7])]
#     qd2 *= x[8]
#     qd3 = @SVector [cos(x[11]), sin(x[11])]
#     qd3 *= x[12]
#     qdd1 = u[ @SVector [1,2] ]
#     qdd2 = u[ @SVector [3,4] ]
#     qdd3 = u[ @SVector [5,6] ]
#     return [qd1; qdd1; qd2; qdd2; qd3; qdd3]
# end
#
# # Discretization info
# tf = 3.0  # final time
# N = 21    # number of knot points
# dt = tf / (N-1) # time step duration
#
# # Define initial and final states (be sure to use Static Vectors!)
# x0 = @SVector [
#                -0.80, -0.05,  0.00, 0.90, #lane1
#                -1.00, -0.05,  0.00, 0.90, #lane2
#                -0.90, -0.30, pi/12, 0.95, #lane4
#                 ]
# xf = @SVector [
#                 1.10, -0.05,  0.00, 0.90, #lane1
#                 0.70, -0.05,  0.00, 0.90, #lane2
#                 0.90, -0.05,  0.00, 0.90, #lane4
#                ]
#
# # Define a quadratic cost
# diag_Q1 = @SVector [ # Player 1 state cost
#     0., 10., 1., 1.,
#     0., 0., 0., 0.,
#     0., 0., 0., 0.]
# diag_Q2 = @SVector [ # Player 2 state cost
#     0., 0., 0., 0.,
#     0., 10., 1., 1.,
#     0., 0., 0., 0.]
# diag_Q3 = @SVector [ # Player 3 state cost
#     0., 0., 0., 0.,
#     0., 0., 0., 0.,
#     0., 10., 1., 1.]
# Q = [0.1*Diagonal(diag_Q1), # Players state costs
#      0.1*Diagonal(diag_Q2),
#      0.1*Diagonal(diag_Q3)]
# Qf = [1.0*Diagonal(diag_Q1),
#       1.0*Diagonal(diag_Q2),
#       1.0*Diagonal(diag_Q3)]
#
# # Players controls costs
# R = [0.1*Diagonal(@SVector ones(length(pu[1]))),
#      0.1*Diagonal(@SVector ones(length(pu[2]))),
#      0.1*Diagonal(@SVector ones(length(pu[3]))),
#     ]
#
# # Players objectives
# obj = [LQRObjective(Q[i],R[i],Qf[i],xf,N) for i=1:p]
#
# # Define the initial trajectory
# xs = SVector{n}(zeros(n))
# us = SVector{m}(zeros(m))
# Z = [KnotPoint(xs,us,dt) for k = 1:N]
# Z[end] = KnotPoint(xs,m)
#
# # Build problem
# actor_radius = 0.08
# actors_radii = [actor_radius for i=1:p]
# actors_types = [:car for i=1:p]
# road_length = 6.20
# road_width = 0.34
# ramp_length = 3.2
# ramp_angle = pi/12
# scenario = MergingScenario(road_length, road_width, ramp_length,
# 	ramp_angle, actors_radii, actors_types)
#
#
# # Create constraint sets
# algames_conSet = ConstraintSet(n,m,N)
# ilqgames_conSet = ConstraintSet(n,m,N)
# con_inds = 1:N # Indices where the constraints will be applied
#
# # Add collision avoidance constraints
# add_collision_avoidance(algames_conSet, actors_radii, px, p, con_inds)
# add_collision_avoidance(ilqgames_conSet, actors_radii, px, p, con_inds)
#
# # Add scenario specific constraints (road boundaries)
# add_scenario_constraints(algames_conSet, scenario, px, con_inds; constraint_type=:constraint)
# add_scenario_constraints(ilqgames_conSet, scenario, px, con_inds; constraint_type=:constraint);
#
# algames_prob = GameProblem(model, obj, algames_conSet, x0, xf, Z, N, tf)
# ilqgames_prob = GameProblem(model, obj, ilqgames_conSet, x0, xf, Z, N, tf);
#
# algames_opts = DirectGamesSolverOptions{T}(
#     iterations=10,
# 	min_steps_per_iteration=1,
# 	record_condition=false,
#     inner_iterations=20,
#     iterations_linesearch=10,
#     log_level=Algames.Logging.Debug)
# algames_solver = DirectGamesSolver(algames_prob, algames_opts)
#
# ilqgames_opts = PenaltyiLQGamesSolverOptions{T}(
#     iterations=200,
#     gradient_norm_tolerance=1e-2,
#     cost_tolerance=1e-4,
#     # iterations_linesearch=5,
#     line_search_lower_bound=0.0,
#     line_search_upper_bound=0.02,
#     log_level=Algames.Logging.Warn)
# ilqgames_solver = PenaltyiLQGamesSolver(ilqgames_prob, ilqgames_opts)
# pen = ones(length(ilqgames_solver.constraints))*100.0
# set_penalty!(ilqgames_solver, pen)
#
# @time solve!(algames_solver);
# @time solve!(ilqgames_solver);
#
# # Define a noise that will be added to the initial state of the system x0.
# # By sampling x0 randomly we generate a class of similar dynamic games.
# # In this case, we add noise to the initial speeds, positions and orientations of the vehicles.
# state_noise = @SVector [ # Uniform noise around x0
#     0.06, 0.06, 2*pi/72, 0.05,
#     0.06, 0.06, 2*pi/72, 0.05,
#     0.06, 0.06, 2*pi/72, 0.05]
#
# opts_monte_carlo = MonteCarloSamplerOptions{n,T}(
#     noise=state_noise, # noise added to the initial state
#     iterations=1000) # number of Monte Carlo samples
#
# algames_sampler = MonteCarloSampler(algames_solver, opts_monte_carlo)
# ilqgames_sampler = MonteCarloSampler(ilqgames_solver, opts_monte_carlo)
#
# monte_carlo_sampling(algames_sampler)
# monte_carlo_sampling(ilqgames_sampler)
#
# # Algames
# # Time required to solve the problem
# visualize_solve_time(algames_sampler; save=false)
#
# # Maximum constraint violation
# visualize_cmax(algames_sampler; save=false)
#
# # Number of iterations (Newton step for Algames, Riccati backwardpass for iLQGames)
# visualize_iterations_total(algames_sampler; save=false)
#
# # Optimality constraint satisfaction
# visualize_optimality_merit(algames_sampler; save=false)
#
#
# # ilQGames
# # Time required to solve the problem
# visualize_solve_time(ilqgames_sampler; save=false)
#
# # Maximum constraint violation
# visualize_cmax(ilqgames_sampler; save=false)
#
# # Number of iterations (Newton step for Algames, Riccati backwardpass for iLQGames)
# visualize_iterations_total(ilqgames_sampler; save=false)

# using BenchmarkTools
# using LinearAlgebra
# using StaticArrays
# using TrajectoryOptimization
# const TO = TrajectoryOptimization
# const AG = Algames
#
# # Define the dynamics model of the game.
# struct InertialUnicycleGame{T} <: AbstractGameModel
#     n::Int  # Number of states
#     m::Int  # Number of controls
#     mp::T
# 	pu::Vector{Vector{Int}} # Indices of the each player's controls
# 	px::Vector{Vector{Int}} # Indices of the each player's x and y positions
#     p::Int  # Number of players
# end
# InertialUnicycleGame() = InertialUnicycleGame(
# 	12,
# 	6,
# 	1.0,
# 	[[1,2],[3,4],[5,6]],
# 	[[1,2],[5,6],[9,10]],
# 	3)
# Base.size(::InertialUnicycleGame) = 12,6,[[1,2],[3,4],[5,6]],3 # n,m,pu,p
#
# # Instantiate dynamics model
# model = InertialUnicycleGame()
# n,m,pu,p = size(model)
# T = Float64
# px = model.px
# function TO.dynamics(model::InertialUnicycleGame, x, u) # Non memory allocating dynamics
#     qd1 = @SVector [cos(x[3]), sin(x[3])]
#     qd1 *= x[4]
#     qd2 = @SVector [cos(x[7]), sin(x[7])]
#     qd2 *= x[8]
#     qd3 = @SVector [cos(x[11]), sin(x[11])]
#     qd3 *= x[12]
#     qdd1 = u[ @SVector [1,2] ]
#     qdd2 = u[ @SVector [3,4] ]
#     qdd3 = u[ @SVector [5,6] ]
#     return [qd1; qdd1; qd2; qdd2; qd3; qdd3]
# end
#
# # Discretization info
# tf = 3.0  # final time
# N = 21    # number of knot points
# dt = tf / (N-1) # time step duration
#
# # Define initial and final states (be sure to use Static Vectors!)
# x0 = @SVector [
#                -0.80, -0.05,  0.00, 0.90, #lane1
#                -1.00, -0.05,  0.00, 0.90, #lane2
#                -0.90, -0.30, pi/12, 0.95, #lane4
#                 ]
# xf = @SVector [
#                 1.10, -0.05,  0.00, 0.90, #lane1
#                 0.70, -0.05,  0.00, 0.90, #lane2
#                 0.90, -0.05,  0.00, 0.90, #lane4
#                ]
#
# # Define a quadratic cost
# diag_Q1 = @SVector [ # Player 1 state cost
#     0., 10., 1., 1.,
#     0., 0., 0., 0.,
#     0., 0., 0., 0.]
# diag_Q2 = @SVector [ # Player 2 state cost
#     0., 0., 0., 0.,
#     0., 10., 1., 1.,
#     0., 0., 0., 0.]
# diag_Q3 = @SVector [ # Player 3 state cost
#     0., 0., 0., 0.,
#     0., 0., 0., 0.,
#     0., 10., 1., 1.]
# Q = [0.1*Diagonal(diag_Q1), # Players state costs
#      0.1*Diagonal(diag_Q2),
#      0.1*Diagonal(diag_Q3)]
# Qf = [1.0*Diagonal(diag_Q1),
#       1.0*Diagonal(diag_Q2),
#       1.0*Diagonal(diag_Q3)]
#
# # Players controls costs
# R = [0.1*Diagonal(@SVector ones(length(pu[1]))),
#      0.1*Diagonal(@SVector ones(length(pu[2]))),
#      0.1*Diagonal(@SVector ones(length(pu[3]))),
#     ]
#
# # Players objectives
# obj = [LQRObjective(Q[i],R[i],Qf[i],xf,N) for i=1:p]
#
# # Define the initial trajectory
# xs = SVector{n}(zeros(n))
# us = SVector{m}(zeros(m))
# Z = [KnotPoint(xs,us,dt) for k = 1:N]
# Z[end] = KnotPoint(xs,m)
#
# # Build problem
# actor_radius = 0.08
# actors_radii = [actor_radius for i=1:p]
# actors_types = [:car for i=1:p]
# road_length = 6.20
# road_width = 0.34
# ramp_length = 3.2
# ramp_angle = pi/12
# scenario = MergingScenario(road_length, road_width, ramp_length,
# 	ramp_angle, actors_radii, actors_types)
#
#
# # Create constraint sets
# algames_conSet = ConstraintSet(n,m,N)
# ilqgames_conSet = ConstraintSet(n,m,N)
# con_inds = 1:N # Indices where the constraints will be applied
#
# # Add collision avoidance constraints
# add_collision_avoidance(algames_conSet, actors_radii, px, p, con_inds)
# add_collision_avoidance(ilqgames_conSet, actors_radii, px, p, con_inds)
#
# # Add scenario specific constraints (road boundaries)
# add_scenario_constraints(algames_conSet, scenario, px, con_inds; constraint_type=:constraint)
# add_scenario_constraints(ilqgames_conSet, scenario, px, con_inds; constraint_type=:constraint);
#
# algames_prob = GameProblem(model, obj, algames_conSet, x0, xf, Z, N, tf)
# ilqgames_prob = GameProblem(model, obj, ilqgames_conSet, x0, xf, Z, N, tf);
#
# algames_opts = DirectGamesSolverOptions{T}(
#     iterations=10,
# 	min_steps_per_iteration=1,
# 	record_condition=false,
#     inner_iterations=20,
#     iterations_linesearch=10,
#     log_level=Algames.Logging.Debug)
# algames_solver = DirectGamesSolver(algames_prob, algames_opts)
#
# ilqgames_opts = PenaltyiLQGamesSolverOptions{T}(
#     iterations=200,
#     gradient_norm_tolerance=1e-2,
#     cost_tolerance=1e-4,
#     # iterations_linesearch=5,
#     line_search_lower_bound=0.0,
#     line_search_upper_bound=0.02,
#     log_level=Algames.Logging.Warn)
# ilqgames_solver = PenaltyiLQGamesSolver(ilqgames_prob, ilqgames_opts)
# pen = ones(length(ilqgames_solver.constraints))*100.0
# set_penalty!(ilqgames_solver, pen)
#
# @time solve!(algames_solver);
# @time solve!(ilqgames_solver);
#
# # Define a noise that will be added to the initial state of the system x0.
# # By sampling x0 randomly we generate a class of similar dynamic games.
# # In this case, we add noise to the initial speeds, positions and orientations of the vehicles.
# state_noise = @SVector [ # Uniform noise around x0
#     0.06, 0.06, 2*pi/72, 0.05,
#     0.06, 0.06, 2*pi/72, 0.05,
#     0.06, 0.06, 2*pi/72, 0.05]
#
# opts_monte_carlo = MonteCarloSamplerOptions{n,T}(
#     noise=state_noise, # noise added to the initial state
#     iterations=1000) # number of Monte Carlo samples
#
# algames_sampler = MonteCarloSampler(algames_solver, opts_monte_carlo)
# ilqgames_sampler = MonteCarloSampler(ilqgames_solver, opts_monte_carlo)
#
# monte_carlo_sampling(algames_sampler)
# monte_carlo_sampling(ilqgames_sampler)
#
# # Algames
# # Time required to solve the problem
# visualize_solve_time(algames_sampler; save=false)
#
# # Maximum constraint violation
# visualize_cmax(algames_sampler; save=false)
#
# # Number of iterations (Newton step for Algames, Riccati backwardpass for iLQGames)
# visualize_iterations_total(algames_sampler; save=false)
#
# # Optimality constraint satisfaction
# visualize_optimality_merit(algames_sampler; save=false)
#
#
# # ilQGames
# # Time required to solve the problem
# visualize_solve_time(ilqgames_sampler; save=false)
#
# # Maximum constraint violation
# visualize_cmax(ilqgames_sampler; save=false)
#
# # Number of iterations (Newton step for Algames, Riccati backwardpass for iLQGames)
# visualize_iterations_total(ilqgames_sampler; save=false)


text:
"""
    RiemannianDistance(M::Manifold)

An implementation of the `Metric` interface from the Distances.jl package that
returns Riemannian distance between points stored in a linearized array of coordinates.
"""
struct RiemannianDistance{TM<:Manifold} <: Metric
    manifold::TM
end

function (dist::RiemannianDistance)(a, b)
    rep_size = representation_size(dist.manifold)
    p = reshape(a, rep_size)
    q = reshape(b, rep_size)
    return distance(dist.manifold, p, q)
end

ext: jl

max_stars_repo_name: JuliaManifolds/ManifoldML

max_issues_repo_name: JuliaManifolds/ManifoldML

max_forks_repo_name: JuliaManifolds/ManifoldML


text:
push!(LOAD_PATH, "./jl")

using Base.Test
using RBT: Iso, Opt, Seq, Rel, Env, EnvRel, Unit, Zero
using DataFrames

int_here = RBT.HerePipe(Int)
println(repr(int_here))

@test int_here(5) == 5
@test_throws Exception int_here("hi there")

any_here = RBT.HerePipe(Any)
println(repr(any_here))

@test any_here(5) == 5
@test any_here("hi there") == "hi there"

unit_here = RBT.HerePipe(Unit)
println(repr(unit_here))

@test unit_here() == nothing

int_unit = RBT.UnitPipe(Int)
println(repr(int_unit))

@test int_unit(5) == nothing
@test_throws Exception int_unit(nothing)

unit_unit = RBT.UnitPipe(Unit)
println(repr(unit_unit))

@test unit_unit() == nothing
@test_throws Exception unit_unit(5)

zero_unit = RBT.UnitPipe(Zero)
println(repr(zero_unit))

@test_throws Exception zero_unit()

unit_zero = RBT.ZeroPipe(Unit)
println(repr(unit_zero))

@test_throws Exception unit_zero()

const5 = RBT.ConstPipe(5)
println(repr(const5))

@test const5() == 5
@test const5(nothing) == 5
@test_throws Exception const5(5)

null = RBT.NullPipe()
println(repr(null))

@test isequal(null(), Nullable())
@test_throws Exception null(5)

empty = RBT.EmptyPipe()
println(repr(empty))

@test empty() == Union{}[]
@test_throws Exception empty(5)

primes = RBT.SetPipe(:primes, Int[2,3,5,7,11], ismonic=true)
println(repr(primes))

@test primes() == Int[2,3,5,7,11]

sqmap = RBT.IsoMapPipe(:sq, Dict(1=>1,2=>4,3=>9,4=>16,5=>25,6=>36,7=>49,8=>64,9=>81,10=>100,11=>121), true, false)
println(repr(sqmap))

@test sqmap(5) == 25

decmap = RBT.OptMapPipe(:dec, Dict(2=>1,3=>2,4=>3,5=>4,6=>5,7=>6,8=>7,9=>8), true, false)
println(repr(decmap))

@test isequal(decmap(5), Nullable(4))
@test isequal(decmap(1), Nullable{Int}())

divmap = RBT.SeqMapPipe(:div, Dict(4=>[2],6=>[2,3],8=>[2,4],9=>[3],10=>[2,5]), false, false, false)
println(repr(divmap))

@test divmap(6) == Int[2,3]
@test divmap(3) == Int[]

println(repr(primes >> sqmap))
println(repr(primes >> decmap))
println(repr(primes >> divmap))

@test (primes >> sqmap)() == [4,9,25,49,121]
@test (primes >> decmap)() == [1,2,4,6]
@test (primes >> divmap)() == []
@test (divmap >> decmap)(6) == [1,2]
@test isequal((decmap >> decmap >> decmap)(4), Nullable(1))

tuple0 = RBT.TuplePipe()
println(repr(tuple0))

@test tuple0() == ()
@test tuple0(5) == ()

tuple3a = RBT.TuplePipe(sqmap, decmap, divmap)
println(repr(tuple3a))

@test isequal(tuple3a(6), (36, Nullable(5), [2,3]))
@test isequal(tuple3a(11), (121, Nullable{Int}(), Int[]))

tuple3a_items = RBT.TuplePipe(RBT.ItemPipe(tuple3a, 3), RBT.ItemPipe(tuple3a, 2), RBT.ItemPipe(tuple3a, 1))
println(repr(tuple3a_items))

@test isequal((tuple3a >> tuple3a_items)(6), ([2,3], Nullable(5), 36))
@test isequal((tuple3a >> tuple3a_items)(11), (Int[], Nullable{Int}(), 121))

tuple3b = RBT.TuplePipe(int_here, sqmap, decmap >> divmap)
println(repr(tuple3b))

@test tuple3b(7) == (7, 49, [2,3])

vector0 = RBT.VectorPipe()
println(repr(vector0))

@test vector0() == []

vector1a = RBT.VectorPipe(sqmap)
println(repr(vector1a))

@test vector1a(6) == [36]

vector1b = RBT.VectorPipe(decmap)
println(repr(vector1b))

@test vector1b(6) == [5]
@test vector1b(1) == []

vector1c = RBT.VectorPipe(divmap)
println(repr(vector1c))

@test vector1c === divmap
@test vector1c(6) == [2,3]
@test vector1c(1) == []

vector4 = RBT.VectorPipe(int_here, sqmap, decmap, divmap)
println(repr(vector4))

@test vector4(6) == [6,36,5,2,3]
@test vector4(1) == [1,1]

@test (primes >> vector4)() == [2,4,1,3,9,2,5,25,4,7,49,6,11,121]

@test RBT.IsoPipe(sqmap) === sqmap

decmap_as_iso = RBT.IsoPipe(decmap)
println(repr(decmap_as_iso))

@test decmap_as_iso(5) == 4
@test_throws Exception decmap_as_iso(1)

divmap_as_iso = RBT.IsoPipe(divmap)
println(repr(divmap_as_iso))

@test divmap_as_iso(4) == 2
@test_throws Exception divmap_as_iso(1)
@test_throws Exception divmap_as_iso(6)

@test RBT.OptPipe(decmap) === decmap

sqmap_as_opt = RBT.OptPipe(sqmap)
println(repr(sqmap_as_opt))

@test isequal(sqmap_as_opt(5), Nullable(25))

divmap_as_opt = RBT.OptPipe(divmap)
println(repr(divmap_as_opt))

@test isequal(divmap_as_opt(4), Nullable(2))
@test isequal(divmap_as_opt(1), Nullable{Int}())
@test_throws Exception divmap_as_opt(6)

@test RBT.SeqPipe(divmap) === divmap

sqmap_as_seq = RBT.SeqPipe(sqmap)
println(repr(sqmap_as_seq))

@test sqmap_as_seq(5) == [25]

decmap_as_seq = RBT.SeqPipe(decmap)
println(repr(decmap_as_seq))

@test decmap_as_seq(5) == [4]
@test decmap_as_seq(1) == Int[]

id_eq_sq = int_here .== sqmap
println(repr(id_eq_sq))

@test id_eq_sq(1) == true
@test id_eq_sq(5) == false

dec_dec_sq_ne_id = (decmap >> decmap >> sqmap) .!= int_here
println(repr(dec_dec_sq_ne_id))

@test isequal(dec_dec_sq_ne_id(4), Nullable(false))
@test isequal(dec_dec_sq_ne_id(6), Nullable(true))
@test isequal(dec_dec_sq_ne_id(2), Nullable{Bool}())

dec_dec_le_div = (decmap >> decmap) .<= divmap
println(repr(dec_dec_le_div))

@test dec_dec_le_div(4) == [true]
@test dec_dec_le_div(6) == [false, false]
@test dec_dec_le_div(10) == Bool[]

between_3_and_6 = (RBT.ConstPipe(Int, 3) .<= int_here) & (int_here .<= RBT.ConstPipe(Int, 6))
println(repr(between_3_and_6))

@test between_3_and_6(5) == true
@test between_3_and_6(1) == false

not_0_or_1 = ~((int_here .== RBT.ConstPipe(Int, 0)) | (int_here .== RBT.ConstPipe(Int, 1)))
println(repr(not_0_or_1))

@test not_0_or_1(1) == false
@test not_0_or_1(5) == true

any_div_between_3_and_6 = RBT.AnyPipe(divmap >> between_3_and_6)
println(repr(any_div_between_3_and_6))

@test any_div_between_3_and_6(1) == false
@test any_div_between_3_and_6(4) == false
@test any_div_between_3_and_6(6) == true
@test any_div_between_3_and_6(10) == true

all_div_between_3_and_6 = RBT.AllPipe(divmap >> between_3_and_6)
println(repr(all_div_between_3_and_6))

@test all_div_between_3_and_6(1) == true
@test all_div_between_3_and_6(4) == false
@test all_div_between_3_and_6(6) == false
@test all_div_between_3_and_6(9) == true

count_primes = RBT.CountPipe(primes)
println(repr(count_primes))

@test count_primes() == 5

max_prime = RBT.IntMaxPipe(primes)
println(repr(max_prime))

@test isequal(max_prime(), Nullable(11))

max_prime_or_0 = RBT.IntMaxPipe(RBT.VectorPipe(primes, RBT.ConstPipe(0)))
println(repr(max_prime_or_0))

@test max_prime_or_0() == 11

min_prime = RBT.IntMinPipe(primes)
println(repr(min_prime))

@test isequal(min_prime(), Nullable(2))

sum_primes = RBT.IntSumPipe(primes)
println(repr(sum_primes))

@test sum_primes() == 28

mean_prime = RBT.IntMeanPipe(primes)
println(repr(mean_prime))

@test isequal(mean_prime(), Nullable(28/5))

first_prime = RBT.FirstPipe(primes)
println(repr(first_prime))

@test isequal(first_prime(), Nullable(2))

last_prime = RBT.FirstPipe(primes, true)
println(repr(last_prime))

@test isequal(last_prime(), Nullable(11))

mixed = sqmap * decmap * divmap
println(repr(mixed))

@test mixed(6) == [(36,5,2), (36,5,3)]

range1 = RBT.RangePipe(RBT.ConstPipe(1), RBT.ConstPipe(10))
println(repr(range1))

@test range1() == [1,2,3,4,5,6,7,8,9,10]

range2 = RBT.RangePipe(RBT.ConstPipe(1), RBT.ConstPipe(2), RBT.ConstPipe(10))
println(repr(range2))

@test range2() == [1,3,5,7,9]

tagged = RBT.PackPipe(:sq => sqmap, :dec => decmap, :div => divmap)
println(repr(tagged))

@test tagged(6) == Pair{Symbol,Int}[:sq => 36, :dec => 5, :div => 2, :div => 3]

untagged = RBT.CasePipe(
    Int, false, :dec => int_here, :div => decmap)
println(repr(untagged))

@test (tagged >> untagged)(6) == [5,1,2]

primes_filter_select =
    RBT.SelectPipe(
        RBT.FilterPipe(
            primes,
            (int_here .>= RBT.ConstPipe(Int, 5)) & (int_here .<= RBT.ConstPipe(Int, 7))),
        int_here, decmap >> divmap)
println(repr(primes_filter_select))

@test primes_filter_select() == [((5, [2]), 5), ((7, [2,3]), 7)]

rev_primes = RBT.SortPipe(primes, true)
println(repr(rev_primes))

@test rev_primes() == [11,7,5,3,2]

primes_by_count_dec_div = RBT.SortByPipe(primes, RBT.CountPipe(decmap >> divmap))
println(repr(primes_by_count_dec_div))

@test primes_by_count_dec_div() == [2,3,11,5,7]

first_prime_by_count_dec_div = RBT.FirstByPipe(primes, RBT.CountPipe(decmap >> divmap))
println(repr(first_prime_by_count_dec_div))

@test isequal(first_prime_by_count_dec_div(), Nullable(7))

first_3_primes = RBT.TakePipe(primes, RBT.ConstPipe(3))
println(repr(first_3_primes))

@test first_3_primes() == [2,3,5]

last_2_primes = RBT.TakePipe(primes, RBT.ConstPipe(3), true)
println(repr(last_2_primes))

@test last_2_primes() == [7,11]

reversed_primes = RBT.ReversePipe(primes)
println(repr(reversed_primes))

@test reversed_primes() == [11,7,5,3,2]

get_prime_5 = RBT.GetPipe(primes, RBT.ConstPipe(5))
println(repr(get_prime_5))

@test isequal(get_prime_5(), Nullable(5))

divs = RBT.ConnectPipe(divmap)
println(repr(divs))

@test divs(8) == [2,4,2]

self_divs = RBT.ConnectPipe(divmap, true)
println(repr(self_divs))

@test self_divs(8) == [8,2,4,2]

depths = range1 >> RBT.DepthPipe(divmap)
@test depths() == [0,0,0,1,0,1,0,2,1,1]

range_sorted_by_div = RBT.SortConnectPipe(range1, divmap)
println(repr(range_sorted_by_div))

@test range_sorted_by_div() == [1,2,4,8,3,6,9,5,10,7]

unique_div = RBT.UniquePipe(range1 >> divmap, true)
println(repr(unique_div))

@test unique_div() == [5,4,3,2]

range1_by_count_div = RBT.GroupPipe(range1, RBT.CountPipe(divmap))
println(repr(range1_by_count_div))

@test range1_by_count_div() == [((0,),[1,2,3,5,7]), ((1,),[4,9]), ((2,),[6,8,10])]

range1_partition_by_count_div =
    RBT.PartitionPipe(
        range1,
        true,
        (RBT.RangePipe(RBT.ConstPipe(-1), RBT.ConstPipe(1)), RBT.CountPipe(divmap), true))
println(repr(range1_partition_by_count_div))

@test isequal(range1_partition_by_count_div(),
    [((Nullable(1),),[4,9]),
     ((Nullable(0),),[1,2,3,5,7]),
     ((Nullable(-1),),[]),
     ((Nullable{Int}(),),[1,2,3,4,5,7,9])])

dict4 = RBT.DictPipe(
    :here => int_here,
    :sq => sqmap,
    :dec => RBT.NullToVoidPipe(decmap),
    :div => divmap)
println(repr(dict4))

@test dict4(6) == Dict{Any,Any}(:here => 6, :sq => 36, :dec => 5, :div => [2,3])
@test dict4(1) == Dict{Any,Any}(:here => 1, :sq => 1, :dec => nothing, :div => [])

df =
    RBT.DataFramePipe(
        RBT.RangePipe(RBT.ConstPipe(1), RBT.ConstPipe(4)),
        :here => int_here, :dec => decmap, :div => divmap)
println(repr(df))

@test isequal(df(), DataFrame(Any[[1,2,3,4], @data([NA,1,2,3]), Vector{Int}[[],[],[],[2]]], [:here,:dec,:div]))

sqr = RBT.FilterPipe(range1, sqmap .== RBT.ParamPipe(Int, :X, Int))
println(repr(sqr))

@test sqr(X=36) == [6]

in_div = RBT.FilterPipe(range1, RBT.InPipe(RBT.ParamPipe(Int, :X, Int), divmap))
println(repr(in_div))

@test in_div(X=2) == [4,6,8,10]

top_div =
    RBT.FilterPipe(
        range1,
        RBT.CountPipe(divmap) .==
            RBT.IntMaxPipe(RBT.RelativePipe(Int, true, true, true) >> RBT.CountPipe(divmap)))
println(repr(top_div))

@test top_div() == [6,8,10]

max_by_div =
    RBT.FilterPipe(
        range1,
        int_here .> RBT.IntMaxPipe(RBT.RelativeByPipe(RBT.CountPipe(divmap), true, false, true)))
println(repr(max_by_div))

@test max_by_div() == [7,9,10]

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using KNITRO
using Compat.Test

#    min  9 - 8x1 - 6x2 - 4x3
#         + 2(x1^2) + 2(x2^2) + (x3^2) + 2(x1*x2) + 2(x1*x3)
#    subject to  c[0]:  x1 + x2 + 2x3 <= 3
#                x1 >= 0
#                x2 >= 0
#                x3 >= 0
#    initpt (0.5, 0.5, 0.5)
#
#    Solution is x1=4/3, x2=7/9, x3=4/9, lambda=2/9  (f* = 1/9)
#
#  The problem comes from Hock and Schittkowski, HS35.

function eval_f(x::Vector{Float64})
    linear_terms = 9.0 - 8.0*x[1] - 6.0*x[2] - 4.0*x[3]
    quad_terms = 2.0*x[1]^2 + 2.0*x[2]^2 + x[3]^2 + 2.0*x[1]*x[2] + 2.0*x[1]*x[3]
    return linear_terms + quad_terms
end

function eval_g(x::Vector{Float64}, cons::Vector{Float64})
    cons[1] = x[1] + x[2] + 2.0*x[3]
end

function eval_grad_f(x::Vector{Float64}, grad::Vector{Float64})
    grad[1] = -8.0 + 4.0*x[1] + 2.0*x[2] + 2.0*x[3]
    grad[2] = -6.0 + 2.0*x[1] + 4.0*x[2]
    grad[3] = -4.0 + 2.0*x[1]            + 2.0*x[3]
end

function eval_jac_g(x::Vector{Float64}, jac::Vector{Float64})
    jac[1] = 1.0
    jac[2] = 1.0
    jac[3] = 2.0
end

objGoal = KTR_OBJGOAL_MINIMIZE
objType = KTR_OBJTYPE_QUADRATIC

n = 3
x_L = zeros(n)
x_U = [KTR_INFBOUND,KTR_INFBOUND,
KTR_INFBOUND]

m = 1
c_Type = [KTR_CONTYPE_LINEAR]
c_L = [-KTR_INFBOUND]
c_U = [3.0]

jac_con = Int32[0,0,0]
jac_var = Int32[0,1,2]

x       = [0.5,0.5,0.5]
lambda  = zeros(n+m)
obj     = [0.0]

kp = createProblem()
setOption(kp, "hessopt", 6)
initializeProblem(kp, objGoal, objType, x_L, x_U, c_Type, c_L, c_U, jac_var, jac_con)
setFuncCallback(kp, eval_f, eval_g)
setGradCallback(kp, eval_grad_f, eval_jac_g)
solveProblem(kp)

# --- test optimal solutions ---
@testset "Test optimal solution" begin
    @test applicationReturnStatus(kp) == :Optimal
    @test isapprox(kp.x, [1.3333333, 0.7777777, 0.4444444], atol=1e-5)
    @test isapprox(kp.obj_val[1], 0.1111111, atol=1e-5)
end

freeProblem(kp)


text:
using BenchmarkTools

# Threads.threadid() tells us the ID number of the thread currently doing a task.
function thread_test()
    Threads.@threads for i in 1:6
        id = Threads.threadid()
        println("Thread ID: ", id, " | i: ", i)
    end
end

# thread_test()

function serialsum(xs)
    r = 0
    for i in xs
        r += i
    end
    return r
end

function threadsum(xs)
    r = zeros(Threads.nthreads())
    Threads.@threads for i in xs
        # Each thread has it's own "sum bucket" to accumulate
        # bits and pieces of xs -- we'll add these buckets together
        # at the end.
        r[Threads.threadid()] += i
    end
    return sum(r)
end

values = randn(10_000, 10_000);
# values = randn(100);

# @time serialsum(values)
# @time threadsum(values)

xs = []
Threads.@threads for i in 1:1000
    push!(xs, i) # Add i to the end of xs
end

println(length(xs))


text:
scipy_optimize = pyimport("scipy.optimize")

function prediction_error(
    mountaincar::MountainCar,
    params::Array{Float64},
    data::Array{MountainCarContTransition},
)
    err = []
    for i = 1:length(data)
        transition = data[i]
        predicted_state =
            step(mountaincar, transition.initial_state, transition.action, params)[1]
        # push!(err, norm(vec(predicted_state) - vec(transition.final_state)))
        push!(err, (predicted_state.position -
                    transition.final_state.position)^2 + (predicted_state.speed -
                                                          transition.final_state.speed)^2)

    end
    err
end

function get_least_squares_fit(
    mountaincar::MountainCar,
    params::MountainCarParameters,
    data::Array{MountainCarContTransition},
)
    fn(p) = prediction_error(mountaincar, p, data)
    scipy_optimize.leastsq(fn, vec(params))[1]
end

function get_least_squares_fit(
    mountaincar::MountainCar,
    params::Array{Float64},
    data::Array{MountainCarContTransition},
)
    get_least_squares_fit(
        mountaincar,
        unvec_params(params),
        data,
    )
end


text:
function bernoulli_iter_redundant(sys, binary_vars::Array{Int,1})

    redundant_iter_sub = Dict()

    # starting from N+2 index as 0th and 1st order terms stay the same
    for iter in sys.iter_all[sys.N+2:end]

        iter_temp = [iter...] # need to transform into array as Tuple is immutable
        for ind in binary_vars
            if iter[ind] > 1
                iter_temp[ind] = 1
            end
        end
        iter_temp = Tuple(iter_temp)

        if iter != iter_temp
            redundant_iter_sub[iter] = iter_temp
        end

    end

    redundant_iter_sub

end

function bernoulli_reduce(sys::MomentEquations, binary_vars::Array{Int,1})

    N = sys.N

    redundant_iter_sub = bernoulli_iter_redundant(sys, binary_vars)
    redundant_iter = keys(redundant_iter_sub)

    redundant_eqs = Int[]
    for ind in binary_vars
        for (i, iter) in enumerate(sys.iter_m)
            if iter[ind] > 1
                push!(redundant_eqs, i+N)
            end
        end
    end

    if sys isa RawMomentEquations
        μ = copy(sys.μ)
        iter_sub = [μ[key] => μ[val] for (key, val) in redundant_iter_sub]
    else
        μ = define_μ(N, sys.q_order)
        μ_redundant_sub = [μ[key] => μ[val] for (key, val) in redundant_iter_sub]

        clean_iter = setdiff(sys.iter_all[N+2:end], redundant_iter)
        central_to_raw = central_to_raw_moments(N, sys.q_order)
        μ_clean_sub = Dict([μ[iter] => central_to_raw[iter] for iter in clean_iter])

        raw_to_central = raw_to_central_moments(N, sys.q_order)
        iter_sub = Dict()
        for iter in redundant_iter
            M_temp = raw_to_central[iter]
            M_temp = simplify(substitute(M_temp, μ_redundant_sub))
            M_temp = substitute(M_temp, μ_clean_sub)
            M_temp = simplify(M_temp)
            iter_sub[sys.M[iter]] = M_temp
        end

    end

    redundant_iter, redundant_eqs, iter_sub

end

"""
    bernoulli_moment_eqs(sys::MomentEquations, binary_vars::Array{Int,1})

Given `MomentEquations` and an array of indices specifying the species which molecule numbers are
binary variables (either 0 or 1), apply identities of Bernoulli variables to remove the redundant
ODEs and return the *cleaned up* `MomentEquations`. See [here](@ref geometric-and-conditional) for
example usage.
"""
function bernoulli_moment_eqs(sys::MomentEquations, binary_vars::Array{Int,1})

    # Construct moment equations removing the redundant ones
    # noting the properties of the Bernoulli variables in the system

    redundant_iter, redundant_eqs, iter_sub = bernoulli_reduce(sys, binary_vars)

    # construct the cleaned moment equations
    clean_eqs = Equation[]
    for (i, eq) in enumerate(sys.odes.eqs)
        if !(i in redundant_eqs)
            clean_rhs = substitute(eq.rhs, iter_sub)
            clean_rhs = expand(clean_rhs)
            push!(clean_eqs, Equation(eq.lhs, clean_rhs))
        end
    end

    # construct a new Raw/Central/MomentEquations system saving only the clean iterators
    clean_iter = setdiff(sys.iter_all, redundant_iter)
    iter_m = filter(x -> 2 <= sum(x) <= sys.m_order, clean_iter)
    iter_q = filter(x -> sys.m_order < sum(x) <= sys.q_order, clean_iter)

    field_values = [getfield(sys, field) for field in fieldnames(typeof(sys))]
    ind_iter_m = findall(x -> x==:iter_m, fieldnames(typeof(sys)))[1]
    ind_iter_q = findall(x -> x==:iter_q, fieldnames(typeof(sys)))[1]
    ind_iter_all = findall(x -> x==:iter_all, fieldnames(typeof(sys)))[1]

    field_values[ind_iter_m] = iter_m       #sys.iter_m
    field_values[ind_iter_q] = iter_q       #sys.iter_q
    field_values[ind_iter_all] = clean_iter #sys.iter_all

    ## fixing ODE system to preserve consistent ordering of parameters
    iv = get_iv(sys.odes)
    ps = get_ps(sys.odes)

    if sys isa RawMomentEquations
        vars = extract_variables(clean_eqs, sys.μ)
    else
        vars = extract_variables(clean_eqs, sys.μ, sys.M)
    end

    odename = Symbol(nameof(sys.odes), "_bernoulli")
    odes = ODESystem(clean_eqs, iv, vars, ps; name=odename)

    new_system = typeof(sys)(odes, field_values[2:end]...)

    new_system

end


text:
##
##  i n t e g r a t i o n . j l  Integration Routines
##

function trapz(x::Vector{Tx}, y::Vector{Ty}) where {Tx<:Number, Ty<:Number}
    # Trapezoidal integration rule
    local n = length(x)
    if (length(y) != n)
        error("Vectors 'x', 'y' must be of same length")
    end
    r = zero(zero(Tx) + zero(Ty))
    if n == 1; return r; end
    for i in 2:n
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    #= correction -h^2/12 * (f'(b) - f'(a))
    ha = x[2] - x[1]
    he = x[end] - x[end-1]
    ra = (y[2] - y[1]) / ha
    re = (y[end] - y[end-1]) / he
    r/2 - ha*he/12 * (re - ra)
    =#
    return r/2
end

function romberg(f::Function, a::Real, b::Real;
                   maxit = 50, tol = 1.0e-15)
    if a == b; return 0.0, 0.0; end
    sgn = 1.0
    if a > b
        a, b = b, a
        sgn = -1.0
    end

    local n::Int = 1; local iter::Int = 0
    I = zeros(maxit+1, maxit+1)

    err = 1.0
    while err > tol && iter < maxit
        iter += 1
        n *= 2
        h = (b - a)/ n
        S = f(a)
        for i = 1:(n-1)
            xi = a + h*i
            S += 2.0*f(xi)
        end
        I[iter+1, 1] = (S + f(b)) * h / 2.0
        for k = 2:(iter+1)
            j = 2 + iter - k
            I[j,k] = (4^(k-1)*I[j+1,k-1] - I[j,k-1]) / (4^(k-1)-1)
        end
        err = abs((I[1,iter+1] - I[2,iter]) / I[1,iter+1])
     end
     if iter >= maxit
         warn("Maximum number of iterations has been reached.")
     end
     if err == 0.0; err = tol; end

     return I[1, iter+1], err
end


text:
using Random
using Distributions
using MultiAgentAllocationTransit
using BenchmarkTools
using Statistics
using JLD2
using Logging
global_logger(SimpleLogger(stderr, Logging.Error))

rng = MersenneTwister(6789)

# EXAMPLES OF ARGUMENTS
# const bb_params_file = "./data/sfmta/sf_bb_params.toml"
# const stop_coords_file = "./data/sfmta/stop_to_coords.json"
# const trips_file = "./data/sfmta/trips.json"
# const drone_params_file = "./data/drone_params.toml"
# const city_tt_outfile = "./data/sfmta/sf_halton_travel_estimates.jld2"
# const n_halton_points = 100

const bb_params_file = ARGS[1]
const stop_coords_file = ARGS[2]
const trips_file = ARGS[3]
const drone_params_file = ARGS[4]
const city_tt_outfile = ARGS[5]
const n_halton_points = parse(Int64, ARGS[6])

const TRANSIT_CAP_RANGE = (3, 5)
const ECBS_WEIGHT = 1.1

function main()
    bb_params = parse_city_params(bb_params_file)
    drone_params = parse_drone_params(drone_params_file)

    tg = load_transit_graph_latlong(stop_coords_file, trips_file, TRANSIT_CAP_RANGE, rng)
    tg, stop_idx_to_trips =
                    transit_graph_preprocessing(tg, MultiAgentAllocationTransit.distance_lat_lon_euclidean, drone_params)

    # Generate halton NN of city
    halton_nn_tree, city_halton_points = generate_city_halton_nn(bb_params; n_points = n_halton_points)

    # Now create a bunch of fake depots and sites
    midway = n_halton_points / 2
    depots = [LatLonCoords((lat = chp[1], lon = chp[2])) for chp in city_halton_points]
    # sites = [LatLonCoords((lat = chp[1], lon = chp[2])) for chp in city_halton_points[midway+1:end]]

    # Create OTG
    dummy_otg = OffTransitGraph(depots = depots, sites = LatLonCoords[])
    state_graph, depot_sites_to_vtx, trip_to_vtx_range = setup_state_graph(tg, dummy_otg)

    env = MAPFTransitEnv(off_transit_graph = dummy_otg, transit_graph = tg, state_graph = state_graph,
                         agent_states = AgentState[], depot_sites_to_vtx = depot_sites_to_vtx, trip_to_vtx_range = trip_to_vtx_range,
                         stop_idx_to_trips = stop_idx_to_trips, aug_trips_fws_dists = zeros(1, 1),
                         drone_params = drone_params, dist_fn = MultiAgentAllocationTransit.distance_lat_lon_euclidean,
                         curr_site_points = [], threshold_global_conflicts = 5)

    println("Computing estimates")
    @time travel_time_estimates = compute_all_pairs_estimates(env, n_halton_points, ECBS_WEIGHT)

    # Save to file
    @save city_tt_outfile halton_nn_tree city_halton_points travel_time_estimates
end

main()


text:
######################################################################
# Suite of macros to automatically generate getters and setters for
# object properties along with the accompanying `hasproperty` and
# `propertynames` specializations.
# -----
# Licensed under MIT License
module GenerateProperties

const Optional{T} = Union{T, Nothing}
const GetterSetterBody = Tuple{Optional{LineNumberNode}, Expr}

struct StructProp{S, F} end
struct StructField{S, F} end

assign(inst, field::Symbol, value) = assign(StructProp{typeof(inst), field}, inst, value)
assign( ::Type{StructProp{ S, F}}, inst::S, value) where {S, F} = assign(StructField{S, F}, inst, value)
assign(T::Type{StructField{S, F}}, inst::S, value) where {S, F} = assign(getfieldtype(T), inst, F, value)
assign(::Type{T},  inst, field::Symbol, value)     where {T}                      = setfield!(inst, field, convert(T, value))
assign(::Type{T1}, inst, field::Symbol, value::T2) where {T1, T2<:T1}             = setfield!(inst, field, value)
assign(::Type{T1}, inst, field::Symbol, value::T2) where {T1<:Number, T2<:Number} = setfield!(inst, field, T1(value))

retrieve(inst, field::Symbol) = retrieve(StructProp{typeof(inst), field}, inst)
retrieve(::Type{StructProp{ S, F}}, inst::S) where {S, F} = retrieve(StructField{S, F}, inst)
retrieve(::Type{StructField{S, F}}, inst::S) where {S, F} = getfield(inst, F)

getparambody(x) = x
getparambody(unionall::UnionAll) = getparambody(unionall.body)

getfieldtype(S::Type, F::Symbol) = getfieldtype(StructField{S, F})
@generated function getfieldtype(::Type{StructField{S, F}}) where {S, F}
    @assert !isa(S, UnionAll) "Cannot reliably use getfieldtype with UnionAll"

    idx = findfirst(field->field==F, fieldnames(S))
    @assert idx !== nothing "Field $F not found in type $S"

    T = S.types[idx]
    :($T)
end

export @generate_properties
macro generate_properties(T, block)
    if !isa(block, Expr) || block.head != :block
        throw(ArgumentError("Second argument to @generate_properties must be a block"))
    end

    result = Expr(:block)
    props  = Set{Symbol}()
    symget = Symbol("@get")
    symset = Symbol("@set")
    symeq  = Symbol("=")

    lastlinenumber = nothing
    for expr ∈ block.args
        if isa(expr, LineNumberNode)
            lastlinenumber = expr
        else
            if expr.head != :macrocall || expr.args[1] ∉ (symget, symset)
                throw(ArgumentError("Every line must be a call to either @get or @set"))
            end

            args = filterlinenumbers(expr.args)
            if args[2].head != symeq throw(ArgumentError("Getter/Setter not an assignment")) end
            prop, body = filterlinenumbers(args[2].args)
            body = replace_self(T, body)
            push!(props, prop)

            if expr.args[1] == symget
                push!(result.args, generate_getter(T, prop, lastlinenumber, body))
            elseif expr.args[1] == symset
                push!(result.args, generate_setter(T, prop, lastlinenumber, body))
            end

            lastlinenumber = nothing
        end
    end

    # Generate propertynames
    push!(result.args, quote
        @generated function Base.propertynames(::$T)
            res = tuple(union($props, fieldnames($T))...)
            :($res)
        end
    end)

    push!(result.args, :(Base.getproperty(self::$T, prop::Symbol) = GenerateProperties.retrieve(self, prop)))
    push!(result.args, :(Base.setproperty!(self::$T, prop::Symbol, value) = GenerateProperties.assign(self, prop, value)))

    esc(result)
end

export @get, @set
macro get(args...) end
macro set(args...) end

filterlinenumbers(exprs) = filter(expr->!isa(expr, LineNumberNode), exprs)

replace_self(_, expr) = expr
function replace_self(T::Union{Symbol, Expr}, expr::Expr)
    @assert isa(T, Symbol) || T.head == :curly
    if expr.head == :(=)
        lhs, rhs = expr.args

        if isa(lhs, Expr) && lhs.head == :. && lhs.args[1] == :self
            prop = lhs.args[2]::QuoteNode
            expr = :(GenerateProperties.assign(self, $rhs))
            insert!(expr.args, 2, structfieldexpr(prop))
        end
    elseif expr.head == :.
        if expr.args[1] == :self
            prop = expr.args[2]::QuoteNode
            expr = :(GenerateProperties.retrieve(self))
            insert!(expr.args, 2, structfieldexpr(prop))
        end
    end
    expr.args = map(sub->replace_self(T, sub), expr.args)
    expr
end

function generate_getter(T, prop::Symbol, linenumber::Optional{LineNumberNode}, body)
    block = Expr(:block)
    if linenumber !== nothing push!(block.args, linenumber) end
    push!(block.args, body)

    fnexpr = :(GenerateProperties.retrieve(self::$T) = $block)
    insert!(fnexpr.args[1].args, 2, argsubtypeexpr(structpropexpr(T, prop)))
    fnexpr
end

function generate_setter(T, prop::Symbol, linenumber::Optional{LineNumberNode}, body)
    block = Expr(:block)
    if linenumber !== nothing push!(block.args, linenumber) end
    push!(block.args, body)

    fnexpr = :(GenerateProperties.assign(self::$T, value) = $block)
    insert!(fnexpr.args[1].args, 2, argsubtypeexpr(structpropexpr(T, prop)))
    fnexpr
end

structpropexpr(T, prop::QuoteNode) = Expr(:curly, :(GenerateProperties.StructProp), :(<:$T), prop)
structpropexpr(T, prop::Symbol)    = structpropexpr(T, QuoteNode(prop))
structfieldexpr(prop::QuoteNode) = Expr(:curly, :(GenerateProperties.StructField), :(typeof(self)), prop)
structfieldexpr(prop::Symbol)    = structfieldexpr(QuoteNode(prop))

function argtypeexpr(type::Expr)
    Expr(:(::), Expr(:curly, :Type, type))
end

function argsubtypeexpr(type::Expr)
    Expr(:(::), Expr(:curly, :Type, Expr(:<:, type)))
end

end # module GenerateProperties


text:
# greedy spectral embedding

# packages
using PyPlot
using Debug

# local
using KPCA

# simple data set:
dim = 2
grouped = 10
n = 4*grouped
n2 = floor(n/2)

#regenerate = false
regenerate = true
if regenerate
	# square with 2 pts on each corner - dictionary should be approximately one from each corner
	X = randn(n,dim) .* 1
	L = zeros(n,1)

	i = 1+0*grouped
	X[i:i+grouped-1,:] += 20
	L[i:i+grouped-1] = 1

	i = 1+1*grouped
	X[i:i+grouped-1,:] -= 20
	L[i:i+grouped-1] = 2

	i = 1+2*grouped
	X[i:i+grouped-1,1] += 5
	X[i:i+grouped-1,2] -= 10
	L[i:i+grouped-1] = 3

	i = 1+3*grouped
	X[i:i+grouped-1,1] -= 5
	X[i:i+grouped-1,2] += 10
	L[i:i+grouped-1] = 4
end

if false #regenerate

	# two gaussian groups of pts similar to figure in paper
	X = randn(n,dim)
	X[1:n2,:] .*= 10
	X[n2:n,1] .+= 20
	X[n2:n,1] .*= 10

	L = ones(n,1)
	L[n2:n] = 2

	figure()
	maxv = maximum(X)
	minv = minimum(X)
	p_i = 1
	for r=1:dim
			for c=1:dim
					subplot(dim,dim,p_i)
					scatter(X[:,c],X[:,r])
					xlim(minv,maxv)
					ylim(minv,maxv)

					p_i += 1
			end
	end

end

for i=1:5 # mix up well
	newind = shuffle([1:n])
	X = X[newind,:]
	L = L[newind]
end

println(mean(X,1))

K = LinearKernel(X)

# test definitions - separated into function to make code isolation easier..

function test_kpca(K)
	# run KPCA to compare:
	Kc, P_kpca, Y_kpca = KernelPCA(K, 0, 2, "none")

	kpca_proj_err = 0
	for i=1:size(K,1)
		# kpca
		ki = K[:,i]
		kpca_proj_err += K[i,i] - (ki' * (K \ ki) )

	end
	# avg proj err
	kpca_proj_err /= size(K,1)
	println("kpca proj err: ", kpca_proj_err)
end

# run the tests:
test_kpca(K)


text:
"""
    SparseMatrixCSR{Bi,Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}

Matrix type for storing sparse matrices
in the Compressed Sparse Row format with `Bi`-based indexing (typically 0 or 1).
The standard way of constructing `SparseMatrixCSR` is through the
[`sparsecsr`](@ref) function.

# Properties

- `m::Int` Number of columns
- `n::Int` Number of rows
- `rowptr::Vector{Ti}` Row `i` (1-based) in `nzval` vector at indices (1-based) `(rowptr[i]+(1-Bi)):(rowptr[i+1]+(1-Bi)-1)`
- `colval::Vector{Ti}` Col indices (`Bi`-based) of stored values
- `nzval::Vector{Tv}`  Stored values, typically non zeros

# Inner constructor

    SparseMatrixCSR{Bi}(
      m::Integer,
      n::Integer,
      rowptr::Vector{Ti},
      colval::Vector{Ti},
      nzval::Vector{Tv}) where {Bi,Tv,Ti<:Integer}

"""
struct SparseMatrixCSR{Bi,Tv,Ti} <: AbstractSparseMatrix{Tv,Ti}
  m::Int
  n::Int
  rowptr::Vector{Ti}
  colval::Vector{Ti}
  nzval::Vector{Tv}
  function SparseMatrixCSR{Bi}(m::Integer, n::Integer, rowptr::Vector{Ti}, colval::Vector{Ti},
    nzval::Vector{Tv}) where {Bi,Tv,Ti<:Integer}
    @noinline throwsz(str, k, r) =
      throw(ArgumentError("$str must be $r, got $k"))
    m < 0 && throwsz("m",">=0",m)
    n < 0 && throwsz("n",">=0", n)
    length(rowptr) != m+1 && throwsz("lengh(rowptr)",m+1,length(rowptr))
    rowptr[end]-Bi != length(nzval) && throwsz("rowptr[end]-Bi",rowptr[m+1]-Bi, length(nzval))
    new{Bi,Tv,Ti}(Int(m), Int(n), rowptr, colval, nzval)
  end
end

"""
    SparseMatrixCSR(a::Transpose{Tv,<:SparseMatrixCSC} where Tv)

Build a 1-based `SparseMatrixCSR` from the lazy transpose of a `SparseMatrixCSC`.
The resulting matrix takes ownership of the internal storage of input matrix.
Modifying the values of one, will mutate also the other.
"""
function SparseMatrixCSR(a::Transpose{Tv,<:SparseMatrixCSC} where Tv)
  at = a.parent
  SparseMatrixCSR{1}(size(a,1),size(a,2),at.colptr,rowvals(at),nonzeros(at))
end

"""
    SparseMatrixCSR{Bi}(a::Transpose{Tv,<:SparseMatrixCSC} where Tv) where Bi

Build a `Bi`-based `SparseMatrixCSR` from the lazy transpose of a `SparseMatrixCSC`.
The resulting matrix takes ownership of the internal storage of input matrix
for any `Bi` and modifies the internal storage when `Bi != 1`.
The input matrix will become unusable in the latter case.
"""
function SparseMatrixCSR{Bi}(a::Transpose{Tv,<:SparseMatrixCSC} where Tv) where Bi
  at = a.parent
  at.colptr .-= getoffset(Bi)
  rowvals(at) .-= getoffset(Bi)
  SparseMatrixCSR{Bi}(size(a,1),size(a,2),at.colptr,rowvals(at),nonzeros(at))
end
SparseMatrixCSR{1}(a::Transpose{Tv,<:SparseMatrixCSC} where Tv) = SparseMatrixCSR(a)

"""
    sparsecsr(args...)
    sparsecsr(::Val{Bi},args...) where Bi

Create  a `SparseMatrixCSR` with `Bi`-based indexing (1 by default)
from the same `args...` as one constructs a `SparseMatrixCSC`
with the [`sparse`](@ref) function.
"""
sparsecsr(I,J,V) = SparseMatrixCSR(transpose(sparse(J,I,V,dimlub(J),dimlub(I))))
sparsecsr(I,J,V,m,n) = SparseMatrixCSR(transpose(sparse(J,I,V,n,m)))
sparsecsr(I,J,V,m,n,combine) = SparseMatrixCSR(transpose(sparse(J,I,V,n,m,combine)))
sparsecsr(::Val{Bi},I,J,V) where Bi = SparseMatrixCSR{Bi}(transpose(sparse(J,I,V,dimlub(J),dimlub(I))))
sparsecsr(::Val{Bi},I,J,V,m,n) where Bi = SparseMatrixCSR{Bi}(transpose(sparse(J,I,V,n,m)))
sparsecsr(::Val{Bi},I,J,V,m,n,combine) where Bi = SparseMatrixCSR{Bi}(transpose(sparse(J,I,V,n,m,combine)))
dimlub(I) = isempty(I) ? 0 : Int(maximum(I))

Base.convert(::Type{T},a::T) where T<:SparseMatrixCSR = a
function Base.convert(
  ::Type{SparseMatrixCSR{Bi,Tv,Ti}},a::SparseMatrixCSR{Bi}) where {Bi,Tv,Ti}
  rowptr = convert(Vector{Ti},a.rowptr)
  colval = convert(Vector{Ti},a.colval)
  nzval = convert(Vector{Tv},a.nzval)
  SparseMatrixCSR{Bi}(a.m,a.n,rowptr,colval,nzval)
end
function Base.convert(
  ::Type{SparseMatrixCSR{Bi,Tv,Ti}},a::SparseMatrixCSR{Bj}) where {Bi,Tv,Ti,Bj}
  rowptr = Vector{Ti}(undef,length(a.rowptr))
  colval = Vector{Ti}(undef,length(a.colval))
  rowptr .= a.rowptr .+ (Bi-Bj)
  colval .= a.colval .+ (Bi-Bj)
  nzval = convert(Vector{Tv},a.nzval)
  SparseMatrixCSR{Bi}(a.m,a.n,rowptr,colval,nzval)
end
function Base.convert(
  ::Type{SparseMatrixCSR{Bi,Tv,Ti}},a::Transpose{Tu,<:SparseMatrixCSC} where Tu) where {Bi,Tv,Ti}
  convert(SparseMatrixCSR{Bi,Tv,Ti},SparseMatrixCSR{Bi}(a))
end
function Base.convert(
  ::Type{SparseMatrixCSR{Bi,Tv,Ti}},a::AbstractMatrix) where {Bi,Tv,Ti}
  at = convert(SparseMatrixCSC{Tv,Ti},transpose(a))
  convert(SparseMatrixCSR{Bi,Tv,Ti},transpose(at))
end

function Base.copy(a::SparseMatrixCSR{Bi}) where Bi
  SparseMatrixCSR{Bi}(a.m,a.n,copy(a.rowptr),copy(a.colval),copy(a.nzval))
end

_copy_and_increment(x) = copy(x) .+ 1

function LinearAlgebra.fillstored!(a::SparseMatrixCSR,v)
  fill!(a.nzval,v)
  a
end

function LinearAlgebra.rmul!(a::SparseMatrixCSR,v::Number)
  rmul!(a.nzval,v)
  a
end

function LinearAlgebra.lu(a::SparseMatrixCSR{0})
  rowptr = _copy_and_increment(a.rowptr)
  colval = _copy_and_increment(a.colval)
  Transpose(lu(SparseMatrixCSC(a.m,a.n,rowptr,colval,a.nzval)))
end

function LinearAlgebra.lu(a::SparseMatrixCSR{1})
  Transpose(lu(SparseMatrixCSC(a.m,a.n,a.rowptr,a.colval,a.nzval)))
end

function LinearAlgebra.lu!(
    translu::Transpose{T,<:SuiteSparse.UMFPACK.UmfpackLU{T}},
    a::SparseMatrixCSR{1}) where {T}
  Transpose(lu!(translu.parent,SparseMatrixCSC(a.m,a.n,a.rowptr,a.colval,a.nzval)))
end

function LinearAlgebra.lu!(
  translu::Transpose{T,<:SuiteSparse.UMFPACK.UmfpackLU{T}},
  a::SparseMatrixCSR{0}) where {T}
  rowptr = _copy_and_increment(a.rowptr)
  colval = _copy_and_increment(a.colval)
  Transpose(lu!(translu.parent,SparseMatrixCSC(a.m,a.n,rowptr,colval,a.nzval)))
end

size(S::SparseMatrixCSR) = (S.m, S.n)
IndexStyle(::Type{<:SparseMatrixCSR}) = IndexCartesian()
function getindex(A::SparseMatrixCSR{Bi,T}, i0::Integer, i1::Integer) where {Bi,T}
  if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
  o = getoffset(A)
  r1 = Int(getrowptr(A)[i0]+o)
  r2 = Int(getrowptr(A)[i0+1]-Bi)
  (r1 > r2) && return zero(T)
  i1o = i1-o
  k = searchsortedfirst(colvals(A), i1o, r1, r2, Base.Order.Forward)
  ((k > r2) || (colvals(A)[k] != i1o)) ? zero(T) : nonzeros(A)[k]
end
function setindex!(A::SparseMatrixCSR{Bi,Tv,Ti}, _v, _i0::Integer, _i1::Integer) where {Bi,Tv,Ti}
  errmsg = "Trying to set an entry outside sparsity pattern"
  v  = convert(Tv,_v)
  i0 = convert(Ti,_i0)
  i1 = convert(Ti,_i1)
  if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
  o = getoffset(A)
  r1 = Int(getrowptr(A)[i0]+o)
  r2 = Int(getrowptr(A)[i0+1]-Bi)
  (r1 > r2) && throw(ArgumentError(errmsg))
  i1o = i1-o
  k = searchsortedfirst(colvals(A), i1o, r1, r2, Base.Order.Forward)
  if ((k > r2) || (colvals(A)[k] != i1o))
    throw(ArgumentError(errmsg))
  end
  A.nzval[k]=v
end

getrowptr(S::SparseMatrixCSR) = S.rowptr
getnzval(S::SparseMatrixCSR) = S.nzval
getcolval(S::SparseMatrixCSR) = S.colval

"""
    getBi(S::SparseMatrixCSR{Bi}) where {Bi}

Return `Bi`.
"""
getBi(S::SparseMatrixCSR{Bi}) where {Bi} = Bi

"""
    getoffset(S::SparseMatrixCSR{Bi}) where {Bi}

Return `1-Bi`. Useful to convert from 1-based to `Bi`-based indexing
(by subtracting the offset).
"""
getoffset(S::SparseMatrixCSR{Bi}) where Bi = getoffset(Bi)
@inline getoffset(Bi::Integer) = 1-Bi

"""
    issparse(S::SparseMatrixCSR)

Returns `true`.
"""
issparse(S::SparseMatrixCSR) = true

"""
    nnz(S::SparseMatrixCSR)

Returns the number of stored (filled) elements in a sparse array.
"""
nnz(S::SparseMatrixCSR) = length(nonzeros(S))

"""
    nonzeros(S::SparseMatrixCSR)

Return a vector (1-based) of the structural nonzero values in sparse array S.
This includes zeros that are explicitly stored in the sparse array.
The returned vector points directly to the internal nonzero storage of S,
and any modifications to the returned vector will mutate S as well.
"""
nonzeros(S::SparseMatrixCSR) = S.nzval

"""
    colvals(S::SparseMatrixCSR{Bi}) where {Bi}

Return a vector of the col indices of `S`. The stored values are indexes to arrays
with `Bi`-based indexing, but the `colvals(S)` array itself is a standard 1-based
Julia `Vector`.
Any modifications to the returned vector will mutate S as well.
Providing access to how the col indices are stored internally
can be useful in conjunction with iterating over structural
nonzero values. See also [`nonzeros`](@ref) and [`nzrange`](@ref).
"""
colvals(S::SparseMatrixCSR) = S.colval

"""
    nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi}

Return the range of indices to the structural nonzero values of a
sparse matrix row. The returned range of indices is always 1-based even for `Bi != 1`.
"""
nzrange(S::SparseMatrixCSR{Bi}, row::Integer) where {Bi} = S.rowptr[row]+getoffset(S):S.rowptr[row+1]-Bi

"""
    findnz(S::SparseMatrixCSR{Bi,Tv,Ti})

Return a tuple `(I, J, V)` where `I` and `J` are the row and column 1-based indices
of the stored ("structurally non-zero") values in sparse matrix A,
and V is a vector of the values. The returned vectors are newly allocated
and are unrelated to the internal storage of matrix `S`.
"""
function findnz(S::SparseMatrixCSR{Bi,Tv,Ti}) where {Bi,Tv,Ti}
  numnz = nnz(S)
  I = Vector{Ti}(undef, numnz)
  J = Vector{Ti}(undef, numnz)
  V = Vector{Tv}(undef, numnz)
  count = 1
  o = getoffset(S)
  @inbounds for row in 1:size(S,1)
    @inbounds for k in nzrange(S,row)
      I[count] = row
      J[count] = S.colval[k]+o
      V[count] = S.nzval[k]
      count += 1
    end
  end
  return (I, J, V)
end

"""
    count(pred, S::SparseMatrixCSR)
    count(S::SparseMatrixCSR)

Count the number of elements in `nonzeros(S)` for which predicate `pred` returns `true`.
If  `pred` not given, it counts the number of `true` values.
"""
count(pred, S::SparseMatrixCSR) = count(pred, nonzeros(S))
count(S::SparseMatrixCSR) = count(i->true, nonzeros(S))

function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector, α::Number, β::Number)
  A.n == size(v, 1) || throw(DimensionMismatch())
  A.m == size(y, 1) || throw(DimensionMismatch())
  if β != 1
    β != 0 ? rmul!(y, β) : fill!(y, zero(eltype(y)))
  end
  o = getoffset(A)
  for row = 1:size(y, 1)
    @inbounds for nz in nzrange(A,row)
      col = A.colval[nz]+o
      y[row] += A.nzval[nz]*v[col]*α
    end
  end
  return y
end

function mul!(y::AbstractVector,A::SparseMatrixCSR,v::AbstractVector)
  A.n == size(v, 1) || throw(DimensionMismatch())
  A.m == size(y, 1) || throw(DimensionMismatch())
  fill!(y, zero(eltype(y)))
  o = getoffset(A)
  for row = 1:size(y, 1)
    @inbounds for nz in nzrange(A,row)
      col = A.colval[nz]+o
      y[row] += A.nzval[nz]*v[col]
    end
  end
  return y
end

*(A::SparseMatrixCSR, v::Vector) = (y = similar(v,size(A,1));mul!(y,A,v))

function show(io::IO, ::MIME"text/plain", S::SparseMatrixCSR)
  xnnz = nnz(S)
  print(io, S.m, "×", S.n, " ", typeof(S), " with ", xnnz, " stored ",
        xnnz == 1 ? "entry" : "entries")
  if xnnz != 0
    print(io, ":")
    show(IOContext(io, :typeinfo => eltype(S)), S)
  end
end
show(io::IO, S::SparseMatrixCSR) = show(convert(IOContext, io), S::SparseMatrixCSR)

function show(io::IOContext, S::SparseMatrixCSR{Bi}) where{Bi}
  nnz(S) == 0 && return show(io, MIME("text/plain"), S)
  ioc = IOContext(io, :compact => true)

  function _format_line(r, col, padr, padc,o)
    print(ioc, "\n  [", rpad(col+o, padr), ", ", lpad(S.colval[r]+o, padc), "]  =  ")
    if isassigned(S.nzval, Int(r))
      show(ioc, S.nzval[r])
    else
      print(ioc, Base.undef_ref_str)
    end
  end

  function _get_cols(from, to)
    idx = eltype(S.rowptr)[]
    c = searchsortedlast(S.rowptr, from)
    for i = from:to
      while i == S.rowptr[c+1]
        c +=1
      end
      push!(idx, c)
    end
    idx
  end

  rows = displaysize(io)[1] - 4 # -4 from [Prompt, header, newline after elements, new prompt]
  if !get(io, :limit, false) || rows >= nnz(S) # Will the whole matrix fit when printed?
    cols = _get_cols(Bi, nnz(S)-getoffset(S))
    padr, padc = ndigits.((maximum(S.colval[1:nnz(S)]), cols[end]))
    _format_line.(1:nnz(S), cols.-getoffset(S), padr, padc,getoffset(S))
  else
    if rows <= 2
      print(io, "\n  \u22ee")
      return
    end
    s1, e1 = 1, div(rows - 1, 2) # -1 accounts for \vdots
    s2, e2 = nnz(S) - (rows - 1 - e1) + 1, nnz(S)
    cols1, cols2 = _get_cols(s1-getoffset(S), e1-getoffset(S)), _get_cols(s2-getoffset(S), e2-getoffset(S))
    padr = ndigits(max(maximum(S.colval[s1:e1]), maximum(S.colval[s2:e2])))
    padc = ndigits(cols2[end])
    _format_line.(s1:e1, cols1, padr, padc)
    print(io, "\n  \u22ee")
    _format_line.(s2:e2, cols2, padr, padc)
  end
  return
end


text:
#! /usr/bin/julia

# Rosetta Code, Order disjoint list items

function order_disjoint{T<:AbstractArray}(m::T, n::T)
    rlen = length(n)
    rdis = zeros(Int, rlen)
    for (i, e) in enumerate(n)
        j = findfirst(m, e)
        while j in rdis && j != 0
            j = findnext(m, e, j+1)
        end
        rdis[i] = j
    end
    if 0 in rdis
        throw(DomainError())
    end
    sort!(rdis)
    p = copy(m)
    p[rdis] = n
    return p
end

testm = {["the", "cat", "sat", "on", "the", "mat"],
         ["the", "cat", "sat", "on", "the", "mat"],
         ["A", "B", "C", "A", "B", "C", "A", "B", "C"],
         ["A", "B", "C", "A", "B", "D", "A", "B", "E"],
         ["A", "B"],
         ["A", "B"],
         ["A", "B", "B", "A"],
         }

testn = {["mat", "cat"],
         ["cat", "mat"],
         ["C", "A", "C", "A"],
         ["E", "A", "D", "A"],
         ["B"],
         ["B", "A"],
         ["B", "A"],
         }

for i in 1:length(testm)
    m = join(testm[i], " ")
    n = join(testn[i], " ")
    p = join(order_disjoint(testm[i], testn[i]), " ")
    println("    (", m, ", ", n, ") => ", p)
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Crazyflie
using Test

# TODO
# No non-instrusive packages for mocking tests
# Might need to use the python unittest library


text:
mutable struct Binding
    name::EXPR
    val::Union{Binding,EXPR,SymbolServer.SymStore,Nothing}
    type::Union{Binding,EXPR,SymbolServer.SymStore,Nothing}
    refs::Vector{Any}
    prev::Union{Binding,SymbolServer.SymStore,Nothing}
    next::Union{Binding,SymbolServer.SymStore,Nothing}
end
Binding(x::EXPR) = Binding(CSTParser.get_name(x), x, nothing, [], nothing, nothing)

function Base.show(io::IO, b::Binding)
    printstyled(io, "Binding(", Expr(b.name),
        b.type === nothing ? "" : ":: ",
        b.refs isa Vector ? "($(length(b.refs)) refs))" : ")", color = :blue)
end

hasbinding(x::EXPR) = hasmeta(x) && hasbinding(x.meta)
bindingof(x) = nothing
bindingof(x::EXPR) = bindingof(x.meta)

hasref(x::EXPR) = hasmeta(x) && hasref(x.meta)
refof(x::EXPR) = hasmeta(x) ? x.meta.ref : nothing

function gotoobjectofref(x::EXPR)
    r = refof(x)
    if r isa SymbolServer.SymStore
        return r
    elseif r isa Binding

    end
end

# Note to self, check consistency of marking self-reference of bindings (i.e.
# for, `function f end` we resolve `f` to itself at this stage.)
function mark_bindings!(x::EXPR, state)
    if hasbinding(x)
        return
    end
    if !hasmeta(x)
        x.meta = Meta()
    end
    if typof(x) === BinaryOpCall
        if kindof(x[2]) === CSTParser.Tokens.EQ
            if CSTParser.is_func_call(x[1])
                name = CSTParser.get_name(x)
                mark_binding!(x)
                mark_sig_args!(x[1])
                if isidentifier(name)
                    setref!(name, bindingof(x))
                end
            elseif typof(x[1]) === CSTParser.Curly
                mark_typealias_bindings!(x)
            else
                mark_binding!(x[1], x)
            end
        elseif kindof(x[2]) === CSTParser.Tokens.ANON_FUNC
            mark_binding!(x[1], x)
        end
    elseif typof(x) === WhereOpCall
        for i = 3:length(x)
            typof(x[i]) === PUNCTUATION && continue
            mark_binding!(x[i])
        end
    elseif typof(x) === CSTParser.For
        markiterbinding!(x[2])
    elseif typof(x) === CSTParser.Generator
        for i = 3:length(x)
            typof(x[i]) === PUNCTUATION && continue
            markiterbinding!(x[i])
        end
    elseif typof(x) === CSTParser.Filter
        for i = 1:length(x) - 2
            typof(x[i]) === PUNCTUATION && continue
            markiterbinding!(x[i])
        end
    elseif typof(x) === CSTParser.Flatten && length(x) === 1 && length(x[1]) >= 3 && length(x[1][1]) >= 3
        for i = 3:length(x[1][1])
            typof(x[1][1][i]) === PUNCTUATION && continue
            markiterbinding!(x[1][1][i])
        end
        for i = 3:length(x[1])
            typof(x[1][i]) === PUNCTUATION && continue
            markiterbinding!(x[1][i])
        end
    elseif typof(x) === CSTParser.Do
        if typof(x[3]) === CSTParser.TupleH
            for i in 1:length(x[3])
                typof(x[3][i]) === PUNCTUATION && continue
                mark_binding!(x[3][i])
            end
        end
    elseif typof(x) === FunctionDef || typof(x) === Macro
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.Function, [], nothing, nothing)
        if isidentifier(name)
            setref!(name, bindingof(x))
        end
        mark_sig_args!(CSTParser.get_sig(x))
    elseif CSTParser.defines_module(x)
        x.meta.binding = Binding(x[2], x, CoreTypes.Module, [], nothing, nothing)
        setref!(x[2], bindingof(x))
    elseif typof(x) === CSTParser.Try && length(x) > 3 && isidentifier(x[4])
        mark_binding!(x[4])
        setref!(x[4], bindingof(x[4]))
    elseif CSTParser.defines_datatype(x)
        name = CSTParser.get_name(x)
        x.meta.binding = Binding(name, x, CoreTypes.DataType, [], nothing, nothing)
        if isidentifier(name)
            setref!(name, bindingof(x))
        end
        mark_parameters(CSTParser.get_sig(x))
        if CSTParser.defines_struct(x) # mark field block
            blocki = typof(x[3]) === CSTParser.Block ? 3 : 4
            for i in 1:length(x[blocki])
                CSTParser.defines_function(x[blocki][i]) && continue
                mark_binding!(x[blocki][i])
            end
        end
    elseif typof(x) === CSTParser.Local
        if length(x) == 2
            if isidentifier(x[2])
                mark_binding!(x[2])
                setref!(x[2], bindingof(x[2]))
            elseif typof(x[2]) === CSTParser.TupleH
                for i = 1:length(x[2])
                    if isidentifier(x[2][i])
                        mark_binding!(x[2][i])
                        setref!(x[2][i], bindingof(x[2][i]))
                    end
                end
            end
        end

    end
end

function mark_binding!(x::EXPR, val = x)
    if typof(x) === CSTParser.Kw || (typof(x) === CSTParser.BinaryOpCall && kindof(x[2]) === CSTParser.Tokens.DECLARATION && typof(x[1]) === CSTParser.TupleH)
        mark_binding!(x[1], x)
    elseif typof(x) === CSTParser.TupleH || typof(x) === Parameters
        for arg in x
            typof(arg) === PUNCTUATION && continue
            mark_binding!(arg, val)
        end
    elseif typof(x) === CSTParser.InvisBrackets
        mark_binding!(CSTParser.rem_invis(x), val)
    elseif !(typof(x) == UnaryOpCall && kindof(x[1]) === CSTParser.Tokens.DECLARATION)
        if !hasmeta(x)
            x.meta = Meta()
        end
        x.meta.binding = Binding(CSTParser.get_name(x), val, nothing, [], nothing, nothing)
    end
    return x
end

function mark_parameters(sig::EXPR)
    signame = CSTParser.rem_where_subtype(sig)
    if typof(signame) === CSTParser.Curly
        for i = 3:length(signame) - 1
            typof(signame[i]) === PUNCTUATION && continue
            mark_binding!(signame[i])
        end
    end
    return sig
end

function markiterbinding!(iter::EXPR)
    if typof(iter) === BinaryOpCall && kindof(iter[2]) in (CSTParser.Tokens.EQ, CSTParser.Tokens.IN, CSTParser.Tokens.ELEMENT_OF)
        mark_binding!(iter[1], iter)
    elseif typof(iter) === CSTParser.Block
        for i = 1:length(iter)
            typof(iter[i]) === PUNCTUATION && continue
            markiterbinding!(iter[i])
        end
    end
    return iter
end

function mark_sig_args!(x::EXPR)
    if typof(x) === Call || typof(x) === CSTParser.TupleH
        if typof(x[1]) === CSTParser.InvisBrackets && typof(x[1][2]) === BinaryOpCall && kindof(x[1][2][2]) === CSTParser.Tokens.DECLARATION
            mark_binding!(x[1][2])
        end
        for i = 2:length(x) - 1
            a = x[i]
            if typof(a) === Parameters
                for j = 1:length(a)
                    aa = a[j]
                    if !(typof(aa) === PUNCTUATION)
                        mark_binding!(aa)
                    end
                end
            elseif !(typof(a) === PUNCTUATION)
                mark_binding!(a)
            end
        end
    elseif typof(x) === WhereOpCall
        for i in 3:length(x)
            if !(typof(x[i]) === PUNCTUATION)
                mark_binding!(x[i])
            end
        end
        mark_sig_args!(x[1])
    elseif typof(x) === BinaryOpCall
        if kindof(x[2]) == CSTParser.Tokens.DECLARATION
            mark_sig_args!(x[1])
        else
            mark_binding!(x[1])
            mark_binding!(x[3])
        end
    elseif typof(x) == UnaryOpCall && typof(x[2]) == CSTParser.InvisBrackets
        mark_binding!(x[2][2])
    end
end

function mark_typealias_bindings!(x::EXPR)
    mark_binding!(x, x)
    setscope!(x, Scope(x))
    for i = 2:length(x[1])
        if isidentifier(x[1][i])
            mark_binding!(x[1][i])
        elseif typof(x[1][i]) === BinaryOpCall && kindof(x[1][i][2]) === CSTParser.Tokens.ISSUBTYPE && isidentifier(x[1][i][1])
            mark_binding!(x[1][i][1])
        end
    end
    return x
end

function _in_func_def(x::EXPR)
    # only called in WhereOpCall
    # check 1st arg contains a call (or op call)
    ex = x[1]
    while true
        if typof(ex) === WhereOpCall || (typof(ex) === BinaryOpCall && kindof(ex[2]) === CSTParser.Tokens.DECLARATION)
            ex = ex[1]
        elseif typof(ex) === Call || (typof(ex) === BinaryOpCall && kindof(ex[2]) !== CSTParser.Tokens.DOT) || typof(ex) == CSTParser.UnaryOpCall
            break
        else
            return false
        end
    end
    # check parent is func def
    ex = x
    while true
        if !(parentof(ex) isa EXPR)
            return false
        elseif typof(parentof(ex)) === WhereOpCall || typof(parentof(ex)) === CSTParser.InvisBrackets
            ex = parentof(ex)
        elseif typof(parentof(ex)) === FunctionDef || CSTParser.is_assignment(parentof(ex))
            return true
        else
            return false
        end
    end
    return false
end

function add_binding(x, state, scope = state.scope)
    if bindingof(x) isa Binding
        bindingof(x).prev = nothing
        bindingof(x).next = nothing
        if typof(bindingof(x).name) === IDENTIFIER
            name = valof(bindingof(x).name)
        elseif typof(bindingof(x).name) === CSTParser.NONSTDIDENTIFIER
            name = valof(bindingof(x).name[2])
        elseif CSTParser.typof(bindingof(x).name) === CSTParser.MacroName
            name = string(Expr(bindingof(x).name))
        elseif CSTParser.typof(bindingof(x).name) === CSTParser.OPERATOR
            name = string(Expr(bindingof(x).name))
        else
            return
        end
        # check for global marker
        if isglobal(name, scope)
            scope = _get_global_scope(state.scope)
        end

        if typof(x) === Macro
            scope.names[string("@", name)] = bindingof(x)
            mn = CSTParser.get_name(x)
            if isidentifier(mn)
                setref!(mn, bindingof(x))
            end
        else
            if scopehasbinding(scope, name)
                bindingof(x).prev = scope.names[name]
                scope.names[name] = bindingof(x)
                bindingof(x).prev.next = bindingof(x)
            else
                # Checks for bindings which overwrite other module's bindings
                if typof(parentof(bindingof(x).name)) === CSTParser.Quotenode && typof(parentof(parentof(bindingof(x).name))) === BinaryOpCall && typof(parentof(parentof(bindingof(x).name))[1]) === IDENTIFIER # Only checks 1 level (e.g. M.name)
                    s1 = scope
                    while true
                        if s1.modules !== nothing
                            if scopehasmodule(s1, Symbol(valof(parentof(parentof(bindingof(x).name))[1]))) # this scope (s1) has a module with matching name
                                mod = getscopemodule(s1, Symbol(valof(parentof(parentof(bindingof(x).name))[1])))
                                if mod isa SymbolServer.ModuleStore && haskey(mod, Symbol(name))
                                    bindingof(x).prev = mod[Symbol(name)]
                                end
                            end
                            break # We've reached a scope that loads modules, no need to keep searching upwards
                        end
                        if parentof(s1) === nothing
                            break
                        else
                            s1 = parentof(s1)
                        end
                    end
                end
                # hoist binding for inner constructor to parent scope
                if (typof(scope.expr) === CSTParser.Struct || typof(scope.expr) === CSTParser.Mutable) && CSTParser.defines_function(x) && parentof(scope) isa Scope
                    return add_binding(x, state, parentof(scope))
                end
                scope.names[name] = bindingof(x)
            end
        end
        infer_type(bindingof(x), scope, state)
    elseif bindingof(x) isa SymbolServer.SymStore
        scope.names[valof(x)] = bindingof(x)
    end
end

isglobal(name, scope) = false
isglobal(name::String, scope) = scopehasbinding(scope, "#globals") && name in scope.names["#globals"].refs

function mark_globals(x::EXPR, state)
    if typof(x) === CSTParser.Global
        if !scopehasbinding(state.scope, "#globals")
            state.scope.names["#globals"] = Binding(EXPR(IDENTIFIER, EXPR[], 0, 0, "#globals", CSTParser.NoKind, false, nothing, nothing), nothing, nothing, [], nothing, nothing)
        end
        for i = 2:length(x)
            if typof(x[i]) === IDENTIFIER && !scopehasbinding(state.scope, valof(x[i]))
                push!(state.scope.names["#globals"].refs, valof(x[i]))
            end
        end

    end
end


text:
# Utility function without constraint
using Optim

kappa = 0.5
eta = 2.0
p1 = 1.0
p2 = 1.0
w = 10.0
e1 = (eta-1)/eta
e2 = eta/(eta-1)
initial_guess = [1.0]

function f_first!(storage, c1)
    storage[1] = -1*e2*
    ((kappa*c1[1]^(e1)+(1-kappa)*((w-p1*c1[1])/p2)^(e1))^(e2-1))*
    (e1*kappa*(c1[1])^((e1)-1)-e1*(1-kappa)*((w-p1*c1[1])/p2)^(e1-1)*(p1/p2))
end

function g_first!(storage, c1)
    storage[1] = -1*(e2)*
    ((kappa*(log(exp(c1[1])))^(e1)+(1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(e1))^((e2)-1))*
    ((e1)*kappa*(log(exp(c1[1])))^((e1)-1)-(e1)*(1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((e1)-1)*(p1/p2))
end

function g_second!(storage, c1)
    storage[1] = -1*(e2)*((e2)-1)*((kappa*log(exp(c1[1]))^(e1))+(1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(e1))^((e2)-2)*
    ((e1)*kappa*log(exp(c1[1]))^((e1)-1)-(e1)*(1-kappa)*(p1/p2)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((e1)-1))*
    ((e1)*kappa*log((exp(c1[1])))^((e1)-1)-(e1)*(1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((e1)-1)*(p1/p2))-
    ((e1)*((e1)-1)*kappa*log(exp(c1[1]))^((e1)-2)+(e1)*(1-kappa)*((e1)-1)*log(exp((w-p1*log(exp(c1[1]))/p2)))^((e1)-2)*(p1/p2)^2)*
    ((kappa*log(exp(c1[1]))^(e1)+(1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^(e1))^((e2)-1))

end

function utility_maximizer(kappa, eta, p1, p2, w, initial_guess, solver)
    if solver == "Newton"
        g(c1) =  -(kappa*log(exp(c1[1])))^e1 - ((1-kappa)*log(exp((w-p1*log(exp(c1[1]))/p2)))^e1)^e2
        output = Optim.minimizer(optimize(g, g_first!, g_second!, initial_guess, Newton()))
    elseif solver == "Conjugate gradient"
        lower = [1e-3]
        upper = [w/p1]
        f(c1) = -1*(kappa*c1[1]^e1+(1-kappa)*((w-p1*c1[1])/p2)^e1)^e2
        output = Optim.minimizer(optimize(f, f_first!, lower, upper, initial_guess, Fminbox(ConjugateGradient())))
    end
    return output
end

utility_maximizer(kappa, eta, p1, p2, w, initial_guess, "Newton")


text:
module SIREN

include("utils.jl")
include("dense.jl")
include("conv.jl")

export(uniform)
export(SIREN_init)
end

max_issues_repo_name: MacKenzieHnC/SIREN

max_forks_repo_name: MacKenzieHnC/SIREN


text:
function readinput(filename)
    map(value -> parse(Int, value), readlines(filename))
end

lines = readinput("d01.in")
length = size(lines, 1)

sum = 2020
P1, P2 = 0, 0

for i in 1:length
    for j in i+1:length
        if lines[i] + lines[j] == sum
            global P1 = lines[i] * lines[j]
        end
        for k in j+1:length
            if lines[i] + lines[j] + lines[k] == sum
                global P2 = lines[i] * lines[j] * lines[k]
            end
        end
    end
end

println("P1: ", P1)
println("P2: ", P2)

max_stars_repo_name: eduellery/adventofcode

max_issues_repo_name: eduellery/adventofcode

max_forks_repo_name: eduellery/adventofcode


text:
@testset "875.koko-eating-bananas.jl" begin
    @test min_eating_speed([3,6,7,11], 8) == 4
    @test min_eating_speed([30,11,23,4,20], 5) == 30
    @test min_eating_speed([30,11,23,4,20], 6) == 23
end


text:
module Schema

using UUIDs
using Dates
using Configurations

export GateInfo, DeviceInfo, NDUV, GateProp,
    Properties, Status, JobStatus, ExpData, ExpResult, Result,
    ExpConfig, Instruction, BooleanFunction, CopyFunction,
    Gate, Barrier, Reset, Measure, Experiment, Qobj,
    UserInfo

abstract type IBMQSchema end

@option struct Subscription <: IBMQSchema
    surveys::Int = 0
    newsletter::Int = 0
    updates::Int = 0
    tips::Int = 0
    updatesInTool::Int = 0
end

@option struct UserInfo <: IBMQSchema
    id::String
    username::String
    email::String
    urls::Dict{String, Any}
    userType::String="Standard"
    firstName::Maybe{String} = nothing
    lastName::Maybe{String} = nothing
    institution::Maybe{String} = nothing
    applications::Vector{String} = String[]
    subscriptions::Subscription = Subscription()
    readOnly::Bool=false
    needsRefill::Bool=false
    iqxPreferences::Dict{String, Any} = Dict{String, Any}()
    loginAccounts::Vector{Any} = []
    ibmQNetwork::Bool=true
    qNetworkRoles::Vector{String} = []
    roles::Vector{Any} = []
    emailVerified::Bool=false
    terms::Dict{String, Any} = Dict{String, Any}()
    canScheduleBackends::Bool=true
    servicesRoles::Vector{String}=[]
end

UserInfo(d::AbstractDict{String}) = from_dict(UserInfo, d)

"""
    struct GateInfo <: IBMQSchema

The gate config data structure has the following keys.

- `name`: Gate name, as it will be referred to in the OpenQASM circuit.
- `parameters`: List of parameters for the gate (empty if no parameters).
- `coupling_map` (optional): List of qubits that the gate applies to, each element of the list is
an n-qubit list where n is the size of the gate (e.g. 1-qubit gate, 2-qubit gate).
- `qasm_def`: OpenQASM definition of the gate in terms of the basis gates [U,CX].
Each unitary gate has an efficient representation in this basis.
- `conditional` (optional): Gate supports conditional operation (true/false). If not
listed then defaults to the backend setting.
- `latency_map` (optional): List for each gate of length n registers that indicates
if the feedback speed to the register is fast (1) or slow (0). In the above example
the u3 gate for qubit 0 has low latency for conditionals to register 0, but qubits
1 and 2 have low latency to both registers 1 and 2.
- `description` (optional): Description of the gate.
"""
@option struct GateInfo <: IBMQSchema
    name::String
    parameters::Maybe{Vector{String}} = nothing
    coupling_map::Maybe{Vector{Vector{Int}}} = nothing
    qasm_def::Maybe{String} = nothing
    conditional::Maybe{Bool} = nothing
    latency_map::Maybe{Vector{Vector{Int}}} = nothing
    description::Maybe{String} = nothing
end

"""
    struct DeviceInfo <: IBMQSchema

The backend will have a method Backend.configuration() which returns the required
backend config data structure. Backends can include additional items to this
structure.

- `backend_name`: Unique (to provider) backend identifier name. This could describe
a setup that goes through several changes, but retains common elements (e.g.,
for a physical device backend this could include the same coupling map and the
physical location, etc.)
backend version: Backend version string in the form `"X.X.X"`. Versions could
indicate, e.g., code changes, equipement upgrades, different cooldowns, new optimizations, etc.
- `n_qubits`: Number of qubits in the backend. Simulator backends return `"-1"`.
- `basis_gates`: List of the available gates on the backend as an array of gate names
(these should match the entries in gates).
- `coupling_map`: Representation of the physical coupling map on the device (the
coupling maps for each gate are defined in gates).
- `gates`: List of the available gates on the backend as a gate config data structure, see [`GateInfo`](@ref).
- `local`: Backend runs locally (`true`) or online (`false`).
- `simulator`: Backend is a simulator (`true`) or an experimental device (`false`).
- `conditional`: Backend supports conditional gates (`true`) or does not (`false`).
Individual gates may also support or not support conditionals.
- `configurable`: Backend (if simulator) is configurable (`true`). If true then there
are user specified configuration parameters (e.g., the topology, noise parameters,
etc.). The data structure for these settings is set by the specific backend.
- `n_registers` (required if conditional is `true`): Specifies the number of registers
slots (i.e. the number of register bits) that are available for conditional operations.
Each register can hold a bit value.
- `register_map` (required if conditional is `true`): Specifies the registers that each
qubit can store measurements. For this example, qubits 0 and 1 can store in
registers 0 and 1, qubit 2 in registers 2 and 3, and qubits 3 and 4 in registers 3
and 4.
- `open_pulse`: OpenPulse experiments are accepted on this backend.

The configuration structure may also have the following optional fields,

- `online_date`: Date that the backend was put online.
- `display_name`: Alternate name for the backend that is more descriptive that can
be used for display purposes.
- `sample_name`: Name of the sample for this given backend (likely blank for a
simulator).
- `description`: String to describe the backend.
- `url`: Internet address to the backend (if applicable).
- `tags`: List of tag strings for the backend that indicate true/false properties, e.g.,
`"credits_required"` (backend requires credits to run). Any absent tag means that
the property is false and new tags can be added.
- `quantum_volume`: quantum volume of given device.
- `max_experiments`: maximum experiment.
- `max_shots`: maximum shot can spawn.
"""
@option struct DeviceInfo <: IBMQSchema
    backend_name::String
    backend_version::VersionNumber
    n_qubits::Int
    memory::Bool = false
    max_shots::Maybe{Int} = nothing
    max_experiments::Maybe{Int} = nothing
    var"local"::Bool = false
    simulator::Bool = false
    conditional::Bool = true
    configurable::Bool = false
    credits_required::Bool = false
    allow_q_object::Bool = false
    allow_object_storage::Bool=false
    n_registers::Maybe{Int} = nothing
    open_pulse::Bool = false
    quantum_volume::Maybe{Int} = nothing
    online_date::Maybe{DateTime} = nothing
    display_name::Maybe{String} = nothing
    sample_name::Maybe{String} = nothing
    description::Maybe{String} = nothing
    url::Maybe{String} = nothing
    tags::Maybe{Vector{String}} = nothing
    register_map::Maybe{Vector{Vector{Int}}} = nothing
    coupling_map::Maybe{Vector{Vector{Int}}} = nothing
    basis_gates::Maybe{Vector{String}} = nothing
    gates::Vector{GateInfo}
end

function Configurations.convert_to_option(::Type{DeviceInfo}, ::Type{DateTime}, s::String)
    DateTime(s, dateformat"yyyy-mm-ddTHH:MM:SSZ")
end

"""
    struct NDUV <: IBMQSchema

Each of the parameters is expressed as a nduv struct (name-date-unit-value structure)
as below,
"""
@option struct NDUV <: IBMQSchema
    name::String
    date::DateTime
    unit::String
    value::Any
end

"""
    struct GateProp <: IBMQSchema

- `qubits`: Qubits involved in the gate.
- `gate`: Gate name, must be one of the gates from `"gates"` in the backend configuration structure.
- `parameters`: List of parameter structures which could generically include `"gate_err"`
(by gate error we mean the 1−Favg for the particular gate) and `"gate_time"`. Note
that each backend may measure gate error using different methodologies, this will
have to be conveyed by the backend over separate channels (e.g. at the URL).
"""
@option struct GateProp <: IBMQSchema
    qubits::Vector{Int}
    gate::String
    parameters::Vector{NDUV}
end

"""
    struct Properties <: IBMQSchema

The backend will have a call `Backend.properties()` which will return a backend props
data structure with backend properties (e.g. calibrations and coherences). Note that
this information is optionally provided by the backend, which will set how often and/or
under what conditions calibrations and characterizations need to be updated.

- `backend_name`, `backend_version`: Backend identifiers (from Backend.configuration())
that specify what backend these results were obtained from.
- `last_update_date`: Date/time of the last run calibration.
- `gates`: List of the qubit gate parameters (as a gate prop structure, see below).
- `qubits`: List of list of qubit parameters (e.g. coherences) which is in order of the
qubits. The qubit parameters could generically include “T1”, “T2”, “readoutErr”
and “frequency”.
- `general`: List of general backend parameters, see [`NDUV`](@ref).
"""
@option struct Properties <: IBMQSchema
    backend_name::String
    backend_version::VersionNumber
    last_update_date::DateTime
    gates::String
    qubits::Vector{Vector{NDUV}}
    general::Vector{NDUV}
end

"""
    struct Status <: IBMQSchema

The backend will have a call Backend.status() which returns status information on
the backend in the backend status structure (only the operational and status msg
fields are required),

- `backend_name`, `backend_version`: Backend identifiers that specify the backend.
- `operational`: Backend is operational (true/false), i.e., currently running jobs.
- `pending_jobs`: Number of jobs in the queue for the backend (if no queue return
0).
- `status_msg`: Status message for the backend. For example, `"The backend is down
for calibration, will be back at 19:00"`.
"""
@option struct Status <: IBMQSchema
    backend_name::String
    backend_version::VersionNumber
    operational::Bool
    pending_jobs::Int
    status_msg::String
end

"""
    struct JobStatus <: IBMQSchema

A call to `Job.status()` returns back a job status data structure of the following
form.

- `job_id`: Backend generated id corresponding to this job (this will only be nonzero
if the job has been successfully initialized and accepted to run on the backend).
- `status`: String value corresponding to the job status (`"ERROR"`,`"QUEUED"`,
`"INITIALIZING"`, `"RUNNING"`, `"CANCELLED"` and `"DONE"`).
- `status_msg`: Backend defined status message.
"""
@option struct JobStatus <: IBMQSchema
    job_id::String
    status::String
    status_msg::String
end

"""
    struct ExpData <: IBMQSchema

The measurement data that is returned in exp data has one of several possible forms:
as a histogram of counts of the memory states, the memory, or (for simulators) the
statevector or unitary matrix.

- `counts`: Histogram of counts in the different memory states. Only states with
non-zero counts are listed as keys. The states are labeled in hex (e.g., a 4 slot
memory `"1010"` (bit string) is decimal 10 and hex `"0xA"`).
- `memory`: State of the classical memory. For OpenQASM (or OpenPulse Level 2)
this is a list of hex strings indicating the state for each shot.
- `statevector` (optional): Final statevector corresponding to evolution of the zero state.
- `unitary` (optional): Final unitary matrix corresponding to the quantum circuit.
- `snapshots` (optional): Snapshots data structure that returns data as dictated by the snapshot
command used by simulators.

## Snapshots

The snapshot is a special command reserved for simulators which allows a "snapshot"
of the simulator state to be recorded.

```json
{
    "name": "snapshot",
    "label": "snap1",
    "type": "state"
}
```

- `name`: `"snapshot"`.
- `label`: Snapshot label which is used to identify the snapshot in the output.
- `type`: Type of snapshot, e.g., `"state"` (take a snapshot of the quantum state).
The types of snapshots offered are defined in a separate specification document
for simulators.
"""
@option struct ExpData <: IBMQSchema
    counts::Dict{String, Int}
    memory::Maybe{Vector{String}} = nothing
    statevector::Maybe{Matrix{ComplexF64}} = nothing
    unitary::Maybe{Matrix{ComplexF64}} = nothing
    snapshots::Maybe{Dict{String, Any}} = nothing
end

"""
    struct ExpResult <: IBMQSchema

Each individual experiment returns an exp result data structure.

- `shots`: If a single integer, then this is the number of shots taken to obtain this
data (``s = n₁``). If the backend allows asynchronous calls to measurement, the
value of ``n₂`` will increase as more data is taken. For backends that return the data
in sections (e.g. for bandwidth reasons) shots is given as a two-element list where
the data is from shot ``n₁`` to shot ``n₂`` (``s = n₂ − n₁``). The next call will give the
data starting at ``n₂ + 1``.
- `status`: Status message for this particular experiment.
- `success`: Success of the experiment (bool).
- `header` (optional): Header structure for the experiment that was passed in with the Qobj.
- `seed` (optional): Experiment seed (for simulator backends).
- `meas_return` (optional): String which determines whether the returned data is
averaged over the shots avg or contains each shot single. This is an OpenPulse
option, but could also apply to snapshots.
- `data`: Generic return experiment data structure exp data that will depend on
the type of experiment (`"QASM"` or `"PULSE"`) and/or the type of backend (e.g.
simulator data). See below.
"""
@option struct ExpResult <: IBMQSchema
    metadata::Maybe{Dict{String, Any}} = nothing
    header::Maybe{Dict{String, Any}} = nothing
    shots::Union{Int, Vector{Int}}
    status::String
    success::Bool
    time_taken::Float64
    seed_simulator::Maybe{Int} = nothing
    seed::Maybe{Int} = nothing
    meas_return::Maybe{String} = nothing
    data::ExpData
end

"""
    struct Result <: IBMQSchema

The results data structure from Job.result().

- `backend_name`, `backend_version`: Backend identifiers that specify what backend these results were obtained from.
- `qobj_id`: User generated id corresponding to the qobj id in the `Qobj`.
- `job_id`: Unique backend job identifier corresponding to these results.
- `date`: Date when the job was run.
- `header` (optional): Header structure for the job that was passed in with the `Qobj`.
- `results`: List of `m` (number of experiments) exp result data structures (defined below).
"""
@option struct Result <: IBMQSchema
    header::Maybe{Dict{String, Any}} = nothing
    metadata::Maybe{Dict{String, Any}} = nothing
    time_taken::Maybe{Float64} = nothing
    qobj_id::String
    job_id::String
    backend_name::String
    backend_version::VersionNumber
    date::DateTime
    status::String
    success::Bool
    results::Vector{ExpResult}
end

function Configurations.convert_to_option(::Type{Result}, ::Type{DateTime}, s::String)
    # NOTE: we only use millisecond since Dates doesn't support 6 digits precision
    # might worth an issue in upstream?
    DateTime(s[1:21], dateformat"yyyy-mm-ddTHH:MM:SS.s")
end

"""
    struct ExpConfig <: IBMQSchema

- `shots`: Number of times to repeat the experiment (for some simulators this may
be limited to 1, e.g., a unitary simulator).
- `memory_slots`: Number of classical memory slots used in this job. Memory slots
are used to record the results of qubit measurements and read out at the end of
an experiment. They cannot be used for feedback (those are the registers).
- `seed` (optional): Randomization seed for simulators.
- `max_credits` (optional): For credit-based backends, the maximum number of
credits that a user is willing to spend on this run (an error will be thrown if the
run required more than max credits).
"""
@option struct ExpConfig <: IBMQSchema
    shots::Int = 1024
    memory_slots::Int
    seed::Maybe{Int} = nothing
    max_credits::Maybe{Int} = nothing
end

"""
    abstract type Instruction <: IBMQSchema end

Abstract type for Qobj schema instructions.
"""
abstract type Instruction <: IBMQSchema end

Configurations.is_option(::Type{T}) where {T <: Instruction} = true

"""
    struct BooleanFunction <: Instruction

A boolean function is a command that takes in a register as an argument and computes
a boolean value that is written back into one of the register slots.

- `name`: `"bfunc"`.
- `mask`: Hex value which is applied as an AND to the register bits. In the given example ,`"0xF"` uses the first 4 bits of the register. The backend may put constraints
on the number of register bits that can be used in this function.
- `relation`: Relational operator for comparing the masked register to the val
(`"="`: equals, `"!="` not equals).
- `val`: Value to which to compare the masked register. In other words, the output
of the function is (register AND mask) relation val. In the above example this
is true when the first 4 register bits are `"0101"`
- `register`: Register slot in which to store the boolean function result. This
register value can then be used to apply conditional commands (see the following
sections).
- `memory` (optional): Memory slot in which to store the boolean function result.
"""
@option struct BooleanFunction <: Instruction
    name::String = "bfunc"
    mask::String
    relation::String
    val::String
    register::Int
    memory::Maybe{Int} = nothing
end

"""
    struct CopyFunction <: Instruction

A copy function is a command that copies a register slot.

- `name`: `"copy"`.
- `register_orig`: Register slot to copy.
- `register_copy`: Register slot(s) to copy to.
"""
@option struct CopyFunction <: Instruction
    name::String = "copy"
    register_orig::Int
    register_copy::Vector{Int}
end

"""
    struct Gate <: Instruction

- `name`: Name of the gate.
- `qubits`: List of qubits to apply the gate.
- `params` (optional): List of parameters for the gate (if the gate has parameters,
such as u1, u2, u3).
- `texparams` (optional): List of parameters for the gate in latex notation.
- `conditional` (optional): Apply the gate if the given register (in this example
register 3) is 1 (true) and conditionals are supported. If left blank then the gate
has no conditional element (i.e. no feedback). By default this is blank.
"""
@option struct Gate <: Instruction
    name::String
    qubits::Vector{Int}
    params::Maybe{Vector{Float64}} = nothing
    texparams::Maybe{Vector{String}} = nothing
    conditional::Maybe{Int} = nothing
end

"""
    struct Barrier <: Instruction

- `name`: `"barrier"`.
- `qubits`: List of qubits on which to apply the barrier.
"""
@option struct Barrier <: Instruction
    name::String = "barrier"
    qubits::Vector{Int}
end

"""
    struct Reset <: Instruction

- `name`: `"reset"`.
- `qubits`: List of qubits to reset.
"""
@option struct Reset <: Instruction
    name::String = "reset"
    qubits::Vector{Int}
end

"""
    struct Measure <: Instruction

- `name`: `"measure"`.
- `qubits`: List of qubits to measure.
- `memory`: List of memory slots in which to store the measurement results (must
be the same length as qubits). Subsequent measurements that write to the same
memory slot will overwrite the previous measurement.
- `register` (optional): List of register slots in which to store the measurement
results (must be the same length as qubits). These can be used for fast feedback
(if allowed). The allowed slots for a qubit may be constrained by the backend
register map.
"""
@option struct Measure <: Instruction
    name::String = "measure"
    qubits::Vector{Int}
    memory::Vector{Int}
    register::Maybe{Vector{Int}} = nothing
end

"""
    struct Experiment <: IBMQSchema

- `header` (optional): User-defined structure that contains metadata on each experiment and
is not used by the backend. The header will be passed through to the result data
structure unchanged. For example, this may contain a fitting parameters for the
experiment. In addition, this header can contain a mapping of backend memory
and backend qubits to OpenQASM registers. This is because an OpenQASM
circuit may contain multiple classical and quantum registers, but Qobj flattens
them into a single memory and single set of qubits.
- `config` (optional): Configuration structure for user settings that can be different in each
experiment. These will override the configuration settings of the whole job. See [`ExpConfig`](@ref).
`instructions`: List of sequence commands that define the experiment. See [`Instruction`](@ref).
"""
@option struct Experiment <: IBMQSchema
    header::Maybe{Dict{String, Any}} = nothing
    # NOTE: the schema specification didn't mention
    # but this can be optional
    config::Maybe{ExpConfig} = nothing
    instructions::Vector{Instruction} = Instruction[]
end

"""
    struct Qobj <: IBMQSchema

Experiments are loaded in through the backend using a Qobj data structure which
encapsulates the user configuration settings and experiment sequences.

- `qobj_id`: User generated run identifier
- `type`: Type of experiment, can be either `"QASM"` for openQASM experiments or `"PULSE"` for OpenPulse experiments.
- `schema_version`: Version of the schema that was used to generate and validate this `Qobj`.
- `experiments`: : List of m experiment sequences to run. Each experiment is an
experiment data structure. Each experiment is run once in the order that they
are specified in this list and then the sequence is repeated until the specified
number of shots has been performed.
- `header` (optional): User-defined structure that contains metadata on the job and is not used
by the backend. The header will be passed through to the result data structure
unchanged. For example, this may contain a description of the full job and/or
the backend that the experiments were compiled for.
- `config`: Configuration settings structure, see also [`ExpConfig`](@ref).
"""
@option struct Qobj <: IBMQSchema
    qobj_id::String = string(uuid1())
    type::String = "QASM"
    schema_version::VersionNumber = v"1.3"
    experiments::Vector{Experiment} = Experiment[]
    header::Maybe{Dict{String, Any}} = nothing
    config::ExpConfig
end

Configurations.to_dict(::Type{T}, x::VersionNumber) where {T <: IBMQSchema} = string(x)

end


text:
using MUMPS
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
include("getDivGrad.jl");
Ns = (8,16,24,32,48,64)
MUMPStime=zeros(length(Ns))
Juliatime=zeros(length(Ns))
nrhs = 1

for i=1:length(Ns)
	println(@sprintf("Solve div-grad on %d x %d x %d grid",Ns[i],Ns[i],Ns[i]));
	A = getDivGrad(Ns[i],Ns[i],Ns[i]);
	n = size(A,1);
    r = ones(n);
	A = A + im*spdiagm(r,0);
	rhs = randn(n,nrhs)+ im*randn(n,nrhs)

	# solve using mumps
	tic()
	x = solveMUMPS(A,rhs,1,0);
	MUMPStime[i]= toc();

	# solve using mumps
	tic()
	# x = solveMUMPS(A,rhs,[],1,1);
	x = A\rhs;
	Juliatime[i]= toc();
end

println(@sprintf("| %4s | %8s | %8s | %8s |","n","MUMPS","Julia","speedup"))

for i=1:length(Ns)
	println(@sprintf("| %d^3  | %1.5f | %1.5f | %1.5f |",Ns[i],MUMPStime[i],Juliatime[i],Juliatime[i]/MUMPStime[i]))
end


text:
buttons = Dict([
    # This button is disabled until the game starts.
    :bRestart =>
        MenuButton(false, UIPixelPos(0,0), UIPixelDims(200, 30), "New Game",
            ()->(resetGame(); buttons[:bNewContinue].callBack();))
    # Note that this text changes to "Continue" after first press.
    :bNewContinue =>
        MenuButton(true, UIPixelPos(0,0), UIPixelDims(200, 30), "New Game",
               ()->(global paused,game_started,buttons;
                    paused[] = false; game_started[] = true;
                    buttons[:bNewContinue].text = "Continue"; # After starting game
                    buttons[:bRestart].enabled = true;        # After starting game
                    ))
    :bSoundToggle =>
        CheckboxButton(true,
            MenuButton(true, UIPixelPos(0,0), UIPixelDims(200, 30), "Sound on/off",
                (enabled)->(toggleAudio(enabled)))
            )
    :bQuit =>
        MenuButton(true, UIPixelPos(0,0), UIPixelDims(200, 30), "Quit",
            ()->(global paused, playing; paused[] = playing[] = false;))

     # Key controls buttons
    :keyALeft =>
        KeyButton(true, UIPixelPos(0,0), UIPixelDims(120, 20), keyDisplayNames[keySettings[:keyP1Collector]],
               ()->(tryChangingKeySettingButton(:keyALeft)))
    :keyARight =>
        KeyButton(true, UIPixelPos(0,0), UIPixelDims(120, 20), keyDisplayNames[keySettings[:keyP1Fighter]],
               ()->(tryChangingKeySettingButton(:keyARight)))
    :keyP2Collector =>
        KeyButton(true, UIPixelPos(0,0), UIPixelDims(120, 20), keyDisplayNames[keySettings[:keyP2Collector]],
               ()->(tryChangingKeySettingButton(:keyP2Collector)))
    :keyP2Fighter =>
        KeyButton(true, UIPixelPos(0,0), UIPixelDims(120, 20), keyDisplayNames[keySettings[:keyP2Fighter]],
               ()->(tryChangingKeySettingButton(:keyP2Fighter)))

    :bResetDefaultKeys =>
        KeyButton(false, UIPixelPos(0,0), UIPixelDims(240, 30), "Reset Default Controls",
               ()->(resetDefaultKeys()))
  ])
paddleAControlsX() = screenCenterX()-260
paddleBControlsX() = screenCenterX()+260
function recenterButtons!()
    global buttons
    buttons[:bRestart].pos     = screenOffsetFromCenter(0,-25)
    buttons[:bNewContinue].pos = screenOffsetFromCenter(0,9)
    buttons[:bSoundToggle].button.pos = screenOffsetFromCenter(0,43)
    buttons[:bQuit].pos        = screenOffsetFromCenter(0,77)
    buttons[:keyALeft].pos    = UIPixelPos(paddleAControlsX(), winHeight[]-147)
    buttons[:keyARight].pos   = UIPixelPos(paddleAControlsX(), winHeight[]-122)
    buttons[:keyP2Collector].pos    = UIPixelPos(paddleBControlsX(), winHeight[]-147)
    buttons[:keyP2Fighter].pos   = UIPixelPos(paddleBControlsX(), winHeight[]-122)
    buttons[:bResetDefaultKeys].pos   = UIPixelPos(screenCenterX(), winHeight[]-102)
end
function toggleAudio(enabled)
    global audioEnabled;
    audioEnabled = enabled;
    if (audioEnabled) SDL2.Mix_ResumeMusic()
    else  SDL2.Mix_PauseMusic()
    end
end