text:
"""
    parameter(sma, ecc)

Get the parameter (semi-latus rectum) of the orbit.
"""
parameter(sma, ecc) = sma * (1 - ecc * ecc)

"""
    periapsis_radius(sma, ecc)

Get the smallest radius along the orbit.
"""
periapsis_radius(sma, ecc) = sma * (1 - ecc)

"""
    apoapsis_radius(sma, ecc)

Get the largest radius along the orbit.
"""
apoapsis_radius(sma, ecc) = sma * (1 + ecc)

"""
    radius(sma, ecc, ta)

Compute the radius of the orbit at the specified true anomaly
"""
radius(sma, ecc, ta) = parameter(sma, ecc) / (1 + ecc * cos(ta))

"""
    true_longitude(raan, aop, ta)

Compute the true longitude
"""
true_longitude(raan, aop, ta) = raan + aop + ta

"""
    isclosed(ecc)

Return true if the two-body orbit is closed, false otherwise.
"""
isclosed(ecc) = ecc < 1

"""
    isopen(ecc)

Return true if the two-body orbit is open, false otherwise.
"""
isopen(ecc) = !isclosed(ecc)

"""
    iscircular(ecc)

Return true if the two-body orbit is circular (``e = 0``), false otherwise.
"""
iscircular(ecc) = ecc == 0

"""
    iselliptical(ecc)

Return true if the two-body orbit is elliptical (``e < 1``), false otherwise.
"""
iselliptical(ecc) = ecc < 1

"""
    isparabolic(ecc)

Return true if the two-body orbit is parabolic (``e = 1``), false otherwise.
"""
isparabolic(ecc) = ecc == 1

"""
    ishyperbolic(ecc)

Return true if the two-body orbit is hyperbolic (``e > 1``), false otherwise.
"""
ishyperbolic(ecc) = ecc > 1

# Defining this allows us to write f(x) rather than f(x) _and_ f(cb::CelestialBody)
_GM(x) = gravitational_parameter(x)
_GM(x::Real) = x

"""
    circular_velocity(gm, r)
    circular_velocity(cb::CelestialBody, r)

Compute the circular velocity of a two-body orbit with circular radius `r`.
"""
circular_velocity(gm, r) = sqrt(_GM(gm) / r)

"""
    mean_motion(gm, sma)
    mean_motion(cb::CelestialBody, sma)

Compute the mean motion of an orbit with semi-major axis, `sma`.
"""
mean_motion(gm, sma) = circular_velocity(_GM(gm), sma) / sma

"""
    period(gm, sma)
    period(cb::CelestialBody, sma)

Compute the period of an orbit with semi-major axis, `sma`.
"""
period(gm, sma) = 2π / mean_motion(_GM(gm), sma)

"""
    angular_momentum(gm, sma, ecc)
    angular_momentum(cb::CelestialBody, sma, ecc)

Compute the angular momentum of the two-body orbit.
"""
angular_momentum(gm, sma, ecc) = sqrt(_GM(gm) * parameter(sma, ecc))

"""
    flight_path_angle(sma, ecc, ta)

Compute the angle between the velocity vector and the local horizon.

A positive flight path angle corresponds to a negative rotation from the local horizon to
the velocity vector around the angular momentum.
"""
function flight_path_angle(sma, ecc, ta)
    r = radius(sma, ecc, ta)
    return acos(sqrt(sma^2 * (1 - ecc^2) / (r * (2sma - r))))
end

include("elements.jl")
include("keplerian.jl")


text:
using Oceananigans
using Oceanostics

grid = RegularRectilinearGrid(size=(4,4,4), extent=(1,1,1))
model = IncompressibleModel(grid=grid)
u, v, w = model.velocities

ke = Oceanostics.KineticEnergy(model, u, v, w)

SPx = Oceanostics.XShearProduction(model, u, v, w, 0, 0, 0)
SPy = Oceanostics.YShearProduction(model, u, v, w, 0, 0, 0)
SPz = Oceanostics.ZShearProduction(model, u, v, w, 0, 0, 0)

#ε_iso = Oceanostics.IsotropicViscousDissipation(model, ν, u, v, w)
#ε_ani = Oceanostics.AnisotropicViscousDissipation(model, ν, ν, ν, u, v, w)

ext: jl


text:
# Autogenerated wrapper script for Libtask_jll for x86_64-apple-darwin
export libtask_julia

