text:
## Mesh convergence tests
using AeroMDAO
using Plots
using DataFrames

## Wing
wing = Wing(foils     = fill(naca4((0,0,1,2)), 2),
            chords    = [1.0, 0.6],
            twists    = [0.0, 0.0],
            spans     = [5.0],
            dihedrals = [5.],
            sweeps      = [5.]);
print_info(wing, "Wing")
S, b, c = projected_area(wing), span(wing), mean_aerodynamic_chord(wing);

## Assembly
ρ       = 1.225
ref     = [0.25, 0., 0.]
V, α, β = 1.0, 5.0, 5.0
Ω       = [0.0, 0.0, 0.0]
fs      = Freestream(V, α, β, Ω)

## Panelling error closure
panel_case(x, y) = solve_case(wing, fs;
                              rho_ref   = ρ,
                              r_ref     = ref,
                              area_ref  = S,
                              span_ref  = b,
                              chord_ref = c,
                              span_num  = y,
                              chord_num = x,
                              viscous   = false,
                              x_tr      = [0.3, 0.3],
                              spacing   = Uniform()
                             );

## Distributions
xs = 10:2:50
ys = xs .* 2# fill(maximum(xs) * 4, length(xs))

num = length(xs)

coeff_data = zeros(num, 9 + 2)
for (i, (x, y)) in enumerate(zip(xs, ys))
    @time nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, camber_panels, horseshoes, Γs =
    panel_case(x, y)

    println("Distribution: nx = $x, ny = $y")
    print_coefficients(nf_coeffs, ff_coeffs)

    coeff_data[i,1]     = x
    coeff_data[i,2]     = y
    coeff_data[i,3:5]   = ff_coeffs
    coeff_data[i,6:end] = nf_coeffs
end

##
data = DataFrame(coeff_data, :auto)
rename!(data, [:chord_num, :span_num, :CD_ff, :CY_ff, :CL_ff, :CD_nf, :CY_nf, :CL_nf, :Cl, :Cm, :Cn])

##
errors = @. abs((data[3:end,3:end] - data[2:end-1,3:end]) / data[2:end-1,3:end])

##
plot(xlabel = "Chordwise Panels", ylabel = "Error", yscale = :log10)

# CDs Errors
plot!(data[2:end-1,"chord_num"], errors[!,"CD_nf"], label = "CD Nearfield")
plot!(data[2:end-1,"chord_num"], errors[!,"CD_ff"], label = "CD Farfield")

# CLs Errors
plot!(data[2:end-1,"chord_num"], errors[!,"CL_nf"], label = "CL Nearfield")
plot!(data[2:end-1,"chord_num"], errors[!,"CL_ff"], label = "CL Farfield")

# CYs Errors
plot!(data[2:end-1,"chord_num"], errors[!,"CY_nf"], label = "CY Nearfield")
plot!(data[2:end-1,"chord_num"], errors[!,"CY_ff"], label = "CY Farfield")

##
CL = plot(data[2:end,"chord_num"], data[2:end,"CL_nf"], marker = :dot, label = "Nearfield", ylabel = "CL")
plot!(data[2:end,"chord_num"], data[2:end,"CL_ff"], marker = :dot, label = "Farfield")

CD = plot(data[2:end,"chord_num"], data[2:end,"CD_nf"], marker = :dot, label = "Nearfield", ylabel = "CD")
plot!(data[2:end,"chord_num"], data[2:end,"CD_ff"], marker = :dot, label = "Farfield")

CY = plot(data[2:end,"chord_num"], data[2:end,"CY_nf"], marker = :dot, label = "Nearfield", ylabel = "CY")
plot!(data[2:end,"chord_num"], data[2:end,"CY_ff"], marker = :dot, label = "Farfield")

Cl = plot(data[2:end,"chord_num"], data[2:end,"Cl"], marker = :dot, label = "Nearfield", ylabel = "Cl")
Cm = plot(data[2:end,"chord_num"], data[2:end,"Cm"], marker = :dot, label = "Nearfield", ylabel = "Cm")
Cn = plot(data[2:end,"chord_num"], data[2:end,"Cn"], marker = :dot, label = "Nearfield", ylabel = "Cn")

plot(CL, CD, CY, Cl, Cm, Cn, layout = (2,3), size = (1280, 720))


text:
struct FeatureSubset{T, R} <: FeatureCreator
    fc::T
    subset::R
end

Base.size(fc::FeatureSubset) = Base.size(fc.fc)
Base.length(fc::FeatureSubset) = Base.size(fc.fc)

(fc::FeatureSubset)(s; kwargs...) =
    fc.fc(s[fc.subset]; kwargs...)

(fc::FeatureSubset)(s,a,s_tp1) = fc.fc(s[fc.subset],a,s_tp1[fc.subset])

struct Sparsify{FC}
    fc::FC
end

function (fc::Sparsify)(s; kwargs...)
    ints = fc.fc(s; kwargs...)
    s = spzeros(size(fc))
    if ints isa AbstractArray
        s[ints] .= 1
    else
        s[ints] = 1
    end
    return s
end

function get_active_action_state_vector(state::SparseVector, action, feature_size, num_actions)
    vec_length = feature_size * num_actions
    new_ind = (state.nzind .- 1) * num_actions .+ action
    active_state_action = sparsevec(new_ind, state.nzval, vec_length)
    return active_state_action
end

mutable struct ValueFeatureProjector <: FeatureCreator
    func::Function
    pf_length::Int
end

function project_features(FP::ValueFeatureProjector, state)
    return FP.func(state, action)
end
(FP::ValueFeatureProjector)(state,action) = project_features(FP, state, action)

Base.size(VFP::ValueFeatureProjector) = VFP.pf_length

mutable struct ActionValueFeatureProjector{T} <: FeatureCreator
    func::T
    num_actions::Int
end

function project_features(FP::ActionValueFeatureProjector,state,action,state_tp1)
    state_action = spzeros(size(FP.func) * FP.num_actions)
    state = FP.func(state,action,state_tp1)
    state_length = size(FP.func)
    state_action[(action-1)*state_length + 1: action*state_length] = state
    return state_action
end

(FP::ActionValueFeatureProjector)(state,action,next_state) = project_features(FP, state, action, next_state)

Base.length(AVFP::ActionValueFeatureProjector) = size(AVFP.func) * AVFP.num_actions
Base.size(AVFP::ActionValueFeatureProjector) = size(AVFP.func) * AVFP.num_actions

struct FeatureProjector{FC} <: FeatureCreator
    fc::FC
    next_obs::Bool
end

Base.size(fc::FeatureProjector) = Base.size(fc.fc)
Base.length(fc::FeatureProjector) = Base.size(fc.fc)

(fc::FeatureProjector)(s, a, s_prime; kwargs...) = if fc.next_obs
    fc.fc(s_prime; kwargs...)
else
    fc.fc(s; kwargs...)
end

ext: jl

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# Usage:
# julia scripts/generate.jl /path/to/dir/containing/ros/repos

include("../src/gen_msgs.jl")

if !isdir("src")
    error("Please run this in the RobotOSData root directory")
end

mkpath("src/ros")
src_std = "$(ARGS[1])/std_msgs"
common_root = "$(ARGS[1])/common_msgs"

if !isdir(joinpath(src_std, "msg"))
    error("Wrong root directory: message files not found")
end

mkpath("src/ros/std_msgs/")
gen_module(:StdMsgs, [src_std], "src/ros/std_msgs/")

# Order these according to mutual dependencies:
src_common = [
    "$common_root/geometry_msgs",
    "$common_root/actionlib_msgs",
    "$common_root/diagnostic_msgs",
    "$common_root/nav_msgs",
    "$common_root/sensor_msgs",
    "$common_root/shape_msgs",
    "$common_root/stereo_msgs",
    "$common_root/trajectory_msgs",
    "$common_root/visualization_msgs",
]

mkpath("src/ros/common_msgs/")
gen_module(:CommonMsgs, src_common, "src/ros/common_msgs/", :(RobotOSData.StdMsgs), :(RobotOSData.CommonMsgs))


text:
"""
    ZerosBridge{T} <: Bridges.Variable.AbstractBridge

Transforms constrained variables in [`MathOptInterface.Zeros`](@ref) to zeros,
which ends up creating no variables in the underlying model.
The bridged variables are therefore similar to parameters with zero values.
Parameters with non-zero value can be created with constrained variables in
[`MOI.EqualTo`](@ref) by combining a [`VectorizeBridge`](@ref) and this bridge.
The functions cannot be unbridged, given a function, we cannot determine, if
the bridged variables were used.
The dual values cannot be determined by the bridge but they can be determined
by the bridged optimizer using [`MathOptInterface.Utilities.get_fallback`](@ref)
if a `CachingOptimizer` is used (since `ConstraintFunction` cannot be got
as functions cannot be unbridged).
"""
struct ZerosBridge{T} <: AbstractBridge
    n::Int # Number of variables
end
function bridge_constrained_variable(
    ::Type{ZerosBridge{T}},
    model::MOI.ModelLike,
    set::MOI.Zeros,
) where {T}
    return ZerosBridge{T}(MOI.dimension(set))
end

function supports_constrained_variable(::Type{<:ZerosBridge}, ::Type{MOI.Zeros})
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:ZerosBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:ZerosBridge})
    return Tuple{DataType,DataType}[]
end

# Attributes, Bridge acting as a model
MOI.get(bridge::ZerosBridge, ::MOI.NumberOfVariables) = 0
function MOI.get(bridge::ZerosBridge, ::MOI.ListOfVariableIndices)
    return MOI.VariableIndex[]
end

# References
function MOI.delete(::MOI.ModelLike, ::ZerosBridge) end

# Attributes, Bridge acting as a constraint

function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::ZerosBridge)
    return MOI.Zeros(bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.ConstraintPrimal,
    bridge::ZerosBridge{T},
) where {T}
    return zeros(T, bridge.n)
end

function MOI.get(
    ::MOI.ModelLike,
    ::MOI.VariablePrimal,
    ::ZerosBridge{T},
    ::IndexInVector,
) where {T}
    return zero(T)
end

function MOIB.bridged_function(::ZerosBridge{T}, ::IndexInVector) where {T}
    return zero(MOI.ScalarAffineFunction{T})
end
function unbridged_map(::ZerosBridge, ::MOI.VariableIndex, ::IndexInVector)
    return nothing
end


text:
using JuMP
using Stat31020
using Base.Test

x0 = ones(2); x1 = x0[1]; x2=x0[2]

# testing:
## -----------------------------------------------------------------------------
## min:  f = sum(x[1]^4 - 2.0*x[1]^2*x[2] + x[1]^2 + x[2]^2 - 2*x[1] for i=I)
## s.t.: -(x[1] + 0.25)^2 + 0.75*x[2] >= 0)
##       -(x[1] + 0.25) + 0.75*x[2]^2 >= 0)
## -----------------------------------------------------------------------------
f, g, h, j, kkt = Stat31020.constrained_wrap(x0, 2; kkt=true, μ=1.0, λ=ones(2))

function h_analytic(x1, x2)
    h = [[(12.0*x1^2 - 4.0*x2 + 2.) (-4.0*x1)];[(-4.0*x1) (2.)]] +
        [[-2. 0.];[0. 0.]] + [[0. 0.];[0. 1.5]]
    return h
end
function j_analytic(x1, x2)
    j = [[(-2.0*(x1 + 0.25)) 0.75];[-1.0 (1.5*x2)]]
    return j
end

