text:
export AstrobeeSE2
export init_traj_straightline, init_traj_geometricplan

mutable struct AstrobeeSE2 <: DynamicsModel
  # state: r v p omega
  x_dim
  u_dim
  clearance

  # Parameters that can be updated
  f::Vector
  A
  B
end

function AstrobeeSE2()
  x_dim = 6
  u_dim = 3
  clearance = 0.05
  AstrobeeSE2(x_dim, u_dim, clearance, [], [], [])
end

function SCPParam(model::AstrobeeSE2, fixed_final_time::Bool)
  convergence_threshold = 0.01

  SCPParam(fixed_final_time, convergence_threshold)
end

function SCPParam_GuSTO(model::AstrobeeSE2)
  Δ0 = 2.
  ω0 = 1.
  ω_max = 1.0e10
  ε = 1.0e-6
  ρ0 = 0.01
  ρ1 = 0.3
  β_succ = 2.
  β_fail = 0.5
  γ_fail = 5.

  SCPParam_GuSTO(Δ0, ω0, ω_max, ε, ρ0, ρ1, β_succ, β_fail, γ_fail)
end

function SCPParam_Mao(model::AstrobeeSE2)
  ρ = [0.;0.25;0.9]
  Δ_u0 = 0.1
  λ = 1.
  α = 2.
  SCPParam_Mao(ρ, Δ_u0, λ, α)
end

function SCPParam_TrajOpt(model::AstrobeeSE2)
  μ0 = 1.
  s0 = 1.
  c = 10.
  τ_plus = 2.
  τ_minus = 0.5
  k = 5.
  ftol = 0.01
  xtol = 0.1
  ctol = 0.01
  max_penalty_iteration = 5
  max_convex_iteration = 5
  max_trust_iteration = 5

  SCPParam_TrajOpt(μ0, s0, c, τ_plus, τ_minus, k, ftol, xtol, ctol,max_penalty_iteration,max_convex_iteration,max_trust_iteration)
end

###############
# Gurobi stuff
###############
function cost_true(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2}) where T
  U, N = traj.U, SCPP.N
  dtp = traj_prev.dt
  Jm = 0
  for k in 1:N-1
    Jm += dtp*norm(U[:,k])^2
  end
  return Jm
end

function cost_true_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  cost_true(traj, traj_prev, SCPP)
end

#############################
# Trajectory Initializations
#############################
function init_traj_straightline(TOP::TrajectoryOptimizationProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  model, x_init, x_goal = TOP.PD.model, TOP.PD.x_init, TOP.PD.x_goal
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess
  N = TOP.N

  X = hcat(range(x_init, stop=x_goal, length=N)...)
  U = zeros(u_dim, N)
  Trajectory(X, U, tf_guess)
end

# TODO(acauligi): Add geometric plan
function int_traj_geometricplan(TOP::TrajectoryOptimizationProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  return Trajectory(TOP)  # Placeholder
end

####################
# Constraint-adding
####################
function initialize_model_params!(SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim = model.x_dim, model.u_dim
  Xp, Up = traj_prev.X, traj_prev.U

  model.f, model.A, model.B = [], A_dyn(Xp[:,1],robot,model), B_dyn(Xp[:,1],robot,model)
  for k = 1:N-1
    push!(model.f, f_dyn(Xp[:,k],Up[:,k],robot,model))
  end

  SCPP.PD.robot.Jcollision = []
  Jcollision = SCPP.PD.robot.Jcollision
  for k = 1:N
    push!(Jcollision, [Eye(2) zeros(2,1); zeros(1,3)])
  end
end

function update_model_params!(SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  Xp, Up, f = traj_prev.X, traj_prev.U, model.f

  for k = 1:N-1
    update_f!(f[k], Xp[:,k], Up[:,k], robot, model)
  end

  xb, Jcollision = SCPP.PD.robot.xb, SCPP.PD.robot.Jcollision
  for k = 1:N
    thetap = traj_prev.X[3,k]
    Jcollision[k][1,3] = -xb[1]*sin(thetap)-xb[2]*cos(thetap)
    Jcollision[k][2,3] =  xb[1]*cos(thetap)-xb[2]*sin(thetap)
  end
end

macro constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  quote
    X, U, Tf = $(esc(traj)).X, $(esc(traj)).U, $(esc(traj)).Tf
    Xp, Up, Tfp, dtp = $(esc(traj_prev)).X, $(esc(traj_prev)).U, $(esc(traj_prev)).Tf, $(esc(traj_prev)).dt
    robot, model, WS, x_init, x_goal = $(esc(SCPP)).PD.robot, $(esc(SCPP)).PD.model, $(esc(SCPP)).WS, $(esc(SCPP)).PD.x_init, $(esc(SCPP)).PD.x_goal
    x_dim, u_dim, N, dh = model.x_dim, model.u_dim, $(esc(SCPP)).N, $(esc(SCPP)).dh
    X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh
  end
end

## Dynamics constraints
function dynamics_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)

  return A_dyn_discrete(X[:,k],dtp,robot,model)*X[:,k] + B_dyn_discrete(X[:,k],dtp,robot,model)*U[:,k] - X[:,k+1]

  # fp, Ap, Bp = get_f(k, model), get_A(k, model), get_B(k, model)
  # if k == N-1
  #   return Tf*fp + Tfp*(Ap*(X[:,k]-Xp[:,k]) + Bp*(U[:,k]-Up[:,k])) - (X[:,k+1]-X[:,k])/dh
  # else
  #   return 0.5*(Tf*(fp + get_f(k+1, model)) + Tfp*(Ap*(X[:,k]-Xp[:,k]) + Bp*(U[:,k]-Up[:,k])) +
  #     Tfp*(Ap*(X[:,k+1]-Xp[:,k+1]) + Bp*(U[:,k+1]-Up[:,k+1]))) - (X[:,k+1]-X[:,k])/dh
  # end
end

# Get current dynamics structures for a time step
get_f(k::Int, model::AstrobeeSE2) = model.f[k]
get_A(k::Int, model::AstrobeeSE2) = model.A
get_B(k::Int, model::AstrobeeSE2) = model.B

# Generate full dynamics structures for a time step
function f_dyn(x::Vector, u::Vector, robot::Robot, model::AstrobeeSE2)
  x_dim = model.x_dim
  f = zeros(x_dim)
  update_f!(f, x, u, robot, model)
  return f
end

function update_f!(f, x::Vector, u::Vector, robot::Robot, model::AstrobeeSE2)
  F, M = u[1:2], u[3]
  f[4:5] = 1/robot.mass*F
  f[6] = robot.Jinv*M
end

function A_dyn(x::Vector, robot::Robot, model::AstrobeeSE2)
  kron([0 1; 0 0], Eye(3))
end

function B_dyn(x::Vector, robot::Robot, model::AstrobeeSE2)
  B = zeros(6,3)
  B[1:6,1:3] = kron([0;1], Eye(3))
  B[4:5,:] = 1/robot.mass * B[4:5,:]
  B[6,:] = robot.Jinv * B[6,:]
  return B
end

# Generate full discrete update version of dynamics matrices for a time step
# TODO(ambyld): Rename these? Initialize these once?
function A_dyn_discrete(x, dt, robot::Robot, model::AstrobeeSE2)
  kron([1 dt; 0 1], Eye(3))
end

function B_dyn_discrete(x, dt, robot::Robot, model::AstrobeeSE2)
  Jzz_inv = robot.Jinv
  B = [0.5*dt^2*Eye(3); dt*Eye(3)]
  B[1:2,:] *= 1/robot.mass  # translational force -> acceleration
  B[4:5,:] *= 1/robot.mass
  B[3,:] *= Jzz_inv      # rotational moment -> acceleration
  B[6,:] *= Jzz_inv
  return B
end

## Convex state inequality constraints
function csi_translational_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  return norm(X[4:5,k]) - robot.hard_limit_vel
end

function csi_angular_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  return norm(X[6,k]) - robot.hard_limit_omega
end

## Convex control inequality constraints
function cci_translational_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  return 1/robot.mass*norm(U[1:2,k]) - robot.hard_limit_accel
end

function cci_angular_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  Jzz_inv = robot.Jinv
  return abs(Jzz_inv*U[3,k]) - robot.hard_limit_α
end

## Nonconvex state inequality constraints
function ncsi_obstacle_avoidance_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)

  rb_idx, env_idx = 1, i
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  r = [X[1:2,k];0]
  dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

  # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
  return clearance - dist
end

## Nonconvex state inequality constraints (convexified)
function ncsi_body_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)

  env_ = WS.btenvironment_keepout
  env_idx = i
  clearance = model.clearance

  # Base
  rb_idx = 1
  r0 = get_workspace_location(traj_prev, SCPP, k)
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r0, env_idx)
  r = get_workspace_location(traj, SCPP, k)

  if dist < SCPP.param.obstacle_toggle_distance
    # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return clearance - (dist + nhat'*(r-r0))
  else
    return 0.
  end
end

function ncsi_arm_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)

  env_ = WS.btenvironment_keepout
  env_idx = i
  clearance = model.clearance

  r0 = get_workspace_location(traj_prev, SCPP, k)

  # Arm
  rb_idx = 2
  xb = robot.xb
  thetap = Xp[3,k]
  R01p = [cos(thetap) -sin(thetap) 0.;
          sin(thetap)  cos(thetap) 0.;
          0. 0. 1.]
  r1 = R01p*xb + r0
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r1, env_idx)

  # return 0.   # toggle off

  if dist < SCPP.param.obstacle_toggle_distance
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return clearance - (dist + nhat'*robot.Jcollision[k]*(X[1:3,k]-r0))
  else
    return 0.
  end
end

## State trust region inequality constraints
function stri_state_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)

  return norm(X[:,k]-Xp[:,k])^2