JLLWrappers.@generate_wrapper_header("Libtask")
JLLWrappers.@declare_library_product(libtask_julia, "@rpath/libtask_julia.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtask_julia,
        "lib/libtask_julia.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SafeTestsets
@safetestset "MA Tests" begin include("MA_tests.jl") end
@safetestset "RI Tests" begin include("RI_tests.jl") end
@safetestset "SEG_LIST Tests" begin include("SEG_LIST_tests.jl") end


text:
#
# C45.jl - Decision Trees for Julia
#
# This file implements decision trees in a fairly straight-forward manner compliant with the
# C4.5 algorithm. The goal is to provide a framework to make decision trees with more than
# two child nodes in cases where there are discrete attributes, unlike the DecisionTree.jl
# package.
#
# This package is inspired by the DecisionTree.jl package, and Prof. Sean Luke's AI lectures
# during the 2013 Fall semester at GMU. This is actually a port of my original Common Lisp
# with extensions taken from the DecisionTree package.

module C45

using DataFrames

export DTLeaf, DTInternal, information, gini

# abstract type for internal node use
abstract DTNode

# Empty node class, for constructors
type EmptyNode <: DTNode
end

# Leaf node
type DTLeaf <: DTNode
	majority::Any
	values:Vector
end

# internal node
type DTInternal <: DTNode
	featid::Integer
	evl::Function # should return the index of the appropriate node given a value in the domain of the feature
	children::Vector{DTNode}

	DTInternal() = new(0,def_choice,Vector(DTNode,1))
end

# A wrapper around DataFrame to add labels to data.
# Note that labels can have NA as a value.
type LabeledData
	data::DataFrame
	labels::Array
end

function def_choice(array::Vector)
	array[rand(1:length(array))]
end

# getColumn - A function to...get a column from a matrix
# Inputs:
#	a - the3 array to get the folumn from
#	col - the column number, index starting at 1
#
# Output:
#	vec - A column vector from a at index col
#
function getColumn(a::Array, col::Integer)
	if 2 < dims(a) || dims(a) < 2
		throw(ArgumentError())
	end

	if col < 1
		throw(ArgumentError())
	end

	rows = size(a,1)
	vec = Array(eltype(a), rows, 1)

	for i in 1:rows
		vec[i] = a[i, col]
	end

	return vec
end

# log_2 - a function to handle the undefined case of log2(0)
# Inputs:
#	x - the input value
# Outputs:
#	Outputs either 0 if x == 0, or log2(x)
function log_2(x)
	if x == 0
		return x
	else
		return log2(x)
	end
end

# countAppearances - a function to see how many times something appears in something else
# Inputs:
#	a - the array object to search through
#	thing - what you are searching for
#
# Output:
#	count - the integer of the number of times thing appeared in a
function countAppearances(a::Array, thing)
	count = 0
	for e in a
		if a == thing
			count += 1
		end
	end
	return count
end

# information - the function to calculate the information content as per information theory
# Inputs:
#	probabilities - the data array, should be a column vector
# Outputs:
#	sum - the information contained in the probability distribution
function information(probabilities::Array)
	sum = 0.0
	for i in probabilities
		sum += (i * log_2(i))
	end
	return (sum * -1)
end

# gini - calculates the GINI coefficient for a list of probabilities
# Inputs:
#	probabilities - a list of probabilities, should sum to 1
# Outputs:
#	the gini coefficient.
function gini(probabilities::Array)
	sum = 0.0
	for p in probabilities
		sum += p^2
	end
	return (1 - sum)
end

# Note that the labels NEED to be in the same order.
#function genProbabilities(feature::Array, samples::Array, labels::Array, vals::Array=unique(samples))
#	data = getColumn(samples, feature[1])
#	if feature[2] == "discrete"
#		counts = map(x -> countAppearance(samples, x), vals)
		pos_counts = Array(eltype(samples), size(samples,1))
#		for s in 1:size(samples,1)

		end
#	else
#
#	end
#end

function chooseFeature(features::Array, samples::Array, evl::Function = information)
	best = Array(Any,2)
	best[2] = 0
	best[1] = 0
	for f in features
		probs = genProbabilities(f, samples)
		if evl(probs) > best[2]
			best[1] = f
			best[2]=evl(probs)
		end
	end
	return best[1]
end

function choosePivot(samples::Array, evl::Function=information,  labels::Array)
	end = size(samples,1) - 1
	pivots = Array(eltype(samples),end)
	for i in 1:end
		pivot[i] = (samples[i] + samples[i+1]) / 2
	end

	best = [NA,0]
	for p in pivots
		probs = genProbabilities(samples. labels)
		res = evl(probs)
		if res > best[2]
			best[1] = p
			best[2] = res
		end
	end
	return best[1]
end


text:
a = create_var(1, 1)
b = eye(3)
c = a * b
test_expr(c, (a,))

d = create_var(5, 3)
f = a * 3.1 + 5.3 * d
test_expr(f, (a, d))

g = rand(4, 5) * d
test_expr(g, (d,))

h = d * rand(3, 4)
test_expr(h, (d,))

j = d / 4
test_expr(j, (d,))

info("All multiply/divide tests passed")


text:
module MRC

using Dates, Statistics, Mmap, CodecZlib, CodecBzip2, CodecXz, TranscodingStreams

export MRCHeader, MRCExtendedHeader, MRCData
export cellangles,
    cellangles!,
    cellsize,
    cellsize!,
    data,
    eachmapcol,
    eachmaprow,
    eachmapsection,
    eachstackunit,
    extendedheader,
    gridsize,
    gridsize!,
    header,
    origin,
    origin!,
    read_mmap,
    start,
    start!,
    updateheader!,
    voxelaxes,
    voxelsize,
    voxelsize!

const COMPRESSIONS = (
    gz=(
        magic=b"\x1f\x8b\b",
        extension=".gz",
        compressor=GzipCompressorStream,
        decompressor=GzipDecompressorStream,
    ),
    bz2=(
        magic=b"BZh",
        extension=".bz2",
        compressor=Bzip2CompressorStream,
        decompressor=Bzip2DecompressorStream,
    ),
    xz=(
        magic=b"\xfd7zXZ\0",
        extension=".xz",
        compressor=XzCompressorStream,
        decompressor=XzDecompressorStream,
    ),
    none=(magic=b"", extension="", compressor=NoopStream, decompressor=NoopStream),
)
const COMPRESSOR_MAGICS = Dict(spec.magic => type for (type, spec) in pairs(COMPRESSIONS))
const COMPRESSOR_EXTENSIONS = Dict(
    spec.extension => type for (type, spec) in pairs(COMPRESSIONS)
)
@enum ByteOrder LittleEndian BigEndian
const MACHINE_STAMP_LITTLE = b"DD\0\0"
const MACHINE_STAMP_LITTLE_ALT = b"DA\0\0"
const MACHINE_STAMP_BIG = b"\x11\x11\0\0"
const HEADER_LENGTH = 1024
const MAP_NAME = "MAP "
const MODE_TO_TYPE = Dict(
    Int32(0) => Int8,
    Int32(1) => Int16,
    Int32(2) => Float32,
    Int32(3) => Complex{Int16},
    Int32(4) => ComplexF32,
    Int32(6) => UInt16,
)
const TYPE_TO_MODE = Dict(reverse(p) for p in MODE_TO_TYPE)
const MRC2014_VERSION = Int32(20140)
@enum MRCDataType Image ImageStack Volume VolumeStack Unknown

include("utils.jl")
include("header.jl")
include("extended_header.jl")
include("data.jl")
include("io.jl")

end # module


text:
export
    SimulationOptions,
    build_environment,
    build_uavs,
    get_joint_actions,
    load_policy,
    ACTIONS

#=
Functions for building the environment from options and default values.
=#

### Simulation settings and constants ###

# These constants do not go in the options class because they
# should essentially never change
const RANGES = [0, 25, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 900, 1000, 1200, 1400, 1600, 1800, 2000, 2500, 3000, 4000, 6000, 8000, 10000, 12800]
const ACTIONS = deg2rad([-12, -6, 0, 6, 12, rad2deg(-1)])
const THETA_DIM = 37
const THETAS = collect(linspace(-pi, pi, THETA_DIM))
const BEARING_DIM = 37
const BEARINGS = collect(linspace(0, 2 * pi, BEARING_DIM))
const SPEED_DIM = 3
const SPEEDS = collect(linspace(35, 45, SPEED_DIM))
const RESPONDINGS = [0, 1]

### Helper functions for building the environment ###

# type for default simulation settings
type SimulationOptions
    # coordination table policy
    coordination_table_one_policy::Bool
    coordination_policy_filepath::String

    # ownship options
    ownship_policy_filepath::String
    is_ownship_policy_joint::Bool
    policy_name::String
    ownship_strategy::Function # {greedy_strategy, coordinated_strategy}
    ownship_action_selection::String # {best*, worst*, average*}
    ownship_range_sensor_std::Float64
    ownship_theta_sensor_std::Float64
    ownship_psi_sensor_std::Float64
    ownship_velocity_sensor_std::Float64
    use_ownship_sensor_noise::Bool

    # intruder options
    intruder_policy_filepath::String
    is_intruder_policy_joint::Bool
    intruder_strategy::Function # {greedy_strategy, coordinated_strategy}
    intruder_action_selection::String # {best*, worst*, average*}
    intruder_range_sensor_std::Float64
    intruder_theta_sensor_std::Float64
    intruder_psi_sensor_std::Float64
    intruder_velocity_sensor_std::Float64
    use_intruder_sensor_noise::Bool

    # environment options
    coordination_message_scheme::Int64
    start_state_generator::Function
    use_noise::Bool
    random_seed::Int64
    gauss_sigma::Float64
    gauss_mu::Float64
    max_steps::Int64
    min_start_range_ratio::Float64
    max_start_range_ratio::Float64
    range_std::Float64
    initial_heading_offset::Float64
    initial_min_separation::Float64
    delayed_pilot_response::Bool

    # dynamics options
    omega::Float64
    dt::Int64
    dt_step::Float64
    pid_b_std::Float64
    pid_b_std_coc::Float64

    function SimulationOptions()
        # optimal policy used for coordination table
        coordination_table_one_policy = true ## currently this is whether or not solve real time...
        table_filepath = joinpath(Pkg.dir("Coord"), "data", "qvalue_tables", "optimal.jld")
        coordination_policy_filepath = table_filepath

        # ownship options
        ownship_policy_filepath = table_filepath
        is_ownship_policy_joint = true # {false = non-joint, true = joint}

        policy_name = "solQ"
        ownship_strategy = greedy_strategy # {greedy_strategy, coordinated_strategy}
        ownship_action_selection = "average" # {best, worst, average}

        # intruder options
        intruder_policy_filepath = table_filepath
        is_intruder_policy_joint = true # {false = non-joint, true = joint}

        intruder_strategy = greedy_strategy
        intruder_action_selection = "average" # {best, worst, average}
        if ownship_strategy == greedy_strategy
            intruder_strategy = greedy_strategy
        end

        # environment options
        # coordination messaging scheme selection
        coordination_message_scheme = 0 # {0, 1, 2, 3, 4}

        # generators: {start_towards_origin_encounter, start_head_on_encounter}
        start_state_generator = start_towards_origin_encounter
        use_noise = true

        # for running on sherlock with multiple subtasks
        # get the id of this subtask to use as the seed
        random_seed = 1

        gauss_sigma = deg2rad(2)
        gauss_mu = 0
        max_steps = 20
        min_start_range_ratio = 5 / 64 # gives min range of 2000
        max_start_range_ratio = 5 / 32 # gives max range of 2000
        range_std = 100.
        initial_heading_offset = pi / 32
        initial_min_separation = 153.
        delayed_pilot_response = false

        # ownship and intruder sensor noise defaults
        ownship_range_sensor_std = 15. # meters
        ownship_theta_sensor_std = deg2rad(1.) # rad
        ownship_psi_sensor_std = deg2rad(1.) # rad
        ownship_velocity_sensor_std = 1. # m/s
        use_ownship_sensor_noise = true
        intruder_range_sensor_std = 15. # meters
        intruder_theta_sensor_std = deg2rad(1.) # rad
        intruder_psi_sensor_std = deg2rad(1.) # rad
        intruder_velocity_sensor_std = 1. # m/s
        use_intruder_sensor_noise = true

        # dynamics options
        omega = .2
        dt = 5
        dt_step = 1.
        pid_b_std = deg2rad(3)
        pid_b_std_coc = deg2rad(5)

        # build the default options
        return new(coordination_table_one_policy, coordination_policy_filepath, ownship_policy_filepath,
            is_ownship_policy_joint, policy_name, ownship_strategy, ownship_action_selection,
            ownship_range_sensor_std, ownship_theta_sensor_std, ownship_psi_sensor_std,
            ownship_velocity_sensor_std, use_ownship_sensor_noise, intruder_policy_filepath,
            is_intruder_policy_joint, intruder_strategy, intruder_action_selection,
            intruder_range_sensor_std, intruder_theta_sensor_std, intruder_psi_sensor_std,
            intruder_velocity_sensor_std, use_intruder_sensor_noise, coordination_message_scheme,
            start_state_generator, use_noise, random_seed, gauss_sigma, gauss_mu, max_steps,
            min_start_range_ratio, max_start_range_ratio, range_std, initial_heading_offset,
            initial_min_separation, delayed_pilot_response, omega, dt, dt_step, pid_b_std, pid_b_std_coc)
    end
end

#=
Description:
Returns an enviornment object

Return Value:
- environment: an environment object
=#
function build_environment()
    # should instead gather options from cmd line
    options = SimulationOptions()

    # build uavs and dynamics
    uavs = build_uavs(options)
    dynamics = build_dynamics(options)

    # build actual environment
    max_steps = options.max_steps
    use_noise = options.use_noise
    gauss_sigma = options.gauss_sigma
    gauss_mu = options.gauss_mu
    coordination_message_scheme = options.coordination_message_scheme
    start_state_generator = options.start_state_generator
    min_start_range_ratio = options.min_start_range_ratio
    max_start_range_ratio = options.max_start_range_ratio
    range_std = options.range_std
    initial_heading_offset = options.initial_heading_offset
    initial_min_separation = options.initial_min_separation
    delayed_pilot_response = options.delayed_pilot_response

    rng = MersenneTwister(options.random_seed)
    state_generator_rng = MersenneTwister(options.random_seed)
    e = Environment(max_steps, use_noise, uavs, gauss_sigma,
                    gauss_mu, dynamics, coordination_message_scheme, start_state_generator,
                    min_start_range_ratio, max_start_range_ratio, range_std,
                    initial_heading_offset, initial_min_separation,
                    delayed_pilot_response, rng, state_generator_rng)
    return e
end

#=
Description:
Make joint actions given list of actions

Parameters:
- actions: list of float valued actions

Return Value:
- array of actions such that each row corresponds to a uav
    and col to a joint action index, with the cell value being the action
=#
function get_joint_actions(actions)
    num_actions = length(actions)
    joint_actions = zeros(2, num_actions^2)
    i_actions = 1
    for ia2 = 1:num_actions
        for ia1 = 1:num_actions
            joint_actions[1, i_actions] = actions[ia1]
            joint_actions[2, i_actions] = actions[ia2]
            i_actions += 1
        end
    end
    return joint_actions
end

#=
Description:
Load in a policy from a filepath

Parameters:
- filepath: filepath to policy
- options: options type with info to build policy

Return Value:
- policy_table type object
=#
function load_policy(filepath::String, is_policy_joint::Bool, options::SimulationOptions)
    if is_policy_joint # joint policy
        q_table = load(filepath, options.policy_name)
        q_grid = RectangleGrid(RANGES, THETAS, BEARINGS, SPEEDS, SPEEDS, RESPONDINGS, RESPONDINGS)
        num_q_table_states = reduce(*, q_grid.cut_counts) + 1
        msg = "the q_table (rows = $(size(q_table, 1))) does not match the
            q_grid (length = $num_q_table_states), which was created using constant values"
        @assert num_q_table_states == size(q_table, 1) msg
        policy = PolicyTable(q_table, q_grid, num_q_table_states)
        return policy
    else # non joint policy
        q_table = load(filepath, options.policy_name)
        q_grid = RectangleGrid(RANGES, THETAS, BEARINGS, SPEEDS, SPEEDS, RESPONDINGS)
        num_q_table_states = reduce(*, q_grid.cut_counts) + 1
        msg = "the q_table (rows = $(size(q_table, 1))) does not match the
            q_grid (length = $num_q_table_states), which was created using constant values"
        @assert num_q_table_states == size(q_table, 1) msg
        policy = PolicyTable(q_table, q_grid, num_q_table_states)
        return policy
    end
end

#=
Description:
Returns a list of uav objects

Return Value:
- uavs: list of uav objects
=#
function build_uavs(options::SimulationOptions)
    # rng - this is for the sensor noise
    m = MersenneTwister(options.random_seed)

    # alias this constant value
    actions = ACTIONS

    # load ownship policy
    policy = load_policy(options.ownship_policy_filepath, options.is_ownship_policy_joint, options)

    # determine the joint actions
    joint_actions = get_joint_actions(actions)

    # build the coordination table
    if options.coordination_table_one_policy
        optimal_policy = load(options.coordination_policy_filepath, options.policy_name)
        c_table = make_coordination_table(optimal_policy, joint_actions)
        c_grid = RectangleGrid(RANGES, THETAS, BEARINGS, SPEEDS, SPEEDS, RESPONDINGS, RESPONDINGS)
        num_coord_states = reduce(*, c_grid.cut_counts)
        coord = CoordinationTable(c_table, c_grid, num_coord_states)
    else
        c_table = load(options.coordination_policy_filepath, options.policy_name)
        c_grid = RectangleGrid(RANGES, THETAS, BEARINGS, SPEEDS, SPEEDS, RESPONDINGS, RESPONDINGS)
        num_coord_states = reduce(*, c_grid.cut_counts)
        coord = CoordinationTable(c_table, c_grid, num_coord_states)
    end

 # create the ownship
    ownship_id = 1
    ownship_uav = UAV(ownship_id, actions, joint_actions, policy, coord, options.is_ownship_policy_joint,
                    options.ownship_strategy, options.ownship_action_selection, m,
                    options.ownship_range_sensor_std, options.ownship_theta_sensor_std,
                    options.ownship_psi_sensor_std, options.ownship_velocity_sensor_std,
                    options.use_ownship_sensor_noise)

    # create the intruder
    if options.intruder_policy_filepath == options.ownship_policy_filepath
        intruder_policy = policy
    else
        intruder_policy = load_policy(options.intruder_policy_filepath, options.is_intruder_policy_joint, options)
    end

    intruder_id = 2
    strategy = options.intruder_strategy
    intruder_uav = UAV(intruder_id, actions, joint_actions, intruder_policy, coord,
                    options.is_intruder_policy_joint, strategy, options.intruder_action_selection,
                    m, options.intruder_range_sensor_std, options.intruder_theta_sensor_std,
                    options.intruder_psi_sensor_std, options.intruder_velocity_sensor_std,
                    options.use_intruder_sensor_noise)

    return [ownship_uav, intruder_uav]
end

#=
Description:
Returns a dynamics object

Return Value:
- dynamics: a dynamics object
=#
function build_dynamics(options::SimulationOptions)
    m = MersenneTwister(options.random_seed)
    omega = options.omega
    use_noise = options.use_noise
    dt = options.dt
    dt_step = options.dt_step
    pid_b_std = options.pid_b_std
    pid_b_std_coc = options.pid_b_std_coc
    d = Dynamics(dt, dt_step, omega, pid_b_std, pid_b_std_coc, use_noise, m)
    return d
end


text:
using Documenter
using DecisionTrees

makedocs(
    modules = [DecisionTrees],
    sitename = "DecisionTrees.jl",
    pages = [
        "Home" => "index.md",
    ],
)

deploydocs(repo = "github.com/laschuet/DecisionTrees.jl.git")


text:
@testset "Automorphisms" begin

   G = PermutationGroup(Int8(4))

   @testset "AutSymbol" begin
      @test_throws MethodError Groups.AutSymbol(:a)
      @test_throws MethodError Groups.AutSymbol(:a, 1)
      f = Groups.AutSymbol(:a, 1, Groups.FlipAut(2))
      @test isa(f, Groups.GSymbol)
      @test isa(f, Groups.AutSymbol)
      @test isa(Groups.perm_autsymbol(Int8.([1,2,3,4])), Groups.AutSymbol)
      @test isa(Groups.rmul_autsymbol(1,2), Groups.AutSymbol)
      @test isa(Groups.lmul_autsymbol(3,4), Groups.AutSymbol)
      @test isa(Groups.flip_autsymbol(3), Groups.AutSymbol)
   end

   a,b,c,d = gens(FreeGroup(4))
   D = NTuple{4,FreeGroupElem}([a,b,c,d])

   @testset "flip_autsymbol correctness" begin
      @test Groups.flip_autsymbol(1)(deepcopy(D)) == (a^-1, b,c,d)
      @test Groups.flip_autsymbol(2)(deepcopy(D)) == (a, b^-1,c,d)
      @test Groups.flip_autsymbol(3)(deepcopy(D)) == (a, b,c^-1,d)
      @test Groups.flip_autsymbol(4)(deepcopy(D)) == (a, b,c,d^-1)
      @test inv(Groups.flip_autsymbol(1))(deepcopy(D)) == (a^-1, b,c,d)
      @test inv(Groups.flip_autsymbol(2))(deepcopy(D)) == (a, b^-1,c,d)
      @test inv(Groups.flip_autsymbol(3))(deepcopy(D)) == (a, b,c^-1,d)
      @test inv(Groups.flip_autsymbol(4))(deepcopy(D)) == (a, b,c,d^-1)
   end

   @testset "perm_autsymbol correctness" begin
      σ = Groups.perm_autsymbol([1,2,3,4])
      @test      σ(deepcopy(D)) == deepcopy(D)
      @test inv(σ)(deepcopy(D)) == deepcopy(D)

      σ = Groups.perm_autsymbol([2,3,4,1])
      @test      σ(deepcopy(D)) == (b, c, d, a)
      @test inv(σ)(deepcopy(D)) == (d, a, b, c)

      σ = Groups.perm_autsymbol([2,1,4,3])
      @test      σ(deepcopy(D)) == (b, a, d, c)
      @test inv(σ)(deepcopy(D)) == (b, a, d, c)

      σ = Groups.perm_autsymbol([2,3,1,4])
      @test      σ(deepcopy(D)) == (b, c, a, d)
      @test inv(σ)(deepcopy(D)) == (c, a, b, d)
   end

   @testset "rmul/lmul_autsymbol correctness" begin
      i,j = 1,2
      r = Groups.rmul_autsymbol(i,j)
      l = Groups.lmul_autsymbol(i,j)
      @test      r(deepcopy(D)) == (a*b,   b, c, d)
      @test inv(r)(deepcopy(D)) == (a*b^-1,b, c, d)
      @test      l(deepcopy(D)) == (b*a,   b, c, d)
      @test inv(l)(deepcopy(D)) == (b^-1*a,b, c, d)

      i,j = 3,1
      r = Groups.rmul_autsymbol(i,j)
      l = Groups.lmul_autsymbol(i,j)
      @test      r(deepcopy(D)) == (a, b, c*a,   d)
      @test inv(r)(deepcopy(D)) == (a, b, c*a^-1,d)
      @test      l(deepcopy(D)) == (a, b, a*c,   d)
      @test inv(l)(deepcopy(D)) == (a, b, a^-1*c,d)

      i,j = 4,3
      r = Groups.rmul_autsymbol(i,j)
      l = Groups.lmul_autsymbol(i,j)
      @test      r(deepcopy(D)) == (a, b, c, d*c)
      @test inv(r)(deepcopy(D)) == (a, b, c, d*c^-1)
      @test      l(deepcopy(D)) == (a, b, c, c*d)
      @test inv(l)(deepcopy(D)) == (a, b, c, c^-1*d)

      i,j = 2,4
      r = Groups.rmul_autsymbol(i,j)
      l = Groups.lmul_autsymbol(i,j)
      @test      r(deepcopy(D)) == (a, b*d,   c, d)
      @test inv(r)(deepcopy(D)) == (a, b*d^-1,c, d)
      @test      l(deepcopy(D)) == (a, d*b,   c, d)
      @test inv(l)(deepcopy(D)) == (a, d^-1*b,c, d)
   end

   @testset "AutGroup/Automorphism constructors" begin

      f = Groups.AutSymbol(:a, 1, Groups.FlipAut(1))
      @test isa(Automorphism{3}(f), Groups.GWord)
      @test isa(Automorphism{3}(f), Automorphism)
      @test isa(AutGroup(FreeGroup(3)), AbstractAlgebra.Group)
      @test isa(AutGroup(FreeGroup(1)), Groups.AbstractFPGroup)
      A = AutGroup(FreeGroup(1))
      @test isa(Groups.gens(A), Vector{Automorphism{1}})
      @test length(Groups.gens(A)) == 1
      A = AutGroup(FreeGroup(1), special=true)
      @test length(Groups.gens(A)) == 0
      A = AutGroup(FreeGroup(2))
      @test length(Groups.gens(A)) == 7
      gens = Groups.gens(A)

      @test isa(A(Groups.rmul_autsymbol(1,2)), Automorphism)
      @test A(Groups.rmul_autsymbol(1,2)) in gens

      @test isa(A(Groups.rmul_autsymbol(2,1)), Automorphism)
      @test A(Groups.rmul_autsymbol(2,1)) in gens

      @test isa(A(Groups.lmul_autsymbol(1,2)), Automorphism)
      @test A(Groups.lmul_autsymbol(1,2)) in gens

      @test isa(A(Groups.lmul_autsymbol(2,1)), Automorphism)
      @test A(Groups.lmul_autsymbol(2,1)) in gens

      @test isa(A(Groups.flip_autsymbol(1)), Automorphism)
      @test A(Groups.flip_autsymbol(1)) in gens

      @test isa(A(Groups.flip_autsymbol(2)), Automorphism)
      @test A(Groups.flip_autsymbol(2)) in gens

      @test isa(A(Groups.perm_autsymbol([2,1])), Automorphism)
      @test A(Groups.perm_autsymbol([2,1])) in gens
   end

   A = AutGroup(FreeGroup(4))

   @testset "eltary functions" begin

      f = Groups.perm_autsymbol([2,3,4,1])
      @test (Groups.change_pow(f, 2)).pow == 1
      @test (Groups.change_pow(f, -2)).pow == 1
      @test (inv(f)).pow == 1

      f = Groups.perm_autsymbol([2,1,4,3])
      @test isa(inv(f), Groups.AutSymbol)

      @test_throws MethodError f*f

      @test A(f)^-1 == A(inv(f))
   end

   @testset "reductions/arithmetic" begin
      f = Groups.perm_autsymbol([2,3,4,1])

      f² = Groups.r_multiply(A(f), [f], reduced=false)
      @test Groups.simplifyperms!(f²) == false
      @test f²^2 == A()

      a = A(Groups.rmul_autsymbol(1,2))*Groups.flip_autsymbol(2)
      b = Groups.flip_autsymbol(2)*A(inv(Groups.rmul_autsymbol(1,2)))
      @test a*b == b*a
      @test a^3 * b^3 == A()
      g,h = Groups.gens(A)[[1,8]] # (g, h) = (ϱ₁₂, ϱ₃₂)

      @test Groups.domain(A) == NTuple{4, FreeGroupElem}(gens(A.objectGroup))

      @test (g*h)(Groups.domain(A)) == (h*g)(Groups.domain(A))
      @test (g*h).savedhash == zero(UInt)
      @test (h*g).savedhash == zero(UInt)
      a = g*h
      b = h*g
      @test hash(a) != zero(UInt)
      @test hash(b) == hash(a)
      @test a.savedhash == b.savedhash
      @test length(unique([a,b])) == 1
      @test length(unique([g*h, h*g])) == 1

      # Not so simple arithmetic: applying starting on the left:
      # ϱ₁₂*ϱ₂₁⁻¹*λ₁₂*ε₂ == σ₂₁₃₄

      g = A(Groups.rmul_autsymbol(1,2))
      x1, x2, x3, x4 = Groups.domain(A)
      @test g(Groups.domain(A)) == (x1*x2, x2, x3, x4)
      g = g*inv(A(Groups.rmul_autsymbol(2,1)))
      @test g(Groups.domain(A)) == (x1*x2, x1^-1, x3, x4)
      g = g*A(Groups.lmul_autsymbol(1,2))
      @test g(Groups.domain(A)) == (x2, x1^-1, x3, x4)
      g = g*A(Groups.flip_autsymbol(2))
      @test g(Groups.domain(A)) == (x2, x1, x3, x4)

      @test g(Groups.domain(A)) == A(Groups.perm_autsymbol([2,1,3,4]))(Groups.domain(A))

      @test g == A(Groups.perm_autsymbol([2,1,3,4]))

      g_im = g(Groups.domain(A))
      @test length(g_im[1]) == 5
      @test length(g_im[2]) == 3
      @test length(g_im[3]) == 1
      @test length(g_im[4]) == 1
      @test length.(Groups.reduce!.(g_im)) == (1,1,1,1)

   end

   @testset "specific Aut(F4) tests" begin
      N = 4
      G = AutGroup(FreeGroup(N))
      S = G.gens
      @test isa(S, Vector{Groups.AutSymbol})
      S = [G(s) for s in unique(S)]
      @test isa(S, Vector{Automorphism{N}})
      @test S == gens(G)
      @test length(S) == 51
      S_inv = [S..., [inv(s) for s in S]...]
      @test length(unique(S_inv)) == 75

      G = AutGroup(FreeGroup(N), special=true)
      S = gens(G)
      S_inv = [G(), S..., [inv(s) for s in S]...]
      S_inv = unique(S_inv)
      B_2 = [i*j for (i,j) in Base.product(S_inv, S_inv)]
      @test length(B_2) == 2401
      @test length(unique(B_2)) == 1777
   end

   @testset "linear_repr tests" begin
      N = 3
      G = AutGroup(FreeGroup(N))
      S = unique([gens(G); inv.(gens(G))])
      R = 3

      @test Groups.linear_repr(G()) isa Matrix{Int}
      @test Groups.linear_repr(G()) == Matrix{Int}(I, N, N)

      M = Matrix{Int}(I, N, N)
      M[1,2] = 1
      ϱ₁₂ = G(Groups.rmul_autsymbol(1,2))
      λ₁₂ = G(Groups.rmul_autsymbol(1,2))

      @test Groups.linear_repr(ϱ₁₂) == M
      @test Groups.linear_repr(λ₁₂) == M

      M[1,2] = -1

      @test Groups.linear_repr(ϱ₁₂^-1) == M
      @test Groups.linear_repr(λ₁₂^-1) == M

      @test Groups.linear_repr(ϱ₁₂*λ₁₂^-1) == Matrix{Int}(I, N, N)
      @test Groups.linear_repr(λ₁₂^-1*ϱ₁₂) == Matrix{Int}(I, N, N)

      M = Matrix{Int}(I, N, N)
      M[2,2] = -1
      ε₂ = G(Groups.flip_autsymbol(2))

      @test Groups.linear_repr(ε₂) == M
      @test Groups.linear_repr(ε₂^2) == Matrix{Int}(I, N, N)

      M = [0 1 0; 0 0 1; 1 0 0]

      σ = G(Groups.perm_autsymbol([2,3,1]))
      @test Groups.linear_repr(σ) == M
      @test Groups.linear_repr(σ^3) == Matrix{Int}(I, 3, 3)
      @test Groups.linear_repr(σ)^3 == Matrix{Int}(I, 3, 3)

      function test_homomorphism(S, r)
         for elts in Iterators.product([[g for g in S] for _ in 1:r]...)
            prod(Groups.linear_repr.(elts)) == Groups.linear_repr(prod(elts)) || error("linear representaton test failed at $elts")
         end
         return 0
      end

      @test test_homomorphism(S, R) == 0
   end
end


text:
export UnKnownRegionError, RegionIterator

struct UnKnownRegionError{R <: AbstractRegion} <: Exception
end

UnKnownRegionError(::Type{R}) where R = UnKnownRegionError{R}()

import Base: showerror

function showerror(io::IO, e::UnKnownRegionError{R}) where R
    print(io, "unknown region $R")
end

mutable struct RegionIterator{
        R <: AbstractRegion,
        B,
        LI<:LatticeIterator
    } <: AbstractHamiltonianIterator

    hamilton::LocalHamiltonian{B, R}
    lattice_it::LI
    RHS::Sites
end

RegionIterator(
    h::LocalHamiltonian{B, Region{0}},
    lattice::AbstractLattice,
    rhs::AbstractSites
) where B = RegionIterator(h, sites(lattice), rhs)

RegionIterator(
    h::LocalHamiltonian{B, Region{K}},
    lattice::AbstractLattice,
    rhs::AbstractSites
) where {B, K} = RegionIterator(h, bonds(lattice, K), rhs)

mutable struct RegionIterState{
        LST, Lv,
        HIT <: AbstractHamiltonianIterator,
        HST
    }

    lattice_st::LST # Lattice State
    lattice_va::Lv # Lattice Value
    hamilton_it::HIT # Hamilton Iter
    hamilton_st::HST # Hamilton State
end

function RegionIterState(itr::RegionIterator{Region{0}})
    lattice_st = start(itr.lattice_it)
    i, lattice_st = next(itr.lattice_it, lattice_st)
    hamilton_it = LHIterator(itr.hamilton, itr.RHS[i])
    hamilton_st = start(hamilton_it)
    return RegionIterState(lattice_st, i, hamilton_it, hamilton_st)
end

function RegionIterState(itr::RegionIterator{Region{K}}) where K
    lattice_st = start(itr.lattice_it)
    (i, j), lattice_st = next(itr.lattice_it, lattice_st)
    hamilton_it = LHIterator(itr.hamilton, itr.RHS[i], itr.RHS[j])
    hamilton_st = start(hamilton_it)
    return RegionIterState(lattice_st, (i, j), hamilton_it, hamilton_st)
end

import Base: start, next, done

# eltype(itr::RegionIterator) = Tuple{}
start(itr::RegionIterator) = RegionIterState(itr)

function next(itr::RegionIterator{Region{K}}, state) where K
    (val, (Si, Sj)), state.hamilton_st = next(state.hamilton_it, state.hamilton_st)
    i, j = state.lattice_va
    LHS = copy(itr.RHS)
    LHS[i] = Si; LHS[j] = Sj
    return (val, LHS), state
end

function done(itr::RegionIterator{Region{K}}, state) where K
    if done(state.hamilton_it, state.hamilton_st)
        if done(itr.lattice_it, state.lattice_st)
            return true
        end

        (i, j), state.lattice_st = next(itr.lattice_it, state.lattice_st)
        state.hamilton_it = LHIterator(itr.hamilton, itr.RHS[i], itr.RHS[j])
        state.hamilton_st = start(state.hamilton_it)
        state.lattice_va = (i, j)
    end
    return false
end

function next(itr::RegionIterator{Region{0}}, state)
    (val, Si), state.hamilton_st = next(state.hamilton_it, state.hamilton_st)
    LHS = copy(itr.RHS)
    LHS[i] = Si
    return (val, LHS), state
end

function done(itr::RegionIterator{Region{0}}, state)
    if done(state.hamilton_it, state.hamilton_st)
        if done(itr.lattice_it, state.lattice)
            return true
        end

        i, state.lattice_st = next(itr.lattice_it, state.lattice_st)
        state.hamilton_it = LHIterator(itr.hamilton, itr.RHS[i])
        state.hamilton_st = start(state.hamilton_it)
        state.lattice_va = i
    end
end


text:
using ModelingToolkit, OrdinaryDiffEq

@parameters t
@connector function HeatPort(;name)
    @variables Q_flow(t) T(t)
    ODESystem(Equation[], t, [Q_flow, T], [], name=name, defaults=[Q_flow=>0.0, T=>273.15]) # how to we update these guess values in the future/require them to be updated
end

# add in additional heat ports for flow and state - can there be some auto reminder of what to connect?

function ModelingToolkit.connect(::Type{HeatPort}, ports...)
    eqs =
    [
        0 ~ sum(port->port.Q_flow, ports)
    ]

    for i in 1:length(ports)-1
        push!(eqs, ports[i].T ~ ports[i+1].T)
    end

    return eqs
end


text:
# Date: July 1, 2018
# Author: Bernhard König
# Title: Modelling the residual of a GLM with a boosted tree

###########

t0 = time_ns()
cd(joinpath(GLOBAL_julia_code_folder, "DecisionTrees.jl"))
# @info("You may want to run 'pkg> instantiate' when you first run this. Use ] to enter the package mode.")
using Distributed
# Distributed.@everywhere using Revise
Distributed.@everywhere import CSV
Distributed.@everywhere import DataFrames
Distributed.@everywhere import DataFrames: DataFrame

Distributed.@everywhere using DecisionTrees

tela = (time_ns() - t0) / 1e9
@show tela # precompilation (first time usage) may take some time

##############################
# Read the data
##############################
datafile = joinpath("data", "freMTPL2", "freMTPL2.csv")
@assert isfile(datafile);
@time fullData = CSV.read(datafile, DataFrame, strict=true);

# add AreaInteger as new variable, i.e. A:F -> 1:6 (as suggested in the paper of Wüthrich, Noll, Salzmann)
areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

##############################
# Some Checks
##############################

# check poisson error versus errors reported in R
# this is a good check to ensure that the data is the same (train/test split, fitted GLM values are the same as in R)

glmFit = fullData[:GLMfitted]
errsGLM = poissonError(fullData[:ClaimNb], fullData[:GLMfitted])

inSampleErrorR = 31.26737970992647  # see R Code
outOfSampleErrorR = 32.171234111851014 # see R Code
valBitArrayIdx = fullData[:trnTest] .== 0;
trnBitArrayIdx = fullData[:trnTest] .== 1;
trnAvgErrorGLM = sum(errsGLM[trnBitArrayIdx]) / sum(trnBitArrayIdx)
valAvgErrorGLM = sum(errsGLM[valBitArrayIdx]) / sum(valBitArrayIdx)
trnAvgErrorGLM - inSampleErrorR / 100 # should be zero
valAvgErrorGLM - outOfSampleErrorR / 100 # should be zero

# consider average model
_avgfit = sum(fullData[:ClaimNb][trnBitArrayIdx]) / sum(fullData[:Exposure][trnBitArrayIdx])
avgfit = ones(size(fullData, 1)) .* _avgfit .* fullData[:Exposure]
errsAVGFit = poissonError(fullData[:ClaimNb], avgfit)
trnAvgErrAVGFit = sum(errsAVGFit[trnBitArrayIdx]) / sum(trnBitArrayIdx)
valAvgErrAVGFit = sum(errsAVGFit[valBitArrayIdx]) / sum(valBitArrayIdx)

# Just to get a 'feeling' of the size of the error: the GLM is 5% better (relative) than the average model.
trnAvgErrAVGFit - trnAvgErrorGLM # about 0.017
valAvgErrAVGFit - valAvgErrorGLM # about 0.017
1 / trnAvgErrAVGFit * trnAvgErrorGLM
1 / valAvgErrAVGFit * valAvgErrorGLM

originalTrnValIndex = deepcopy(fullData[:trnTest])

##############################
# Prepare the data
##############################

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]