h_a = h_analytic(x1, x2)
j_a = j_analytic(x1, x2)
z = zeros(2, 2)
kkt_a = [h_a j_a'; j_a z]

display(full(kkt))
println();
display(full(kkt_a))
println();
@test full(kkt) == full(kkt_a)


text:
using Test
using DelimitedFiles, LinearAlgebra
import MLJBase
const Mlj = MLJBase
using StableRNGs
#rng = StableRNG(123)
using BetaML.Trees

println("*** Additional testing for the Testing Decision trees/Random Forest algorithms...")

println("Testing MLJ interface for Trees models....")
X, y                           = Mlj.@load_boston
model_dtr                      = DecisionTreeRegressor()
regressor_dtr                  = Mlj.machine(model_dtr, X, y)
Mlj.evaluate!(regressor_dtr, resampling=Mlj.CV(), measure=Mlj.rms, verbosity=0)

model_rfr                      = RandomForestRegressor()
regressor_rfr                  = Mlj.machine(model_rfr, X, y)
Mlj.evaluate!(regressor_rfr, resampling=Mlj.CV(), measure=Mlj.rms, verbosity=0)

X, y                           = Mlj.@load_iris
model_dtc                      = DecisionTreeClassifier()
regressor_dtc                  = Mlj.machine(model_dtc, X, y)
Mlj.evaluate!(regressor_dtc, resampling=Mlj.CV(), measure=Mlj.LogLoss())

model_rfc                      = RandomForestClassifier(maxFeatures=3)
regressor_rfc                  = Mlj.machine(model_rfc, X, y)
Mlj.evaluate!(regressor_rfc, resampling=Mlj.CV(), measure=Mlj.LogLoss())

#=
using MLJ
X, y                           = Mlj.@load_boston
MLJ.models(MLJ.matching(X,y))
Tree = @load DecisionTreeRegressor pkg=BetaML
tree = Tree()
=#


text:
using Test
using Flux
using FeedbackNets

@testset "FeedbackNets" begin

include("splitter_tests.jl")
include("merger_tests.jl")
include("feedbackchain_tests.jl")
include("feedbacktree_tests.jl")

end # @testset "FeedbackNets"


text:
#=Currently to test Assoc please use include()
=#
include("../src/D4M.jl");

A = Assoc("a|b|c|d|","v,x,y,z,",[1,2,3,4]);
a = A["a,b,","x,v,"];
B = Assoc("v,x,y,z,","a+b+c+d+",[5,6,7,8]);
#currently Dump works

#=
ToDos
writecsv and readcsv functions - utilize Julia csv functions.
=#
println("Constructing histogram of A on column x and z..");
C = transpose( nocol(A[:,"v,"]) & nocol(A[:,"z,"])) *A;
print(C);


text:
__precompile__()

module Quaternions

  import Base: +, -, *, /, ^, ==
  import Base: abs, abs2, angle, conj, cos, exp, inv, isreal, isfinite, isinf, iszero, isnan, log, real, sin, sqrt
  import Base: convert, promote_rule, float
  import Base: rand, randn
  import LinearAlgebra: lyap, norm, normalize, sylvester
  using LinearAlgebra: cross, dot
  using Random

  Base.@irrational INV_SQRT_EIGHT 0.3535533905932737622004 sqrt(big(0.125))

  include("Quaternion.jl")
  include("Octonion.jl")
  include("DualQuaternion.jl")

  export Quaternion,
    QuaternionF16,
    QuaternionF32,
    QuaternionF64,
    Octonion,
    OctonionF16,
    OctonionF32,
    OctonionF64,
    DualQuaternion,
    DualQuaternionF16,
    DualQuaternionF32,
    DualQuaternionF64
  export quat
  export octo
  export imag_part
  export dualquat
  export angleaxis
  export angle
  export axis
  export linpol
  export normalize
  export normalizea
  export dconj
  export quatrand
  export nquatrand
  export octorand
  export dualquatrand
  export ndualquatrand
  export qrotation
  export rotationmatrix
  export slerp
end


text:
using RvSpectML
using DataFrames, CSV

function ccf_total(order_list_timeseries::AbstractChunkListTimeseries, line_list_df::DataFrame, pipeline::PipelinePlan; recalc::Bool = false,
                  output_fn_suffix::String = "", range_no_mask_change::Real=RvSpectMLBase.max_bc, ccf_mid_velocity::Real=0.0,
                  mask_scale_factor::Real=1, mask_type::Symbol = :tophat, use_expr::Bool = false )
    if need_to(pipeline,:ccf_total) || recalc
      if verbose println("# Computing CCF.")  end
      @assert !need_to(pipeline,:extract_orders)
      @assert !need_to(pipeline,:clean_line_list_tellurics)
      if mask_type == :tophat
        mask_shape = CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor*mask_scale_factor)
      else
        @error("Requested mask shape (" * string(mask_type) * " not avaliable.")
      end

      line_list = CCF.BasicLineList(line_list_df.lambda, line_list_df.weight)
      ccf_plan = CCF.BasicCCFPlan(mask_shape = mask_shape, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=range_no_mask_change)
      v_grid = CCF.calc_ccf_v_grid(ccf_plan)
      if use_expr
        @time ccfs = CCF.calc_ccf_chunklist_timeseries_expr(order_list_timeseries, ccf_plan)
      else
        @time ccfs = CCF.calc_ccf_chunklist_timeseries(order_list_timeseries, ccf_plan)
      end
      #mask_shape_expr = CCF.TopHatCCFMask(order_list_timeseries.inst, scale_factor=tophap_ccf_mask_scale_factor*11.0)
      #mask_shape_expr = CCF.GaussianCCFMask(order_list_timeseries.inst, scale_factor=9)
      # Warning:  CCF with a shape other than a tophat is still experimental
      #ccf_plan_expr = CCF.BasicCCFPlan(mask_shape = mask_shape_expr, line_list=line_list, midpoint=ccf_mid_velocity, range_no_mask_change=22e3)
      #@time ccfs_expr = CCF.calc_ccf_chunklist_timeseries_expr(order_list_timeseries, ccf_plan_expr) #, verbose=true)
      #println("# Ratio of max(ccfs_expr)/max(ccfs) = ", mean(maximum(ccfs_expr,dims=1)./maximum(ccfs,dims=1)) )
      #=
      mask_shape_expr2 = RvSpectML.CCF.CosCCFMask(order_list_timeseries.inst, scale_factor=18)   # Why does such a large value help so much?
      ccf_plan_expr2 = RvSpectML.CCF.BasicCCFPlan(mask_shape = mask_shape_expr2, line_list=line_list, midpoint=ccf_mid_velocity)
      @time ccfs_expr2 = RvSpectML.CCF.calc_ccf_chunklist_timeseries_expr(order_list_timeseries, ccf_plan_expr2) #, verbose=true)
      =#
      if save_data(pipeline, :ccf_total)
         CSV.write(joinpath(output_dir,target_subdir * "_ccfs" * output_fn_suffix * ".csv"),Tables.table(ccfs',header=Symbol.(v_grid)))
         #CSV.write(joinpath(output_dir,target_subdir * "_ccfs_expr.csv"),Tables.table(ccfs_expr',header=Symbol.(v_grid)))
      end
      set_cache!(pipeline, :ccf_total, (ccfs=ccfs, v_grid=v_grid) )
      dont_need_to!(pipeline,:ccf_total)
    end

    if has_cache(pipeline,:ccf_total) return read_cache(pipeline,:ccf_total)
    else   @error("Invalid pipeline state.")          end
end


text:
using Compose
using Test

Compose.main()


text:
function simulate!(s::Simulator{T};
	clock_time_noise=1.0e-3,
	verbose=false) where T
    status = false

    p = s.policy
	w = s.dist
	traj = s.traj

	N = length(traj.u)

	# state
	t = 1
	x = [traj.q[t+1]; (traj.q[t+1] - traj.q[t]) ./ s.h; traj.γ[t]]

	# last_computed_control = p.ref_traj.u[1] / N_sample
	clock_time = 0.0
    for t = 1:N
		# println("timestep $t / $N, time $time")
		policy_time = @elapsed traj.u[t] = exec_policy(p, x, clock_time)
		traj.u[t] .*= s.h # transform force from policy to impulse for RoboDojo simulator
		s.opts.record && (s.stats.policy_time[t] = policy_time)

        # disturbances
        traj.w[t] .= disturbances(w, traj.q[t+1], t)

        # step
        status = RoboDojo.step!(s, t, verbose=verbose)
		x .= [traj.q[t+2]; (traj.q[t+2] - traj.q[t+1]) ./ s.h; traj.γ[t]]
		clock_time += s.h + clock_time_noise * rand(1)[1]

        !status && break
    end

    return status
end

function exec_policy(p::CIMPC{T,NQ,NU,NW,NC}, x::Vector{T}, t::T) where {T,NQ,NU,NW,NC}
	# reset
	if t ≈ 0.0
		# p.next_time_update = 0.0
		p.q0 .= p.ref_traj.q[1]
		# p.altitude .= 0.0
		p.buffer_time = 0.0
		# set_trajectory!(p.traj, p.ref_traj)
		set_implicit_trajectory!(p.im_traj, p.im_traj_cache)
		reset_window!(p.window)
	end

	policy_time = @elapsed begin
		if p.buffer_time <= 0.0

			# window
			idx_nearest = findnearest(p.times_reference, t % (p.ref_traj.h * (p.ref_traj.H - 1)))[1]
			for i = 1:(p.H + 2)
				p.window[i] = i + idx_nearest - 1 > p.ref_traj.H ? i - p.ref_traj.H + idx_nearest - 1 : i + idx_nearest - 1
			end
			set_window!(p.traj, p.ref_traj, p.window)

			# optimize
			q1 = x[1:NQ]
			q0 = x[1:NQ] - x[NQ .+ (1:NQ)] .* p.traj.h
			newton_solve!(p.newton, p.s, q0, q1,
				p.window, p.im_traj, p.traj, warm_start = t > 0.0)
		end
	end

	# update buffer time
	(p.buffer_time <= 0.0) && (p.buffer_time = policy_time)
	p.buffer_time -= p.traj.h

	# control
	p.u .= p.newton.traj.u[1]

	# add gains
	if p.opts.gains
		K1 = p.K_traj[p.window[1]]

		q1 = x[1:NQ]
		q0 = x[1:NQ] - x[NQ .+ (1:NQ)] .* p.traj.h
		p.u .+= K1 * ([p.traj.q[p.window[1]]; p.traj.q[p.window[2]]] - [q0; q1])
	end

	p.u ./= p.traj.h

	return p.u
end

# horizon_reference = 11
# horizon_mpc = 3
# timestep = 0.1
# time_reference = timestep * (horizon_reference - 1)
# times_reference = range(0, stop=time_reference, length=horizon_reference)

# t = time_reference - 0.1

function findnearest(a,x)
	length(a) > 0 || return 0:-1
	r = searchsorted(a,x)
	length(r) > 0 && return r
	last(r) < 1 && return searchsorted(a,a[first(r)])
	first(r) > length(a) && return searchsorted(a,a[last(r)])
	x-a[last(r)] < a[first(r)]-x && return searchsorted(a,a[last(r)])
	x-a[last(r)] > a[first(r)]-x && return searchsorted(a,a[first(r)])
	return first(searchsorted(a,a[last(r)])):last(searchsorted(a,a[first(r)]))
end

# idx_nearest = findnearest(times_reference, t % time_reference)[1]
# idx_window_wrap = [i + idx_nearest - 1 > horizon_reference ? i - horizon_reference + idx_nearest - 1 : i + idx_nearest - 1 for i in 1:(horizon_mpc + 2)]

# ###
# # window
# t = 0.2
# p.window .= 0
# idx_nearest = findnearest(p.times_reference, t % (p.ref_traj.h * (p.ref_traj.H - 1)))[1]
# for i = 1:(p.H + 2)
# 	p.window[i] = i + idx_nearest - 1 > p.ref_traj.H ? i - p.ref_traj.H + idx_nearest - 1 : i + idx_nearest - 1
# end
# p.window
# @show idx_nearest
# @show length(p.traj.q)
# @show length(p.traj.u)
# @show length(p.ref_traj.q)
# @show length(p.ref_traj.u)
# @show length(p.window)
# @show p.window
# set_window!(p.traj, p.ref_traj, p.window)


text:
#Functions recognized by Mathenize
sym = [ :sqrt, :+, :-, :/, :^, :tan, :*,
        :sin, :cos, :sincos, :sind, :cosd,
        :tand, :sinh, :cosh, :tanh, :log,
        :pi, :π, :\, :fma, :muladd, :inv,
        :div, :÷, :fld, :cld, :mod, :eps,
        :rem2pi, :mod2pi, :divrem, :fldmod,
        :fld1, :mod1, ://, :rationalize,
        :numerator, :denominator, :<<, :>>,
        :>>>, :bitrotate, :cmp, :complex,
        :~, :&, :|, :asin, :acos, :atan,
        :asind, :acosd, :atand, :sec, :csc,
        :cot, :secd, :cscd, :cotd, :asec,
        :acsc, :acot, :asecd, :acscd, :acotd,
        :sech, :csch, :coth, :asinh, :acosh,
        :atanh, :atan, :asech, :acsch, :acoth,
        :sinc, :cosc, :deg2rad, :hypot, :log2,
        :log10, :log1p, :frexp, :exp, :exp2,
        :exp10, :ldexp, :modf, :expm1, :round,
        :ceil, :floor, :trunc, :unsafe_trunc,
        :min, :max, :minmax, :clamp, :clamp!,
        :abs, :abs2, :copysign, :sign, :signbit,
        :flipsign, :isqrt, :cbrt, :real, :imag,
        :reim, :conj, :angle, :cis, :cispi,
        :binomial, :factorial, :gcd, :lcm,
        :gcdx, :ispow2, :nextpow, :prevpow,
        :nextprod, :invmod, :powermod, :ndigits,
        :widemul, :evalpoly, :@evalpoly, :im,
        :vcat, :hcat, :sum]

#Iterate through expression to verify all items in it are valid symbols and operations.
function subtasking(math, tasks, sbtask, LOG_INFO, print_info)

    success = false
    for ñ in sbtask[1:tasks]

        if hastask(ñ)

            for i in 1:length(ñ.args)

                if ñ.args[i] in sym
                    success = true
                    continue
                elseif ñ.args[i] isa Number
                    success = true
                    continue
                elseif ispermitted(ñ.args[i], LOG_INFO)
                    success = true
                    continue
                else
                    push!(LOG_INFO, "       └ $(ñ) -> $(typeof(ñ)) ||| $(ñ.args[1]) -> $(typeof(ñ)) was not a recognized subtask")
                    unknownmath(ñ, LOG_INFO, print_info)
                end
            end

        elseif ispermitted(ñ, LOG_INFO)
            success = true
        else
            push!(LOG_INFO, "       └ $(ñ) was not a recognized task")
            unknownmath(ñ, LOG_INFO, print_info)
        end
    end

    push!(LOG_INFO, "    Subtasks were all read\n(...) Checking if expression was parsed and can be computed")
    if success
        push!(LOG_INFO, "! Expression was parsed and computed succesfully\n")
        return Core.eval(Base.Math, math)
    else
        push!(LOG_INFO, "! Expression was not parsed succesfully\n")
        return nothing
    end

end

#Check if value is a valid math operation, such as a mathematical function, number, vector, or matrix.
function ispermitted(tsk, LOG_INFO)

    #Adding log info.
    push!(LOG_INFO, "    └ -> $(tsk) is being checked, its type is: $(typeof(tsk))")

    if tsk in sym

        #tsk is in Array named sym.
        push!(LOG_INFO, "        └ $(tsk) is permitted that belongs to: $(tsk)")
        return true

    elseif tsk isa Number

        push!(LOG_INFO, "        └ $(tsk) is permitted that belongs to: $(tsk)")
        return true

    elseif tsk isa Expr && hasproperty(tsk, :head) && tsk.head == :row

        #NOT SAFE, NOT PROVEN TO BE EVAL SAFE, DO NOT USE WITHOUT PRECAUTION.
        push!(LOG_INFO, "        └ $(tsk) is a valid row that belongs to: $(tsk)")
        return true

    elseif tsk isa Expr && hasproperty(tsk, :head) && tsk.head == :hcat || tsk.head == :vcat || tsk.head == :vect

        push!(LOG_INFO, "        └ $(tsk) is a valid matrix or vector that belongs to: $(tsk)")
        return true

    elseif tsk isa Expr && tsk.args[1] in sym

        #Checking if expression's first argument is in sym.
        push!(LOG_INFO,"        └ $(tsk.args[1]) is a valid expression found in sym that belongs to: $(tsk)")
        return true

    else

        #Not found.
        push!(LOG_INFO, "       └ $(tsk) was not recognized as permitted")
        return false

    end
end

#Check if argument contains subarguments.
function hastask(sb)
    return hasproperty(sb, :args) ? true : false
end

#Error when given input contains an unknown operation.
function unknownmath(ñ, LOG_INFO, print_info::Bool)

    items = "and is an empty value."
    if !print_info
        message_info = "Check the log using calculate(math::String, true)"
    else
        message_info = ""
        @info join(LOG_INFO, "\n")
    end
    if hastask(ñ) items = "that contains $(ñ.args)" end
    error("Error in Mathenize syntax. $(message_info)\n└ ->$(ñ) is not recognized as a valid math operation. \n └ The input given is a $(typeof(ñ)) $(items)")
    return nothing

end


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(256, board_shim)
brainflow.release_session(board_shim)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
using Stage
using Base.Test

@expect 1 == 1
try @expect 1 == 2
catch e
  @test e.msg == "test failed 1 == 2"
end

@expect log(e) == 1.0
try @expect log(e) == 2.0
catch e
  @test e.msg == "test failed log(e) == 2.0"
end

@expect log(e) < 2.0
try @expect log(e) >= 2.0
catch e
  @test e.msg == "test failed log(e) >= 2.0"
end

@expect abs(1.0 - 0.9) < 0.1


text:
using Pkg

Pkg.update()

Pkg.add("Plots")
Pkg.add("DataFrames")
Pkg.add("CSV")

Pkg.build("Plots")
Pkg.build("DataFrames")
Pkg.build("CSV")

max_stars_repo_name: yaozhenghangma/Julia-notebook

max_issues_repo_name: yaozhenghangma/Julia-notebook

max_forks_repo_name: yaozhenghangma/Julia-notebook


text:
using SpecialFunctions

"""
    Humlicek (1982) rational approximation:  region I.
"""
function humlicek1(z::Complex)
    FT = eltype(real(z));
    #t = z.imag-1j*z.real;   w = t * recSqrtPi / (0.5 + t*t)
    w = 1im*FT(1/sqrt(pi)) *z / (z*z-FT(0.5))
end

"""
    Humlicek (1982) rational approximation:  region II.
"""
function humlicek2(z::Complex)
    FT = eltype(real(z));
    # this implementation corresponds to the fortran code
    t = imag(z)-1im*real(z)
    u = t * t
    w = (t * (FT(1.410474) + u*FT(1/sqrt(pi))))/ (3/4 + (u *(3+u)))
end

"""
    Humlicek (1982) rational approximation:  region II.
"""
function Humlicek2(z::Complex)
    FT = eltype(real(z));
    # this implementation corresponds to Eq. (12) of the manuscript
    zz = z*z
    w = 1im* (z * (zz*FT(1/sqrt(pi))-FT(1.410474)))/ (3/4 + zz*(zz-3))
end

"""
    Humlicek (1982) rational approximation:  region III.
"""
function humlicek3(z::Complex)
    FT = eltype(real(z));
    t = imag(z)-1im*real(z)
    w =  (FT(16.4955) + t * (FT(20.20933) + t * (FT(11.96482) + t * (FT(3.778987) + FT(0.5642236)*t)))) /
         (FT(16.4955) + t * (FT(38.82363) + t * (FT(39.27121) + t * (FT(21.69274) + t * (FT(6.699398) + t)))))
end

"""
    Humlicek (1982) rational approximation:  region IV.
"""
function humlicek4(z::Complex)
    FT = eltype(real(z));
    t = imag(z)-1im*real(z)
    u = t*t
    nom = t*(FT(36183.31)-u*(FT(3321.99)-u*(FT(1540.787)-u*(FT(219.031)-u*(FT(35.7668)-u*(FT(1.320522)-u*FT(.56419)))))))
    den = FT(32066.6)-u*(FT(24322.8)-u*(FT(9022.23)-u*(FT(2186.18)-u*(FT(364.219)-u*(FT(61.5704)-u*(FT(1.84144)-u))))))
    w  = exp(u) - nom/den
end

"""
    Humlicek (1982) complex probability function:  w4 rational approximations.
"""
function humlicek(z::Complex)
    s = abs(real(z)) + imag(z)
    # Choices depending on regions:
    if s>15.0
        w = humlicek1(z)
    elseif s>5.5
        w = humlicek2(z)
    else
        if imag(z)>=0.195*abs(real(z))-0.176
            w = humlicek3(z)
        else
            w = humlicek4(z)
        end
    end
    return w
end

"""
    Humlicek (1979) complex probability function:  rational approximation for y>0.85 OR |x|<18.1*y+1.65.
"""
function cpf12a(z::Complex)
    FT = eltype(real(z));
    cr  = FT(1.5);
    crr = FT(2.25);
    #a = similar(z)
    ct1 = FT(.3142403762544);  ct2 = FT(.9477883912402);  ct3 = FT(1.5976826351526);
    ct4 = FT(2.2795070805011); ct5 = FT(3.0206370251209); ct6 = FT(3.88972489786978);

    ca1 = FT(-1.393236997981977);  ca2 = FT(-0.2311524061886763);  ca3 = FT(+0.1553514656420944);
    ca4 = FT(-0.006218366236965554); ca5 = FT(9.190829861057117e-5); ca6 = FT(+6.275259577e-7);

    cb1 = FT(1.011728045548831);  cb2 = FT(-0.7519714696746353);  cb3 = FT(0.01255772699323164);
    cb4 = FT(0.01002200814515897); cb5 = FT(-2.420681348155727e-4); cb6 = FT(5.008480613664576e-7);

    x = real(z)
    y = imag(z)
    ry = cr + y
    ryy = ry * ry
    wk =  (ca1*(x-ct1) + cb1*ry) / ((x-ct1)^2 + ryy) - (ca1*(x+ct1) - cb1*ry) / ((x+ct1)^2 + ryy) +
          (ca2*(x-ct2) + cb2*ry) / ((x-ct2)^2 + ryy) - (ca2*(x+ct2) - cb2*ry) / ((x+ct2)^2 + ryy) +
          (ca3*(x-ct3) + cb3*ry) / ((x-ct3)^2 + ryy) - (ca3*(x+ct3) - cb3*ry) / ((x+ct3)^2 + ryy) +
          (ca4*(x-ct4) + cb4*ry) / ((x-ct4)^2 + ryy) - (ca4*(x+ct4) - cb4*ry) / ((x+ct4)^2 + ryy) +
          (ca5*(x-ct5) + cb5*ry) / ((x-ct5)^2 + ryy) - (ca5*(x+ct5) - cb5*ry) / ((x+ct5)^2 + ryy) +
          (ca6*(x-ct6) + cb6*ry) / ((x-ct6)^2 + ryy) - (ca6*(x+ct6) - cb6*ry) / ((x+ct6)^2 + ryy);
    wl =  (cb1*(x-ct1) - ca1*ry) / ((x-ct1)^2 + ryy) + (cb1*(x+ct1) + ca1*ry) / ((x+ct1)^2 + ryy) +
          (cb2*(x-ct2) - ca2*ry) / ((x-ct2)^2 + ryy) + (cb2*(x+ct2) + ca2*ry) / ((x+ct2)^2 + ryy) +
          (cb3*(x-ct3) - ca3*ry) / ((x-ct3)^2 + ryy) + (cb3*(x+ct3) + ca3*ry) / ((x+ct3)^2 + ryy) +
          (cb4*(x-ct4) - ca4*ry) / ((x-ct4)^2 + ryy) + (cb4*(x+ct4) + ca4*ry) / ((x+ct4)^2 + ryy) +
          (cb5*(x-ct5) - ca5*ry) / ((x-ct5)^2 + ryy) + (cb5*(x+ct5) + ca5*ry) / ((x+ct5)^2 + ryy) +
          (cb6*(x-ct6) - ca6*ry) / ((x-ct6)^2 + ryy) + (cb6*(x+ct6) + ca6*ry) / ((x+ct6)^2 + ryy);
    a = wk + im * wl
    #println(a)
    return a#a   # wk, wl
end

"""
    Humlicek (1979) complex probability function:  rational approximation for y<0.85 AND |x|>18.1*y+1.65.
"""
function cpf12b(z::Complex)
    FT = eltype(real(z));
    ct1 = FT(.3142403762544);  ct2 = FT(.9477883912402);  ct3 = FT(1.5976826351526);
    ct4 = FT(2.2795070805011); ct5 = FT(3.0206370251209); ct6 = FT(3.88972489786978);

    ca1 = FT(-1.393236997981977);  ca2 = FT(-0.2311524061886763);  ca3 = FT(+0.1553514656420944);
    ca4 = FT(-0.006218366236965554); ca5 = FT(9.190829861057117e-5); ca6 = FT(+6.275259577e-7);

    cb1 = FT(1.011728045548831);  cb2 = FT(-0.7519714696746353);  cb3 = FT(0.01255772699323164);
    cb4 = FT(0.01002200814515897); cb5 = FT(-2.420681348155727e-4); cb6 = FT(5.008480613664576e-7);
    cr  = FT(1.5);
    crr = FT(2.25);

    x = real(z);
    y = imag(z);
    ry   = cr+y      # yy0   = y+1.5
    y2r  = y +2cr  # y2y0  = y+3.0
    rry  = cr*ry     # y0yy0 = 1.5*(y+1.5)
    ryry = ry^2     # yy0sq = (y+1.5)^2
    wk =  ( cb1*((x-ct1)^2-rry) - ca1*(x-ct1)*y2r ) / (((x-ct1)^2+crr)*((x-ct1)^2+ryry)) +
          ( cb1*((x+ct1)^2-rry) + ca1*(x+ct1)*y2r ) / (((x+ct1)^2+crr)*((x+ct1)^2+ryry)) +
          ( cb2*((x-ct2)^2-rry) - ca2*(x-ct2)*y2r ) / (((x-ct2)^2+crr)*((x-ct2)^2+ryry)) +
          ( cb2*((x+ct2)^2-rry) + ca2*(x+ct2)*y2r ) / (((x+ct2)^2+crr)*((x+ct2)^2+ryry)) +
          ( cb3*((x-ct3)^2-rry) - ca3*(x-ct3)*y2r ) / (((x-ct3)^2+crr)*((x-ct3)^2+ryry)) +
          ( cb3*((x+ct3)^2-rry) + ca3*(x+ct3)*y2r ) / (((x+ct3)^2+crr)*((x+ct3)^2+ryry)) +
          ( cb4*((x-ct4)^2-rry) - ca4*(x-ct4)*y2r ) / (((x-ct4)^2+crr)*((x-ct4)^2+ryry)) +
          ( cb4*((x+ct4)^2-rry) + ca4*(x+ct4)*y2r ) / (((x+ct4)^2+crr)*((x+ct4)^2+ryry)) +
          ( cb5*((x-ct5)^2-rry) - ca5*(x-ct5)*y2r ) / (((x-ct5)^2+crr)*((x-ct5)^2+ryry)) +
          ( cb5*((x+ct5)^2-rry) + ca5*(x+ct5)*y2r ) / (((x+ct5)^2+crr)*((x+ct5)^2+ryry)) +
          ( cb6*((x-ct6)^2-rry) - ca6*(x-ct6)*y2r ) / (((x-ct6)^2+crr)*((x-ct6)^2+ryry)) +
          ( cb6*((x+ct6)^2-rry) + ca6*(x+ct6)*y2r ) / (((x+ct6)^2+crr)*((x+ct6)^2+ryry))
    wl =  (cb1*(x-ct1) - ca1*ry) / ((x-ct1)^2 + ryry) + (cb1*(x+ct1) + ca1*ry) / ((x+ct1)^2 + ryry) +
          (cb2*(x-ct2) - ca2*ry) / ((x-ct2)^2 + ryry) + (cb2*(x+ct2) + ca2*ry) / ((x+ct2)^2 + ryry) +
          (cb3*(x-ct3) - ca3*ry) / ((x-ct3)^2 + ryry) + (cb3*(x+ct3) + ca3*ry) / ((x+ct3)^2 + ryry) +
          (cb4*(x-ct4) - ca4*ry) / ((x-ct4)^2 + ryry) + (cb4*(x+ct4) + ca4*ry) / ((x+ct4)^2 + ryry) +
          (cb5*(x-ct5) - ca5*ry) / ((x-ct5)^2 + ryry) + (cb5*(x+ct5) + ca5*ry) / ((x+ct5)^2 + ryry) +
          (cb6*(x-ct6) - ca6*ry) / ((x-ct6)^2 + ryry) + (cb6*(x+ct6) + ca6*ry) / ((x+ct6)^2 + ryry)
    a = exp(-x * x) + y * wk + im * wl   # wk, wl
    return a
end

"""
    Complex error function --- J.A.C. Weideman: SIAM J. Num. Anal. 31 1497-1518 (1994); equation (38.I) and table I.
"""
function weideman32a(z::Complex)
    FT = eltype(real(z));
    L32 = FT(sqrt(32/sqrt(2)));
    a1=  FT(2.5722534081245696e+00);  a2 = FT(2.2635372999002676e+00);  a3 = FT(1.8256696296324824e+00);  a4 = FT(1.3455441692345453e+00);
    a5=  FT(9.0192548936480144e-01);  a6 = FT(5.4601397206393498e-01);  a7 = FT(2.9544451071508926e-01);  a8 = FT(1.4060716226893769e-01);
    a9=  FT(5.7304403529837900e-02);  a10= FT(1.9006155784845689e-02);  a11= FT(4.5195411053501429e-03);  a12= FT(3.9259136070122748e-04);
    a13= FT(-2.4532980269928922e-04); a14= FT(-1.3075449254548613e-04); a15= FT(-2.1409619200870880e-05); a16= FT(6.8210319440412389e-06);
    a17= FT(4.4015317319048931e-06);  a18= FT(4.2558331390536872e-07);  a19= FT(-4.1840763666294341e-07); a20= FT(-1.4813078891201116e-07);
    a21= FT(2.2930439569075392e-08);  a22= FT(2.3797557105844622e-08);  a23= FT(8.1248960947953431e-10);  a24= FT(-3.2080150458594088e-09);
    a25= FT(-5.2310170266050247e-10); a26= FT(4.1537465934749353e-10);  a27= FT(1.1658312885903929e-10);  a28= FT(-5.5441820344468828e-11);
    a29= FT(-2.1542618451370239e-11); a30= FT(8.0314997274316680e-12);  a31= FT(3.7424975634801558e-12);  a32= FT(-1.3031797863050087e-12);

    iz  = 1im*real(z) - imag(z)
    lpiz    = L32 + iz  # wL - y + x*complex(0.,1.)
    lmiz    = L32 - iz  # wL + y - x*complex(0.,1.)
    recLmiZ = 1 / lmiz
    Z       = lpiz * recLmiZ
    polynom = a1+(a2+(a3+(a4+(a5+(a6+(a7+(a8+(a9+(a10+(a11+(a12+(a13+(a14+(a15+(a16+(a17+(a18+(a19+(a20+(a21+(a22+(a23+(a24+(a25+(a26+(a27+(a28+(a29+(a30+(a31+a32*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z)*Z
    w       = (FT(1/sqrt(pi))  +  2 * polynom * recLmiZ)  *  recLmiZ
    return w
end

"""
    Humlicek (1979) complex probability function w(z): a single rational approximation.
"""
function w(::CPF12ErrorFunction, z::Complex)
    FT = eltype(real(z));
    # Needs to be double checked, implemented slightly differently from Schreier
    if (abs(real(z))<FT(18.1)*imag(z)+FT(1.65)) || (imag(z)>FT(0.85))
        return cpf12a(z)
    else
        return cpf12b(z)
    end
end

"""
    Complex error function w(z) using Humlicek's and Weideman's rational approximation:
        |x|+y>15:  Humlicek (JQSRT, 1982) rational approximation for region I;
        else:      J.A.C. Weideman (SIAM-NA 1994); equation (38.I) and table I.
"""
function w(::HumlicekWeidemann32VoigtErrorFunction, z::Complex)

    FT = eltype(real(z));
    if abs(real(z))+imag(z)>15
        w = 1im*FT(1/sqrt(pi)) * z / (z*z-FT(0.5))  # Humlicek (1982) approx 1 for s>15
    else
        w = weideman32a(z)
    end
    return w
end

"""
 Complex error function w(z) using Humlicek's and Weideman's rational approximation:
        |x|+y>8:  Humlicek (JQSRT, 1982) rational approximation for region II;
        else:      J.A.C. Weideman (SIAM-NA 1994); equation (38.I) and table I.
"""
function w(::HumlicekWeidemann32SDErrorFunction, z::Complex)
    FT = eltype(real(z));
    if abs(real(z))+imag(z)>=8
        w = humlicek2(z)  # Humlicek (1982) approx 1 for s>8
    else
        w = weideman32a(z)
    end
    return w
end

"""
    Complex error function w(z) using Humlicek's and Weideman's rational approximation:
        |x|+y>15:  Humlicek (JQSRT, 1982) rational approximation for region I;
        else:      using erfc(-iz) from Special Functions
"""
function w(::ErfcHumliErrorFunctionVoigt, z::Complex)
    FT = eltype(real(z));
    if abs(real(z))+imag(z)>15
        w = 1im*FT(1/sqrt(pi)) * z / (z*z-FT(0.5))  # Humlicek (1982) approx 1 for s>15
    else
        w = erfcx(-1im*z)
    end
    return w
end

"""
 Complex error function w(z) using Humlicek's and Weideman's rational approximation:
        |x|+y>8:  Humlicek (JQSRT, 1982) rational approximation for region II;
        else:     using erfc(-iz) from Special Functions
"""
function w(::ErfcHumliErrorFunctionSD, z::Complex)
    FT = eltype(real(z));
    if abs(real(z))+imag(z)>=8
        w = humlicek2(z)  # Humlicek (1982) approx 1 for s>8
    else
        w = erfcx(-1im*z)
    end
    return w
end

"""
 Complex error function w(z) using erfc(-iz) from Special Functions
"""
function w(::ErfcErrorFunction, z::Complex)
    erfcx(-im*z)
end


text:
using PDESolver
using SimpleODEMod

mesh, sbp, eqn, opts = run_simpleode(ARGS[1])


text:
#
# Bayes Net
#
lualatex_available() = try success(`lualatex -v`) catch; false end
Base.showable(::MIME"image/svg+xml", bn::BayesNet) = true

plot(dag::DAG, nodelabel) = gplot(dag,
								  nodelabel=nodelabel,
								  layout=stressmajorize_layout,
								  nodefillc="lightgray",
								  edgestrokec="black",
								  EDGELINEWIDTH=0.3) # GraphPlot (default plotting)

# called at runtime (replaces plot with TikzGraphs, if loaded)
function __init__()
	@require TikzGraphs="b4f28e30-c73f-5eaf-a395-8a9db949a742" begin
		if lualatex_available()
			plot(dag::DAG, nodelabel) = TikzGraphs.plot(dag, nodelabel)
		end
	end
end

function plot(bn::BayesNet)
	if !isempty(names(bn))
		dag = bn.dag
		nodelabel = AbstractString[string(s) for s in names(bn)] # NOTE: sometimes the same var shows up twice
	else
		dag = DiGraph(1)
		nodelabel = ["Empty Graph"]
	end

	plot(dag, nodelabel)
end

function Base.show(f::IO, a::MIME"image/svg+xml", bn::BayesNet)
 	show(f, a, plot(bn))
end


text:
## Drag polar analysis case
using AeroMDAO
using DataFrames
using Plots

## Lifting surfaces

## Surfaces

# Wing
wing = Wing(foils     = fill(naca4(0,0,1,2), 2),
            chords    = [1.0, 0.6],
            twists    = [0.0, 0.0],
            spans     = [8.0],
            dihedrals = [0.],
            sweeps    = [5.]);

# Horizontal tail
htail = Wing(foils     = fill(naca4(0,0,1,2), 2),
             chords    = [0.7, 0.42],
             twists    = [0.0, 0.0],
             spans     = [2.5],
             dihedrals = [0.],
             sweeps    = [6.39],
             position  = [4., 0, 0.],
             angle     = 0.,
             axis      = [0., 1., 0.])

# Vertical tail
vtail = HalfWing(foils     = fill(naca4(0,0,0,9), 2),
                 chords    = [0.7, 0.42],
                 twists    = [0.0, 0.0],
                 spans     = [1.0],
                 dihedrals = [0.],
                 sweeps    = [7.97],
                 position  = [4., 0, 0],
                 angle     = 90.,
                 axis      = [1., 0., 0.])

# Print info
print_info(wing, "Wing")
print_info(htail, "Horizontal Tail")
print_info(vtail, "Vertical Tail")

## WingMesh type
wing_mesh  = WingMesh(wing, [12], 6,
                      span_spacing = Cosine()
                     )
htail_mesh = WingMesh(htail, [12], 6,
                      span_spacing = Cosine()
                     )
vtail_mesh = WingMesh(vtail, [12], 6,
                      span_spacing = Cosine()
                     )

aircraft = ComponentArray(
                          wing  = make_horseshoes(wing_mesh),
                          htail = make_horseshoes(htail_mesh),
                          vtail = make_horseshoes(vtail_mesh)
                         );

## Case
fs  = Freestream(alpha = 0.0,
                 beta  = 0.0,
                 omega = [0., 0., 0.]);

ref = References(speed     = 1.0,
                 density   = 1.225,
                 viscosity = 1.5e-5,
                 area      = projected_area(wing),
                 span      = span(wing),
                 chord     = mean_aerodynamic_chord(wing),
                 location  = mean_aerodynamic_center(wing))

## Evaluate cases
using Setfield

αs = -5:0.5:5
results = permutedims(combinedimsview(
    map(αs) do α
        fst = @set fs.alpha = deg2rad(α)
        sys = solve_case(aircraft, fst, ref)
        [ α; farfield(sys); nearfield(sys) ]
    end
))

## Generate DataFrame
data = DataFrame(results, :auto)
rename!(data, [:α, :CD_ff, :CY_ff, :CL_ff, :CD_nf, :CY_nf, :CL_nf, :Cl, :Cm, :Cn])

##
# pyplot(dpi = 300)
plot(xlabel = "CD", ylabel = "CL", title = "Drag Polar")
plot!(data[!,"CD_ff"], data[!,"CL_ff"], label = :none, marker = :dot)


text:
using Distributions
using StatsBase
using Base.Test

# Weights

wk = [0.20; 0.50; 0.25; 0.05]

wk = wk / sum(wk)

# Perform resampling

res = Int64[]

for i = 1:100000

    append!(res, resample(wk))

end

# Check results

nr = counts(res, 1:maximum(res))

@test all(wk .== round(nr / sum(nr), 2))


text:
export removeZeros, complex2real, real2complex

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

function removeZeros( A::Array )
# Make array shorter by removing all zero values.
A = A[A .!= 0]
return A
end # function removeZeros

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

function complex2real( a::Array{ComplexF64} )
# Convert from [a1 a2 ...] to [ ar1 ai1 ar2 ai2 ... ]
n = length(a)
b = Array{Float64}(undef,n*2)

j = 1
for i = 1:n
   b[j], b[j+1] = reim(a[i])
   j += 2
end  # i

return b
end # function complex2real

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

function real2complex( b::Array{Float64} )
# Convert from [ br1 bi1 br2 bi2 ... ] to [b1 b2 ...]
n = length(b)
a = Array{ComplexF64}(undef,div(n,2))

j = 1
for i = 1:2:n
   a[j] = complex( b[i], b[i+1] )
   j += 1
end  # i

return a
end  # function real2complex


text:
# ClimateMachine driver configurations
#
# Contains helper functions to establish simulation configurations to be
# used with the ClimateMachine driver. Currently:
# - AtmosLESConfiguration
# - AtmosGCMConfiguration
# - OceanBoxGCMConfiguration
# - SingleStackConfiguration
#
# User-customized configurations can use these as templates.

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius

abstract type AbstractSolverType end

struct ExplicitSolverType <: AbstractSolverType
    solver_method::Function
    ExplicitSolverType(; solver_method = LSRK54CarpenterKennedy) =
        new(solver_method)
end

struct IMEXSolverType <: AbstractSolverType
    linear_model::Type
    linear_solver::Type
    solver_method::Function
    function IMEXSolverType(;
        # FIXME: this is Atmos-specific
        linear_model = AtmosAcousticGravityLinearModel,
        linear_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
    )
        return new(linear_model, linear_solver, solver_method)
    end
end

struct MultirateSolverType <: AbstractSolverType
    linear_model::Type
    solver_method::Type
    slow_method::Function
    fast_method::Function
    timestep_ratio::Int
    function MultirateSolverType(;
        # FIXME: this is Atmos-specific
        linear_model = AtmosAcousticGravityLinearModel,
        solver_method = MultirateRungeKutta,
        slow_method = LSRK54CarpenterKennedy,
        fast_method = LSRK54CarpenterKennedy,
        timestep_ratio = 100,
    )
        return new(
            linear_model,
            solver_method,
            slow_method,
            fast_method,
            timestep_ratio,
        )
    end
end

DefaultSolverType = IMEXSolverType

abstract type ConfigSpecificInfo end
struct AtmosLESSpecificInfo <: ConfigSpecificInfo end
struct AtmosGCMSpecificInfo{FT} <: ConfigSpecificInfo
    domain_height::FT
    nelem_vert::Int
    nelem_horz::Int
end
struct OceanBoxGCMSpecificInfo <: ConfigSpecificInfo end
struct SingleStackSpecificInfo <: ConfigSpecificInfo end

"""
    ClimateMachine.DriverConfiguration

Collects all parameters necessary to set up a ClimateMachine simulation.
"""
struct DriverConfiguration{FT}
    config_type::ClimateMachineConfigType

    name::String
    N::Int
    array_type
    solver_type::AbstractSolverType
    #
    # Model details
    bl::BalanceLaw
    #
    # execution details
    mpicomm::MPI.Comm
    #
    # mesh details
    grid::DiscontinuousSpectralElementGrid
    #
    # DGModel details
    numerical_flux_first_order::NumericalFluxFirstOrder
    numerical_flux_second_order::NumericalFluxSecondOrder
    numerical_flux_gradient::NumericalFluxGradient
    #
    # configuration-specific info
    config_info::ConfigSpecificInfo

    function DriverConfiguration(
        config_type,
        name::String,
        N::Int,
        FT,
        array_type,
        solver_type::AbstractSolverType,
        bl::BalanceLaw,
        mpicomm::MPI.Comm,
        grid::DiscontinuousSpectralElementGrid,
        numerical_flux_first_order::NumericalFluxFirstOrder,
        numerical_flux_second_order::NumericalFluxSecondOrder,
        numerical_flux_gradient::NumericalFluxGradient,
        config_info::ConfigSpecificInfo,
    )
        return new{FT}(
            config_type,
            name,
            N,
            array_type,
            solver_type,
            bl,
            mpicomm,
            grid,
            numerical_flux_first_order,
            numerical_flux_second_order,
            numerical_flux_gradient,
            config_info,
        )
    end
end

function print_model_info(model)
    msg = "Model composition\n"
    for key in fieldnames(typeof(model))
        msg =
            msg * @sprintf(
                "    %s = %s\n",
                string(key),
                string((getproperty(model, key)))
            )
    end
    @info msg
end

function AtmosLESConfiguration(
    name::String,
    N::Int,
    (Δx, Δy, Δz)::NTuple{3, FT},
    xmax::FT,
    ymax::FT,
    zmax::FT,
    param_set::AbstractParameterSet,
    init_LES!;
    xmin = zero(FT),
    ymin = zero(FT),
    zmin = zero(FT),
    array_type = ClimateMachine.array_type(),
    solver_type = IMEXSolverType(linear_solver = SingleColumnLU),
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        init_state_conservative = init_LES!,
    ),
    mpicomm = MPI.COMM_WORLD,
    boundary = ((0, 0), (0, 0), (1, 2)),
    periodicity = (true, true, false),
    meshwarp = (x...) -> identity(x),
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
) where {FT <: AbstractFloat}

    print_model_info(model)

    brickrange = (
        grid1d(xmin, xmax, elemsize = Δx * N),
        grid1d(ymin, ymax, elemsize = Δy * N),
        grid1d(zmin, zmax, elemsize = Δz * N),
    )
    topology = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = N,
        meshwarp = meshwarp,
    )

    @info @sprintf(
        """
Establishing Atmos LES configuration for %s
    precision        = %s
    polynomial order = %d
    domain           = %.2f m x%.2f m x%.2f m
    resolution       = %dx%dx%d
    MPI ranks        = %d
    min(Δ_horz)      = %.2f m
    min(Δ_vert)      = %.2f m""",
        name,
        FT,
        N,
        xmax,
        ymax,
        zmax,
        Δx,
        Δy,
        Δz,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        AtmosLESConfigType(),
        name,
        N,
        FT,
        array_type,
        solver_type,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        AtmosLESSpecificInfo(),
    )
end

function AtmosGCMConfiguration(
    name::String,
    N::Int,
    (nelem_horz, nelem_vert)::NTuple{2, Int},
    domain_height::FT,
    param_set::AbstractParameterSet,
    init_GCM!;
    array_type = ClimateMachine.array_type(),
    solver_type = DefaultSolverType(),
    model = AtmosModel{FT}(
        AtmosGCMConfigType,
        param_set;
        init_state_conservative = init_GCM!,
    ),
    mpicomm = MPI.COMM_WORLD,
    meshwarp::Function = cubedshellwarp,
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
) where {FT <: AbstractFloat}

    print_model_info(model)

    _planet_radius::FT = planet_radius(param_set)
    vert_range = grid1d(
        _planet_radius,
        FT(_planet_radius + domain_height),
        nelem = nelem_vert,
    )

    topology = StackedCubedSphereTopology(mpicomm, nelem_horz, vert_range)

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = N,
        meshwarp = meshwarp,
    )

    @info @sprintf(
        """
Establishing Atmos GCM configuration for %s
    precision        = %s
    polynomial order = %d
    #horiz elems     = %d
    #vert elems      = %d
    domain height    = %.2e m
    MPI ranks        = %d
    min(Δ_horz)      = %.2f m
    min(Δ_vert)      = %.2f m""",
        name,
        FT,
        N,
        nelem_horz,
        nelem_vert,
        domain_height,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        AtmosGCMConfigType(),
        name,
        N,
        FT,
        array_type,
        solver_type,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        AtmosGCMSpecificInfo(domain_height, nelem_vert, nelem_horz),
    )
end

function OceanBoxGCMConfiguration(
    name::String,
    N::Int,
    (Nˣ, Nʸ, Nᶻ)::NTuple{3, Int},
    model::HydrostaticBoussinesqModel;
    FT = Float64,
    array_type = ClimateMachine.array_type(),
    solver_type = ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    ),
    mpicomm = MPI.COMM_WORLD,
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
    periodicity = (false, false, false),
    boundary = ((1, 1), (1, 1), (2, 3)),
)

    brickrange = (
        range(FT(0); length = Nˣ + 1, stop = model.problem.Lˣ),
        range(FT(0); length = Nʸ + 1, stop = model.problem.Lʸ),
        range(FT(-model.problem.H); length = Nᶻ + 1, stop = 0),
    )

    topology = StackedBrickTopology(
        mpicomm,
        brickrange;
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = N,
    )

    return DriverConfiguration(
        OceanBoxGCMConfigType(),
        name,
        N,
        FT,
        array_type,
        solver_type,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        OceanBoxGCMSpecificInfo(),
    )
end

function SingleStackConfiguration(
    name::String,
    N::Int,
    nelem_vert::Int,
    zmax::FT,
    param_set::AbstractParameterSet,
    model::BalanceLaw;
    zmin = zero(FT),
    array_type = ClimateMachine.array_type(),
    solver_type = ExplicitSolverType(),
    mpicomm = MPI.COMM_WORLD,
    boundary = ((0, 0), (0, 0), (1, 2)),
    periodicity = (true, true, false),
    meshwarp = (x...) -> identity(x),
    numerical_flux_first_order = RusanovNumericalFlux(),
    numerical_flux_second_order = CentralNumericalFluxSecondOrder(),
    numerical_flux_gradient = CentralNumericalFluxGradient(),
) where {FT <: AbstractFloat}

    print_model_info(model)

    xmin, xmax = zero(FT), one(FT)
    ymin, ymax = zero(FT), one(FT)
    brickrange = (
        grid1d(xmin, xmax, nelem = 1),
        grid1d(ymin, ymax, nelem = 1),
        grid1d(zmin, zmax, nelem = nelem_vert),
    )
    topology = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = periodicity,
        boundary = boundary,
    )

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = array_type,
        polynomialorder = N,
        meshwarp = meshwarp,
    )

    @info @sprintf(
        """
Establishing single stack configuration for %s
    precision        = %s
    polynomial order = %d
    domain           = %.2f m x%.2f m x%.2f m
    #vert elems      = %d
    MPI ranks        = %d
    min(Δ_horz)      = %.2f m
    min(Δ_vert)      = %.2f m""",
        name,
        FT,
        N,
        xmax,
        ymax,
        zmax,
        nelem_vert,
        MPI.Comm_size(mpicomm),
        min_node_distance(grid, HorizontalDirection()),
        min_node_distance(grid, VerticalDirection())
    )

    return DriverConfiguration(
        SingleStackConfigType(),
        name,
        N,
        FT,
        array_type,
        solver_type,
        model,
        mpicomm,
        grid,
        numerical_flux_first_order,
        numerical_flux_second_order,
        numerical_flux_gradient,
        SingleStackSpecificInfo(),
    )
end


text:
module HyperFlows

# We need a function where input ranges from 0 to 1, and output is an x and y coordinate, which we can integrate over

struct LimFunc{A<:Real,B}
    startx::A
    endx::A
    func::B
end

function geomConstructor(LimFuncs::LimFunc...)
    totalLength = sum(map(x -> abs(x.startx-x.endx), LimFuncs ))
    return function(t)
        t = t*totalLength
        runningTotal = 0
        for func in LimFuncs
            runningTotal = runningTotal + abs(func.startx - func.endx)
            if t < runningTotal
                x = func.startx + (t - runningTotal + abs(func.startx - func.endx )) / (func.endx-func.startx)
                return x
            end
        end
    end#returns funtion that maps a value t to x and y values
end

end # module

# scratch code
first = LimFunc(2,3,x->x-1)
second = LimFunc(3,2,x->-2*x+8)
third = LimFunc(2,1,x->2*x)
forth = LimFunc(1,2,x->-x+3)

g = geomConstructor(first, second, third, forth)
g(0.74)


text:
"""
`process_sessions`
"""

function process_sessions(DataIndex::DataFrames.AbstractDataFrame)
    c=0
    b=0
    pokes = DataFrame()
    bouts = DataFrame()
    streaks = DataFrame()
    for i=1:size(DataIndex,1)
        path = DataIndex[i,:Bhv_Path]
        session = DataIndex[i,:Session]
        filetosave = DataIndex[i,:Preprocessed_Path]
        if ~isfile(filetosave)
            pokes_data = process_pokes(path)
            FileIO.save(filetosave,pokes_data)
            bouts_data = process_bouts(pokes_data)
            streaks_data = process_streaks(pokes_data)
            b=b+1
        else
            pokes_data = FileIO.load(filetosave)|> DataFrame
            booleans=[:Reward,:Stim,:Wall,:Correct,:Stim_Day,:LastPoke]#columns to convert to Bool removed :Side,:SideHigh
            for x in booleans
                pokes_data[!,x] = eltype(pokes_data[!,x]) == Bool ? pokes_data[!,x] : occursin.(pokes_data[!,x],"true")
            end
            bouts_data = process_bouts(pokes_data)
            streaks_data = process_streaks(pokes_data)
            c=c+1
        end
        if isempty(pokes)
            pokes = pokes_data
            bouts = bouts_data
            streaks = streaks_data
        else
            try
                append!(pokes, pokes_data)
                append!(bouts,bouts_data)
                append!(streaks, streaks_data)
            catch
                println(DataIndex[i,:Bhv_Path])
                append!(pokes, pokes_data[:, names(pokes)])
                append!(bouts, bouts_data[:, names(pokes)])
                append!(streaks, streaks_data[:, names(streaks)])
            end
        end
    end
    println("Existing file = ",c," Preprocessed = ",b)
    return pokes, bouts, streaks
end

"""
`create_exp_dataframes`
"""

function create_exp_dataframes(DataIndex::DataFrames.AbstractDataFrame)
    exp_dir = DataIndex[1,:Saving_Path]
    pokes, bouts, streaks = process_sessions(DataIndex)
    exp_calendar = by(pokes,:MouseID) do dd
        Flipping.create_exp_calendar(dd,:Day)
    end
    protocol_calendar = by(pokes,:MouseID) do dd
        Flipping.create_exp_calendar(dd,:Day,:Protocol)
    end
    if !any(protocol_calendar[:,:Flexi])
        select!(protocol_calendar,DataFrames.Not([:Manipulation,:Flexi]))
    end
    pokes = add_exp_calendar(pokes,exp_calendar,protocol_calendar)
    pokes = Flipping.check_fiberlocation(pokes,exp_dir)
    filetosave = joinpath(exp_dir,"pokes"*splitdir(exp_dir)[end]*".csv")
    CSVFiles.save(filetosave,pokes)

    bouts = add_exp_calendar(bouts,exp_calendar,protocol_calendar)
    bouts = Flipping.check_fiberlocation(bouts,exp_dir)
    filetosave = joinpath(exp_dir,"bouts"*splitdir(exp_dir)[end]*".csv")
    CSVFiles.save(filetosave,bouts)

    streaks = add_exp_calendar(streaks,exp_calendar,protocol_calendar)
    streaks = Flipping.check_fiberlocation(streaks,exp_dir)
    filetosave = joinpath(exp_dir,"streaks"*splitdir(exp_dir)[end]*".jld2")
    @save filetosave streaks
    simple = DataFrames.select(streaks,DataFrames.Not(:PokeSequence))
    filetosave = joinpath(exp_dir,"streaks"*splitdir(exp_dir)[end]*".csv")
    CSVFiles.save(filetosave,simple)
    return pokes, bouts, streaks, DataIndex
end

function create_exp_dataframes(Directory_path::String,Exp_type::String,Exp_name::String, Mice_suffix::String)
    DataIndex = Flipping.find_behavior(Directory_path, Exp_type, Exp_name,Mice_suffix)
    create_exp_dataframes(DataIndex)
end

function create_exp_dataframes(Raw_data_dir::String)
    DataIndex = Flipping.find_behavior(Raw_data_dir)
    create_exp_dataframes(DataIndex)
end


text:
# The SLEEF.jl package is licensed under the MIT "Expat" License:

# > Copyright (c) 2016: Mustafa Mohamad and other contributors:
# >
# > https://github.com/musm/SLEEF.jl/graphs/contributors
# >
# > Permission is hereby granted, free of charge, to any person obtaining a copy
# > of this software and associated documentation files (the "Software"), to deal
# > in the Software without restriction, including without limitation the rights
# > to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# > copies of the Software, and to permit persons to whom the Software is
# > furnished to do so, subject to the following conditions:
# >
# > The above copyright notice and this permission notice shall be included in all
# > copies or substantial portions of the Software.
# >
# > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# > IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# > FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# > AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# > LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# > OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# > SOFTWARE.
# >

# SLEEF.jl includes ported code from the following project

#     - [SLEEF](https://github.com/shibatch/SLEEF) [public domain] Author Naoki Shibata

using Base.Math: significand_bits

isnzero(x::T) where {T <: AbstractFloat} = signbit(x)
ispzero(x::T) where {T <: AbstractFloat} = !signbit(x)

# function cmpdenorm(x::Tx, y::Ty) where {Tx <: AbstractFloat, Ty <: AbstractFloat}
#     sizeof(Tx) < sizeof(Ty) ? y = Tx(y) : x = Ty(x) # cast larger type to smaller type
#     (isnan(x) && isnan(y)) && return true
#     (isnan(x) || isnan(y)) && return false
#     (isinf(x) != isinf(y)) && return false
#     (x == Tx(Inf)  && y == Ty(Inf))  && return true
#     (x == Tx(-Inf) && y == Ty(-Inf)) && return true
#     if y == 0
#         (ispzero(x) && ispzero(y)) && return true
#         (isnzero(x) && isnzero(y)) && return true
#         return false
#     end
#     (!isnan(x) && !isnan(y) && !isinf(x) && !isinf(y)) && return sign(x) == sign(y)
#     return false
# end

# the following compares the ulp between x and y.
# First it promotes them to the larger of the two types x,y
infh(::Type{Float64}) = 1e300
infh(::Type{Float32}) = 1e37
function countulp(::Type{T}, __x, __y) where {T}
    _x, _y = promote(__x, __y)
    x, y = convert(T, _x), convert(T, _y) # Cast to smaller type

    ulpc = convert(T, abs(_x - _y) / ulp(y))
    nanulp = VectorizationBase.ifelse(isnan(x) ⊻ isnan(y), T(10000), T(0))
    infulp = VectorizationBase.ifelse((sign(x) == sign(y)) & (abs(y) > infh(T)), T(0), T(10001))

    ulpc = VectorizationBase.ifelse(isinf(x), infulp, VectorizationBase.ifelse(isfinite(y), ulpc, T(10003)))
    ulpc = VectorizationBase.ifelse(isnan(x) | isnan(y), nanulp, ulpc)
    ulpc = VectorizationBase.ifelse(iszero(y), VectorizationBase.ifelse(iszero(x), T(0), T(10002)), ulpc)
    return ulpc
end

DENORMAL_MIN(::Type{Float64}) = 2.0^-1074
DENORMAL_MIN(::Type{Float32}) = 2f0^-149

function ulp(x::Union{<:VectorizationBase.AbstractSIMD{<:Any,T},T}) where {T<:AbstractFloat}
    e = exponent(x)
    ulpc = max(VectorizationBase.vscalef(T(1.0), e - significand_bits(T)), DENORMAL_MIN(T))
    ulpc = VectorizationBase.ifelse(x == T(0.0), DENORMAL_MIN(T), ulpc)
    return ulpc
end

countulp(x::T, y::T) where {T <: AbstractFloat} = countulp(T, x, y)
countulp(x::VectorizationBase.AbstractSIMD{W,T}, y::VectorizationBase.AbstractSIMD{W,T}) where {W,T <: AbstractFloat} = countulp(T, x, y)

# test the accuracy of a function where fun_table is a Dict mapping the function you want
# to test to a reference function
# xx is an array of values (which may be tuples for multiple arugment functions)
# tol is the acceptable tolerance to test against
function test_acc(f1, f2, T, xx, tol, ::StaticInt{W} = pick_vector_width(T); debug = false, tol_debug = 5) where {W}
    @testset "accuracy $(f1)" begin

        reference = map(f2 ∘ big, xx)
        comp = similar(xx)
        i = 0
        spc = VectorizationBase.zstridedpointer(comp); spx = VectorizationBase.zstridedpointer(xx);
        GC.@preserve xx comp begin
            while i < length(xx)
                vstore!(spc, f1(vload(spx, (MM{W}(i),))), (MM{W}(i),))
                i += W
            end
        end
        rmax = 0.0
        rmean = 0.0
        xmax = map(zero, first(xx))
        for i ∈ eachindex(xx)
            q = comp[i]
            c = reference[i]
            u = countulp(T, q, c)
            rmax = max(rmax, u)
            xmax = rmax == u ? xx[i] : xmax
            rmean += u
            if debug && u > tol_debug
                @show f1, q, f2, T(c), x, T(c)
            end
        end
        rmean = rmean / length(xx)

        t = @test trunc(rmax, digits=1) <= tol

        fmtxloc = isa(xmax, Tuple) ? join(xmax, ", ") : string(xmax)
        println(
            rpad(f1, 18, " "), ": max ", rmax,
            rpad(" at x = " * fmtxloc, 40, " "),  ": mean ", rmean
        )

    end
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 04911b34-229d-11eb-1658-3d8632a4e006
using Distributions, StatsBase, Plots, OTSIN

# ╔═╡ bc3b1416-2355-11eb-16b4-c1e191b1dc94
md"""
# OT solutions arise from local choices

In this notebook, I will show that the optimal transportation model arises from local choises of species interactions. It suffices to do local choices of species based on the **utility** and Metropolis heuristic. This leads to behaviour for which the optimal transportation problem is an exact solution.
"""

# ╔═╡ 8547775c-2359-11eb-1a41-673a988bfdca
md"First, we define a utility matrix."

# ╔═╡ 6144ef0e-2333-11eb-17a9-0392d858553e
begin
	Ula = zeros(30, 25)
	Ula[1:10, 1:12] .= 2.0
	Ula[9:20, 10:20] .= 1.5
	Ula[22:30, 20:end] .= 1.0
	Ula[11:end, 1:10] .= 0.5
	heatmap(Ula)
end

# ╔═╡ d4fe37fa-2335-11eb-31f6-bfda1afcd852
Usm = [3 3 3 1 0;
	   3 2.5 2 1 0.5;
	   2.5 2 2 0 0;
       1  1 0 0 1;
       1 0 0 0.5 1]

# ╔═╡ 60672f2a-229d-11eb-1f7f-a59e191a8a65
#M = 10*[ones(15, 10) zeros(15, 10); 0.1ones(10, 12) 0.5ones(10, 8)];
M = Usm

# ╔═╡ 3a16e734-229d-11eb-3187-d721246abe92
n, m = size(M)

# ╔═╡ bd664532-229d-11eb-19aa-1d5847ad8687
a = ones(n) / n

# ╔═╡ c99500bc-229d-11eb-1459-fb202f664fba
b = range(10, 1, length=m) |> x -> x/sum(x)

# ╔═╡ eac665fa-229d-11eb-1404-73687f0ea99e
multinomr(n, p) = rand(Multinomial(n, p))

# ╔═╡ 51c9802a-229e-11eb-0800-d180f7c63640
sample(1:n, Weights(a), 2)

# ╔═╡ 214d7c6c-229e-11eb-2b54-cf6068a5ff69

function simulated_annealing!(P, M, λ; nsteps=10000)
	n, m = size(P)
	n_interactions = sum(P)
	rowweigths = Weights(sum(P, dims=2)[:])
	for t in 1:nsteps
		i1, i2 = sample(1:n, rowweigths, 2)
		j1 = sample(1:m, Weights(P[i1,:]))
		j2 = sample(1:m, Weights(P[i2,:]))
		ΔM = (M[i1,j2] + M[i2, j1]) - (M[i1,j1] + M[i2,j2])
		if ΔM > 0.0 || rand() ≤ exp(ΔM * λ)
			P[i1,j1] -= 1
			P[i2,j2] -= 1
			P[i1,j2] += 1
			P[i2,j1] += 1
		end
	end
	return P
end

# ╔═╡ b85f7bca-229f-11eb-1e15-c7cb10437450
n_interactions = 200

# ╔═╡ c0c50460-229f-11eb-3a58-7d43a336d5e3
a_sampled = multinomr(n_interactions, a)

# ╔═╡ ce10d4b6-229f-11eb-06bd-3720e2a99287
b_sampled = multinomr(n_interactions, b)

# ╔═╡ d9bc7310-229f-11eb-37fe-3f30f4814257
function findP(a_sampled, b_sampled, n_interactions)
	n, m = length(a_sampled), length(b_sampled)
	P = zeros(Int, n, m)
	while n_interactions > 0
		i, j = rand(1:n), rand(1:m)
		if sum(P[i,:]) < a_sampled[i] && sum(P[:,j]) < b_sampled[j]
			P[i,j] += 1
			n_interactions -= 1
		end
	end
	return P
end

# ╔═╡ 7fb496e4-22a0-11eb-1080-ebf1eb4f42bc
P = findP(a_sampled, b_sampled, n_interactions)

# ╔═╡ fd192064-22a0-11eb-1fcc-b11f49ad44ab
simulated_annealing!(P, M, 10; nsteps=1000_000)

# ╔═╡ 2f8e074e-22a1-11eb-136d-2b83877dd889
marginalsheatmap(P)

# ╔═╡ 6896da6a-22a2-11eb-1fcc-839ef4100dc1
Pot = optimaltransport(M, a, b, λ=10, maxitter=100)

# ╔═╡ 806ed552-22a2-11eb-203a-419eee9307a3
marginalsheatmap(Pot)

# ╔═╡ 295134bc-22a3-11eb-3e28-fb55a02936e8
λs = (10).^(-2:0.2:1.5)

# ╔═╡ 76315a66-22ab-11eb-0a31-adee0035db0d
n_experiments = length(λs)

# ╔═╡ 844fa4c2-22ab-11eb-1186-f1b3644cc72f
n_repetitions = 50

# ╔═╡ de07a940-22ab-11eb-0113-83f986874b02
numb_interactions = [100, 1000, 10_000]

# ╔═╡ bf4dc6c6-236e-11eb-17b5-63f843e8cce7
minheatmap(P; kwargs...) = heatmap(P, legend=false, seriescolor=cgrad([:white, :blue]),showaxis=true, xticks=(1:size(P,1), ('a':'z')[1:size(P,1)]),
				yticks=(1:size(P,2), ('A':'Z')[1:size(P,2)]); kwargs...)

# ╔═╡ b6270016-229d-11eb-1d5a-1b124847f727
pM = minheatmap(M, color=cgrad([:white, :red]))

# ╔═╡ f8797b2a-27ea-11eb-2be7-d76ac2a24e4c
savefig(pM, "../figures/M_sa_sim.svg")

# ╔═╡ c1c25be2-236e-11eb-0b1e-27ce2c4a7416
minheatmap(P)

# ╔═╡ 7fba27a0-22a5-11eb-1a5f-85c0b7fe545c
function entropy(P)
	total = sum(P)
	logtotal = log(total)
	h = 0.0
	for Pij in P
		if Pij > 0
			Pij /= total
			h += -Pij * log(Pij)
		end
	end
	return h
end

# ╔═╡ f51dd178-22a7-11eb-286d-756e77b64941
average_utity(P, M) = sum(P .* M) / sum(P)

# ╔═╡ e29aaa28-22b4-11eb-097b-2379ac1d3ee8
average_utity(P, M)

# ╔═╡ f4b75eca-22a2-11eb-3454-b7cf0c58ae0b
begin
	av_utilities = zeros(n_experiments, length(numb_interactions));
for (ind_int, n_int) in enumerate(numb_interactions)
	for rep in 1:n_repetitions
		P = findP(multinomr(n_int, a), multinomr(n_int, b), n_int)
		for (i, λ) in enumerate(λs)
			simulated_annealing!(P, M, λ; nsteps=100000)
			av_utilities[i, ind_int] += average_utity(P, M)
		end

	end

end
av_utilities ./= n_repetitions
end

# ╔═╡ 2dc31e48-2372-11eb-00c0-a3fdd9d08e56
begin
	utilities = zeros(n_experiments, n_repetitions)
	for rep in 1:n_repetitions
		P = findP(multinomr(n_interactions, a), multinomr(n_interactions, b), n_interactions)
		for (i, λ) in enumerate(λs)
			simulated_annealing!(P, M, λ; nsteps=100_000)
			utilities[i, rep] += average_utity(P, M)
		end
	end
	utilities
end

# ╔═╡ 3e1c6016-2374-11eb-3b85-3b419741ec61
μ_ut = mean(utilities, dims=2)

# ╔═╡ 459210c0-2374-11eb-1fab-a123f5b4fe1f
σ_ut = std(utilities, dims=2)

# ╔═╡ 73810eaa-22af-11eb-2417-6d34b137efac
ot_utilites_theo = [average_utity(optimaltransport(M, a, b, λ=λ), M) for λ in λs]

# ╔═╡ 846a5720-22ad-11eb-1adf-d5f29a11edd2
begin
	plambda = plot(λs, [μ_ut μ_ut], fillrange=[μ_ut+1σ_ut μ_ut-1σ_ut], fillalpha=0.4, color=:orange, label=["sampled utitily" ""], lw=4, xscale=:log, xlabel="lambda", ylabel="average utility")
	#plot!(λs, μ_ut, label = "", , color=:orange)
	plot!(λs, ot_utilites_theo, label="OT bound", ls=:dash, lw=4, color=:black)
	plot!(legend=:topleft)
end

# ╔═╡ e91cf2a0-236c-11eb-13b7-99f464bd8c2c
begin
	l = @layout [a{0.60h}; b c d; e f g]
	p_sum = plot(plambda,
		minheatmap(optimaltransport(M, a, b, λ=0.05), title="lambda=0.05"),
		minheatmap(optimaltransport(M, a, b, λ=0.5), title="lambda=0.5"),
		minheatmap(optimaltransport(M, a, b, λ=5), title="lambda=5"),
		minheatmap(simulated_annealing!(P, M, 0.05, nsteps=1000_000), seriescolor=cgrad([:white, :green])),
		minheatmap(simulated_annealing!(P, M, 0.5, nsteps=1000_000), seriescolor=cgrad([:white, :green])),
		minheatmap(simulated_annealing!(P, M, 5, nsteps=1000_000), seriescolor=cgrad([:white, :green])),
		layout=l, size=(600, 800))
	savefig("../figures/OT_MH_approx.svg")
	p_sum
end

# ╔═╡ 594a8728-2374-11eb-01bc-e9bf600f89d3
begin
	pband = plot(λs, av_utilities,
			labels = reshape(["n interactions=$nintr" for nintr in numb_interactions], 1, :), lw=2, xscale=:log, xlabel="lambda", ylabel="average utility")
	plot!(λs, ot_utilites_theo, label="OT bound", ls=:dash, lw=4, color=:black)
	plot!(legend=:topleft)
end

# ╔═╡ 2184453a-22a8-11eb-301f-9553f16a0c4c
sum(Pot)

# ╔═╡ ac9d6a72-22b5-11eb-2060-21112685115f
Ps = reshape(rand(Multinomial(10, vec(Pot))), n, m)

# ╔═╡ 28543ba0-22b6-11eb-3457-937945c42c47
average_utity(Ps, M)

# ╔═╡ Cell order:
# ╠═bc3b1416-2355-11eb-16b4-c1e191b1dc94
# ╠═04911b34-229d-11eb-1658-3d8632a4e006
# ╟─8547775c-2359-11eb-1a41-673a988bfdca
# ╠═3a16e734-229d-11eb-3187-d721246abe92
# ╠═6144ef0e-2333-11eb-17a9-0392d858553e
# ╠═b6270016-229d-11eb-1d5a-1b124847f727
# ╠═f8797b2a-27ea-11eb-2be7-d76ac2a24e4c
# ╠═d4fe37fa-2335-11eb-31f6-bfda1afcd852
# ╠═60672f2a-229d-11eb-1f7f-a59e191a8a65
# ╠═bd664532-229d-11eb-19aa-1d5847ad8687
# ╠═c99500bc-229d-11eb-1459-fb202f664fba
# ╠═eac665fa-229d-11eb-1404-73687f0ea99e
# ╠═51c9802a-229e-11eb-0800-d180f7c63640
# ╠═214d7c6c-229e-11eb-2b54-cf6068a5ff69
# ╠═b85f7bca-229f-11eb-1e15-c7cb10437450
# ╠═c0c50460-229f-11eb-3a58-7d43a336d5e3
# ╠═ce10d4b6-229f-11eb-06bd-3720e2a99287
# ╠═d9bc7310-229f-11eb-37fe-3f30f4814257
# ╠═7fb496e4-22a0-11eb-1080-ebf1eb4f42bc
# ╠═fd192064-22a0-11eb-1fcc-b11f49ad44ab
# ╠═2f8e074e-22a1-11eb-136d-2b83877dd889
# ╠═e29aaa28-22b4-11eb-097b-2379ac1d3ee8
# ╠═6896da6a-22a2-11eb-1fcc-839ef4100dc1
# ╠═806ed552-22a2-11eb-203a-419eee9307a3
# ╠═295134bc-22a3-11eb-3e28-fb55a02936e8
# ╠═76315a66-22ab-11eb-0a31-adee0035db0d
# ╠═844fa4c2-22ab-11eb-1186-f1b3644cc72f
# ╠═de07a940-22ab-11eb-0113-83f986874b02
# ╠═f4b75eca-22a2-11eb-3454-b7cf0c58ae0b
# ╠═2dc31e48-2372-11eb-00c0-a3fdd9d08e56
# ╠═3e1c6016-2374-11eb-3b85-3b419741ec61
# ╠═459210c0-2374-11eb-1fab-a123f5b4fe1f
# ╠═73810eaa-22af-11eb-2417-6d34b137efac
# ╠═846a5720-22ad-11eb-1adf-d5f29a11edd2
# ╠═594a8728-2374-11eb-01bc-e9bf600f89d3
# ╠═e91cf2a0-236c-11eb-13b7-99f464bd8c2c
# ╠═bf4dc6c6-236e-11eb-17b5-63f843e8cce7
# ╠═c1c25be2-236e-11eb-0b1e-27ce2c4a7416
# ╠═7fba27a0-22a5-11eb-1a5f-85c0b7fe545c
# ╠═f51dd178-22a7-11eb-286d-756e77b64941
# ╠═2184453a-22a8-11eb-301f-9553f16a0c4c
# ╠═ac9d6a72-22b5-11eb-2060-21112685115f
# ╠═28543ba0-22b6-11eb-3457-937945c42c47

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
using Fresa
println("Running Fresa tests")


text:
# pure-julia (currently) implementation of basic position logic:
# - a Position() instance represents a particular "arrangement"
# - move!(...) updates the position to include the new move
# - groups of stones and spaces are tracked
# - the implementation is written in a style that should be easy to move
#   to a GPU (OpenCL)
# - additional info useful for presenting to a neural net (liek flood
#   fill distances) is included

# https://en.wikipedia.org/wiki/Komidashi
komi = Float32(7.5)

# --- data structures

"""a point on the board - empty, or a stone."""
@enum Point empty=0 black=1 white=-1

"""if the point is a stone, give a stone of the other colour.  if
empty, return empty."""
other(t::Point) = t == empty ? empty : (t == black ? white : black)

"""a row of 19 trits fits nicely in 32 bits binary."""
typealias Row UInt32

"""a row of n zeros (signed trits)."""
emptyrow(n) = Row(sum([3^(n-1) for n in 1:n]))

"""a trit-map of the entire board."""
@auto_hash_equals immutable Board{N}
    row::Vector{Row}
    Board() = new(fill(emptyrow(N), N))
    Board(b::Board) = new(copy(b.row))
end
# https://groups.google.com/forum/#!topic/julia-users/-n3xsEHxwAI
Board{N}(b::Board{N}) = Board{N}(b)

"""identify points common to a group and provide basic stats."""
@auto_hash_equals immutable Groups{N}
    index::Array{UInt8, 2}
    # do we need to care about overflow here?
    size::Array{UInt8, 1}
    lives::Array{UInt8, 1}
    Groups() = new(zeros(UInt8, N, N), zeros(UInt8, 255), zeros(UInt8, 255))
    Groups(g::Groups{N}) = new(copy(g.index), copy(g.size), copy(g.lives))
end
Groups{N}(g::Groups{N}) = Groups{N}(g)

"""a map of cartesian signed distances to the nearest non-empty point,
where the sign indicates whether the nearest point is black or white
(via the Point enumeration).

distances are offset by 1, so a value of +/-1 indicates that the point
contains a stone.  a value of zero is used initially (replaced on
first move by flood_to_point!) and as a marker for deleted groups
(replaced by flood_dead_group!).

this contains almost all the information about the state of play (the
+/-1 values duplicate Board).  it may be used as input to a neural
net."""
@auto_hash_equals immutable Flood{N}
    distance::Array{Int8, 2}
    Flood() = new(zeros(Int8, N, N))
    Flood(f::Flood{N}) = new(copy(f.distance))
end
Flood{N}(g::Flood{N}) = Flood{N}(g)

@auto_hash_equals immutable Space{N}
    border::Array{Int8, 2}
    index::Array{UInt8, 2}
    Space() = new(zeros(Int8, N, N), zeros(UInt8, N, N))
    Space(s::Space{N}) = new(copy(s.border), copy(s.index))
end
Space{N}(g::Space{N}) = Space{N}(g)

@auto_hash_equals type Details
    spaces::Int
    prisoners::Int
    Details() = new(0, 0)
    Details(d::Details) = new(d.spaces, d.prisoners)
end

@auto_hash_equals type Stats
    prev::Nullable{Tuple{Int,Int}}
    score::Float32
    moves::Int
    owned::Int
    stones::Int
    colour::Dict{Point, Details}
    Stats() = new(Nullable{Tuple{Int,Int}}(), -komi, 0, 0, 0, Dict{Point, Details}(black=>Details(), white=>Details()))
    Stats(s::Stats) = new(s.prev, s.score, s.moves, s.owned, s.stones, Dict([x=>Details(s.colour[x]) for x in (black, white)]))
end

"""a single position (implicitly, in a search tree).  combines Board,
Groups, Flood, Space and Stats."""
@auto_hash_equals immutable Position{N}
    board::Board{N}
    groups::Groups{N}
    flood::Flood{N}
    space::Space{N}
    stats::Stats
    Position() = new(Board{N}(), Groups{N}(), Flood{N}(), Space{N}(), Stats())
    Position(p::Position) =
        new(Board(p.board), Groups(p.groups), Flood(p.flood), Space(p.space), Stats(p.stats))
end
Position{N}(g::Position{N}) = Position{N}(g)
Position() = Position{19}()

# --- points

"""extract the point at a given (x, y)."""
point(r::Row, x) = Point(mod(div(r, 3^(x-1)), 3)-1)

"""extract the point at a given (x, y)."""
point(b::Board, x, y) = point(b.row[y], x)

"""extract the point at a given (x, y)."""
point(p::Position, x, y) = point(p.board, x, y)

# --- display

"""the location of board markers (dots)."""
markers(n) = Set([Integer(floor(n/4)), (n+1)/2, n+1-Integer(floor(n/4))])

function fmtpoint(r::Row, y, n)
    m = markers(n)
    fmt(x) = (y in m && x in m ? "O.X" : "O X")[2+Int(point(r, x))]
    # extra space aligns with other components that use 2 characters
    # per point
    string(" ", join(map(fmt, 1:n), " "))
end
fmtpoint{N}(b::Board{N}) = map(y -> fmtpoint(b.row[y], y, N), N:-1:1)

Base.print(io::IO, b::Board) = print(io, join(fmtpoint(b), "\n"))

function fmtindex(index::Array{UInt8, 2}, y, n)
    m = markers(n)
    function fmt(x)
        if index[x, y] > 0
            @sprintf("%02x", index[x, y])
        elseif y in m && x in m
            "__"
        else
            "  "
        end
    end
    join(map(fmt, 1:n))
end
fmtindex{N}(g::Groups{N}) = map(y -> fmtindex(g.index, y, N), N:-1:1)

stats_per_row = 4
function fmtstats(g::Groups, row)
    function fmt(col)
        i = col + (row-1) * stats_per_row
        @sprintf("%02x:%2d/%2d", i, g.lives[i], g.size[i])
    end
    join(map(fmt, 1:stats_per_row), ", ")
end
fmtstats{N}(g::Groups{N}) = map(r -> fmtstats(g, r), 1:N)

function Base.print(io::IO, g::Groups)
    print(io, join(fmtindex(g), "\n"))
    print(io, "\n")
    print(io, join(fmtstats(g), "\n"))
end

function fmtdistance{N}(f::Flood{N}, y)
    function fmt(x)
        d = f.distance[x, y]
        if abs(d) < 10
            @sprintf("%2d", d)
        elseif d > 0
            string(" ", Char(Int('a') + d - 10))
        else
            string(" ", Char(Int('A') - d - 10))
        end
    end
    join(map(fmt, 1:1N))
end
fmtdistance{N}(f::Flood{N}) = map(y -> fmtdistance(f, y), N:-1:1)

Base.print(io::IO, f::Flood) = print(io, join(fmtdistance(f), "\n"))

fmtindex{N}(s::Space{N}) = map(y -> fmtindex(s.index, y, N), N:-1:1)

function fmtborder{N}(s::Space{N}, y)
    fmt(x) = " XO."[s.border[x, y]+1]
    string(" ", join(map(fmt, 1:N), " "))
end
fmtborder{N}(s::Space{N}) = map(y -> fmtborder(s, y), N:-1:1)

function Base.print(io::IO, s::Space)
    print(io, join(fmtborder(s), "\n"))
    print(io, "\n")
    print(io, join(fmtindex(s), "\n"))
end

Base.print(io::IO, s::Stats) =
@printf(io, " %5.1f  [black  pr=%-2d sp=%-3d] [white  pr=%-2d sp=%-3d]",
        s.score,
        s.colour[black].prisoners, s.colour[black].spaces,
        s.colour[white].prisoners, s.colour[white].spaces)

function Base.print(io::IO, p::Position)
    print(io,
          join(map(x -> join(x, "  "),
                   zip(fmtpoint(p.board),
                       fmtindex(p.groups))), "\n"))
    print(io, "\n\n")
    print(io,
          join(map(x -> join(x, "  "),
                   zip(fmtdistance(p.flood),
                       fmtstats(p.groups))), "\n"))
    print(io, "\n\n")
    print(io,
          join(map(x -> join(x, "  "),
                   zip(fmtborder(p.space),
                       fmtindex(p.space))), "\n"))
    print(io, "\n\n")
    print(io, p.stats)
end

# --- update state on move

function k_lowest_unused(k, a::Array{UInt8, 2}, n)
    available = collect(UInt8, 1:255)
    result = zeros(UInt8, k)
    function extract(result, value)
        if value > 0
            for i in 1:k
                if result[i] == 0
                    result[i] = value
                    break
                end
            end
        end
        result
    end
    @forall_fold i j n extract result available begin
        if a[i, j] > 0
            available[a[i, j]] = 0
        end
    end
end

"""merge the group at (x,y) with newgroup."""
function merge_group!{N}(g::Groups{N}, newgroup, x, y)
    oldgroup = g.index[x, y]
    if newgroup != oldgroup
        g.size[newgroup], g.size[oldgroup] = g.size[newgroup] + g.size[oldgroup], 0
        @forall i j N begin
            if g.index[i, j] == oldgroup
                g.index[i, j] = newgroup
            end
        end
    end
end

"""delete group g: remove it from Groups.index, sets Flood.distance to
zero (these are later replaced in flood_dead_group!) and Space.index
to the next available value (patched up in fix_space!)."""
function delete_group!{N}(p::Position{N}, t, g)
    space = k_lowest_unused(1, p.space.index, N)[1]
    p.stats.colour[t].prisoners = p.stats.colour[t].prisoners + p.groups.size[g]
    p.groups.size[g] = 0
    @forall i j N begin
        if p.groups.index[i, j] == g
            # remove stone
            move!(p.board, empty, i, j)
            # erase group index
            p.groups.index[i, j] = 0
            p.space.index[i, j] = space
            # set flood to special value (see flood_dead_group!)
            p.flood.distance[i, j] = 0
        end
    end
end

"""replace Flood.distance values that have been zeroed."""
function flood_dead_group!{N}(f::Flood{N}, b::Board{N}, t::Point)
    while ! @forall_fold i j N (a,b) -> (a || b!=0) false f.distance begin
        if f.distance[i, j] == 0
           # find closest neighbour
           mind = 256
           @forneighbours i j N ii jj begin
               d = f.distance[ii, jj]
               if d != 0
                   mind = min(mind, abs(d))
               end
           end
           # if we have a neighbour, add one to that
           if mind < 256
               f.distance[i, j] = Int(t) * (mind+1)
           end
       end
    end
        # while zero distances exist (repeat above)
    end
end

"""update Flood.distance when a stone is played.  this is done by
moving from each point towards the new stone.  if the path reaches the
stone without meeting any other, and the distance is less than the
current value, then it is updated.  i am not completely sure this is
exact, but it appears to work ok."""
function flood_to_point!{N}(f::Flood{N}, b::Board{N}, t::Point, x, y)
    @forall i j N begin
        if point(b, i, j) == empty
            ii, jj, d = i, j, 1
            while true
                di, dj, d = x-ii, y-jj, d+1
                if abs(di) >= abs(dj)
                    ii = ii + sign(di)
                else
                    jj = jj + sign(dj)
                end
                if point(b, ii, jj) != empty
                    break
                end
            end
            if ii == x && jj == y
                prev = f.distance[i, j]
                if prev == 0 || abs(d) < abs(prev)
                    f.distance[i, j] = d * Int(t)
                end
            end
        end
    end
end

"""re-calculate (from scratch) the number of lives associated with
each group."""
function calculate_lives!{N}(g::Groups{N})
    # the general idea here is that for each space, we see whether it
    # contributes a life to any group (by lookng at neighbouting
    # groups).  then we add those lives up.  this avoids "double
    # counting" spaces that touch more than one stone.
    for k in 1:255
        g.lives[k] = 0
    end
    # local vars
    temp = zeros(UInt8, N, N, 4)   # neighbouring groups
    n = zeros(UInt8, N, N)         # number of neighbouring groups (<=4)
    @forall i j N begin
        # for spaces
        if g.index[i, j] == 0
            @forneighbours i j N ii jj begin
                group = g.index[ii, jj]
                # if group is not already known
                if group != 0 && !any([group == temp[i, j, k] for k in 1:n[i, j]])
                    n[i, j] = n[i, j] + 1
                    temp[i, j, n[i, j]] = group
                end
            end
            # thread group mem with memory guard
            for k in 1:n[i, j]
                g.lives[temp[i, j, k]] = g.lives[temp[i, j, k]] + 1
            end
        end
    end
end

"""the Space.border array uses this mask to indicate the border
colour(s).  it is chosen so that updating adjacent values using or
works ok, but needs translating at the net level."""
function border_mask(t::Point)
    if t == empty
        return 0
    elseif t == black
        return 1
    else
        return 2
    end
end

"""iteratively update space wherever the border is (initially) zero
until the index and border are consistent acrosss the group (lowest
index, combined border)."""
function fix_space!{N}(s::Space{N}, b::Board{N})

    my_active = zeros(Bool, N, N)
    my_nsew = zeros(UInt8, N, N)   # bitmask of stones (to ignore in round 2)
    my_index = zeros(UInt8, N, N)
    my_border = zeros(UInt8, N, N)
    my_k = zeros(UInt8, N, N)

    # first round, set border from stones
    @forall i j N begin
        my_active[i, j] = s.border[i, j] == 0 && point(b, i, j) == empty
        if my_active[i, j]
            my_index[i, j] = s.index[i, j]
            @forneighbours i j N ii jj begin
                t = point(b, ii, jj)
                if t != empty
                    my_border[i, j] = my_border[i, j] | border_mask(t)
                    my_nsew[i, j] = my_nsew[i, j] | 2^my_k[i, j]
                end
                my_k[i, j] = my_k[i, j] + 1
            end
            s.border[i, j] = my_border[i, j]
        end
    end

    # second round - update until consistent
    our_changed = true
    while our_changed
        our_changed = false
        my_k = zeros(UInt8, N, N)

        @forall i j N begin
            my_changed = zeros(Bool, N, N)
            if my_active[i, j]
                @forneighbours i j N ii jj begin
                    if my_nsew[i, j] & 2^my_k[i, j] == 0
                        index = s.index[ii, jj]
                        if index != my_index[i, j] && index > 0 && (my_index[i, j] == 0 || index < my_index[i, j])
                            my_changed[i, j] = true
                            my_index[i, j] = index
                        end
                        border = s.border[ii, jj]
                        if !(border == 0 || border & my_border[i, j] == border)
                            my_changed[i, j] = true
                            my_border[i, j] = my_border[i, j] | border
                        end
                    end
                    my_k[i, j] = my_k[i, j] + 1
                end
                if my_changed[i, j]
                    s.index[i, j] = my_index[i, j]
                    s.border[i, j] = my_border[i, j]
                    our_changed = true
                end
            end
        end
    end

end

"""add space indices around a new move.  these are then resolved by
fix_space!()."""
function index_new_space!{N}(s::Space{N}, b::Board{N}, x, y)
    available = k_lowest_unused(4, s.index, N)
    k = 1
    @forneighbours x y N xx yy begin
        if point(b, xx, yy) == empty
            s.index[xx, yy] = available[k]
            k = k+1
        end
    end
end

function stats!{N}(z::Stats, s::Space{N}, x, y)
    function f(z, b)
        if b == 1
            z.colour[black].spaces = z.colour[black].spaces + 1
        elseif b == 2
            z.colour[white].spaces = z.colour[white].spaces + 1
        end
        z
    end
    z.colour[black].spaces, z.colour[white].spaces = 0, 0
    foldl(f, z, s.border)
    z.stones = count(x -> x == 0, s.border)
    z.prev = (x, y)
    stats!(z)
end

function stats!(z::Stats)
    b, w = z.colour[black], z.colour[white]
    z.score = b.prisoners + b.spaces - (w.prisoners + w.spaces) - komi
    z.owned = b.spaces + w.spaces
    z.moves += 1
end

type IllegalMove <: Exception end

function assert_empty(p, x, y)
    if point(p, x, y) != empty
        throw(IllegalMove())
    end
end

function assert_alive(p, x, y)
    if p.groups.lives[p.groups.index[x, y]] == 0
        throw(IllegalMove())
    end
end

"""update the row to include the given Point."""
function move(row::Row, t::Point, x)
    k = 3^(x-1)
    l, r = divrem(row, k)
    l = 3 * div(l, 3) + Row(Int(t)+1)
    l * k + r
end

"""update the board to include the given point."""
move!(b::Board, t::Point, x, y) = b.row[y] = move(b.row[y], t, x)

"""update Flood.distance to include the given point (only - full
flood update requires flood_dead_group! and flood_to_point!)."""
move!(f::Flood, t::Point, x, y) = f.distance[x, y] = Int8(t)

function move!{N}(s::Space{N}, x, y)
    index = s.index[x, y]
    @forall i j N begin
        if s.index[i, j] == index
            s.index[i, j] = 0
            s.border[i, j] = 0
        end
    end
end

"""update the position, given a new point.  this includes all
processing."""
function move!{N}(p::Position{N}, t::Point, x, y)
    assert_empty(p, x, y)
    move!(p.board, t, x, y)
    move!(p.flood, t, x, y)
    move!(p.space, x, y)
    flood_to_point!(p.flood, p.board, t, x, y)
    index_new_space!(p.space, p.board, x, y)
    newgroup = k_lowest_unused(1, p.groups.index, N)[1]
    p.groups.index[x, y] = newgroup
    p.groups.size[newgroup] = 1
    @forneighbours x y N xx yy begin
        tt = point(p.board, xx, yy)
        if tt == t
            merge_group!(p.groups, newgroup, xx, yy)
        elseif tt == other(t)
            g = p.groups.index[xx, yy]
            # if it had a single life, and we just filled a space next
            # door, then it must have died
            if p.groups.lives[g] == 1
                delete_group!(p, t, g)
            end
        end
    end
    flood_dead_group!(p.flood, p.board, t)
    calculate_lives!(p.groups)
    fix_space!(p.space, p.board)
    assert_alive(p, x, y)
    stats!(p.stats, p.space, x, y)
    p   # support call with new instance
end

function pass!(t, p::Position)
    p.stats.prev = Nullable{Tuple{Int,Int}}()
    p.stats.colour[other(t)].prisoners += 1
    stats!(p.stats)
    p
end

# --- query state

function valid{N}(h::Vector{Position{N}}, t::Point, x, y)
    if point(h[end], x, y) != empty
        return false
    end
    if ko(h, t, x, y)
        return false
    end
    @forneighbours x y N xx yy begin
        tt = point(h[end], xx, yy)
        if tt == empty
            return true
        elseif tt == t
            if h[end].groups.lives[h[end].groups.index[xx, yy]] > 1
                return true
            end
        else
            if h[end].groups.lives[h[end].groups.index[xx, yy]] == 1
                return true
            end
        end
    end
    false
end

# test for ko
function ko{N}(h::Vector{Position{N}}, t::Point, x, y)
    # we need at least two previous positions (actually more, but
    # perhaps we only store two)
    if length(h) < 2
        return false
    end
    # for a ko situation, the previous move must have captured a
    # single piece.
    opponent = other(t)
    if (h[end].stats.colour[opponent].prisoners -
        h[end-1].stats.colour[opponent].prisoners) != 1
        return false
    end
    if (h[end].stats.colour[opponent].spaces -
        h[end-1].stats.colour[opponent].spaces) != 1
        return false
    end
    # and there must be a possible group that we capture now, with a
    # single stone, that has just moved there.
    candidate = false
    @forneighbours x y N xx yy begin
        if point(h[end], xx, yy) == opponent
            g = h[end].groups.index[xx, yy]
            if h[end].groups.lives[g] == 1
                if h[end].groups.size[g] == 1 && point(h[end-1], xx, yy) == empty && !candidate
                    candidate = true
                else
                    # we killed a non-ko group, so this isn't ko
                    return false
                end
            end
        end
    end
    return candidate
end


text:
# This file was generated, do not modify it. # hide
combine(groupby(iris, :Species), :PetalLength => mean)


text:
promptname(level, name) = "$level|$name> "
function RunDebugger(frame, repl = Base.active_repl, terminal = Base.active_repl.t; initial_continue=false)
    state = DebuggerState(frame, repl, terminal)

    # Setup debug panel
    panel = LineEdit.Prompt(promptname(state.level, "debug");
        prompt_prefix = "\e[38;5;166m",
        prompt_suffix = Base.text_colors[:normal],
        on_enter = s->true)

    panel.hist = REPL.REPLHistoryProvider(Dict{Symbol,Any}(:debug => panel))
    REPL.history_reset_state(panel.hist)

    search_prompt, skeymap = LineEdit.setup_search_keymap(panel.hist)
    search_prompt.complete = REPL.LatexCompletions()

    state.main_mode = panel

    panel.on_done = (s,buf,ok)->begin
        line = String(take!(buf))
        old_level = state.level
        if !ok || strip(line) == "q"
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if isempty(strip(line)) && length(panel.hist.history) > 0
            command = panel.hist.history[end]
        else
            command = strip(line)
        end
        do_print_status = true
        cmd1 = split(command,' ')[1]
        do_print_status = try
            execute_command(state, Val{Symbol(cmd1)}(), command)
        catch err
            # No point showing an internal stacktrace
            Base.display_error(Base.pipe_writer(terminal), err, [])
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if old_level != state.level
            panel.prompt = promptname(state.level, "debug")
        end
        LineEdit.reset_state(s)
        if state.frame === nothing
            LineEdit.transition(s, :abort)
            LineEdit.reset_state(s)
            return false
        end
        if do_print_status
            print_status(Base.pipe_writer(terminal), active_frame(state))
        end
        return true
    end

    key = '`'
    repl_switch = Dict{Any,Any}(
        key => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                prompt = julia_prompt(state)
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, prompt) do
                    LineEdit.state(s, prompt).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s,key)
            end
        end
    )

    state.standard_keymap = Dict{Any,Any}[skeymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap([repl_switch;state.standard_keymap])

    if initial_continue
        execute_command(state, Val(:c), "c")
        state.frame === nothing && return state.overall_result
    end
    print_status(Base.pipe_writer(terminal), active_frame(state))
    REPL.run_interface(terminal, LineEdit.ModalInterface([panel,search_prompt]))

    return state.overall_result
end

function julia_prompt(state::DebuggerState)
    # Return early if this has already been called on the state
    isassigned(state.julia_prompt) && return state.julia_prompt[]

    julia_prompt = LineEdit.Prompt(promptname(state.level, "julia");
        # Copy colors from the prompt object
        prompt_prefix = state.repl.prompt_color,
        prompt_suffix = (state.repl.envcolors ? Base.input_color : state.repl.input_color),
        complete = DebugCompletionProvider(state),
        on_enter = REPL.return_callback)
    julia_prompt.hist = state.main_mode.hist
    julia_prompt.hist.mode_mapping[:julia] = julia_prompt

    julia_prompt.on_done = (s,buf,ok)->begin
        if !ok
            LineEdit.transition(s, :abort)
            return false
        end
        xbuf = copy(buf)
        command = String(take!(buf))
        @static if VERSION >= v"1.2.0-DEV.253"
            response = _eval_code(active_frame(state), command)
            REPL.print_response(state.repl, response, true, true)
        else
            ok, result = _eval_code(active_frame(state), command)
            REPL.print_response(state.repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
        end
        println(state.terminal)
        LineEdit.reset_state(s)
    end
    julia_prompt.keymap_dict = LineEdit.keymap([REPL.mode_keymap(state.main_mode); state.standard_keymap])
    state.julia_prompt[] = julia_prompt
    return julia_prompt
end

@static if VERSION >= v"1.2.0-DEV.253"
    function _eval_code(frame::Frame, code::AbstractString)
        try
            return eval_code(frame, code), false
        catch
            return Base.catch_stack(), true
        end
    end
else
    function _eval_code(frame::Frame, code::AbstractString)
        try
            return true, eval_code(frame, code)
        catch err
            return false, (err, catch_backtrace())
        end
    end
end

function eval_code(frame::Frame, command::AbstractString)
    expr = Base.parse_input_line(command)
    isexpr(expr, :toplevel) && (expr = expr.args[end])
    # see https://github.com/JuliaLang/julia/issues/31255 for the Symbol("") check
    vars = filter(v -> v.name != Symbol(""), JuliaInterpreter.locals(frame))
    res = gensym()
    eval_expr = Expr(:let,
        Expr(:block, map(x->Expr(:(=), x...), [(v.name, v.value) for v in vars])...),
        Expr(:block,
            Expr(:(=), res, expr),
            Expr(:tuple, res, Expr(:tuple, [v.name for v in vars]...))
        ))
    eval_res, res = Core.eval(moduleof(frame), eval_expr)
    j = 1
    for (i, v) in enumerate(vars)
        if v.isparam
            frame.framedata.sparams[j] = res[i]
            j += 1
        else
            frame.framedata.locals[frame.framedata.last_reference[v.name]] = Some{Any}(res[i])
        end
    end
    eval_res
end

# Completions

mutable struct DebugCompletionProvider <: REPL.CompletionProvider
    state::DebuggerState
end

function LineEdit.complete_line(c::DebugCompletionProvider, s)
    partial = REPL.beforecursor(s.input_buffer)
    full = LineEdit.input_string(s)
    ret, range, should_complete = completions(c, full, lastindex(partial))
    return unique!(map(REPLCompletions.completion_text, ret)), partial[range], should_complete
end

function completions(c::DebugCompletionProvider, full, partial)
    mod = moduleof(c.state.frame)
    ret, range, should_complete = REPLCompletions.completions(full, partial, mod)

    # TODO Add local variable completions
    return ret, range, should_complete
end