end

## Convex state inequality constraints
function ctri_control_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  return norm(U[:,k]-Up[:,k])
end

function get_workspace_location(traj, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, k::Int, i::Int=0) where {T,E}
  return [traj.X[1:2,k]; 0]
end

## Constructing full list of constraint functions
function SCPConstraints(SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  model = SCPP.PD.model
  x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N
  WS = SCPP.WS

  SCPC = SCPConstraints()

  ## Dynamics constraints
  for k = 1:N-1
    push!(SCPC.dynamics, (dynamics_constraints, k, 0))
  end

  ## Convex state equality constraints
  # Init and goal (add init first for convenience of getting dual)
  for i = 1:x_dim
    push!(SCPC.convex_state_eq, (cse_init_constraints, 0, i))
  end
  for i = 1:x_dim
    push!(SCPC.convex_state_eq, (cse_goal_constraints, 0, i))
  end

  ## Convex state inequality constraints
  for k = 1:N
    push!(SCPC.convex_state_ineq, (csi_translational_velocity_bound, k, 0))
    push!(SCPC.convex_state_ineq, (csi_angular_velocity_bound, k, 0))
  end

  ## Nonconvex state equality constraints (convexified)
  nothing

  ## Nonconvex state inequality constraints (convexified)
  env_ = WS.btenvironment_keepout
  for k = 1:N, i = 1:length(env_.convex_env_components)
    push!(SCPC.nonconvex_state_convexified_ineq, (ncsi_body_obstacle_avoidance_constraints_convexified, k, i))
    # push!(SCPC.nonconvex_state_convexified_ineq, (ncsi_arm_obstacle_avoidance_constraints_convexified, k, i))
  end

  ## Convex control equality constraints
  nothing

  ## Convex control inequality constraints
  for k = 1:N-1
    push!(SCPC.convex_control_ineq, (cci_translational_accel_bound, k, 0))
    push!(SCPC.convex_control_ineq, (cci_angular_accel_bound, k, 0))
  end

  ## State trust region ineqality constraints
  for k = 1:N
    push!(SCPC.state_trust_region_ineq, (stri_state_trust_region, k, 0))
  end

  ## Constrol trust region inequality constraints
  for k = 1:N-1
    push!(SCPC.control_trust_region_ineq, (ctri_control_trust_region, k, 0))
  end

  return SCPC
end

# TODO: Generalize this? Need to make A always a vector
function trust_region_ratio_gusto(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  fp, Ap = model.f, model.A
  num,den = 0, 0
  env_ = WS.btenvironment_keepout

  for k in 1:N-1
    linearized = fp[k] + Ap*(X[:,k]-Xp[:,k])
    num += norm(f_dyn(X[:,k],U[:,k],robot,model) - linearized)
    den += norm(linearized)
  end

  clearance = model.clearance

  # TODO(ambyld): Update for arm
  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx, body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        linearized = clearance - (dist + nhat'*(r-r0))

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

        num += abs((clearance-dist) - linearized)
        den += abs(linearized)
      end
    end
  end
  return num/den
end

function trust_region_ratio_trajopt(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  fp = model.f
  num, den = 0, 0
  env_ = WS.btenvironment_keepout
  dt = traj.dt

  for k in 1:N-1
    ϕ_old = norm(fp[k] - (Xp[:,k]-Xp[:,k])/dtp, 1)
    ϕ_new = norm(f_dyn(X[:,k],U[:,k],robot,model) - (X[:,k]-X[:,k])/dt, 1)
    ϕ_hat_new = norm(dynamics_constraints(traj,traj_prev,SCPP,k,0), 1)
    num += (ϕ_old-ϕ_new)
    den += (ϕ_old-ϕ_hat_new)
  end

  clearance = model.clearance

  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx,body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        ϕ_old = clearance-dist

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        ϕ_new = clearance-dist
        ϕ_hat_new = clearance - (dist + nhat'*(r-r0))

        num += (ϕ_old-ϕ_new)
        den += (ϕ_old-ϕ_hat_new)
      end
    end
  end

  return num/den
end

function trust_region_ratio_mao(traj, traj_prev::Trajectory, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,x_goal,x_dim,u_dim,N,dh = @constraint_abbrev_astrobeeSE2(traj, traj_prev, SCPP)
  num,den = 0, 0

  cost_true_prev = cost_true(traj_prev, traj_prev, SCPP)
  cost_true_new = cost_true(traj, traj, SCPP)
  for k in 1:N-1
    cost_true_prev += norm(dynamics_constraints(traj_prev, traj_prev, SCPP, k, 0), Inf)
    cost_true_new += norm(dynamics_constraints(traj, traj, SCPP, k, 0), Inf)
  end
  cost_linearized_new = cost_true(traj, traj, SCPP)
  return (cost_true_prev-cost_true_new)/(cost_true_prev-cost_linearized_new)
end

function get_dual_cvx(prob::Convex.Problem, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}, solver) where {T,E}
	if solver == "Mosek"
		return -MathProgBase.getdual(prob.model)[1:SCPP.PD.model.x_dim]
	else
		return []
	end
end

function get_dual_jump(SCPS::SCPSolution, SCPP::SCPProblem{Astrobee2D{T}, AstrobeeSE2, E}) where {T,E}
	@show -MathProgBase.getconstrduals(SCPS.solver_model.internalModel)[1:SCPP.PD.model.x_dim]
end

max_stars_repo_name: StanfordASL/GuSTO

max_issues_repo_name: StanfordASL/GuSTO

max_forks_repo_name: StanfordASL/GuSTO


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

using JSON
using Libdl
using LinearAlgebra

################################################################################
# First try to detect and load existing libmxnet
################################################################################
libmxnet_detected = false
libmxnet_curr_ver = get(ENV, "MXNET_COMMIT", "master")
curr_win = "20180211"  # v1.1.0

if haskey(ENV, "MXNET_HOME")
  MXNET_HOME = ENV["MXNET_HOME"]
  @info("MXNET_HOME environment detected: $MXNET_HOME")
  @info("Trying to load existing libmxnet...")
  # In case of macOS, if user build libmxnet from source and set the MXNET_HOME,
  # the output is still named as `libmxnet.so`.
  lib = Libdl.find_library(["libmxnet.$(Libdl.dlext)", "libmxnet.so"])
  if !isempty(lib)
    @info("Existing libmxnet detected at $lib, skip building...")
    libmxnet_detected = true
  else
    error("Failed to load existing libmxnet")
  end
end

ext: jl


text:
using LinearAlgebra: sqrt
using Printf: Threads
using Plots: Threads
using ProgressMeter
using Printf
include("PolygonBase.jl")
include("ESBase.jl")
include("silhouette&pieces.jl")

# X[1つめのピースのx座標, 1つめのピースのy座標, １つめの回転角Θ, 2つめのx座標, ...]
function loss_poly(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    tmp = LibGEOS.intersection(unioned, MYPolygon2LibGEOS(silhouette))
    A = LibGEOS.area(tmp) / LibGEOS.area(MYPolygon2LibGEOS(silhouette))

    # 各2pieceの共通部分を求める
    inter = 1.0 - LibGEOS.area(unioned) / sum_pieces
    outer = 1.0 - LibGEOS.area(tmp)/sum_pieces

    # 頂点の誤差成分 (シルエットベース 緩い)
    Δv = 0.0
    for v_s in 1:silhouette.n
        Δmin = 100.0
        for j in 1:length(pieces)
            tmp_p = move(pieces[j], X[3j - 2], X[3j - 1])
            rotate!(tmp_p, X[3j] * 2π)
            for v_p in 1:tmp_p.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - tmp_p.vertexes[:, v_p]).^2))
            end
        end
        Δv += Δmin / silhouette.n
        # Δv = max(Δv, Δmin)
    end

    # 頂点の誤差成分 (ピースベース 厳しい)
    """
    Δv = 0.0
    for i in 1:length(pieces)
        p = move(pieces[i], X[3i - 2], X[3i - 1])
        rotate!(p, X[3i] * 2π)
        Δmax = 0.0
        # 頂点ごとに誤差を見る
        for v_p in 1:p.n
            Δmin = 100.0
            # シルエットについて見る
            for v_s in 1:silhouette.n
                Δmin = min(Δmin, sum((silhouette.vertexes[:,v_s] - p.vertexes[:, v_p]).^2))
            end
            # その他のピースについて見る
            for j in 1:length(pieces)
                if j == i continue end
                q = move(pieces[j], X[3j - 2], X[3j - 1])
                rotate!(p, X[3j] * 2π)
                for v_q in 1:q.n
                    Δmin = min(Δmin, sum((q.vertexes[:,v_q] - p.vertexes[:, v_p]).^2))
                end
            end
            Δmax = max(Δmax, Δmin)
        end
        # Δv += 2tanh(2Δmin) / (π*silhouette.n)
        # Δv = max(Δv, Δmin)
        Δv += Δmax
    end
    """

    # 辺の誤差成分
    E = 0.0
    for i in 1:length(pieces)
        P = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(P, X[3i] * 2π)

        cosθ = 0.0

        # シルエットの各頂点と距離が最小な点を探す
        min_d = 100.0
        for p in 0:P.n-1, s in 0:silhouette.n-1
            d = sum((P.vertexes[:, p + 1] - silhouette.vertexes[:, s + 1]).^2)
            # 最小の大きさの角度を探す
            if d < min_d
                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ1 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p, P.n) + 1] - P.vertexes[:, mod(p-1, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ2 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s, silhouette.n) + 1] - silhouette.vertexes[:, mod(s-1, silhouette.n) + 1]
                cosθ3 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                eP = P.vertexes[:, mod(p+1, P.n) + 1] - P.vertexes[:, mod(p, P.n) + 1]
                eS = silhouette.vertexes[:, mod(s+1, silhouette.n) + 1] - silhouette.vertexes[:, mod(s, silhouette.n) + 1]
                cosθ4 = abs(eP' * eS / (norm(eP) * norm(eS) ))

                cosθ = max(cosθ1, cosθ2, cosθ3, cosθ4)
            end
        end
        E += cosθ / length(pieces)
    end

    return -loss_args[1] * A + loss_args[2] * Δv + loss_args[3] * outer + loss_args[4] * (E - 1)
end

function loss_poly2(X::Array, loss_args::Array{Float64, 1})
    global pices, silhouette

    # silhouetteとpieceの共通部分の面積を求める
    tmp_p = move(pieces[1], X[1], X[2])
    rotate!(tmp_p, X[3] * 2π)
    unioned = MYPolygon2LibGEOS(tmp_p)
    sum_pieces = LibGEOS.area(MYPolygon2LibGEOS(tmp_p))

    for i in 2:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        unioned = LibGEOS.union(unioned, MYPolygon2LibGEOS(tmp_p))
        sum_pieces += MYPolygon2LibGEOS(tmp_p) |> LibGEOS.area
    end
    F = LibGEOS.area(LibGEOS.intersection(unioned, silhouette)) / silhouette.area

    # 辺の誤差成分
    E = 0.0
    for P in pieces
        minimum_sinθ = 1
        for i in 1:P.n
            p1 = P.vertexes[:, i]
            p2 = i < n ? P.vertexes[:, i+1] : P.vertexes[:, 1]

            p21 = p2 - p1
            norm_p21 = norm(p21)

            # シルエットを見る
            for j in 1:silhouette.n
                s1 = silhouette.vertexes[:, j]
                s2 = j < silhouette.n ? silhouette.vertexes[:, j+1] : silhouette.vertexes[:, 1]
                sinθ = cross(p21, s2-s1)/(norm_p21*norm(s2-s1))
                minimum_sinθ = min(sinθ, minimum_sinθ)
            end
            # 他のピースを見る
            for Q in pieces
                if P == Q
                    continue
                end
                for j in 1:Q.n
                    q1 = Q.vertexes[j, :]
                    q2 = j < Q.n ? Q.vertexes[:, j+1] : Q.vertexes[:, 1]
                    sinθ = cross(p21, q2-q1)/(norm_p21*norm(q2-q1))
                    minimum_sinθ = min(sinθ, minimum_sinθ)
                end
            end
        end
        E += minimum_sinθ / length(pieces) * 2P.area
    end

    # 頂点の誤差成分
    V = 0.0
    for P in pieces
        v = 0.0
        for i in 1:P.n
            p = P.vertexes[:, i]

            min_d = 100.0
            # シルエットについて見る
            for j in 1:silhouette.n
                ps = norm(p - silhouette.vertexes[:, j])
                min_d = min(min_d, ps)
            end

            # 他のピースについて見る
            for Q in pieces
                if Q == P
                    continue
                end
                for j in 1:Q.n
                    pq = norm(p - Q.vertexes[:, j])
                    min_d = min(min_d, pq)
                end
            end
            v += min_d / P.n
        end
        V += 2 * P.area * v / length(pieces)
    end

    return loss_args[1] * F + loss_args[2] * E + loss_args[3] * V
end

# 面積ベースのピース批評
# シルエットから飛び出たやつしか扱えない
function reviewer_1(pieces::Array{MYPolygon{Float64}, 1}, X, silhouette::MYPolygon, threshold)
    # silhouetteと各pieceの共通部分の面積を求める
    score = zeros(Bool, length(pieces))

    for i in 1:length(pieces)
        tmp_p = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(tmp_p, X[3i] * 2π)
        GEOS_p = MYPolygon2LibGEOS(tmp_p)
        score[i] = threshold < intersectionArea(tmp_p, silhouette) / LibGEOS.area(GEOS_p)
    end
    return score
end

# 頂点ベースの批評
# おそらく、あらゆるケースを扱えるが計算量がエグそう
function reviewer_2(pieces::Array{MYPolygon{Float64}, 1}, X, silhouette::MYPolygon, threshold::T) where T <: Number
    # 各ピースの各頂点の中でもっともシルエットや他のピースの辺から離れているものを探す
    score = trues(length(pieces))

    for i in 1:length(pieces)
        optied_piece = move(pieces[i], X[3i-2], X[3i-1])
        rotate!(optied_piece, X[3i] * 2π)
        for p in 1:optied_piece.n
            P = optied_piece.vertexes[:, p]
            vertex_check = false

            # シルエットについて
            for s in 1:silhouette.n
                A = silhouette.vertexes[:, s]
                if s < silhouette.n
                    B = silhouette.vertexes[:, s+1]
                else
                    B = silhouette.vertexes[:, 1]
                end

                if distance_seg2point(A, B, P) < threshold
                    vertex_check = true
                    break
                end
            end

            # 他のピースについて
            for j in 1:length(pieces)
                if vertex_check
                    break
                elseif j == i
                    continue
                end
                optied_qiece = move(pieces[j], X[3j-2], X[3j-1])
                rotate!(optied_qiece, X[3j] * 2π)
                for q in 1:pieces[j].n
                    A = optied_qiece.vertexes[:, q]
                    if q < pieces[j].n
                        B = optied_qiece.vertexes[:, q+1]
                    else
                        B = optied_qiece.vertexes[:, 1]
                    end

                    if distance_seg2point(A, B, P) < threshold
                        vertex_check = true
                        break
                    end
                end
                if vertex_check
                    break
                end
            end

            # 一つでもダメな頂点があったらそのピースは終わり
            if vertex_check
            else
                score[i] = false
                break
            end
        end
    end
    return score
end

# p_Taiikusuwari は図形の頂点指定ミスでエラー起こる
puzzles = [p_house2, p_house4, p_arrow2, p_trapezoid2_1, p_trapezoid2_2, p_tri_above_square, p_fish2, p_small_stand, p_triforce3, p_rect_s, p_hurt3, p_hexagon_m1, p_hexagon_m2, p_hexagon5, p_square_3, p_square_4, p_square_7, p_butterfly_s, p_butterfly, p_pesudo_butterfly, p_butterfly4, p_note, p_diamond, p_antena, p_pencile4, p_crown, p_snake, p_faucet, p_crab_scissor, p_candy7, p_flame7, p_hurt7, p_katana]
loss_args = [100.0, 260.0, 0.0, 60.0]

# 初期化
max_gen = 128
sample_num = 128
λ = 18
recycling_durability = 4
recycle_threshold = 0.1

# 本体
rm("outputs", force=true, recursive=true)
mkdir("outputs")

for puzzle in puzzles
    global puzzle_name, silhouette, pieces = puzzle

    fpath = "outputs\\outputs_$puzzle_name"
    rm(fpath, force=true, recursive=true)
    mkdir(fpath)

    best_fitnesses = Array{Float64, 1}()
    p = Progress(sample_num, desc="$puzzle_name: ")
    for t in 1:sample_num
        cmaes = init_CMAES(zeros(3 * length(pieces)), 0.8, λ)
        fitnesses = zeros(cmaes.dim)
        X = zeros((cmaes.dim, cmaes.λ))
        Y = zeros((0, cmaes.λ))
        rng = MersenneTwister(t)

        for recycle in 1:recycling_durability
            for gen in 1:max_gen
                # 個体生成
                X = samplePopulation(cmaes, rng=rng)

                # 回転角のパラメータを -π ~ π までに正規化
                for j in 1:size(X)[2]
                    for i in 3:3:size(X)[1]
                        X[i, j] = mod(X[i, j], 1.0)
                    end
                end

                # 評価値
                loss(x) = loss_poly(x, loss_args)
                fitnesses = get_fitness(vcat(X, Y), loss)

                # 更新
                update!(cmaes, X, fitnesses, gen)
            end
            # 正解っぽいピースと不正解っぽいピースに分ける
            Z = vcat(X, Y)
            best_arg = argmin(fitnesses)
            score = reviewer_2(pieces, Z[:, best_arg], silhouette, recycle_threshold)

            new_pieces = Array{MYPolygon{Float64}, 1}()
            y = zeros(Float64, 0)
            x = zeros(Float64, 0)
            bad_num = 0
            for i in 1:length(pieces)
                if score[i]
                    push!(new_pieces, pieces[i])
                    push!(y, Z[3i-2:3i, best_arg]...)
                else
                    pushfirst!(new_pieces, pieces[i])
                    pushfirst!(x, Z[3i-2:3i, best_arg]...)
                    bad_num += 1
                end
            end
            if bad_num == 0 break end

            global pieces = deepcopy(new_pieces)
            Y = repeat(y, 1, λ)
            X = repeat(x, 1, λ)
            cmaes = init_CMAES(zeros(3 * bad_num), 1.0, λ)
        end
        push!(best_fitnesses, maximum(fitnesses))

        # 表示用
        best_arg = argmin(fitnesses)
        best_pieces = Array{MYPolygon, 1}()
        Z = vcat(X, Y)
        for i in 1:length(pieces)
            todisplay = move(pieces[i], Z[3i-2, best_arg], Z[3i-1, best_arg])
            rotate!(todisplay, Z[3i, best_arg] * 2π)
            push!(best_pieces, todisplay)
        end
        display(silhouette, best_pieces...)
        savefig(@sprintf "outputs\\outputs_%s\\%02.0f_%d_results" puzzle_name -best_fitnesses[end] t)

        next!(p)
    end
end

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@testset "Gherkin combinators  " begin
    @testset "Block text" begin
        @testset "Empty; OK" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                \"\"\"
            """)
            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == ""
        end

        @testset "Empty, then Quux; OK" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                \"\"\"
                Quux
            """)
            # Act
            p = Sequence{String}(BlockText(), Line("Quux"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["", "Quux"]
        end

        @testset "Foo; OK" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo
                \"\"\"
            """)
            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo"
        end

        @testset "Foo Bar Baz; OK" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo
                Bar
                Baz
                \"\"\"
            """)
            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo\nBar\nBaz"
        end

        @testset "Foo Bar Baz, then Quux; OK" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo
                Bar
                Baz
                \"\"\"
                Quux
            """)

            # Act
            p = Sequence{String}(BlockText(), Line("Quux"))
            result = p(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["Foo\nBar\nBaz", "Quux"]
        end

        @testset "Foo, empty line, then Baz; Empty line is included" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo

                Baz
                \"\"\"
            """)

            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo\n\nBaz"
        end

        @testset "Foo, comment, then Baz; Comment line is included" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo
                # Comment line
                Baz
                \"\"\"
            """)

            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo\n# Comment line\nBaz"
        end

        @testset "Foo, comment and empty, then Baz; Comment line and empty are included" begin
            # Arrange
            input = ParserInput("""
                \"\"\"
                Foo
                # Comment line

                Baz
                \"\"\"
            """)

            # Act
            p = BlockText()
            result = p(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "Foo\n# Comment line\n\nBaz"
        end
    end

    @testset "KeywordParser" begin
        @testset "Scenario:, Scenario:; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario:
            """)

            # Act
            parser = KeywordParser("Scenario:")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Keyword}
            @test result.value.keyword == "Scenario:"
            @test result.value.rest == ""
        end

        @testset "Given, Given; OK" begin
            # Arrange
            input = ParserInput("""
                Given
            """)

            # Act
            parser = KeywordParser("Given")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Keyword}
            @test result.value.keyword == "Given"
            @test result.value.rest == ""
        end

        @testset "Scenario:, Scenario Outline:; Not OK" begin
            # Arrange
            input = ParserInput("""
                Scenario Outline:
            """)

            # Act
            parser = KeywordParser("Scenario:")
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Keyword}
        end

        @testset "Scenario: then Scenario Outline:, Scenario:, Scenario Outline:; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario:
                Scenario Outline:
            """)

            # Act
            parser = Sequence{Keyword}(KeywordParser("Scenario:"), KeywordParser("Scenario Outline:"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{Keyword}}
            @test result.value[1].keyword == "Scenario:"
            @test result.value[1].rest == ""
            @test result.value[2].keyword == "Scenario Outline:"
            @test result.value[2].rest == ""
        end

        @testset "Scenario:, Scenario: Some description; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some description
            """)

            # Act
            parser = KeywordParser("Scenario:")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Keyword}
            @test result.value.keyword == "Scenario:"
            @test result.value.rest == "Some description"
        end

        @testset "Scenario:, Scenario:Some description, without a space; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario:Some description
            """)

            # Act
            parser = KeywordParser("Scenario:")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Keyword}
            @test result.value.keyword == "Scenario:"
            @test result.value.rest == "Some description"
        end

        @testset "Given; Given on the description as well; Only first given is removed" begin
            # Arrange
            input = ParserInput("""
                Given Given
            """)

            # Act
            parser = KeywordParser("Given")
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Keyword}
            @test result.value.keyword == "Given"
            @test result.value.rest == "Given"
        end
    end

    @testset "GivenParser" begin
        @testset "Given some precondition; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
            """)

            # Act
            parser = GivenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Given}
            @test result.value == Given("some precondition")
        end

        @testset "Given some other precondition; OK" begin
            # Arrange
            input = ParserInput("""
                Given some other precondition
            """)

            # Act
            parser = GivenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Given}
            @test result.value == Given("some other precondition")
        end

        @testset "When some action; Not OK" begin
            # Arrange
            input = ParserInput("""
                When some action
            """)

            # Act
            parser = GivenParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Given}
        end

        @testset "Followed by block text; OK and the text is present" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                    \"\"\"
                    Some block text.
                    On two lines.
                    \"\"\"
            """)

            # Act
            parser = GivenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Given}
            @test result.value isa Given
            @test result.value.text == "some precondition"
            @test result.value.block_text == "Some block text.\nOn two lines."
        end

        @testset "Givennospace; Not OK" begin
            # Arrange
            input = ParserInput("""
                Givennospace some precondition
            """)

            # Act
            parser = GivenParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Given}
        end
    end

    @testset "WhenParser" begin
        @testset "When some action; OK" begin
            # Arrange
            input = ParserInput("""
                When some action
            """)

            # Act
            parser = WhenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{When}
            @test result.value == When("some action")
        end

        @testset "When some other action; OK" begin
            # Arrange
            input = ParserInput("""
                When some other action
            """)

            # Act
            parser = WhenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{When}
            @test result.value == When("some other action")
        end

        @testset "Given some precondition; Not OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
            """)

            # Act
            parser = WhenParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{When}
        end

        @testset "Followed by block text; OK and the text is present" begin
            # Arrange
            input = ParserInput("""
                When some action
                    \"\"\"
                    Some block text.
                    On two lines.
                    \"\"\"
            """)

            # Act
            parser = WhenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{When}
            @test result.value isa When
            @test result.value.text == "some action"
            @test result.value.block_text == "Some block text.\nOn two lines."
        end

        @testset "Whennospace; Not OK" begin
            # Arrange
            input = ParserInput("""
                Whennospace some action
            """)

            # Act
            parser = WhenParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{When}
        end
    end

    @testset "ThenParser" begin
        # Tests for Given and When demonstrate correct behavior
        # and the design of the parser is such that this step will
        # have the same behavior, so I'm merely demonstrating the existence
        # of a ThenParser, not fully testing it.
        @testset "Then some postcondition; OK" begin
            # Arrange
            input = ParserInput("""
                Then some postcondition
            """)

            # Act
            parser = ThenParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Then}
            @test result.value == Then("some postcondition")
        end

        @testset "Thennospace; Not OK" begin
            # Arrange
            input = ParserInput("""
                Thennospace some postcondition
            """)

            # Act
            parser = ThenParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Then}
        end
    end

    @testset "Steps parser" begin
        @testset "Two givens; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                Given some other precondition
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] == Given("some precondition")
            @test result.value[2] == Given("some other precondition")
        end

        @testset "Two givens separated by a block text; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                    \"\"\"
                    Block text line 1.
                    Block text line 2.
                    \"\"\"
                Given some other precondition
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] isa Given
            @test result.value[1].text == "some precondition"
            @test result.value[1].block_text == "Block text line 1.\nBlock text line 2."
            @test result.value[2] == Given("some other precondition")
        end

        @testset "Two givens follow by a block text; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                Given some other precondition
                    \"\"\"
                    Block text line 1.
                    Block text line 2.
                    \"\"\"
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] == Given("some precondition")
            @test result.value[2] isa Given
            @test result.value[2].text == "some other precondition"
            @test result.value[2].block_text == "Block text line 1.\nBlock text line 2."
        end

        @testset "No givens; OK" begin
            # Arrange
            input = ParserInput("""
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value == []
        end

        @testset "Given then When; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                 When some action
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] == Given("some precondition")
            @test result.value[2] == When("some action")
        end

        @testset "When, Then; OK" begin
            # Arrange
            input = ParserInput("""
                 When some action
                 Then some postcondition
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] == When("some action")
            @test result.value[2] == Then("some postcondition")
        end

        @testset "Given, When, Then; OK" begin
            # Arrange
            input = ParserInput("""
                Given some precondition
                 When some action
                 Then some postcondition
            """)

            # Act
            parser = StepsParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{ScenarioStep}}
            @test result.value[1] == Given("some precondition")
            @test result.value[2] == When("some action")
            @test result.value[3] == Then("some postcondition")
        end
    end

    @testset "BackgroundParser" begin
        @testset "Empty Background, no description; OK" begin
            # Arrange
            input = ParserInput("""
                Background:
            """)

            # Act
            parser = BackgroundParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Background}
            @test result.value.description == ""
            @test result.value.steps == []
        end

        @testset "Empty Background, Some description; OK" begin
            # Arrange
            input = ParserInput("""
                Background: Some description
            """)

            # Act
            parser = BackgroundParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Background}
            @test result.value.description == "Some description"
            @test result.value.steps == []
        end

        @testset "Scenario:; Not OK" begin
            # Arrange
            input = ParserInput("""
                Scenario:
            """)

            # Act
            parser = BackgroundParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Background}
        end

        @testset "Given/When/Then; OK" begin
            # Arrange
            input = ParserInput("""
                Background: Some new description
                    Given some precondition
                     When some action
                     Then some precondition
            """)

            # Act
            parser = BackgroundParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Background}
            @test result.value.description == "Some new description"
            @test result.value.steps == [Given("some precondition"), When("some action"), Then("some precondition")]
        end
    end

    @testset "FeatureParser" begin
        @testset "Empty feature; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
        end

        @testset "Feature with two scenarios any steps; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario

                    Scenario: Some other scenario
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[2].description == "Some other scenario"
        end

        @testset "Feature with a scenario; Default background" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.background.description == ""
            @test result.value.background.steps == []
        end

        @testset "Feature and Rule with one scenario; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: A feature description

                    Rule: Some rule description

                        Scenario: Some scenario
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "A feature description"
            @test result.value.scenarios[1].description == "Some rule description"
            rule = result.value.scenarios[1]
            @test rule.scenarios[1].description == "Some scenario"
        end

        @testset "Feature with a background; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Background: Some background
                        Given some precondition
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.background.description == "Some background"
            @test result.value.background.steps == [Given("some precondition")]
            @test result.value.scenarios == []
        end

        @testset "Feature with a background, then a Scenario; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Background: Some background
                        Given some precondition

                    Scenario: Some scenario
                        When some action
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.background.description == "Some background"
            @test result.value.background.steps == [Given("some precondition")]
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[1].steps == [When("some action")]
        end

        @testset "Feature with a scenario outline; Scenario outline is the only step" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario Outline: Some scenario
                        Given some precondition <Foo>

                        Examples:
                            | Foo |
                            | bar |
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            outline = result.value.scenarios[1]
            @test outline.description == "Some scenario"
            @test outline.placeholders == ["Foo"]
            @test outline.examples == [["bar"]]
        end
    end

    @testset "FeatureFileParser" begin
        @testset "Feature with two scenarios any steps; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario

                    Scenario: Some other scenario
            """)

            # Act
            parser = FeatureFileParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[2].description == "Some other scenario"
        end

        @testset "Feature, then unallowed new Feature; Not OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    Scenario: Some scenario

                    Scenario: Some other scenario

                    Feature: Not allowed here
            """)

            # Act
            parser = FeatureFileParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Feature}
        end
    end

    @testset "DataTable" begin
        @testset "One row, one column; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Foo"]]
        end

        @testset "One row, one column; OK" begin
            # Arrange
            input = ParserInput("""
                | Bar |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Bar"]]
        end

        @testset "Table row then line Baz; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo |
                Baz
            """)

            # Act
            parser = Sequence{Union{DataTable, String}}(DataTableParser(), Line("Baz"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{Union{DataTable, String}}}
            @test result.value[1] == [["Foo"]]
            @test result.value[2] == "Baz"
        end

        @testset "No pipes on the line; Not OK" begin
            # Arrange
            input = ParserInput("""
                Baz
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{DataTable}
        end

        @testset "Two columns; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo | Bar |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Foo", "Bar"]]
        end

        @testset "Three columns; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo | Bar | Baz |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Foo", "Bar", "Baz"]]
        end

        @testset "Two rows; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo |
                | Bar |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Foo"], ["Bar"]]
        end

        @testset "EOF; Not OK" begin
            # Arrange
            input = ParserInput("""
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{DataTable}
        end

        @testset "Many columns and rows; OK" begin
            # Arrange
            input = ParserInput("""
                | Foo | Bar |
                | Baz | Quux |
            """)

            # Act
            parser = DataTableParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{DataTable}
            @test result.value == [["Foo", "Bar"], ["Baz", "Quux"]]
        end
    end

    @testset "TagsParser" begin
        @testset "AnyLine; @tag; OK" begin
            # Arrange
            input = ParserInput("@tag")

            # Act
            parser = Experimental.AnyLine()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{String}
            @test result.value == "@tag"
        end

        @testset "Splitter; @tag1 and @tag2; OK" begin
            # Arrange
            input = ParserInput("@tag1 @tag2")

            # Act
            parser = Experimental.Splitter(Experimental.AnyLine(), isspace)
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2"]
        end

        @testset "Validator; All strings begin with an @; OK" begin
            # Arrange
            input = ParserInput("@tag1 @tag2")

            # Act
            inner = Experimental.Splitter(Experimental.AnyLine(), isspace)
            parser = Experimental.Validator{String}(inner, x -> startswith(x, "@"))
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2"]
        end

        @testset "Validator; Not all strings begin with an @; Not OK" begin
            # Arrange
            input = ParserInput("@tag1 tag2")

            # Act
            inner = Experimental.Splitter(Experimental.AnyLine(), isspace)
            parser = Experimental.Validator{String}(inner, x -> startswith(x, "@"))
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "@tag; OK" begin
            # Arrange
            input = ParserInput("""
                @tag
            """)

            # Act
            parser = TagParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag"]
        end

        @testset "tag; Not OK" begin
            # Arrange
            input = ParserInput("""
                tag
            """)

            # Act
            parser = TagParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "EOF; Not OK" begin
            # Arrange
            input = ParserInput("")

            # Act
            parser = TagParser()
            result = parser(input)

            # Assert
            @test result isa BadParseResult{Vector{String}}
        end

        @testset "@tag followed by Scenario:; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1
                Scenario: Some scenario
            """)

            # Act
            parser = Sequence{Union{Vector{String}, Keyword}}(
                        TagParser(),
                        KeywordParser("Scenario:")
            )
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{Union{Vector{String}, Keyword}}}
            @test result.value[1] == ["@tag1"]
            @test result.value[2] == Keyword("Scenario:", "Some scenario")
        end

        @testset "@tag1 @tag2; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1 @tag2
            """)

            # Act
            parser = TagParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2"]
        end

        @testset "@tag1 @tag2 with multiple spaces; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1     @tag2
            """)

            # Act
            parser = TagParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2"]
        end

        @testset "@tag1 then @tag2; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1
                @tag2
            """)

            # Act
            parser = TagLinesParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2"]
        end

        @testset "@tag1 @tag2, then @tag3; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1 @tag2
                @tag3
            """)

            # Act
            parser = TagLinesParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Vector{String}}
            @test result.value == ["@tag1", "@tag2", "@tag3"]
        end
    end

    @testset "Tags" begin
        @testset "Empty Scenario, with tags; OK" begin
            # Arrange
            input = ParserInput("""
                @tag1
                @tag2
                Scenario: Some description
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some description"
            @test result.value.tags == ["@tag1", "@tag2"]
        end

        @testset "Empty feature with tags; OK" begin
            # Arrange
            input = ParserInput("""
                @sometag
                @othertag
                Feature: Some feature
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.tags == ["@sometag", "@othertag"]
        end

        @testset "Empty Rule with tags; OK" begin
            # Arrange
            input = ParserInput("""
                @sometag
                @othertag
                Rule: Some rule
            """)

            # Act
            parser = RuleParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Rule}
            @test result.value.tags == ["@sometag", "@othertag"]
        end
    end

    @testset "Long descriptions parser" begin
        @testset "Description Foo; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some new description
                    Foo
                    Given some precondition
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some new description"
            @test result.value.long_description == "Foo"
        end

        @testset "Description Bar; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some new description
                    Bar
                    Given some precondition
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some new description"
            @test result.value.long_description == "Bar"
        end

        @testset "Description Bar then When; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some new description
                    Bar
                    When some action
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some new description"
            @test result.value.long_description == "Bar"
        end

        @testset "Description Bar then Then; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some new description
                    Bar
                    Then some postcondition
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some new description"
            @test result.value.long_description == "Bar"
        end
    end

    @testset "Long descriptions" begin
        @testset "Description on a feature, no Background; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    This is a description.
                    On two lines.

                    Scenario: Some scenario
                        When some action
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.header.long_description == ["This is a description.\nOn two lines."]
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[1].steps == [When("some action")]
        end

        @testset "Description on a feature, with Background; OK" begin
            # Arrange
            input = ParserInput("""
                Feature: Some feature

                    This is a description.
                    On two lines.

                    Background: Some background
                        Given some precondition

                    Scenario: Some scenario
                        When some action
            """)

            # Act
            parser = FeatureParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Feature}
            @test result.value.header.description == "Some feature"
            @test result.value.header.long_description == ["This is a description.\nOn two lines."]
            @test result.value.background.description == "Some background"
            @test result.value.background.steps == [Given("some precondition")]
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[1].steps == [When("some action")]
        end

        @testset "Scenario with a description; OK" begin
            # Arrange
            input = ParserInput("""
                Scenario: Some new description
                    Foo
                    Bar
                    Baz
                    Given some precondition
                    Given some other precondition
            """)

            # Act
            parser = ScenarioParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Scenario}
            @test result.value.description == "Some new description"
            @test result.value.long_description == "Foo\nBar\nBaz"
        end

        @testset "Background with a description; OK" begin
            # Arrange
            input = ParserInput("""
                Background: Some background description

                    This is a description.
                    On two lines.

                    Given some precondition
                    Given some other precondition
            """)

            # Act
            parser = BackgroundParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Background}
            @test result.value.description == "Some background description"
            @test result.value.long_description == "This is a description.\nOn two lines."
        end

        @testset "Description on a Rule; OK" begin
            # Arrange
            input = ParserInput("""
                Rule: Some rule

                    This is a description.
                    On two lines.

                    Scenario: Some scenario
                        When some action
            """)

            # Act
            parser = RuleParser()
            result = parser(input)

            # Assert
            @test result isa OKParseResult{Rule}
            @test result.value.description == "Some rule"
            @test result.value.longdescription == "This is a description.\nOn two lines."
            @test result.value.scenarios[1].description == "Some scenario"
            @test result.value.scenarios[1].steps == [When("some action")]
        end

        # TODO Scenario Outline
    end
end


text:
module EvolutionaryModelingTools

using RecordedArrays
using Random
using MacroTools

export @cfunc, @ufunc, @reaction
export Reaction, gillespie, gillespie!

include("tools.jl")
include("sample.jl")
include("model.jl")

end # module


text:
calcMidpoints(edges::AbstractVector) = Float64[0.5 * (edges[i] + edges[i+1]) for i in 1:length(edges)-1]

"Make histogram-like bins of data"
function binData(data, nbins)
  lo, hi = extrema(data)
  edges = collect(linspace(lo, hi, nbins+1))
  midpoints = calcMidpoints(edges)
  buckets = Int[max(2, min(searchsortedfirst(edges, x), length(edges)))-1 for x in data]
  counts = zeros(Int, length(midpoints))
  for b in buckets
    counts[b] += 1
  end
  edges, midpoints, buckets, counts
end

"""
A hacky replacement for a histogram when the backend doesn't support histograms directly.
Convert it into a bar chart with the appropriate x/y values.
"""
function histogramHack(; kw...)
  d = Dict(kw)

  # we assume that the y kwarg is set with the data to be binned, and nbins is also defined
  edges, midpoints, buckets, counts = binData(d[:y], d[:nbins])
  d[:x] = midpoints
  d[:y] = float(counts)
  d[:linetype] = :bar
  d[:fillrange] = d[:fillrange] == nothing ? 0.0 : d[:fillrange]
  d
end

"""
A hacky replacement for a bar graph when the backend doesn't support bars directly.
Convert it into a line chart with fillrange set.
"""
function barHack(; kw...)
  d = Dict(kw)
  midpoints = d[:x]
  heights = d[:y]
  fillrange = d[:fillrange] == nothing ? 0.0 : d[:fillrange]

  # estimate the edges
  dists = diff(midpoints) * 0.5
  edges = zeros(length(midpoints)+1)
  for i in 1:length(edges)
    if i == 1
      edge = midpoints[1] - dists[1]
    elseif i == length(edges)
      edge = midpoints[i-1] + dists[i-2]
    else
      edge = midpoints[i-1] + dists[i-1]
    end
    edges[i] = edge
  end

  x = Float64[]
  y = Float64[]
  for i in 1:length(heights)
    e1, e2 = edges[i:i+1]
    append!(x, [e1, e1, e2, e2])
    append!(y, [fillrange, heights[i], heights[i], fillrange])
  end

  d[:x] = x
  d[:y] = y
  d[:linetype] = :path
  d[:fillrange] = fillrange
  d
end

"""
A hacky replacement for a sticks graph when the backend doesn't support sticks directly.
Convert it into a line chart that traces the sticks, and a scatter that sets markers at the points.
"""
function sticksHack(; kw...)
  dLine = Dict(kw)
  dScatter = copy(dLine)

  # these are the line vertices
  x = Float64[]
  y = Float64[]
  fillrange = dLine[:fillrange] == nothing ? 0.0 : dLine[:fillrange]

  # calculate the vertices
  yScatter = dScatter[:y]
  for (i,xi) in enumerate(dScatter[:x])
    yi = yScatter[i]
    for j in 1:3 push!(x, xi) end
    append!(y, [fillrange, yScatter[i], fillrange])
  end

  # change the line args
  dLine[:x] = x
  dLine[:y] = y
  dLine[:linetype] = :path
  dLine[:markershape] = :none
  dLine[:fillrange] = nothing

  # change the scatter args
  dScatter[:linetype] = :none

  dLine, dScatter
end

function regressionXY(x, y)
  # regress
  β, α = convert(Matrix{Float64}, [x ones(length(x))]) \ convert(Vector{Float64}, y)

  # make a line segment
  regx = [minimum(x), maximum(x)]
  regy = β * regx + α
  regx, regy
end

# ------------------------------------------------------------------------------------

nop() = nothing

get_mod(v::AVec, idx::Int) = v[mod1(idx, length(v))]
get_mod(v::AMat, idx::Int) = size(v,1) == 1 ? v[1, mod1(idx, size(v,2))] : v[:, mod1(idx, size(v,2))]
get_mod(v, idx::Int) = v

makevec(v::AVec) = v
makevec{T}(v::T) = T[v]

"duplicate a single value, or pass the 2-tuple through"
maketuple(x::Real) = (x,x)
maketuple{T,S}(x::@compat(Tuple{T,S})) = x

unzip{T,S}(v::AVec{@compat(Tuple{T,S})}) = [vi[1] for vi in v], [vi[2] for vi in v]

# given 2-element lims and a vector of data x, widen lims to account for the extrema of x
function _expand_limits(lims, x)
  try
    e1, e2 = extrema(x)
    lims[1] = min(lims[1], e1)
    lims[2] = max(lims[2], e2)
  # catch err
  #   warn(err)
  end
  nothing
end

# if the type exists in a list, replace the first occurence.  otherwise add it to the end
function addOrReplace(v::AbstractVector, t::DataType, args...; kw...)
    for (i,vi) in enumerate(v)
        if isa(vi, t)
            v[i] = t(args...; kw...)
            return
        end
    end
    push!(v, t(args...; kw...))
    return
end

function replaceType(vec, val)
  filter!(x -> !isa(x, typeof(val)), vec)
  push!(vec, val)
end

function replaceAliases!(d::Dict, aliases::Dict)
  ks = collect(keys(d))
  for k in ks
    if haskey(aliases, k)
      d[aliases[k]] = d[k]
      delete!(d, k)
    end
  end
end

createSegments(z) = collect(repmat(z',2,1))[2:end]

Base.first(c::Colorant) = c
Base.first(x::Symbol) = x

sortedkeys(d::Dict) = sort(collect(keys(d)))

function fakedata(sz...)
  y = zeros(sz...)
  for r in 2:size(y,1)
    y[r,:] = 0.95 * y[r-1,:] + randn(size(y,2))'
  end
  y
end

isijulia() = isdefined(Main, :IJulia) && Main.IJulia.inited

istuple(::Tuple) = true
istuple(::Any) = false
isvector(::AVec) = true
isvector(::Any) = false
ismatrix(::AMat) = true
ismatrix(::Any) = false
isscalar(::Real) = true
isscalar(::Any) = false

# ticksType{T<:Real,S<:Real}(ticks::@compat(Tuple{T,S})) = :limits
ticksType{T<:Real}(ticks::AVec{T}) = :ticks
ticksType{T<:AVec,S<:AVec}(ticks::@compat(Tuple{T,S})) = :ticks_and_labels
ticksType(ticks) = :invalid

limsType{T<:Real,S<:Real}(lims::@compat(Tuple{T,S})) = :limits
limsType(lims::Symbol) = lims == :auto ? :auto : :invalid
limsType(lims) = :invalid

Base.convert{T<:Real}(::Type{Vector{T}}, rng::Range{T}) = T[x for x in rng]
Base.convert{T<:Real,S<:Real}(::Type{Vector{T}}, rng::Range{S}) = T[x for x in rng]

Base.merge(a::AbstractVector, b::AbstractVector) = sort(unique(vcat(a,b)))

# ---------------------------------------------------------------

"""
Allows temporary setting of backend and defaults for Plots. Settings apply only for the `do` block.  Example:
```
with(:gadfly, size=(400,400), type=:hist) do
  plot(rand(10))
  plot(rand(10))
end
```
"""
function with(f::Function, args...; kw...)

  # dict to store old and new keyword args for anything that changes
  newdefs = Dict(kw)
  olddefs = Dict()
  for k in keys(newdefs)
    olddefs[k] = default(k)
  end

  # save the backend
  if CURRENT_BACKEND.sym == :none
    pickDefaultBackend()
  end
  oldbackend = CURRENT_BACKEND.sym

  for arg in args

    # change backend?
    if arg in backends()
      backend(arg)
    end

    # # TODO: generalize this strategy to allow args as much as possible
    # #       as in:  with(:gadfly, :scatter, :legend, :grid) do; ...; end
    # # TODO: can we generalize this enough to also do something similar in the plot commands??

    # k = :linetype
    # if arg in _allTypes
    #   olddefs[k] = default(k)
    #   newdefs[k] = arg
    # elseif haskey(_typeAliases, arg)
    #   olddefs[k] = default(k)
    #   newdefs[k] = _typeAliases[arg]
    # end

    k = :legend
    if arg in (k, :leg)
      olddefs[k] = default(k)
      newdefs[k] = true
    end

    k = :grid
    if arg == k
      olddefs[k] = default(k)
      newdefs[k] = true
    end
  end

  # display(olddefs)
  # display(newdefs)

  # now set all those defaults
  default(; newdefs...)

  # call the function
  ret = f()

  # put the defaults back
  default(; olddefs...)

  # revert the backend
  if CURRENT_BACKEND.sym != oldbackend
    backend(oldbackend)
  end

  # return the result of the function
  ret
end

# ---------------------------------------------------------------

type DebugMode
  on::Bool
end
const _debugMode = DebugMode(false)

function debugplots(on = true)
  _debugMode.on = on
end

debugshow(x) = show(x)
debugshow(x::AbstractArray) = print(summary(x))

function dumpdict(d::Dict, prefix = "", alwaysshow = false)
  _debugMode.on || alwaysshow || return
  println()
  if prefix != ""
    println(prefix, ":")
  end
  for k in sort(collect(keys(d)))
    @printf("%14s: ", k)
    debugshow(d[k])
    println()
  end
  println()
end

function dumpcallstack()
  error()  # well... you wanted the stacktrace, didn't you?!?
end

# ---------------------------------------------------------------

# push/append/clear/set the underlying plot data
# NOTE: backends should implement the getindex and setindex! methods to get/set the x/y data objects

# index versions
function Base.push!(plt::Plot, i::Integer, x::Real, y::Real)
  xdata, ydata = plt[i]
  plt[i] = (extendSeriesData(xdata, x), extendSeriesData(ydata, y))
  plt
end
function Base.push!(plt::Plot, i::Integer, y::Real)
  xdata, ydata = plt[i]
  # if !isa(xdata, UnitRange)
  #   error("Expected x is a UnitRange since you're trying to push a y value only. typeof(x) = $(typeof(xdata))")
  # end
  plt[i] = (extendSeriesByOne(xdata), extendSeriesData(ydata, y))
  plt
end

Base.push!(plt::Plot, y::Real) = push!(plt, 1, y)

# update all at once
function Base.push!(plt::Plot, x::AVec, y::AVec)
  nx = length(x)
  ny = length(y)
  for i in 1:plt.n
    push!(plt, i, x[mod1(i,nx)], y[mod1(i,ny)])
  end
  plt
end

function Base.push!(plt::Plot, x::Real, y::AVec)
  push!(plt, [x], y)
end

function Base.push!(plt::Plot, y::AVec)
  ny = length(y)
  for i in 1:plt.n
    push!(plt, i, y[mod1(i,ny)])
  end
  plt
end

# append to index
function Base.append!(plt::Plot, i::Integer, x::AVec, y::AVec)
  @assert length(x) == length(y)
  xdata, ydata = plt[i]
  plt[i] = (extendSeriesData(xdata, x), extendSeriesData(ydata, y))
  plt
end

function Base.append!(plt::Plot, i::Integer, y::AVec)
  xdata, ydata = plt[i]
  if !isa(xdata, UnitRange{Int})
    error("Expected x is a UnitRange since you're trying to push a y value only")
  end
  plt[i] = (extendSeriesByOne(xdata, length(y)), extendSeriesData(ydata, y))
  plt
end

# used in updating an existing series

extendSeriesByOne(v::UnitRange{Int}, n::Int = 1) = isempty(v) ? (1:n) : (minimum(v):maximum(v)+n)
extendSeriesByOne(v::AVec, n::Integer = 1) = isempty(v) ? (1:n) : vcat(v, (1:n) + maximum(v))
extendSeriesData{T}(v::Range{T}, z::Real) = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::Range{T}, z::AVec) = extendSeriesData(float(collect(v)), z)
extendSeriesData{T}(v::AVec{T}, z::Real) = (push!(v, convert(T, z)); v)
extendSeriesData{T}(v::AVec{T}, z::AVec) = (append!(v, convert(Vector{T}, z)); v)

# ---------------------------------------------------------------

function supportGraph(allvals, func)
  vals = reverse(sort(allvals))
  bs = sort(backends())
  x = ASCIIString[]
  y = ASCIIString[]
  for val in vals
    for b in bs
        supported = func(Plots.backendInstance(b))
        if val in supported
            push!(x, string(b))
            push!(y, string(val))
        end
      end
  end
  n = length(vals)

  scatter(x,y,
          m=:rect,
          ms=10,
          size=(300,100+18*n),
          # xticks=(collect(1:length(bs)), bs),
          leg=false
         )
end

supportGraphArgs() = supportGraph(_allArgs, supportedArgs)
supportGraphTypes() = supportGraph(_allTypes, supportedTypes)
supportGraphStyles() = supportGraph(_allStyles, supportedStyles)
supportGraphMarkers() = supportGraph(_allMarkers, supportedMarkers)
supportGraphScales() = supportGraph(_allScales, supportedScales)
supportGraphAxes() = supportGraph(_allAxes, supportedAxes)

function dumpSupportGraphs()
  for func in (supportGraphArgs, supportGraphTypes, supportGraphStyles,
               supportGraphMarkers, supportGraphScales, supportGraphAxes)
    plt = func()
    png(IMG_DIR * "/supported/$(string(func))")
  end
end

# ---------------------------------------------------------------

# Some conversion functions
# note: I borrowed these conversion constants from Compose.jl's Measure
const PX_PER_INCH = 100
const DPI = PX_PER_INCH
const MM_PER_INCH = 25.4
const MM_PER_PX = MM_PER_INCH / PX_PER_INCH
inch2px(inches::Real) = float(inches * PX_PER_INCH)
px2inch(px::Real) = float(px / PX_PER_INCH)
inch2mm(inches::Real) = float(inches * MM_PER_INCH)
mm2inch(mm::Real) = float(mm / MM_PER_INCH)
px2mm(px::Real) = float(px * MM_PER_PX)
mm2px(mm::Real) = float(px / MM_PER_PX)


text:
plan_clenshawcurtis(μ) = length(μ) > 1 ? FFTW.plan_r2r!(μ, FFTW.REDFT00) : fill!(similar(μ),1)'

"""
Compute nodes of the Clenshaw—Curtis quadrature rule.
"""
clenshawcurtisnodes(::Type{T}, N::Int) where T = chebyshevpoints(T, N; kind = 2)

"""
Compute weights of the Clenshaw—Curtis quadrature rule with modified Chebyshev moments of the first kind ``\\mu``.
"""
clenshawcurtisweights(μ::Vector) = clenshawcurtisweights!(copy(μ))
clenshawcurtisweights!(μ::Vector) = clenshawcurtisweights!(μ, plan_clenshawcurtis(μ))
function clenshawcurtisweights!(μ::Vector{T}, plan) where T
    N = length(μ)
    rmul!(μ, inv(N-one(T)))
    plan*μ
    μ[1] *= half(T); μ[N] *= half(T)
    return μ
end


text:
using DataStructures: CircularBuffer

"""
    FixedWindowAssociativeOp{T,Op,Op!}(window)
    FixedWindowAssociativeOp{T,Op}(window)

State necessary for accumulation over a rolling window of fixed size.

# Fields
- `window_state::WindowedAssociativeOp{T,Op}`: The underlying general-window state.
- `remaining_window::Int`: How much of the window remains to be filled. Initially this will
    be set to the window size, and will then reduce for every value added until it reaches
    zero.
"""
mutable struct FixedWindowAssociativeOp{T,Op,Op!}
    window_state::WindowedAssociativeOp{T,Op,Op!,CircularBuffer{T}}
    remaining_window::Int

    function FixedWindowAssociativeOp{T,Op,Op!}(
        window::Integer
    ) where {T,Op,Op!}
        window < 1 && throw(ArgumentError("Got window $window, but it must be positive."))
        window_state = WindowedAssociativeOp{T,Op,Op!,CircularBuffer{T}}(
            CircularBuffer{T}(window - 1), CircularBuffer{T}(window)
        )
        return new(window_state, window)
    end
end

function FixedWindowAssociativeOp{T,Op}(window) where {T,Op}
    return FixedWindowAssociativeOp{T,Op,Op}(window)
end

"""
    update_state!(state::FixedWindowAssociativeOp, value) -> state

Add the specified `value` to the `state`. Drop a value from the window iff the window is
full.

# Returns
- `::FixedWindowAssociativeOp`: The instance `state` that was passed in.
"""
function update_state!(state::FixedWindowAssociativeOp, value)
    num_dropped_from_window = if state.remaining_window > 0
        state.remaining_window -= 1
        0
    else
        1
    end

    # With @inbounds, we assert that num_dropped_from_window will never exceed the size of
    # the window.
    @inbounds update_state!(state.window_state, value, num_dropped_from_window)
    return state
end

window_value(state::FixedWindowAssociativeOp) = window_value(state.window_state)
window_size(state::FixedWindowAssociativeOp)::Int = window_size(state.window_state)

"""
    window_full(state::FixedWindowAssociativeOp)::Bool

Returns true iff the given `state` has a full window.
"""
window_full(state::FixedWindowAssociativeOp)::Bool = state.remaining_window == 0


text:
module ACP_MC

include("../Misc.jl")
include("BackwardKernel.jl")
using .Misc, .BackwardKernel

using Distributions, SpecialFunctions
import Base.Iterators: product

export Particle, smc_weight

mutable struct Particle{ℜ <: Real}
    S_R::Array{ℜ,1}
    S_I::Array{Array{ℜ,2}}
    allocations
    function Particle(X::Array{ℜ}, R::Ƶ) where {ℜ<:Real, Ƶ<:Int}
        return new{ℜ}(zeros(ℜ,R), [zeros(ℜ,R,Iₙ) for Iₙ ∈ size(X)], [Tuple(zeros(length(size(X))+2)) for i in 1:sum(X)])
    end
    function Particle(S_R::Array{ℜ,1}, S_I::Array{Array{ℜ,2}}) where {ℜ<:Real}
        return new{ℜ}(S_R, S_I, [Tuple(zeros(length(S_I)+2)) for i in 1:sum(S_R)])
    end
end

function smc_weight(X::Array{ℜ,N}, R::Ƶ, P::Ƶ=1; a=1.0, b=a/nansum(X), debug=true) where {ℜ<:Real, Ƶ<:Integer,N}
    X_full = deepcopy(X)
    X = sum(X_full, N)
    I = size(X)

    α_I = (a/R) ./ I
    α_R = a/R

    T = Ƶ(round(sum(X)))
    ESS = zeros(T)

    Π::Array{Particle} = [Particle(X, R) for p=1:P]
    Π₂::Array{Particle} = [Particle(X, R) for p=1:P]

    log_Z = a*log(b) - (a+T)*log(b + 1) + lgamma(a + T) - lgamma(a) - sum(lgamma, X .+ 1)

    log_w, W, cum_W = fill(log_Z,P), fill(1.0/P,P), zeros(P)
    log_q, q = zeros(R), zeros(R)
    r = 1
    for (t,i) in enumerate(EventQueue(X_full))
        z = i[end]
        i = i[1:end-1]
        for (p,π) ∈ enumerate(Π)

            log_q .=  sum(n -> log.(π.S_I[n][:,i[n]] .+ α_I[n]),1:(N-1)) - (N-1-1)*log.(π.S_R .+ α_R) .- log(t+a-1)
            log_ν = logsumexp(log_q)
            log_q .-= log_ν

            q .= exp.(log_q)
            r = rand(Categorical(q))

            π.S_R[r] += 1
            for n ∈ 1:(N-1)
                π.S_I[n][r,i[n]] += 1
            end

            log_w[p] += log_ν
            π.allocations[t] = Tuple([i...,r,z])
        end

        log_Z = logmeanexp(log_w)
        W .= exp.(log_w .- logsumexp(log_w))

        ESS[t] = 1.0/sum(W .* W)

        if ESS[t] < P/2
            cum_W .= cumsum(W)
            u = rand()/P

            for p ∈ 1:P
                Π₂[p].S_R .= Π[p].S_R
                for n ∈ 1:(N-1)
                    Π₂[p].S_I[n] .= Π[p].S_I[n]
                end
            end

            p₂ = 1
            for p ∈ 1:P # systematic resampling
                while cum_W[p₂] < u
                    p₂ += 1
                end
                Π[p].S_R .= Π₂[p₂].S_R
                for n ∈ 1:(N-1)
                    Π[p].S_I[n] .= Π₂[p₂].S_I[n]
                end
                u += 1.0/P
            end
            log_w .= log_Z
        end
    end
    p = rand(Categorical(W))
    return log_Z, ESS, Π[p]
end
end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
using Composite
using Base.Test

twice  = (+) ∘ dup
@test twice(  3 ) == 6

square = (*) ∘ dup
@test square( 3 ) == 9

twice_square = twice  ∘ square
@test twice_square( 3 ) == 18

square_twice = square ∘ twice
@test square_twice( 3 ) == 36

twice_and_square = twice ^ square
@test twice_and_square( 3 ) == ( 6, 9 )

twice_first_and_square_second  = twice % square
@test twice_first_and_square_second( 3, 4 ) == (6, 16)

add3 = (+) < 3
@test add3() == 3
@test add3( 4 ) == 7

gt3 = (>) > 3
@test ! gt3( 2 )
@test ! gt3( 3 )
@test gt3( 4 )

max_stars_repo_name: sabjohnso/Composite

max_issues_repo_name: sabjohnso/Composite

max_forks_repo_name: sabjohnso/Composite


text:
struct KeldyshTimeGF{T, scalar} <: AbstractTimeGF{T, scalar}
  grid::KeldyshTimeGrid
  gtr::AntiHermitianStorage{T,scalar}
  les::AntiHermitianStorage{T,scalar}
  ξ::GFSignEnum
end

norbitals(G::KeldyshTimeGF) = G.gtr.norb

function KeldyshTimeGF(::Type{T}, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  nt = grid.nt

  gtr = AntiHermitianStorage(T, nt, norb, scalar)
  les = AntiHermitianStorage(T, nt, norb, scalar)

  KeldyshTimeGF(grid, gtr, les, ξ)
end
KeldyshTimeGF(grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = KeldyshTimeGF(ComplexF64, grid, norb, ξ, scalar)

@inline function Base.getindex(G::KeldyshTimeGF, k, l, t1::TimeGridPoint, t2::TimeGridPoint, greater=true)
  greater = t1 == t2 ? greater : heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx

  return greater ? G.gtr[k,l,i,j] : G.les[k,l,i,j]
end

function Base.setindex!(G::KeldyshTimeGF, v, k, l, t1::TimeGridPoint, t2::TimeGridPoint)
  greater = heaviside(t1.bpoint, t2.bpoint)

  i = t1.ridx
  j = t2.ridx

  return greater ? G.gtr[k,l,i,j] = v : G.les[k,l,i,j] = v
end

function Base.similar(G::T, ξ::GFSignEnum) where T <: KeldyshTimeGF
  T(G.grid, similar(G.gtr), similar(G.les), ξ)
end

function Base.zero(G::T, ξ::GFSignEnum) where T <: KeldyshTimeGF
  T(G.grid, zero(G.gtr), zero(G.les), ξ)
end

function Base.similar(G::T) where T <: KeldyshTimeGF
  similar(G, G.ξ)
end

function Base.zero(G::T) where T <: KeldyshTimeGF
  zero(G, G.ξ)
end

function TimeDomain(G::KeldyshTimeGF)
  grid = G.grid
  tplus = grid[forward_branch]
  tminus = reverse(grid[backward_branch])

  nt = grid.nt

  points = Tuple{TimeGridPoint, TimeGridPoint}[]

  for i in 1:nt
    for j in 1:i
      push!(points, (tminus[j], tminus[i])) # greater
      push!(points, (tplus[j], tminus[i])) # lesser
    end
  end

  return TimeDomain(points)
end

function KeldyshTimeGF(f::Function, ::Type{T}, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) where T <: Number
  G = KeldyshTimeGF(T, grid, norb, ξ, scalar)
  D = TimeDomain(G)

  for (t1, t2) in D.points
    G[t1, t2] = f(t1, t2)
  end

  return G
end

KeldyshTimeGF(f::Function, grid::KeldyshTimeGrid, norb=1, ξ::GFSignEnum=fermionic, scalar=false) = KeldyshTimeGF(f, ComplexF64, grid, norb, ξ, scalar)

function KeldyshTimeGF(dos::AbstractDOS, β, grid::KeldyshTimeGrid)
  # result is time invariant so do calculation for time invariant case and then copy over
  G = TimeInvariantKeldyshTimeGF(dos, β, grid)
  KeldyshTimeGF(grid,1,fermionic,true) do t1, t2
    G[t1,t2]
  end
end


text:
function ground(args...)
    background("white")
    sethue("black")
end

@testset "Javis Viewer" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    l1 = @JLayer 20:60 100 100 Point(0, 0) begin
        obj = Object((args...) -> circle(O, 25, :fill))
        act!(obj, Action(1:20, appear(:fade)))
    end

    render(vid; pathname = "")

    action_list = [back, star_obj]

    viewer_win, frame_dims, r_slide, tbox, canvas, actions, total_frames, video =
        Javis._javis_viewer(vid, 100, action_list, false)
    visible(viewer_win, false)

    @test get_gtk_property(viewer_win, :title, String) == "Javis Viewer"

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    second_frame = Javis.get_javis_frame(video, actions, curr_frame, layers = [l1])
    @test Reactive.value(r_slide) == 2

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame, layers = [l1])
    @test Reactive.value(r_slide) == 1

    @test first_frame != second_frame

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    last_frame = Javis.get_javis_frame(video, actions, curr_frame, layers = [l1])
    @test curr_frame == total_frames

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame, layers = [l1])
    @test curr_frame == 1

    @test last_frame != first_frame
end

@testset "Livestreaming" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    conf_local = setup_stream(:local, address = "0.0.0.0", port = 8081)
    @test conf_local isa Javis.StreamConfig
    @test conf_local.livestreamto == :local
    @test conf_local.protocol == "udp"
    @test conf_local.address == "0.0.0.0"
    @test conf_local.port == 8081

    conf_twitch_err = setup_stream(:twitch)
    conf_twitch = setup_stream(:twitch, twitch_key = "foo")
    @test conf_twitch_err isa Javis.StreamConfig
    @test conf_twitch_err.livestreamto == :twitch
    @test isempty(conf_twitch_err.twitch_key)
    @test conf_twitch.twitch_key == "foo"

    render(vid, pathname = "stream_local.gif", streamconfig = conf_local)

    # errors with macos; a good test to have
    # test_local = run(pipeline(`lsof -i -P -n`, `grep ffmpeg`))
    # @test test_local isa Base.ProcessChain
    # @test test_local.processes isa Vector{Base.Process}

    cancel_stream()
    @test_throws ProcessFailedException run(
        pipeline(
            `ps aux`,
            pipeline(`grep ffmpeg`, pipeline(`grep stream_loop`, `awk '{print $2}'`)),
        ),
    )

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    @test_throws ErrorException render(
        vid,
        pathname = "stream_twitch.gif",
        streamconfig = conf_twitch_err,
    )
    rm("stream_twitch.gif")
end


text:
module ISE

export ise

"""
    ```julia
    ise(error)
    ```
    Integral Squared Error (ISE) fromm error vector.
    Inputs:
    error -> error vector.
"""
function ise(error::Vector{T}) where T
    sum((error.^2))
end

"""
    ```julia
    ise(error)
    ```
    Integral Squared Error (ISE) from vector of error vector.
    Inputs:
    error -> vector of error vector.
"""
function ise(error::Vector{Vector{T}}) where T
    map(x->sum((x.^2)), error)
end

"""
    ```julia
    ise(x, t, desired)
    ```
    Integral Squared Error (ISE) from error betwwen desired value and controler output.
    Inputs:
    x -> controler output vector;
    t -> time vector;
    desired -> desired value function.
"""
function ise(x::Vector{T}, t::Vector{T}, desired::Function) where T
    x = x .- map(desired, t)
    sum((x.^2))
end

end