# Let us now model the RESIDUAL of the GLM.
# thus we model observed/glmfitted, if we can model that ratio successfully, our model will improve.
# we have numerator=ClaimNb, denominator=GLMfitted
# thus a residual of 1.5 would mean that the GLM fit is 50% too low
# dfprepped is an intermediary dataframe which is generally not needed
dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, keycol="IDpol", trnvalcol="trnTest", numcol="ClaimNb", denomcol="GLMfitted", weightcol="Exposure", independent_vars=selected_explanatory_vars);

##############################
# Define model SETTINGS
##############################

# let us consider a larger validation data (say 30%) set for the following
# resample_trnvalidx!(dtmtable,.7)

updateSettingsMod!(sett,minWeight=-0.02,model_type="boosted_tree",iterations=10,learningRate=0.05,subsampling_features_prop=.7,calculatePoissonError=false)
# sett.scorebandsstartingpoints=[1,50,100,150,200,300,400,500,600,700,850,900,950]

# resM is the resulting Model
resultingFiles, resM = dtm(dtmtable, sett)

predDF = predict(resM, dtmtable.features)

estimatedResiduals = predDF[:SmoothedEstimate];

# let us see if this improved the error
# we will simply multiply the GLM fit, with the fitted residual
glmFitImprovedByBoosting = fullData[:GLMfitted] .* estimatedResiduals

errsGLMBoosted = poissonError(fullData[:ClaimNb], glmFitImprovedByBoosting)
trnAvgErrorGLMBoosted = sum(errsGLMBoosted[trnBitArrayIdx]) / sum(trnBitArrayIdx)
valAvgErrorGLMBoosted = sum(errsGLMBoosted[valBitArrayIdx]) / sum(valBitArrayIdx)

valAvgErrorGLM - valAvgErrorGLMBoosted
trnAvgErrorGLM - trnAvgErrorGLMBoosted

# let us change the score bands a little bit as there is not enough differentiation
sett.scorebandsstartingpoints = [1,100,250,500,750,900]

# resM is the resulting Model
resultingFiles, resM = dtm(dtmtable, sett)


text:
using MultivariatePolynomials
using JuMP

struct ConeLyap{T, P<:AbstractPolynomial{T}, S}
    p::P
    Q::Matrix{S}
    b::Vector{S}
    c::Vector{Float64}
    H::Matrix{Float64}
    #L::Matrix{JuMP.Variable}
    #λinv::Vector{JuMP.Variable}
end

using RecipesBase

export Ellipsoid

struct Ellipsoid{T}
    Q::Matrix{T}
    c::Vector{T}
end

struct LiftedEllipsoid{T}
    P::Matrix{T}
end

@recipe function f(ell::Ellipsoid)
    @assert Base.LinAlg.checksquare(ell.Q) == 2
    αs = linspace(0, 2π, 64)
    ps = [[cos(α), sin(α)] for α in αs]
    r = [sqrt(dot(p, ell.Q * p)) for p in ps]
    seriestype --> :shape
    legend --> false
    ell.c[1] .+ cos.(αs) ./ r, ell.c[2] .+ sin.(αs) ./ r
end


text:
# Generates, compiles and runs a C file (gen_consts.c) which prints the necessary type
# declarations and constants in Julia syntax.

# For the "handle" constants (MPI_Info, MPI_Op, etc.) we print the Fortran values, since
# for some implementations (OpenMPI), the C versions are only available at link-time (so
# may not be actually constant).

# It may be possible to work around this by using `cglobal` instead, but we can't know
# what symbols the constants resolve to without a C compiler.

@info "Unknown MPI ABI: building constants file"

MPI_LIBRARY_PATH = dirname(dlpath(libmpi))

MPI_PATH = get(ENV, "JULIA_MPI_PATH", "")
MPI_INCLUDE_PATH = get(ENV, "JULIA_MPI_INCLUDE_PATH", MPI_PATH == "" ? "" : joinpath(MPI_PATH,"include"))
mpicc = get(ENV, "JULIA_MPICC", MPI_PATH == "" ? "mpicc" : joinpath(MPI_PATH, "bin", "mpicc"))

if haskey(ENV, "JULIA_MPI_CFLAGS")
    CFLAGS = Base.shell_split(ENV["JULIA_MPI_CFLAGS"])
else
    lname, = split(basename(libmpi),'.')
    if startswith(lname, "lib")
        lname = lname[4:end]
    end
    CFLAGS = `-l$(lname) -L$(MPI_LIBRARY_PATH)`
    if MPI_INCLUDE_PATH != ""
        CFLAGS = `$CFLAGS -I$(MPI_INCLUDE_PATH)`
    end
end

# build gen_consts.c

MPI_op_consts = [
    :MPI_OP_NULL,
    :MPI_BAND,
    :MPI_BOR,
    :MPI_BXOR,
    :MPI_LAND,
    :MPI_LOR,
    :MPI_LXOR,
    :MPI_MAX,
    :MPI_MIN,
    :MPI_SUM,
    :MPI_PROD,
    :MPI_MAXLOC,
    :MPI_MINLOC,
    :MPI_REPLACE,
]

if MPI_VERSION >= v"3.0"
    push!(MPI_op_consts, :MPI_NO_OP)
end

MPI_datatype_consts = [
    :MPI_DATATYPE_NULL,
    :MPI_PACKED,
    :MPI_CHAR,
    :MPI_SIGNED_CHAR,
    :MPI_UNSIGNED_CHAR,
    :MPI_WCHAR,
    :MPI_BYTE,
    :MPI_SHORT,
    :MPI_UNSIGNED_SHORT,
    :MPI_INT,
    :MPI_UNSIGNED,
    :MPI_LONG,
    :MPI_UNSIGNED_LONG,
    :MPI_FLOAT,
    :MPI_DOUBLE,
]

if MPI_VERSION >= v"2.2"
    append!(MPI_datatype_consts, [
        :MPI_INT8_T,
        :MPI_INT16_T,
        :MPI_INT32_T,
        :MPI_INT64_T,
        :MPI_UINT8_T,
        :MPI_UINT16_T,
        :MPI_UINT32_T,
        :MPI_UINT64_T,
        :MPI_C_FLOAT_COMPLEX,
        :MPI_C_DOUBLE_COMPLEX,
    ])
end

MPI_handle = [
    :MPI_Info => [
        :MPI_INFO_NULL,
    ],
    :MPI_Win => [
        :MPI_WIN_NULL,
    ],
    :MPI_Comm => [
        :MPI_COMM_NULL,
        :MPI_COMM_SELF,
        :MPI_COMM_WORLD,
    ],
    :MPI_Request => [
        :MPI_REQUEST_NULL,
    ],
    :MPI_File => [
        :MPI_FILE_NULL,
    ],
    :MPI_Op => MPI_op_consts,
    :MPI_Datatype => MPI_datatype_consts
]

MPI_Cints = [
    :MPI_THREAD_SINGLE,
    :MPI_THREAD_FUNNELED,
    :MPI_THREAD_SERIALIZED,
    :MPI_THREAD_MULTIPLE,
    :MPI_PROC_NULL,
    :MPI_ANY_SOURCE,
    :MPI_ANY_TAG,
    :MPI_SUCCESS,
    :MPI_UNDEFINED,
    :MPI_LOCK_EXCLUSIVE,
    :MPI_LOCK_SHARED,
    :MPI_MAX_INFO_KEY,
    :MPI_MAX_INFO_VAL,
    :MPI_TAG_UB,
    :MPI_COMM_TYPE_SHARED,
    :MPI_ORDER_C,
    :MPI_ORDER_FORTRAN,
    :MPI_UNIVERSE_SIZE,
    :MPI_MAX_ERROR_STRING,
    :MPI_SUCCESS,
    :MPI_MODE_RDONLY,
    :MPI_MODE_RDWR,
    :MPI_MODE_WRONLY,
    :MPI_MODE_CREATE,
    :MPI_MODE_EXCL,
    :MPI_MODE_DELETE_ON_CLOSE,
    :MPI_MODE_UNIQUE_OPEN,
    :MPI_MODE_SEQUENTIAL,
    :MPI_MODE_APPEND,
    :MPI_SEEK_SET,
    :MPI_SEEK_CUR,
    :MPI_SEEK_END,
    :MPI_IDENT,
    :MPI_CONGRUENT,
    :MPI_SIMILAR,
    :MPI_UNEQUAL,
]

MPI_pointers = [
    :MPI_BOTTOM,
    :MPI_IN_PLACE,
    :MPI_STATUS_IGNORE,
    :MPI_STATUSES_IGNORE,
    ]

const libptr = dlopen_e(libmpi)

open("gen_consts.c","w") do f
    print(f,"""
#include <stdio.h>
#include <stddef.h>
#include <inttypes.h>
#include "mpi.h"

int main(int argc, char *argv[]) {
    MPI_Init(&argc, &argv);

    FILE *fptr;
    fptr = fopen("consts.jl", "a");

""")

    println(f,"  fprintf(fptr, \"const MPI_Aint = Int%d\\n\", 8*(int)sizeof(MPI_Aint));")
    println(f,"  fprintf(fptr, \"const MPI_Offset = Int%d\\n\", 8*(int)sizeof(MPI_Offset));")
    println(f,"  fprintf(fptr, \"const MPI_Count = Int%d\\n\", 8*(int)sizeof(MPI_Count));")

    println(f,"  fprintf(fptr, \"const MPI_Status_size = %d\\n\", (int)sizeof(MPI_Status));")
    println(f,"  fprintf(fptr, \"const MPI_Status_Source_offset = %d\\n\", (int)offsetof(MPI_Status, MPI_SOURCE));")
    println(f,"  fprintf(fptr, \"const MPI_Status_Tag_offset = %d\\n\", (int)offsetof(MPI_Status, MPI_TAG));")
    println(f,"  fprintf(fptr, \"const MPI_Status_Error_offset = %d\\n\", (int)offsetof(MPI_Status, MPI_ERROR));")

    for (T,constants) in MPI_handle
        T_f2c = T == :MPI_Datatype ? :MPI_Type_f2c : Symbol(T, :_f2c)
        T_c2f = T == :MPI_Datatype ? :MPI_Type_c2f : Symbol(T, :_c2f)
        if Libdl.dlsym_e(libptr, T_f2c) == C_NULL
            println(f,"  fprintf(fptr, \"const $T = Cint\\n\");")
            for constant in constants
                println(f,"  fprintf(fptr, \"const $constant = Cint(%i)\\n\", $constant);")
            end
        else
            println(f,"  fprintf(fptr, \"const $T = Ptr{Cvoid}\\n\");")
            println(f,"  fprintf(fptr, \"$(Symbol(T,:_f2c))(c::Cint) = ccall(:$T_f2c,$T,(Cint,),c)\\n\");")
            for constant in constants
                println(f,"  fprintf(fptr, \"const $constant = Cint(%i)\\n\", $T_c2f($constant));")
            end
        end
    end

    for constant in MPI_Cints
        println(f,"  fprintf(fptr, \"const $constant = Cint(%i)\\n\", (int)($constant));")
    end

    for constant in MPI_pointers
        if Sys.WORD_SIZE == 32
            println(f,"  fprintf(fptr, \"const $constant = reinterpret(SentinelPtr, 0x%08\" PRIx32 \")\\n\", (uint32_t)($constant));")
        else
            println(f,"  fprintf(fptr, \"const $constant = reinterpret(SentinelPtr, 0x%016\" PRIx64 \")\\n\", (uint64_t)($constant));")
        end
    end

    print(f,"""
  fclose(fptr);
  MPI_Finalize();
  return 0;
}
""")
end

# compile gen_consts.c
run(`$mpicc gen_consts.c -o gen_consts $CFLAGS`)

open("consts.jl", write=true) do f
    println(f, "# This file is automatically generated")
    println(f, "# Do not edit")
    println(f)
end

# execute gen_consts
run(`$mpiexec_cmd -n 1 ./gen_consts`)


text:
using NearestNeighbors
using Distances: chebyshev
using SpecialFunctions: digamma
using LsqFit: curve_fit
using StatsBase: autocor

export estimate_delay
export estimate_dimension
export mutinfo, mutinfo_delaycurve
export estimate_dimension, stochastic_indicator

#####################################################################################
#                                Mutual Information                                 #
#####################################################################################
"""
    mutinfo(k, X1, X2[, ..., Xm]) -> MI

Calculate the mutual information `MI` of the given vectors
`X1, X2, ...`, using `k` nearest-neighbors.

The method follows the second algorithm ``I^{(2)}`` outlined by Kraskov in [1].

## References
[1] : A. Kraskov *et al.*, [Phys. Rev. E **69**, pp 066138 (2004)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.69.066138)

## Performance Notes
This functin gets very slow for large `k`.

See also [`estimate_delay`](@ref) and [`mutinfo_delaycurve`](@ref).
"""
function mutinfo(k, Xm::Vararg{<:AbstractVector,M}) where M
    @assert M > 1
    @assert (size.(Xm,1) .== size(Xm[1],1)) |> prod
    k += 1
    N = size(Xm[1],1)
    invN = 1/N

    d = Dataset(Xm...)
    tree = KDTree(d.data, Chebyshev())

    n_x_m = zeros(M)

    Xm_sp = zeros(Int, N, M)
    Xm_revsp = zeros(Int, N, M)
    for m in 1:M
        Xm_sp[:,m] .= sortperm(Xm[m]; alg=QuickSort)
        Xm_revsp[:,m] .= sortperm(Xm_sp[:,m]; alg=QuickSort)
    end

    I = digamma(k) - (M-1)/k + (M-1)*digamma(N)

    nns = (x = knn(tree, d.data, k)[1]; [ind[1] for ind in x])

    I_itr = zeros(M)
    # Makes more sense computationally to loop over N rather than M
    for i in 1:N
        ϵ = abs.(d[nns[i]] - d[i])./2

        for m in 1:M # this loop takes 8% of time
            hb = lb = Xm_revsp[i,m]
            while abs(Xm[m][Xm_sp[hb,m]] - Xm[m][i]) <= ϵ[m] && hb < N
                hb += 1
            end
            while abs(Xm[m][Xm_sp[lb,m]] - Xm[m][i]) <= ϵ[m] && lb > 1
                lb -= 1
            end
            n_x_m[m] = hb - lb
        end

        I_itr .+= digamma.(n_x_m)
    end

    I_itr .*= invN

    I -= sum(I_itr)

    return max(0, I)
end

"""
    mutinfo_delaycurve(x; maxtau=100, k=1)

Return the [`mutinfo`](@ref) between `x` and itself for delays of `1:maxtau`.
"""
function mutinfo_delaycurve(X::AbstractVector; maxtau=100, k=1)
    I = zeros(maxtau)

    @views for τ in 1:maxtau
        I[τ] = mutinfo(k, X[1:end-τ],X[τ+1:end])
    end

    return I
end

#####################################################################################
#                               Estimate Delay Times                                #
#####################################################################################
"""
    localextrema(y) -> max_ind, min_ind
Find the local extrema of given array `y`, by scanning point-by-point. Return the
indices of the maxima (`max_ind`) and the indices of the minima (`min_ind`).
"""
function localextrema end
@inbounds function localextrema(y)
    l = length(y)
    i = 1
    maxargs = Int[]
    minargs = Int[]
    if y[1] > y[2]
        push!(maxargs, 1)
    elseif y[1] < y[2]
        push!(minargs, 1)
    end

    for i in 2:l-1
        left = i-1
        right = i+1
        if  y[left] < y[i] > y[right]
            push!(maxargs, i)
        elseif y[left] > y[i] < y[right]
            push!(minargs, i)
        end
    end

    if y[l] > y[l-1]
        push!(maxargs, l)
    elseif y[l] < y[l-1]
        push!(minargs, l)
    end
    return maxargs, minargs
end

function exponential_decay_extrema(c::AbstractVector)
    ac = abs.(c)
    ma, mi = localextrema(ac)
    # ma start from 1 but correlation is expected to start from x=0
    ydat = ac[ma]; xdat = ma .- 1
    # Do curve fit from LsqFit
    model(x, p) = @. exp(-x/p[1])
    decay = curve_fit(model, xdat, ydat, [1.0]).param[1]
    return decay
end

function exponential_decay(c::AbstractVector)
    # Do curve fit from LsqFit
    model(x, p) = @. exp(-x/p[1])
    decay = curve_fit(model, 0:length(c)-1, abs.(c), [1.0]).param[1]
    return decay
end

"""
    estimate_delay(s, method::String) -> τ

Estimate an optimal delay to be used in [`reconstruct`](@ref).
Return the exponential decay time `τ` rounded to an integer.

The `method` can be one of the following:

* `"first_zero"` : find first delay at which the auto-correlation function becomes 0.
* `"first_min"` : return delay of first minimum of the auto-correlation function.
* `"exp_decay"` : perform an exponential fit to the `abs.(c)` with `c` the
  auto-correlation function of `s`.
* `"mutual_inf"` : return the first minimum of the mutual information function
  (see [`mutinfo_delaycurve`](@ref)).
  This option also has the following keyword arguments:
    * `maxtau::Integer=100` : stops the delay calculations after the given `maxtau`.
    * `k::Integer=1` : the number of nearest-neighbors to include.

*WARNING* - `"mutual_inf"` has bad results with data from maps. In addition
it is much slower than the other alternatives.
"""
function estimate_delay(x::AbstractVector, method::String; maxtau=100, k=1)
    method ∈ ("first_zero", "first_min", "exp_decay", "mutual_inf") ||
        throw(ArgumentError("Unknown method"))

    if method=="first_zero"
        c = autocor(x, 0:length(x)÷10; demean=true)
        i = 1
        # Find 0 crossing:
        while c[i] > 0
            i+= 1
            i == length(c) && break
        end
        return i

    elseif method=="first_min"
        c = autocor(x, 0:length(x)÷10, demean=true)
        i = 1
        # Find min crossing:
        while  c[i+1] < c[i]
            i+= 1
            i == length(c)-1 && break
        end
        return i
    elseif method=="exp_decay"
        c = autocor(x, demean=true)
        # Find exponential fit:
        τ = exponential_decay(c)
        return round(Int,τ)
    elseif method=="mutual_inf"
        m = mutinfo(k, x, x)
        L = length(x)
        for i=1:maxtau
            n = mutinfo(k, view(x, 1:L-i), view(x, 1+i:L))
            n > m && return i
            m = n
        end
    end
end

#####################################################################################
#                                Estimate Dimension                                 #
#####################################################################################
"""
    estimate_dimension(s::AbstractVector, τ:Int, Ds = 1:5) -> E₁s

Compute a quantity that can estimate an optimal amount of
temporal neighbors `D` to be used in [`reconstruct`](@ref).

## Description
Given the scalar timeseries `s` and the embedding delay `τ` compute the
values of `E₁` for each `D ∈ Ds`, according to Cao's Method (eq. 3 of [1]).
Please be aware that in **DynamicalSystems.jl** `D` stands for the amount of temporal
neighbors, *not* the dimension as in [1]. The dimension is `D+1`.

Return the vector of all computed `E₁`s. To estimate a good value for `D` from this,
find `D` for which the value `E₁` saturates at some value around 1.

*Note: This method does not work for datasets with perfectly periodic signals.*

## References

[1] : Liangyue Cao, [Physica D, pp. 43-50 (1997)](https://www.sciencedirect.com/science/article/pii/S0167278997001188?via%3Dihub)
"""
function estimate_dimension(s::AbstractVector{T}, τ::Int, Ds = 1:5) where {T}
    E1s = zeros(T, length(Ds))
    aafter = zero(T)
    aprev = _average_a(s, Ds[1], τ)
    for (i, D) ∈ enumerate(Ds)
        aafter = _average_a(s, D+1, τ)
        E1s[i] = aafter/aprev
        aprev = aafter
    end
    return E1s
end
# then use function `saturation_point(Ds, E1s)` from ChaosTools

function _average_a(s::AbstractVector{T},D,τ) where T
    #Sum over all a(i,d) of the Ddim Reconstructed space, equation (2)
    R1 = reconstruct(s,D+1,τ)
    R2 = reconstruct(s[1:end-τ],D,τ)
    tree2 = KDTree(R2)
    nind = (x = knn(tree2, R2.data, 2)[1]; [ind[1] for ind in x])
    e=0.
    for (i,j) in enumerate(nind)
        δ = norm(R2[i]-R2[j], Inf)
        #If R2[i] and R2[j] are still identical, choose the next nearest neighbor
        if δ == 0.
            j = knn(tree2, R2[i], 3, true)[1][end]
            δ = norm(R2[i]-R2[j], Inf)
        end
        e += norm(R1[i]-R1[j], Inf) / δ
    end
    return e / length(R1)
end

function dimension_indicator(s,D,τ) #this is E1, equation (3) of Cao
    return _average_a(s,D+1,τ)/_average_a(s,D,τ)
end

"""
    stochastic_indicator(s::AbstractVector, τ:Int, Ds = 1:4) -> E₂s

Compute an estimator for apparent randomness in a reconstruction with `Ds` temporal
neighbors.

## Description
Given the scalar timeseries `s` and the embedding delay `τ` compute the
values of `E₂` for each `D ∈ Ds`, according to Cao's Method (eq. 5 of [1]).

Use this function to confirm that the
input signal is not random and validate the results of [`estimate_dimension`](@ref).
In the case of random signals, it should be `E₂ ≈ 1 ∀ D`.

## References

[1] : Liangyue Cao, [Physica D, pp. 43-50 (1997)](https://www.sciencedirect.com/science/article/pii/S0167278997001188?via%3Dihub)
"""
function stochastic_indicator(s::AbstractVector{T},τ, Ds=1:4) where T # E2, equation (5)
    #This function tries to tell the difference between deterministic
    #and stochastic signals
    #Calculate E* for Dimension D+1
    E2s = Float64[]
    for D ∈ Ds
        R1 = reconstruct(s,D+1,τ)
        tree1 = KDTree(R1[1:end-1-τ])
        method = FixedMassNeighborhood(2)

        Es1 = 0.
        nind = (x = neighborhood(R1[1:end-τ], tree1, method); [ind[1] for ind in x])
        for  (i,j) ∈ enumerate(nind)
            Es1 += abs(R1[i+τ][end] - R1[j+τ][end]) / length(R1)
        end

        #Calculate E* for Dimension D
        R2 = reconstruct(s,D,τ)
        tree2 = KDTree(R2[1:end-1-τ])
        Es2 = 0.
        nind = (x = neighborhood(R2[1:end-τ], tree2, method); [ind[1] for ind in x])
        for  (i,j) ∈ enumerate(nind)
            Es2 += abs(R2[i+τ][end] - R2[j+τ][end]) / length(R2)
        end
        push!(E2s, Es1/Es2)
    end
    return E2s
end


text:
abstract type AbstractPooledResource{T} end

mutable struct PooledResource{T} <: AbstractPooledResource{T}
    resource::T
    ref_count::Int
    dispose::Function
    function PooledResource(r::T, dispose::Function) where T
        new{T}(r, 1, dispose)
    end
end

resource(r) = r.resource
ref_count(r) = r.ref_count

function release!(r)
    if r.ref_count < 1
        throw(ArgumentError("Cannot release a resource with ref_count zero!"))
    end
    @debug "Releasing $(r)"
    r.ref_count -= 1
    if r.ref_count == 0
        @debug "Disposing $(r)"
        r.dispose(r)
    end
end

function retain!(r)
    if r.ref_count < 1
        throw(ArgumentError("Cannot retain a resource with ref_count zero!"))
    end
    @debug "Retaining $(r)"
    r.ref_count += 1
end


text:
using Base.Test
include("print_dashed_line.jl")

include(joinpath("..", "src", "MolecularGeometry.jl"))
using MolecularGeometry
import MolecularGeometry.gage

include(joinpath("function","runtests.jl"))


text:
abstract Spline{T} <: ScalarFunc{T}

# cubic 1D spline
type Spline1D{T,V} <: Spline{T}
    x::V
    y::V
    y2::V
    Spline1D(x::AbstractVector{T}, y::AbstractVector{T}, y2::AbstractVector{T}) = new(x,y,y2)
end

Spline1D{T}(x::AbstractVector{T}, y::AbstractVector{T}, y2::AbstractVector{T}) =
    Spline1D{eltype(x),typeof(x)}(x,y,y2)

typealias Spline1Df Spline1D{Float64,Vector{Float64}}

@fastmath function Spline1D(x::AbstractVector{Float64}, y::AbstractVector{Float64})
    @assert length(x) == length(y)
    local sig,p
    y2 = similar(y)
    n = length(x)
    u = Array{Float64}(n-1)
    y2[1] = u[1] = 0.0
    @inbounds for i in 2:n-1
        sig = (x[i]-x[i-1])/(x[i+1]-x[i-1])
        p = sig*y2[i-1]+2.0
        y2[i]=(sig-1.0)/p
        u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1])
        u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p
    end
    qn=un=0.0
    y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0)
    @inbounds for k in n-1:-1:1  # This is the backsubstitution loop of the tridiagonal
        y2[k]=y2[k]*y2[k+1]+u[k];
    end
    Spline1D(x,y,y2)
end

function Spline1D(x::AbstractVector{Float64}, y::AbstractVector{Float64}, num::Int)
    n = length(x)
    ind = num < n ? round(Int,linspace(1,n,num)) : Colon()
    Spline1D(x[ind],y[ind])
end

Base.show(io::IO, f::Spline1D) = write(io, "Spline1D $(length(f.x))-points on interval [$(f.x[1]),$(f.x[end])]")

@fastmath function locate(xx::AbstractVector{Float64}, x::Float64)
    local k
    klo,khi = 1,length(xx)
    while khi-klo > 1
        k = (khi+klo) >> 1  # mean of khi and klo
        if xx[k] > x
            khi = k
        else
            klo = k
        end
    end
    klo,khi
end

@fastmath function Base.call(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    a*f.y[klo] + b*f.y[khi] + ((a*a*a-a)*f.y2[klo]+(b*b*b-b)*f.y2[khi])*(h*h)/6.
end

@vectorize2 Base.call Spline1D Float64 Float64

# function Base.call(f::Spline1D, dst::AbstractVector{Float64}, x::AbstractVector{Float64})
#     @assert length(dst) == length(x)
#     for i in 1:length(x)
#         dst[i] = f(x[i])
#     end
# end

# function Base.call(f::Spline1D, x::AbstractVector{Float64})
#     dst = similar(x)
#     f(dst,x)
#     dst
# end

@fastmath function deriv1(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    (f.y[khi] - f.y[klo])/h + ((1.-3.*a*a)*f.y2[klo] + (3*b*b-1.)*f.y2[khi])*h/6.
end

@vectorize2 deriv1 Spline1D Float64 Float64

# function deriv1(f::Spline1D, dst::AbstractVector{Float64}, x::AbstractVector{Float64})
#     @assert length(dst) == length(x)
#     @inbounds for i in 1:length(x)
#         dst[i] = deriv1(f,x[i])
#     end
# end

# function deriv1(f::Spline1D, x::AbstractVector{Float64})
#     dst = similar(x)
#     deriv1(f,dst,x)
#     dst
# end

@fastmath function deriv2(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    a*f.y2[klo] + b*f.y2[khi]
end

@vectorize2 deriv2 Spline1D Float64 Float64

# function deriv2(f::Spline1D, dst::AbstractVector{Float64}, x::AbstractVector{Float64})
#     @assert length(dst) == length(x)
#     @inbounds for i in 1:length(x)
#         dst[i] = deriv2(f,x[i])
#     end
# end

# function deriv2(f::Spline1D, x::AbstractVector{Float64})
#     dst = similar(x)
#     deriv2(f,dst,x)
#     dst
# end

@fastmath function eval_deriv1(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    a*f.y[klo] + b*f.y[khi] + ((a*a*a-a)*f.y2[klo]+(b*b*b-b)*f.y2[khi])*(h*h)/6.,
        (f.y[khi] - f.y[klo])/h + ((1.-3.*a*a)*f.y2[klo] + (3*b*b-1.)*f.y2[khi])*h/6.
end

@fastmath function eval_deriv2(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    a*f.y[klo] + b*f.y[khi] + ((a*a*a-a)*f.y2[klo]+(b*b*b-b)*f.y2[khi])*(h*h)/6.,
        a*f.y2[klo] + b*f.y2[khi]
end

@fastmath function eval_deriv12(f::Spline1D, x::Float64)
    klo,khi = locate(f.x,x)
    h = f.x[khi] - f.x[klo]  # @assert h != 0.0
    a = (f.x[khi]-x)/h
    b = (x-f.x[klo])/h
    a*f.y[klo] + b*f.y[khi] + ((a*a*a-a)*f.y2[klo]+(b*b*b-b)*f.y2[khi])*(h*h)/6.,
        (f.y[khi] - f.y[klo])/h + ((1.-3.*a*a)*f.y2[klo] + (3*b*b-1.)*f.y2[khi])*h/6.,
            a*f.y2[klo] + b*f.y2[khi]
end

precompile(Spline1D, (Vector{Float64}, Vector{Float64}))
precompile(locate, (Vector{Float64}, Float64))
precompile(Base.call, (Spline1D, Float64))
precompile(Base.call, (Spline1D, Vector{Float64}, Float64))
precompile(Base.call, (Spline1D, Vector{Float64}))
precompile(deriv1, (Spline1D, Float64))
precompile(deriv1, (Spline1D, Vector{Float64}, Float64))
precompile(deriv1, (Spline1D, Vector{Float64}))
precompile(deriv2, (Spline1D, Float64))
precompile(deriv2, (Spline1D, Vector{Float64}, Float64))
precompile(deriv2, (Spline1D, Vector{Float64}))
precompile(eval_deriv1, (Spline1D, Float64))
precompile(eval_deriv2, (Spline1D, Float64))
precompile(eval_deriv12, (Spline1D, Float64))

max_stars_repo_name: johncsnyder/Element


text:
@testset "17.letter-combinations-of-a-phone-number.jl" begin
    @test letter_combinations("23") ==
        ["ad", "ae", "af", "bd", "be", "bf", "cd", "ce", "cf"]
    @test isempty(letter_combinations(""))
    @test letter_combinations("8") == ["t", "u", "v"]
end


text:
using Pkg
Pkg.activate(pwd())

# # Linear regression revisited
# ## Hypothesis testing

using Random
using Statistics
using LinearAlgebra
using Plots

function plot_histogram(xs, f; kwargs...)
    plt = histogram(xs;
        label="Sampled density",
        xlabel = "x",
        ylabel = "pdf(x)",
        nbins = 85,
        normalize = :pdf,
        opacity = 0.5,
        kwargs...
    )

    plot!(plt, range(minimum(xs), maximum(xs); length=100), f;
        label="True density",
        line=(4,:black),
    )

    return plt
end

#+

Random.seed!(666)
n = 1000
xs = randn(n)

# ### Exercise:
# Use the ``t``-test to verify whether the samples were generated from a distribution with
# zero mean.
#
# **Hints:**
# - The Student's distribution is invoked by `TDist()`.
# - The probability ``\mathbb P(T\le t)`` equals to the
# [distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function)
# ``F(t)``, which can be called by `cdf`.
#
# ---
# ### Solution:

# ---

using HypothesisTests

OneSampleTTest(xs)

# ## Linear models

using RDatasets

wages = dataset("plm", "Snmesp")
wages.W = 2. .^ (wages.W)

#+

X = Matrix(wages[:, [:N, :Y, :I, :K, :F]])
X = hcat(ones(size(X,1)), X)
y = wages[:, :W]

w0 = (X'*X) \ (X'*y)

#+

using GLM

model = lm(@formula(W ~ 1 + N + Y + I + K + F), wages)

# ### Exercise:
# Check that the solution computed by hand and by `lm` are the same.
#
# Then remove the feature with the highest ``p``-value and observe whether there was any
# performance drop. The performance is usually evaluated by the
# [coeffient of determination](https://en.wikipedia.org/wiki/Coefficient_of_determination)
# denoted by ``R^2\in[0,1]``. Its higher values indicate a better model.
#
# **Hint**: Use functions `coef` and `r2`.
#
# ---
# ### Solution:

# ---

y_hat = predict(model)

plot_histogram(y_hat, x -> pdf(Normal(mean(y_hat), std(y_hat)), x)) |> display
plot_histogram(y_hat, x -> pdf(fit(Normal, y_hat), x)) |> display

#+

test_normality = ExactOneSampleKSTest(y_hat, Normal(mean(y_hat), std(y_hat)))

# ## Generalized linear models

model = glm(@formula(W ~ 1 + N + Y + I + K + F), wages, InverseGaussian(), SqrtLink())

# ### Exercise:
# Create the scatter plot of predictions and labels. Do not use the `predict` function.
#
# ---
# ### Solution:


text:
#using ExoplanetsSysSim
#using StatsFuns
#using JLD
#using DataFrames
#import ExoplanetsSysSim.StellarTable.df
import Compat: UTF8String, ASCIIString

function calc_distance_model_ks(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  d1 = calc_distance_num_planets(summary1,summary2,sim_param,verbose=verbose)
  d2 = calc_distance_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d3 = calc_distance_num_planets_binned(summary1,summary2,sim_param,verbose=verbose)
  #d4 = calc_distance_mean_std_log_period_depth(summary1,summary2,sim_param,verbose=verbose)
  d5 = calc_distance_ks_period(summary1,summary2,sim_param,verbose=verbose)
  d6 = calc_distance_ks_depth(summary1,summary2,sim_param,verbose=verbose)
  d7 = calc_distance_ks_period_ratios(summary1,summary2,sim_param,verbose=verbose)
  d8 = calc_distance_ks_duration_ratios(summary1,summary2,sim_param,verbose=verbose)
  return vcat(d1, d2, d5, d6, d7, d8)
end

function calc_distance_model_kl(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  d1 = calc_distance_kl_num_planets(summary1,summary2,sim_param,verbose=verbose)
  d2 = calc_distance_kl_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d2 = calc_distance_hellinger_num_n_tranet_systems(summary1,summary2,sim_param,verbose=verbose)
  #d3 = calc_distance_num_planets_binned(summary1,summary2,sim_param,verbose=verbose)
  #d4 = calc_distance_mean_std_log_period_depth(summary1,summary2,sim_param,verbose=verbose)
  d5 = calc_distance_kl_period(summary1,summary2,sim_param,verbose=verbose)
  d6 = calc_distance_kl_depth(summary1,summary2,sim_param,verbose=verbose)
  d7 = calc_distance_kl_period_ratios(summary1,summary2,sim_param,verbose=verbose)
  d8 = calc_distance_kl_duration_ratios(summary1,summary2,sim_param,verbose=verbose)
  return vcat(d1, d2, d5, d6, d7, d8)
end
calc_distance_model(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false) = calc_distance_model_kl(summary1,summary2,sim_param,verbose=verbose)

function calc_distance_num_planets(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = haskey(summary1.stat,"num_tranets") ? summary1.stat["num_tranets"] : summary1.stat["expected planets detected"]
    np2 = haskey(summary2.stat,"num_tranets") ? summary2.stat["num_tranets"] : summary2.stat["expected planets detected"]
    #println("np1 = ",np1,", np2 = ",np2)

    dist_np = dist_L1_abs(np1/summary1.stat["num targets"], np2/summary2.stat["num targets"])
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", d = ",dist_np)
  return dist_np
end

function calc_distance_num_planets_binned(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = haskey(summary1.stat,"planets table") ? summary1.stat["planets table"] : summary1.stat["expected planets table"]
    np2 = haskey(summary2.stat,"planets table") ? summary2.stat["planets table"] : summary2.stat["expected planets table"]
    #println("np1 = ",np1,", np2 = ",np2)

    dist_np_bin = zeros(length(np1))
    for n in 1:length(np1)
      dist_np_bin[n] = dist_L1_abs(np1[n]/summary1.stat["num targets"], np2[n]/summary2.stat["num targets"])
      #println("True # [Bin ", n,"] = ",np1[n],", Expected # [Bin ", n,"] = ",np2[n])
    end
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", dist = ",dist_np_bin)
  return dist_np_bin
end

function calc_distance_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    d[n] = dist_L1_abs(summary1.stat["num n-tranet systems"][n]/summary1.stat["num targets"], summary2.stat["num n-tranet systems"][n]/summary2.stat["num targets"])
  end
  return d
end

function calc_distance_mean_std_log_period_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    d1 = dist_L1_abs(summary1.stat["mean log10 P"], summary2.stat["mean log10 P"])
    d2 = dist_L1_abs(summary1.stat["mean log10 depth"], summary2.stat["mean log10 depth"])
    d3 = dist_L1_abs(summary1.stat["std log10 P"], summary2.stat["std log10 P"])
    d4 = dist_L1_abs(summary1.stat["std log10 depth"], summary2.stat["std log10 depth"])
    return [d1, d2, d3, d4]
end

# compute supremum of differences between empirical cdfs.
# Borrowed from JuliaStats/HypothesisTests.jl
function ksstats{T<:Real, S<:Real}(x::AbstractVector{T}, y::AbstractVector{S})
    n_x, n_y = length(x), length(y)
    sort_idx = sortperm([x; y])
    pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)
end

##### Function to compute the KS distance for discrete arrays (i.e. of the planet multiplicities in the systems):
function ksstats_ints{T<:Int, S<:Int}(x::AbstractVector{T}, y::AbstractVector{S})
    n_x, n_y = length(x), length(y)
    xy_max = max(maximum(x), maximum(y)) #maximum of both x and y
    x_counts = [sum(x .== n) for n in 1:xy_max]
    y_counts = [sum(y .== n) for n in 1:xy_max]
    pdf_diffs = x_counts./n_x .- y_counts./n_y
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)
end

##### Function to compute the 2-sample Anderson-Darling (AD) distance between two empirical CDFs (continuous distributions):
function ADstats{T<:Real, S<:Real}(x::AbstractVector{T}, y::AbstractVector{S})
    #This function computes the AD distance according to A. N. Pettitt (1976) Eq. (1.2)
    n, m = length(x), length(y)
    N = n + m
    sort_idx = sortperm([x; y]) #array of indices that would sort the combined array
    M_i_diffs = [ones(n); zeros(m)][sort_idx]
    M_i_array = cumsum(M_i_diffs)[1:end-1] #array of M_i except for last element, i.e. from i=1 to i=N-1
    i_array = 1:(N-1) #array of i from i=1 to i=N-1
    AD_dist = (1./(n*m))*sum(((M_i_array*N - n*i_array).^2.)./(i_array.*(N - i_array))) #AD distance
    return AD_dist
end

function ADstats_mod{T<:Real, S<:Real}(x::AbstractVector{T}, y::AbstractVector{S})
    #This function is the same as 'ADstats' except without the factor of 'nm/N' before the integral
    n, m = length(x), length(y)
    N = n + m
    sort_idx = sortperm([x; y]) #array of indices that would sort the combined array
    M_i_diffs = [ones(n); zeros(m)][sort_idx]
    M_i_array = cumsum(M_i_diffs)[1:end-1] #array of M_i except for last element, i.e. from i=1 to i=N-1
    i_array = 1:(N-1) #array of i from i=1 to i=N-1
    AD_dist = (N/((n*m)^2.))*sum(((M_i_array*N - n*i_array).^2.)./(i_array.*(N - i_array))) #AD distance
    return AD_dist
end

##### Function to compute the Cressie-Read Power Divergence (CRPD) statistic for the distributions of observed planet multiciplities:
function CRPDstats{T<:Int, S<:Int}(En::AbstractVector{T}, On::AbstractVector{S})
    #En and On must be arrays of the total numbers of systems with 1,2,3,... observed planets, in the simulated (i.e. expected) and the actual (i.e. observed Kepler) data, respectively
    E_array = En./sum(En) #normalized numbers (fractions) of simulated systems with 1,2,3,... observed planets
    O_array = On./sum(On) #normalized numbers (fractions) of actual Kepler systems with 1,2,3,... observed planets
    rho = 0.
    for (i,E_i) in enumerate(E_array)
        if En[i] != 0
            rho += O_array[i]*((O_array[i]/E_array[i])^(2./3.) - 1.)
        end
    end
    rho = (9./5.)*rho
    return rho
end

function calc_distance_ks_period(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["P list"]
  samp2 = summary2.stat["P list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["depth list"]
  samp2 = summary2.stat["depth list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_period_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["period_ratio_list"]
  samp2 = summary2.stat["period_ratio_list"]
  return ksstats(samp1,samp2)[5]
end

function calc_distance_ks_duration_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = summary1.stat["duration_ratio_list"]
  samp2 = summary2.stat["duration_ratio_list"]
  return ksstats(samp1,samp2)[5]
end

# Function for Relative Entropy / K-L divergence
include("kde.jl")

function calc_distance_kl_num_planets(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
    np1 = (haskey(summary1.stat,"num_tranets") ? summary1.stat["num_tranets"] : summary1.stat["expected planets detected"])
    np2 = (haskey(summary2.stat,"num_tranets") ? summary2.stat["num_tranets"] : summary2.stat["expected planets detected"])
    ntarg1 = summary1.stat["num targets"]
    ntarg2 = summary2.stat["num targets"]
    mintarg = min(ntarg1,ntarg2)
    alpha1 = 1+np1
    beta1 = 1+ntarg1
    alpha2 = 1+np2
    beta2 = 1+ntarg2
    #= Approximate Estimated Rate Posterior as a Normal Distribution
    mu1 = alpha1/beta1
    sigma1 = mu1 /sqrt(alpha1)
    mu2 = alpha2/beta2
    sigma2 = mu2 /sqrt(alpha1)
    dist = 0.5*( (sigma2/sigma1)^2 - 1 + 2*log(sigma1/sigma2) + (mu2-mu1)^2/sigma1^2 )  # Normal approximation, substitute denom for difference in rates
    =#
    #= Approximate Estimated Rate Posterior as a Gamma Distribution
    dist = abs( (alpha2-alpha1)*digamma(alpha1) - lgamma(alpha1) + lgamma(alpha2)  ) # Gamma approximation, same beta, why need abs?
    dist += abs( (alpha1-alpha2)*digamma(alpha2) - lgamma(alpha2) + lgamma(alpha1)  )  # Gamma approximation, same beta, why need abs?
    print("# a1=",alpha1, " a2=",alpha2, " digamma=",digamma(alpha1), " ", digamma(alpha2)," gamma term=",lgamma(alpha1)-lgamma(alpha2))
    if ntarg1!=ntarg2
       print("# ntarg=",ntarg1," ",ntarg2," dist(before betas)=", dist)
       dist += alpha2*(log(beta1)-log(beta2)) + alpha1*(beta2-beta1)/beta1  # Gammma approximation, beta terms
       dist += alpha1*(log(beta2)-log(beta1)) + alpha2*(beta1-beta2)/beta2  # Gammma approximation, beta terms
       println(" dist (after betas)=",dist)
    end
    =#
    ## #= Approximate Estimated Rate Posterior as a Poisson Distribution
    #dist = abs( alpha2-alpha1+alpha1*log(alpha1/alpha2) )
    #dist += abs( alpha2-alpha1+alpha1*log(alpha1/alpha2) )
    rate1 = alpha1/beta1*mintarg
    rate2 = alpha2/beta2*mintarg
    dist = (rate1-rate2)*log(rate1/rate2)
    ## =#
    dist /= 2
    return dist
end

function calc_distance_kl_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  #= categorical distribution
  f1sum = sum(summary1.stat["num n-tranet systems"]) # summary1.stat["num targets"]
  f2sum = sum(summary2.stat["num n-tranet systems"]) # summary2.stat["num targets"]
  if !(f1sum>0 && f2sum>0) return 0.0  end
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    f1 = summary1.stat["num n-tranet systems"][n]/f1sum
    f2 = summary2.stat["num n-tranet systems"][n]/f2sum
    m = (f1+f2)/2
    if m>zero(m)
    if f1>zero(f1)
       d[n] += 0.5*f1*log(f1/m)
    end
    if f2>zero(f2)
       d[n] += 0.5*f2*log(f2/m)
    end
#   else
#       d += Inf
    end
  end
  =#
  # Poisson distributions for each
    ntarg1 = summary1.stat["num targets"]
    ntarg2 = summary2.stat["num targets"]
    #mintarg = min(ntarg1,ntarg2)
  d = zeros(max_tranets_in_sys)
  for n in 1:max_tranets_in_sys
    f1 = (1+summary1.stat["num n-tranet systems"][n])/(1+ntarg1) # *(1+mintarg)
    f2 = (1+summary2.stat["num n-tranet systems"][n])/(1+ntarg2) # *(1+mintarg)
    d[n] = (f1-f2)*log(f1/f2)
  end
  return d
end

function calc_distance_hellinger_num_n_tranet_systems(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
  f1sum = sum(summary1.stat["num n-tranet systems"]) # summary1.stat["num targets"]
  f2sum = sum(summary2.stat["num n-tranet systems"]) # summary2.stat["num targets"]
  d = 1
  for n in 1:max_tranets_in_sys
    f1 = summary1.stat["num n-tranet systems"][n]/f1sum
    f2 = summary2.stat["num n-tranet systems"][n]/f2sum
    d -= sqrt(f1*f2)
  end
  #d = sqrt(d)
  return d
end

function calc_distance_kl_period(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  #println("# P list 1: n=",length(summary1.stat["P list"])," min=",minimum(summary1.stat["P list"]), " max=",maximum(summary1.stat["P list"]))
  #println("# P list 2: n=",length(summary2.stat["P list"])," min=",minimum(summary2.stat["P list"]), " max=",maximum(summary2.stat["P list"]))
  samp1 = log.(summary1.stat["P list"] )
  samp2 = log.(summary2.stat["P list"] )
  calc_kl_distance_ab(samp1,samp2,log(0.5),log(320.) )
end

function calc_distance_kl_depth(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  samp1 = log.(summary1.stat["depth list"] )
  samp2 = log.(summary2.stat["depth list"] )
  calc_kl_distance_ab(samp1,samp2,log(0.000025),log(0.025) )
end

function calc_distance_kl_period_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  min_ratios_to_compute_distances = 3
  distance_when_not_enough_ratios = 10.0
  samp1 = summary1.stat["period_ratio_list"]
  samp2 = summary2.stat["period_ratio_list"]
  if length(samp1)<min_ratios_to_compute_distances || length(samp2)<min_ratios_to_compute_distances
     return distance_when_not_enough_ratios
  end
  calc_kl_distance_ab(samp1,samp2,0.0,1.0)
end

function calc_distance_kl_duration_ratios(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, sim_param::SimParam ; verbose::Bool = false)
  min_ratios_to_compute_distances = 3
  distance_when_not_enough_ratios = 10.0
  samp1 = log.( summary1.stat["duration_ratio_list"] )
  samp2 = log.( summary2.stat["duration_ratio_list"] )
  if length(samp1)<min_ratios_to_compute_distances || length(samp2)<min_ratios_to_compute_distances
     return distance_when_not_enough_ratios
  end
  calc_kl_distance_ab(samp1,samp2,-log(10.0),log(10.0) )
end

## test functions in this file, requires model, summary stats and functions from ExoplanetsSysSim
function test_distance()
  sim_param = setup_sim_param_model()
  cat_phys1 = generate_kepler_physical_catalog(sim_param)
  cat_ref = observe_kepler_targets_single_obs(cat_phys1,sim_param)
  summary_stat_ref = calc_summary_stats_model(cat_ref,sim_param)
  cat_phys2 = generate_kepler_physical_catalog(sim_param)
  cat_obs = observe_kepler_targets_single_obs(cat_phys2,sim_param)
  summary_stat_obs = calc_summary_stats_model(cat_obs,sim_param)
  dist = calc_distance_model(summary_stat_ref,summary_stat_obs,sim_param)
end


text:
## LINEAR LEARNING NETWORKS (FOR INTERNAL USE ONLY)

# constructs and returns an unsupervised node, and its machine, for a
# given model and input. the node is static if model is a function
# instead of bona fide model. Example: `node_(pca, X)` (dynamic) or
# `node_(MLJ.matrix, X)` (static).
function node_(model, X)
    if model isa Model
        mach = machine(model, X)
        return transform(mach, X), mach
    else
        n = node(model, X)
        return n, n.machine
    end
end

# `models_and_functions` can include both functions (static
# operatations) and bona fide models. If `ys == nothing` the learning
# network is assumed to be unsupervised. Otherwise: If `target ===
# nothing`, then no target transform is applied; if `target !==
# nothing` and `inverse === nothing`, then corresponding `transform`
# and `inverse_transform` are applied; if neither `target` nor
# `inverse` are `nothing`, then both `target` and `inverse` are
# assumed to be StaticTransformations, to be applied
# respectively to `ys` and to the output of the `models_and_functions`
# pipeline. Note that target inversion is applied to the output of the
# *last* nodal machine in the pipeline, corresponding to the last
# element of `models_and_functions`.

# No checks whatsoever are performed. Returns a learning network.
function linear_learning_network(Xs, ys, target, inverse, models_and_functions...)

    n = length(models_and_functions)

    if ys !== nothing && target !== nothing
        yt, target_machine = node_(target, ys)
    else
        yt  = ys
    end

    nodes = Vector(undef, n + 1)
    nodes[1] = Xs

    for i = 2:(n + 1)
        m = models_and_functions[i-1]
        if m isa Supervised
            supervised_machine = machine(m, nodes[i-1], yt)
            nodes[i] = predict(supervised_machine, nodes[i-1])
       else
            nodes[i] = node_(m, nodes[i-1]) |> first
        end
    end

    if target === nothing
        terminal_node=nodes[end]
    else
        if inverse === nothing
            terminal_node = inverse_transform(target_machine, nodes[end])
        else
            terminal_node = node_(inverse, nodes[end]) |> first
        end
    end

    return terminal_node

end

## PREPROCESSING

function eval_and_reassign(modl, ex)
    s = gensym()
    evaluated = modl.eval(ex)
    modl.eval(:($s = $evaluated))
    return s, evaluated
end

pipe_alert(message) = throw(ArgumentError("@pipeline error.\n"*
                                     string(message)))
pipe_alert(k::Int) = throw(ArgumentError("@pipeline error $k. "))

# does expression processing, syntax and semantic
# checks. is_probabilistic is `true`, `false` or `missing`.
function pipeline_preprocess(modl, ex, is_probabilistic::Union{Missing,Bool})

    ex isa Expr || pipe_alert(2)
    length(ex.args) > 1 || pipe_alert(3)
    ex.head == :call || pipe_alert(4)
    pipetype_ = ex.args[1]
    pipetype_ isa Symbol || pipe_alert(5)

    fieldnames_ = []               # for `from_network`
    models_ = []                   # for `from_network`
    models_and_functions_ = []     # for `linear_learning_network`
    models_and_functions  = []     # for `linear_learning_network`
    for ex in ex.args[2:end]
        if ex isa Expr
            if ex.head == :kw
                variable_ = ex.args[1]
                variable_ isa Symbol || pipe_alert(8)
                value_, value = eval_and_reassign(modl, ex.args[2])
                if variable_ == :target
                    if value isa Function
                        value_, value =
                            eval_and_reassign(modl,
                                              :(MLJ.StaticTransformer($value)))
                    end
                    value isa Unsupervised ||
                        pipe_alert("Got $value where a function or "*
                                   "Unsupervised model instance was "*
                                   "expected. ")
                    target_ = value_
                    target = value
                    push!(fieldnames_, :target)
                elseif variable_ == :inverse
                    if value isa Function
                        value_, value =
                            eval_and_reassign(modl,
                                              :(MLJ.StaticTransformer($value)))
                    else
                        pipe_alert(10)
                    end
                    inverse_ = value_
                    inverse = value
                    push!(fieldnames_, :inverse)
                else
                    value isa Model || pipe_alert(11)
                    push!(models_and_functions_, value_)
                    push!(models_and_functions, value)
                    push!(fieldnames_, variable_)
                end
                push!(models_, value_)
            else
                f = modl.eval(ex)
                f isa Function ||
                    pipe_alert("Perhaps a missing name, as `name=$f'?")
                push!(models_and_functions_, ex)
                push!(models_and_functions, f)
            end
        else
            f = modl.eval(ex)
            f isa Function || pipe_alert(7)
            push!(models_and_functions_, ex)
            push!(models_and_functions, f)
        end
    end

    (@isdefined target)  || (target = nothing; target_ = :nothing)
    (@isdefined inverse) || (inverse = nothing; inverse_ = :nothing)
    inverse !== nothing && target === nothing &&
        pipe_alert("You have specified `inverse=...` but no `target`.")

    supervised(m) = !(m isa Unsupervised) && !(m isa Function)

    supervised_components = filter(supervised, models_and_functions)

    length(supervised_components) < 2 ||
        pipe_alert("More than one component of the pipeline is a "*
                   "supervised model .")

    is_supervised  =
        length(supervised_components) == 1

    # `kind_` is defined in composites.jl:
    kind = kind_(is_supervised, is_probabilistic)
    ismissing(kind) &&
        pipe_alert("Network has no supervised components and so "*
                  "`is_probabilistic=true` "*
                  "declaration is not allowed. ")

    target isa StaticTransformer && inverse == nothing &&
        pipe_alert("It appears `target` is a function. "*
                   "You must therefore specify `inverse=...` .")

    target == nothing || is_supervised ||
        pipe_alert("`target=...` has been specified but no "*
                   "supervised components have been specified. ")

    return (pipetype_, fieldnames_, models_,
            models_and_functions_, target_, inverse_, kind)

end

function pipeline_preprocess(modl, ex, kw_ex)
    kw_ex isa Expr || pipe_alert(10)
    kw_ex.head == :(=) || pipe_alert(11)
    kw_ex.args[1] == :is_probabilistic ||
        pipe_alert("Unrecognized keywork `$(kw_ex.args[1])`.")
    value_ = kw_ex.args[2]
    if value_ == :missing
        value = missing
    else
        value = value_
        value isa Bool ||
            pipe_alert("`is_probabilistic` can only be `true` or `false`.")
    end
    return pipeline_preprocess(modl, ex, value)
end

function pipeline_(modl, ex, kw_ex)

    (pipetype_, fieldnames_, models_, models_and_functions_,
     target_, inverse_, kind) =
         pipeline_preprocess(modl, ex, kw_ex)

    if kind === :UnsupervisedNetwork
        ys_ = :nothing
    else
        ys_ = :(source(nothing, kind=:target))
    end

    N_ex = quote

        MLJ.linear_learning_network(source(nothing), $ys_,
                                    $target_, $inverse_,
                                    $(models_and_functions_...))
    end

    from_network_(modl, pipetype_, fieldnames_, models_, N_ex, kind)

    return pipetype_

end

pipeline_(modl, ex) = pipeline_(modl, ex, :(is_probabilistic=missing))

"""
    @pipeline NewPipeType(fld1=model1, fld2=model2, ...)
    @pipeline NewPipeType(fld1=model1, fld2=model2, ...) is_probabilistic=false

Create a new pipeline model type `NewPipeType` that composes the types of
the specified models `model1`, `model2`, ... . The models are composed
in the specified order, meaning the input(s) of the pipeline goes to
`model1`, whose output is sent to `model2`, and so forth.

At most one of the models may be a supervised model, in which case
`NewPipeType` is supervised. Otherwise it is unsupervised.

The new model type `NewPipeType` has hyperparameters (fields) named
`:fld1`, `:fld2`, ..., whose default values for an automatically
generated keyword constructor are deep copies of `model1`, `model2`,
... .

*Important.* If the overall pipeline is supervised and makes probabilistic
predictions, then one must declare `is_probabilistic=true`. In the
deterministic case the keyword argument can be omitted.

Static (unlearned) transformations - that is, ordinary functions - may
also be inserted in the pipeline as shown in the following example
(the classifier is probabilistic but the pipeline itself is
deterministic):

    @pipeline MyPipe(X -> coerce(X, :age=>Continuous),
                     hot=OneHotEncoder(),
                     cnst=ConstantClassifier(),
                     yhat -> mode.(yhat))

### Return value

An instance of the new type, with default
hyperparameters (see above), is returned.

### Target transformation and inverse transformation

A learned target transformation (such as standardization) can also be
specified, using the keyword `target`, provided the transformer
provides an `inverse_transform` method:

    @load KNNRegressor
    @pipeline MyPipe(hot=OneHotEncoder(),
                     knn=KNNRegressor(),
                     target=UnivariateTransformer())

A static transformation can be specified instead, but then
an `inverse` must also be given:

    @load KNNRegressor
    @pipeline MyPipe(hot=OneHotEncoder(),
                     knn=KNNRegressor(),
                     target = v -> log.(v),
                     inverse = v -> exp.(v))

*Important.* While the supervised model in a pipeline providing a
 target transformation can appear anywhere in the pipeline (as in
 `ConstantClassifier` example above), the inverse operation is always
 performed on the output of the *final* model or static
 transformation in the pipeline.

See also: [@from_network](@ref)

"""
macro pipeline(exs...)
    pipetype_ = pipeline_(__module__, exs...)

    esc(quote
        $pipetype_()
        end)

end


text:
# Autogenerated wrapper script for czmq_jll for i686-w64-mingw32
export libczmq

using ZeroMQ_jll
JLLWrappers.@generate_wrapper_header("czmq")
JLLWrappers.@declare_library_product(libczmq, "libczmq.dll")
function __init__()
    JLLWrappers.@generate_init_header(ZeroMQ_jll)
    JLLWrappers.@init_library_product(
        libczmq,
        "bin\\libczmq.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
##Reineke SDI chart

using RecipesBase

@userplot sdi_chart
@recipe function f(dat::sdi_chart ;maxsdi=450)

diarng = 1:1:20
#tparng = [1,1200] #unusued for now, TODO: figure out how to autoscale plot?
maxline() = maxsdi ./ (diarng ./ 10.0).^1.605

  # plot attributes
  xlabel --> "Trees Per Acre"
  ylabel --> "Quadratic Mean Diameter"
  xticks --> (0:50:450)
  yticks --> (0:20:200)
  xscale --> :log10
  yscale --> :log10

  legend := false
  primary := false

#max SDI line
@series begin
    linecolor := :black
    maxline()
end

#competition induced mortality SDI line
@series begin
    linecolor := :orange
    maxline() * 0.55
end

#crown closure SDI line
@series begin
    linecolor := :red
    maxline() * 0.35
end

 # the point and the annotations
 @series begin
    primary := true
    seriestype := :scatter
    #annotations := (100:38:442, 20:10:110, ["$i%" for i in 20:10:110])

    a,b = dat.args
    [a], [b]
 end
end


text:
__precompile__(true)
module Roots

using LinearAlgebra

if VERSION >= v"0.7-"
    using Printf
else
    using Missings
end

using ForwardDiff
using Compat: @nospecialize

export fzero,
       fzeros,
       newton, halley,  # deprecate these 4?
       secant_method, steffensen,
       D

export find_zero,
       Order0, Order1, Order2, Order5, Order8, Order16
export Bisection, FalsePosition

## load in files
include("adiff.jl")
include("find_zero.jl")
include("bracketing.jl")
include("derivative_free.jl")
include("newton.jl")

## Main functions are
## * fzero(f, ...) to find _a_ zero of f
## * fzeros(f, ...) to attempt to find all zeros of f

"""

Find zero of a function using an iterative algorithm

* `f`: a scalar function or callable object
* `x0`: an initial guess, finite valued.

Keyword arguments:

* `ftol`: tolerance for a guess `abs(f(x)) < ftol`
* `xtol`: stop if `abs(xold - xnew) <= xtol + max(1, |xnew|)*xtolrel`
* `xtolrel`: see `xtol`
* `maxeval`: maximum number of steps
* `verbose`: Boolean. Set `true` to trace algorithm
* `order`: Can specify order of algorithm. 0 is most robust, also 1, 2, 5, 8, 16.
* `kwargs...` passed on to different algorithms. There are `maxfneval` when `order` is 1,2,5,8, or 16 and `beta` for orders 2,5,8,16,

This is a polyalgorithm redirecting different algorithms based on the value of `order`.

"""
function fzero(f, x0::Real; kwargs...)
    x = float(x0)
    isinf(x) && throw(ConvergenceFailed("An initial value must be finite"))
    derivative_free(f, x; kwargs...)
end

"""
Find zero of a function within a bracket

Uses a modified bisection method for non `big` arguments

Arguments:

* `f` A scalar function or callable object
* `a` left endpont of interval
* `b` right endpont of interval
* `xtol` optional additional tolerance on sub-bracket size.

For a bracket to be valid, it must be that `f(a)*f(b) < 0`.

For `Float64` values, the answer is such that `f(prevfloat(x)) *
f(nextfloat(x)) < 0` unless a non-zero value of `xtol` is specified in
which case, it stops when the sub-bracketing produces an bracket with
length less than `max(xtol, abs(x1)*xtolrel)`.

For `Big` values, which defaults to the algorithm of Alefeld, Potra, and Shi, a
default tolerance is used for the sub-bracket length that can be enlarged with `xtol`.

If `a==-Inf` it is replaced with `nextfloat(a)`; if `b==Inf` it is replaced with `prevfloat(b)`.

Example:

    `fzero(sin, 3, 4)` # find pi
    `fzero(sin, [big(3), 4]) find pi with more digits
"""
function fzero(f, a::Real, b::Real; kwargs...)
    bracket = adjust_bracket((a, b))
    a0, b0 = a < b ? promote(float(a), b) : promote(float(b), a)
    (a0 == -Inf) && (a0 = nextfloat(a0))
    (b0 == Inf) && (b0 = prevfloat(b0))
    find_zero(f, bracket, Bisection(); kwargs...)
end

"""
Find a zero with bracket specified via `[a,b]`, as `fzero(sin, [3,4])`.
"""
function fzero(f, bracket::Vector{T}; kwargs...)  where {T <: Real}
    fzero(f, float(bracket[1]), float(bracket[2]); kwargs...)
end

"""
Find a zero within a bracket with an initial guess to *possibly* speed things along.
"""
function fzero(f, x0::Real, bracket::Vector{T}; kwargs...)  where {T <: Real}

    a, b = adjust_bracket(bracket)

    try
        ex = a42a(f, a, b, float(x0); kwargs...)
    catch ex
        if isa(ex, StateConverged)
            return(ex.x0)
        else
            rethrow(ex)
        end
    end
end

"""
Find zero using Newton's method.
"""
fzero(f::Function, fp::Function, x0::Real; kwargs...) = newton(f, fp, float(x0); kwargs...)

## fzeros
"""

`fzeros(f, a, b)`

Attempt to find all zeros of `f` within an interval `[a,b]`.

Simple algorithm that splits `[a,b]` into subintervals and checks each
for a root.  For bracketing subintervals, bisection is
used. Otherwise, a derivative-free method is used. If there are a
large number of zeros found relative to the number of subintervals, the
number of subintervals is increased and the process is re-run.

There are possible issues with close-by zeros and zeros which do not
cross the origin (non-simple zeros). Answers should be confirmed
graphically, if possible.

"""
function fzeros(f, a::Real, b::Real; kwargs...)
    find_zeros(f, float(a), float(b); kwargs...)
end
fzeros(f, bracket::Vector{T}; kwargs...) where {T <: Real} = fzeros(f, bracket[1], bracket[2]; kwargs...)

end